# 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).

"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…. # Course on risk, insurance, and uncertainty The course on risk and insurance in Luminy, starts at 10.30 on Friday (here) instead of Thursday (I switched with Patrice Bertail). The slides can be found here, Then, it will be time to leave Marseille, # Modélisation des tables de mortalité En cours d’actuariat, nous commencerons à parler des tables de mortalité prospectives. Côté théorie, je peux renvoyer vers quelques slides de Pierre Devolder, ici. Je renvoie ici pour un court papier de vulgarisation sur le sujet. Sur l’ajustement de lois de mortalité prospectives, je renvoie au très beau mémoire en ligne ici (ou pour des choses très proches en anglais). Pour l’approche statique, je renvoie également ici pour une méthode simple d’ajustement de loi de Makeham. Pour les aspects dynamique, il y a des documents techniques ici ou . Sinon pour les données que nous utiliserons en TD, elles sont ici (pour l’exposition, i.e. le nombre de vivants) et (pour les décès). Il faut retravailler un peu les données, à cause d’une valeur “110+”, mais rien de plus simple, > DEATH=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/R/FR-Deaths_1x1.txt",header=TRUE) > EXPOSURE=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/R/FR-Exposures_1x1.txt",header=TRUE) > DEATH$Age=as.numeric(as.character(DEATH$Age)) Message d'avis : NAs introduits lors de la conversion automatique > DEATH$Age[is.na(DEATH$Age)]=110 > EXPOSURE$Age=as.numeric(as.character(EXPOSURE$Age)) Message d'avis : NAs introduits lors de la conversion automatique > EXPOSURE$Age[is.na(EXPOSURE$Age)]=110 Coté code, je renvoie à un précédant billet sur le sujet, ici, et pour le package (non officiel) de démographie développé par Rob Hyndman. Enfin, à la fin du cours, on parlera un peu de modèles de durée. J’essayerai de taper quelques choses dans les jours à venir, mais côté TD, on utilisera ce que fait Frédéric Planchet (qui connaît fort bien le sujet), et qui se trouve en ligne ici. Les données sont en ligne ici, avec un peu de code R pour retraiter les donnéesici et , et un peu de code pour des modèles classiques, comme du modèle de Cox (ici) ou du Kaplan Meier (). Les codes R utilisés en TD sont les suivants (avec en plus les graphiques obtenus), > D=DEATH[DEATH$Year==an,]
> E=EXPOSURE[EXPOSURE$Year==an,] > MU = D[,3:5]/E[,3:5] > plot(0:110,log(MU[,1]),type="l",col="red") > lines(0:110,log(MU[,2]),col="blue") Avec ci-dessus les log des taux de décès instantannés pour l’année 1986. > PH=matrix(NA,111,111) # ligne X, colonne H > PF=PH > for(x in 0:110){ + PH[x+1,1:(111-x)]=exp(-cumsum(MU[(x+1):111,2])) + PF[x+1,1:(111-x)]=exp(-cumsum(MU[(x+1):111,1])) + } > x=0 > plot(1:111,PH[x+1,],ylim=c(0,1),type="l",col="blue") > lines(1:111,PF[x+1,],col="red") pour la fonction de survie à la naissance (ci-dessus) ou à 50 ans (ci-dessous), > somme=function(x){sum(x,na.rm=TRUE)} > EH=EF=rep(NA,111) > EH=apply(PH,1,somme) > EF=apply(PF,1,somme) > plot(0:110,EH,type="l",col="blue") > lines(0:110,EF,col="red") Ce qui correspond aux espérances de vie résiduelles. > for(k in 1:3){ + Z=MU[,k] + Z[is.nan(Z)==TRUE]=NA + MU[,k]=Z + } > AGE=unique(DEATH$Age)
> ANNEE=unique(DEATH$Year) > MUF=matrix(MU[,1],length(AGE),length(ANNEE)) > MUH=matrix(MU[,2],length(AGE),length(ANNEE)) > persp(AGE[1:100],ANNEE,log(MUH[1:100,]), + theta=-30,col="light green",shade=TRUE) > persp(AGE[1:100],ANNEE,log(MUH[1:100,107:1]), + theta=-30,col="light green",shade=TRUE) > image(AGE,ANNEE,log(MUH)) On a retourné le graph en mettant les années à l’envers, mais on peut aussi regarder la surface pour les femmes, ainsi que l’allure des courbes de niveau > persp(AGE[1:100],ANNEE,log(MUF[1:100,107:1]), + theta=-30,col="light green",shade=TRUE) > image(AGE,ANNEE,log(MUF)) Mais nous n’avons ici qu’une lecture statique des tables, et nous ne suivons pas vraiment un individu: en même temps que la personne vieilli, le temps passe…. Considérons une personne née en 1900, et regardons (ex post) l’évolution de sa probabilité de décès (en regardant non pas la probabilité de décéder à 18 d’une personne de 18 ans en 1900, mais en tenant compte du fait que la personne a fêté ses 18 ans en 1918). > MUFP=MUHP=rep(NA,111) # MU PROSPECTIF > for(i in 1:111){ + MUHP[i]=MUH[i,colonne+i] + MUFP[i]=MUF[i,colonne+i] + } > plot(0:110,log(MUHP),col="purple",type="l") > D=DEATH[DEATH$Year==annee.naissance,]
> E=EXPOSURE[EXPOSURE$Year==annee.naissance,] > MU = D[,3:5]/E[,3:5] > lines(0:110,log(MU[,2]),col="blue") On peut aussi l’ajustement par un modèle à la Lee-Carter (via la fonction fournie par LifeMetrics) > source("https://perso.univ-rennes1.fr/arthur.charpentier/fitModels.r") > DEATH=DEATH[DEATH$Age<90,] # virer age>90
> EXPOSURE=EXPOSURE[EXPOSURE$Age<90,] # virer age>90 > XV=unique(DEATH$Age)
> YV=unique(DEATH$Year) > ETF=t(matrix(EXPOSURE[,3],length(XV),length(YV))) > DTF=t(matrix(DEATH[,3],length(XV),length(YV))) > ETH=t(matrix(EXPOSURE[,4],length(XV),length(YV))) > DTH=t(matrix(DEATH[,4],length(XV),length(YV))) > WA=matrix(1,length(YV),length(XV)) > LCF = fit701(xv=XV,yv=YV,etx=ETF,dtx=DTF,wa=WA) -1166100 -> -263940.2 ->-204175.2 -188375.9 -> -185822.4 ->-184403.5 -183223.3 -> -183069.8 ->-182350.3 -181701.1 -> -181600.6 ->-181145.5 -180729.6 -> -180671.5 ->-180371.8 -180101.9 -> -180068.6 ->-179871.2 > LCH = fit701(xv=XV,yv=YV,etx=ETH,dtx=DTH,wa=WA) -2564403 -> -1028776 ->-953629.9 -881334.2 -> -868279.3 ->-834413.6 -804946.8 -> -802817.2 ->-786343 -772775 -> -772007.5 ->-764501.7 -758653.8 -> -758399 ->-755212 -752931 -> -752851 ->-751620.2 -750811.1 -> -750786.7 ->-750351.4 > plot(LCF$x,LCF$beta1,type="l",col="red") > lines(LCH$x,LCH$beta1,col="blue") > plot(LCF$x,LCF$beta2,type="l",col="red") > lines(LCH$x,LCH$beta2,col="blue") > plot(LCF$y,LCF$kappa2,type="l",col="red") > lines(LCH$y,LCH$kappa2,col="blue") > MUF = DTF/ETF > MUH = DTH/ETH > MUFMODEL = MUF; MUHMODEL = MUH > for(i in 1:nrow(MUF)){ + for(j in 1:ncol(MUF)){ + MUFMODEL[i,j] = exp(LCF$beta1[j]
+ +LCF$beta2[j]*LCF$kappa2[i])
+ MUHMODEL[i,j] = exp(LCH$beta1[j] + +LCH$beta2[j]*LCH$kappa2[i]) + }} soit, pour avoir un dessin, > persp(LCH$y,LCH$x,log(MUH), + theta=45,col="light green",shade=TRUE) > persp(LCH$y,LCH\$x,log(MUHMODEL),