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 http://freakonometrics.blog.free.fr/public/maths/ETS-NN-3.png, 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 BIC
941.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)
> ETS
ETS(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 BIC
937.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: Y
ARIMA(0,1,0) with drift
 
Call: Arima(x = Y, order = c(0, 1, 0), include.drift = TRUE)
 
Coefficients:
drift
-2.0337
s.e. 1.1904
 
sigma^2 estimated as 138.9: log likelihood = -380.8
AIC = 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: Y
ARIMA(0,1,1) with drift
 
Call: auto.arima(x = Y, allowdrift = TRUE)
 
Coefficients:
ma1 drift
-0.3868 -2.0139
s.e. 0.0970 0.6894
 
sigma^2 estimated as 122.3: log likelihood = -374.65
AIC = 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=2000
x=40
s=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=.035
m=70
h=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 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...



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.