Freakonometrics

Following the previous post on life contingencies and actuarial models in life insurance, I upload additional material for the short course at the 6th R/Rmetrics Meielisalp Workshop & Summer School on Computational Finance and Financial Engineering organized by ETH Zürich, https://www.rmetrics.org/. The second part of the talk (on Actuarial models with R) will be dedicated to longevity and mortality. A complete set of slides can be downloaded from the blog, but again, only some part will be presented.

As mentioned earlier, the codes are from a book on actuarial science in R, written with Christophe Dutang and Vincent Goulet (so far in French) that should appear, some day... The code used in the slides above can be downloaded from here, and datasets are the following,

> DEATH <- read.table(
+ "http://freakonometrics.free.fr/Deces-France.txt",
+ header=TRUE)
> EXPO  <- read.table(
+ "http://freakonometrics.free.fr/Exposures-France.txt",
+ header=TRUE,skip=2)

For additional resources, I will use Rob Hyndman's package on demography, Heather Turner and David Firth's package on generalized nonlinear models (e.g. the slides of the short course Heather gave in Rennes at the UseR! conference in 2009), as well as functions developed by JPMorgan's LifeMetrics (functions are  fully documented in the LifeMetrics Technical Document). All those functions can be obtained using

> library(demography)
> library(gnm)
> source("http://freakonometrics.free.fr/fitModels.R")

In the previous post (here) discussing forecasts of actuarial quantities, I did not mention much how to forecast the temporal component in the Lee-Carter model. Actually, many things can be done. Consider here some exponential smoothing techniques (see here for more details). For instance, with a simple model, with an additive error, no trend and no seasonal component,

or equivalently

Then , and if we want a confidence interval

> (ETS.ANN=ets(Y,model="ANN"))
ETS(A,N,N)
 
Call:
ets(y = Y, model = "ANN")
 
Smoothing parameters:
alpha = 0.7039
 
Initial states:
l = 75.0358
 
sigma:  11.4136
 
AIC     AICc      BIC
941.0089 941.1339 946.1991
> plot(forecast(ETS.ANN,h=100))

Here, a natural idea will probably be to take into account a linear trend on the series, i.e.

where

then the forecast we make is  i.e.

ETS.AAN=ets(Y,model="AAN")
plot(forecast(ETS.AAN,h=100))

It is also possible to ask for an "optimal" exponential smoothing model (without any specification)

> ETS
ETS(A,A,N)
 
Call:
 ets(y = Y)
 
  Smoothing parameters:
    alpha = 0.6107
    beta  = 1e-04
 
  Initial states:
    l = 74.5622
    b = -1.9812
 
  sigma:  11.0094
 
     AIC     AICc      BIC
937.8695 938.2950 948.2500
> plot(forecast(ETS,h=100))

Another idea can be to consider some autoregressive integrated moving average models (ARIMA), here with a linear trend

For instance, if we want only the integrated component, i.e.

then the code is

> ARIMA010T=Arima(Y,order=c(0,1,0),include.drift=TRUE)
Series: Y
ARIMA(0,1,0) with drift
 
Call: Arima(x = Y, order = c(0, 1, 0), include.drift = TRUE)
 
Coefficients:
drift
-2.0337
s.e.   1.1904
 
sigma^2 estimated as 138.9:  log likelihood = -380.8
AIC = 765.6   AICc = 765.73   BIC = 770.77
> plot(forecast(ARIMA010T,h=100))

But note that any kind of ARIMA model can be considered, e.g. integrated with an autoregressive component (here of order 1)
or with also a moving average component (again of order 1)
Actually, note that, once again, we can ask for an automatic selection of the model,

> (autoARIMA=auto.arima(Y,allowdrift=TRUE))
Series: Y
ARIMA(0,1,1) with drift
 
Call: auto.arima(x = Y, allowdrift = TRUE)
 
Coefficients:
          ma1    drift
      -0.3868  -2.0139
s.e.   0.0970   0.6894
 
sigma^2 estimated as 122.3:  log likelihood = -374.65
AIC = 755.29   AICc = 755.55   BIC = 763.05
> plot(forecast(autoARIMA,h=100))

Finally, it is possible to use also also some specific functions of the package, for instance to consider a random walk with a drift,

RWF=rwf(Y,h=100,drift=TRUE)
plot(RWF)

or Holt model (with a trend)

HOLT=holt(Y,h=100,damped=TRUE)
plot(HOLT)

