**Freakonometrics - Tag - R-english**, and kindly contributed to R-bloggers)

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-A

base=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)

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])

On the other hand, the Lee-Carter (age-year) component is

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)

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

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$kt

K2=LC3$kt[99]+LC3F$kt.f$mean

K=-c(K1,K2)/50

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])

}}

It is also possible to compute residual life expectancy , including dynamics

x=40

E=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)

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=40

t=2000

r=.035

m=70

h=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)

**leave a comment**for the author, please follow the link and comment on his blog:

**Freakonometrics - Tag - R-english**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...