nathan — May 6, 2014, 1:45 PM
# STANDARD SELECT AND ULTIMATE SURVIVAL MODEL
# (makeham's law, ux = A + B*c^x, with a 2 year select period
# see DHW Actuarial Mathematics for Life Contingent Risks
# if tp[x], tp[x]+1, tpx are overidden, the survival model will be changed (feel free to do this)
# note that if tp[x] = tp[x]+1 = tpx, then there is no select period
# specific insurance and annuity products can be created by combining implemented functions
#### Variables ####
A=0.00022
B=2.7e-06
c=1.124
x=20
w=131
m=12
i=0.05
l20 = 100000
#### Main Functions ####
#' @title Annuity Due
#' @description present value of n-year annuity due
#' @param i the interest rate
#' @param n the number of years
#' @details uses sum of geometric series formula for computation
#' @examples adueanglen(0.05,5)
#' @export
adueanglen <- function(i,n) {
return( (1-v(i,n))/i)
}
#' @title Annuity Immediate
#' @description present value of n-year annuity immediate
#' @param i the interest rate
#' @param n the number of years
#' @details uses sum of geometric series formula for computation. differs by (1+i) factor from adueanglen(i,n) by which it is smaller since payments are made later at the end of the year
#' @examples aanglen(0.05,5)
#' @export
aanglen <- function(i,n) {
return( (1-v(i,n))/i)
}
#' @title Annuity Due (mthly)
#' @description present value of n-year annuity due (mthly)
#' @param i the interest rate
#' @param n the number of periods
#' @param m the number of times interest is compounded per periods
#' @details only different from Annuity Due formula is denominator
#' @examples adueanglen(0.05,5,12)
#' @export
adueanglenupperm <- function(i,n,m) {
return((1-v(i,n))/dupperm(i,m))
}
#' @title Whole life annuity due (select life aged x)
#' @description Whole life annuity due is a contingent cash flow. Only X param is required (uses Standard Select and Ulimate Survival Model)
#' @param x the age of the select life
#' @details Contingent whole life annuity - payments are made at beginning of year and stop when x dies. survivorship and interest is accounted for
#' @examples adueselectx(4)
#' @references DHW Actuarial Mathematics for Life Contingent Risks (2nd)
#' @export
adueselectx <- function(x) {
return (( 1 - Aselectx(x))/ d(i))
}
#' @title Whole life annuity due (select life aged x+1)
#' @description n-year annuity due
#' @export
adueselectxplus1 <- function(x) {
return((1-Aselectxplus1(x))/d(i))
}
#' @title Term life annuity due (select life aged x+1)
#' @description n-year annuity due
#' @export
adueselectxplus1term <- function(x,n) {
return((1-Aselectxplus1endowment(x,n))/d(i))
}
#' @title Term life annuity due (select life aged x+1)
#' @description n-year annuity due
#' @export
adueselectxplus1termupperm <- function(x,n,m) {
return((1-Aselectxplus1endowmentupperm(x,n,m))/dupperm(i,m))
}
#' @title Whole life guaranteed annuity due (select life aged x)
#' @description n-year annuity due
#' @export
adueselectxguaranteed <- function(x,n) {
if (n == 1) {
return(adueanglen(i, 1) + nEselectx(1, x) * adueselectxplus1(x))
}
else {
return(adueanglen(i, n) + nEselectx(n, x) * aduex(x + n))
}
}
#' @title Whole life guaranteed annuity due mthly (select life aged x)
#' @description n-year annuity due
#' @export
adueselectxguaranteedupperm <- function(x, n, m) {
if (n == 1) {
return(adueanglenupperm(i, 1, m) + nEselectx(n, x) * adueselectxplus1upperm(x, m))
} else {
return(adueanglenupperm(i, n, m) + nEselectx(n, x) * aduexupperm(x + n, m))
}
}
#' @title Whole life annuity due mthly (select life aged x)
#' @description n-year annuity due
#' @export
adueselectxupperm <- function(x, m) {
return((1 - Aselectxupperm(x, m)) / dupperm(i, m))
}
#' @title Whole life annuity due mthly (select life aged x+1)
#' @description n-year annuity due
#' @export
adueselectxplus1upperm <- function(x, m) {
return ((1 - Aselectxupperm(x, m)) / dupperm(i, m))
}
#' @title Term life annuity due (select life aged x)
#' @description n-year annuity due
#' @export
adueselectxterm <- function(x, n) {
return ((1 - Aselectxendowment(x, n)) / d(i))
}
#' @title Term life annuity due mthly (select life aged x)
#' @description n-year annuity due
#' @export
adueselectxtermupperm <- function(x, n, m) {
return((1 - Aselectxendowmentupperm(x, n, m)) / dupperm(i, m))
}
#' @title Whole life annuity due
#' @description n-year annuity due
#' @export
aduex <- function(x) {
return((1 - Ax(x)) / d(i))
}
#' @title Term life annuity due
#' @description n-year annuity due
#' @export
aduexterm <- function(x, n) {
return((1 - Axendowment(x, n)) / d(i))
}
#' @title Term life annuity due mthly
#' @description n-year annuity due
#' @export
aduextermupperm <- function(x, n, m) {
return((1 - Axendowmentupperm(x, n, m)) / dupperm(i, m))
}
#' @title Endowment insurance
#' @description n-year annuity due
#' @export
Axendowment <- function(x, n, Moment= 1) {
return(Axterm(x, n, Moment) + nEx(n, x, Moment))
}
#' @title Endowment insurance mthly
#' @description n-year annuity due
#' @export
Axendowmentupperm <- function(x, n, m, Moment= 1) {
return(Axtermupperm(x, n, m, Moment) + nEx(n, x, Moment))
}
#' @title Term insurance
#' @description n-year annuity due
#' @export
Axterm<-function(x, n, Moment= 1) {
#whole life insurance minus a deferred whole life insurance
return(Ax(x, Moment) - ndeferredAx(n, x, Moment))
}
#' @title Term insurance mthly
#' @description n-year annuity due
#' @export
Axtermupperm <- function(x, n, m, Moment= 1) {
# whole life insurance minus a deferred whole life insurance
return(Axupperm(x, m, Moment) - ndeferredAxupperm(n, x, m, Moment))
}
#' @title Whole life guaranteed annuity
#' @description n-year annuity due
#' @export
aduexguaranteed <- function(x, n){
return(adueanglen(i, n) + nEx(n, x) * aduex(x + n))
}
#' @title Whole life guaranteed annuity mthly
#' @description n-year annuity due
#' @export
aduexguaranteedupperm <- function(x, n, m) {
return(adueanglenupperm(i, n, m) + nEx(n, x) * aduexupperm(x + n, m))
}
#' @title Whole life annuity mthly
#' @description n-year annuity due
#' @export
aduexupperm <- function(x, m) {
return((1 - Axupperm(x, m)) / dupperm(i, m))
}
#' @title Whole life insurance (select life aged x)
#' @description n-year annuity due
#' @export
Aselectx <- function(x, Moment= 1) {
return(v(i, 1, Moment) * tqselectx(1, x) + v(i, 2, Moment) * tqselectxplus1(1, x) * tpselectx(1, x) + v(i, 2, Moment) * tpselectx(2, x) * Ax(x + 2, Moment))
}
#' @title Whole life insurance mthly (select life aged x)
#' @description n-year annuity due
#' @export
Aselectxupperm <- function(x, m,Moment= 1) {
if (Moment == 2) {
istar = exp(2 * log(1 + i)) - 1
istarupperm = iupperm(istar, m)
return(Aselectx(x, Moment) * istar / istarupperm)
} else {
return(Aselectx(x) * i / iupperm(i, m))
}
}
#' @title Whole life insurance (select life aged x+1)
#' @description n-year annuity due
#' @export
Aselectxplus1 <- function(x, Moment= 1) {
return(v(i, 1, Moment) * tqselectxplus1(1, x) + v(i, 1, Moment) * tpselectxplus1(1, x) * Ax(x + 2, Moment))
}
#' @title Whole life insurance mthly (select life aged x+1)
#' @description n-year annuity due
#' @export
Aselectxplus1upperm <- function(x, m , Moment= 1) {
if (Moment == 2) {
istar = exp(2 * log(1 + i)) - 1
istarupperm = iupperm(istar, m)
return(Aselectxplus1(x, Moment) * istar / istarupperm)
} else {
return(Aselectxplus1(x) * i / iupperm(i, m))
}
}
#' @title Endowment insurance (select life aged x)
#' @description n-year annuity due
#' @export
Aselectxendowment <- function(x, n, Moment= 1) {
return(selectxterm(x, n, Moment) + nEselectx(n, x, Moment))
}
#' @title Endowment insurance mthly (select life aged x)
#' @description n-year annuity due
#' @export
Aselectxendowmentupperm <- function(x, n, m , Moment = 1) {
return(Aselectxtermupperm(x, n, m, Moment) + nEselectx(n, x, Moment))
}
#' @title Term insurance (select life aged x)
#' @description n-year annuity due
#' @export
Aselectxterm <- function(x, n, Moment = 1) {
return(Aselectx(x, Moment) - ndeferredAselectx(n, x, Moment))
}
#' @title Term insurance (select life aged x+1)
#' @description n-year annuity due
#' @export
Aselectxplus1term <- function(x, n, Moment= 1) {
return(Aselectxplus1(x, Moment) - ndeferredAselectxplus1(n, x, Moment))
}
#' @title Term insurance mthly (select life aged x+1)
#' @description n-year annuity due
#' @export
Aselectxplus1termupperm <- function(x, n, m, Moment = 1) {
return(Aselectxplus1upperm(x, m, Moment) - ndeferredAselectxplus1upperm(n, x, m, Moment))
}
#' @title Deferred whole life insurance (select life aged x+1)
#' @description n-year annuity due
#' @export
ndeferredAselectxplus1 <- function(n, x, Moment= 1) {
if (n == 0) {
return(nEselectx(n, x, Moment) * Aselectxplus1(x, Moment))
} else {
return(nEselectx(n, x, Moment) * Ax(x + n, Moment))
}
}
#' @title Deferred whole life insurance mthly (select life aged x+1)
#' @description n-year annuity due
#' @export
ndeferredAselectxplus1upperm <- function(n, x, m, Moment= 1) {
if (n == 0) {
return(nEselectx(n, x, Moment) * Aselectxplus1upperm(x, Moment))
} else {
return(nEselectx(n, x, Moment) * Axupperm(x + n, Moment))
}
}
#' @title Endowment insurance (select life aged x+1)
#' @description n-year annuity due
#' @export
Aselectxplus1endowment <- function(x, n,Moment= 1) {
return(Aselectxplus1term(x, n, Moment) + nEselectxplus1(n, x, Moment))
}
#' @title Endowment insurance mthly (select life aged x+1)
#' @description n-year annuity due
#' @export
Aselectxplus1endowmentupperm <- function(x, n, m, Moment = 1) {
return(Aselectxplus1termupperm(x, n, m, Moment) + nEselectxplus1(n, x, Moment))
}
#' @title Term insurance mtly (select life aged x)
#' @description n-year annuity due
#' @export
Aselectxtermupperm <- function(x, n, m, Moment = 1) {
return(Aselectxupperm(x, m, Moment) - ndeferredAselectxupperm(n, x, m, Moment))
}
#' @title Whole life insurance
#' @description n-year annuity due
#' @export
Ax <- function(x, Moment = 1) {
Total = 0
for (k in 0:(w - x - 1)) {
Total = Total + udeferredtqx(k, 1, x) * v(i, k + 1, Moment)
}
return(Total)
}
#' @title Whole life insurance mthly
#' @description n-year annuity due
#' @export
Axupperm <- function(x, m, Moment= 1) {
if (Moment == 2){
istar = exp(2 * log(1 + i)) - 1
istarupperm = iupperm(istar, m)
return(Ax(x, Moment) * istar / istarupperm)
} else {
return(Ax(x, Moment) * i / iupperm(i, m))
}
}
#' @title discount rate of interest
#' @description n-year annuity due
#' @export
d<-function(i) {
return(i / (1 + i))
}
#' @title nominal discount rate of interest mthly
#' @description n-year annuity due
#' @export
dupperm <- function(i, m) {
return( m * (1 - (1 - d(i)) ^ (1 / m)))
}
#' @title nominal interest rate mthly
#' @description n-year annuity due
#' @export
iupperm <- function(i, m) {
return( m * ((1 + i) ^ (1 / m) - 1))
}
#' @title Deferred Whole life annuity (select life aged x)
#' @description n-year annuity due
#' @export
ndeferredadueselectx <- function(n, x) {
if (n == 0) {
return(nEselectx(n, x) * adueselectx(x))
} else if(n==1){
return(nEselectx(1, x) * adueselectxplus1(x))
} else {
return(nEselectx(n, x) * aduex(x + n))
}
}
#' @title Deferred Whole life annuity mthly (select life aged x)
#' @description n-year annuity due
#' @export
ndeferredadueselectxupperm <- function(n, x, m ) {
if (n== 0) {
return(nEselectx(n, x) * adueselectxupperm(x, m))
} else if(n==1) {
return(nEselectx(1, x) * adueselectxplus1upperm(x, m))
} else {
return(nEselectx(n, x) * aduexupperm(x + n, m))
}
}
#' @title Deferred Whole life annuity
#' @description n-year annuity due
#' @export
ndeferredaduex <- function(n, x) {
return(nEx(n, x) * aduex(x + n))
}
#' @title Deferred Whole life annuity mthly
#' @description n-year annuity due
#' @export
ndeferredaduexupperm <- function(n, x , m) {
return(nEx(n, x) * aduexupperm(x + n, m))
}
#' @title Deferred Whole life insurance (select life aged x)
#' @description n-year annuity due
#' @export
ndeferredAselectx <- function(n, x, Moment= 1) {
if (n == 0) {
return(nEselectx(n, x, Moment) * Aselectx(x, Moment))
} else if (n == 1) {
return(nEselectx(1, x, Moment) * Aselectxplus1(x, Moment))
} else {
return(nEselectx(n, x, Moment) * Ax(x + n, Moment))
}
}
#' @title Deferred whole life insurance mthly (select life aged x)
#' @description n-year annuity due
#' @export
ndeferredAselectxupperm <- function(n, x, m, Moment= 1) {
if (n == 0) {
return(nEselectx(n, x, Moment) * Aselectxupperm(x, m, Moment))
} else if (n == 1) {
return(nEselectx(1, x, Moment) * Aselectxplus1upperm(x, m, Moment))
} else {
return(nEselectx(n, x, Moment) * Axupperm(x + n, m, Moment))
}
}
#' @title Deferred Whole life insurance
#' @description n-year annuity due
#' @export
ndeferredAx <- function(n, x, Moment= 1) {
return(nEx(n, x, Moment) * Ax(x + n, Moment))
}
#' @title Deferred Whole life insurance (mthly)
#' @description n-year annuity due
#' @export
ndeferredAxupperm <- function(n, x, m, Moment = 1) {
return(nEx(n, x, Moment) * Axupperm(x + n, m, Moment))
}
#' @title Actuarial Discount factor (select life aged x+1)
#' @description n-year annuity due
#' @export
nEselectxplus1 <- function(n, x, Moment= 1) {
return(v(i, n, Moment) * tpselectxplus1(n, x))
}
#' @title Actuarial Discount factor (select life aged x)
#' @description n-year annuity due
#' @export
nEselectx <- function(n, x, Moment= 1) {
return(v(i, n, Moment) * tpselectx(n, x))
}
#' @title Actuarial Discount factor
#' @description n-year annuity due
#' @export
nEx <- function(n, x,Moment= 1) {
return(v(i, n, Moment) * tpx(n, x))
}
#' @title Survival Function (select life aged x)
#' @description Based on the Standard and Ultimate Survival Model (Makeham's Law with )
#' @export
tpselectx <- function(t, age=x) {
# commented this out so that Aselectx plots will work, but this may cause problems
# 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(tpselectx(2,age)*tpx(t-2,age+2))
}
}
#' @title Survival Function (select life aged x+1)
#' @description n-year annuity due
#' @export
tpselectxplus1 <- function(t, age=x) {
# commented this out so that Aselectx plots will work, but this may cause problems
# if ((age + t) >= (w - 1)) {
# 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(tpselectxplus1(1,age)*tpx(t-1,age+1))
}
}
#' @title Survival Function
#' @description n-year annuity due
#' @export
tpx <- function(t, age=x) {
ifelse(((age + t) >= w), 0,exp(-((A * t) + (B / log(c)) * (c ^ age) * (c ^ t - 1))))
}
#' @title Cumulative Distribution Function (select life aged x)
#' @description n-year annuity due
#' @export
tqselectx <- function(t, age=x) {
return(1 - tpselectx(t, age))
}
#' @title Cumulative Distribution Function (select life aged x+1)
#' @description n-year annuity due
#' @export
tqselectxplus1 <- function(t, age=x) {
return(1 - tpselectxplus1(t, age))
}
#' @title Cumulative Distribution Function
#' @description n-year annuity due
#' @export
tqx <- function(t, age=x) {
return(1 - tpx(t, age))
}
#' @title Deferred Whole life guaranteed annuity due (select life aged x)
#' @description n-year annuity due
#' @export
udeferredadueselectxguaranteed <- function(u, x, n) {
if ((u + n) == 0) {
return(adueanglen(i, 1) + nEselectx(1, x) * adueselectx(x))
} else if (u+n==1) {
return(adueanglen(i, 1) + nEselectx(1, x) * adueselectxplus1(x))
} else {
return(v(i, u) * tpselectx(u, x) * adueanglen(i, n) + nEselectx(u + n, x) * aduex(x + u + n))
}
}
#' @title Deferred Whole life guaranteed annuity due mthly (select life aged x)
#' @description n-year annuity due
#' @export
udeferredadueselectxguaranteedupperm <- function(u, x, n, m) {
if ((u + n) == 0) {
return(adueanglenupperm(i, 1, m) + nEselectx(n, x) * adueselectxupperm(x, m))
} else if ((u + n) == 1){
return(adueanglenupperm(i, 1, m) + nEselectx(n, x) * adueselectxplus1upperm(x, m))
} else {
return(v(i, u) * tpselectx(u, x) * adueanglenupperm(i, n, m) + nEselectx(u + n, x) * aduexupperm(x + u + n, m))
}
}
#' @title Deferred whole life guaranteed annuity due (select life aged x+1)
#' @description n-year annuity due
#' @export
udeferredadueselectxplus1guaranteed <- function(u, x, n) {
return(v(i, u) * tpselectxplus1(u, x) * adueanglen(i, n) + nEselectxplus1(u + n, x) * aduex(x + u + n))
}
#' @title Deferred whole life guaranteed annuity due mthly (select life aged x+1)
#' @description n-year annuity due
#' @export
udeferredadueselectxplus1guaranteedupperm <- function(u, x , n, m) {
return(v(i, u) * tpselectxplus1(u, x) * adueanglenupperm(i, n, m) + nEselectxplus1(u + n, x) * aduexupperm(x + u + n, m))
}
#' @title Deferred whole life guaranteed annuity due
#' @description n-year annuity due
#' @export
udeferredaduexguaranteed <- function (u, x, n) {
return(v(i, u) * adueanglen(i, n) * tpx(u, x) + nEx(u + n, x) * aduex(x + u + n))
}
#' @title Deferred whole life guaranteed annuity due mthly
#' @description n-year annuity due
#' @export
udeferredaduexguaranteedupperm <- function(u, x, n, m) {
return(v(i, u) * tpx(u, x) * adueanglenupperm(i, n, m) + nEx(u + n, x) * aduexupperm(x + u + n, m))
}
#' @title Deferred Cumulative Distribution Function
#' @description n-year annuity due
#' @export
udeferredtqx <- function(u, t, x) {
return(tpx(u, x) - tpx(u + t, x))
}
#' @title Deferred Cumulative Distribution Function (select life aged x)
#' @description n-year annuity due
#' @export
udeferredtqselectx <- function(u, t, x) {
return(tpselectx(u, x) - tpselectx(u + t, x))
}
#' @title Present Value Factor
#' @description The present value (discount) factor for a cash flow
#' @param i the interest rate
#' @param n the number of periods to discount
#' @export
v <- function(i, n, Moment= 1) {
return((1 + i) ^ -(n * Moment))
}
#### Life Table Functions ####
#' @title Create Life Table
#' @description creates a table like the one shown in appendix D of DHW (2nd)
#' @param x The age to start the 100000 individuals at
#' @param w The limiting age (omega)
#' @references DHW Actuarial Mathematics for Life Contingent Risks (2nd)
#' @examples TableD.1(20,131)
#' @export
TableD.1 <- function(x, w) {
age <- numeric(0)
lselectx <- numeric(0)
lselectxplus1 <- numeric(0)
lxplus2 <- numeric(0)
ageplus2 <- numeric(0)
for (j in (x-2) : (w-2) ) {
lxplus2 = c(lxplus2, l20*tpx( j-(x-2), x) )
ageplus2 = c(ageplus2, j+2)
if(j<x) {
age=c(age, 0)
lselectx=c(lselectx, 0)
lselectxplus1 = c(lselectxplus1, 0)
} else {
age = c(age,j)
lselectx = c(lselectx, lxplus2[j-(x-3)]/tpselectx(2,j) )
lselectxplus1 = c(lselectxplus1, lxplus2[j-(x-3)]/tpselectxplus1(1,j))
}
}
table <- data.frame(age, lselectx, lselectxplus1, lxplus2, ageplus2)
names(table) = c("x", "l[x]", "l[x]+1", "lx+2", "x+2")
return(table[0:(101-x),])
}
#' @title Create Select Insurance Table
#' @description creates a table like the one shown in appendix D of DHW (2nd)
#' @param x The age to start the 100000 individuals at
#' @param w The limiting age (omega)
#' @references DHW Actuarial Mathematics for Life Contingent Risks (2nd)
#' @examples TableD.2(20,131)
#' @export
TableD.2 <- function(x,w) {
age <- numeric(0)
aselectxcol = Aselectxcolumn(x,w,i,1)
aselectx2col = Aselectxcolumn(x,w,i,2)
adueselectxcol = adueselectxcolumn(x,w,i,1)
adueselectxplus1col = adueselectxplus1column(x,w,i,1)
Eselectx5 <- numeric(0)
Eselectx10 <- numeric(0)
Eselectx20 <- numeric(0)
for (j in (x:w)) {
age = c(age,j)
Eselectx5 = c(Eselectx5, nEselectx(5,j))
Eselectx10 = c(Eselectx10, nEselectx(10,j))
Eselectx20 = c(Eselectx20, nEselectx(20,j))
}
table <- data.frame(age, adueselectxcol, adueselectxplus1col, aselectxcol, aselectx2col, Eselectx5, Eselectx10, Eselectx20,age)
#table = formatDataFrame(table, 5)
names(table) = c("x", '\"a[x]', '\"a[x]+1', "A[x]", "2A[x]", "5E[x]", "10E[x]", "20E[x]","x")
return(table[0:(81-x),])
}
#' @title Create Insurance Table
#' @description creates a table like the one shown in appendix D of DHW (2nd)
#' @param x The age to start the 100000 individuals at
#' @param w The limiting age (omega)
#' @references DHW Actuarial Mathematics for Life Contingent Risks (2nd)
#' @examples TableD.3(20,131)
#' @export
TableD.3 <- function(x,w) {
age <- numeric(0)
aduexcol = aduexcolumn(x,w,i,1)
axcol = Axcolumn(x,w,i,1)
ax2col = Axcolumn(x,w,i,2)
Ex5 <- numeric(0)
Ex10 <- numeric(0)
Ex20 <- numeric(0)
for (j in (x:w)) {
age = c(age,j)
Ex5 = c(Ex5, nEx(5,j))
Ex10 = c(Ex10, nEx(10,j))
Ex20 = c(Ex20, nEx(20,j))
}
table <- data.frame(age, aduexcol, axcol, ax2col, Ex5, Ex10, Ex20,age)
#table = formatDataFrame(table,5)
names(table) = c("x", '\"ax' ,"Ax", "2Ax", "5Ex", "10Ex", "20Ex","x")
return(table[0:(81-x),])
}
#### Helper Methods - Don't Export ####
Axcolumn <- function(x,w,i,moment=1) {
axcol <- numeric(0); temp=0
for (k in 0:(w-x)) {
if(k==0) temp = v(i,1,moment)
else temp = axcol[k]*tpx(1,(w-k))*v(i,1,moment) + v(i,1,moment)*tqx(1,(w-k))
axcol = c(axcol, temp)
}
return(rev(axcol))
}
Aselectxcolumn <- function(x,w,i,moment=1) {
aselectxcol <- numeric(0); temp = 0
axplus2col = Axcolumn(x+2, w+2, i, moment)
for (k in 0:(w-x)) {
temp = (v(i,1,moment)*tqselectx(1,x+k) + v(i, 2, moment)*tqselectxplus1(1,x+k)*tpselectx(1,x+k)
+ v(i,2,moment)*tpselectx(2,x+k)*axplus2col[k+1])
aselectxcol = c(aselectxcol, temp)
}
return(aselectxcol)
}
Aselectxplus1column <- function(age=x,w,i,moment=1) {
aselectxplus1col <- numeric(0); temp = 0
for(k in age:w) {
aselectxplus1col = c(aselectxplus1col, Aselectxplus1(k))
}
return(aselectxplus1col)
}
aduexcolumn <- function(x,w,i,moment=1) {
aduexcol <- numeric(0); temp=0
axcol = Axcolumn(x,w,i,moment)
for (k in 0:(w-x)) {
aduexcol = c(aduexcol, (1-axcol[k])/d(i))
}
aduexcol = c(aduexcol,0)
return(aduexcol)
}
adueselectxcolumn <- function(x,w,i,moment=1) {
aduexcol <- numeric(0); temp=0
aselectxcol = Aselectxcolumn(x,w,i,moment)
for (k in 0:(w-x)) {
aduexcol = c(aduexcol, (1-aselectxcol[k])/d(i))
}
aduexcol = c(aduexcol,0)
return(aduexcol)
}
adueselectxplus1column <- function(age=x, w, i, moment=1) {
aduexcol <- numeric(0); temp=0
aselectxplus1col = Aselectxplus1column(age,w,i,moment)
for (k in 0:(w-age)) {
aduexcol = c(aduexcol, (1-aselectxplus1col[k])/d(i))
}
aduexcol = c(aduexcol,0)
return(aduexcol)
}
#### Data to Export ####
#' @title TableD.1 (DHW 2nd)
#' @description A life table with l[x], l[x]+1 and lx+2 columns for x=20 and w=131
#' @details A different table can be created by using the TableD.1(x,w) function.
#' @export
table1 = TableD.1(20,131)
#' @title TableD.2 (DHW 2nd)
#' @description A life table with adue[x], adue[x]+1, A[x], 2A[x], 5E[x], 10E[x], 20E[x] columns.
#' @details A different table can be created by using the TableD.2(x,w) function. Note that 2A[x] corresponds to A[x] at twice the force of interest (2*delta)
#' @export
table2 = TableD.2(20,131)
#' @title TableD.3 (DHW 2nd)
#' @description A life table with aduex, Ax, 2Ax, 5Ex, 10Ex, 20Ex columns.
#' @details A different table can be created by using the TableD.3(x,w) function. Note that 2Ax corresponds to Ax at twice the force of interest (2*delta)
#' @export
table3 = TableD.3(20,131)
#### Simulation - most of these methods shouldn't be exported ####
# used by simulationTable
kdeferredqx <- function(k,x) {
return(tpx(k,x) - tpx(k+1,x))
}
# used by simulationTable
kdeferredselectqx <- function(k,x) {
return(tpselectx(k,x) - tpselectx(k+1,x))
}
# 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])
}
#' @title Create Simulation Table
#' @description k, v^k, P(Z=v^k, k|qx, using select probabilities
#' @param x the age of the individuals (select age) for the simulation
#' @param w the limiting age (120)
#' @details slight error on the last few cdf values due to rounding to 5 decimal places
#' @examples simulationTable(50,131)
#' @export
simulationTable <- function(x,w) {
# method that formats at truncates vectors in the table
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))
}
# use select probabilities
pdfZ <- function(k,x) {
if(0<k & k<(w-1)) return(kdeferredselectqx(k-1,x))
return(0)
}
kdeferredqx_ <- numeric(0)
kx <- numeric(0)
v_k <- numeric(0)
zv_k <- numeric(0)
truncateAt = w-x
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]"
return(table)
}
#' @title Run Simulation
#' @description run simulation used in Assignment 4 of Acma 320
#' @param n the sample size (default 50)
#' @param age the age (x) for the simulation (default 50)
#' @examples runSimulation(); runSimulation(n=500)
#' @export
runSimulation <- function(n=50, age=50) {
theoreticalZ = simtable[,5]
cdfK = simtable[,3]
sdfZ = simtable[,6]
MeanZ = Aselectx(age)
SdZ = sqrt(Aselectx(age,2)-Aselectx(age,1)^2)
# P(Z>=A[age])
prob = sdf(MeanZ,sdfZ)
cat("\n"); cat("P(Z>= A[age]): ", prob); cat("\n")
# cdfK is a variable declared outside the function
# w is also accessible since it was declared outside the function
simKx = kNums(n,age,w,cdfK)
# default interest rate is 0.05
simZ = zNums(simKx, 0.05)
hist(simZ, main=paste("Histogram of Z (n = ", n, ")"), xlab ="Z", col="pink")
simMeanZ = mean(simZ)
simSdZ = sd(simZ)
cat("n is: ", n)
cat("\nMeanZ: ", MeanZ)
cat("\nsimMeanZ: ", simMeanZ)
cat("\nSdZ: ", SdZ)
cat("\nsimSdZ: ", simSdZ)
simCdfZ = ecdf(simZ)
cat("\nP(Z>= A[age]: ~", (1-simCdfZ(MeanZ)))
cat("\n")
}
#' @title Simulation Data Table
#' @description Provides the data used in the example simulation
#' @details This table was created using x = 50, and generating random K(x) values based on the underlying survival model. Other tables can also be created using the same approach. Note that K(x) values were generated using random uniform distribution transformation
#' @export
simtable = simulationTable(50,131)
#' @title Simulation Z plot
#' @description Plot simulated z values
#' @examples simZplot()
#' @export
simZplot <- function() {
plot(simtable[,5], main="Theoretical Distribution of Z", ylab="P(Z)=z)", xlab="Z", col="blue")
}
#### Curtate, Complete Expectations, Variance Formulas ####
Rx <- function(x) {
if (x >= 0 & x <= 39) {
return(0.999)
} else if (x >= 40 & x <= 59) {
return(0.998)
} else if (x >= 60 & x <= 79) {
return(0.995)
} else if (x >= 80 & x <= w) {
return(0.99)
} else {
return(0)
}
}
# use survival function
ncurtate <- function(n) {
total = 0
for (k in 0:(n - 1)) {
total = total + k * kdeferredqx(k,x)
}
total = total + n * tpx(n)
return(total)
}
# get answer from integrating the survival function
ncomplete <- function(n) {
return(midpoint(tpx, 0, n, 1000))
}
# use survival function; solve by letting tpx = 0.5
# or use life table and find where radix/2 has been surpassed
median <- function() {
i = 1
# account for fact that they all may die at the end
if (lxcol[w-x] > (radix/2)) {
return(w - 1)
} else {
while (lxcol[i] > radix/2) {
i = i + 1
}
return(i + x - 1)
}
}
# calculate this based on a life table
mode <- function() {
# mode of age at death is when the most people die calculate the mode based on some age x
max = 0
age = 0
# process each entry in dx
for (i in 1:(w - x)) {
if (dxcol[i] > max) {
max = dxcol[i]
age = i
}
}
# age corresponds to additional age, so add this to x and subtract 1 for the floor
return(age + x - 1)
}
# get answer from integrating the survival function
stdevTx <- function() {
Etx = midpoint(tpx, 0, (w - x), 1000)
intFunction <- function(t) {
return(t*tpx(t))
}
Etx_squared = 2 * midpoint(intFunction, 0, (w - x), 1000)
return(sqrt((Etx_squared - Etx^2)))
}
stdevKx <- function() {
Ekx <- function() {
total = 0
for (k in 2:(w - x)) {
total = total + lxcol[k]
}
total = total/lxcol[1]
return(total)
}
Ekx_squared <- function() {
total = 0
for (k in 2:(w - x)) {
total = total + lxcol[k] * (2 * (k - 1) - 1)
}
total = total/lxcol[1]
return(total)
}
return(sqrt(Ekx_squared() - Ekx()^2))
}
# adjust the life.table provided by s years
adjust <- function(lifetable, s) {
lx <- numeric(0)
dx <- numeric(0)
qx <- numeric(0)
px <- numeric(0)
agecol = lifetable[,1]
for (i in x:(w - 1)) {
if (i == x) {
lx = c(lx, radix * tpx(i - x, x))
# adjust qx column by appropriate amount
qx = c(qx, (lifetable["qx"][((i - x) + 1), 1] * Rx(x)^s))
dx = c(dx, (qx[(i - x) + 1] * lx[(i - x) + 1]))
px = c(px, (1 - qx[(i - x) + 1]))
} else if (x < i & i < (w - 2)) {
lx = c(lx, (lx[(i - x)] - dx[(i - x)]))
qx = c(qx, (lifetable["qx"][((i - x) + 1), 1] * Rx(i)^s))
dx = c(dx, (qx[(i - x) + 1] * lx[(i - x) + 1]))
px = c(px, (1 - qx[(i - x) + 1]))
# doesnt apply the Rx adjustment to the (lw-1) row to ensure that all the
# survivors die by the limiting age
} else {
lx = c(lx, (lx[(i - x)] - dx[(i - x)]))
qx = c(qx, (lifetable["qx"][((i - x) + 1), 1] * Rx(i)^0))
dx = c(dx, (qx[(i - x) + 1] * lx[(i - x) + 1]))
px = c(px, (1 - qx[(i - x) + 1]))
}
}
return(data.frame(agecol,lx, dx, qx, px))
}
# results
# table1 = createLifeTable(0)
# table2 = createLifeTable(x)
# table3 = adjust(table2, 10)
# table4 = adjust(table2, 25)
#
# lxcol = table2[,2]
# dxcol = table2[,3]
# qxcol = table2[,4]
# pxcol = table2[,5]
#
# plot(tpx, 0, (w-x), ylab = "tpx", xlab = "T(x)", main = "Survival Function")
# plot(seq(0:(w-x-1))-1, dxcol, type = "l", main = "Curve of Deaths", xlab = "T(x)", ylab = "Deaths")
#### Assignment 7 ####
# insurance1epv <- function(age=x,n,m=1) {
# return((1-Aselectxterm(age,n,m))/d(i,m))
# }
#
# insurance1var <- function(age=x,n,m=1) {
# return( (Aselectxterm(age,n,m,moment=2) - Aselectxterm(age,n,m)^2)/(d(i,m)^2))
# }
#
# insurance2epv <- function(age=x,n,m=1) {
# return (nEselectx(n,age)*(1-(Ax(age+n)*i/iupperm(i,m)))/d(i,m))
# }
#
# insurance2var <- function(age=x,n,m=1) {
# j=exp(2*log(1+i))-1
# return(tpselectx(n,age)*v(i,n,moment=2)*
# ((((1-Ax(age+n)*i/iupperm(i,m))/d(i,m))^2*(1-tpselectx(n,age)))
# + (( Ax(age+n,moment=2)*j/iupperm(j,m) - (Ax(age+n)*i/iupperm(i,m))^2)/d(i,m)^2)))
# }
#
# getPercentile <- function(EPV,SD,percentile,size=numberContracts) {
# zvalue = qnorm(percentile)
# return (zvalue*sqrt(size)*SD + size*EPV)
# }
#
# EPV = c(insurance1epv(x,n,m),insurance2epv(x,n,m))
# SD = c(sqrt(insurance1var(x,n,m)), sqrt(insurance2var(x,n,m)))
#
# numberContracts = 100
# p70 = c(getPercentile(EPV[1],SD[1],0.70),getPercentile(EPV[2],SD[2],0.70))
# p95 = c(getPercentile(EPV[1],SD[1],0.95),getPercentile(EPV[2],SD[2],0.95))
# results <- data.frame(EPV, SD, p70,p95)
# cat(paste("NUMBER OF CONTRACTS = ", numberContracts))