# Forecasting Techniques for Demographics

Tomorrow, we will discuss in our cours forecasting tools for demographics. But first, we will see basic static tools. Before playing with longitudinal dataset, let us use “standard” life tables. Some (French) datasets are available on the INED website, but let us use the most popular ones, the TV8890 and TD8890. Those two tables can be downloaded from wikipedia,

url="https://fr.wikipedia.org/wiki/Table_de_mortalit%C3%A9" download.file(url,"mortalite.html") library(XML) tables=readHTMLTable("mortalite.html")

The code to get the dataset is the following (here, we should be careful since there is a space in the number, as a separator for thousands)

TV8890_0=tables[[2]] a1=as.numeric(as.character(TV8890_0[,1])) a2=as.numeric(as.character(TV8890_0[,3])) espace=function(x) gsub("[[:space:]]", "", x) b1=espace(as.character(TV8890_0[,2])) b2=espace(as.character(TV8890_0[,4])) TV8890=data.frame(x=c(a1,a2),lx=as.numeric(c(b1,b2))

One can also read the second life table (note that there is a typo in wikipedia since here, the two tables are exactly the same)

TD8890_0=tables[[1]] a1=as.numeric(as.character(TD8890_0[,1])) a2=as.numeric(as.character(TD8890_0[,3])) b1=espace(as.character(TD8890_0[,2])) b2=espace(as.character(TD8890_0[,4])) TD8890=data.frame(x=c(a1,a2),lx=as.numeric(c(b1,b2)))

It is possible to use that survival function to compute some sort to life expectancy at birth

sum(TV8890$lx)/100000-1 [1] 72.01518 One can visualize the survival probability (up to a 100,000 scaling constant) plot(TV8890,type="l") or the death probability, i.e. the probability to die at some specific age $x$, given that you did reach age $x$, also called the force of mortality n=nrow(TV8890) px=(TV8890$lx[1:(n-1)]-TV8890$lx[2:n])/ TV8890$lx[1:(n-1)] x=TV8890$x[1:(n-1)] plot(x,px,type="l",xlab="age") A more popular visualization is obtained with a log scale for the probability plot(x,px,type="l",log="y") Finally, we can compute the density of the age at death pbx=TV8890$lx[1:(n-1)]*px/100000 plot(x,pbx,type="l")

that can also be used to compute life expectancy

sum(x*pbx) [1] 72.01518

That is for the static case. For longitudinal tables, we can use those from the Human Mortality Database. For instance, for France, we can get information of the number of death at age $x$ during year $t$, as well as the exposure (number of people alive). For France, there are datasets available here

url="http://freakonometrics.free.fr/FranceDeaths_1x1.txt" download.file(url,"FRD.txt") url="http://freakonometrics.free.fr/FranceExposures_1x1.txt" download.file(url,"FRE.txt")

url="http://freakonometrics.free.fr/CanadaDeaths_1x1.txt" download.file(url,"CAD.txt") url="http://freakonometrics.free.fr/CanadaExposures_1x1.txt" download.file(url,"CAE.txt")

The following code can be use to read those files.

FRD=read.table("FRD.txt",skip = 3,header=TRUE) tail(FRD) Year Age Female Male Total 22195 2015 105 297.96 35.86 333.82 22196 2015 106 182.95 20.39 203.34 22197 2015 107 104.87 10.50 115.37 22198 2015 108 57.27 5.07 62.34 22199 2015 109 31.93 2.59 34.52 22200 2015 110+ 33.03 1.61 34.64

# Proportion of people alive in 1945 that are still alive

In demography, we like to use life tables to estimate the probability that someone born in 1945 (say) is still alive nowadays.  But another interesting quantity might be the probability that someone alive in 1945 is still alive nowadays.

The main difference is that we do not know when that person, alive in 1945, was born. Someone who was old in 1945 is very unlikely still alive in 2017. To compute those probabilities, we can use datasets from http://www.mortality.org/hmd/. More precisely, we need both death and birth data. I assume that datasets (text files) were downloaded (it is necessary to register – for free – to get the data).

D=read.table("FRDeaths_1x1.txt",skip=1,header=TRUE) B=read.table("FRBirths.txt",skip=1,header=TRUE)

In the death dataset, there is a “110+” for people older than 110 years. For convenience, let us cap our observations at 110 years old,

D$Age=as.numeric(as.character(D$Age)) D$Age[is.na(D$Age)]=110

Consider now a first function that will return, for people born in 1930 (say) two informations

• the number of people (here, let us consider women only) born in 1930 (from the birth database)
• the number of death of people of age 0 in 1930, people of age 1 in 1931, people of age 2 in 1932, etc…

The code is simple

+   date2=as.Date(paste(data$date_death_y,"-",data$date_death_m,"-",data$date_death_d,sep=""),"%Y-%m-%d") + idx=which(!(is.na(date1)|is.na(date2))) + date1=date1[idx] + date2=date2[idx] + itg=try(age<-age_years(date1,date2),silent=TRUE) + if(inherits(itg, "try-error")) age=trunc((date2-date1)/365.25) + w=weekdays(date2) + T=table(age,w) + Tab=matrix(0,106,7) + for(i in 1:nrow(T)) if(as.numeric(rownames(T)[i])<106) Tab[as.numeric(rownames(T)[i]),]=T[i,] + return(Tab) + } + D <- lapply( seq_len(nrow( go_through)),count_birthday) + T=D[[1]] + for(s in 2:length(D)) T=T+D[[s]] + return(T) + } If we run that function on the three files > D1=TABLE_AGE_DAY("ssdm1") |========================================| 100% > D2=TABLE_AGE_DAY("ssdm2") |========================================| 100% > D3=TABLE_AGE_DAY("ssdm3") |========================================| 100% we can visualize not percentages, as on the figure above, but counts > D=D1+D2+D3 > colnames(D)= c("Sun","Thu","Mon","Tue","Wed","Sat","Fri") > D=D1[, c("Sun","Mon","Tue","Wed","Thu","Fri","Sat")] and we have here (I remove the Saturday to get a better output) > D[,1:6] Sun Mon Tue Wed Thu Fri [1,] 2843 2888 2943 3020 2979 3038 [2,] 2007 1866 1918 1974 1990 2137 [3,] 1613 1507 1532 1530 1515 1613 [4,] 1322 1256 1263 1259 1207 1330 [5,] 1155 1061 1092 1128 1112 1171 [6,] 1067 985 950 1082 1009 1055 [7,] 1129 901 915 954 941 1044 [8,] 1026 927 944 935 911 1005 [9,] 1029 1012 871 908 939 998 [10,] 1093 1011 974 958 928 1018 [11,] 1106 1031 1019 1036 1087 1122 [12,] 1289 1219 1176 1215 1141 1292 [13,] 1618 1455 1487 1484 1466 1633 [14,] 2121 2000 1900 1941 1845 2138 [15,] 2949 2647 2519 2499 2524 2748 [16,] 4488 3885 3798 3828 3747 4267 [17,] 5709 4612 4520 4422 4443 5005 [18,] 7280 5618 5400 5271 5344 5986 [19,] 8086 6172 5833 5820 6004 6628 [20,] 8389 6507 6166 6055 6430 6955 [21,] 8794 7038 6794 6628 6841 7572 [22,] 8578 6528 6512 6472 6757 7342 [23,] 8345 6750 6483 6469 6714 7338 [24,] 8361 6859 6589 6623 6854 7369 [25,] 8398 6974 6892 6766 6964 7613 [26,] 8432 7210 7012 7175 7343 7801 [27,] 8757 7641 7526 7352 7674 7950 [28,] 9190 8041 7843 7851 7940 8268 [29,] 9495 8409 8555 8400 8469 8934 [30,] 9876 9041 9015 9166 9106 9641 [31,] 10567 9952 9506 9634 9770 10212 [32,] 11417 10428 10402 10275 10455 11169 [33,] 11992 11306 11124 11095 11243 11749 [34,] 12665 12327 11760 12025 12137 12443 [35,] 13629 13135 13179 13037 12968 13724 [36,] 14560 14009 13927 13822 14105 14436 [37,] 15660 14990 15013 15009 15101 15700 [38,] 16749 16504 16148 16091 15912 16863 [39,] 17815 17760 17519 17144 17553 17943 [40,] 19366 19057 18918 18517 18760 19604 [41,] 20770 20458 20154 20339 20349 21238 [42,] 21962 22194 22020 21499 21690 22347 [43,] 23803 23922 23701 23681 23437 24227 [44,] 25685 26133 25559 25209 25287 26115 [45,] 27506 28110 27363 27042 27272 28228 [46,] 29366 29744 29555 29245 29678 30444 [47,] 31444 32193 31817 31504 31753 32302 [48,] 33452 34719 33529 33954 33441 34618 [49,] 36186 37150 36005 36064 36226 37138 [50,] 38401 39244 38813 38465 38506 39884 [51,] 40331 41830 41168 41110 40937 42014 [52,] 43181 44351 43975 43949 43579 44734 [53,] 45307 47134 46522 46149 46089 47286 [54,] 47996 49441 49139 48678 48629 49903 [55,] 50635 52424 51757 51433 51477 52550 [56,] 53509 55337 54556 54482 54406 55906 [57,] 55703 58482 58016 57400 57097 58758 [58,] 59016 61453 60652 61024 60557 62473 [59,] 62475 65651 64169 63824 63829 65592 [60,] 66621 69185 68885 68217 68752 69963 [61,] 69759 73144 72421 71784 71745 73414 [62,] 80346 84253 83044 83177 82416 83833 [63,] 86851 90059 89002 88985 89245 90334 [64,] 91839 95465 94602 93985 94154 96195 [65,] 98461 102846 101348 101328 101306 103170 [66,] 104569 108722 107768 107711 107729 109350 [67,] 111230 115477 114418 114743 113935 116356 [68,] 116999 122053 120727 120342 119782 122926 [69,] 123695 128339 127184 126822 126639 129037 [70,] 129956 136123 134555 135120 133842 137390 [71,] 137984 142964 141316 142855 141419 143620 [72,] 145132 150708 148407 149345 149448 151910 [73,] 152877 157993 155861 156349 155924 158725 [74,] 159109 164652 162722 163499 163157 165744 [75,] 165848 172121 170730 170482 170585 173431 [76,] 172457 179036 177185 177328 177392 180215 [77,] 179936 185015 183223 183932 183237 186663 [78,] 185900 191053 189986 189730 189639 193038 [79,] 191498 196694 194246 194810 195246 197812 [80,] 195505 201289 199684 199561 198968 203226 [81,] 199031 204927 202204 202622 202951 205792 [82,] 201589 207928 204929 204001 204396 208224 [83,] 201665 206743 205194 204676 205256 207980 [84,] 200965 205653 203422 202393 203422 206012 [85,] 197445 202692 199498 199730 200075 201728 [86,] 192324 195961 193589 194754 193800 196102 [87,] 183732 188063 185153 186104 186021 188176 [88,] 174258 177474 175822 176078 176761 177449 [89,] 163180 166706 162810 164367 164281 166436 [90,] 149169 151738 150148 150212 150535 152435 [91,] 134218 136866 134959 134922 135027 136381 [92,] 118936 121106 119591 119509 119793 120998 [93,] 102734 104955 102944 102865 103345 104776 [94,] 87418 88885 88023 86963 87546 87872 [95,] 72023 72698 72151 71579 71530 72287 [96,] 56985 58238 57478 57319 57163 57615 [97,] 44447 45058 44607 44469 43888 44868 [98,] 33457 34132 33022 33409 33454 33642 [99,] 24070 24317 24305 24089 24020 24383 [100,] 17165 17295 16755 17115 16957 17207 [101,] 11799 12125 11709 11816 11824 11719 [102,] 7714 7741 7959 7691 7648 7633 [103,] 5024 5012 4822 4792 4882 4916 [104,] 2987 3101 2978 3049 3093 2906 [105,] 1781 1894 1811 1756 1734 1834 So clearly, for young people, the number of deaths is rather small… And to visualize it, as above, we can use > P=D/apply(D,1,sum)*100 > range(P) [1] 12.34857 17.59386 > dP=trunc((P-min(P))/(max(P)+.01-min(P))*11) > library(RColorBrewer) > CLR=rev(brewer.pal(name="RdYlBu", 11)) > plot(0:1,0:1,ylim=c(55,110),xlim=c(-1,7)) > for(i in 1:106){ + for(j in 1:7){ + rect(j-1,108-i,j,107-i,col=CLR[dP[i,j]]) + }} > text(rep(-.5,106),107.5-1:106,0:105,cex=.4) As above, we observe a strong difference among weekdays for the date of death for young people (below 30) which disappear after (even if there is still a sunday effect) # Men set to live as long as women by 2030? A few months ago, in Men set to live as long as women, figures show, it was mentioned that (in the U.K.) the gap between male and female life expectancy is closing and men could catch up by 2030, according to an adviser for the Office for National Statistics. (the slides are available online http://cass.city.ac.uk/…). # Your Life in Weeks This week, I discovered a picture on http://waitbutwhy.com/, which represent a (so-called) typical human life, in weeks, I found that interesting. But the first problem is that I don’t understand the limit, below: 90 years, that’s not the average life length. That’s not what you should expect to live when you get born. The second problem is that it cannot be as static as it might seem, when you look at the picture. I mean, life expectancy at age 0 is not the same as life expectancy at age 30, or 50. So I did try to make an animated graph, using prospective life tables. Here a code to generate life tables, at different period, for a French population (I distinguish, here male and female) library(demography) france.LC1 <- lca(fr.mort,adjust="e0",series="female",years=c(1900,2100)) france.fcast <- forecast(france.LC1,h=100) L2 <- lifetable(france.fcast) ex2=L2$ex
L1=lifetable(fr.mort,series="female")
ex1=L1$ex exF=cbind(ex1,ex2) france.LC1 <- lca(fr.mort,adjust="e0",series="male",years=c(1900,2100)) france.fcast <- forecast(france.LC1,h=100) L2 <- lifetable(france.fcast) ex2=L2$ex
L1=lifetable(fr.mort,series="male")
ex1=L1$ex exM=cbind(ex1,ex2) Y=colnames(exF) Based on those lifetables, we can extract remaining life expectancy, at various ages (say, for instance 50, 51, 52, etc), for someone born on some given year (say 1950). Based on those expected remaining lifetimes, we can plot picture=function(yearborn=1950,age=50){ k=which(Y==yearborn) M=diag(exM[,k+0:100]) F=diag(exF[,k+0:100]) par(mfrow=c(1,2)) va=0:(52*100-1) plot(va%%52,va%/%52,cex=.6,pch=15,col=c("light yellow","light blue","white")[1+ (va>=age*52)*1+(va>(age+M[age+1])*52)*1],ylim=c(100,0),axes=FALSE,xlab="Week", ylab="Age",main=paste("Man, born on ",yearborn, ", age ",age,sep="")) axis(1) axis(2) plot(va%%52,va%/%52,cex=.6,pch=15,col=c("light yellow","pink","white")[1+ (va>=age*52)*1+(va>(age+F[age+1])*52)*1],ylim=c(100,0),axes=FALSE,xlab="Week", ylab="Age",main=paste("Woman, born on ",yearborn, ", age ",age,sep="")) axis(1) axis(2)} For instance, if we want the graph above, for someone age 30, born in 1980, we use picture(1980,30) Now, if we run a code to get an animated gif, we can get, for someone born in 1950, and for someone born in 2000 Now, if I could get historical datasets, with the average time spent in schools, ages of retirement, etc, I guess I could add it on the graph. But that’s another story… # Smoothing mortality rates This morning, I was working with Julie, a student of mine, coming from Rennes, on mortality tables. Actually, we work on genealogical datasets from a small region in Québec, and we can observe a lot of volatiliy. If I borrow one of her graph, we get something like Since we have some missing data, we wanted to use some Generalized Nonlinear Models. So let us see how to get a smooth estimator of the mortality surface. We will write some code that we can use on our data later on (the dataset we have has been obtained after signing a lot of official documents, and I guess I cannot upload it here, even partially). DEATH <- read.table( "http://freakonometrics.free.fr/Deces-France.txt", header=TRUE) EXPO <- read.table( "http://freakonometrics.free.fr/Exposures-France.txt", header=TRUE,skip=2) library(gnm) D=DEATH$Male
E=EXPO$Male A=as.numeric(as.character(DEATH$Age))
Y=DEATH$Year I=(A<100) base=data.frame(D=D,E=E,Y=Y,A=A) subbase=base[I,] subbase=subbase[!is.na(subbase$A),]

