# Statistique de l’assurance STT6705V, partie 12 bis

December 7, 2010
By

(This article was first published on Freakonometrics - Tag - R-english, and kindly contributed to R-bloggers)

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,

To leave a comment for the author, please follow the link and comment on their blog: Freakonometrics - Tag - R-english.

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

If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Tags: , , , , , , , , , , , , , , ,

Comments are closed.

## Sponsors

Contact us if you wish to help support R-bloggers, and place your banner here.

# Never miss an update! Subscribe to R-bloggers to receive e-mails with the latest R posts.(You will not see this message again.)

Click here to close (This popup will not appear again)