## 14.02.18 ## 12.10-12.15
#### 12.10
# a)--------------------------------------------------------------------------
#compute estimates of Japanese men (first column) and Us men (second column) mortalities
#using ql = (l - nl+1)/nl, nl is survivors up to age l years
nn = matrix(scan("japusmale.txt"),byrow=T,ncol=2)
l_e = nrow(nn) -1
ql = 1-nn[1:l_e +1,]/nn[1:l_e,]
list(q=q,nn=nn,l_e=l_e)
# Japanese male mortalities in 2008 in column 1 and US males (2007) in
# column 2 read from file.
#
# On output estimated mortalities are stored as the columns of "q" (matrix)
# and the number of individuals under risk as the columns of "nn" (matrix).
# "l_e" is the maximum age.
#b)
#Plot the estimated mortalities jointly against age
par(mfrow=c(1,1))
ll = 1:l_e -1
#Japanese
plot(ll,ql[,1], xlab="age", ylab="q", main="mortality risk")
#US
lines(ll,ql[,2],lwd=2)
legend("topleft",c("US males 2007","Japanese males 2008"),lty=c(1,NA),pch = c(NA,1))
#text(40,0.4,"Points: Japanese males 2008")
#text(40,0.38,"Solid: US males 2007")
# from the plot it looks like the US males lived slightly shorter than the Japanese males
# in 2008. This is because the line is slightly higer than the plotted points, which means
# they had a higher mortality rate.
#c)
#plot of the logatithms
#Japanese
plot(ll,log(ql[,1]), xlab="age", ylab="log(q)", main="log mortality risk")
#US
lines(ll,log(ql[,2]))
legend("topleft",c("US males 2007","Japanese males 2008"),lty=c(1,NA),pch = c(NA,1))
# Again looking at the log of the data. The US males seam to live shorter than the Japanese
# There is a clear differensce up tp the age of 60 then the data become more and more equal
#--------------------------------------------------------------------------
#### 12.11
#redo the exercise, but for female
#dir()
nn = matrix(scan("japusfem.txt"),byrow=T,ncol=2)
l_e = nrow(nn)-1
ql = 1-nn[1:l_e +1,]/nn[1:l_e,]
list(q=q,nn=nn,l_e=l_e)
# Japanese female mortalities in 2008 in column 1 and US females (2007)
# in column 2 read from file.
#
# On output estimated mortalities are stored as the columns of "q" (matrix)
# and the number of individuals under risk as the columns of "nn" (matrix).
# "l_e" is the maximum age.
#plot
ll = 1:l_e -1
#Japanese
plot(ll,ql[,1], xlab="age", ylab="mortalities", main="mortality risk")
#US
lines(ll,ql[,2])
legend("topleft",c("US females 2007","Japanese females 2008"),lty=c(1,NA),pch = c(NA,1))
#plot of the logatithms
#Japanese
plot(ll,log(ql[,1]), xlab="age", ylab="log mortalities",main="survivalprobabilities per year on Log Scale")
#US
lines(ll,log(ql[,2]))
legend("topleft",c("US females 2007","Japanese females 2008"),lty=c(1,NA),pch = c(NA,1))
# We see the same trend for females as for males, where the US females in 2007 have a higher
# probability of dying at all ages, compared to the Japanese females in 2008.
# Looking at the log of the probability of mortality per year for the females
# we see the same trend as before. The US females in 2007 have a higer log probability
# than the Japanese females in 2008. There is a clearer difference between the two
# groups of females than the two groups of males. Where they reach a higher age before meeting
# the same mortality rate.
#--------------------------------------------------------------------------
#### 12.12
# a)
#using data from 12.10
# Life table computed "K" periods ahead for every initial age.
# The input is the mortality vector "q".
#
# On output "kp" (matrix) contains the life table as columns for
# each initial age up to maximum age "l_e".
ex1212 = function(q,K){
l_e=length(q)
p=c(1-q,rep(0,l_e+1))
kp=matrix(1,K+1,l_e+1)
for(l in 0:l_e+1){
kp[1:K+1,l] = cumprod(p[l:(l+K-1)])
}
list(kp=kp,l_e=l_e)
}
#"kp" (matrix) contains the life table as columns for
# each initial age up to maximum age "l_e"
# Life table computed "K" periods ahead for every initial age.
# The input is the mortality vector "q".
# b)
nn = matrix(scan("japusmale.txt"),byrow=T,ncol=2)
l_e = nrow(nn)-1
ql = 1-nn[1:l_e +1,]/nn[1:l_e,]
q1 = ql[,1]
ex1212(q1,50)
# The life table "K" periods ahead plotted for Japanese males when "q"
# is the non-parameteric estimate from 12.10 a).
#plot when K = 50, l_e=110 against l = 10,30,...,100 jointly for Japanese males
ind = 1+1:10*10
K=50
matplot(0:K,ex1212(q1,50)$kp[,ind],"l",xlab="Time ahead (years)",ylab="Survival probability")
title("Survival Japanese males: Initial age is 10,20,...,100 years")
#which line belongs to which group of the given initial ages?
#c)
# The life table "K" periods ahead plotted for US males when "q"
# is the non-parameteric estimate from 12.10 a).
#plot when K = 50, l_e=110 against l = 10,30,...,100 jointly for US males
q2 = ql[,2]
ind = 1+1:10*10
matplot(0:K,ex1212(q2,50)$kp[,ind],"l",xlab="Time ahead (years)",ylab="Survival probability")
title("Survival US males: Initial age is 10,20,...,100 years")
#--------------------------------------------------------------------------
#### 12.13
# a) use p and kp from exercise 12.12
# The probability of dying after "k" periods computed for "k" up to
# "K" every intial age in the probabiliy vector "q".
#
# On output "kq" (matrix) contains the probabilities as columns for
# each initial age up to the maximum age "l_e".
ex1213=function(p,K,kp,l_e){
kq=matrix(0,K+1,l_e+1)
for(l in 0:l_e+1){
kq[1:K+1,l]=kp[1:K,l]*(1-p[1:K+(l-1)])
}
list(kq=kq,l_e=l_e)
}
p1=c(1-q1,rep(0,ex1212(q1,100)$l_e+1))
o = ex1213(p1,100,ex1212(q1,100)$kp,ex1212(q1,100)$l_e)
o
# b)
K=100
ind=1+(1:10)*min(floor(o$l_e/10),10)
matplot(0:K,o$kq[,ind],"l",xlab="Number of years until dying",
ylab="Probability of dying")
text(0.5*K,0.5*max(o$kq[,ind]),"Intial age 10,20,...,100 years,")
text(0.5*K,0.45*max(o$kq[,ind]), "one for each curve")
title("Mortalities over time for Japanse males 2008")
# c)
p2=c(1-q2,rep(0,ex1212(q2,100)$l_e+1))
o = ex1213(p2,100,ex1212(q2,100)$kp,ex1212(q2,100)$l_e)
o
K=100
ind=1+(1:10)*min(floor(o$l_e/10),10)
matplot(0:K,o$kq[,ind],"l",xlab="Number of years until dying",
ylab="Probability of dying")
text(0.5*K,0.5*max(o$kq[,ind]),"Intial age 10,20,...,100 years,")
text(0.5*K,0.45*max(o$kq[,ind]), "one for each curve")
title("Mortalities over time for US males 2007")
#-------------------------------------------------------------------------
#### 12.14
# a)
nn = matrix(scan("japusmale.txt"),byrow=T,ncol=2)
l_e = nrow(nn)-1
ql = 1-nn[1:l_e +1,]/nn[1:l_e,]
list(q=q,nn=nn,l_e=l_e)
ll=1:l_e-1
length(ll)
z = NA
z = nn[1:110,1]-nn[2:111,1]
#z1 = c(0,z)
length(z)
length(nn[1:110,1])
theta_hat= gompmakfit1(ll,n=nn[1:110,1],z)$theta
qgm=1-exp(-theta_hat[1]-theta_hat[2]*(exp(theta_hat[3])-1)*exp(theta_hat[3]*ll)/theta_hat[3])
# Now the plot:
plot(ll,ql[,1],xlab="Age (year)",ylab="Mortalities")
lines(ll,qgm)
# Text to plot.
legend("topleft",c("Data","Gomperz-Makeham estimation"),pch=c(1,NA),lty=c(NA,1))
title("Japanese males 2008")
list(theta_hat=theta_hat)
#3.268162e-04 1.636826e-05 1.021868e-01
# b) the same, but for US males instead of Japanese males
ll=1:l_e-1
length(ll)
z = NA
z = nn[1:110,2]-nn[2:111,2]
#z1 = c(0,z)
length(z)
length(nn[1:110,2])
theta_hat= gompmakfit1(ll,n=nn[1:110,2],z)$theta
qgm=1-exp(-theta_hat[1]-theta_hat[2]*(exp(theta_hat[3])-1)*exp(theta_hat[3]*ll)/theta_hat[3])
# Now the plot:
plot(ll,ql[,2],xlab="Age (year)",ylab="Mortalities")
lines(ll,qgm)
# Text to plot.
legend("topleft",c("Data","Gomperz-Makeham estimation"),pch=c(1,NA),lty=c(NA,1))
title("US males 2007")
list(theta_hat=theta_hat)
#7.815133e-04 3.945333e-05 9.266762e-02
#c) the same, but for Japanese females.
nn = matrix(scan("japusfem.txt"),byrow=T,ncol=2)
l_e = nrow(nn)-1
ql = 1-nn[1:l_e +1,]/nn[1:l_e,]
list(q=q,nn=nn,l_e=l_e)
ll=1:l_e-1
length(ll)
z = NA
z = nn[1:110,1]-nn[2:111,1]
length(z)
length(nn[1:110,1])
theta_hat= gompmakfit1(ll,n=nn[1:110,1],z)$theta
qgm=1-exp(-theta_hat[1]-theta_hat[2]*(exp(theta_hat[3])-1)*exp(theta_hat[3]*ll)/theta_hat[3])
# Now the plot:
plot(ll,ql[,1],xlab="Age (year)",ylab="Mortalities")
lines(ll,qgm)
# Text to plot.
legend("topleft",c("Data","Gomperz-Makeham estimation"),pch=c(1,NA),lty=c(NA,1))
title("Japanese females 2008")
list(theta_hat=theta_hat)
#3.364727e-04 1.817205e-06 1.213192e-01
#d) the same but for US females
ll=1:l_e-1
length(ll)
z = NA
z = nn[1:110,2]-nn[2:111,2]
length(z)
length(nn[1:110,2])
theta_hat= gompmakfit1(ll,n=nn[1:110,2],z)$theta
qgm=1-exp(-theta_hat[1]-theta_hat[2]*(exp(theta_hat[3])-1)*exp(theta_hat[3]*ll)/theta_hat[3])
# Now the plot:
plot(ll,ql[,2],xlab="Age (year)",ylab="Mortalities")
lines(ll,qgm)
# Text to plot.
legend("topleft",c("Data","Gomperz-Makeham estimation"),pch=c(1,NA),lty=c(NA,1))
title("US females 2007")
list(theta_hat=theta_hat)
#4.674124e-04 2.598666e-05 9.753583e-02
par(mfrow=c(2,2))
##---------------------------------------------------------------------------------------
#### 12.15
# a)
# The probability of dying after "k" periods uo to "K" computed for every
# intial age up to "l_e". The model is Gomperz-Makeham with parameters
# in "theta" (vector).
#
# On output "kq" (matrix) contains the probabilities as columns for
# each initial age up to the maximum "l_e".
ex1215=function(l_e,theta,K) {
ll=0:(2*l_e)
hh=exp(theta[3]*ll)
kp=exp(-theta[1]*ll-(theta[2]/theta[3])*(hh-1)%o%hh)
kq=matrix(0,K+1,l_e+1)
for (l in 0:l_e+1){
kq[1:K+1,l]=kp[1:K,l]*(1-kp[2,l:(l+K-1)])
}
list(kq=kq)
}
# Numerical application of exercise 12.15 a) with parameters those for US
# males 2007.
kq=ex1215(l_e=110,theta=c(0.000778,0.0000374,0.0928),K=100)$kq
# After computing the mortaties, then the plot
ind=1+(1:10)*min(floor(l_e/10),10)
matplot(0:K,kq[,ind],"l",xlab="Number of years until dying",
ylab="Mortalities")
text(0.5*K,0.5*max(kq[,ind]),"Intial age 10,20,...,100 years,")
text(0.5*K,0.45*max(kq[,ind]), "one for each curve")
title("Mortalities over long for US male 2007 Gomperz-Makeham")
## 21.02.18 ## 12.16, 12.18, 12.21-12.24 ##
#### 12.16
#a)
l_1=18;l_2=100;theta=c(0.000778,0.0000374,0.0928)
h=1
E=rep(0,l_2+1)
for (l in l_1:l_2){
E[l+1]=gompmakmean(theta,l,h)
}
plot(l_1:l_2,E[l_1:l_2+1],"l",xlab="Age (years)",ylab="Remaining life (years)")
title("Remaining number of years")
# Computes and plots the expected number of years remaining against age
# under the Gomperz-Makeham model with parameters "theta" (vector). The
# plot applies to age between "l_1" and "l_2". Time increment "h=1" for the
# numerical implementation.
#
# R-function used: gompmakmean.
#b)
l_1=18;l_2=100;theta=c(0.000778,0.0000374,0.0928)
E=rep(0,l_2+1)
# h = 1 (as in a))
h=1
for (l in l_1:l_2){
E[l+1]=gompmakmean(theta,l,h)
}
plot(l_1:l_2,E[l_1:l_2+1],"l",xlab="Age (years)",ylab="Remaining life (years)")
# h = 0.1
h=0.1
for (l in l_1:l_2){
E[l+1]=gompmakmean(theta,l,h)
}
lines(l_1:l_2,E[l_1:l_2+1],"l",col="blue")
# h = 0.01
h=0.01
for (l in l_1:l_2){
E[l+1]=gompmakmean(theta,l,h)
}
lines(l_1:l_2,E[l_1:l_2+1],"l",col="red",lty=2)
legend("topright",c("h=1","h=0.1","h=0.01"),col=c("black","blue","red"),pch="-")
text(l_1+0.7*(l_2-l_1),0.75*max(E),"Different time increments")
text(l_1+0.7*(l_2-l_1),0.70*max(E),"for the numerical integration")
title("Remaining number of years")
#c)
l_1=18;l_2=100;theta=c(0.000778,0.0000374,0.0928)
E=rep(0,l_2+1)
# h = 1
h=1
for (l in l_1:l_2){
E[l+1]=gompmakmean(theta,l,h)+l
}
plot(l_1:l_2,E[l_1:l_2+1],"l",xlab="Age (years)",ylab="Total life (years)")
# h = 0.1
h=0.1
for (l in l_1:l_2){
E[l+1]=gompmakmean(theta,l,h)+l
}
lines(l_1:l_2,E[l_1:l_2+1],"l",col="blue")
# h = 0.01
h=0.01
for (l in l_1:l_2){
E[l+1]=gompmakmean(theta,l,h)+l
}
lines(l_1:l_2,E[l_1:l_2+1],"l",col="red",lty=2)
text(l_1+0.6*(l_2-l_1),0.95*max(E),"Different time increments")
text(l_1+0.54*(l_2-l_1),0.935*max(E),"for the numerical integration")
title("Total number of years")
legend("topleft",c("h=1","h=0.1","h=0.01"),col=c("black","blue","red"),pch="-")
help(text)
#------------------------------------------------------------------------------------
#### 12.18 by hand
#### 12.21 by hand
#### 12.22
#a)
ex1222 = function(kp,d,K_1){
#kp = lifetable
#d = discount
#K_1 = max annuity
M=d**(0:(K_1-1))*kp[1:K_1,]
a=apply(M,2,cumsum)
list(a=a)
#Output is a matrix "a" with coeffiecients for each initial age stored
#as columns.
}
#b)
nn = matrix(scan("japusmale.txt"),byrow=T,ncol=2)
l_e = nrow(nn)-1
ql = 1-nn[1:l_e +1,]/nn[1:l_e,]
q2 = ql[,2]
o = ex1212(q2,50)
K_1=50
#After the life table for USmales, now the annutites:
a=ex1222(o$kp,d=0.97,K_1=50)$a
par(mfrow=c(1,1))
matplot(1:K_1,a[,l_1:l_2+1],"l",xlab="Years ahead",ylab="Coefficients")
title("Mortality adjusted annuity coeffcients for US males")
#c)
#Each age group has a curve from age 25 to 70.
#highest curve is for 25 year olds
#lowest curve is for 70 year olds
#increasing curves in terms of duration
#capital PI -> pension: s = PI/a(lo:K)
#higher pension at an older age -> same sum over fewer years
#curve flattens out when individual dies
#the coefficients is the quantity for payments in advance when alive.
#see p.444 in curriculum book.
#### 12.23
#a) by hand
#b) by hand
#c)
# Computes and prints annuities up to different expiries defined in "K"
# (vector) with mortalities NOT taken into account. "d" is the discount.
d=0.98;K=1:10*5
K
a_k=round((1-d**K)/(1-d),digits=3)
a_K=cbind(K,a_k)
print(a_K)
#for a fixed premium PI the pension recieved at each period gets smaller the more periods we wish
#to sign. This can be seen as the aK values increase the more periods ahead we want.
#d)
d=0.98;l=40;K=1:10*5
K_1=max(K)
# Life table for US males:
o=ex1212(q2,K_1+1)
# Mortality adjusted annuities:
M=d**(0:(K_1-1))*o$kp[1:K_1,l]
a1=matrix(cumsum(M),K_1,1)
a1=round(a1[K,],digits=3)
b=cbind(K,a1)
print(b)
# K a
# [1,] 5 4.782
# [2,] 10 9.040
# [3,] 15 12.800
# [4,] 20 16.082
# [5,] 25 18.905
# [6,] 30 21.277
# [7,] 35 23.199
# [8,] 40 24.670
# [9,] 45 25.692
# [10,] 50 26.294
c = cbind(b,a_k)
print(c)
# K a1 a_k
# [1,] 5 4.782 4.804
# [2,] 10 9.040 9.146
# [3,] 15 12.800 13.072
# [4,] 20 16.082 16.620
# [5,] 25 18.905 19.827
# [6,] 30 21.277 22.726
# [7,] 35 23.199 25.346
# [8,] 40 24.670 27.715
# [9,] 45 25.692 29.856
# [10,] 50 26.294 31.792
# a_k column is the highest, meaning that when we take mortalities into account
# the pension s received at each period, keeping PI fixed will be higher as the
# US males live shorter (a1 < a_k) at age 40 than whats expected when not concidering mortalities and
# only using discount values for K periods ahead.
#e)
d=0.98;l=20;K=1:10*5
K_1=max(K)
# Life table for US males:
#q2
o=ex1212(q2,K_1+1)
# Mortality adjusted annuities:
M=d**(0:(K_1-1))*o$kp[1:K_1,l]
a2=matrix(cumsum(M),K_1,1)
a2=round(a2[K,],digits=3)
x=cbind(c,a2)
print(x)
d=0.98;l=60;K=1:10*5
K_1=max(K)
# Life table for US males:
#q2
o=ex1212(q2,K_1+1)
# Mortality adjusted annuities:
M=d**(0:(K_1-1))*o$kp[1:K_1,l]
a3=matrix(cumsum(M),K_1,1)
a3=round(a3[K,],digits=3)
x2=cbind(x,a3)
print(x2)
colnames(x2) = c("K","l=40","age ind.","l=20","l=60")
print(x2)
compare = cbind(x2[,1],x2[,3],x2[,4],x2[,2],x2[,5])
colnames(compare) = c("K","age ind.","l=20","l=40","l=60")
# K age ind. l=20 l=40 l=60
# [1,] 5 4.804 4.791 4.782 4.699
# [2,] 10 9.146 9.091 9.040 8.647
# [3,] 15 13.072 12.950 12.800 11.847
# [4,] 20 16.620 16.411 16.082 14.294
# [5,] 25 19.827 19.508 18.905 15.995
# [6,] 30 22.726 22.266 21.277 16.998
# [7,] 35 25.346 24.701 23.199 17.450
# [8,] 40 27.715 26.827 24.670 17.581
# [9,] 45 29.856 28.655 25.692 17.603
# [10,]50 31.792 30.191 26.294 17.604
# columns: increasing
# rows: decreasing
# the higher the a the lower the pension for a fixed PI
# make sense that the higher age group get a higher pension (smaller ak) for any period ahead
# as it's a smaller chance that they will be alive to receive the pension-> greater chance that the
# insurance company earn money
#### 12.24
#a)
#re-read table of US males, just in case of an overwrite.
nn = matrix(scan("japusmale.txt"),byrow=T,ncol=2)
l_e = nrow(nn)-1
ql = 1-nn[1:l_e +1,]/nn[1:l_e,]
q2 = ql[,2]
ex1224=function(q2,l_0,s,l_r,d){
# Computes the equivalence premium for a pension scheme entered
# at age "l_0" and paid for in advance each period up to retirement
# age "l_r". The pension is "s" (in advance) and the discount "d".
#
# "pi" is the premium.
#Life table and annuities for US males
l_e=length(q2)
kp=ex1212(q2,l_e+1)$kp
a=ex1222(kp,d,l_e)$a
#Premium
pi=s*d**(l_r-l_0)*kp[l_r-l_0+1,l_0+1]*a[l_e,l_r+1]/a[l_r-l_0,l_0+1]
list(pi=pi)
}
#test program for entry age 30 and retirement age 65, with pension s=1 and discount 0.98
ex1224(q2,l_0=30,s=1,l_r =65,0.98)
# $pi
# [1] 0.2457182
#b)
# Computes and plot equivalence premia against the retirement
# age in "llr" (vector) for a pension scheme. The pension is "s" and the
# discount "d". All payments are in advance. Joint plots are generated
# for individuals entering the scheme at age between "l_1" and "l_2".
l_e=length(q2)
kp=ex1212(q2,l_e+1)$kp
a=ex1222(kp,d,l_e)$a
# After life table and annuities for US males, then premia lowest age.
s=1;d=0.97
llr=55:75;l_1=20;l_2=45
pi=s*d**(llr-l_1)*kp[llr-l_1+1,l_1+1]*a[110,llr+1]/a[llr-l_1,l_1+1]
plot(llr,pi,"l",ylim=c(0,1.5),xlab="Retirement age (years)")
#Premia at higher age
for (l_0 in (l_1+1):l_2) {
pi=s*d**(llr-l_0)*kp[llr-l_0+1,l_0+1]*a[110,llr+1]/a[llr-l_0,l_0+1]
lines(llr,pi)
}
# the below can e used to check bottom and top lines in plot can be used
#bottom line
pi=s*d**(llr-20)*kp[llr-20+1,20+1]*a[110,llr+1]/a[llr-20,20+1]
lines(llr,pi,col="red")
#top line
pi=s*d**(llr-45)*kp[llr-45+1,45+1]*a[110,llr+1]/a[llr-45,45+1]
lines(llr,pi,col="blue")
text(mean(llr),0.75*max(pi),"Initial age of indvidual")
text(mean(llr),0.70*max(pi),"varied from 20 to 45")
title("Equivalence premium against retirement age")
legend("topright",c("lo=20","lo=45"),pch="-",col=c("red","blue"))
# Annual pension s in advance purchased at age lo lasting to end of life le=110
#retirement age lr varies from 55 to 75 years old
#starting age l0 varies from 20 to 45 years old
#funding has to be higher when we start at an early age: longer time to support/get pension.
#lowest curve is l0=20 as we have the longest support/contribution time.
#highest curve is l0=45 as we have the shortest support/contribution time.
## 28.02.18 ## 12.25-12.28 + 12.31-12.32 ##
#### 12.25
# a) make function to compute coefficients for given variables.
ex1225=function(kq,d,K_1){
# kq = table of mortialities (computed from exercise 12.13)
# d = discount
# K = maximum number of periods we want to run for.
# A table of A-coeffients for term insurance is computed given the above.
# Output is matrix of coefficients "A" with one column per age.
#
# Warning: If e1213a is used to set up "kq", it must be called with K> K_1.
M=d**(1:K_1)*kq[1:K_1+1,]
A=apply(M,2,cumsum)
list(A=A)
}
#b)
#data from US males
nn = matrix(scan("japusmale.txt"),byrow=T,ncol=2)
l_e = nrow(nn)-1
ql = 1-nn[1:l_e +1,]/nn[1:l_e,]
q2 = ql[,2]
#use exercise 12.12 to calculate p2 for K=50 periods
p2=c(1-q2,rep(0,ex1212(q2,50)$l_e+1))
#use exrcise 12.13 to compute "kq" (matrix) containing the probabilities as columns for
#each initial age up to the maximum age "l_e".
o = ex1213(p2,50,ex1212(q2,50)$kp,ex1212(q2,50)$l_e)
#use function made in a) to compute coefficients:
A = ex1225(o$kq,d=0.97,K_1=50)$A
#minimum age
l_1 = 25
#maximum age
l_2 = 55
#plot
matplot(1:K_1, A[,l_1:l_2+1],"l",
xlab="Duration of scheme (years)",ylab="A-coeffcient")
title("A-coefficients for US males")
#c interpreter the plot/pattern you see
#### 12.26
#a) make a general program to compute PI using equation 12.28 from the book
ex1226_PI=function(s,a,A,K_1){
#s = sum released
#a = annuity coefficient in contributing period (matrix)
#A = annuity coefficient in "receiving period" (matrix)
#K_1 = maximum of periods we want to run for
# On output "pi" (vector) contains equivalance premia columnnwise for each age for "K"
# between 1 and "K_1".
pi=s*A[1:K_1,]/a[1:K_1,]
list(pi=pi)
}
#b)
# Numerical application of e1226a for mortalities for US males 2007
# with equivalence premia plotted against the duration of the contracts
# up to "K_1" with the age of the policy holder varied between "l_1" and
# "l_2". "s" is released if the policy holder dies. The discount is "d".
kp=ex1212(q2,K_1+1)$kp
kq = ex1213(p2,50,ex1212(q2,50)$kp,ex1212(q2,50)$l_e)$kq
#The annuitiy coeficients
a=ex1222(kp,d,K_1)$a
A=ex1225(kq,d,K_1)$A
# The premia
pi = ex1226_PI(s=1,a=a,A=A,K_1=50)$pi
#Plotting
matplot(1:K_1,pi[,l_1:l_2+1],"l",xlab="Duration of contract (years)",ylab="Premium per period")
text(8,0.022,"Age policy")
text(8,0.0205,"holder varied")
title("Equivalence premium for term insurance")
#### 12.27
#a)
ex1227=function(q,s,l_r,d){
#q = mortality vector, age of infuviduals vary as in q
#s = pension starting at l_r
#l_r = pension age
#d = discount
# Output is a tabulation "pi" of one-time premia against age when the
# contract is drawn up.
#The life table
l_e=length(q)
kp=ex1212(q,l_e)$kp
# After the life table has been computed, the indicator matrix is set up.
I=matrix(0,l_e+1,l_e+1)
I[row(I)+col(I)>l_r+1]=1
# Finally a tabulation of the one-time premium.
ll=0:l_e
M=s*(d**ll)*kp*I
pi=apply(M,2,sum)
list(pi=pi)
}
nn = matrix(scan("japusmale.txt"),byrow=T,ncol=2)
l_e = nrow(nn) -1
ql = 1-nn[1:l_e +1,]/nn[1:l_e,]
q2= ql[,2]
#b)
l_1= 20
l_2=109
PI_1 = ex1227(q=q2,s=1,l_r=60,d=0.97)$pi
plot(l_1:l_2,PI_1[l_1:l_2+1],"l",xlab="Age (years)",ylab="Premium")
title("One-time premium for pension against age")
#c)
PI_2 = ex1227(q=q2,s=1,l_r=60,d=0.96)$pi
lines(l_1:l_2,PI_2[l_1:l_2+1],col="blue",lty=2)
legend("topright",c("d=0.97","d=0.96"),col=c("black","blue"),pch = "-")
#### 12.28
#a) general program to tabulate PI_lo
ex1228=function(l_e,theta,s,l_r,d){
#l_e = maximum age of individuals
#theta = a vector on length 3 with the Gomperz-Makeham parameters
#s = pension s starting at retirement age
#l_r = retirement age
#d = discount
#
# Output is a tabulation "pi" of all one-time premias for individuals
# of age 0 to l_e.
ll=0:l_e
hh=exp(theta[3]*ll)
kp=exp(-theta[1]*ll-(theta[2]/theta[3])*(hh-1)%o%hh)
I=matrix(0,l_e+1,l_e+1)
I[row(I)+col(I)>l_r+1]=1
ll=0:l_e
M=s*(d**ll)*kp*I
pi=apply(M,2,sum)
list(pi=pi)
}
#b) run the program for certain values of the given parameters
l_1=20;l_2=109
pi=ex1228(l_2+1,theta=c(0.000780,0.0000374,0.0928),s=1,l_r=60,d=0.97)$pi
# plot of G-M one-time premia estimates.
plot(l_1:l_2,pi[l_1:l_2+1],"l",xlab="Age (years)",ylab="Premium")
title("One-time premium against age")
#One-time premia based on unsmoothed mortalities for US males from 12.27.
lines(l_1:l_2,PI_1[l_1:l_2+1], col = "red",lty=2)
#### 12.31 by hand
#### 12.32
# a)
#given parameter values:
l=20;K=40;M=1;theta=c(0.000775,0.0000377,0.0927);psi=c(0.005,0.0000759,0.0875)
zeta=c(0.005,0,0)
#R-function used to compute kp: disableprob. (can be foud on STK4500 page)
kp=disableprob(theta,psi,zeta,l,K,M)$kp
#Plot of probability of an induvidual being in state active, disabled or dead
#K years ahead given age l:
matplot(0:K,t(kp[,1,]),"l",xlab="Years ahead",ylab="")
title("Fraction active, disabled and dead")
text(35,0.8,"active")
text(35,0.25,"disabled")
text(36,0.12,"dead")
# b) same as we did in a), but now with M=12
M=12
kp=disableprob(theta,psi,zeta,l,K,M)$kp
lines(0:K,kp[1,1,],col="blue")
lines(0:K,kp[2,1,],col="blue")
lines(0:K,kp[3,1,],col="blue")
#take a look at kp matrix if we got time.
## 07.03.18 ## 12.33-12.37 ## (only 12.33 and 12.34 are R tasks, rest by hand) ##
## 12.33
# a)
l_1=20;l_2=40
K=30
theta=c(0.000780,0.0000374,0.0928)
psi=c(0.005,0.00001,0.12)
zeta=c(0.0005,0,0)
kp=disableprob(theta,psi,zeta,l_1,K,M=12)$kp
plot(0:K,kp[1,1,],"l",xlab="Years ahead",ylab="",ylim=c(0.0,1))
#code below shows what line belongs to initial age l_1
#plot(0:K,kp[1,1,],"l",xlab="Years ahead",ylab="",ylim=c(0.0,1),col="red")
# Adding the other transition probabilities
for (l in (l_1+1):l_2){
kp=disableprob(theta,psi,zeta,l,K,M=12)$kp
lines(0:K,kp[1,1,])
}
title("Fraction active when initial age is varied")
#probability of staying active plotted against years ahead for induviduals from
#age l_1 to l_2 using the function disableporb.
#b) we do the same computations as above, but plot kP(2|1)
plot(0:K,kp[2,1,],"l",xlab="Years ahead",ylab="",ylim=c(0,0.5))
for (l in (l_1+1):l_2){
kp=disableprob(theta,psi,zeta,l,K,M=12)$kp
lines(0:K,kp[2,1,])
}
#-------------------------------------------------------------------------------------
## 12.34
#a) redo 12.32 b): (just copy paste mix of what we did in 12.32)
l=20;K=40
#NB: the values of thetha, psi ans zeta are not the same as in the book, done to see the differences clearer.
#we get the same result using the numbers given in the book but the fidderence is not as clear.
theta=c(0.000775,0.0000377,0.0927)
psi=c(0.005,0.0000759,0.0875)
zeta=c(0.005,0,0)
M=12
kp=disableprob(theta,psi,zeta,l,K,M)$kp
matplot(0:K,t(kp[,1,]),"l",xlab="Years ahead",ylab="")
#b)
zeta=c(0,0,0)
kp=disableprob(theta,psi,zeta,l,K,M)$kp
lines(0:K,kp[1,1,],col = "blue")
lines(0:K,kp[2,1,],col = "blue")
lines(0:K,kp[3,1,],col = "blue")
#c) compare the two evaluations
#-------------------------------------------------------------------------------------