Freakonometrics

To content | To menu | To search

Teaching › actuariat 10/11 STT6705V

Entries feed - Comments feed

Wednesday, January 26 2011

Ecrêter les coûts de sinistres en tarification (2)

Pour poursuivre le précédant billet (en ligne ici), revenons sur le principe même de la tarification: calculer un cout moyen de sinistre. En utilisant la formule des probabilités totales, on peut écrire
http://freakonometrics.blog.free.fr/public/maths/ecret01.png
En conditionnant par le fait d'atteindre, ou pas, le seuil définissant les grands risques, i.e. http://freakonometrics.blog.free.fr/public/maths/ecret2.png - avec http://freakonometrics.blog.free.fr/public/maths/ecret03.png -on obtient alors
http://freakonometrics.blog.free.fr/public/maths/ecret04.png
Maintenant si on raffine un peu, en tarification, on cherche un coût moyen conditionnellement à des variables tarifaire, et donc, on peut écrire - en notant que la probabilité http://freakonometrics.blog.free.fr/public/maths/ecret07.png utilisée importait peu, et que l'on aurait la même relation sous http://freakonometrics.blog.free.fr/public/maths/media08.png
http://freakonometrics.blog.free.fr/public/maths/ecret10.png
Le premier terme correspond aux sinistres 'normaux'. Pour le second, on notera que
http://freakonometrics.blog.free.fr/public/maths/ecret13.png
Autrement dit, on peut être tenté par ne plus distinguer par classe pour le coût moyen des très très gros sinistres. Mais on répartira proportionnellement à la fréquence des gros sinistres sinistres. C'est ce que nous avions fait ici en répartissant parmi tout le monde, en notant aussi que
http://freakonometrics.blog.free.fr/public/maths/ecret14.png
On peut toutefois aller un peu plus loin ici, en modélisant la probabilité conditionnelle
ce qui se fait en faisant - par exemple - une régression logistique,
http://freakonometrics.blog.free.fr/public/maths/ecret-probit-error.png
En reprenant la base précédant, on peut essayer de modéliser la probabilité d'être touché par un très gros sinistre, puis on l'utilise
base2=base
base2$Y[1]=500000+base$Y[1]
seuil=100000
I=base2$Y>seuil
base7=data.frame(base6,I)
regI=glm(I~bs(age,2)+color+sex,data=base7,family=binomial)

On note que seuls les véhicules de 5 ans interviennent ici (la spline ne lisse pas vraiment car on a un seul sinistre, et des données discrètes)
YpI=predict(regI,newdata=data.frame(age=0:20,color="red",sex="H"),
type="response")
reg=glm(Y~bs(age)+color+sex,data=base,family=poisson(link="log"))
Yinf=predict(reg,newdata=data.frame(age=0:20,color="red",sex="H"),
type="response")
Ysup=mean(500000)
Yp=Yinf*(1-YpI)+Ysup*YpI
plot(0:20,Yp,col="red",type="l",ylim=c(1000,2000))
En l'occurrence, ici on va très largement faire contribuer les conducteurs hommes conduisant une voiture rouge de 5 ans, car leur coût moyen devient 11157.862 (soit 10 fois plus que les autres classes de risque)
> Yp
0 1 2 3 4 5 6 7
1732.516 1684.248 1645.529 1614.871 1591.009 11157.862 1559.079 1548.988
Autrement dit, l'écrêtement n'a pas vraiment fonctionné. Toutefois, on peut aussi enlever la fonction spline pour lisser davantage,
base2=base
base2$Y[1]=500000+base$Y[1]
seuil=100000
I=base2$Y>seuil
base7=data.frame(base6,I) regI=glm(I~age+color+sex,data=base7,family=binomial)

et cette fois le cout moyen devient

Ici, sans splines, on a une fonction décroissante de l'âge du véhicule pour la probabilité d'avoir un très gros sinistre: on va pénaliser surtout les véhicules neufs. Si c'est assez peu concluant sur cet exemple, je renvois au polycopié de cours, en ligne ici, pour des mises en œuvre qui marchent davantage. Encore une fois, l'idée est de mutualiser ces très gros sinistres sur des classes voisines de celles qui ont été touchées. Mais ça n'est pas forcément simple...

Ecrêter les coûts de sinistres en tarification (1)

