e1210a=function()
{
# 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.
nn=matrix(scan("japusmale.txt"),byrow=T,ncol=2)
l_e=nrow(nn)-1
q=1-nn[1:l_e+1,]/nn[1:l_e,]
list(q=q,nn=nn,l_e=l_e)
}
e1210b=function()
{
# The mortalities in e1210a plotted againt age.
q=e1210a()$q
#Plotting
ll=1:nrow(q)-1
plot(ll,q[,1],xlab="Age (years)",ylab="",ylim=c(0,0.55))
lines(ll,q[,2],'l')
text(40,0.4,"Points: Japanese males 2008")
text(40,0.38,"Solid: US males 2007")
title("Survival probabilities per year")
}
e1210c=function()
{
# The mortalities in e1210a on log-scale plotted againt age.
q=e1210a()$q
#Plotting
ll=1:nrow(q)-1
plot(ll,log(q[,1]),xlab="Age (years)",ylab="Log")
lines(ll,log(q[,2]),'l')
text(40,-2,"Points: Japanese males 2008")
text(40,-2.6,"Solid: US males 2007")
title("Survival probabilities per year on LOG scale")
}
e1211a=function()
{
# 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.
nn=matrix(scan("japusfem.txt"),byrow=T,ncol=2)
l_e=nrow(nn)-1
q=1-nn[1:l_e+1,]/nn[1:l_e,]
list(q=q,nn=nn,l_e=l_e)
}
e1211b=function()
{
# The mortalities in e1211a plotted againt age.
q=e1211a()$q
#Plotting
ll=1:nrow(q)-1
plot(ll,q[,1],xlab="Age (years)",ylab="",ylim=c(0,0.55))
lines(ll,q[,2],'l')
text(40,0.4,"Points: Japanese females 2008")
text(40,0.38,"Solid: US females 2007")
title("Survival probabilities per year")
}
e1211c=function()
{
# The mortalities in e1211a on log-scale plotted againt age.
q=e1211a()$q
#Plotting
ll=1:nrow(q)-1
plot(ll,log(q[,1]),xlab="Age (years)",ylab="Log")
lines(ll,log(q[,2]),'l')
text(40,-2,"Points: Japanese females 2008")
text(40,-2.6,"Solid: US females 2007")
title("Survival probabilities per year on LOG scale")
}
e1212a=function(q,K)
{
# 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".
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)])
#kp is the life table
list(kp=kp,l_e=l_e)
}
e1212b=function(K=50)
{
# The life table "K" periods ahead plotted for Japanese males when "q"
# is the non-parameteric estimate from e1210a.
q=e1210a()$q[,1]
o=e1212a(q,K)
#Plotting
ind=1+(1:10)*min(floor(o$l_e/10),10)
matplot(0:K,o$kp[,ind],"l",xlab="Time ahead (years)",ylab="Survival probability")
title("Survival Japanese males: Initial age is 10,20,...,100 years")
}
e1212c=function(K=50)
{
# The life table "K" periods ahead plotted for for US males when "q"
# is the non-parameteric estimate from e1210a.
q=e1210a()$q[,2]
o=e1212a(q,K)
#Plotting
ind=1+(1:10)*min(floor(o$l_e/10),10)
matplot(0:K,o$kp[,ind],"l",xlab="Time ahead (years)",
ylab="Survival probability")
text(0.6*K,0.5*max(o$kp),"p")
title("Survival US males: Initial age is 10,20,...,100 years")
}
e1213a=function(q,K)
{
# 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".
#
o=e1212a(q,K)
p=c(1-q,rep(0,o$l_e+1))
kq=matrix(0,K+1,o$l_e+1)
for (l in 0:o$l_e+1)kq[1:K+1,l]=o$kp[1:K,l]*(1-p[1:K+(l-1)])
list(kq=kq,l_e=o$l_e)
}
e1213b=function(K=100)
{
# The mortalities in e1213a plotted for Japanese males when "q"
# is the non-parametric estimate from e1210a.
q=e1210a()$q[,1]
o=e1213a(q,K)
#Plotting
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 long for Japanse males 2008")
}
e1213c=function(K=100)
{
# The mortalities in e1213a plotted for US males when "q"
# is the non-parametric estimate from e1210a.
q=e1210a()$q[,2]
o=e1213a(q,K)
#Plotting
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 long for US males 2007")
}
e1214a=function()
{
# The Gomperz-Makeham model is fitted Japanese males 2008 and the
# fit checked through a plot.
#
# Output "theta_hat" is the fitted parameters.
o=e1210a()
# After the data and the un-smoothed probabilities,
# then the Gomperz-Makeham fitting.
ll=1:o$l_e-1
theta_hat=gompmakfit(o$nn[,1])$theta_hat
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,o$q[,1],xlab="Age (year)",ylab="Mortalities")
lines(ll,qgm)
# Text to plot.
text(40,0.4,"Dots: Data")
text(40,0.37,"Solid: Gomperz-Makeham")
title("Japanese males 2008")
list(theta_hat=theta_hat)
}
e1214b=function()
{
# The Gomperz-Makeham model is fitted US males 2007 and the
# fit checked through a plot.
#
# Output "theta_hat" is the fitted parameters.
o=e1210a()
ll=1:o$l_e-1
# After the data and the un-smoothed probabilities,
# then the Gomperz-Makeham fitting.
theta_hat=gompmakfit(o$nn[,2])$theta_hat
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,o$q[,2],xlab="Age (year)",ylab="Mortalities")
lines(ll,qgm)
# Text to plot
text(40,0.4,"Dots: Data")
text(40,0.37,"Solid: Gomperz-Makeham")
title("US males 2007")
list(theta_hat=theta_hat)
}
e1214c=function()
{
# The Gomperz-Makeham model is fitted Japanese females 2008 and the
# fit checked through a plot.
#
# Output "theta_hat" is the fitted parameters.
o=e1211a()
ll=1:o$l_e-1
# After the data and the un-smoothed probabilities,
# then the Gomperz-Makeham fitting.
theta_hat=gompmakfit(o$nn[,1])$theta_hat
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,o$q[,1],xlab="Age (year)",ylab="Mortalities")
lines(ll,qgm)
# Text to plot
text(40,0.4,"Dots: Data")
text(40,0.37,"Solid: Gomperz-Makeham")
title("Japanese females 2008")
list(theta_hat=theta_hat)
}
e1214d=function()
{
# The Gomperz-Makeham model is fitted US females 2007 and the
# fit checked through a plot.
#
# Output "theta_hat" is the fitted parameters.
o=e1211a()
ll=1:o$l_e-1
# After the data and the un-smoothed probabilities,
# then the Gomperz-Makeham fitting.
theta_hat=gompmakfit(o$nn[,2])$theta_hat
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,o$q[,2],xlab="Age (year)",ylab="Mortalities")
lines(ll,qgm)
# Text to plot
text(40,0.4,"Dots: Data")
text(40,0.37,"Solid: Gomperz-Makeham")
title("US females 2007")
list(theta_hat=theta_hat)
}
e1215a=function(l_e,theta,K)
{
# 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".
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)
}
e1215b=function(K=100,l_e=110,theta=c(0.000778,0.0000374,0.0928))
{
# Numerical application of e1215a with parameters those for US
# males 2007.
kq=e1215a(l_e,theta,K)$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")
}
e1216a=function(l_1=18,l_2=100,theta=c(0.000778,0.0000374,0.0928))
{
# 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.
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")
}
e1216b=function(l_1=18,l_2=100,theta=c(0.000778,0.0000374,0.0928))
{
# 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 increments for the
# the numerical implementation are "h=1", "h=0.1", "h=0.01".
#
# R-function used: gompmakmean.
E=rep(0,l_2+1)
#First h
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)")
#Second h
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")
#Third h
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")
text(l_1+0.7*(l_2-l_1),0.75*max(E),"Different increments")
text(l_1+0.7*(l_2-l_1),0.70*max(E),"for the numerical integration")
title("Remaining number of years")
}
e1216c=function(l_1=18,l_2=100,theta=c(0.000778,0.0000374,0.0928))
{
# Computes and plots the expected total number of years against age
# under the Gomperz-Makeham model with parameters "theta" (vector). The
# plot applies to age between "l_1" and "l_2". Time increments for the
# the numerical implementation are "h=1", "h=0.1", "h=0.01".
#
# R-function used: gompmakmean.
E=rep(0,l_2+1)
#First h
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)")
#Second h
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")
#Third h
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")
text(l_1+0.3*(l_2-l_1),0.95*max(E),"Different increments")
text(l_1+0.3*(l_2-l_1),0.935*max(E),"for the numerical integration")
title("Total number of years")
}
e1217c=function(theta,epsilon)
{
# Computes the percentiles of how long people live under a
# Gomperz-Makeham model with parameters "theta" (vector). "epsilon"
# (vector) defines which percentiles are computed.
#
# Output is "epsilon" and the percentiles "y_epsilon" (vector).
y_epsilon=epsilon
fe=function(y,e)-theta[1]*y-theta[2]/theta[3]*(exp(theta[3]*y)-1)-log(1-e)
for (i in 1:25){y_epsilon[i]=uniroot(fe,e=epsilon[i],c(0.01,120))[[1]]}
list(epsilon=epsilon,y_epsilon=y_epsilon)
}
e1217d=function(theta=c(0.000775,0.0000377,0.0928),epsilon=1:25*0.04-0.02)
{
# Numerical application of e1217c.
o=e1217c(theta,epsilon)
plot(o$epsilon,o$y_epsilon,ylab="Percentile (years)",xlab="Probabilities")
title("Percentiles under Gomperz-Makeham")
}
e1222a=function(kp,d,K_1)
# A table of mortality adjusted annuties (in advance) is computed for a
# life table "kp" and the discount "d". Duration "K" of the annuity varies
# from 0 to "K_1" and the individuals are between 0 and "l_e" of age.
#
# Output is a matrix "a" with coeffiecients for each initial age stored
# as columns.
#
# Warning: If e1212a is used to set up "kp", it must be called with K> K_1.
{
M=d**(0:(K_1-1))*kp[1:K_1,]
a=apply(M,2,cumsum)
list(a=a)
}
e1222b=function(d=0.97,K_1=50,l_1=25,l_2=70)
{
# Computes and plot mortality adjusted annuities (in advance) up to "K_1"
# for US males 2007 between "l_1" and "l_2" years of age. "d" is the discount.
q=e1210a()$q[,2]
o=e1212a(q,K_1+1)
#After the life table for USmales, now the annutites:
a=e1222a(o$kp,d,K_1)$a
matplot(1:K_1,a[,l_1:l_2+1],"l",xlab="Years ahead",ylab="Coefficients")
title("Mortality adjusted annuity coeffcients for US males")
}
e1223c=function(d=0.98,K=1:10*5)
{
# Computes and prints annuities up to different expiries defined in "K"
# (vector) with mortalities NOT taken into account. "d" is the discount.
a_K=round((1-d**K)/(1-d),digits=3)
a_K=cbind(K,a_K)
print(a_K)
}
e1223d=function(d=0.98,l=40,K=1:10*5)
{
# Computes and prints annuities (in advance) up to different expiries
# defined in "K" (vector). Mortalities are those for a US male 2007 of
# age "l". "d" is the discount.
K_1=max(K)
# Life table for US males first
q=e1210a()$q[,2]
o=e1212a(q,K_1+1)
# Mortality adjusted annuities here
M=d**(0:(K_1-1))*o$kp[1:K_1,l]
a=matrix(cumsum(M),K_1,1)
a=round(a[K,],digits=3)
b=cbind(K,a)
print(b)
}
e1224a=function(q,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(q)
kp=e1212a(q,l_e+1)$kp
a=e1222a(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)
}
e1224b=function(l_1=20,l_2=45,llr=55:75,s=1,d=0.97)
{
# 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".
q=e1210a()$q[,2]
kp=e1212a(q,111)$kp
a=e1222a(kp,d,110)$a
# After life table and annuities for US males, then premia lowest age.
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)
}
text(mean(llr),0.9*max(pi),"Initial age of indvidual varied")
title("Equivalence premium against retirement age")
}
e1225a=function(kq,d,K_1)
{
# A table of A-coeffients for term insurance is computed given a table
# of mortalities "kq"; consult e1213a. The discount is "d" and the
# the payment stream terminates at K which runs from 0 to "K_1". The
# individuals are between 0 and "l_e" of age.
#
# 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)
}
e1225b=function(d=0.97,K_1=50,l_1=25,l_2=55)
{
# Numerical application of e1225a for US males with A-coefficients
# plotted against duration of the scheme when age of the individuals
# are between "l_1" and "l_2".
q=e1210a()$q[,2]
o=e1213a(q,K_1+1)
A=e1225a(o$kq,d,K_1)$A
# After the A-coefficients, then the 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")
}
e1226a=function(s,a,A,K_1)
{
# Equivalence premia for term insurance are computed when contributions
# are in equal amounts in advance up to termination "K" of the contract
# or if the policy holder dies in which case a one-time payment "s" is
# is released. Annuity coefficents "a" and "A" (both matrices) are needed.
#
# On output "pi" (vector) contains 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)
}
e1226b=function(s=1,d=0.97,K_1=50,l_1=25,l_2=55)
{
# 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".
q=e1210a()$q[,2]
kp=e1212a(q,K_1+1)$kp
kq=e1213a(q,K_1+1)$kq
#The annuitiy coeficients
a=e1222a(kp,d,K_1)$a
A=e1225a(kq,d,K_1)$A
# The premia
pi=s*A[1:K_1,]/a[1:K_1,]
#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.021,"holder varied")
title("Equivalence premium for term insurance")
}
e1227a=function(q,s,l_r,d)
{
# Computes for a given mortality vector "q" one-time premia for pensions
# "s" starting at "l_r". The discount is "d". The age of the individuals
# are those in "q".
#
# 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=e1212a(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 tthe one-time premium.
ll=0:l_e
M=s*(d**ll)*kp*I
pi=apply(M,2,sum)
list(pi=pi)
}
e1227b=function(s=1,l_r=60,d=0.97,l_1=20,l_2=109)
{
# Numerical application of e1227a. One-time premia for
# individuals between "l_1" and "l_2" are plotted against age
# when the mortalities are those for US males 2007. The pension is
# "s", retirement age "l_r" and discount "d".
q=e1210a()$q[,2]
pi=e1227a(q,s,l_r,d)$pi
# After having omputed the one-time premia, then the plot.
plot(l_1:l_2,pi[l_1:l_2+1],"l",xlab="Age (years)",ylab="Premium")
title("One-time premium against age")
}
e1228a=function(l_e,theta,s,l_r,d)
{
# Computes for a given Gomperz-Makeham model with parameters
# "theta" (vector) one-time premia for pensions "s" starting at "l_r".
# The discount is "d". Age of the individuals are from 0 to "l_e".
#
# Output is a tabulation "pi" of all these one-time premia.
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)
}
e1228b=function(s=1,l_r=60,d=0.97,l_1=20,l_2=109,
theta=c(0.000780,0.0000374,0.0928))
{
# Numerical application of the Gomperz-Makeham model for
# US males 2007 with the one-time premia plotted for
# individuals between "l_1" and "l_2". Pension, retirement age
# and discount are "s", "l_r" and "d".
pi=e1228a(l_2+1,theta,s,l_r,d)$pi
# After one-time prmeia, then the plot.
plot(l_1:l_2,pi[l_1:l_2+1],"l",xlab="Age (years)",ylab="Premium")
title("One-time premium against age")
}
e1228c=function(s=1,l_r=60,d=0.97,l_1=20,l_2=109,
theta=c(0.000775,0.0000377,0.0927))
# Numerical application similar to e1228b, but now wih the
# Gomperz-Makeham mortalities for US male 2007 compared to
# results under the unsmoothed mortalities.
{
#One-time premia based on unsmoothed mortalities for US males
e1227b(s,l_r,d,l_1,l_2)
# After premia bsewd on the unsmoothed mortalities,
# then the Gomperz-Makeham version.
pi=e1228a(l_2+1,theta,s,l_r,d)$pi
lines(l_1:l_2,pi[l_1:l_2+1])
}
e1230a=function(kp,kpt,K)
{
# Computes the transition probabilities for a widow scheme given
# life tables "kp" for males and "kpt" for females. Both spouses
# are alive at time 0.
#
# Output are the probability vectors up to time "K" for both being alive
# ("p_00"), he dead/she alive ("p_10"), she dead/he alive ("p_01") and
# both dead ("p_11").
pm=kp[0:K+1,]
pf=kpt[0:K+1,]
p00=pm*pf
p10=(1-pm)*pf
p01=pm*(1-pf)
p11=(1-pm)*(1-pf)
list(p00=p00,p10=p10,p01=p01,p11=p11)
}
e1230b=function(l=30,K=50)
{
# Numerical application of e1230a for US males/females 2007. The
# four transition probabilities are plotted against time up to "K"
# when both spouses are of age "l" at the beginning.
om=e1210a()
of=e1211a()
# After the mortalities from fil, the transition probabilities for
# each gender:
okpm=e1212a(om$q[,2],K)
okpf=e1212a(of$q[,2],K)
# Now transition probabilities for couples:
o=e1230a(okpm$kp,okpf$kp,K)
# Plotting
plot(0:K,o$p00[,l+1],"l",ylim=c(0,1),xlab="Years ahead",ylab="Probabilities")
lines(0:K,o$p10[,l+1])
lines(0:K,o$p01[,l+1])
lines(0:K,o$p11[,l+1])
#text(10,0.7,"Which curve to which")
#text(10,0.66,"set of transitions?")
title("Transition probabilities for couples")
}
e1232a=function(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))
{
# Numerical application of R-function disableprob; consult it for
# details. The probability of state active, disabled, dead for an
# individual of age l are plotted against time up to "K". "M" is the
# technical parameter of the algorithm.
#
# R-function used: disableprob.
kp=disableprob(theta,psi,zeta,l,K,M)$kp
#Plotting
matplot(0:K,t(kp[,1,]),"l",xlab="Years ahead",ylab="")
title("Fraction active, disabled and dead")
}
e1232b=function(l=20,K=40)
{
# This is similar to e1232a, but now the results for M=1 amd M=12 are
# plotted jointly
#
# R-function used: disableprob.
# Starting for M=1
e1232a(l=l,K=K,M=1)
# Redoing for M=12
M=12
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,K,M)$kp
lines(0:K,kp[1,1,])
lines(0:K,kp[2,1,])
lines(0:K,kp[3,1,])
}
e1233a=function(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))
{
#
# Numerical application of R-function disableprob for
# transition probabilities for a disability scheme; consult
# disableprob for details. The probability of remaining in
# state active are plotted against time for individuals between
# "l_1" and "l_2".
#
# R-function used: disableprob.
# First set of transition probabilities
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))
title("Fraction active when initial age is varied")
# 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,])}
}
e1233b=function(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))
{
#
# Numerical application of R-function disableprob for
# transition probabilities for a disability scheme; consult
# disableprob for details. The probability of being in
# state disabled are plotted against time for individuals between
# "l_1" and "l_2" and in state active at the beginning.
#
# R-function used: disableprob.
# First set of transition probabilities
kp=disableprob(theta,psi,zeta,l_1,K,M=12)$kp
plot(0:K,kp[2,1,],"l",xlab="Years ahead",ylab="",ylim=c(0,0.5))
title("Fraction disabled when initial age is varied")
# 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[2,1,])}
}
e1234=function(l=20,K=40,theta=c(0.000780,0.0000374,0.0928),
psi=c(0.005,0.0000759,0.0875),zeta=c(0.005,0,0))
{
#
# Numerical application of R-function disableprob for the
# transition probabilities for a disability scheme; consult
# disableprob for details. The probability of being in
# states active, disabled and dead are plotted against time up to "K"
# for individuals of age "l". The results in e1232a are compared
# with those obtained under rehabilitation rate zero.
#
# R-function used: disableprob
# Starting as in Exercise 12.32a
e1232a(l=l,K=K,theta=theta,psi=psi,zeta=zeta)
# Redoing for new zeta
kp=disableprob(theta,psi,zeta=c(0,0,0),l,K,1)$kp
lines(0:K,kp[1,1,])
lines(0:K,kp[2,1,])
lines(0:K,kp[3,1,])
}
e1240a=function(q,s,l_r,d,N_sum,l_1,l_2,gamma,l_m)
{
# The present value of a pension portfolio is computed. The age profile
# between "l_1" and "l_2" is defined by "gamma" and "l_m"; see Exercise
# 12.40, and "N_sum" is the total number of individuals. Other input
# quantities are the vector of mortalities "q", the pension "s" per
# individual, the retirement age "l_r" and the discount "d".
#
# Output is present value "PV".
# First the age distribution of the portfolio.
N=exp(-gamma*abs(l_1:l_2-l_m))
N=(N_sum/sum(N))*N
# Then one-time premia.
pi=e1227a(q,s,l_r,d)$pi
# Finally the present value.
PV=sum(pi[l_1:l_2+1]*N)
list(PV=PV)
}
e1240b=function(s=1,l_r=60,N_sum=10000,l_1=20,l_2=100,d=0.97,gamma=0.1,
l_m=c(30,50,70))
{
# Numerical application of e1240a with the present value computed for
# three age profile as "l_m" is varied. Mortalities are those for US
# males 2007.
q=e1210a()$q[,2]
# After the mortalities have been read frm file, the present values
# are computed and printed.
PV=l_m
for (i in 1:3)PV[i]=round(e1240a(q,s,l_r,d,N_sum,l_1,l_2,gamma,l_m[i])$PV,
digits=0)
print(rbind(l_m,PV))
}
e1241a=function(q,s,l_r,N_sum,l_1,l_2,gamma,l_m)
{
# The liability sequence for the pension portfolio in e1240a is
# evaluated. The age profile between "l_1" and "l_2" is defined by
# "gamma" and "l_m"; see Exercise 12.40, and "N_sum" is the total
# number of individuals. Other input quantities are the mortality
# vector "q", the pension per individual "s" and the retirement age "l_r".
#
# On output "X" (vector) contains the liability sequence.
# First the life table.
l_e=length(q)
kp=e1212a(q,l_e)$kp
# Then the age distribution of portfolio
N=exp(-gamma*abs(l_1:l_2-l_m))
N=(N_sum/sum(N))*N
N=c(rep(0,l_1),N,rep(0,l_e-l_2))
# Finally the portfolio liabilities.
I=matrix(0,l_e+1,l_e+1)
I[row(I)+col(I)>l_r+1]=1
X=s*(kp*I)%*%N
list(X=X)
}
e1241b=function(s=1,l_r=60,N_sum=10000,l_1=20,l_2=100,gamma=0.1,
l_m=c(30,50,70),K=80)
{
# Numerical application of e1241a. The liabilities are plotted
# "K" years ahead for three different age distributions. Mortalities
# are those for US males 2007.
q=e1210a()$q[,2]
X=e1241a(q,s,l_r,N_sum,l_1,l_2,gamma,l_m[3])$X
plot(0:K,X[0:K+1],"l",xlab="Years ahead",ylab="Expense")
X=e1241a(q,s,l_r,N_sum,l_1,l_2,gamma,l_m[2])$X
lines(0:K,X[0:K+1])
X=e1241a(q,s,l_r,N_sum,l_1,l_2,gamma,l_m[1])$X
lines(0:K,X[0:K+1])
title("Portfolio liabilities for three different age profiles")
}
e1242b=function(s=1,l_r=60,N_sum=10000,l_1=20,l_2=100,d=0.97,gamma=0.1,
l_m=c(30,50,70))
{
# Another numerical application of e1241a where the present values
# of the liabilities in e1241b are tabulated for different age
# profiles in "l_m" (vector). The discount is "d".
q=e1210a()$q[,2]
# After the mortalities, the present values are computed and printed.
PV=l_m
for ( i in 1:3){
X=e1241a(q,s,l_r,N_sum,l_1,l_2,gamma,l_m[i])$X
kk=length(X)-1
PV[i]=round(sum(d**(0:kk)*X),digit=0)}
print(rbind(l_m,PV))
}
e1243a=function(q,pi,l_r,N_sum,l_1,l_2,gamma,l_m)
{
# Future premium income for the the pension portfolio in e1240a is
# evaluated. The age profile between "l_1" and "l_2" is defined by
# "gamma" and "l_m"; see Exercise 12.40, and "N_sum" is the total
# number of individuals. Other input quantities are the mortality
# vector "q", the premium per individual per period "pi" and the
# retirement age "l_r".
#
# Output is "Pi" (vector) containing future premium income.
# First the life table:
l_e=length(q)
kp=e1212a(q,l_e)$kp
# Then the age distribution of portfolio:
N=exp(-gamma*abs(l_1:l_2-l_m))
N=(N_sum/sum(N))*N
N=c(rep(0,l_1),N,rep(0,l_e-l_2))
# Finally the premium income:
I=matrix(0,l_e+1,l_e+1)
I[row(I)+col(I)