LAT4072_S7
Mauricio Rodriguez Abreu
10/2/2022
Multiple decrements example
setwd("C:/Users/23043/Dropbox/UDLAP/Cursos/2022 Primavera/Pensiones y SS/Datos")
data<-read.csv("MDM.csv")
head(data)## Age Inv Death DInv DNInv Ces Ret
## 1 0 0 0 0 0 0 0
## 2 1 0 0 0 0 0 0
## 3 2 0 0 0 0 0 0
## 4 3 0 0 0 0 0 0
## 5 4 0 0 0 0 0 0
## 6 5 0 0 0 0 0 0We can start by estimating the probabilities, remember that for four cases:
data$q1<-data$Inv*(1-0.5*(data$Death+data$Ces+data$Ret)+
(1/3)*(data$Death*data$Ces+data$Death*data$Ret+
data$Ces*data$Ret)-
(1/4)*(data$Death*data$Ces*data$Ret))
data$q2<-data$Death*(1-0.5*(data$Inv+data$Ces+data$Ret)+
(1/3)*(data$Inv*data$Ces+data$Inv*data$Ret+
data$Ces*data$Ret)-
(1/4)*(data$Inv*data$Ces*data$Ret))
data$q3<-data$Ces*(1-0.5*(data$Death+data$Inv+data$Ret)+
(1/3)*(data$Death*data$Inv+data$Death*data$Ret+
data$Inv*data$Ret)-
(1/4)*(data$Death*data$Inv*data$Ret))
data$q4<-data$Ret*(1-0.5*(data$Death+data$Ces+data$Inv)+
(1/3)*(data$Death*data$Ces+data$Death*data$Inv+
data$Ces*data$Inv)-
(1/4)*(data$Death*data$Ces*data$Inv))For \(l_x\), lets first estimate the total probability of leaving and the probability of surviving. Then, we can create a new dataset for our estimations:
data$qxt<-data$q1+data$q2+data$q3+data$q4
data$px<-1-data$qxt
data2<-data[which(data$Age>=60 & data$Age<=70),]Now, we can start with a \(l_{60}=50,000\) and work our estimations from there:
data2$lx<-50000
for (i in 2:11) {
data2$lx[i]<-data2$lx[i-1]*data2$px[i-1]
}And just for the sake of it
library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.1.2ggplot(data2,aes(x=Age,y=lx))+geom_line()
Now, we can estimate the population leaving for each cause:
data2$d1<-data2$lx*data2$q1
data2$d2<-data2$lx*data2$q2
data2$d3<-data2$lx*data2$q3
data2$d4<-data2$lx*data2$q4In order to get the total, we can use the function summaryBy from the doBy library:
library(doBy)## Warning: package 'doBy' was built under R version 4.1.2ans1<-summaryBy(d1+d2+d3+d4~1,data=data2,FUN=sum)And finally, the survivors are:
50000-ans1$d1.sum-ans1$d2.sum-ans1$d3.sum-ans1$d4.sum## [1] 279.6286The second problem is more complicated. We have to estimate annuities conditioned to something happening.
remove(ans1,data2)
dataM<-data
dataF<-read.csv("MDF.csv")
remove(data)First, lets define all possible annuities.But first, we have to define all life tables we need:
library(lifecontingencies)## Package: lifecontingencies
## Authors: Giorgio Alfredo Spedicato [aut, cre]
## (<https://orcid.org/0000-0002-0315-8888>),
## Christophe Dutang [ctb] (<https://orcid.org/0000-0001-6732-1501>),
## Reinhold Kainhofer [ctb] (<https://orcid.org/0000-0002-7895-1311>),
## Kevin J Owens [ctb],
## Ernesto Schirmacher [ctb],
## Gian Paolo Clemente [ctb] (<https://orcid.org/0000-0001-6795-4595>),
## Ivan Williams [ctb]
## Version: 1.3.7
## Date: 2021-03-21 22:00:02 UTC
## BugReport: https://github.com/spedygiorgio/lifecontingencies/issuesltM<-probs2lifetable(dataM$Death,radix=1, type="qx")
ltF<-probs2lifetable(dataF$Death,radix=1, type="qx")
ltMi<- probs2lifetable(dataM$DInv,radix=1,type="qx")
dataMb<-dataM[which(dataM$Age>=60 & dataM$Age<=70),]We can now, write functions to call when we need to estimate the benefits:
# Regular pension
axy<-function(r) {
if (r<65) {0} else {
12*15000*axn(ltM,x=r,i=0.03,k=12,payment="due")+
12*10000*(axn(ltF,x=r-2,i=0.03,k=12,payment="due")-
axyzn(tablesList=list(ltM,ltF),x=c(r,r-2),i=0.03,k=12,
status="joint",payment="due"))
}
}
# Early retirement
axy_er<- function(r) {
if(r<60) {0} else {
if (r==60) {
250000+0.75*12*15000*axn(ltM,x=r,i=0.03,k=12,payment="due")+
0.75*12*10000*(axn(ltF,x=r-2,i=0.03,k=12,payment="due")-
axyzn(tablesList=list(ltM,ltF),x=c(r,r-2),i=0.03,k=12,
status="joint",payment="due"))
} else {
if (r==61) {
200000+0.8*12*15000*axn(ltM,x=r,i=0.03,k=12,payment="due")+
0.8*12*10000*(axn(ltF,x=r-2,i=0.03,k=12,payment="due")-
axyzn(tablesList=list(ltM,ltF),x=c(r,r-2),i=0.03,k=12,
status="joint",payment="due"))
} else {
if (r==62) {
150000+0.85*12*15000*axn(ltM,x=r,i=0.03,k=12,payment="due")+
0.85*12*10000*(axn(ltF,x=r-2,i=0.03,k=12,payment="due")-
axyzn(tablesList=list(ltM,ltF),x=c(r,r-2),i=0.03,k=12,
status="joint",payment="due"))
} else {
if (r==63) {
100000+0.90*12*15000*axn(ltM,x=r,i=0.03,k=12,payment="due")+
0.9*12*10000*(axn(ltF,x=r-2,i=0.03,k=12,payment="due")-
axyzn(tablesList=list(ltM,ltF),x=c(r,r-2),i=0.03,k=12,
status="joint",payment="due"))
} else {
if (r==64){
50000+0.95*12*15000*axn(ltM,x=r,i=0.03,k=12,payment="due")+
0.95*12*10000*(axn(ltF,x=r-2,i=0.03,k=12,payment="due")-
axyzn(tablesList=list(ltM,ltF),x=c(r,r-2),i=0.03,k=12,
status="joint",payment="due"))
} else {0}}}}}}}
# pension for the wife
ay<- function(r) {
12*15000*axn(ltF,x=r-2,i=0.03,k=12,payment="due")
}
# For the disability before retirement
axi_y<-function(r) {
if (r<2) {0} else {
12*20000*(axn(ltMi,x=r,i=0.03,k=12,payment="due")+
0.75*(axn(ltF,x=r-2,i=0.03,k=12,payment="due")-
axyzn(tablesList=list(ltMi,ltF),x=c(r,r-2),i=0.03,k=12,
status="joint",payment="due")))
}
}Now, we have to estimate the benefit for each age. B1: Disability B2: Death B3: Early retirement B4: Retirement
B1<-c(60:70)
b<-rep(0,11)
B1<-cbind(B1,b)
for(i in 1:11){
B1[i,2]<-axi_y(B1[i,1])
}
B1<-B1[,2]
dataMb<-cbind(dataMb,B1)
B1<-c(60:70)
b<-rep(0,11)
B1<-cbind(B1,b)
for(i in 1:11){
B1[i,2]<-ay(B1[i,1])
}
B2<-B1[,2]
dataMb<-cbind(dataMb,B2)
B1<-c(60:70)
b<-rep(0,11)
B1<-cbind(B1,b)
for(i in 1:11){
B1[i,2]<-axy_er(B1[i,1])
}
B3<-B1[,2]
dataMb<-cbind(dataMb,B3)
B1<-c(60:70)
b<-rep(0,11)
B1<-cbind(B1,b)
for(i in 1:11){
B1[i,2]<-axy(B1[i,1])
}
B4<-B1[,2]
dataMb<-cbind(dataMb,B4)We can see on average how much money it would cost with Pr=1
summary(dataMb$B1)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4097616 4319679 4525649 4512787 4714940 4889336summary(dataMb$B2)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 3348008 3538488 3716172 3705452 3880015 4030685summary(dataMb$B3)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 1476098 3214970 3393428summary(dataMb$B4)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 3067339 1779097 3264421 3449157Now, for each age, we can multiply each benefit by age and estimate the expected value for the benefits we promised
dataMb$APV1<-dataMb$q1*dataMb$B1
dataMb$APV2<-dataMb$q2*dataMb$B2
dataMb$APV3<-dataMb$q3*dataMb$B3
dataMb$APV4<-dataMb$q4*dataMb$B4
dataMb$APVT<-dataMb$APV1+dataMb$APV2+dataMb$APV3+dataMb$APV4
summaryBy(APVT+APV1+APV2+APV3+APV4~1,data=dataMb,FUN=sum)## APVT.sum APV1.sum APV2.sum APV3.sum APV4.sum
## 1 13106220 82512.63 244768.5 5820610 6958329