Considérons la base de données (simulée) suivante
set.seed(1)
n=10000
age=sample(0:20,size=n,replace=TRUE)
sex=sample(c("H","F"),size=n,replace=TRUE)
color=sample(c("blue","red","white","green","grey"),
size=n,replace=TRUE)
param=log(1000*exp(-age/100))
Y=rlnorm(n,param,1)
base=data.frame(Y,age,color,sex)
On dipose de coût de sinistres, qui sont fonction de l'âge (du véhicule), mais on dispose aussi de la couleur du véhicule et du sexe du conducteur principal. Supposons que le premier sinistre observé soit incroyablement élevé (comme cela peut arriver en pratique).
base2=base
base2$Y[1]=500000+base$Y[1]
Ce sinistre est associé à un véhicule rouge, de 5 ans, conduit par un homme,
> head(base2)
Y age color sex
1 500768.0434 5 red H
2 838.0116 7 blue F
3 557.3347 12 blue F
Une tarification standard donne
reg=glm(Y~bs(age)+color+sex,data=base2,family=poisson(link="log"))
Yp=predict(reg,newdata=data.frame(age=0:20,color="red",sex="H")
,type="response")
plot(0:20,Yp,col="red",type="l",ylim=c(1000,2000))
Yp=predict(reg,newdata=data.frame(age=0:20,color="red",sex="F")
,type="response")
lines(0:20,Yp,col="red",type="l",lty=2)
Yp=predict(reg,newdata=data.frame(age=0:20,color="blue",sex="H")
,type="response")
lines(0:20,Yp,col="blue",type="l")
Yp=predict(reg,newdata=data.frame(age=0:20,color="blue",sex="F")
,type="response")
lines(0:20,Yp,col="blue",type="l",lty=2)
abline(h=mean(base2$Y),col="grey")
abline(h=mean(Y),col="grey",lwd=.5,lty=2)
soit, graphiquement

où les traits plein sont des hommes, et les traits pointillés sont des femmes, et la couleur désigne la couleur du véhicule. Ici, on observe que le pic de la prédiction de cout moyen est maximal pour les véhicules rouges, pour les hommes, et pour les véhicules de 4-5 ans. Le trait horizontal est la charge moyenne.
L'écrêtement consiste a dire que l'on peut écrêter ce(s) gros sinistre(s). La charge extrême est alors partagée entre certains assurés. Par exemple, le cas extrême consiste à répartir uniformément la surcharge de ce premier sinistre, entre tout le monde.
base3=base
base3$Y=500000/10000+base$Y
reg=glm(Y~bs(age)+color+sex,data=base3,family=poisson(link="log"))
soit graphiquement

Ici, tout le monde participera au financement de ce sinistre extrême. On note que seul l'âge intervient alors véritablement (ce qui est effectivement le cas dans ma simulation). Le sinistre extrême est alors complètement mutualisé.
Toutefois, de la même façon que l'on peut légitimer la segmentation du tarif, on peut se dire que cette mutualisation n'est pas "juste". On peut alors mutualiser uniquement parmi les hommes, par exemple,
base4=base
base4$Y[base4$sex=="H"]=500000/sum(base4$sex=="H")+base$Y
reg=glm(Y~bs(age)+color+sex,data=base4,family=poisson(link="log"))
ce qui donnerait

i.e. un cout moyen toujours plus élevé pour les hommes (qui subissent tous - de manière uniforme) la surprime. Mais on pourrait répartir la charge sur crête entre les hommes conduisant une voiture rouge,
base5=base
base5$Y[(base5$sex=="H")&(base5$color=="red")]=
500000/sum((base5$sex=="H")&(base5$color=="red"))+
base$Y[(base5$sex=="H")&(base5$color=="red")]
reg=glm(Y~bs(age)+color+sex,data=base5,family=poisson(link="log"))
ce qui donnerait ici

un cout moyen plus élevé pour les voitures rouges, et les conducteurs hommes.

Voilà pour les méthodes les plus simples... Mais comme je l'ai fait en cours, on peut aussi utiliser un modèle (binomial) permettant de modéliser les caractéristiques des sinistres les plus importants. à suivre donc...

Sunday, January 2 2011

Statistique de l'assurance STT6705V, questions/réponses

Séance de questions/réponses mercredi après midi avec les élèves de Rennes. Je serais en ligne de 16h30 à 19h30 pour répondre aux éventuelles questions, via ma webcam. Je risque d'avoir un peu de retard au début car je serais en cours juste avant... Petite précision: je me contenterais de répondre aux questions (plutôt de programmation dans le cadre des projets), il ne s'agit pas de refaire le cours.

Tuesday, December 14 2010

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 quantitative point 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 ?

Tuesday, December 7 2010

Statistique de l'assurance STT6705V, partie 12 bis

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,

Thursday, December 2 2010

Statistique de l'assurance STT6705V, les projets

