assign4.R

Nathan — Apr 7, 2014, 12:19 PM

rm(list=ls()); cat("\014"); options(warn=-1)

graphics.off(); plotDims=c(2,2); par(mfrow=plotDims)

#### vars ####
A = 0.00022
B = 2.5e-05
c = 1.1
l20 = 1e+05
x = 20
i = 0.05
w = 120

writeToExcel = FALSE
if(writeToExcel==TRUE){tryCatch({library(xlsx);library(rJava)},error=function(e){install.packages("xlsx")})}

#### functions ####
fix <- function(vector, k, truncate) {
  decimals <- function(x,k) {
    format(round(x,k), nsmall=k)
  }
  newvector <- numeric(0); vector = as.numeric(vector)

  for (j in 1:length(vector)) newvector[j] = decimals(vector[j],k)

  newvector = newvector[1:truncate]
  return(as.numeric(newvector))
}

display <- function(message, num) {
  cat(paste(message, num, sep="")); cat("\n")
}

uxplust <- function(t,...){ 
  age = if(length(list(...))==1) list(...)[[1]] else x
  return(A+(B*c^age)*c^t)
}

uxpluss <- function(t,...) {
  if(t<0 || t>2) return(uxplust(t,...))
  else return(0.9^(2-t)*uxplust(t,...))
}

fxt <- function(t,...) {
  return(tpx(t,...)*uxplust(t,...))
}

tpx <- function(t,...) {
  age = if(length(list(...))==1) list(...)[[1]] else x

  # account for the fact that no one lives beyond age w
  if( (age+t)>= w) return(0)
  return(exp(-(A*t + B/log(c)*(c^age)*(c^t -1))))
}

tpselectx <- function(t,age=x) {
  if( (age+t)>= w) return(0)
  if (0<= t & t<=2) {
    return(exp(0.9^(2-t)* ((A*(1 - 0.9^t)/log(0.9)) + (B*c^age*(c^t - 0.9^t)/log(0.9/c)))))
  } else {
    #return(tpx(t,age))
    return(tpselectx(2,age)*tpx(t-2,age+2))
  }
}


tpselectxplus1 <- function(t,...) {
  age = if(length(list(...))==1) list(...)[[1]] else x

  # account for the act that no one lives beyond age w
  if ( (age+t)>=w) return(0)
  if(0<=t & t<=1) {
    return(exp(0.9^(1-t)*( (A/log(0.9)*(1-0.9^t)) + (B*c^(age+1)/log(0.9/c)*(c^t - 0.9^t)))))
  } else {
    return(tpx(t,...))
  }
}

tqx <- function(t,...) {
  return(1-tpx(t,...))
}

tqselectx <- function(t,...) {
  return(1-tpselectx(t,...))
}

tqselectxplus1 <- function(t,...) {
  return(1-tpselectxplus1(t,...))
}

kdeferredqx <- function(k,...) {
  return(tpx(k,...) - tpx(k+1,...))
}

kdeferredselectqx <- function(k,...) {
  return(tpselectx(k,...) - tpselectx(k+1,...))
}

v <- function(i, n, ...) {
  moment = if(length(list(...)) == 1) list(...)[[1]] else 1
  delta = log(1+i)*moment
  return(exp(-delta*n))
}

Ax <- function(x,...) {
  total = 0
  # w-x-1 since w-x would make A120 non-zero
  # sum from k = 0 to w-x-1 of v^k+1 * k|qx 
  #however, (w-x-1) = (120-20-1) = 99
  for (k in (0:(w-x-1))) {
    total = total + v(i, (k+1),...)*kdeferredqx(k,x)
  }
  return(total)
}

Axcolumn <- function(x,w,i,...) {
  axcol <- numeric(0); temp = 0
  # w-x-1 since w-x would make A120 non-zero
  for (k in 0:(w-x)) {
    if(k==0) temp = v(i,1,...)
    else temp = axcol[k]*tpx(1,(w-k))*v(i,1,...) + v(i,1,...)*tqx(1,(w-k))
    axcol = c(axcol, temp)
  }

  # reverse the vector due to backwards recursion
  return(rev(axcol))
}

Aselectx <- function(x,...) {
  return(v(i,1,...)*tqselectx(1,x) + v(i,2,...)*tqselectxplus1(1,x)*tpselectx(1,x) + 
           v(i,2,...)*tpselectx(2,x)*(Ax(x+2,...)))
}

Aselectxcolumn <- function(x,w,i,...) {
  aselectxcol <- numeric(0); temp = 0
  axplus2col = Axcolumn(x+2, w+2, i, ...)

  for (k in 0:(w-x)) {
    # temp = Aselectx(k+x,...)
    temp = (v(i,1,...)*tqselectx(1,x+k) + v(i, 2, ...)*tqselectxplus1(1,x+k)*tpselectx(1,x+k) 
            + v(i,2,...)*tpselectx(2,x+k)*axplus2col[k+1])
    aselectxcol = c(aselectxcol, temp)
  }
  return(aselectxcol)
}

