# Reinterpreting Lee-Carter Mortality Model

Last week, while I was giving my crash course on R for insurance, we’ve been discussing possible extensions of Lee & Carter (1992) model. If we look at the seminal paper, the model is defined as follows

# Smoothing mortality rates

This morning, I was working with Julie, a student of mine, coming from Rennes, on mortality tables. Actually, we work on genealogical datasets from a small region in Québec, and we can observe a lot of volatiliy. If I borrow one of her graph, we get something like

Since we have some missing data, we wanted to use some Generalized Nonlinear Models. So let us see how to get a smooth estimator of the mortality surface.  We will write some code that we can use on our data later on (the dataset we have has been obtained after signing a lot of official documents, and I guess I cannot upload it here, even partially).

DEATH <- read.table(
"http://freakonometrics.free.fr/Deces-France.txt",
"http://freakonometrics.free.fr/Exposures-France.txt",
library(gnm)
D=DEATH$Male E=EXPO$Male
A=as.numeric(as.character(DEATH$Age)) Y=DEATH$Year
I=(A<100)
base=data.frame(D=D,E=E,Y=Y,A=A)
subbase=base[I,]
subbase=subbase[!is.na(subbase$A),] The first idea can be to use a Poisson model, where the mortality rate is a smooth function of the age and the year, something like $D_{x,t}\sim\mathcal{P}(E_{x,t}\cdot \exp[{\color{blue}s(x,t)}])$that can be estimated using library(mgcv) regbsp=gam(D~s(A,Y,bs="cr")+offset(log(E)),data=subbase,family=quasipoisson) predmodel=function(a,y) predict(regbsp,newdata=data.frame(A=a,Y=y,E=1)) vX=trunc(seq(0,99,length=41)) vY=trunc(seq(1900,2005,length=41)) vZ=outer(vX,vY,predmodel) persp(vZ,theta=-30,col="green",shade=TRUE,xlab="Ages (0-100)", ylab="Years (1900-2005)",zlab="Mortality rate (log)") The mortality surface is here It is also possible to extract the average value of the years, which is the interpretation of the $a_x$ coefficient in the Lee-Carter model, predAx=function(a) mean(predict(regbsp,newdata=data.frame(A=a, Y=seq(min(subbase$Y),max(subbase$Y)),E=1))) plot(seq(0,99),Vectorize(predAx)(seq(0,99)),col="red",lwd=3,type="l") We have the following smoothed mortality rate Recall that the Lee-Carter model is $D_{x,t}\sim\mathcal{P}(E_{x,t}\cdot \exp[{\color{blue}a_x+b_x\cdot k_t}])$ where parameter estimates can be obtained using regnp=gnm(D~factor(A)+Mult(factor(A),factor(Y))+offset(log(E)), data=subbase,family=quasipoisson) predmodel=function(a,y) predict(regnp,newdata=data.frame(A=a,Y=y,E=1)) vZ=outer(vX,vY,predmodel) persp(vZ,theta=-30,col="green",shade=TRUE,xlab="Ages (0-100)", ylab="Years (1900-2005)",zlab="Mortality rate (log)") The (crude) mortality surface is with the following $a_x$ coefficients. plot(seq(1,99),coefficients(regnp)[2:100],col="red",lwd=3,type="l") Here we have a lot of coefficients, and unfortunately, on a smaller dataset, we have much more variability. Can we smooth our Lee-Carter model ? To get something which looks like $D_{x,t}\sim\mathcal{P}(E_{x,t}\cdot \exp[{\color{blue}s_a(x)+s_b(x)\cdot s_k(t)}])$ Actually, we can, and the code is rather simple library(splines) knotsA=c(20,40,60,80) knotsY=c(1920,1945,1980,2000) regsp=gnm(D~bs(subbase$A,knots=knotsA,Boundary.knots=range(subbase$A),degre=3)+ Mult(bs(subbase$A,knots=knotsA,Boundary.knots=range(subbase$A),degre=3), bs(subbase$Y,knots=knotsY,Boundary.knots=range(subbase$Y),degre=3))+ offset(log(E)),data=subbase, family=quasipoisson) BpA=bs(seq(0,99),knots=knotsA,Boundary.knots=range(subbase$A),degre=3)
BpY=bs(seq(min(subbase$Y),max(subbase$Y)),knots=knotsY,Boundary.knots= range(subbase$Y),degre=3) predmodel=function(a,y) predict(regsp,newdata=data.frame(A=a,Y=y,E=1)) v Z=outer(vX,vY,predmodel) persp(vZ,theta=-30,col="green",shade=TRUE,xlab="Ages (0-100)", ylab="Years (1900-2005)",zlab="Mortality rate (log)") The mortality surface is now and again, it is possible to extract the average mortality rate, as a function of the age, over the years, BpA=bs(seq(0,99),knots=knotsA,Boundary.knots=range(subbase$A),degre=3)
Ax=BpA%*%coefficients(regsp)[2:8]
plot(seq(0,99),Ax,col="red",lwd=3,type="l")

We can then play with the smoothing parameters of the spline functions, and see the impact on the mortality surface

knotsA=seq(5,95,by=5)
knotsY=seq(1910,2000,by=10)
regsp=gnm(D~bs(A,knots=knotsA,Boundary.knots=range(subbase$A),degre=3)+ Mult(bs(A,knots=knotsA,Boundary.knots=range(subbase$A),degre=3),
bs(Y,knots=knotsY,Boundary.knots=range(subbase$Y),degre=3)) +offset(log(E)),data=subbase,family=quasipoisson) predmodel=function(a,y) predict(regsp,newdata=data.frame(A=a,Y=y,E=1)) vZ=outer(vX,vY,predmodel) persp(vZ,theta=-30,col="green",shade=TRUE,xlab="Ages (0-100)", ylab="Years (1900-2005)",zlab="Mortality rate (log)") We now have to use those functions our our small data sample ! That should be fun…. # Combien de temps profite-t-on de ses grands parents ? Ce Hier matin, je suis tombé un peu par hasard sur deux graphiques de l’INSEE (en France) avec l’age moyen des mères à l’accouchement, en fonction de rang de naissance de l’enfant, avec tout d’abord 1905-1965, puis 1960-2000 Ces graphiques sont passionnants en soi – comme en ont témoigné pas mal de followers sur Twitter – mais ils m’ont fait m’interroger. En particulier, sur la croissance observée depuis 30 ans, qui me faisait penser à la tendance croissante observée sur les durées de vie. On n’a – malheureusement – pas accès au données complètes sur le site, mais on peut trouver d’autres donnéesintéressantes (en l’occurrence l’age moyen à la naissance). > agenaissance=read.table("http://freakonometrics.blog.free.fr/ public/data/agenaissance.csv",header=TRUE,sep=",") > agenaissance$Age=as.character(agenaissance$AGE) > agenaissance$AGE=as.numeric(substr(agenaissance$Age,1,2))+ + as.numeric(substr(agenaissance$Age,4,4))/10
> plot(agenaissance$ANNEE+.5,agenaissance$AGE,
+ type="l",lwd=2,col="blue")

Visuellement, on retrouve la courbe en bleu foncée sur les graphiques ci-dessus,

On peut alors aller en cran plus loin, en se demandant non pas quel était l’âge moyen de la mère, mais de la grand-mère (au sens la mère de la mère)

> agenaissance$NAIS.MERE=(agenaissance$ANNEE+.5)-
+ agenaissance$AGE > w=(trunc(agenaissance$NAIS.MERE-.5))
> rownames(agenaissance)=agenaissance$ANNEE > a1=agenaissance[as.character(w),]$NAIS.MERE
> a2=agenaissance[as.character(w+1),]$NAIS.MERE > p=agenaissance$NAIS.MERE-(w+.5)
> agenaissance$NAIS.GRD.MERE=(1-p)*a1+p*a2 > agenaissance$age.GRD.MERE=agenaissance$ANNEE+.5- + agenaissance$NAIS.GRD.MERE
> tail(agenaissance)
ANNEE  AGE   Age NAIS.MERE NAIS.GRD.MERE age.GRD.MERE
2000  2000 30.3 30,3     1970.2       1942.87        57.63
2001  2001 30.4 30,4     1971.1       1943.80        57.70
2002  2002 30.4 30,4     1972.1       1944.92        57.58
2003  2003 30.5 30,5     1973.0       1945.95        57.55
2004  2004 30.5 30,5     1974.0       1947.05        57.45
2005  2005 30.6 30,6     1974.9       1948.04        57.46
> plot(agenaissance$ANNEE+.5,agenaissance$age.GRD.MERE,
+ type="l",lwd=2,col="red")

Là encore, on peut visualiser l’âge de la grand-mère maternelle à la naissance

A partir de là, on peut se demander combien de temps on profite de ses grands-parents (ou tout du moins ici de sa grand mère maternelle), en se basant sur les calculs d’espérance de vie résiduelle. En utilisant le modèle de Lee-Carter pour modéliser les taux de décès annuel, et en extrapolant sur le siècle en cours, on peut extrapoler les espérances de vie résiduelles.

> Deces <- read.table("http://freakonometrics.free.fr/
> Deces$Age <- as.numeric(as.character(Deces$Age))
> Deces$Age[is.na(Deces$Age)] <- 110
> Expo$Age <- as.numeric(as.character(Expo$Age))
> Expo$Age[is.na(Expo$Age)] <- 110
>  library(forecast)
>  library(demography)
>  YEAR <- unique(Deces$Year);nC=length(YEAR) > AGE <- unique(Deces$Age);nL=length(AGE)
>  MUF  <- matrix(Deces$Female/Expo$Female,nL,nC)
>  POPF <- matrix(Expo$Female,nL,nC) > BASEF <- demogdata(data=MUF, pop=POPF,ages=AGE, + years=YEAR, type="mortality", + label="France", name="Femmes", lambda=1) > LCF <- lca(BASEF) > LCFf<-forecast(LCF,h=100) > A <- LCF$ax
> B <- LCF$bx > K1 <- LCF$kt
> K2 <- K1[length(K1)]+LCFf$kt.f$mean
> K <- c(K1,K2)
> 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]) }}
> esp.vie = function(xentier,T){
+ s <- seq(0,99-xentier-1)
+ MUd <- MU[xentier+1+s,T+s-1898]
+ Pxt <- cumprod(exp(-diag(MUd)))
+ ext <- sum(Pxt)
+ return(ext) }
> EVIE = function(x,T){
+ x1 <- trunc(x)
+ x2 <- x1+1
+ return((1-(x-x1))*esp.vie(x1,T)+(x-x1)*esp.vie(x2,T)) }
> agenaissance$EV=NA > for(i in 1:100){ + t <- 2006-i + agenaissance$EV[agenaissance$ANNEE==t]= + EVIE(x=agenaissance$age.GRD.MERE[
+ agenaissance$ANNEE==t],t) } > tail(agenaissance) ANNEE AGE Age NAIS.MERE NAIS.GRD.MERE age.GRD.MERE EV 2000 30.3 30,3 1970.2 1942.87 57.63 29.13876 2001 30.4 30,4 1971.1 1943.80 57.70 29.17047 2002 30.4 30,4 1972.1 1944.92 57.58 29.39027 2003 30.5 30,5 1973.0 1945.95 57.55 29.52041 2004 30.5 30,5 1974.0 1947.05 57.45 29.72511 2005 30.6 30,6 1974.9 1948.04 57.46 29.80398 Autrement dit, sur la dernière ligne, l’espérance de vie (résiduelle) pour une femme de 57.46 ans en 2005 était d’environ 29.80 ans. On peut alors visualiser non seulement l’âge moyen de sa grand-mère à la naissance, mais son espérance de vie résiduelle, > plot(agenaissance$ANNEE+.5,agenaissance$EV, + type="l",lwd=2,col="purple") On note que depuis 30 ans, en France, la durée (moyenne) pendant laquelle les petits-enfants vont profiter de leur grands parents s’est stabilisé à une trentaine d’années. On peut aussi continuer, et remonter d’un cran (en refaisant tourner le code avec quelques modifications): on a alors l’âge (moyen) de son arrière grand mère à la naissance, et la durée de vie (résiduelle) des arrière grand mères On manque ici de données, mais il semble que l’on profite – en moyenne – environ 5 ans de son arrière grand mère. Maintenant on peut aussi s’interroger sur les limites de cette étude rapide. En particulier, de même qu’il existe une corrélation forte entre les durées de vie de conjoints (e.g. broken heart syndrom de Jagger & Sutton (1991)), on peut se demander si la naissance d’enfants et de petits-enfants a un impact sur la durée de vie résiduelle d’une personne (ou si on peut supposer l’indépendance comme on l’a fait ici). # Is it that stupid to make extremely long term forecast when studying mortality ? I received recently a comment by FCA (here) who raised an important question, about forecast in dynamic mortality models. (S)he mentioned that from his(her) point of view, the econometric models I considered were “good to predict for the next, say, 3 or 4 years. Not for the next 50 years…”. Which was the message I tried to stress last year in a conference about retirement in France (here). But from a quantitativepoint of view, how inconsistent were forecasts made 35 years ago, or 60 years ago ? Consider here the Lee Carter model, obtained on the periods 1816-1950 (in black below), 1816-1975 (in red) and 1816-2000 (in blue), unfortunately, it is difficult to compare ‘s since we have identifiability problems here. Nevertheless, we if consider affine transformation so that ‘s are equal in 1900 and 1950 (say), we obtain On that graph, we considered an ETS (AAN) forecast. If we do not consider the entire series for forecasting, but only observations following WWI (1945), we obtain For sketches of the R code, T=1980 base0=data.frame(D,E,A,Y,a=as.factor(A), y=as.factor(Y)) base=base0[base0$Y<=T,]
LC2=gnm(D~a+Mult(a,y),offset=log(E),family=
poisson,data=base)
A=LC2$coefficients[1]+LC2$coefficients[2:110]
B=LC2$coefficients[111:220] K0=LC2$coefficients[221:length(LC2$coefficients)] Y=as.numeric(K0) K1=c(K0,forecast(ets(Y,model="AAN"),h=240)$mean)
K2=c(K0,forecast(auto.arima(Y,allowdrift=TRUE),h=240)$mean) MU=matrix(NA,length(A),length(K1)) MU1=MU2=MU for(i in 1:length(A)){ for(j in 1:length(K1)){ MU1[i,j]=exp(A[i]+B[i]*K1[j]) MU2[i,j]=exp(A[i]+B[i]*K2[j]) }} x=40 s=seq(0,109-x-1) t=2000 Pxt1=cumprod(exp(-diag(MU1[x+1+s,t+s-base1$Year[1]-1])))
Pxt2=cumprod(exp(-diag(MU2[x+1+s,t+s-base1$Year[1]-1]))) r=.035 m=70 h=seq(0,39) V1=1/(1+r)^(m-x+h)*Pxt1[m-x+h] V2=1/(1+r)^(m-x+h)*Pxt2[m-x+h] M=cbind(V1,V2) apply(M,2,sum) Actually, it is not that bad…. even if it is only a qualitative intuition. Again, I am not a demographer, and my interest is more on actuarial science… so if we look at the estimation of annuities (still the same insurance contract, as here) for some insured of age 40 in 2000, we get the following graph (where forecasts ‘s were obtained on the complete series, i.e. from 1816 until the year we consider), (here it means that in 1900, I had to forecast mortality for someone of age 40 in 2000… so we had to forecast mortality with a 150 year horizon). Obviously, even if we are able to forecast improvement of mortality rates, it is not enough since it looks like, each year, improvement are alway higher than what what expected. Note that if we run it twice (since there might be problem with initial values in the econometric procedure) we obtain something similar, So, the output is consistent. And if we change the way we predict future values, e.g. on focusing only on the past 50 years, i.e. K1=c(K0,forecast(ets(Y[(length(Y)-50):length(Y)], model="AAN"),h=240)$mean)
K2=c(K0,forecast(auto.arima(Y[(length(Y)-50):length(Y)],
allowdrift=TRUE),h=240)\$mean)

we obtain the following graph for the annuity associated to an insurance contract sold in 2000,

so that relative changes compared with 1980 are (in %)

Hence, over a bit more than 25 years, we underestimated annuities of 25%. We if start to take into account possible investments, it is not so bad, I think….  don’t you think ?