Les projets sont à rendre pour la rentrée des vacances de Noël, pour le lundi 10 janvier (au plus tard) pour les étudiants de l'Université de Montréal, et le lundi 31 janvier (au plus tard) pour les étudiants de l'Université de Rennes1. J'attends trois rapports (un rapport par projet) envoyé par adresse électronique à l'adresse indiquée par courriel (le message devra mettre en copie l'autre membre du binôme pour les projets faits à deux). Merci de préciser les noms de tous les membres du binôme sur la première page du rapport.

  • pour le rapport de tarification (ici pour les bases de données) j'avais donné beaucoup d'indications ici, mais pour faire rapide, pas plus de 15-20 pages. Le projet doit se terminer par plusieurs propositions de primes (basés sur plusieurs modèles) pour un ou plusieurs assurés "types" (indiqués sur le billet mentionné auparavant),
  • pour le rapport de provisionnement (ici pour les bases de données) j'attends un rapport d'une dizaine de pages. Une première partie détaillera sur quel triangle vous allez utiliser les méthodes vues en cours (je vous autorise à ne pas prendre en compte - et donc à exclure - des années trop atypiques, à condition de le justifier un minimum). Une seconde partie présentera l'utilisation de la méthode Mack, avec un best estimate et le mse du montant total de provision. Une troisième partie présentera une ou plusieurs mises en œuvres de méthodes économétriques avec des simulations. Enfin une dernière partie présentera pour tous les modèles considérés un best estimate, mais aussi un quantile à 95% et à 99.5%.
  • pour le rapport de mortalité (ici pour les bases de données) j'attends un rapport d'une dizaine de pages, présentant sommairement dans une première partie la mortalité dans le pays étudié (surface de taux de mortalité et pics apparents); la seconde partie proposera une modélisation du taux de mortalité et une projection pour les 90 ans à venir (le modèle attendu sera celui de Lee-Carter mais tout autre complément sera le bienvenu*); une troisième partie présentera l'évolution de l'espérance de vie entre 1950 et 2000; et enfin, dans une dernière partie, le calcul de la valeur actuelle (à la signature du contrat, en 2000) probable du contrat suivant sera demandé (pour une personne qui avait 45 ans en 2000): la personne veut toucher 10,000 tous les ans à partir de ses 65 ans (à terme échu) tant qu'elle est en vie, avec une garantie décès d'un montant de 50,000 si la personne décède avant 65 ans, et 20,000 si elle décède après 65 ans (strictement). Le taux d'actualisation sera fixé à 3%.
Sinon je ne peux que vous encourager à revoir ce qui a été dit en cours... La règle est que le nombre de pages indiqué est une borne supérieure: merci d'être concis et d'illustrer autant que possible par des graphiques et des tableaux. Les sorties informatiques ne sont pas demandées.

* je mettrais en ligne rapidement un billet détaillant un peu les différentes méthodes de prédiction du coefficient temporel.

Statistique de l'assurance STT6705V, partie 12

The final course (since courses end this week in Montréal) can be watched here and there. The drawings from the course can be downloaded here (including last week's). First, to come back on last week's course , we considered Lee-carter model, i.e.

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)
Here, parameters have the following shape.
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])
The first component is rather similar to what we had before
On the other hand, the Lee-Carter (age-year) component is

and the cohort effect (age-cohort) is

Note that the year-component captures wars and chocks that can affect anyone some given years, and that gains on life expectancy is more a cohort effect.
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 http://freakonometrics.blog.free.fr/public/maths/viekt.png 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)
while it is
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

plot(ANNEE,-LC2$coefficients[199:297]*50)
lines(ANNEE,LC3$kt,col="red")

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])
}}
Once we have that matrix, we simply have to work on diagonal (i.e. cohorts) to calculate anything. For instance, to get the dynamic version of http://freakonometrics.blog.free.fr/public/maths/viehpx.png, the code is (here for someone of age 40 in 2000)
t=2000
x=40
s=seq(0,99-x-1)
MUd=MU[x+1+s,t+s-1898]
Pxt=cumprod(exp(-diag(MUd)))
It is also possible to compute residual life expectancy http://freakonometrics.blog.free.fr/public/maths/viex.png, 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)
Note that this has been computed given the maximum lifetime we had in the dataset (here 99 years old). Hence, we assume here that no one can live more than 100 year (which will impact life expectancy and is a rather strong assumption, especially in 2050 - which might explain the concave shape of the curve on the right part).
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)
(in order to compare the value of such a contract, from 1900 up to 2050).

Monday, November 29 2010

Statistique de l'assurance STT6705V, partie 11

Last course will be uploaded soon (the links will be here and there). The R code considered is given below. First, we had to work a little bit on the datasets,

