question5.R

Nathan — Mar 7, 2014, 10:45 PM

#### Variables ####
x = 50
delta = 0.05

#### Functions ####
integrate <- function(f, a, b, n=1000) {
  l = a + ((b - a)/n) * 0.5
  A = 0

  for (i in 0:(n - 1)) {
    A = A + f(l) * ((b - a)/n)
    l = l + ((b - a)/n)
  }

  return(A)
}

v <- function(delta, n) {
  return(exp(-delta*n))
}

ux <- function(x,t) {
  if(x<= 65) return (0.008)
  else return(0.020)
}

tpx <- function(t,x) {
  return(exp(-ux(x)*t))
}

pdf <- function(t,x) {
  return(tpx(t,x)*ux(x))
}

# 15-year term 
Z1 <- function(t) {
  if (t<=15) return(v(delta,t))
  else return(0)
}

# 15-year deferred whole life
Z2 <- function(t) {
  if (t>=15) return(v(delta,t))
  else return (0)
}

# whole life
Z3 <- function(t) {
  return(v(delta,t))
}

# 15 year pure endowment
Z4 <- function(t) {
  if(t>=15) return(v(delta,15))
  else return(0)
}

# annually increasing wholfe life
Z5 <- function(t) {
  if(t<=15) {
    return(floor(t+1)*v(delta,t))
  } else {
    return(15*v(delta,t))
  }
}

EZ1 <- function(moment=1) {
  delta <<- delta*moment 

  f <- function(t,x=50) return(pdf(t,x)*Z1(t))
  ans = round(integrate(f,0,1000),5)

  delta <<- delta/moment
  return(ans)
}

EZ2 <- function(moment=1) {
  delta <<- delta*moment 

  f <- function(t,x=50) return(pdf(t,x)*Z2(t))
  ans = round(integrate(f,0,1000),5)

  delta <<- delta/moment
  return(ans)
}

EZ3 <- function(moment=1) {
  delta <<- delta*moment 

  f <- function(t,x=50) return(pdf(t,x)*Z3(t))
  ans = round(integrate(f,0,1000),5)

  delta <<- delta/moment
  return(ans)
}

EZ4 <- function(moment=1) {
  delta <<- delta*moment 

  # f <- function(t,x=50) return(pdf(t,x)*Z4(t))
  # ans = round( integrate(f,0,1000),5)
  ans = round(tpx(15,50)*v(delta,15),5)

  delta <<- delta/moment
  return(ans)
}

EZ5 <- function(moment=1) {
  temp <- delta; delta <<- delta*moment 

  f <- function(t,x=50) return(pdf(t,x)*Z5(t))
  ans = round( (integrate(f,0,15) + 15*EZ2()), 5)
  #ans = round ( (integrate(f,0,1000)),5)
  delta <<- temp
  return(ans)
}

# Variance formula for Z1, Z2, Z3, and Z4
VZ <- function(f) {
  return(f(moment=2)- (f(moment=1))^2) 
}

VZ5 <- function(f) {
  f <- function(t,x=50) return(pdf(t,x)*Z5(t))

  additionalTerm <- function(t,x=50) {
    intFunction <- function(t,x=50) {
      return (ux(x)*floor(t+1)^2 * exp((-2*delta + ux(x))*t))
    }
    return(integrate(intFunction,0,15) - integrate(f,0,15)^2)
  }

  return( 225*VZ(EZ2) - 2*(15*EZ2()*integrate(f,0,15)) + additionalTerm())
}

q5 <- function() {
  Zi = c(1,2,3,4,5)
  EZi = c(EZ1(),EZ2(),EZ3(),EZ4(),EZ5())
  VZi = c(VZ(EZ1),VZ(EZ2),VZ(EZ3),VZ(EZ4),VZ5())
  return(data.frame(Zi,EZi,VZi))
}

ans = q5()
print(q5())

  Zi     EZi     VZi
1  1 0.08013 0.05297
2  2 0.05778 0.01131
3  3 0.13791 0.05502
4  4 0.41895 0.02238
5  5 1.42216 4.99252