Category Archives: Demography

Job for life ? Bishop of Rome ?

The job of Bishop of Rome – i.e. the Pope – is considered to be a life-long commitment. I mean, it usually was. There have been 266 popes since 32 A.D. (according to http://oce.catholic.com/…): almost all popes have served until their death. But that does not mean that they were in the job for long… One can easily extract the data from the website,

> L2=scan("http://oce.catholic.com/index.php?title=List_of_Popes",what="character")
Read 4485 items
> index=which(L2=="</td><td>Reigned")
> X=L2[index+1]
> Y=strsplit(X,split="-")

But one should work a little bit because sometimes, there are inconsistencies, e.g. 911-913 and then 913-14, so we need some more lines. Further, we can extract from this file the years popes started to reign, the year it ended, and the length, using those functions

> diffyears=function(x){
+ s=NA
+ if(sum(substr(x,1,1)=="c")>0){x[substr(x,1,1)=="c"]=substr(x[substr(x,1,1)=="c"],3,nchar(x[substr(x,1,1)=="c"]))}
+ if(length(x)==1){s=1}
+ if(length(x)==2){s=diff(as.numeric(x))}
+ return(s)}
> whichyearsbeg=function(x){
+ s=NA
+ if(sum(substr(x,1,1)=="c")>0){x[substr(x,1,1)=="c"]=substr(x[substr(x,1,1)=="c"],3,nchar(x[substr(x,1,1)=="c"]))}
+ if(length(x)==1){s=as.numeric(x)}
+ if(length(x)==2){s=as.numeric(x)[1]}
+ return(s)}
> whichyearsend=function(x){
+ s=NA
+ if(sum(substr(x,1,1)=="c")>0){x[substr(x,1,1)=="c"]=substr(x[substr(x,1,1)=="c"],3,nchar(x[substr(x,1,1)=="c"]))}
+ if(length(x)==1){s=as.numeric(x)}
+ if(length(x)==2){s=as.numeric(x)[2]}
+ return(s)}

On our file, we have

> Years=unlist(lapply(Y,whichyearsbeg))
> YearsB=c(Years[1:91],752,Years[92:length(Years)])
> YearsB[187]=1276
> Years=unlist(lapply(Y,whichyearsend))
> YearsE=c(Years[1:91],752,Years[92:length(Years)])
> YearsE[187]=1276
> YearsE[266]=2013
> YearsE[122]=914 
> W=unlist(lapply(Y,diffyears))
> W=c(W[1:91],1,W[92:length(W)])
> W[W==-899]=1
> which(is.na(W))
[1] 187 266
> W[187]=1
> W[266]=2013-2005

If we plot it, we have the following graph,

> plot(YearsB,W,type="h")

and if we look at the average length, we have the following graph,

> n=200
> YEARS = seq(0,2000,length=n)
> Z=rep(NA,n)
> for(i in 2:(n-1)){
+ index=which((YearsB>YEARS[i]-50)&(YearsE<YEARS[i]+50))
+ Z[i] = mean(W[index])}
> plot(YEARS,Z,type="l",ylim=c(0,30))
> n=50
> YEARS = seq(0,2000,length=n)
> Z=rep(NA,n)
> for(i in 2:(n-1)){
+ index=which((YearsB>YEARS[i]-50)&(YearsE<YEARS[i]+50))
+ Z[i] = mean(W[index])}
> lines(YEARS,Z,type="l",col="grey")

which does not reflect mortality improvements that have been observed over two millenniums. It might related to the fact that the average age at time of election has  increased over time (for instance, Benedict XVI was elected at 78 – one of the oldest to be elected). Actually, serving a bit more than 7 years is almost the median,

> mean(W>=7.5)
[1] 0.424812

(42% of the Popes did stay at least 7 years in charge) or we can look at the histogram,

> hist(W,breaks=0:35)

Unfortunately, I could not find more detailed database (including the years of birth for instance) to start a life-table of Popes.

 

Les députés sont-ils à l’image de la population

Beaucoup de choses ont été écrites sur le fait que les députés ne sont pas vraiment le reflet de la population, que ce soit en terme de profession, de sexe, d’origine, d’age, etc. La liste pourrait être longue. Il y a plusieurs mois, j’avais commencé à regarder le profil des députés, par age. En effet, le site http://assemblee-nationale.fr/ permet d’accéder à des données sur tous les députés, depuis la Révolution. Y compris leur date de naissance. En croisant ces données avec des données de population, par exemple via http://www.mortality.org/, on peut comparer la répartition des ages des députés, avec la répartition des ages de la population.

Pour les amateurs, le code pour récupérer les données (ou au moins les dates de naissance des députés) ressemble à

N=2002
URL=paste("http://www.assemblee-nationale.fr/
sycomore/result.asp?radio_dept=tous_departements&
regle_nom=est&Nom=&departement=&
choixdate=intervalle&D%C3%A9butMin=01%2F01%2F",
N,"&FinMin=31%2F12%2F",N,"&Dateau=&legislature=",
s,"&choixordre=chrono&Rechercher=
Lancer+la+recherche",sep="")
HTML=scan(URL,what="character")

k=which(HTML=="class=\"titre\">Né")
vHTML=HTML[k:length(HTML)]
vk=which(substr(vHTML,1,7)==">&nbsp;")
liste=vHTML[vk]
naissance=liste[seq(1,length(liste),by=2)]
NAISSANCE=as.Date(substr(naissance,8,17),
"%d/%m/%Y")

Maintenant, pour être tout à fait honnête, je ne suis pas certain de ce qui est vraiment renvoyé, et j’ai des doutes que cela correspondent réellement à la requête faite. En effet, même si je demande à avoir la liste des députés après l’élection, j’ai trop de monde… mais peut-être est-ce du aux décès éventuels, et il est possible que l’ensemble des députés qui ont siégé pendant la mandature apparaissent dans le résultat de la requête.

Sur la figure suivante on voit, sur plusieurs élections depuis plus de 100 ans, comment les deux distributions se déforment, avec en rouge la distribution de l’age des députés, et en bleu, la distribution de la population française, dans son ensemble (population de plus de 18 ans)

Si on veut tout suivre sur un graphique, au lieu de se regarder une animation, on peut représenter les différents quantiles (10%, 25%, 75% et 90%, retenus sur la population de plus de 18 ans, et l’age médian, au centre), avec la population française l’année de l’élection,

et l’ensemble des élus au parlement,

Si on veut faciliter la comparaison, on peut se contenter de visualiser l’évolution des ages moyens,

ou encore, du ratio (en % de différence) entre l’age moyen des députés, et celui de l’ensemble de la population.

Sur ce graphique, on voit que depuis 30 ans, l’age moyen des députés croit plus vite que celui de la population: la population français vieilli, mais moins que ses députés… La gérontocratie perdure donc en France. En espérant que cela ne débouche pas sur le clash générationnel que l’on semble observer ces temps-ci au Québec…

Date of death, birthday and Elvis Presley

10 days ago, a study published on http://www.annalsofepidemiology.org/ mentioned that “Death has a preference for birthdays” (as claimed in the title). The conclusion of the paper is that, in general, birthdays do not evoke a postponement mechanism but appear to end up in a lethal way more frequently than expected (“anniversary reaction”). Well, this is not new, and several previous articles have mentioned that point, e.g. Angermeyer et al. (1987).

I found the idea interesting since in demography, there is a large literature trying to extrapolate death rates from discrete to continuous time. Extrapolation are usually extremely smooth. But none of them integrate that aspect of mortality precisely on the birthday. The problem is that it is rather difficult to say something since datasets with individual observations are rare, online.

But yesterday, @coulmont sent me a tweet mentioning a website. I do not know if this is legal (even if some explanations are given), but I will mention courtesy of http://ssdmf.info/. It is a so-called Social Security Death Master File, containing individual informations about deaths in the US, as well as geographic information (as described on http://www.ssa.gov/), for people having a social security number.

With R, it is possible to work on those files (even they are huge, with tens of millions observations). For instance, we can check who is inside.

> elvis=scan("ssdm2",skip=22371720,n=1,what="character",sep=",")
> elvis
[1] " 409522002PRESLEY         ELVIS     0800197701081935  "

If you believe that Elvis is dead, you might agree that this database can be accurate (or at least, not too bad). And further, we can see here how to read the result: Elvis was born on January 8, 1935 (8 last digits), and died on August 16, 1977 (8 digits before). Obviously here, there are some problems with the dataset (we do not have the day of the death of Elvis). So here, we remove all the observations that do not give us proper dates. Then, the idea is to assume that the person died in 2000 (or any year since the point is to focus on days and months). Then, we count the number of days between the day of death and the birthday in 2001 (that would have been after) and the one in 2000 (that was either before or after the death), so that we can derive the number of days after the birthday,

dates=substr(base,66,81)
death=as.Date(substr(dates,1,8),"%m%d%Y")
birth=as.Date(substr(dates,9,16),"%m%d%Y")
indice=is.na(death)|is.na(birth)
mean(indice)
mdeath=substr(dates,1,2)
ddeath=substr(dates,3,4)
mbirth=substr(dates,9,10)
dbirth=substr(dates,11,12)
indice=which(ddeath!="00")
birth1=as.Date(paste(mbirth[indice],
dbirth[indice],"2000",sep=""),"%m%d%Y")
birth2=as.Date(paste(mbirth[indice],
dbirth[indice],"2001",sep=""),"%m%d%Y")
death=as.Date(paste(mdeath[indice],ddeath[indice],
"2000",sep=""),"%m%d%Y")
k=length(indice)
diffday=cbind((as.numeric(death-birth1))[1:k],
(as.numeric(death-birth2))[1:k])
DIFF=apply(diffday,1,function(x) {min(x[x>=0])})

What we have here is the number of days following the previous birthday. If we look at the distribution of that number of days, we obtain

counts=table(DIFF)
plot(as.numeric(names(counts)),
as.numeric(counts))
counts["0"]/(mean(counts[100:200]))
> counts["0"]/(mean(counts[100:200]))
0
1.121261

Thus, the death excess on the day of birth was around 12%, which is rather close to the one obtained from the Swiss mortality statistics 1969–2008 (in Ajdacic-Gross et al. (2012)). Note that here, we just play with a small subset of the entire dataset,

That database is probably extremely interesting, except that it suffers a huge selection bias, since only dead people are in that database. So it might be useless if we wish to study life expectancy of people named Bill versus people named Georges (that was something I wanted to investigate initially). But we’ll see what else we can do with it (since Ewen have been able to write some code to go through that huge dataset).

La pyramide des âges, au Canada

L’autre jour, Statistique Canada publiait une jolie animation interactive permettant de visualiser la déformation de la “pyramide” des âges (évoquée par exemple surhttp://www.lapresse.ca/). Bon, j’ai des données assez proches, tirées dehttp://www.mortality.org/).

> pop=read.table(
+ "http://freakonometrics.blog.free.fr/PopulationCanada.txt",
+ header=TRUE,skip=2)
> pop$Age=as.numeric(as.character(pop$Age))
> pop$Year=as.numeric(as.character(pop$Year))
> pop=pop[is.na(pop$Age)==FALSE,]

On peut faire assez facilement des pyramides des ages, par année, avec le code suivant,

> library(plotrix)
> seuils=seq(0,110,by=10)
> pop$tranche=cut(pop$Age,seuils, right = FALSE)
> an=1955
> base=pop[pop$Year==an,]
> women=aggregate(base$Female,
+ list(base$tranche),sum)[,2]
> men=aggregate(base$Female,
+ list(base$tranche),sum)[,2]
> nom=as.character(unique(pop$tranche))
> pyramid.plot(men/sum(men)*100,
+ women/sum(women)*100,labels=nom,gap=2,
+  lxcol=c("blue","blue","purple","purple","purple",
+ "purple","red","red","red","red","red"),
+  rxcol=c("blue","blue","purple","purple","purple",
+ "purple","red","red","red","red","red"))

http://freakonometrics.blog.free.fr/public/perso6/pyramide-ages-canada2.gif

Au lieu de faire des barres horizontales, on peut cumuler, par tranche d’age, et regarder l’évolution dans le temps des proportions “moins de 20 ans” (bleu), “entre 20 et 60 ans” (mauve) et “plus de 60 ans” (rouge).

> seuils=c(0,20,60,110)
> pop$tranche2=cut(pop$Age,seuils, right = FALSE)
> YEAR=1921:2008
> totaux=matrix(NA,length(YEAR),3)
> for(i in 1:length(YEAR)){
+ base=pop[pop$Year==YEAR[i],]
+ totaux[i,]=aggregate(base$Total,
+ list(base$tranche2),sum)[,2]
+ }
> stotaux=apply(totaux,1,sum)
> X=totaux[,1]/stotaux
> Y=X+totaux[,2]/stotaux

Si on regarde bien la pyramide, on s’aperçoit qu’après guerre, on a une génération importante, qui se déplace ensuite progressivement vers le haut, les baby-boomers. Pour visualiser encore davantage cette génération, on peut l’isoler sur le graphique ci-dessous, en suivant la cohorte née entre 1945 et 1950.

Mais on reviendra un peu plus en détails très bientôt sur cette génération (peut-être plutôt sur des données françaises cette fois afin de limiter les erreurs d’interprétation). A suivre…

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

Births and week-ends, in France

This week, I have seen on the internet (sorry, I cannot find proper references) the graph produced here on the right: which birthday is most likely ? The fact that I have no further information is important, since I do not know in which country such a graph was obtained. At least, I know it should not be France…

In France, I have already mentioned that there is a strong week-end effect: nowadays, there is 25% less deliveries during week-ends than during the week. Calot (1981) observed already that there were less deliveries on Sundays. This has been confirmed more recently, e.g. in http://www.lepoint.fr/ or http://www.prepabl.fr/, with a significant difference between week days, and week-ends. Here  is the number of birth per day, over 40 years, with in blue the average trend during the week, and in red, during week-ends,

naissance=read.table(
"http://freakonometrics.free.fr/naissanceFR2.txt")
attach(naissance)
date=as.Date(date)
plot(date, nbre,cex=.5)
t2=as.POSIXlt(date)
jour=t2$wday
X=naissance$date
Y=naissance$nbre
J=jour
df=data.frame(X,Y,J)
library(splines)
regs=lm(Y~bs(X,df=20),data=df[jour%in%c(0,6),])
Yp=predict(regs,newdata=df)
lines(X,Yp,col="red",lwd=3)
regs=lm(Y~bs(X,df=20),data=df[jour%in%1:5,])
Yp=predict(regs,newdata=df)
lines(X,Yp,col="blue",lwd=3)

If we look at the evolution of the ratio week-ends over weeks days, we have the following graph

t2=as.POSIXlt(date)
jour=t2$wday
jour=jour[1:(1982*7)]
nbre2=jour
for(i in 1:1982){
taux=sum(nbre[6:7+7*(i-1)])/
sum(nbre[1:5+7*(i-1)])/2*5
nbre2[1:5+7*(i-1)]=nbre[1:5+7*(i-1)]*taux
nbre2[6:7+7*(i-1)]=nbre[6:7+7*(i-1)]
nbre2[1:7+7*(i-1)]=
mean(nbre[1:7+7*(i-1)])/mean(nbre2[1:7+7*(i-1)])*
nbre2[1:7+7*(i-1)]
}
nbretaux=jour
for(i in 1:1982){
taux=sum(nbre[6:7+7*(i-1)])/
sum(nbre[1:5+7*(i-1)])/2*5
nbretaux[1:7+7*(i-1)]=taux
}
plot(date[1:length(nbre2)],nbretaux)
X= date[1:length(nbre2)]
Y=nbretaux
library(splines)
reg=lm(Y~bs(X,df=20))
Yp=predict(reg)
lines(X,Yp,col="red",lwd=3)

In the beginning of the 70’s, during week-ends, there were 5% less deliveries, but 25% less around 2000. It is then possible to produce the same kind of graphs as the one above, per year of birth. And here, we clearly observe the importance of the week end effect (maybe also because of color choice)

naissance=read.csv(
"http://freakonometrics.free.fr/naissanceFR.csv",
sep=";")
M=as.matrix(naissance[,3:ncol(naissance)])
BIRTH=as.vector(t(M))
YEAR=rep(1968:2005,each=12*31)
MONTH=rep(rep(1:12,each=31),38)
DAY=rep(1:31,12*38)
X=NA
for(y in 1968:2005){
sbase=base[YEAR==y,]
X=c(X,sbase$BIRTH/sum(sbase$BIRTH,
na.rm=TRUE))
}
base=data.frame(YEAR,MONTH,DAY,
BIRTH,BIRTHDAYPROB=X[-1])

m1=min(base$BIRTHDAYPROB,na.rm=TRUE)
m2=max(base$BIRTHDAYPROB,na.rm=TRUE)
y=1980
colr=rev(heat.colors(100))
sbase=base[YEAR==y,]
plot(0:1,0:1,col="white",xlim=c(-1,12),
ylim=c(-31,1),axes=FALSE,xlab=
paste("Naissance en",y,sep=" "),ylab="")
for(x in 1:nrow(sbase)){
a=sbase$MONTH[x];b=sbase$DAY[x]
polygon(c(a-.9,a-.9,a-.1,a-.1),-c(b-.9,b-.1,
b-.1,b-.9),col=colr[(sbase$BIRTHDAYPROB[x]-m1)/
(m2-m1)*100],border=NA)
}
text((1:12)-.5,.5,c("J","F","M","A","M","J","J",
"A","S","O","N","D"),cex=.7)
text(-.5,-(1:31)+.5,1:31,cex=.7)

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-France.txt",header=TRUE)
> Expo  <- read.table("http://freakonometrics.free.fr/
Exposures-France.txt",header=TRUE,skip=2)
> 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).

Les garçons ont-ils plus de sœurs que les filles ?

Ce qui est pénible avec les paradoxes, en tant que scientifique, c’est qu’ils nous obligent à remettre en cause notre intuition, sans cesse. Car mine de rien, on fonctionne tous à l’intuition, c’est un peu comme ça qu’on monte un programme de recherche: on a l’intuition que l’on pourra obtenir des résultats, mais tout reste à prouver.

Prenons la question suivante “les garçons on-t-ils plus de sœurs que les filles ?” en moyenne. N’importe qui avec un peu de bon sens répondrait “oui“, l’idée naturelle est que la fille ne peut pas se compter elle même comme sœur, donc les filles devraient en compter un peu moins… Et comme toujours la formulation est importante. On ne demande pas, par exemple “les garçons ont-t-ils plus chance d’avoir au moins une sœur que les filles ?” (auquel cas la réponse serait différente).

Maintenant réfléchissons un peu, et faisons un petit modèle où on suppose que la taille de la fratrie ne dépend pas des sexes de la fratrie. Et comptons, pour mieux voir ce qui se passe. Et commençons par un cas simple, par exemple avec trois enfants,

  • (G,G,G) – personne n’a de sœur
  • (G,G,F), (G,F,G) ou (F,G,G) – la fille n’a pas de sœur, par contre, il y a deux garçons avec une sœur
  • (G,F,F), (F,G,F) ou (F,F,G) – dans chaque cas, il y a deux filles avec une sœur, et un garçon avec deux sœurs,
  • (F,F,F) – ce qui fait trois filles avec chacune deux sœurs.

Si on récapitule pour les garçons, on a 3 cas avec 2 garçons ayant 1 sœur et 3 cas avec 1 garçon ayant 2 sœurs, soit (à un facteur de normalisation près) un nombre espéré de sœurs de 12. Et pour les filles, on a 3 cas avec 2 filles ayant 1 sœur et 1 cas avec 3 filles ayant 2 sœurs, soit (au même facteur de normalisation près) un nombre espéré de sœur de 12.
Posons http://freakonometrics.blog.free.fr/public/perso4/soeur07.gif la taille de la fratrie, que l’on va supposer fixé. Essayons de formaliser les calculs qu’on vient de faire. S’il y a http://freakonometrics.blog.free.fr/public/perso4/soeur08.gif filles, et http://freakonometrics.blog.free.fr/public/perso4/soeur11.gif garçons,

  • http://freakonometrics.blog.free.fr/public/perso4/soeur08.gif filles et http://freakonometrics.blog.free.fr/public/perso4/soeur11.gif garçons – chacun des http://freakonometrics.blog.free.fr/public/perso4/soeur11.gif garçons a http://freakonometrics.blog.free.fr/public/perso4/soeur08.gif sœurs, alors que chacune des http://freakonometrics.blog.free.fr/public/perso4/soeur08.gif filles a http://freakonometrics.blog.free.fr/public/perso4/soeur13.gif sœurs,

Bref, dans le calcul de l’espérance pour les garçons, on doit avoir des http://freakonometrics.blog.free.fr/public/perso4/soeur12.gif et pour les filles, http://freakonometrics.blog.free.fr/public/perso4/soeur14.gif. Et en faisant intervenir le nombre de combinaisons possibles de http://freakonometrics.blog.free.fr/public/perso4/soeur08.gif filles, et http://freakonometrics.blog.free.fr/public/perso4/soeur11.gif garçons, et par qui était notre facteur de proportionnalité de tout à l’heure (i.e. le nombre de possibilités d’agencement des sexes dans les familles avec http://freakonometrics.blog.free.fr/public/perso4/soeur07.gif enfants, quand on se focalise sur le sexe de tous, sauf un, qui est choisi au hasard, et qui est la personne qui compte ses soeurs), et en notant http://freakonometrics.blog.free.fr/public/perso4/soeur06.gif l’espérance du nombre de sœurs chez les garçons, et http://freakonometrics.blog.free.fr/public/perso4/soeur05.gif chez les filles, on a (merci à Pierre Antoine qui m’a fait remarqué une erreur dans mes calculs initiaux)

alors que

Un peu de calculs permet de montrer que ces deux quantités sont en fait égales… et valent

ce qui a du sens car il reste http://freakonometrics.blog.free.fr/public/perso4/soeur07.gif-1 personnes dans la fratrie (je peux renvoyer au livre de Ruma Falk, Understanding Probability and Statistics: A Book of Problems, page 46, qui passe par les fonctions génératrices pour une démonstration propre). Amusant, non ? En tous les cas suffisamment pour attirer mon attention alors que j’essaye de récupérer d’un soucis de santé, en dormant dans le lit de ma (grande) fille au milieu des peluches et des poupées…

Valeur de l’euro de rente

Il y a quelques années, j’avais été sollicité pour rédiger les Annexes de l’ouvrageévaluation du préjudice corporel, de Max le Roy. Comme l’indique le descriptif du livre, j’ai un peu la pression car pas mal de jugements rendus sont basés sur mes calculs, “La résolution 75-7 du comité des ministres du Conseil de l’Europe est intervenue pour faciliter l’harmonisation des législations et des jurisprudences en matière de réparation des dommages en cas de lésion corporelle ou de décès. Elle a proposé un certain nombre de principes qui, bien que non obligatoires pour les États membres, sont très largement suivis par les juridictions françaises. Ce sont ces principes et leur application en droit français qui sont exposés et commentés dans cet ouvrage qui rend compte, notamment, de la très importante réforme apportée en matière de recours subrogatoire des caisses de sécurité sociale contre les tiers par la loi n° 2006-1640 du 21 décembre 2006. Le barème de capitalisation de rentes tient compte des plus récentes tables de mortalité et prend pour base un taux d’intérêt de 5 % conforme aux données économiques et financières actuelles. Accompagné d’un barème fonctionnel indicatif des incapacités en droit commun, de tableaux de jurisprudence et de formules de jugement, l’ouvrage fait le point sur toutes les questions auxquelles sont confrontés les magistrats, avocats, experts, médecins, assureurs.”

Avec la mise à jour des tables de mortalité en France, j’ai du remonter mes manches et reprogrammer tous les calculs. Sauf que cette fois, j’ai proposé un petit code R pour générer à nouveau ces tables. Le but est de calculer le prix de mille euros de rente viagère ou temporaire limitée. Une petite option permet de dire si on veut la nouvelle table, ou l’ancienne (ce qui m’a permis de vérifier que je ne disais pas des bêtises).

LF=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/TV8890.csv",
header=TRUE,sep=";")
LH=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/TD8890.csv",
header=TRUE,sep=";") 
PRIX=function(annee=2011,age,sex="HOM",taux=0.04,
              duree,C=1000,old=FALSE){
if(old==FALSE){
file  = paste("http://freakonometrics.blog.free.fr/public/data/",
              sex,"-table-SPLx.csv",sep="")
B     = read.table(file,header=TRUE,sep=",")
an    = annee-age; if(an>2005){an=2005}
nom   = paste("X",an,sep="")
L     = B[,nom]}
if(old==TRUE){
if(sex=="HOM"){L = c(LH$Lx,rep(0,20))}
if(sex=="FEM"){L = c(LF$Lx,rep(0,20))}}
Q     = L[(age+1):length(L)]/L[(age+1)]
actualisation = (1+taux)^(0:min(duree,120-age))
prixsup = sum(Q[2:(min(duree,120-age)+1)]/
          actualisation[2:(min(duree,120-age)+1)] )
prixinf = sum(Q[1:(min(duree,120-age))]/
          actualisation[1:(min(duree,120-age))] )
return(C*c(prixsup,prixinf))}

Les tables TGH-05 et TGF-05 sont issues d’une population de rentiers et on peut consulter les tables ici (pour les 38 pages du Journal Officiel). Ces tables ont été construites sur une population différente des tables TH et TF (construites sur l’ensemble de la population française). La méthodologie est décrite ici (merci d’ailleurs à Fred qui a approuvé mon utilisation de ces tables).

Par exemple, on peut retrouver la valeur figurant dans les anciennes Annexes (ici),

> PRIX(age=50,sex="HOM",duree=150,taux=.05,old=TRUE)
[1] 13452.15 14452.15

En bouclant, on peut générer les tables données en Annexe.

AGES=c(150,65,60,55,25,20,18,16)
tableH=matrix(0,106,length(AGES))
for(i in 1:length(AGES)){
for(a in 0:min(c(AGES[i]-1,105))){
tableH[a+1,i]=PRIX(age=a,sex="FEM",duree=
        AGES[i]-a,taux=.04,old=FALSE)[1]
}}
write(tableH,"table-femme-nouvelle-04.csv",sep=";")

(on utilise la première valeur var on veut des rentes versées en fin d’année, en cas de vie). Avec un taux de 4%, on a pour les hommes la table ici, et pour les femmes la table . Si quelqu’un voit une coquille, merci de me signaler avant que le livre ne parte à l’éditeur…

Du sex-ratio en France

En novembre dernier, Baptiste @Coulmont m’avait envoyé un courriel correspondant à ce qu’il a mis en ligne sur son blog (ici) sur l’utilisation du fichier des prénoms français pour analyser le sex ratio à la naissance en France. J’attendais qu’il publie ses commentaires avant de mettre les miens (car le graphique qu’il a mis en ligne ce matin m’avait fait m’interroger).
Pour faire simple, la base de prénoms inclue un sexe. Mais si on regarde le rapport du nombre de garçon sur le nombre de filles, à la naissance, on a le graphique ci-dessous,

dat=read.table("nat2004.csv",sep=";",header=TRUE)
naissancesm=rep(NA,105)
for (i in 1900:2004) {
naissancesm[i-1899]=sum(dat[dat$annais==i&dat$sexe==1,"nombre"],
na.rm=TRUE)
}
naissancesf=rep(NA,105)
for (i in 1900:2004) {
naissancesf[i-1899]=sum(dat[dat$annais==i&dat$sexe==2,"nombre"],
na.rm=TRUE)
}
plot(1900:2004,naissancesm/naissancesf,col="red")

La tendance du début est très surprenante. Par exemple,  si on reprend les chiffres donnés par Pierre Simon Laplace sur les naissances à Paris entre 1750 et 1800 (mentionné hier, ici) on est déjà sur ratio de l’ordre de 1.05 (que l’on retrouve sur la fin de notre graphique, mais pas le début).
> 393386/377555
[1] 1.041930
 
> 251527/241945
[1] 1.039604

Donc il n’y a pas de raison d’avoir cette diff1rence. La conclusion semble être qu’il y a un soucis sur la base des prénoms. En effet, dans un rapport de Anouch Chahnazarian (ici), on retrouve l’évolution suivante,

que l’on peut rapprocher des données d’Éric Brian & Marie Jaisson (ici pour les données et quelques pages) utilisées dans leur ouvrage Le sexisme de la première heure, hasard et sociologie, qui mesure ici la fréquence de garçons à la naissance

Si l’on compare ces dernières données (via le fichier ici), calculé sur les données de l’INSEE et de l’INED, on retrouve un niveau très proche de celui que l’on a sur le fichier des prénoms, avec toutefois un biais constamment négatif.

b=read.table("http://freakonometrics.blog.free.fr/public/data/sex-ratio.txt")
X=b$V1
Y=b$V2/(100-b$V2)
lines(X,Y,col="blue")

Attention donc à la variable de sexe dans la base surtout avant guerre (en espérant que ce soit le seul soucis). Et je renvoie au blog de Baptiste qui publie toujours des choses très amusantes sur les prénoms, ici (je n’ai plus trop eu le temps de travailler dessus depuis les 5 billets en ligne ).

C’était si vieux que ça 77 ans ?

Lors de l’ouverture des cadeaux de Noël, et que nous discutions de Tintin avec les enfants, il y a eu un long débat (sans fin) sur “pourquoi on ne peut plus lire Tintin après 77 ans” ? Tintin est ce reporter inventé au début du XXème siècle par Hergé, dont le nom a donné naissance au Journal de Tintin en septembre 1946 (ou plutôt à Kuifje, dans la version flamande, qui a été traduite en français en 1948), dont le sous-titre était « Le journal des jeunes de 7 à 77 ans ». Bon, autant le 7 ans peut se comprendre (même si ma fille prétendait que quand même, elle avait le droit de lire Tintin si elle voulait), mais le 77 ans a donné lieu a de grandes discussions… Surtout que d’un air docte, j’ai du affirmer “en 1948, après la guerre, 77 ans c’était très très très vieux” (genre Mathusalem quoi, ou Dumbledore, pour reprendre une référence qu’ils connaissent) “alors que maintenant...

Ah, si j’avais su… Pour rappel, en 1948, en France, l’espérance de vie à la naissance était de 74 ans, environ, et la fonction de survie (à la naissance) ressemblait à ça,

autrement dit, la probabilité qu’une personne née en 1948 atteigne 77 ans était de 2 chances sur 3. En fait, si on regarde pour les personnes en vie en 1948 (disons nées après 1900 pour faire simple), la probabilité d’atteindre 77 ans en fonction de leur année de naissance ressemblait à ça,

ou de manière duale, le quantile à 66.67% est (à la naissance)

Si on compare à aujourd’hui, un 77 ans en 1948 correspondrait à un 85 ans en 2010 (en terme d’age atteint avec la même probabilité). Et pas 99 ans comme cela est dit parfois sur les boites de jeux…. Moralité, faut jamais croire son père quand il veut faire son intéressant les soirs de Noël ! (mais les autres jours, si… non mais).

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 http://freakonometrics.blog.free.fr/public/maths/viekt.png‘s since we have identifiability problems here. Nevertheless, we if consider affine transformation so that  http://freakonometrics.blog.free.fr/public/maths/viekt.png‘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 http://freakonometrics.blog.free.fr/public/maths/viekt.png‘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 ?

 

Davantage de régularité dans les naissances ?

L’autre jour, l’INED avait fait beaucoup de bruit en publiant une étude expliquant que le pic des naissance classiquement observé en mai semblait avoir disparu (ici ou ). On va essayer de vérifier ce point, à partir de toutes les naissances observées en France entre 1968 et 2005. Ça sera le billet sexe de l’été…

  • Le nombre quotidien de naissances en France

La série brute ressemble à ça,

On pourrait à la rigueur extraire les cycles directement sur cette série, mais ça ne serait pas forcément convainquant (même si visuellement on retrouve ce qui était évoqué),

mais la première chose que l’on observe est que la différence entre le nombre de naissances le week end (en rouge) et en semaine (en bleu) s’est fortement creusée dans les années 80-90.

On peut regarder le ratio du nombre de naissance le week end, par rapport aux jours de la semaine précédant le week end,

et on peut alors “corriger” ou “lisser” la série, de manière à maintenir la moyenne hebdomadaire, mais en transférant quelques naissances de la semaine au week end.

Sur cette série lissée, on peut extraire une tendance de long terme.

C’est sur le “bruit” autour de cette tendance que l’on peut essayer d’étudier un éventuel cycle annuel. Notons que l’on retrouve les tendances sur la série brute

  • Série corrigée de la tendance de long terme

La série suivante permet de voir l’évolution du “cycle” (série brute à laquelle on retranche la tendance observée ci-dessus) à la fin des années 60 (en rouge), des années 70 (en orange), des années 80 (en jaune), voire des années 90 (en presque blanc),

Sur les années 68-78, on a la tendance suivante

alors que pour les années 95-2005, cela donne

Au final, les deux tendances sont les suivantes, avec en rouge les années 70 et en bleules années 2000.

Effectivement, le pic classiquement observé au mois de mai semble avoir disparu, et la série aujourd’hui semble beaucoup plus lisse. Et comme l’INED, on note qu’il y a effectivement un mini-pic au mois de septembre.

  • Des naissances au comportement sexuel

Tout comme l’INED, on peut relier ça à l’activité sexuelle des français, que l’on peut visualiser sur le graphique ci-dessous (en décalant tout simplement de 9 mois),

Bref, l’activité sexuelle semblait s’intensifier pendant les vacances scolaires, en particulier début août (mais aussi plus généralement pendant les vacances scolaires, que ce soit Pâques, les vacances d’été, ou les vacances de Noël d’ailleurs). Ou peut être pourrait-on conclure à la rigueur que l’été, on sort moins couvert…. car on observe les comportements sexuels associés à une naissance 9 mois plus tard (ce qui ne doit pas correspondre à l’unique activité sexuelle).
Les français semblent avoir intensifié leur activité sexuelle pendant une période plus fraîche (novembre-mars), au détriment de la période estivale. Un effet du réchauffement climatique (il fait trop chaud pour batifoler l’été) ? Un effet des 35 heures et des RTT (et des vacances qui sont plus réparties dans l’année qu’avant) ?

  • L’été s’ra chaud ?

En fait, au lieu de penser en terme de saison, on peut penser en terme de température (j’ai pris la température quotidienne minimale à Paris, c’est à dire la nuit… ).

c’est à dire en zoomant

On peut aussi chercher à étudier davantage ce qui se passait il y a 30 ans et maintenant (comme dans la partie précédente). Il y a 30 ans (courbe rouge) on observait un (très) fort pic estival. Ce pic a disparu aujourd’hui (courbe bleue), avec un comportement comparable, qu’il fasse 0 ou 20 degrés dehors.

  • Quid de l’impact des bourses ?

Mais la recherche de la variable qui pourrait expliquer cette saisonnalité ou ces pics est sans fin. Beaucoup de monde évoque la bonne santé économique pour expliquer les booms de natalité…. ça marche peut être à long terme, mais à court terme ?
En fait, si on croisse l’activité sexuelle avec le niveau de la bourse, on peut presque retrouver ce genre de comportement,

une bourse en très forte hausse semble provoquer un (léger) pic de natalité 9 mois plus tard, alors qu’un effondrement de la bourse se traduit plutôt par une chute (légère certes) 9 mois plus tard…

Les abaques comme outils de pédagogie statistique

Comme je l’avais mentionné à propos d’un vieux papier paru dans Le Canard, c’est souvent compliqué de communiquer sur les retraites. C’est un sujet complexe et assez rapidement technique. Historiquement, les abaques sont les tables de calculs,

Mais il faut admettre que retrouve un peu n’importe quoi dans ce terme d’abaque…. Certaines sont évoquées ici ou , voire pour quelque chose de très joli. On notera qu’il existe une idée sous-jacente commune, à savoir construire un outil graphique permettant de visualiser comment 3 paramètres peuvent être liés (généralement, ils sont liés par une fonction complexe).

  • Abaques et retraites

Parmi les outils que l’on retrouve dans les documents du COR, j’ai été surpris de voir ressortir l’utilisation des abaques (par exemple ici). On nous en fait même la promotion en vidéo,

Pour la petite histoire, j’avais évoqué son utilisation dans un graphique ici sur les problèmes d’échantillonnage. C’est effectivement un outil graphique intéressant, à condition de savoir le décoder…
L’idée est simple, à savoir lier les trois paramètres essentiels(ou supposés comme tels) du bilan d’un système de retraite, à savoir

  • le niveau des prélèvements,
  • le niveau de l’âge effectif de départ en retraite
  • le niveau des pensions.

La difficulté est de trouver une forme graphique permettant de relier ces trois niveaux (la dimension 3 étant souvent trop abstraite). Par exemple dans les sondages, pour construire un intervalle de confiance, les abaques sont utilisées (comme ici) pour relier

  • le niveau des intervalles de confiance (5%-95% ou 10%-90%)
  • le niveau de probabilité attendu dans le sondage
  • la taille de l’échantillon interrogé

On peut d’ailleurs trouver d’autres exemple en statistiques ici ou . Mais pour revenir au problème des retraite, sur le graphique ci-dessous (appendix 1 ici), en abscisse on met le taux de remplacement (c’est à dire le ratio de la pension moyenne par rapport au revenu moyen) et en ordonnée on met le niveau d’un hausse du taux de cotisation (en points de cotisation). chaque ligne représente un âge moyen de départ en retraite effectif.

On se fixe alors des scénarios d’équilibre et on regarde la manière dont se déplacent les équilibres. J’avais prévu toute une analyse technique, mais je viens de me rendre compte qu’Antoine a mis en ligne un billet remarquable sur le sujet sur http://www.ecopublix.eu/2010/…, que je ne pourrais jamais égaler en clarté. Donc je vais proposer une autre construction d’abaque afin de montrer son intérêt pédagogique.

  • Comprendre l’actualisation

Sur le site http://images.math.cnrs.fr/…, Xavier Caruso explique son émerveillement devant la difficulté de faire des calculs d’actualisation. En un sens, ça me rassure de voir que les calculs actuariels ne paraissent pas triviaux à tout le monde (et surtout plein de contre-intuitions). Xavier essaye d’expliquer la différence entre

  • ouvrir un compte rémunéré à 4,5% sur lequel on dépose 500 euros chaque mois pendant 10 ans
  • ouvrir un compte rémunéré à 4,5% sur lequel on dépose 60000 euros initialement sans plus jamais l’alimenter ensuite.

Si effectivement la somme des 120 versements de 500 euros correspond aux 60000 euros, la valeur intégrant les intérêts n’est pas du tout la même au bout de 10 ans, car dans le premier car, comme le disent les banquiers, “les intérêts se capitalisent“. Bon, si Xavier fait des maths (et donc des développements limités par avoir des formules fermés), on va plutôt faire ici un peu de calculs numériques. Tout d’abord, pour visualiser l’évolution des montants disponibles sur les comptes, rien de plus simple,

> tps=1/12*(1:120)
> plot(tps,60000*(1+0.045)^tps,ylim=c(0,100000))
> lines(tps,cumsum(500*(1+0.045)^(tps)))

les montants au bout de 10 ans étant

> sum(500*(1+0.08525)^(1/12*(1:120)))
[1] 93180.33
> sum(500*(1+0.045)^(1/12*(1:120)))
[1] 75514.32

Ce qui confirme l’intuition que nous avions. Comme le note Xavier, il faudrait un taux d’intérêt deux fois plus grand pour le placement mensuel, la preuve

> sum(500*(1+0.09)^(1/12*(1:120)))
[1] 95543
> sum(500*(1+0.08525)^(1/12*(1:120)))
[1] 93180.33

Xavier obtient cette valeur en effectuant des développements limités. Il résume ça à l’aide du dessin

On peut aller plus loin, en notant qu’à taux sur le placement mensuel donné, ainsi que la maturité, le taux sur le placement bloqué qui rapporte la même valeur à échéance est donné par

 

>1))/(500*12*10))^(1/10)-1
[1] 0.04762238

