# Freakonometrics

## Tag - demography

Saturday, June 23 2012

## Actuarial models with R, Meielisalp

I will be giving a short course in Switzerland next week, 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 long version of the slides for the course on Actuarial models with R can be found online with the #Rmetrics tag, and the short version will be uploaded soon. There will be some practicals, based on Swiss mortality table for the part on demography. The datasets can be uploaded using the following code,

```DEATH=read.table(
"http://freakonometrics.free.fr/DeathsSwitzerland.txt",
"http://freakonometrics.free.fr/ExposuresSwitzerland.txt",
DEATH\$Age=as.numeric(as.character(DEATH\$Age))
DEATH=DEATH[-which(is.na(DEATH\$Age)),]
EXPOSURE\$Age=as.numeric(as.character(EXPOSURE\$Age))
EXPOSURE=EXPOSURE[-which(is.na(EXPOSURE\$Age)),]```
• based on those datasets, plot the log mortality rates for male and female,

• for those two datasets, plot the log mortality rates in 1900 and 1950, respectively
• for those two datasets, plot the log mortality rates for the cohort born on 1900 and 1950, respectively
• on the total dataset (male and female), fit a Lee-Carter model. Plot the age coefficients

• plot the time coefficients and propose a forecast for that series of estimators.

• plot the residuals obtained from the regression

• using those estimates, and the forecasts, project the log-mortality rates

• extrapolate the survival function of an insured aged 40 in 2000, and compare it with the one obtained on the longitudinal dataset.

• based on those survival functions, compute actuarial present values for several quantities, e.g. whole life annuities for some insured aged 40, and whole life insurances, and compare those values from 1900 till 2040 (on the graphs below, titles were inverted).

Then, we will briefly mention payment triangles. We will work on the triangle used on http://rworkingparty.wikidot.com/ that can be downloaded below,
```OthLiabData=read.csv(
"http://www.casact.org/research/reserve_data/othliab_pos.csv",
library(plyr)
OL=SumData=ddply(OthLiabData,.(AccidentYear,DevelopmentYear,
DevelopmentLag),summarise,IncurLoss=sum(IncurLoss_H1-BulkLoss_H1),
CumPaidLoss=sum(CumPaidLoss_H1), EarnedPremDIR=
sum(EarnedPremDIR_H1))
LossTri <- as.triangle(OL, origin="AccidentYear",
dev="DevelopmentLag", value="IncurLoss")
Year=as.triangle(OL, origin="AccidentYear",
dev="DevelopmentLag", value="DevelopmentYear")
TRIANGLE=LossTri
TRIANGLE[Year>1997]=NA```
Here, the triangle looks like that
```> TRIANGLE
dev
origin      1      2      3      4      5      6      7      8      9     10
1988 128747 195938 241180 283447 297402 308815 314126 317027 319135 319559
1989 135147 208767 270979 304488 330066 339871 344742 347800 353245     NA
1990 152400 238665 297495 348826 359413 364865 372436 372163     NA     NA
1991 151812 266245 357430 400405 423172 442329 460713     NA     NA     NA
1992 163737 269170 347469 381251 424810 451221     NA     NA     NA     NA
1993 187756 358573 431410 476674 504667     NA     NA     NA     NA     NA
1994 210590 351270 486947 581599     NA     NA     NA     NA     NA     NA
1995 213141 351363 444272     NA     NA     NA     NA     NA     NA     NA
1996 237162 378987     NA     NA     NA     NA     NA     NA     NA     NA
1997 220509     NA     NA     NA     NA     NA     NA     NA     NA     NA```
• suggest an estimation for the amount of reserves, all years.
• using a Poisson regression, propose a VaR with level 99.5% for future payments, for all claims that already occurred.

Monday, June 4 2012

## Longevity and mortality dynamics with R

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",
+ "http://freakonometrics.free.fr/Exposures-France.txt",

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")```

Sunday, December 18 2011

## 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ées inté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).

Tuesday, December 7 2010