And if all that is not enough, it is also possible to go further by changing the size of the series we use to fit the model. A question that naturally arises is the treatment of wars in our model: shouldn't we assume that the forecast should be based only on the last 50 years (and exclude wars from our analysis) ? In that case, for instance, the exponential smoothing technique gives
while the ARIMA procedure returns
And finally, with the Holt technique, we have
So, it looks like we have a lot of techniques that can be used to provide a forecast for the temporal component in the Lee-Carter model,
All those estimators can be used to estimate annuities of insurance contracts (as here),

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)
 
A=LC2$coefficients[1]+LC2$coefficients[2:99]
B=LC2$coefficients[100:198]
K0=LC2$coefficients[199:297]
Y=as.numeric(K01)
K1=c(K0,forecast(ets(Y),h=120)$mean)
K2=c(K0,forecast(auto.arima(Y,allowdrift=TRUE),h=120)$mean)
K3=c(K0,rwf(Y,h=120,drift=TRUE)$mean)
K4=c(K0,holt(Y,h=120,drift=TRUE)$mean)
K5=c(K0,forecast(ets(Y[50:99]),h=120)$mean)
K6=c(K0,forecast(auto.arima(Y[50:99],allowdrift=TRUE),h=120)$mean)
K7=c(K0,rwf(Y[50:99],h=120,drift=TRUE)$mean)
K8=c(K0,holt(Y[50:99],h=120,drift=TRUE)$mean)
 
MU=matrix(NA,length(A),length(K1))
MU1=MU2=MU3=MU4=MU5=MU6=MU7=MU8=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]*K5[j])
MU3[i,j]=exp(A[i]+B[i]*K3[j])
MU4[i,j]=exp(A[i]+B[i]*K4[j])
MU5[i,j]=exp(A[i]+B[i]*K5[j])
MU6[i,j]=exp(A[i]+B[i]*K6[j])
MU7[i,j]=exp(A[i]+B[i]*K7[j])
MU8[i,j]=exp(A[i]+B[i]*K8[j])
}}
 
t=2000
x=40
s=seq(0,98-x-1)
 
Pxt1=cumprod(exp(-diag(MU1[x+1+s,t+s-1898])))
Pxt2=cumprod(exp(-diag(MU2[x+1+s,t+s-1898])))
Pxt3=cumprod(exp(-diag(MU3[x+1+s,t+s-1898])))
Pxt4=cumprod(exp(-diag(MU4[x+1+s,t+s-1898])))
Pxt5=cumprod(exp(-diag(MU5[x+1+s,t+s-1898])))
Pxt6=cumprod(exp(-diag(MU6[x+1+s,t+s-1898])))
Pxt7=cumprod(exp(-diag(MU7[x+1+s,t+s-1898])))
Pxt8=cumprod(exp(-diag(MU8[x+1+s,t+s-1898])))
 
r=.035
m=70
h=seq(0,21)
V1=1/(1+r)^(m-x+h)*Pxt1[m-x+h]
V2=1/(1+r)^(m-x+h)*Pxt2[m-x+h]
V3=1/(1+r)^(m-x+h)*Pxt3[m-x+h]
V4=1/(1+r)^(m-x+h)*Pxt4[m-x+h]
V5=1/(1+r)^(m-x+h)*Pxt5[m-x+h]
V6=1/(1+r)^(m-x+h)*Pxt6[m-x+h]
V7=1/(1+r)^(m-x+h)*Pxt7[m-x+h]
V8=1/(1+r)^(m-x+h)*Pxt8[m-x+h]

Hence, here the difference is significant,

> M=cbind(V1,V2,V3,V4,V5,V6,V7,V8)
> apply(M,2,sum)
V1       V2       V3       V4       V5       V6       V7       V8
4.389372 4.632793 4.406465 4.389372 4.632793 4.468934 4.482064 4.632793

or graphically,

Last course will be uploaded soon (the links will be here and there). The R code considered is given below. First, we had to work a little bit on the datasets,

(this last line is here to play a little bit with the 3d mortality surface). We first used the Lee-Carter function proposed by JPMorgan in LifeMetrics,
Then we considered nonlinear Poisson regression,
As mentioned during the course, this technique is great... but it is sentive to initial values in the optimization procedure. For instance, consider the following loops,
Here are representation of the first component in the Lee-Carter model,
Hence, the estimation depends on the choice of initial values... even if the shape remains unchanged...
Then, we finnished with Rob Hyndman's package,
We concluded with a short discussion about errors of the Lee-Carter model (on French mortality)
The graphs can be seen below (as a function of time, and a function of ages)
his Wednesday, we will discuss how we can use those estimators (or other ones) in life insurance...