(dans le cas d’un placement rapportant 9%). En fait, si on regarde le montant dont on dispose sur les comptes à chaque date, on obtient précisément les valeurs suivantes, en fonction du temps,

on retrouve effectivement qu’à 4,5%, les deux placements ne donnent pas le même montant à échéance (ici 10 ans). Alors que si le placement alimenté tous les mois rapportait 8,525%, on aurait exactement la même somme au bout de 10 ans

ou si on compare avec un placement à 9% (le double, ce qui est obtenu par développement limité), on obtient

Pour trouver la valeur exacte du multiplicateur en fonction du taux offert sur le compte bloqué et la maturité, la fonction est tout simplement

> taux=4.5/100
> T=10
> f=function(k){(1*12*T)*(1+taux)^T-sum(1*(1+k*taux)^(1/12*(1:(12*T)))) }
> uniroot(f,interval=c(0,20))$root
[1] 1.894332

L’idée d’utiliser une abaque peut s’avérer intéressante, car on essaye ici de relier

  • le taux de rendement de l’argent bloqué
  • le taux proposé pour le placement alimenté tous les mois
  • la durée envisagée du placement

Je peux mettre en abscisse le taux du placement bloqué. En ordonné, je met le multiplicateur utilisé pour le placement alimenté tous les mois. A maturité donné (par exemple 10 ans en rouge), on peut regarder la valeur du multiplicateur pour que les valeurs à échéance coïncident,