Aselectxplots <- function(i1,i2,i3,xmin,xmax,ymin,ymax,legend) {
  temp = i
  i <<- i1; plot( Aselectx, xmin, xmax, xlab="Age at Selection",
                  ylab = expression("A"["[x]"]), col="red", 
                  main = expression(bold(paste("A"["[x]"], " with varying interest rates"))),
                  ylim=c(ymin,ymax))
  i <<- i2; plot( Aselectx, xmin, xmax, col="green", add=TRUE)
  i <<- i3; plot( Aselectx, xmin, xmax, col="blue", add=TRUE)
  if (legend==TRUE) {
    legend(xmin, ymax, c(paste(i1*100,"%"), paste(i2*100,"%"), paste(i3*100,"%")),
           lty=c(1,1,1), lwd=c(2,2,2), col=c("red", "green", "blue"), cex=0.7)
  }
  i <<- temp   
}

# use select probabilities
pdfZ <- function(k,x) {
  if(0<k & k<(w-1)) return(kdeferredselectqx(k-1,x))
  return(0)
}

# x, l[x], l[x]+1, lx+2, x+2
createtable1 <- function(x, w, truncateTable) {
  age <- numeric(0)
  age_plus_2 <- numeric(0)
  l_x <- numeric(0)
  l_x_1 <- numeric(0)
  lx_2 <- numeric(0)
  truncateAt = truncateTable-x+1

  for (j in ((x-2):(w-2))) {
    lx_2 = c(lx_2, l20*tpx((j-(x-2)),x))
    age_plus_2 = c(age_plus_2, (j+2))
    # account for blanks in first two rows
    if(j<x) { l_x=c(l_x,"-"); l_x_1 = c(l_x_1, "-"); age=c(age,"-") }
    else {
      l_x = c(l_x, lx_2[(j-(x-3))]/tpselectx(2,j))
      l_x_1 = c(l_x_1, lx_2[(j-(x-3))]/tpselectxplus1(1,j))
      age = c(age,j)
    }
  }
  lx_2 = fix(lx_2, 2, truncateAt)
  age_plus_2 = fix(age_plus_2, 0, truncateAt)
  age = fix(age, 0, truncateAt)
  l_x = fix(l_x,2, truncateAt)
  l_x_1 = fix(l_x_1, 2, truncateAt)
  table <- data.frame(age, l_x, l_x_1, lx_2, age_plus_2)

  names(table)[names(table)=="age"] = "x"
  names(table)[names(table)=="age_plus_2"] = "x+2"
  names(table)[names(table)=="l_x"] = "l[x]"
  names(table)[names(table)=="l_x_1"] = "l[x]+1"
  names(table)[names(table)=="lx_2"] = "lx+2"
  return(table)
}

# x, A[x], 2A[x]
createtable2 <- function(x,w, truncateTable) {
  age <- numeric(0)
  aselectxcol = Aselectxcolumn(x,w,i,1)
  aselectx2col = Aselectxcolumn(x,w,i,2)
  truncateAt = truncateTable-x+1

  for (i in (x:w)) age = c(age,i)

  aselectxcol = fix(aselectxcol, 5, truncateAt)
  aselectx2col = fix(aselectx2col, 5, truncateAt)
  age = fix(age, 0, truncateAt)
  table <- data.frame(age, aselectxcol, aselectx2col)

  names(table)[names(table)=="age"] = "x"
  names(table)[names(table)=="aselectxcol"] = "A[x]"
  names(table)[names(table)=="aselectx2col"] = "2A[x]"

  return(table)
}

# x, Ax, 2Ax 
createtable3 <- function(x,w,truncateTable) {
  age <- numeric(0)

  axcol = Axcolumn(x,w,i,1)
  ax2col = Axcolumn(x,w,i,2)
  # if table goes from 20 to 80 then truncate at index (80-20)+1
  # this is because in R, array indexices start from 1 not 0
  truncateAt = truncateTable-x+1

  for (i in (x:w)) age = c(age,i)

  axcol = fix(axcol, 5, truncateAt)
  ax2col = fix(ax2col, 5, truncateAt)
  age = fix(age,0, truncateAt)
  table <- data.frame(age, axcol, ax2col)
  names(table)[names(table)=="age"] = "x"
  names(table)[names(table)=="axcol"] = "Ax"
  names(table)[names(table)=="ax2col"] = "2Ax"
  return(table)
}