The first idea can be to use a Poisson model, where the mortality rate is a smooth function of the age and the year, something like

$D_{x,t}\sim\mathcal{P}(E_{x,t}\cdot \exp[{\color{blue}s(x,t)}])$that can be estimated using

library(mgcv)
regbsp=gam(D~s(A,Y,bs="cr")+offset(log(E)),data=subbase,family=quasipoisson)
predmodel=function(a,y) predict(regbsp,newdata=data.frame(A=a,Y=y,E=1))
vX=trunc(seq(0,99,length=41))
vY=trunc(seq(1900,2005,length=41))
vZ=outer(vX,vY,predmodel)
ylab="Years (1900-2005)",zlab="Mortality rate (log)")

The mortality surface is here

It is also possible to extract the average value of the years, which is the interpretation of the $a_x$ coefficient in the Lee-Carter model,

predAx=function(a) mean(predict(regbsp,newdata=data.frame(A=a,
Y=seq(min(subbase$Y),max(subbase$Y)),E=1)))
plot(seq(0,99),Vectorize(predAx)(seq(0,99)),col="red",lwd=3,type="l")

We have the following smoothed mortality rate

Recall that the Lee-Carter model is

$D_{x,t}\sim\mathcal{P}(E_{x,t}\cdot \exp[{\color{blue}a_x+b_x\cdot k_t}])$