Bon, je vais faire un rapide billet pour répondre à un commentaire qui avait été posté ici, qui remettait en cause mes conclusions (ce qui ne me gêne pas, loin de là, c'est comme ça qu'on va faire avancer les choses), et surtout demandait des informations complémentaires. J'avais un peu forcé le trait en affirmant que "une personne de 70 ans vivra à peine plus longtemps - en moyenne - qu’une personne de 70 ans en 1950" (c'est peut être le cas pour les centenaires, mais je n'ai pas de données), mais le but de mon précédant billet était de montrer que les gains sont beaucoup plus faibles que pour les plus jeunes.
J'ai décidé de refaire l'étude, de manière plus claire peut être, et plus complète. Et en utilisant des données de l'INED, qui est une source a priori fiable de données.

• Les données (détaillées) 

Comme je le disais, tout d'abord, les données: je prends ici les données de l'INED (et pas de http://mortalityorg). On peut trouver en ligne (au format xls) des données de population par âge et par sexe, ici (populations par âge, de 0 à 100 ans, au premier janvier, de 1899 à 1998, selon le territoire couvert par la statistique des décès, ensemble des deux sexes, avec depuis 1946, les données INSEE et avant 1946, des reconstitutions intercensitaires (Vallin, 1973)). Et là, on a les données du nombre de décès par âge et par sexe (décès par âge et par génération, de 1899 à 1997, ensemble des deux sexes, avec pour les moins de 100 ans, depuis 1907, décès hors pertes militaires (SGF, puis INSEE), pour les années de guerre : décès (INSEE + estimations des pertes militaires (Vallin, 1973)) et entre 1899 et 1906 : décès par âge de l'INSEE, répartition entre les triangles de Lexis (Vallin, 1973). Pour les centenaires, jusqu'en 1967 : SGF, INSEE + répartitions estimées par les auteurs, et depuis 1968 : INSEE (extrait d'enregistrements individuels fournis par l'INSEE dans le cadre d'un avenant à la convention INED-INSEE)). Voilà pour les données, ou presque....
Par paresse, et pour pouvoir partager les données plus facilement, j'ai mis en ligne des fichiers csv sur ma page.
Les données (brutes) contiennent des colonnes sans intérêt, que je me suis permis de virer. Sinon dans la base des populations, certaines lignes ont été supprimées car en 1914 il y avait deux populations, sur territoire d'avant 1914 (sans l'Alsace-Lorraine: Moselle, Bas-Rhin, Haut-Rhin), et le territoire sans opérations militaires (manquent, outre l'Alsace-Lorraine: Aisne, Ardennes, Marne, Meurthe-et-Moselle, Meuse, Nord, Oise, Pas-de-Calais, Somme et Vosges). C'est cette dernière qui est retenues dans la base des décès c'est donc cette dernière que l'on utilisera pour calculer l'exposition.
> tabB=read.table("http://perso.univ-rennes1.fr/arthur.charpentier/tabB.csv",
+      sep=";",header=FALSE)
> ANNEE=tabB[,1]
> BASEB=tabB[,seq(2,246,by=2)]
> BASEB=BASEB[,1:100]
> tabC=read.table("http://perso.univ-rennes1.fr/arthur.charpentier/tabC.csv",
+      sep=";",header=FALSE)
> BASEC=tabC[,2:101]
> BASEC=BASEC[-c(16,23,43,48,51,53),] 
> BASEC=BASEC[1:nrow(BASEB),]
> AGE=0:99

• Les taux de décès (population entière)

Commençons par calculer, et tracer les taux de décès, qui sont les ratios entre les nombres de décès et la population totale.
Mon premier point (pour répondre au commentaire) est que si l'on regarde un peu la surface, les gains en terme de taux de mortalité instantanés sont plus importants à la naissance (courbe rouge) qu'à des âges beaucoup plus élevés, comme 90 ans (en bleu) voire à 65 ans (en vert),

• Un peu de prospective

Comme on le voit, sur nos données, on est limité car on n'a qu'un siècle de données, et c'est donc délicat de suivre les espérances de vie.
Le plus simple est alors de reprendre ce que j'avais fait (en ligne ici), utilisons les fonctions de Rob Hyndman, en ligne sur son blog (là), mais sur nos données (la dernière fois j'avais utilisé les siennes).
>   library(demography)
Le chargement a nécessité le package : forecast
Le chargement a nécessité le package : fracdiff
This is forecast 2.05
Le chargement a nécessité le package : rainbow
Le chargement a nécessité le package : hdrcde
Le chargement a nécessité le package : locfit
Le chargement a nécessité le package : akima
Le chargement a nécessité le package : lattice
locfit 1.5-6     2010-01-20
Le chargement a nécessité le package : ash
Le chargement a nécessité le package : ks
Le chargement a nécessité le package : KernSmooth
KernSmooth 2.23 loaded
Copyright M. P. Wand 1997-2009
Le chargement a nécessité le package : mvtnorm
hdrcde 2.13 loaded
Attachement du package : 'hdrcde'
Le chargement a nécessité le package : MASS
Le chargement a nécessité le package : pcaPP
Le chargement a nécessité le package : cluster
Le chargement a nécessité le package : ftsa
This is demography 1.03
Il y a eu 13 avis (utilisez warnings() pour les visionner)
Commençons par mettre les données au bon format. Pour des raisons numériques, je suis obligé de me débarrasser des âges élevés pour lesquels j'ai pu avoir une population nulle, ou pour lesquels j'ai des valeurs manquantes. J'ai donc exclu les âges supérieurs à 90 ans,
> donnees = demogdata(data=t(as.matrix(BASEB[,1:90]))/t(as.matrix(BASEC[,1:90])),
+                     pop= t(as.matrix(BASEC[,1:90])),
+                     ages=AGE[1:90],
+                     years=ANNEE,
+                     type="mortality", label="France", name="total", lambda=0)
> donnees 
Mortality data for France
    Series: total
    Years: 1899 - 1997
    Ages:  0 - 89