## Statistique de l'assurance STT6705V, partie 12 bis

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      BIC941.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)
`> ETSETS(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      BIC937.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: YARIMA(0,1,0) with drift Call: Arima(x = Y, order = c(0, 1, 0), include.drift = TRUE) Coefficients:drift-2.0337s.e.   1.1904 sigma^2 estimated as 138.9:  log likelihood = -380.8AIC = 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: YARIMA(0,1,1) with drift Call: auto.arima(x = Y, allowdrift = TRUE) Coefficients:          ma1    drift      -0.3868  -2.0139s.e.   0.0970   0.6894 sigma^2 estimated as 122.3:  log likelihood = -374.65AIC = 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=2000x=40s=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=.035m=70h=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,

Thursday, December 2 2010

## Statistique de l'assurance STT6705V, partie 12

The final course (since courses end this week in Montréal) can be watched here and there. The drawings from the course can be downloaded here (including last week's). First, to come back on last week's course , we considered Lee-carter model, i.e.

Note that it is possible to go a bit further, and to consider something that can be seen as a second order development

or even more generally

This idea was implicitly in the initial paper published in the 80's. Hence, I mentioned in the course that (weighted) least square techniques can be considered to estimate parameters,

and the the first order condition is simply

(which was the interpretation that it is simply the average (over years) of mortality rate). Given that estimator, it is possible to write

and then to use the first components of the singular value decomposition to derive estimators.About the singular value decomposition (SVD, with a nice graph stolen from wikipedia, here, on the right - actually the French article on SVD is a bit more punchy, from Снайперская винтовка Драгунова, there). Note that SVD is closely related to PCA...
Thus, instead of focusing only on the first principal components, it is possible to go up to order 2.
But it is also possible to add, instead of a time component, a cohort-based component, i.e.

The estimation is rather simple, as shown below
`D=as.vector(BASEB)E=as.vector(BASEC)A=rep(AGE,each=length(ANNEE))Y=rep(ANNEE,length(AGE))C=Y-Abase=data.frame(D,E,A,Y,C,a=as.factor(A),y=as.factor(Y),c=as.factor(C))LC4=gnm(D~a+Mult(a,y)+Mult(a,c),offset=log(E),family=poisson,data=base)`
Here, parameters have the following shape.
`plot(AGE[-1],LC2\$coefficients[2:90])plot(AGE,LC2\$coefficients[91:180])plot(ANNEE,LC2\$coefficients[181:279])plot(AGE,LC2\$coefficients[280:369])plot(1808:1997,LC2\$coefficients[370:559])`
The first component is rather similar to what we had before
On the other hand, the Lee-Carter (age-year) component is and the cohort effect (age-cohort) is Note that the year-component captures wars and chocks that can affect anyone some given years, and that gains on life expectancy is more a cohort effect.
Finally, we discussed actuarial applications. But first, recall that in actuarial literature (without any dynamics), it is standard to defined

and

The mortality rate is then

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)`
while it is
`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\$ktK2=LC3\$kt[99]+LC3F\$kt.f\$meanK=-c(K1,K2)/50MU=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])}}`
Once we have that matrix, we simply have to work on diagonal (i.e. cohorts) to calculate anything. For instance, to get the dynamic version of , the code is (here for someone of age 40 in 2000)
`t=2000x=40s=seq(0,99-x-1)MUd=MU[x+1+s,t+s-1898]Pxt=cumprod(exp(-diag(MUd)))`
It is also possible to compute residual life expectancy , including dynamics
`x=40E=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)`
Note that this has been computed given the maximum lifetime we had in the dataset (here 99 years old). Hence, we assume here that no one can live more than 100 year (which will impact life expectancy and is a rather strong assumption, especially in 2050 - which might explain the concave shape of the curve on the right part).
It is also possible to compute annuities. For instance, consider a 40 years old insured, in 2000, willing to get 1 every year after age 70, until he dies, i.e.
The annuity for such a contract is
`x=40t=2000r=.035m=70h=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)`
(in order to compare the value of such a contract, from 1900 up to 2050).

Monday, November 29 2010

## Statistique de l'assurance STT6705V, partie 11

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,