where parameter estimates can be obtained using

regnp=gnm(D~factor(A)+Mult(factor(A),factor(Y))+offset(log(E)),
data=subbase,family=quasipoisson)
predmodel=function(a,y) predict(regnp,newdata=data.frame(A=a,Y=y,E=1))
vZ=outer(vX,vY,predmodel)
ylab="Years (1900-2005)",zlab="Mortality rate (log)")

The (crude) mortality surface is

with the following $a_x$ coefficients.

plot(seq(1,99),coefficients(regnp)[2:100],col="red",lwd=3,type="l")

Here we have a lot of coefficients, and unfortunately, on a smaller dataset, we have much more variability. Can we smooth our Lee-Carter model ? To get something which looks like

$D_{x,t}\sim\mathcal{P}(E_{x,t}\cdot \exp[{\color{blue}s_a(x)+s_b(x)\cdot s_k(t)}])$

Actually, we can, and the code is rather simple

library(splines)
knotsA=c(20,40,60,80)
knotsY=c(1920,1945,1980,2000)
regsp=gnm(D~bs(subbase$A,knots=knotsA,Boundary.knots=range(subbase$A),degre=3)+
Mult(bs(subbase$A,knots=knotsA,Boundary.knots=range(subbase$A),degre=3),
bs(subbase$Y,knots=knotsY,Boundary.knots=range(subbase$Y),degre=3))+
offset(log(E)),data=subbase, family=quasipoisson)
BpA=bs(seq(0,99),knots=knotsA,Boundary.knots=range(subbase$A),degre=3) BpY=bs(seq(min(subbase$Y),max(subbase$Y)),knots=knotsY,Boundary.knots= range(subbase$Y),degre=3)
predmodel=function(a,y)
predict(regsp,newdata=data.frame(A=a,Y=y,E=1)) v
Z=outer(vX,vY,predmodel)
ylab="Years (1900-2005)",zlab="Mortality rate (log)")

