#setwd("C:/Users/23043/Dropbox/UDLAP/Cursos/2022 Primavera/Pensiones y SS/Presentaciones")
setwd("~/Dropbox/UDLAP/Cursos/2022 Primavera/Pensiones y SS/Presentaciones")
lt<-read.csv("ILT1.csv")

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.6
## Date:     2019-03-05 23:50:03 UTC
## BugReport: http://github.com/spedygiorgio/lifecontingencies/issues
lt1<-probs2lifetable(lt$qx,radix=100000,type="qx",name="Life Table")

We will also need the salaray scale.

sc<-read.csv("SalaryScale.csv")

First, we will estimate the accumulation of benefits. We will not estimate the flat unit or career average since those are basic…

Now, assume we have a benefit of 1.5% of the 3-year average final salary for a person who starts working at age 30. The formula is: \[ B_r=k(r-y) \frac{1}{n}\sum_{t=r-n}^{r-1} s_t \] Assumign the person makes $1 at age 30…

average<- mean(sc[43:45,3])
k<-0.015
y<-30
r<-65
Br<-k*(r-y)*average
Br
## [1] 4.884775

Or using the second approach: \[ B_r=k(r-y) \frac{1}{n}(S_r - S_{r-n}) \]

sx<-sc[11:45,3]
sx[36]<-0
Sx<-rep(0,36)
for (i in 1:36) {
  Sx[i]<-sum(sx[1:i-1])
}
Br<-k*(r-y)*(1/3)*(Sx[36]-Sx[33])
Br
## [1] 4.884775

Now, for any age x we have \[ B_x=k(x-y) \frac{1}{n}(S_x - S_{x-n}) \]

x<-c(30:65)
Bx<-rep(0,36)
for (i in 1:36) {
  if(i==1) {
    Bx[i]<-k*(x[i]-y)*(1/3)*(Sx[i]-Sx[i])
  }
  else 
  {
    if(i==2) {
      Bx[i]<-k*(x[i]-y)*(1/1)*(Sx[i]-Sx[i-1])
    }
    else {
      if(i==3) {
        Bx[i]<-k*(x[i]-y)*(1/2)*(Sx[i]-Sx[i-2])
      }
      else {
        Bx[i]<-k*(x[i]-y)*(1/3)*(Sx[i]-Sx[i-3])
      }
    }
  }
}

Now, we can have estimate bx and create a dataset.

bx<-rep(0,36)
for (i in 1:35){
  bx[i]<-Bx[i+1]-Bx[i]
}
Accrual<-cbind(x,sx,Sx,Bx,bx)

Now, lets explore the modifications:

  1. Constant dollar modification:

\[ ^{CD}b_x=\frac{B_r}{r-y} \] \[ ^{CD}B_x=\frac{B_r}{r-y}(x-y) \quad for \quad y \leq x \leq r \]

Accrual<-as.data.frame(Accrual)
#COnstant dollar modification

Accrual$CDBx<-Br/(r-y)*(Accrual$x-y)
Accrual$CDbx<-Br/(r-y)

Now the constant percent modification:

\[ ^{CP}b_x=\frac{B_r}{S_r}s_x \] \[ ^{CP}B_x=\frac{B_r}{S_r}S_x \]

Accrual$CPbx<-(Br/Sx[36])*Accrual$sx
Accrual$CPBx<-(Br/Sx[36])*Accrual$Sx

Now, we can compare they way the benefits accumulate throughout the years:

library(ggplot2)
## Registered S3 methods overwritten by 'tibble':
##   method     from  
##   format.tbl pillar
##   print.tbl  pillar
A64<-Accrual[which(Accrual$x<=64),]

ggplot(A64, aes(x=x))+
  geom_line(aes(y=bx),color="red")+
  geom_line(aes(y=CDbx),color="blue")+
  geom_line(aes(y=CPbx),color="green")

ggplot(Accrual, aes(x=x))+
  geom_line(aes(y=Bx),color="red")+
  geom_line(aes(y=CDBx),color="blue")+
  geom_line(aes(y=CPBx),color="green")

Section for Actuarial Liabilities

Plan termination liability

\[ PTL_x=\begin{cases} B_x \cdot {_{r-x}p_x^{(m)}} \cdot v^{r-x} \cdot \ddot{a}_r & x\leq r \\ B_r \cdot \ddot{a}_x& x \geq r \end{cases} \]

Since we only have individuals younger that 65 in our dataset, we dont need the second part. Estimating this:

Accrual$PTL<-Accrual$Bx*pxt(lt1,x=Accrual$x,t=r-Accrual$x)*
  (1.035^(-(r-Accrual$x)))*axn(lt1,x=r,i=0.035,payment="due")

And we can explore the behavior of this liability:

ggplot(Accrual,aes(x=x,y=PTL))+geom_line()

Now, for the PVFB at age x we use the Benefit at retirement:

Accrual$PVFB<-Br**pxt(lt1,x=Accrual$x,t=r-Accrual$x)*
  (1.035^(-(r-Accrual$x)))*axn(lt1,x=r,i=0.035,payment="due")

We need to find the value for k under different methods:

  1. Accrued Benefit method
Accrual$ABk<-Accrual$Bx/Br
  1. Benefit prorate methods b.1) Constant dollar method:
Accrual$BP1k<-(Accrual$x-y)/(r-y)

b.2) Constant percenr method

Accrual$BP2k<-Accrual$Sx/Accrual$Sx[36]
  1. Cost Prorate Methods c.1) Constant Dollar
Accrual$ax1<- axn(lt1,x=y,n=Accrual$x-y,i=0.035,payment="due")
Accrual$ar1<- axn(lt1,x=y,n=r-y,i=0.035,payment="due")
Accrual$CP1k<-Accrual$ax1/Accrual$ar1

c.2) Constant percent

# Pendiente

Now, we can estimate the different actuarial liabilities:

Accrual$ALAB<-Accrual$PVFB*Accrual$ABk
Accrual$ALBP1<-Accrual$PVFB*Accrual$BP1k
Accrual$ALBP2<-Accrual$PVFB*Accrual$BP2k
Accrual$ALCP1<-Accrual$PVFB*Accrual$CP1k

And just plot them…

ggplot(Accrual, aes(x=x))+
  geom_line(aes(y=ALAB),color="red")+
  geom_line(aes(y=ALBP1),color="blue")+
  geom_line(aes(y=ALBP2),color="green")+
  geom_line(aes(y=ALCP1),color="purple")