tabB=read.table("http://perso.univ-rennes1.fr/arthur.charpentier/tabB.csv", 
sep=";",header=FALSE)
ANNEE=tabB[,1]
BASEB=tabB[,seq(2,246,by=2)]
BASEB=as.matrix(BASEB[,1:100])
AGE=0:ncol(BASEB)
tabC=read.table("http://perso.univ-rennes1.fr/arthur.charpentier/tabC.csv",
sep=";",header=FALSE)
an=tabC[,1]
an=an[-c(16,23,43,48,51,53,106)]
BASEC=tabC[,2:101]
BASEC=as.matrix(BASEC[-c(16,23,43,48,51,53,106),])  
BASEB=BASEB[,1:90]
BASEC=BASEC[,1:90]
AGE=AGE[1:90]
MU=as.matrix(log(BASEB/BASEC))
persp(ANNEE,AGE,MU,theta=70,shade=TRUE,
col="green")
library(rgl)
persp3d(ANNEE,AGE,MU,col="light blue")
(this last line is here to play a little bit with the 3d mortality surface). We first used the Lee-Carter function proposed by JPMorgan in LifeMetrics,
source("http://perso.univ-rennes1.fr/arthur.charpentier/fitModels.r")
x=AGE
y=ANNEE
d=BASEB
e=BASEC
w=matrix(1,nrow(d),ncol(e))
LC1=fit701(x,y,e,d,w)
plot(AGE,LC1$beta1)
plot(ANNEE,LC1$kappa2)
plot(AGE,LC1$beta2)
Then we considered nonlinear Poisson regression,
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)
plot(AGE[-1],LC2$coefficients[1]+LC2$coefficients[2:90])
lines(AGE,LC1$beta1,col="blue")
plot(ANNEE,LC2$coefficients[181:279])
plot(ANNEE,LC1$kappa2,col="red")
plot(AGE,LC1$beta2)
As mentioned during the course, this technique is great... but it is sentive to initial values in the optimization procedure. For instance, consider the following loops,
plot(AGE[-1],LC2$coefficients[1]+LC2$coefficients[2:90],
type="l",col="blue",xlab="",ylab="")
for(s in 1:200){
LC2=gnm(D~a+Mult(a,y),offset=log(E),
family=poisson,data=base)
lines(AGE[-1],LC2$coefficients[1]+LC2$coefficients[2:90],col="blue")
}
Here are representation of the first component in the Lee-Carter model,
XXXXHence, the estimation depends on the choice of initial values... even if the shape remains unchanged...
Then, we finnished with Rob Hyndman's package,
library(demography)
base=demogdata(data=t(exp(MU)),pop=t(BASEC),
ages=AGE,years=ANNEE,type="mortality",
label="France",name="Total",lambda=0)
LC3=lca(base)
LC3F=forecast(LC3,100)
plot(LC3$year,LC3$kt,xlim=c(1900,2100),ylim=c(-300,150))
lines(LC3F$year,LC3$kt[99]+LC3F$kt.f$mean)
lines(LC3F$year,LC3$kt[99]+LC3F$kt.f$lower,lty=2)
lines(LC3F$year,LC3$kt[99]+LC3F$kt.f$upper,lty=2)
We concluded with a short discussion about errors of the Lee-Carter model (on French mortality)
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))
RES=residuals(LC2,"pearson")
base$res=RES
plot(base$A,base$res)
couleur=heat.colors(100)
plot(base$A,base$res,col=couleur[base$Y-1898])
plot(base$Y,base$res,col=couleur[base$A+1])
The graphs can be seen below (as a function of time, and a function of ages)
his Wednesday, we will discuss how we can use those estimators (or other ones) in life insurance...

Tuesday, November 23 2010

Statistique de l'assurance STT6705V, partie 11

Pour récupérer les données, et définir un taux de mortalité en fonction de l'année et de l'âge, le code ressemble à ça,

tabB=read.table("http://perso.univ-rennes1.fr/arthur.charpentier/tabB.csv", 
sep=";",header=FALSE)
ANNEE=tabB[,1]
BASEB=tabB[,seq(2,246,by=2)]
BASEB=as.matrix(BASEB[,1:100])
tabC=read.table("http://perso.univ-rennes1.fr/arthur.charpentier/tabC.csv",
sep=";",header=FALSE)
an=tabC[,1]
an=an[-c(16,23,43,48,51,53,106)]
BASEC=tabC[,2:101]
BASEC=as.matrix(BASEC[-c(16,23,43,48,51,53,106),])
AGE=0:(ncol(BASEC)-1)
MU=as.matrix(log(BASEB/BASEC))
persp(ANNEE,AGE,MU,theta=30)

Pour la modélisation des tables de mortalité, nous allons utiliser,

  • le package demography de Rob Hyndman (ici) qui propose par exemple d'implémenter le modèle de Lee & Carter.
  • les fonctions de LifeMetrics de JP Morgan (ici), en particulier dans le code source fitModel.r dont le détail technique traîne en ligne.
  • les fonctions du package gnm de modèles Poissonniens non-linéaires.