The mortality surface is now

and again, it is possible to extract the average mortality rate, as a function of the age, over the years,

BpA=bs(seq(0,99),knots=knotsA,Boundary.knots=range(subbase$A),degre=3) Ax=BpA%*%coefficients(regsp)[2:8] plot(seq(0,99),Ax,col="red",lwd=3,type="l") We can then play with the smoothing parameters of the spline functions, and see the impact on the mortality surface knotsA=seq(5,95,by=5) knotsY=seq(1910,2000,by=10) regsp=gnm(D~bs(A,knots=knotsA,Boundary.knots=range(subbase$A),degre=3)+
Mult(bs(A,knots=knotsA,Boundary.knots=range(subbase$A),degre=3), bs(Y,knots=knotsY,Boundary.knots=range(subbase$Y),degre=3))
+offset(log(E)),data=subbase,family=quasipoisson)
predmodel=function(a,y) predict(regsp,newdata=data.frame(A=a,Y=y,E=1))
vZ=outer(vX,vY,predmodel)
ylab="Years (1900-2005)",zlab="Mortality rate (log)")

We now have to use those functions our our small data sample ! That should be fun….

# How old is the oldest person you know?

Last week, we had a discussion with some colleagues about the fact that – in order to prepare for the SOA exams – we did not have time (so far) to mention results on extreme values in our actuarial program. I did gave an introduction in my nonlife actuarial models class, but it was only an introduction, in three hours, in order to illustrate reinsurance pricing. And I told my students that if they wanted to know more about extreme values, they should start a master program in actuarial science and finance, since I will give a course on extremes (and copulas) next winter.

