ultimate_survival_model.R

global = list(A=0.00022, B=2.7e-06, c=1.124, d=2, x=20, w=131, radix=1e+05, i=0.05) # d is the select period

# i: interest rate, m: compounding
d <- function(i=global$i, m=1) {
  m*(1-(1-i/(1+i))^(1/m))
}

# i: interest rate, n: length, delta: force of interest
v <- function(i=global$i, n=1, delta=log(1+i)) {
  exp(-delta*n)
}

# s: select remaining, x: age, d: total select, A: makeham param, B: makeham param, c: makeham param
uxs <- function(s, x=global$x, d=global$d, A=global$A, B=global$B, c=global$c) { 
  0.9^(pmax(d-s,0))*(A+B*c^(x+s))
}

# t: survival period, x: age, s: remaining select, d: total select, w: limiting age
tpx <- function(t, x=global$x, s=0, d=global$d, w=global$w) { 
  unlist(ifelse((x+t)>=w, 0, lapply(t, function(t) exp(-integrate(function(t) uxs(t,ifelse(s!=d, x+s,x),s), 0, 
                        min(t,d))$value-integrate(function(t) uxs(t,ifelse(s!=d, x+s+d, x+d),0), 0, max(0, t-d))$value))))
}

# t: survival period, x: age, s: remaining select, d: total select, w: limiting age
tqx <- function(t,x=global$x, s=0, d=global$d, w=global$w)  {
  1 - tpx(t,x,s,d,w) 
}

# u: deferral, t: years after deferral, x: age, s: select, d: total select, w: limiting age
udeferredtqx <- function(u, t=1, x=global$x, s=0, d=global$d, w=global$w) {
  unlist(lapply(u, function(u) tpx(u, x, s, d, w) - tpx(u + t, x, s, d, w)))
}

# x: age, s: remaining select, d: total select, w: limiting age, i: interest rate, m: compounding, n: term, c: continuous, e: endowment, mt: moment
Ax <- function(x=global$x, s=0, d=global$d, w=global$w, i=global$i, m=1, n=w-x, c=0, e=0, mt=1) {
  Ax. <- function(x=global$x, s=0, d=global$d, w=global$w, i=global$i, m=1, n=w-x, c=0, e=0) {
    tm=ifelse(c, unlist(lapply(x, function(x) integrate(function(t) tpx(t,x,s,d,w)*uxs(t+d-s,x-d+s,d)*v(i,t,delta=log(1+i)*mt), 0, n)$value)),
         unlist(lapply(x, function(x) sum(udeferredtqx(seq(0,m*n-1)/m,1/m,x,s,d)*v(i,(seq(0,m*n-1)+1)/m, delta=log(1+i)*mt)))))
    ifelse(e,tm+tpx(n,x,s,d,w)*v(i,n,delta=log(1+i)*mt),tm)
  }
  unlist(lapply(n, function(n) Ax.(x,s,d,w,i,m,n,c,e)))
}

# x: age, s: remaining select, d: total select, w: limiting age, i: interest rate
ax <- function(x=global$x, s=0, d=global$d, w=global$w, i=global$i) {
  (1 - Ax(x,s,d,w,i))/d(i)
} 

# x: age, w: limiting age, radix: lx, d: total select
createLifeTable <- function(x=global$x, w=global$w, radix=global$radix, d=global$d) {
  lt = data.frame(c(rep(NA,d),x:(w-d-1)), lapply(0:(d-1), function(y) c(rep(NA,d), 
                  tail(tpx(0:(w-x-1),x)*radix,-d)/unlist(lapply(x:(w-x-1), 
                  function(x) tpx(d-y,x,d-y))))), tpx(0:(w-x-1),x)*radix, x:(w-1))
  names(lt) = c("x", unlist(lapply(0:(d-1), function(x) paste("l[x]+",x,sep=""))), paste("lx+",d,sep=""),"x+2"); lt
}