On retrouve directement que le taux équivalent pour un placement alimenté tous les mois doit être 1,89 fois le taux offert sur le placement bloqué si on se fixe une échéance de 10 ans pour égaliser les placements. Le code est simplement le suivant

> tx=seq(0.005,.05,by=.001)
> tps=seq(5,25,by=.25)
> M=matrix(NA,length(tps),length(tx))
> for(i in 1:length(tps)){
+ for(j in 1:length(tx)){
+ T=tps[i]
+ taux=tx[j]
+ f=function(k){
+ (1*12*T)*(1+taux)^T-sum(1*(1+k*taux)^(1/12*(1:(12*T))))
+ }
+ M[i,j]=uniroot(f,interval=c(0,2))$root
+ }}
> contour(tps,tx,M,lwd=2,
+ xlab="Durée du placement (en années)",ylab="Taux (du placement bloqué)")
> abline(v=seq(5,25,by=1),lty=2,col="grey")
> abline(h=seq(0,.05,by=.005),lty=2,col="grey")

Bref, l’actualisation est un sujet délicat (même ce n’est que des calculs de sommes et de puissances). Et faire des petits dessins permet souvent de mieux comprendre…

  1. sum(500*(1+0.09)^(1/12*(1:120 []

Les tables de mortalité

Chose promise, chose dûe, un court billet expliquant les principales tables utilisées,

  • Les tables TD et TV 88-90

Ces tables datent un peu, et si je continue à en parler en cours, c’est parce qu’elles sont simples à utiliser (et pour continuer à me faire croire que je n’ai pas vieilli depuis mes études). Cette table est d’ailleurs tellement sérieuse qu’on la retrouve dans la loi (ici), dans un arrêté d’avril 1993. La table dite TD 88-90 (pour Décès) a été établie par l’INSEE suite aux observations réalisées entre 1988 et 1990 sur une population d’hommes. Elle était appliquée pour le calcul des primes des contrats d’assurance décès. La table dite TV 88-90 (pour Vie) a été établie par l’INSEE suite aux observations réalisées entre 1988 et 1990 sur une population de femmes. Elle était appliquée pour le calcul des primes des contrats d’assurance en cas de vie. Ces tables peuvent se récupérer à l’aide des codes suivants,

> TD=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/TD8890.csv",
+ sep=";",header=TRUE)
> TV=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/TV8890.csv",
+ sep=";",header=TRUE)

Ces tables ont été remplacées par les tables dites TH et TF, respectivement.

  • Les tables TH et TF 00-02

Ces tables ont été établies à partir des données INSEE de la population française entre 2000 et 2002 et ont été lissées. Ce sont des tables générationnelles, qui nécessitent un correctif d’âge pour tenir compte des écarts de mortalité entre générations. Elles sont applicables à partir du 1er janvier 2006. L’institut des actuaires a proposé une “notice d’utilisation” en ligne ici, et Cimon en avait parlé sur son blog ().

  • Les tables TPRV 95

La table TPRV 93 (pour Table Prospective de Rente Viagère) est un extrait de la table dite plancher pour la tarification des contrats de rente viagère. Elle a été publiée par l’arrêté du 28 juillet 1993 (ici sans les annexes), et correspond à une table prospective qui retrace la mortalité des générations 1887 à 1993 (les tables prospectives sont au programme de Master 2).

La TPRV 93 représente la table complète de la génération 1950. La table est en ligne ici (en csv) lisible sous R avec le code suivant,

> TD=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/TPRV.csv",
+ sep=";",header=TRUE)

Mais ce n’est que la première étape. On utilise ensuite un décallage d’âge.

On a parlé de ce point en TD (mais sur les taux de mortalité), c’est aussi ce qui s’appelle l’hypothèse de Rueff, et traduit la diminution des taux de mortalité par un rajeunissement (i.e. un gain sur les âges), i.e.

Notons qu’en fait,  dépend du niveau des taux d’intérêt. L’ensemble des tables de décallage sont en ligne ici.

  • Les tables TGH et TGF 05

Il s’agit ici, là encore, d’une population de rentiers, comme pour la (ici pour les 38 pages du code). Aussi, ces tables ont été construites sur une population différente des tables TH et TF (construites sur l’ensemble de la population française). La méthodologie est décrite ici.