But actually, extreme values are everywhere ! For instance, there is a Prudential TV commercial where has people place large, round stickers on a number line to represent the age of the oldest person they know. This forms some kind of histogram. The message is to have Prudential prepare you to have adequate money for all these years. And actually, anyone can add his or her own sticker at the Prudential website.

Patrick Honner, on his blog (http://mrhonner.com/…), did mention this interesting representation. But this idea is not new, as mentioned in a post, published three years ago. In 1932, Emil Gumbel gave a talk in France on the “âge limite“. And as he wrote it “on peut donc supposer que la distribution de l’âge limite – c’est à dire la probabilité que cet âge ait une valeur donnée – soit Gaussienne“. In 1932 (not aware of Fisher and Tippett work, he thought that the limiting distribution for a maximum would be Gaussian). But a few years after, he read about Fisher’s work, and observed also that “la distribution d’une valeur extrêmes peut être représentée pour un nombre suffisant d’observations par la formule doublement exponentielle, pourvu que la distribution initiale se comporte asymptotiquement comme une exponentielle. La formule devient rigoureuse si la distribution initiale est exponentielle“, as he wrote in 1935. And in 1937, he wrote a paper on “les centennaires” that can also be related to the work of Bortkiewicz on rare events. One should also mention one of the most important paper in extreme value theory, published in 1974 by Balkema and de Haan, on Residual Life Time at Great Age.

Because in this experiment, the question is “How Old is the Oldest Person You Know?“, so it is the distribution of a maximum. And from Fisher-Tippett theorem, if we assume that the age is bounded (and that there exists some finite upper limit), then the limiting distribution for the maxima (or to be more rigorous, a affine transformation of the maxima) should be Weibull distribution. And this is what it looks like

> plot(-x,dweibull(x,2.25,4),type="l",lwd=2)

As an actuary, the only thing I know about demography, is the distribution of the age of death. For instance, consider the following French life table

> alive <- read.table(
+ "https://perso.univ-rennes1.fr/arthur.charpentier/TV8890.csv",
+ sep=";",header=TRUE)$Lx > nb= -diff(alive) > ages=0:110 > plot(ages,nb,type="h") This is the distribution of the age of the death in a given population. Which is not the same as the distribution mentioned above! What we look for is the following: given that someone is alive, what could be the distribution of his-her age ? Actually, if we assume that the yearly number of birth is constant with time (as well as death probability), then we can compute easily to number of people of age $x$ : we take everyone born (exactly) $x$ years ago, and remove all those who died at at $x$, $x-1$, etc. So the function should be > probadeath=nb/sum(nb) > nbx=function(x) 1-sum(probadeath[1:(x+1)]) > surv=Vectorize(nbx)(ages) > distrage=surv/sum(surv) which looks like But this assumption of constant number of birth is not that relevent. And actually, what we need is the distribution of the age within a population… This is a population pyramid, actually. The French one can be downloaded from http://www.insee.fr/fr/ppp/bases-de-donnees/…. > population <- read.table("popinsee2007.csv",sep=";",header=TRUE)$POPTOT07
> ages=0:107
> plot(ages,population/sum(population),type="h")

(the red line being the one obtained previously, using some natality assumptions). Now, let us use this population to generate acquaintances.

> agemax=function(nsim=1000,size=20){
+ agemax=rep(NA,nsim)
+ for(i in 1:nsim){
+ X=sample(ages,prob=population/sum(population),size=size,replace=TRUE)
+ agemax[i]=max(X)}
+ return(agemax)}

Here, we assume that everyone knows 20 other people, randomly chosen in the entire population, then we return the age of the oldest. And we do that for 1,000 people. Here is the distribution, we obtain

> XS=agemax(10000,20)
> plot(table(XS)/length(XS),type="h",xlim=c(0,108))

where the red line is a Weibull distribution (a transformed one, actually, since in extremely value theory, the distance to the upper bound of the distribution has a Weibull density),

> library(MASS)
> fit=fitdistr(108-XS,dweibull,list(shape=1,scale=1))
> lines(ages,dweibull(108-ages,fit$estimate[1],fit$estimate[2]),col="red")

Which is quite close to the distribution obtained in the commercial, don’t you think ? But still, it should be possible to be more accurate, since people should think of their parents, or grandparents. So I guess it could be possible to build a more accurate algorithm, to get something closer to the distribution obtained on the Prudential website. But first, let us wait to have more stickers, more observations… and then I’ll be back to play with it !

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