Ensuite, on peut reprendre le code que j'avais fait. Notons qu'on utiliser ici une modélisation à la Lee et Carter, ce qui a été utilisée par l'INED pour extrapoler certaines données dans la table.
>   france.LC2 <- lca(donnees,adjust="dt",series="total")
>   france.fcast <- forecast(france.LC2)
>   L2 <- lifetable(france.fcast)
>   ex2=L2$ex
>   L1=lifetable(donnees,series="total")
>   ex1=L1$ex
Pour l'espérance de vie à la naissance, et son évolution dans le temps, si on regarde les gains moyens année après année, on voit que l'on gagne, en moyenne, 1/4 d'espérance de vie,  soit 3 mois par an (comme dans mon précédant billet),
>   age=0
>   ex=c(ex1[age+1,],ex2[age+1,])
>   plot(1899:2047,ex,col="blue")
>   I=(1950:2000)-1898
>   y=ex[I]
>   x=1950:2000
>   lm(y~x)
Call:
lm(formula = y ~ x)
Coefficients:
(Intercept)            x 
  -421.0384       0.2549 
>   abline(lm(y~x),col="red")
>   points(x,y,pch=19,col="red")

On notera qu'en utilisant les autres méthodes d'estimation, j'obtiens la même pente. En revanche, si je regarde l'espérance de vie résiduelle à 65 ans, elle n'est plus que de 0,16 années, soit moins de deux mois chaque années,
>   age=65
>   ex=c(ex1[age+1,],ex2[age+1,])
>   plot(1899:2047,ex,col="blue")
>   I=(1950:2000)-1898
>   y=ex[I]
>   x=1950:2000
>   lm(y~x)
Call:
lm(formula = y ~ x)
Coefficients:
(Intercept)            x 
  -295.6812       0.1612 
>   abline(lm(y~x),col="red")
>   points(x,y,pch=19,col="red")

Dit autrement, entre un bébé qui nait en 1975 et un bébé qui nait en 1995, il y a un gain d'environ 6 ans d'espérance de vie, passant de 81,8 ans à 87,7, alors que pour une personne qui atteignait les 65 ans en 1975, il avait 22 ans à vivre encore, en moyenne, contre seulement 26,5 pour une personne de 65 ans en 1995. Donc je peux maintenir ce que j'affirmais la dernière fois, sur d'autres données.... En fait, si on veut être plus précis, on peut précisément regarder l'évolution de la pente, correspondant au gain moyen d'espérance de vie chaque année,

Bref, de considérables progrès ont été effectués en terme de gains d'espérance de vie, mais beaucoup plus sur les jeunes que sur les personnes plus âgées...