# k, v^k, P(Z=v^k, k|qx
# use select probabilities
# slight error on the last few cdf values due to rounding to 5 decimal places
createtable4 <- function(x,w,truncateTable) {
  kdeferredqx_ <- numeric(0)
  kx <- numeric(0)
  v_k <- numeric(0)
  zv_k <- numeric(0)
  truncateAt = truncateTable-x+1

  for (k in 0:(w-x)) {
    kdeferredqx_ = c(kdeferredqx_, kdeferredselectqx(k,x))
    kx = c(kx,k)
    v_k = c(v_k, v(i,k))
    # probability that Z = v(i,k)
    zv_k = c(zv_k, pdfZ(k,x))
  }

  kx = fix(kx, 0, truncateAt)
  kdeferredqx_ = fix(kdeferredqx_, 5, truncateAt)
  v_k = fix(v_k, 5, truncateAt)
  zv_k = fix(zv_k, 5, truncateAt)

  cdf_k <- numeric(0); sdf_Z <- numeric(0)
  total = 0
  # find cdf of k
  for (j in 1:(length(kdeferredqx_))) {
    total = total+kdeferredqx_[j]
    cdf_k = c(cdf_k, total)
    sdf_Z = c(sdf_Z, total-kdeferredqx_[j])
  }
  cdf_k = fix(cdf_k, 5, truncateAt); sdf_Z = fix(sdf_Z, 5, truncateAt)
  table = data.frame(kx, kdeferredqx_, cdf_k, v_k, zv_k, sdf_Z)

  names(table)[names(table)=="kdeferredqx_"] = "P[K(50)=k]"
  names(table)[names(table)=="kx"] = "K(x)"
  names(table)[names(table)=="v_k"] = "v^k"
  names(table)[names(table)=="zv_k"] = "P[Z=v^k]"
  names(table)[names(table)=="cdf_k"] = "P[K(x)<=k]"
  names(table)[names(table)=="sdf_Z"] = "P[Z>=v^k]"

  theoreticalZ <<- as.numeric(zv_k)
  cdfK <<- as.numeric(cdf_k)
  sdfZ <<- as.numeric(sdf_Z)
  return(table)
}

# draw n observations of K(50) 
# x: age, w: limiting age
kNums <- function(n, w,x, cdf) {
  # use cdf vector to transform from runif to K(x)
  rand <- function(n) return(runif(n,0,1))
  nums = rand(n)

  transform <- function(num) {
    k = 0
    if(num>cdf[length(cdf)]) num = runif(1,0,1)
    while(num>cdf[k+1] & k<69) k=k+1
    return(k)
  }

  knums <- numeric(0)
  for(num in nums) knums = c(knums, transform(num))
  return(knums)
}

# draw n observations of Z corresponding to K(50)
zNums <- function(nums,...) {
  j = if(length(list(...))==1) list(...)[[1]] else 0.05
  newnums <- numeric(0)
  for (k in 1:length(nums)) newnums = c(newnums, v(j, nums[k]+1))
  return(newnums)
}

# 1 - cdf(Z)
sdf <- function(z, sdfZ) {
  k = 0
  while(v(i,k+1)>z & k<200) k=k+1
  return(sdfZ[k+1])
}

#### output ####

# for the first (3 tables and 2 plots) x = 20, and is 50 otherwise
table1 = createtable1(x,w,w-1)  # goes up to 119
table2 = createtable2(x,w,w-3) 
table3 = createtable3(x,w,w-1)
table4 = createtable4(50,w,w-1)

Aselectxplots(0.04,0.05,0.06,x,117,0,1,TRUE) # plot goes up to 117
plot(theoreticalZ, main="Theoretical Distribution of Z", ylab="P(Z)=z)", xlab="Z", col="blue")

# n = 50, n = 500
testCases = c(50,500)
MeanZ = Aselectx(50); SdZ = sqrt(Aselectx(50,2) - Aselectx(50,1)^2)

# P(Z>=A[50])
prob = sdf(MeanZ,sdfZ)
cat("\n"); display("P(Z>= A[50]): ", prob); cat("\n")

P(Z>= A[50]): 0.34213

for (j in 1:length(testCases)) {
  simKx = kNums(testCases[j],50,w,cdfK)
  simZ = zNums(simKx, 0.05)
  hist(simZ, main=paste("Histogram of Z (n = ", testCases[j], ")"), xlab ="Z", col="pink")
  simMeanZ = mean(simZ); simSdZ = sd(simZ)
  display("n is: ", testCases[j])
  display("MeanZ: ", MeanZ); display("simMeanZ: ", simMeanZ)
  display("SdZ: ", SdZ); display("simSdZ: ", simSdZ)
  simCdfZ = ecdf(simZ)
  display("P(Z>= A[50]: ~", (1-simCdfZ(MeanZ))); cat("\n")
}

n is: 50
MeanZ: 0.246181759757277
simMeanZ: 0.225124598643679
SdZ: 0.160387268174681
simSdZ: 0.142766170803151
P(Z>= A[50]: ~0.3

n is: 500
MeanZ: 0.246181759757277
simMeanZ: 0.233912423572489
SdZ: 0.160387268174681
simSdZ: 0.162658173920116
P(Z>= A[50]: ~0.294