`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=as.matrix(BASEB[,1:100])AGE=0:ncol(BASEB)tabC=read.table("http://perso.univ-rennes1.fr/arthur.charpentier/tabC.csv", sep=";",header=FALSE) an=tabC[,1]an=an[-c(16,23,43,48,51,53,106)]BASEC=tabC[,2:101] BASEC=as.matrix(BASEC[-c(16,23,43,48,51,53,106),])  BASEB=BASEB[,1:90]BASEC=BASEC[,1:90]AGE=AGE[1:90]MU=as.matrix(log(BASEB/BASEC))persp(ANNEE,AGE,MU,theta=70,shade=TRUE,col="green")library(rgl)persp3d(ANNEE,AGE,MU,col="light blue")`
(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,
`source("http://perso.univ-rennes1.fr/arthur.charpentier/fitModels.r")x=AGEy=ANNEEd=BASEBe=BASECw=matrix(1,nrow(d),ncol(e))LC1=fit701(x,y,e,d,w)plot(AGE,LC1\$beta1)plot(ANNEE,LC1\$kappa2)plot(AGE,LC1\$beta2)`
Then we considered nonlinear Poisson regression,
`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)plot(AGE[-1],LC2\$coefficients[1]+LC2\$coefficients[2:90])lines(AGE,LC1\$beta1,col="blue")plot(ANNEE,LC2\$coefficients[181:279])plot(ANNEE,LC1\$kappa2,col="red")plot(AGE,LC1\$beta2)`
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,
`plot(AGE[-1],LC2\$coefficients[1]+LC2\$coefficients[2:90],type="l",col="blue",xlab="",ylab="")for(s in 1:200){LC2=gnm(D~a+Mult(a,y),offset=log(E),family=poisson,data=base)lines(AGE[-1],LC2\$coefficients[1]+LC2\$coefficients[2:90],col="blue")}`
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,
`library(demography)base=demogdata(data=t(exp(MU)),pop=t(BASEC),ages=AGE,years=ANNEE,type="mortality",label="France",name="Total",lambda=0)LC3=lca(base)LC3F=forecast(LC3,100)plot(LC3\$year,LC3\$kt,xlim=c(1900,2100),ylim=c(-300,150))lines(LC3F\$year,LC3\$kt[99]+LC3F\$kt.f\$mean)lines(LC3F\$year,LC3\$kt[99]+LC3F\$kt.f\$lower,lty=2)lines(LC3F\$year,LC3\$kt[99]+LC3F\$kt.f\$upper,lty=2)`
We concluded with a short discussion about errors of the Lee-Carter model (on French mortality)
`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))RES=residuals(LC2,"pearson")base\$res=RESplot(base\$A,base\$res)couleur=heat.colors(100)plot(base\$A,base\$res,col=couleur[base\$Y-1898])plot(base\$Y,base\$res,col=couleur[base\$A+1])`
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...

Tuesday, November 23 2010

## Statistique de l'assurance STT6705V, partie 11

Pour récupérer les données, et définir un taux de mortalité en fonction de l'année et de l'âge, le code ressemble à ça,

`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=as.matrix(BASEB[,1:100])tabC=read.table("http://perso.univ-rennes1.fr/arthur.charpentier/tabC.csv", sep=";",header=FALSE) an=tabC[,1]an=an[-c(16,23,43,48,51,53,106)]BASEC=tabC[,2:101] BASEC=as.matrix(BASEC[-c(16,23,43,48,51,53,106),]) AGE=0:(ncol(BASEC)-1)MU=as.matrix(log(BASEB/BASEC))persp(ANNEE,AGE,MU,theta=30)`

Pour la modélisation des tables de mortalité, nous allons utiliser,

• le package demography de Rob Hyndman (ici) qui propose par exemple d'implémenter le modèle de Lee & Carter.
• les fonctions de LifeMetrics de JP Morgan (ici), en particulier dans le code source fitModel.r dont le détail technique traîne en ligne.
• les fonctions du package gnm de modèles Poissonniens non-linéaires.
Les codes R sont les suivants,
`library(demography)BASEF = demogdata(data=MUF, pop=POPF,        ages=AGE, years=YEAR, type="mortality",        label="France", name="Femmes", lambda=1)LCF1 = lca(BASEF)source("http://perso.univ-rennes1.fr/arthur.charpentier/fitModels.r")LCF2 = fit701(xv=XV,yv=YV,etx=ETF,dtx=DTF,wa=WA)library(gnm)LCF3 = gnm(D~factor(Age)+       Mult(Exp(factor(Age)),factor(Year)),       data=base,offset=log(E),family=quasipoisson)`

Tuesday, September 14 2010

## Espérance de vie à la naissance, ou à 65 ans ?

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 , 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.
> ANNEE=tabB[,1]
> BASEB=tabB[,seq(2,246,by=2)]
> BASEB=BASEB[,1:100]
> 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 (), 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
Le chargement a nécessité le package : mvtnorm
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.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...

Thursday, January 29 2009

## Actuariat, démographie et R

Un certain nombre de fonctions sous R permettent de faire de la projection de tables de mortalité.

• le package demography de Rob Hyndman (ici) qui propose par exemple d'implémenter le modèle de Lee & Carter.
• les fonctions de LifeMetrics de JP Morgan (ici), en particulier dans le code source fitModel.r dont le détail technique traîne en ligne.
• plus récemment, l'Institute of Actuaries a lancé le CMI project:"The CMI’s Mortality Projections Working Party has evaluated the P-Spline and Lee-Carter methodologies for projecting mortality. [...] The CMI would like to express its thanks to James Kirkby and Dr Iain Currie of Heriot Watt University for developing this software. The software requires the statistical package R to be installed together with the R(D)COM interface which allows Excel to call R. Please note that the software will not function on Excel 97 or earlier versions as these do not contain some of the VBA commands required by the software. In order to obtain the software from the CMI, please complete the CMI Software Request Form."
à suivre donc...

[04/05/09] Comme cela m'a été fait remarqué, le lien vers le CMI project ne fonctionne plus. Je mets en ligne ici le fichier qu'il était alors possible de télécharger, et que quelqu'un m'a gentiment envoyé.