Les codes R sont les suivants,
library(demography)
BASEF = demogdata(data=MUF, pop=POPF,
ages=AGE, years=YEAR, type="mortality",
label="France", name="Femmes", lambda=1)
LCF1 = lca(BASEF)
source("http://perso.univ-rennes1.fr/arthur.charpentier/fitModels.r")
LCF2 = fit701(xv=XV,yv=YV,etx=ETF,dtx=DTF,wa=WA)
library(gnm)
LCF3 = gnm(D~factor(Age)+
Mult(Exp(factor(Age)),factor(Year)),
data=base,offset=log(E),family=quasipoisson)

Thursday, November 11 2010

Statistique de l'assurance STT6705V, partie 10, bis

Les liens pour le cours d'hier sont ici et . Sinon, pour ceux qui veulent mettre en œuvre une technique MCMC pour étudier le montant de provision, je mets un code ici (correspondant à celui détaillé dans les notes de cours, inspiré d'un travail fait par Nathalie Balson il y a quelques années maintenant). Pour quelques compléments théoriques, je peux renvoyer ici. Le code utilisé en cours est en ligne ici, et les petits dessins .
Sinon pour les bases de données pour le dernier devoir maison, le code est donné en dessous. Le numéro de la base est le même que pour les bases de tarification que vous aviez choisi... Par exemple, pour le groupe qui avait choisi la base 10,

