Statistique de l’assurance STT6705V, partie 12

i.e.

Thus,

Complete expectation of life is then

or its discrete version

Hence, it is possible to write

which can be extended in the dynamic framework as follows

A natural estimator of that quantity is

i.e., with Lee-Carter model

All those quantities can be computed quite simply.

But first, I have to work a little bit on the dataset, at least to be able to predict mortality up to 100 years old (with the demography package, we have to stop at 90). One idea can be to mix two estimation techniques: the nonlinear Poisson regression to get and up to 99 years old, and then to use Rob Hyndman’s package to estimate and predict the component. But first, we have to check that  ‘s and ‘s with the two techniques are not too different. The code for the first model is

BASEB=BASEB[,1:99]
BASEC=BASEC[,1:99]
AGE=AGE[1:99]
library(gnm)
D=as.vector(BASEB)
E=as.vector(BASEC)
A=rep(AGE,each=length(ANNEE))
Y=rep(ANNEE,length(AGE))
base=data.frame(D,E,A,Y,a=as.factor(A),
y=as.factor(Y))
LC2=gnm(D~a+Mult(a,y)+offset=log(E),
family=poisson,data=base)

BASEB=BASEB[,1:90]
BASEC=BASEC[,1:90]
AGE=AGE[1:90]
library(demography)
MU=as.matrix(BASEB/BASEC)
base=demogdata(data=t(MU),pop=t(BASEC),
ages=AGE,years=ANNEE,type="mortality",
label="France",name="Total",lambda=0)
LC3=lca(base,interpolate=TRUE)
LC3F=forecast(LC3,150)

for the second one. Then, we can compare predictions for  ‘s e.g. (with output from econometric regression in blue, and Rob’s package in red)

plot(AGE,LC3$ax,lwd=2,col="red",xlab="",ylab="")
lines(1:98,LC2$coefficients[1]+LC2$coefficients[2:99],
lwd=2,col="blue")

The second problem is that we have to use a linear transformation, since in the econometric package, we do not use the same constraints as in Rob’s package. Here, it was

plot(ANNEE,-LC2$coefficients[199:297]*50)
lines(ANNEE,LC3$kt,col="red")

Then we can compute predicted for all ages and years,

A=LC2$coefficients[1]+LC2$coefficients[2:99]
B=LC2$coefficients[100:190]
K1=LC3$kt
K2=LC3$kt[99]+LC3F$kt.f$mean
K=-c(K1,K2)/50
MU=matrix(NA,length(A),length(K))
for(i in 1:length(A)){
for(j in 1:length(K)){
MU[i,j]=exp(A[i]+B[i]*K[j])
}}

t=2000
x=40
s=seq(0,99-x-1)
MUd=MU[x+1+s,t+s-1898]
Pxt=cumprod(exp(-diag(MUd)))

x=40
E=rep(NA,150)
for(t in 1900:2040){
s=seq(0,90-x-1)
MUd=MU[x+1+s,t+s-1898]
Pxt=cumprod(exp(-diag(MUd)))
ext=sum(Pxt)
E[t-1899]=ext}
plot(1900:2049,E)

x=40
t=2000
r=.035
m=70
h=seq(0,21)
V=1/(1+r)^(m-x+h)*Pxt[m-x+h]
sum(V)
VV=rep(NA,150)

for(t in 1900:2040){
s=seq(0,90-x-1)
MUd=MU[x+1+s,t+s-1898]
Pxt=cumprod(exp(-diag(MUd)))
h=seq(0,30)
V=1/(1+r)^(m-x+h)*Pxt[m-x+h]
VV[t-1899]=sum(V,na.rm=TRUE)}
plot(1900:2049,VV)