k=10
loc="http://perso.univ-rennes1.fr/arthur.charpentier/6705V/"
nom=paste(loc,"base",k,"Exposures_1x1.txt",sep="")
base1=read.table(nom,skip=2,header=TRUE)
libelle1=read.table(nom,nrows=1,sep=";")
nom=paste(loc,"base",k,"Deaths_1x1.txt",sep="")
base2=read.table(nom,skip=2,header=TRUE)
(je donnerais plus d'information sur ce qui est attendu dans le rapport par la suite). Pour rappel, il n'y aura pas cours la semaine prochaine....

Tuesday, November 9 2010

Bornhuetter-Ferguson and claims reserving, STT6705V part 10

In several posts on this blog, I have mentioned Chain Ladder, Mack's view of Chain Ladder, and the Overdispersed Poisson model for incremental payments. But there is also one important model used when modeling IBNR: the so-called Bornhuetter-Ferguson model (from that paper). The slides for tomorrow can be downloaded here,


Here is also a link to the slides presented by William Panning in 2006. And here and there, it is possible to find papers (or slides) discussing uncertainty in Bornhuetter-Ferguson's model.
Reserving techniques will probably take us half of the course. We will finally discuss mortality models. The dataset we will use can be downloaded from here (with population data here and death data there). I have converted them in csv files, and the R code to use them is the following,
tabB=read.table("http://perso.univ-rennes1.fr/arthur.charpentier/tabB.csv", 
sep=";",header=FALSE)
ANNEE=tabB[,1]
BASEB=tabB[,seq(2,246,by=2)]
BASEB=BASEB[,1:100]
tabC=read.table("http://perso.univ-rennes1.fr/arthur.charpentier/tabC.csv",
sep=";",header=FALSE)
BASEC=tabC[,2:101]
BASEC=BASEC[-c(16,23,43,48,51,53),]
BASEC=BASEC[1:nrow(BASEB),] AGE=0:99


(I will upload soon the code for claims reserving techniques...)

Thursday, November 4 2010

Splines: opening the (black) box...

Splines in regression is something which looks like a black box (or maybe like some dishes you get when you travel away from home: it tastes good, but you don't what's inside... even if you might have some clues, you never know for sure*). With splines, it is the same: there are knots, then we consider polynomial interpolations on parts between knots, and we make sure that there is no discontinuity (on the prediction, but on the derivative as well).
That sounds nice, but when you look at the output of the regression... you got figures, but you barely see how to interpret them... So let us have a look at the box, and I mean what is inside that box...
So, consider the following dataset, with the following spline regression,
> library(splines)
> K=c(14,20)
> plot(cars)
> reg=lm(dist~bs(speed,knots=c(K),
+ degree=2),data=cars)
> u=seq(4,25,by=.1)
> B=data.frame(speed=u)
> Y=predict(reg,newdata=B)
> lines(u,Y,lwd=2,col="red")
I.e. we have the following (nice) picture,

But if we look at the output of the regression, we get this
> summary(reg)
Call:
lm(formula = dist ~ bs(speed, knots = c(K), degree = 2), data = cars)
Residuals:
    Min      1Q  Median      3Q     Max
-21.848  -8.702  -1.667   6.508  42.514
Coefficients:
                                     Estimate Std. Error t value Pr(>|t|)   
(Intercept)                             5.991      9.692   0.618 0.539596   
bs(speed, knots = c(K), degree = 2)1    8.087     15.278   0.529 0.599174   
bs(speed, knots = c(K), degree = 2)2   45.540     10.735   4.242 0.000109 ***
bs(speed, knots = c(K), degree = 2)3   49.789     15.704   3.170 0.002738 **
bs(speed, knots = c(K), degree = 2)4   95.550     13.651   7.000 1.02e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 15 on 45 degrees of freedom
Multiple R-squared: 0.6888,     Adjusted R-squared: 0.6611
F-statistic: 24.89 on 4 and 45 DF,  p-value: 6.49e-11
So, what can we do with those numbers ?
First, assume know that we consider only one knot (we have to start somewhere), and we consider a b-spline interpolation of degree 1 (i.e. linear by parts).
> K=c(14)
> reg=lm(dist~bs(speed,knots=c(K),degree=1),data=cars)
> summary(reg)
Call:
lm(formula = dist ~ bs(speed, knots = c(K), degree = 1), data = cars)
Coefficients:
                                     Estimate Std. Error t value Pr(>|t|)   
(Intercept)                             3.290      7.746   0.425   0.6730   
bs(speed, knots = c(K), degree = 1)1   31.483      9.670   3.256   0.0021 **
bs(speed, knots = c(K), degree = 1)2   80.525      9.038   8.910 1.16e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 15.41 on 47 degrees of freedom
Multiple R-squared: 0.657,      Adjusted R-squared: 0.6424
F-statistic: 45.01 on 2 and 47 DF,  p-value: 1.203e-11

Recall that b-splines work like that:given knots http://perso.univ-rennes1.fr/arthur.charpentier/latex/spl01.png (we define splines on the unit interval), then
http://perso.univ-rennes1.fr/arthur.charpentier/latex/spl02.png
for all http://perso.univ-rennes1.fr/arthur.charpentier/latex/spl03.png, while
http://perso.univ-rennes1.fr/arthur.charpentier/latex/spl04.png
for all http://perso.univ-rennes1.fr/arthur.charpentier/latex/spl05.png. The code is something like that
> B=function(x,j,n,K){
+ b=0
+ a1=a2=0
+ if(((K[j+n+1]>K[j+1])&(j+n<=length(K))&(n>0))==TRUE){a2=(K[j+n+1]-x)/
+          (K[j+n+1]-K[j+1])*B(x,j+1,n-1,K) }
+ if(((K[j+n]>K[j])&(n>0))==TRUE){a1=(x-K[j])/
+          (K[j+n]-K[j])*B(x,j,n-1,K)}
+ if(n==0){ b=((x>K[j])&(x<=K[j+1]))*1 }
+ if(n>0){  b=a1+a2}
+ return(b)
+ }
So, for instance, for splines of degree 1, we have
> u=seq(0,1,by=.01)
> plot(u,B(u,1,1,c(0,.4,1,1)),lwd=2,col="red",type="l",ylim=c(0,1))
> lines(u,B(u,2,1,c(0,.4,1,1)),lwd=2,col="blue")
> abline(v=c(0,.4,1),lty=2)

and for splines of degree 2, the basis is
> u=seq(0,1,by=.01)
> plot(u,B(u,1,2,c(0,0,.4,1,1,1)),lwd=2,col="red",type="l",ylim=c(0,1))
> lines(u,B(u,2,2,c(0,0,.4,1,1,1)),lwd=2,col="blue")
> lines(u,B(u,3,2,c(0,0,.4,1,1,1)),lwd=2,col="green")
> abline(v=c(0,.4,1),lty=2)

...etc. Note that I need to duplicate sometimes starting and end points (but it should be possible to fix it in the function).
So, how do we use that, here ? Actually, there are two steps:
  1. we get from http://perso.univ-rennes1.fr/arthur.charpentier/latex/spl06.png to the unit interval (using a simple affine transformation)
  2. we consider http://perso.univ-rennes1.fr/arthur.charpentier/latex/sp07.png
Here, based on the graph above (with the basis function), note that we can use
> u0=seq(0,1,by=.01)
> v=reg$coefficients[2]*u0+reg$coefficients[1]
> x1=seq(min(cars$speed),K,length=length(u0))
> lines(x1,v,col="green",lwd=2)
> u0=seq(0,1,by=.01)
> v=(reg$coefficients[3]-reg$coefficients[2])*u0+ reg$coefficients[1]+reg$coefficients[2]
> x2=seq(K,max(cars$speed),length=length(u0))
> lines(x2,v,col="blue",lwd=2)

which gives us exactly the graph we obtained previously. But we can also consider
> plot(cars)
> reg=lm(dist~bs(speed,knots=c(K),degree=1),data=cars)
> k=(K-min(cars$speed))/(max(cars$speed)-min(cars$speed))
> u0=seq(0,1,by=.01)
> v=reg$coefficients[1]+
+   reg$coefficients[2]*B(u0,1,1,c(0,k,1,1))+
+   reg$coefficients[3]*B(u0,2,1,c(0,k,1,1))
> lines(min(cars$speed)+u0*(max(cars$speed)-min(cars$speed)),
+ v,col="purple",lwd=2)
> abline(v=K,lty=2,col="red")

So, we should be able to try with two knots (but we keep it linear, so far). Hence
> K=c(14,20)
> plot(cars)
> reg=lm(dist~bs(speed,knots=c(K),degree=1),data=cars)
> u=seq(4,25,by=.1)
> B=data.frame(speed=u)
> Y=predict(reg,newdata=B)
> lines(u,Y,lwd=2,col="red")
> abline(v=K,lty=2,col="red")
First, we can plot our basis functions, with two knots,
> u=seq(0,1,by=.01)
> plot(u,B(u,1,1,c(0,.4,.7,1)),lwd=2,col="red",type="l",ylim=c(0,1))
> lines(u,B(u,2,1,c(0,.4,.7,1,1)),lwd=2,col="blue")
> lines(u,B(u,3,1,c(0,.4,.7,1,1)),lwd=2,col="green")
> abline(v=c(0,.4,.7,1),lty=2)

so we can use those functions here, as we did before,
> plot(cars)
> reg=lm(dist~bs(speed,knots=c(K),degree=1),data=cars)
> k=(K-min(cars$speed))/(max(cars$speed)-min(cars$speed))
> u0=seq(0,1,by=.01)
> v=reg$coefficients[1]+
+   reg$coefficients[2]*B(u0,1,1,c(0,k,1,1))+
+   reg$coefficients[3]*B(u0,2,1,c(0,k,1,1))+
+   reg$coefficients[4]*B(u0,3,1,c(0,k,1,1))
> lines(min(cars$speed)+u0*(max(cars$speed)-min(cars$speed)),
+ v,col="red",lwd=2)
> abline(v=K,lty=2,col="red")

Great, it looks promising.... Let us look finally at the case we have two knots, and some quadratic splines. Here, with two knots, the basis is
> u=seq(0,1,by=.01)
> plot(u,B(u,1,2,c(0,0,.4,.7,1,1,1)),lwd=2,col="red",type="l",ylim=c(0,1))
> lines(u,B(u,2,2,c(0,0,.4,.7,1,1,1)),lwd=2,col="blue")
> lines(u,B(u,3,2,c(0,0,.4,.7,1,1,1)),lwd=2,col="green")
> lines(u,B(u,4,2,c(0,0,.4,.7,1,1,1)),lwd=2,col="orange")
> abline(v=c(0,.4,.7,1),lty=2)

so if we just rewrite our previous function, we have
> plot(cars)
> reg=lm(dist~bs(speed,knots=c(K),degree=2),data=cars)
> k=(K-min(cars$speed))/(max(cars$speed)-min(cars$speed))
> u0=seq(0,1,by=.01)
> v=reg$coefficients[1]+
+   reg$coefficients[2]*B(u0,1,2,c(0,0,k,1,1,1))+
+   reg$coefficients[3]*B(u0,2,2,c(0,0,k,1,1,1))+
+   reg$coefficients[4]*B(u0,3,2,c(0,0,k,1,1,1))+
+   reg$coefficients[5]*B(u0,4,2,c(0,0,k,1,1,1))
> lines(min(cars$speed)+u0*(max(cars$speed)-min(cars$speed)),
+ v,col="purple",lwd=2)
> abline(v=K,lty=2,col="red")

And just one final comment: how do I optimally choose my knots ?... A simple idea can be to consider all possible knots, and consider the model which gives us the residuals with the smaller variance,
> vk=seq(.05,.95,by=.05)
> SSR=matrix(NA,length(vk))
> for(i in 1:(length(vk))){
+ k=vk[i]
+ K=min(cars$speed)+k*(max(cars$speed)-min(cars$speed))
+ reg=lm(dist~bs(speed,knots=c(K),degree=2),data=cars)
+ SSR[i]=sum(residuals(reg)^2)
+ }
> plot(vk,SSR,type="b",col="blue")

Here, the best model is obtained when we split 3/4-1/4...

* and I know what I am talking about: I have been living in China for more than a year....

Statistique de l'assurance STT6705V, partie 9, ter

Suite à des petits soucis légaux, je ne peux pas proposer de triangles par trimestres, ou sur plusieurs régions, comme promis. Les données pour le second devoir-maison sont en ligne ici. Les données ne sont pas très jolies, car il manque la partie supérieure gauche des triangles... mais en adaptant les techniques vues en cours, c'est largement faisable. Il y a les triangles de paiements, mais aussi des estimations de charge (par les gestionnaires de sinistres, i.e. charge dossier/dossier).

Les binômes restent les mêmes que pour le premier devoir, il s'agit de me proposer un montant de provision, et une estimation du quantile à 95% du montant de provisions. En cas de problème, vous avez mon adresse électronique.... Et je reviendrais ultérieurement sur la forme de ce que j'attends. La correspondance est simple, et fonction des bases utilisées en tarification

  • bases 1 et 2 : triangle onglet 1
  • bases 3 et 4 : triangle onglet 2
  • bases 5 et 7 : triangle onglet 3
  • bases 8 et 9 : triangle onglet 4
  • bases 10 et 11 : triangle onglet 5
  • bases 17 et 20 : triangle onglet 6
  • bases 21 et 22 : triangle onglet 7
  • bases 23 : triangle onglet 8

Attention, les triangles sont présentés ici à l'envers par rapport au cours (les années de survenance sont en colonne).

Wednesday, November 3 2010

Statistique de l'assurance STT6705V, partie 9, bis

Les liens pour le cours de ce matin sont en ligne ici (en un bloc de deux heures). Sinon, comme cela semblait en intéresser certains, j'ai mis en ligne quelques versions des dessins que j'ai pu  faire au tableau... cours4, cours5, cours6, cours7, cours8 et cours9. Ca ressemble aux dessins de mon fils, donc je ne les mets pas pour une quelconque valeur artistique, mais essentiellement parce que ces derniers temps, j'ai davantage utilisé le tableau blanc que les transparents...

My residuals look weird... aren't they ?

Since I got the same question twice, let us look at it quickly....  Some students show me a graph (from a Poisson regression) which looks like that,

and they asked "isn't it weird ?", i.e."residuals are null or positive... this is not what we should have, right ?". Actually, residuals are always centered in a glm regression (if you keep the intercept). Let us generate a dataset,
> n=500
> set.seed(1)
> X=runif(n)
> E=rnorm(n)*.3
> XB=-5+2*X+E
> lambda=exp(XB)
> Y=rpois(n,lambda)
> reg=glm(Y~X,family=poisson)
> plot(predict(reg),residuals(reg,"pearson"))
If we look carefully, residuals in the lower part, they are strictly negative... And since, on average, we should have zero, then we are very close to 0.
> mean(residuals(reg,"pearson"))
[1] -0.003381202

Actually, the true value is either
> u=seq(-10,2,by=.1)
> lines(u,(0-exp(u))/sqrt(exp(u)),lwd=.5,col="red")
or
> lines(u,(1-exp(u))/sqrt(exp(u)),lwd=.5,col="blue")

So what happens ? Here, the value of the parameter (the expected value of our Poisson distribution) is very small
> mean(Y)
[1] 0.026

A crude estimator for the glm would be to consider a regression only on the intercept. In that case, the error would be either 0-0.026, or 1-0.026, i.e. for 97.5%, we have -0.026 (which is small, close to 0, but negative), while for 2.5% we have 0.974 (which is huge compared to the previous value). Well, on that graph, we plot Pearson's residuals, i..e. we divide the difference by the expected value... so we have those two curves....
Actually, there are only two clouds since the probability to have 1 (when lambda is close to 2.5%) is small, and extremely small to exceed 2...
On the graph below, the average value of Y increases (and then decreases): initially, almost all the points are on the red curve (i.e. we observe 0, and predict something extremely small), then we have more and more points on the blue curve (for more and more individuals we have 1), then we start to see some on the green curve (i.e. 2), etc. Note that areas where we observe clouds remain unchanged (i.e. the colored curves).
> P=seq(-5,-2,by=.05)
> P=c(P,rev(P))
> for(i in 1:length(P)){
+ set.seed(1)
+ XB=P[i]+2*X+E
+ lambda=exp(XB)
+ Y=rpois(n,lambda)
+ reg=glm(Y~X,family=poisson)
+ u=seq(-10,2,by=.1)
+ plot(predict(reg),residuals(reg,"pearson"),ylim=c(-1,5),xlim=c(-6,0))
+ abline(h=0)
+ lines(u,(0-exp(u))/sqrt(exp(u)),lwd=.5,col="red")
+ lines(u,(1-exp(u))/sqrt(exp(u)),lwd=.5,col="blue")
+ lines(u,(2-exp(u))/sqrt(exp(u)),lwd=.5,col="green")
+ rlines(u,(3-exp(u))/sqrt(exp(u)),lwd=.5,col="orange")
+ lines(u,(3-exp(u))/sqrt(exp(u)),lwd=.5,col="orange")
+ lines(u,(4-exp(u))/sqrt(exp(u)),lwd=.5,col="purple")
+ points(predict(reg),residuals(reg,"pearson"))
+ }


So the graph is perfectly normal since we have a discrete distribution (while we predict something continuous).

- page 1 of 3