# Overdispersed Poisson et bootstrap

Pour le dernier cours sur les méthodes de provisionnement, on s’est arrête aux méthodes par simulation. Reprenons là où on en était resté au dernier billet où on avait vu qu’en faisant une régression de Poisson sur les incréments, on obtenait exactement le même montant que la méthode Chain Ladder,

> Y
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 3209 1163   39   17    7   21
[2,] 3367 1292   37   24   10   NA
[3,] 3871 1474   53   22   NA   NA
[4,] 4239 1678  103   NA   NA   NA
[5,] 4929 1865   NA   NA   NA   NA
[6,] 5217   NA   NA   NA   NA   NA

> y=as.vector(as.matrix(Y))
> base=data.frame(y,ai=rep(2000:2005,n),bj=rep(0:(n-1),each=n))
> reg2=glm(y~as.factor(ai)+as.factor(bj),data=base,family=poisson)
> summary(reg2)

Call:
glm(formula = y ~ as.factor(ai) + as.factor(bj), family = poisson,
data = base)

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)        8.05697    0.01551 519.426  < 2e-16 ***
as.factor(ai)2001  0.06440    0.02090   3.081  0.00206 **
as.factor(ai)2002  0.20242    0.02025   9.995  < 2e-16 ***
as.factor(ai)2003  0.31175    0.01980  15.744  < 2e-16 ***
as.factor(ai)2004  0.44407    0.01933  22.971  < 2e-16 ***
as.factor(ai)2005  0.50271    0.02079  24.179  < 2e-16 ***
as.factor(bj)1    -0.96513    0.01359 -70.994  < 2e-16 ***
as.factor(bj)2    -4.14853    0.06613 -62.729  < 2e-16 ***
as.factor(bj)3    -5.10499    0.12632 -40.413  < 2e-16 ***
as.factor(bj)4    -5.94962    0.24279 -24.505  < 2e-16 ***
as.factor(bj)5    -5.01244    0.21877 -22.912  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 46695.269  on 20  degrees of freedom
Residual deviance:    30.214  on 10  degrees of freedom
(15 observations deleted due to missingness)
AIC: 209.52

Number of Fisher Scoring iterations: 4

> base$py2=predict(reg2,newdata=base,type="response") > round(matrix(base$py2,n,n),1)
[,1]   [,2] [,3] [,4] [,5] [,6]
[1,] 3155.7 1202.1 49.8 19.1  8.2 21.0
[2,] 3365.6 1282.1 53.1 20.4  8.8 22.4
[3,] 3863.7 1471.8 61.0 23.4 10.1 25.7
[4,] 4310.1 1641.9 68.0 26.1 11.2 28.7
[5,] 4919.9 1874.1 77.7 29.8 12.8 32.7
[6,] 5217.0 1987.3 82.4 31.6 13.6 34.7

> sum(base$py2[is.na(base$y)])
[1] 2426.985

Le plus intéressant est peut être de noter que la loi de Poisson présente probablement trop peu de variance…

> reg2b=glm(y~as.factor(ai)+as.factor(bj),data=base,family=quasipoisson)
> summary(reg2)

Call:
glm(formula = y ~ as.factor(ai) + as.factor(bj), family = quasipoisson,
data = base)

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)        8.05697    0.02769 290.995  < 2e-16 ***
as.factor(ai)2001  0.06440    0.03731   1.726 0.115054
as.factor(ai)2002  0.20242    0.03615   5.599 0.000228 ***
as.factor(ai)2003  0.31175    0.03535   8.820 4.96e-06 ***
as.factor(ai)2004  0.44407    0.03451  12.869 1.51e-07 ***
as.factor(ai)2005  0.50271    0.03711  13.546 9.28e-08 ***
as.factor(bj)1    -0.96513    0.02427 -39.772 2.41e-12 ***
as.factor(bj)2    -4.14853    0.11805 -35.142 8.26e-12 ***
as.factor(bj)3    -5.10499    0.22548 -22.641 6.36e-10 ***
as.factor(bj)4    -5.94962    0.43338 -13.728 8.17e-08 ***
as.factor(bj)5    -5.01244    0.39050 -12.836 1.55e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasipoisson family taken to be 3.18623)

Null deviance: 46695.269  on 20  degrees of freedom
Residual deviance:    30.214  on 10  degrees of freedom
(15 observations deleted due to missingness)
AIC: NA

Number of Fisher Scoring iterations: 4

Mais on en reparlera dans un instant. Ensuite, on avait commencé à regarder erreurs commises, sur la partie supérieure du triangle. Classiquement, par construction, les résidus de Pearson sont de la forme

$\varepsilon_i=\frac{Y_i-\widehat{Y}_i}{\sqrt{\text{Var}(Y_i)}}$

On avait vu dans le cours de tarification que la variance au dénominateur pouvait être remplacé par le prévision, puisque dans un modèle de Poisson, l’espérance et la variance sont identiques. Donc on considérait

$\varepsilon_i=\frac{Y_i-\widehat{Y}_i}{\sqrt{\widehat{Y}_i}}$

> base$erreur=(base$y-base$py2)/sqrt(base$py2)
> round(matrix(base$erreur,n,n),1) [,1] [,2] [,3] [,4] [,5] [,6] [1,] 0.9 -1.1 -1.5 -0.5 -0.4 0 [2,] 0.0 0.3 -2.2 0.8 0.4 NA [3,] 0.1 0.1 -1.0 -0.3 NA NA [4,] -1.1 0.9 4.2 NA NA NA [5,] 0.1 -0.2 NA NA NA NA [6,] 0.0 NA NA NA NA NA Le soucis est que si $\widehat{Y}_i$ est – asymptotiquement – un bon estimateur pour $\text{Var}(Y_i)$, ce n’est pas le cas à distance finie, car on a alors un estimateur biaisé pour la variance, et donc la variance des résidus n’a que peu de chances d’être de variance unitaire. Aussi, il convient de corriger l’estimateur de la variance, et on pose alors $\varepsilon_i=\sqrt{\frac{n}{n-k}}\cdot\frac{Y_i-\widehat{Y}_i}{\sqrt{\widehat{Y}_i}}$ qui sont alors les résidus de Pearson tel qu’on doit les utiliser. > E=base$erreur[is.na(base$y)==FALSE]*sqrt(21/(21-11)) > E [1] 1.374976e+00 3.485024e-02 1.693203e-01 -1.569329e+00 1.887862e-01 [6] -1.459787e-13 -1.634646e+00 4.018940e-01 8.216186e-02 1.292578e+00 [11] -3.058764e-01 -2.221573e+00 -3.207593e+00 -1.484151e+00 6.140566e+00 [16] -7.100321e-01 1.149049e+00 -4.307387e-01 -6.196386e-01 6.000048e-01 [21] -8.987734e-15 > boxplot(E,horizontal=TRUE) En rééchantillonnant dans ces résidus, on peut alors générer un pseudo triangle. Pour des raisons de simplicités, on va générer un peu rectangle, et se restreindre à la partie supérieure, > Eb=sample(E,size=36,replace=TRUE) > Yb=base$py2+Eb*sqrt(base$py2) > Ybna=Yb > Ybna[is.na(base$y)]=NA
> Tb=matrix(Ybna,n,n)
> round(matrix(Tb,n,n),1)
[,1]   [,2] [,3] [,4] [,5] [,6]
[1,] 3115.8 1145.4 58.9 46.0  6.4 26.9
[2,] 3179.5 1323.2 54.5 21.3 12.2   NA
[3,] 4245.4 1448.1 61.0  7.9   NA   NA
[4,] 4312.4 1581.7 68.7   NA   NA   NA
[5,] 4948.1 1923.9   NA   NA   NA   NA
[6,] 4985.3     NA   NA   NA   NA   NA

Cette fois, on a un nouveau triangle ! on va alors pouvoir faire plusieurs choses,

1. compléter le triangle para la méthode Chain Ladder, c’est à dire calculer les montants moyens que l’on pense payer dans les années futures
2. générer des scénarios de paiements pour les années futurs, en générant des paiements suivant des lois de Poisson (centrées sur les montants moyens que l’on vient de calculer)
3. générer des scénarios de paiements avec des lois présentant plus de variance que la loi de Poisson. Idéalement, on voudrait simuler des lois qusi-Poisson, mais ce ne sont pas de vraies lois. Par contre, on peut se rappeler que dans ce cas, la loi Gamma devrait donner une bonne approximation.

Pour ce dernier point, on va utiliser le code suivant, pour générer des quasi lois,

> rqpois = function(n, lambda, phi, roundvalue = TRUE) {
+ b = phi
+ a = lambda/phi
+ r = rgamma(n, shape = a, scale = b)
+ if(roundvalue){r=round(r)}
+ return(r)
+ }

Je renvois aux diverses notes de cours pour plus de détails sur la justification, ou à un vieux billet. On va alors faire une petite fonction, qui va soit somme les paiements moyens futurs, soit sommer des générations de scénarios de paiements, à partir d’un triangle,

> CL=function(Tri){
+ y=as.vector(as.matrix(Tri))
+ base=data.frame(y,ai=rep(2000:2005,n),bj=rep(0:(n-1),each=n))
+ reg=glm(y~as.factor(ai)+as.factor(bj),data=base,family=quasipoisson)
+ py2=predict(reg,newdata=base,type="response")
+ pys=rpois(36,py2)
+ pysq=rqpois(36,py2,phi=3.18623)
+ return(list(
+ cl=sum(py2[is.na(base$y)]), + sc=sum(pys[is.na(base$y)]),
+ scq=sum(pysq[is.na(base$y)]))) + } Reste alors à générer des paquets de triangles. Toutefois, il est possible de générer des triangles avec des incréments négatifs. Pour faire simple, on mettra des valeurs nulles quand on a un paiement négatif. L’impact sur les quantiles sera alors (a priori) négligeable. > for(s in 1:1000){ + Eb=sample(E,size=36,replace=TRUE)*sqrt(21/(21-11)) + Yb=base$py2+Eb*sqrt(base$py2) + Yb=pmax(Yb,0) + scY=rpois(36,Yb) + Ybna=Yb + Ybna[is.na(base$y)]=NA
+ Tb=matrix(Ybna,6,6)
+ C=CL(Tb)
+ VCL[s]=C$cl + VR[s]=C$sc
+ VRq[s]=C$scq + } Si on regarde la distribution du best estimate, on obtient > hist(VCL,proba=TRUE,col="light blue",border="white",ylim=c(0,0.003)) > boxplot(VCL,horizontal=TRUE,at=.0023,boxwex = 0.0006,add=TRUE,col="light green") > D=density(VCL) > lines(D) > I=which(D$x<=quantile(VCL,.05))
> polygon(c(D$x[I],rev(D$x[I])),c(D$y[I],rep(0,length(I))),col="blue",border=NA) > I=which(D$x>=quantile(VCL,.95))
> polygon(c(D$x[I],rev(D$x[I])),c(D$y[I],rep(0,length(I))),col="blue",border=NA) Mais on peut aussi visualiser des scénarios basés sur des lois de Poisson (équidispersé) ou des scénarios de lois quasiPoisson (surdispersées), ci-dessous Dans ce dernier cas, on peut en déduire le quantile à 99% des paiements à venir. > quantile(VRq,.99) 99% 2855.01 Il faut donc augmenter le montant de provisions de l’ordre 15% pour s’assurer que la compagnie pourra satisfaire ses engagements dans 99% des cas, > quantile(VRq,.99)-2426.985 99% 428.025 # Reserving with negative increments in triangles A few months ago, I did published a post on negative values in triangles, and how to deal with them, when using a Poisson regression (the post was published in French). The idea was to use a translation technique: 1. Fit a model not on $Y_i$‘s but on $Y_i^{(k)}=Y_i+k$, for some $k\geq 0$, 2. Use that model to make predictions, and then translate those predictions, $\widehat{Y}_i^{(k)}-k$ This is what was done to get the following graph, where a Poisson regression was fitted. Black points are $Y_i$‘s while blue points are $\widehat{Y}_i^{(k)}$‘s, for some $k\geq 0$. We fit a model to get the blue prediction, and then translate it to get the red prediction (on the $Y_i$‘s). In this example, there were no negative values, but it is possible to use it get a better understanding on the impact of this technique. The prediction, here, is the red line. And clearly, the value of $k$ has an impact on the prediction (since we do not consider, here, a linear model: with a linear model, translating has not impact at all, except on the intercept). The alternative mentioned in the previous post was to use this technique on several $k$‘s, and them interpolate 1. For a given $k$, fit a model not on $Y_i$‘s but on $Y_i^{(k)}=Y_i+k$, use that model to make predictions, and then translate those predictions, $\widehat{Y}_i^{(k)}-k$. 2. Do it for several $k$‘s. 3. Use it to extrapolate when $k$ is $0$ (which is the case we are interested in). In the context of loss reserving, the idea is extremely simple. Consider a triangle with incremental payments > source("https://perso.univ-rennes1.fr/arthur.charpentier/bases.R") > Y=T=PAID > n=ncol(T) > Y[,2:n]=T[,2:n]-T[,1:(n-1)] > Y [,1] [,2] [,3] [,4] [,5] [,6] [1,] 3209 1163 39 17 7 21 [2,] 3367 1292 37 24 10 NA [3,] 3871 1474 53 22 NA NA [4,] 4239 1678 103 NA NA NA [5,] 4929 1865 NA NA NA NA [6,] 5217 NA NA NA NA NA Now, we do not have negative values, here, but we can still see is translation techniques can be used. The benchmark is the Poisson regression, since we can run it : > y=as.vector(as.matrix(Y)) > base=data.frame(y,ai=rep(2000:2005,n),bj=rep(0:(n-1),each=n)) > reg=glm(y~as.factor(ai)+as.factor(bj),data=base,family=poisson) Here, the amount is reserve is the sum of predicted values in the lower part of the triangle, > py=predict(reg,newdata=base,type="response") > sum(py[is.na(base$y)])
[1] 2426.985

which is exactly Chain Ladder’s estimate.

Now, let us use a translation technique to compute the amount of reserves. The code will be

> decal=function(k){
+ reg=glm(y+k~as.factor(ai)+as.factor(bj),data=base,family=poisson)
+ py=predict(reg,newdata=base,type="response")
+ return(sum(py[is.na(base$y)]-k)) For instance, if we translate of +5, we would get > decal(5) [1] 2454.713 while a translation of +10 would return > decal(10) [1] 2482.29 Clearly, translations do have an impact on the estimation. Here, just to check, if we do not translate, we do have Chain Ladder’s estimate, > decal(0) [1] 2426.985 The idea mentioned in the previous post was to try several translations, and then extrapolate, to get the value in 0. Here, translations will give the following estimates > K=10:20 > (V=Vectorize(decal)(K)) [1] 2482.290 2487.788 2493.279 2498.765 2504.245 2509.719 2515.187 2520.649 [9] 2526.106 2531.557 2537.001 We can plot those values, and run a regression > plot(K,V,xlim=c(0,20),ylim=c(2425,2540)) > abline(h=decal(0),col="red",lty=2) the dotted horizontal line is Chain Ladder. Now, let us extrapolate > b=data.frame(K=K,D=V) > rk=lm(D~K,data=b) > predict(rk,newdata=data.frame(K=0)) 1 2427.623 On has to admit that it is not that bad. But yesterday evening, Karim asked me why I did use a linear regression, for my extrapolation. And to be honest, I do not know. I mean, the only answer might be that points are almost on a straight line. So the first time I saw it, I was exited, and I ran a linear regression. Now, let us see if we can do better. Because here, we do use a translation of +10 or +20 (which might be rather small). What if we use much larger values ? (because we might have large negative incremental values). With the following code, we try, each time 11 consecutive values, the smallest one going from 0 to 50, > hausse=1:50; res=rep(NA,50) > for(k in hausse){ + VK=k:(10+k) + b=data.frame(K=VK,D=Vectorize(decal)(VK)) + rk=lm(D~K,data=b) + res[k]=predict(rk,newdata=data.frame(K=0)) + } > plot(hausse,res,type="l",col="red",ylim=c(2422,2440)) > abline(rk,col="blue") Here, we compute reserves when extrapolations were done after 11 translations, from $k$ to $k+10$. With different values of $k$. The case where $k$ is ten was the one mentioned above, > res[hausse==10] [1] 2427.623 Actually, it might also be possible to consider not 11 translations, but 26, from $k$ to $k+25$. Here, we get > hausse=1:50; res=rep(NA,50) > for(k in hausse){ + VK=k:(25+k) + b=data.frame(K=VK,D=Vectorize(decal)(VK)) + rk=lm(D~K,data=b) + res[k]=predict(rk,newdata=data.frame(K=0)) + } > lines(hausse,res,type="l",col="blue",lty=2) We now have the dotted line Here, it is getting worst. So let us keep here 11 translations. Perhaps, we can try something different. For instance a Poisson regression, with a log like (i.e. we consider an exponential extrapolation), > hausse=1:50; res=rep(NA,50) > for(k in hausse){ + VK=k:(10+k) + b=data.frame(K=VK,D=Vectorize(decal)(VK)) + rk=glm(D~K,data=b,family=poisson) + res[k]=predict(rk,newdata=data.frame(K=0),type="response") + } > lines(hausse,res,type="l",col="purple") The purple line will be a Poisson model, with a log link. Perhaps we can try another link function, like a quadratic one > hausse=1:50; res=rep(NA,50) > for(k in hausse){ + VK=k:(10+k) + b=data.frame(K=VK,D=Vectorize(decal)(VK)) + rk=glm(D~K,data=b,family=poisson(link= + power(lambda = 2))) + res[k]=predict(rk,newdata=data.frame(K=0),type="response") + } > lines(hausse,res,type="l",col="orange") That would be the orange line, Here, we need a link function between identity (the linear model, the blue line) and the quadratic one (the orange one), for instance a power function 3/2, > hausse=1:50; res=rep(NA,50) > for(k in hausse){ + VK=k:(10+k) + b=data.frame(K=VK,D=Vectorize(decal)(VK)) + rk=glm(D~K,data=b,family=poisson(link= + power(lambda = 1.5))) + res[k]=predict(rk,newdata=data.frame(K=0),type="response") + } > lines(hausse,res,type="l",col="green") Here, it looks like we can use that model for any kind of translation, from +10 till +50, even +100 ! But I do not have any intuition about the use of this power function… # Examen intra, éléments de correction L’énoncé de l’examen intra est en pdf ici et comme annoncé par courriel, la correction de l’intra est dans le pdf en ligne. Comme personne ne semble en désaccord avec les réponses proposées, les notes seront mises en ligne très bientôt. Concertant les questions 18 et 19 quelques compléments d’explications (que je n’avais pas tapé dans le pdf). On avait vu que l’estimateur du maximum de vraisemblance pour une régression de Poisson était asymptotiquement Gaussien, $\widehat{\boldsymbol{\beta}}_{P}\sim\mathcal{N}(\boldsymbol{\beta},V_\infty(\widehat{\boldsymbol{\beta}}_{P}))$ (asymptotiquement) avec $V_\infty(\widehat{\boldsymbol{\beta}}_{P})=\left(\sum_{i=1}^n \widehat Y_i \boldsymbol{X}_i\boldsymbol{X}_i'\right)^{-1}$ Quand on a une régression de type binomiale négative, si on note de manière très générale $\omega_i=\text{Var}(Y_i|\boldsymbol{X}_i)$ (on avait vu en cours qu’il existait plusieurs spécifications possibles pour cette variance conditionnelle). Dans ce cas, $\widehat{\boldsymbol{\beta}}_{BN}\sim\mathcal{N}(\boldsymbol{\beta},V_\infty(\widehat{\boldsymbol{\beta}}_{BN}))$ avec $V_\infty(\widehat{\boldsymbol{\beta}}_{P})=\left(\sum_{i=1}^n \widehat Y_i \boldsymbol{X}_i\boldsymbol{X}_i'\right)^{-1}\left[\sum_{i=1}^n \omega_i \boldsymbol{X}_i\boldsymbol{X}_i\right]\left(\sum_{i=1}^n \widehat Y_i \boldsymbol{X}_i\boldsymbol{X}_i'\right)^{-1}$ Bref, tout dépend fondamentalement de la spécification de la variance conditionnelle. Sous R, c’est la régression binomiale négative de type 1 qui est considérée, i.e. $\omega_i=\text{Var}(Y_i|\boldsymbol{X}_i)=\phi\cdot \mathbb{E}(Y_i|\boldsymbol{X}_i)=\phi \cdot \widehat{Y}_i$ On toujours une relation de la forme $\widehat{\boldsymbol{\beta}}_{QP}\sim\mathcal{N}(\boldsymbol{\beta},V_\infty(\widehat{\boldsymbol{\beta}}_{QP}))$ avec (en simplifiant un peu) $V_\infty(\widehat{\boldsymbol{\beta}}_{QP})=\phi\cdot\left(\sum_{i=1}^n \widehat Y_i \boldsymbol{X}_i\boldsymbol{X}_i'\right)^{-1}$ aussi, on a $V_\infty(\widehat{\boldsymbol{\beta}}_{QP})=\phi\cdot V_\infty(\widehat{\boldsymbol{\beta}}_{P})$ Mais comme annoncé en cours, des points étaient données pour ceux qui se contentaient d’affirmer que la variance des estimateurs était plus grande s’il y avait sur-dispersion. # Multiple (smoothed) regression and portfolio exposure Wednesday, in class, we’ve seen how to visualize a multiple regression model (with two continuous explanatory variables). Here, the goal is to predict the average cost of an insurance claim, using some covariates, e.g. the age of the driver, and the age of the car (recall that losses here are liability losses). The prediction obtained from a (standard) generalized linear model, with a log-link > reg1=glm(cout~ageconducteur+agevehicule,data=base,family=Gamma(link="log")) The code to visualize the predicted average cost is the following: first, we have to compute predictions for specific values, > pred=function(x,y){ + predict(reg,newdata=data.frame(ageconducteur=x, + agevehicule=y),type="response") Then, we use this function to compute values on a grid, > X=seq(20,80,by=5) > Y=0:20 > Z=outer(X,Y,p) > image(X,Y,Z,col=rev(heat.colors(101))) > contour(X,Y,Z,add=TRUE, + levels=c(1400,1800,2000,2200,2400,2600,2800,3000,3200,4000,5000)) If we use factors, and not continuous variates (cut versions of those two variates), > reg2=glm(cout~cut(ageconducteur,breaks=c(0,22,35,55,80,100))* + cut(agevehicule,breaks=c(-1,1,3,5,10,100)), + data=base,family=Gamma(link="log")) (note that we consider the Cartesian product, so values are computed for each product of factors, age of the driver and age of the car) we obtain Obviously, we’re missing something here: the most expensive class with one model is the cheapeast for the other one! Of course, it might come from our classes (that were chosen a bit randomly), but it might be interesting to use nonlinear functions of the ages. So, let us use splines to smooth those two variables, > reg3=glm(cout~bs(ageconducteur)+bs(agevehicule),data=base, + family=Gamma(link="log")) With additive smoothed functions, we obtained a symmetric graph (due to the additive property) while with a bivariate spline > library(mgcv) + reg4=gam(cout~s(ageconducteur,agevehicule),data=base, + family=Gamma(link="log")) (for some odd reasons, I could not use – easily – a bivariate spline in the Generalized Linear Model, but it did work considering a Generalized Additive Model – which is, by no means additive now). We can identify here some regions where the average cost can be extremely expensive… But, as mentioned wednesday, one should keep in mind that some parts of the square above are not reached. More precisely, the distribution of the portfolio, as a function of those two covariates is the following Thus, the proportion of young drivers driving a brand new car, and the proportion of old drivers driving a very old car is rather small… If the goal is to find niches, one should look at the prediction more carefully, but if the goal is to make that everyone gets an insurance cover, maybe we should allow that some drivers are under-priced (especially when are rare in the portfolio). And one should keep in mind that average costs are extremely sensitive to large losses, as discussed previously http://freakonometrics.hypotheses.org/3490 (and in class) In the univariate case, I have migrated an old post, we I tried to reproduce (in R and in French) some standard graphs in the insurance industry: it is always interesting to visualize not only the prediction obtained from our models, but also the size of each class in the portfolio, The post is online here http://freakonometrics.hypotheses.org/1224 # Données pour la régression logistique, et de Poisson Pour le cours de mercredi, deux petites bases, pour se pratiquer à modéliser des variable 0/1 ou une variable de comptage, > base = read.table("http://freakonometrics.free.fr/base-glm-act2040.txt", + header=TRUE) ou encore > base = read.table("http://freakonometrics.free.fr/base-pratique-act2040.txt", + header=TRUE) Sinon, une base plus complète pour faire de la tarification, > BASEN=read.table("http://freakonometrics.free.fr/baseN.txt",header=TRUE,sep=";") > BASEY=read.table("http://freakonometrics.free.fr/baseY.txt",header=TRUE,sep=";") > head(BASEN) ageconducteur agepermis sexeconducteur situationfamiliale habitation zone 1 57 39 F Celiba peri-urbain 8 2 54 35 H Celiba urbain 3 3 51 32 F Celiba urbain 1 4 53 35 H Marie rural 4 5 61 43 H Marie urbain 8 6 60 29 F Marie peri-urbain 1 agevehicule proprietaire payment marque poids usage 1 12 locataire Annuel AUTRES 8.>3500kg PROMENADE 2 20 sans mrp Semestriel PEUGEOT 4.3100-3199kg PROMENADE 3 4 sans mrp Annuel RAPIDO 1.<2700kg PROMENADE 4 1 sans mrp Annuel AUTRES 3.3000-3099kg PROMENADE 5 1 proprietaire Annuel FIAT 6.3300-3399kg PROMENADE 6 10 sans mrp Mensuel FIAT 8.>3500kg PROMENADE exposition nombre voiture 1 1 0 Monospace 2 1 0 Berline 3 1 0 sans avp 4 1 0 sans avp 5 1 1 Monospace 6 1 0 sans avp Parmi les variables, la description (sommaire) est la suivante, • ageconducteur: âge du conducteur principal du véhicule • agepermis: ancienneté du permis de conduire du conducteur principal du véhicule • sexeconducteur: sexe du conducteur principal (H ou F) • situationfamiliale: situation familiale du conducteur principal (“Celiba”, “Marie” ou “Veuf/Div”) • habitation: zone d’habitation du conducteur principal (“peri-urbain”, “rural” ou “urbain” ) • zone: zone d’habitation (allant de 1 à 8) • agevehicule: age du véhicule • proprietaire: si le conducteur principal possède un contrat Habitation, son statut (“locataire” ou “proprietaire”) Sinon “sans mrp” • payment:type de fractionnement de la prime d’assurance automobile (“Annuel”, “Mensuel” ou “Semestriel”) • marque: marque du véhicule > levels(BASEN[,10]) [1] "ADRIA" "AUTOSTAR" "AUTRES" "BURSTNER MOBIL" [5] "CHALLENGER" "CHAUSSON" "CITROEN" "FIAT" [9] "FORD" "HYMERMOBIL" "MERCEDES" "PEUGEOT" [13] "PILOTE" "RAPIDO" "RENAULT" "VOLKSWAGEN" • poids: classe de poids du véhicule > levels(BASEN[,11]) [1] "1.<2700kg" "2.2700-2999kg""3.3000-3099kg""4.3100-3199kg" [5] "5.3200-3299kg""6.3300-3399kg""7.3400-3499kg""8.>3500kg" • usage: utilisation du véhicule principal (“PROMENADE” ou “TOUS_DEPLACEMENTS”) • exposition: exposition, en années • nombre: nombre d’accident responsabilité civile du conducteur principal, pendant l’année passée • cout: cout du sinistre • voiture: type de véhicule > levels(BASEN[,15]) [1] "Berline" "Break" "Buggy" [4] "Cabriolet" "Combispace" "Coup\xe9" [7] "Coup\xe9 Cabriolet" "Jeep" "Minibus" [10] "Minispace" "Monospace" "sans avp" La variable d’intérêt est ici le nombre d’accident, > table(BASEN$nombre)

0     1
60155  3264

La base est un peu particulière – on en parlera en classe – les assurés ayant eu 0 ou 1 accident dans l’année.

# Further readings on GLMs and ratemaking

Some articles found in Actuarial journal, on ratemarking,

and in the CAS forums, and Astin conference papers

# Exposure with binomial responses

Last week, we’ve seen how to take into account the exposure to compute nonparametric estimators of several quantities (empirical means, and empirical variances) incorporating exposure. Let us see what can be done if we want to model a binomial response. The model here is the following: ,

• the number of claims $N_i$ on the period $[0,1]$ is unobserved
• the number of claims $Y_i$ on $[0,E_i]$ is observed (as well as $E_i$)

that can be visualize below

Consider the case where the variable of interest is not the number of claims, but simply the indicator of the occurrence of a claim. Then we wish to model the event $\{N=0\}$ versus $\{N>0\}$, interpreted as non-occurrence and occurrence. Given the fact that we can only observe $\{Y=0\}$ versus $\{Y>0\}$. Having an inclusion is not enough to derive a model. Actually, with a Poisson process model, we can get easily that

$\mathbb{P}(Y=0) = \mathbb{P}(N=0)^E$

With words, it means that the probability of not having a claim in the first six months of the year is the square root of not have a claim over a year. Which makes sense. Assume that the probability of not having a claim can be explained by some covariates, denoted $\boldsymbol{X}$, through some link function (using the GLM terminology),

$\mathbb{P}(N=0|\boldsymbol{X})=h(\boldsymbol{X}^{\text{\sffamily T}}\boldsymbol{\beta})$

Now, since we do observe $Y$ – and not $N$ – we have

$\mathbb{P}(Y=0|\boldsymbol{X},E)=h(\boldsymbol{X}^{\text{\sffamily T}}\boldsymbol{\beta})^E$

The dataset we will use is always the same

> sinistre=read.table("http://freakonometrics.free.fr/sinistreACT2040.txt",
> sinistres=sinistre[sinistre$garantie=="1RC",] > sinistres=sinistres[sinistres$cout>0,]
> T=table(sinistres$nocontrat) > T1=as.numeric(names(T)) > T2=as.numeric(T) > nombre1 = data.frame(nocontrat=T1,nbre=T2) > I = contrat$nocontrat%in%T1
> T1= contrat$nocontrat[I==FALSE] > nombre2 = data.frame(nocontrat=T1,nbre=0) > nombre=rbind(nombre1,nombre2) > sinistres = merge(contrat,nombre) > sinistres$nonsin = (sinistres$nbre==0) The first model we can consider is based on the standard logistic approach, i.e. $\mathbb{P}(Y=0|\boldsymbol{X},E)=\left(\frac{\exp(\boldsymbol{X}^{\text{\sffamily T}}\boldsymbol{\beta})}{1+\exp(\boldsymbol{X}^{\text{\sffamily T}}\boldsymbol{\beta})}\right)^E$ That’s nice, but difficult to handle with standard functions. Nevertheless, it is always possible to compute numerically the maximum likelihood estimator of $\boldymbol{\beta}$ given $(Y_i,\boldsymbol{X}_i,E_i)$. > Y=sinistres$nonsin
> X=cbind(1,sinistres$ageconducteur) > E=sinistres$exposition
> logL = function(beta){
+ 	pi=(exp(X%*%beta)/(1+exp(X%*%beta)))^E
+ 	-sum(log(dbinom(Y,size=1,prob=pi)))
+ }
> optim(fn=logL,par=c(-0.0001,-.001),
+ method="BFGS")
$par [1] 2.14420560 0.01040707$value
[1] 7604.073
$counts function gradient 42 10$convergence
[1] 0
$message NULL > parametres=optim(fn=logL,par=c(-0.0001,-.001), + method="BFGS")$par

Now, let us look at alternatives, based on standard regression models. For instance a binomial-log model. Because the exposure appears as a power of the annual probability, everything would be fine if $h$ was the exponential function (or $h^{-1}$ was the log link function), since

$\mathbb{P}(Y=0|\boldsymbol{X},E)=\exp(E+\boldsymbol{X}^{\text{\sffamily T}}\boldsymbol{\beta})$

Now, if we try to code it, it starts quickly to be problematic,

> reg=glm(nonsin~ageconducteur+offset(exposition),
Error: no valid set of coefficients has been found: please supply starting values

I tried (almost) everything I could, but I could not get rid of that error message,

> startglm=c(0,-.001)
> names(startglm)=c("(Intercept)","ageconducteur")
> etaglm=rep(-.01,nrow(sinistresI))
> etaglm[sinistresI$nonsin==0]=-10 > muglm=exp(etaglm) > reg=glm(nonsin~ageconducteur+offset(exposition), + data=sinistresI,family=binomial(link="log"), + control = glm.control(epsilon=1e-5,trace=TRUE,maxit=50), + start=startglm, + etastart=etaglm,mustart=muglm) Deviance = NaN Iterations - 1 Error: no valid set of coefficients has been found: please supply starting values So I decided to give up. Almost. Actually, the problem comes from the fact that $\mathbb{P}(Y=0)$ is closed to 1. I guess everything would be nicer if we could work with probability close to 0. Which is possible, since $\mathbb{P}(Y>0)=1-\mathbb{P}(Y=0) = 1-[1-\mathbb{P}(N>0)]^E$ where $\mathbb{P}(N>0)$ is close to 0. So we can use Taylor’s expansion, $\mathbb{P}(Y>0)\sim1-1+E\cdot \mathbb{P}(N>0)]=E\cdot \mathbb{P}(N>0)]$ Here, the exposure does no longer appears as a power of the probability, but appears multiplicatively. Of course, there are higher order terms. But let us forget them (so far). If – one more time – we consider a log link function, then we can incorporate the exposure, or to be more specific, the logarithm of the exposure. > regopp=glm((1-nonsin)~ageconducteur+offset(log(exposition)), + data=sinistresI,family=binomial(link="log")) which now works perfectly. Now, to see a final model, perhaps we should get back to our Poisson regression model since we do have a model for the probability that $\mathbb{P}(Y=\cdot)$. > regpois=glm(nbre~ageconducteur+offset(log(exposition)), + data=sinistres,family=poisson(link="log")) We can now compare those three models. Perhaps, we should also include the prediction without any explanatory variable. For the second model (actually, it does run without any explanatory variable), we run > regreff=glm((1-nonsin)~1+offset(log(exposition)), + data=sinistres,family=binomial(link="log")) so that the prediction is here > exp(coefficients(regreff)) (Intercept) 0.06776376 This value is comparable with the logistic regression, > logL2 = function(beta){ + pi=(exp(beta)/(1+exp(beta)))^E + -sum(log(dbinom(Y,size=1,prob=pi)))} > param=optim(fn=logL2,par=.01,method="BFGS")$par
> 1-exp(param)/(1+exp(param))
[1] 0.06747777

But is quite different from the Poisson model,

> exp(coefficients(glm(nbre~1+offset(log(exposition)),
(Intercept)
0.07279295

Let us produce a graph, to compare those models,

> age=18:100
> yml1=exp(parametres[1]+parametres[2]*age)/(1+exp(parametres[1]+parametres[2]*age))
> plot(age,1-yml1,type="l",col="purple")
> yp=predict(regpois,newdata=data.frame(ageconducteur=age,
+ exposition=1),type="response")
> yp1=1-exp(-yp)
> ydl=predict(regopp,newdata=data.frame(ageconducteur=age,
+ exposition=1),type="response")
> plot(age,ydl,type="l",col="red")
> lines(age,yp1,type="l",col="blue")
> lines(age,1-yml1,type="l",col="purple")
> abline(h=exp(coefficients(regreff)),lty=2)

Observe here that the three models are quite different. Actually, with two models, it is possible to run more complex regression, e.g. with splines, to visualize the impact of the age on the probability of having – or not – a car accident. If we compare the Poisson regression (still in red) and the log-binomial model, with Taylor’s expansion, we get

The next step is to see how to incorporate the exposure in a tree. But that’s another story…

# Natura non facit saltus

(see John Wilkins’ article on the – interesting – history of that phrase http://scienceblogs.com/evolvingthoughts/…). We will see, this week in class, several smoothing techniques, for insurance ratemaking. As a starting point, assume that we do not want to use segmentation techniques: everyone will pay exactly the same price.

• no segmentation of the premium

And that price should be related to the pure premium, which is proportional to the frequency (or the annualized frequency, as discussed previously), since

$\mathbb{E}_{\mathbb{P}}\left(\sum_{i=1}^N Y_i\right)=\mathbb{E}_{\mathbb{P}}(N) \cdot \mathbb{E}_{\mathbb{P}}(Y_i)$

The probability measure is mentioned here just to recall that we can use any measure. Even $\mathbb{P}_{\boldsymbol{X}}$ (based on some covariates). Without any covariate, the expected frequency should be

> regglm0=glm(nbre~1+offset(log(exposition)),data=sinistres,family=poisson)
> summary(regglm0)

Call:
glm(formula = nbre ~ 1 + offset(log(exposition)), family = poisson,
data = sinistres)

Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.5033  -0.3719  -0.2588  -0.1376  13.2700

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  -2.6201     0.0228  -114.9   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 12680  on 49999  degrees of freedom
Residual deviance: 12680  on 49999  degrees of freedom
AIC: 16353

Number of Fisher Scoring iterations: 6
> exp(coefficients(regglm0))
(Intercept)
0.07279295

Thus, if we do not want to take into account potential heterogeneity, we should assume that $N\sim\mathcal{P}(\lambda)$ where $\lambda$ is closed to 7.28%. Yes, as mentioned in class, it is rather common to see $\lambda$ as a percentage, i.e. a probability, since

$\mathbb{P}(N\neq 0)=1-e^{-\lambda}\approx \lambda$

i.e. $\lambda$ can be interpreted as the probability of not have a claim (see also the law of small numbers). Let us visualize this as a function of the age of the driver,

> a=18:100
> yp=predict(regglm0,newdata=data.frame(ageconducteur=a,exposition=1),type="response",se.fit=TRUE)
> yp0=yp$fit > yp1=yp$fit+2*yp$se.fit > yp2=yp$fit-2*yp$se.fit > plot(a,yp0,type="l",ylim=c(.03,.12)) > abline(v=40,col="grey") > lines(a,yp1,lty=2) > lines(a,yp2,lty=2) > k=23 > points(a[k],yp0[k],pch=3,lwd=3,col="red") > segments(a[k],yp1[k],a[k],yp2[k],col="red",lwd=3) We do predict the same frequency for all drivers, e.g. for some drive aged 40, > cat("Frequency =",yp0[k]," confidence interval",yp1[k],yp2[k]) Frequency = 0.07279295 confidence interval 0.07611196 0.06947393 Let us now consider the case where we try to take into account heterogeneity, e.g. by age, • The (standard) Poisson regression The idea of the (log-)Poisson regression is to assume that instead of having $N\sim\mathcal{P}(\lambda)$, we should have $N|\boldsymbol{X}\sim\mathcal{P}(\lambda_{\boldsymbol{X}})$, where $\lambda_{\boldsymbol{X}}=\exp(\beta_0+\beta_1 \boldsymbol{X}_1+\cdots+\beta_k\boldsymbol{X}_k)$ in a very general setting. Here, let us consider only one explanatory variable, i.e. $\lambda_{X}=\exp(\beta_0+\beta_1 {X})$ Here, we have > yp=predict(regglm1,newdata=data.frame(ageconducteur=a,exposition=1), + type="response",se.fit=TRUE) > yp0=yp$fit
> yp1=yp$fit+2*yp$se.fit
> yp2=yp$fit-2*yp$se.fit
> plot(a,yp0,type="l",ylim=c(.03,.12))
> abline(v=40,col="grey")
> lines(a,yp1,lty=2)
> lines(a,yp2,lty=2)
> points(a[k],yp0[k],pch=3,lwd=3,col="red")
> segments(a[k],yp1[k],a[k],yp2[k],col="red",lwd=3)

i.e. the prediction for the annualized claim frequency for our 40 year old driver is now 7.74% (which is slightly higher than what we had before, 7.28%)

> cat("Frequency =",yp0[k]," confidence interval",yp1[k],yp2[k])
Frequency = 0.07740574  confidence interval 0.08117512 0.07363636

It is possible to compute not the expected frequency , but the ratio $\mathbb{E}(N|X)/\mathbb{E}(N)$.

Above the horizontal blue line, the premium will be higher than the one obtained without segmentation, and (of course) lower below. Here, drivers younger than 44 year old will pay more, while driver older than 44 year old will be less. We have discussed, in the introduction, the necessity of segmentation. If we consider two companies, one segmenting, while the other one has a flat rate, then older drivers will go to the first company (since insurance is cheaper) while younger ones will go to the second one (again, it is cheaper). The problem is that the second company implicitly hopes that older drivers will compensate the risk. But since they’re gone, insurance will be too cheap, and the company will loose money (if not goes bankrupt). So companies have to use segmentation techniques to survive. Now, the problem is that we cannot be sure that this exponential decay of the premium is the proper way the premium should evolve as a function of the age. An alternative can be to use nonparametric techniques to visualize to true influence of the age on claims frequency.

• A pure nonparametric model

A first model can be to consider a premium, per age. This can be done considering the age of the driver as a factor in the regression,

> regglm2=glm(nbre~as.factor(ageconducteur)+offset(log(exposition)),
+ data=sinistres,family=poisson)
> yp=predict(regglm2,newdata=data.frame(ageconducteur=a0,exposition=1),
+ type="response",se.fit=TRUE)
> yp0=yp$fit > yp1=yp$fit+2*yp$se.fit > yp2=yp$fit-2*yp$se.fit > plot(a0,yp0,type="l",ylim=c(.03,.12)) > abline(v=40,col="grey") Here, the forecast for our 40 year old driver is slightly lower than be previous one, but the confidence interval is much larger (since we focus on a very small subclass of the portfolio: drivers aged exactly 40) Frequency = 0.06686658 confidence interval 0.08750205 0.0462311 Here, we consider too small classes, and the premium is too erratic: the premium will decrease of 20% from age 40 to 41, and then increase of 50% from age 41 to 42, > diff(log(yp0[23:25])) 24 25 -0.2330241 0.5223478 There is no chance that the company will keep the insured with this strategy. This discontinuity of the premium is clearly an important issue here. • Using age classes An alternative can be to consider age classes, from very young drivers to senior drivers. > level1=seq(15,105,by=5) > regglmc1=glm(nbre~cut(ageconducteur,level1)+offset(log(exposition)), + data=sinistres,family=poisson) > summary(regglmc1) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -1.6036 0.1741 -9.212 < 2e-16 *** cut(ageconducteur, level1)(20,25] -0.4200 0.1948 -2.157 0.0310 * cut(ageconducteur, level1)(25,30] -0.9378 0.1903 -4.927 8.33e-07 *** cut(ageconducteur, level1)(30,35] -1.0030 0.1869 -5.367 8.02e-08 *** cut(ageconducteur, level1)(35,40] -1.0779 0.1866 -5.776 7.65e-09 *** cut(ageconducteur, level1)(40,45] -1.0264 0.1858 -5.526 3.28e-08 *** cut(ageconducteur, level1)(45,50] -0.9978 0.1856 -5.377 7.58e-08 *** cut(ageconducteur, level1)(50,55] -1.0137 0.1855 -5.464 4.65e-08 *** cut(ageconducteur, level1)(55,60] -1.2036 0.1939 -6.207 5.40e-10 *** cut(ageconducteur, level1)(60,65] -1.1411 0.2008 -5.684 1.31e-08 *** cut(ageconducteur, level1)(65,70] -1.2114 0.2085 -5.811 6.22e-09 *** cut(ageconducteur, level1)(70,75] -1.3285 0.2210 -6.012 1.83e-09 *** cut(ageconducteur, level1)(75,80] -0.9814 0.2271 -4.321 1.55e-05 *** cut(ageconducteur, level1)(80,85] -1.4782 0.3371 -4.385 1.16e-05 *** cut(ageconducteur, level1)(85,90] -1.2120 0.5294 -2.289 0.0221 * cut(ageconducteur, level1)(90,95] -0.9728 1.0150 -0.958 0.3379 cut(ageconducteur, level1)(95,100] -11.4694 144.2817 -0.079 0.9366 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > yp=predict(regglmc1,newdata=data.frame(ageconducteur=a,exposition=1), + type="response",se.fit=TRUE) > yp0=yp$fit
> yp1=yp$fit+2*yp$se.fit
> yp2=yp$fit-2*yp$se.fit
> plot(a,yp0,ylim=c(.03,.12),type="s")
> abline(v=40,col="grey")
> lines(a,yp1,lty=2,type="s")
> lines(a,yp2,lty=2,type="s")

Here we obtain the following predictions,

and for our 40 year old driver, the frequency is now 6.84%.

Frequency = 0.0684573  confidence interval 0.07766717 0.05924742

But our classes were defined arbitrarily here. Perhaps should we consider other classes, to see if the prediction is sensitive to the cutting values,

> level2=level1-2
> regglmc2=glm(nbre~cut(ageconducteur,level2)+offset(log(exposition)),
+ data=sinistres,family=poisson)

which yields the following values for our 40 year old driver,

Frequency = 0.07050614  confidence interval 0.07980422 0.06120807

So here, we did not remove the discontinuity problem. An idea here can be to consider moving regions: if the goal is to predict the frequency for a 40 year old driver, perhaps the class should be (somehow) centered around 40. And center the interval around 35 for drivers aged 35. Etc.

• Moving average

Thus, it is natural to consider some local regressions, where only drivers aged almost 40 should be considered. This almost concept is related to the bandwidth. For instance, drivers between 35 and 45 can be considered as being almost40. In practice we can either consider a subset function, or we can use weights in the regressions

> value=40
> h=5
> sinistres$omega=(abs(sinistres$ageconducteur-value)<=h)*1
> regglmomega=glm(nbre~ageconducteur+offset(log(exposition)),
+ data=sinistres,family=poisson,weights=omega)

To see what’s going on, let us consider an animated plot, where the age of interest is changing,

Here, for our 40 year old drive, we get

Frequency = 0.06913391  confidence interval 0.07535564 0.06291218

We do obtain a curve that can be interpreted as a local regression. But here, we do not take into account that 35 is not as close to 40 as 39 could be. An here, 34 is assumed to be very far away from 40. Clearly, we could improve that technique: kernel functions can considered, i.e. the closer to 40, the larger the weight.

> value=40
> h=5
> sinistres$omega=dnorm(abs(sinistres$ageconducteur-value)/h)
> regglmomega=glm(nbre~ageconducteur+offset(log(exposition)),
+ data=sinistres,family=poisson,weights=omega)

which can be plotted below

Here, our prediction for our 40 year old drive is

Frequency = 0.07040464  confidence interval 0.07981521 0.06099408

This is the idea of kernel regression techniques. But as explained in the slides, other non parametric techniques can be considered, like spline functions.

• Smoothing with splines

In R, it is simple to use spline function (somehow much more simple than kernel smoothers)

> library(splines)
> regglmbs=glm(nbre~bs(ageconducteur)+offset(log(exposition)),
+ data=sinistres,family=poisson)

The prediction for our 40 year old driver is now

Frequency = 0.06928169  confidence interval 0.07397124 0.06459215

Note that this techniques is related to another class of models, the so-called Generalized Additive Models, i.e. GAMs.

> library(mgcv)
> reggam=gam(nbre~s(ageconducteur)+offset(log(exposition)),
+ data=sinistres,family=poisson)

The prediction is extremely close to the one we obtained above (the main differences being observed for very old drivers)

Frequency = 0.06912683  confidence interval 0.07501663 0.06323702
• Comparison of the different models

Somehow, one way or another, all those models are valid. So perhaps we should compare them,

On the graph above, we can visualize the upper and the lower bound of the prediction, for the 9 models. The horizontal line is the predicted value without taking into account heterogeneity. It is possible to consider relative values, with respect to this value,

# Fréquence de sinistres, et surdispersion

Je continue à mettre en ligne les transparents qui serviront de support pour le cours ACT2040. Dans cette dernière partie sur la modélisation de la fréquence de sinistre, on parlera de surdispersion. Les transparents sont en ligne ici,

Sinon, parmi les références complémentaires, je peux suggérer plusieurs documents rédigés par des praticiens, comme Meyers (2009) http://casact.org/education/…, Isamail & Jemain (2009) http://casact.org/pubs/… ou encore le document très intéressant (et critique) de Schmid (2011) http://casact.org/education/…. Les plus motivés pourront aussi survoler les section 2.3 et 2.4. du livre Denuit et al. (2007), en ligne sur http://books.google.ca/…

# Introduction aux modèles linéaires généralisés

J’ai un peu d’avance dans le cours. Je vais mettre en ligne les transparents pour la semaine prochaine (normalement), où nous aborderons la classe des modèles linéaires généralisés. Les transparents sont en ligne ici.

Je n’ai pas mis de section sur lesGeneralized Additive Models, on se contentera de la section sur le lissage évoquée à la fin des transparents sur la modélisation de la fréquence. Afin de légitimer les méthodes de lissage (sur l’âge de l’assuré en particulier), je renvoie vers un graphique produit il y a plusieurs années par un cabinet de conseil, qui notait que la forme de la fonction de lissage, liant l’âge à la fréquence est identique, dans tous les pays,

Mais je pense que je ferais un billet dédié au lissage, dans la problématique de la tarification en assurance IARD.

# Régression de Poisson, et biais minimal

Lors du prochain cours d’actuariat, nous allons finir les arbres de régression, et introduire la régression de Poisson. Les transparents sont en ligne ici,

Je vais présenter la régression en Poisson, en faisant un parallèle avec la régression logistique, la session suivante portera sur la généralisation obtenue avec les modèles linéaires généralisés. Sur la régression de Poisson, je suggère de lire Frees (2010) chapitre 12 (p 343-361), Greene (2012), section 18.3 (p 802-828) ou encore de Jong Heller (2008) chapitre 6. Sur les méthodes de biais minimal, de Jong Heller (2008), section 1.3 et l’article de Sholom Feldblum, http://www.casact.org/…. Sur le passage de ces dernières méthodes (introduites par Robert Bailey dans les années 60, http://www.casact.org/… et http://www.casact.org/…), je recommande la lecture de l’article de Ben Zehnwirth, Ratemaking From Bailey and Simon (1960) to Generalized Linear Regression Models, en ligne sur http://www.casact.org/…

Comme annoncé au premier cours, j’essaye de mettre en ligne les transparents au fur et à mesure, mais j’avais pris l’habitude d’écrire au tableau ces dernières années. Il faut donc que je tape tout. Pour le devoir un courriel sera envoyé d’ici la fin de semaine à tous les groupes qui se sont inscrits.

# Régression logistique et arbres

Pour le cours de mercredi prochain, la base utilisée sera une base tirée du livre de Jed Frees, http://instruction.bus.wisc.edu/jfrees/…

> baseavocat=read.table("http://freakonometrics.free.fr/AutoBI.csv",header=TRUE,sep=",")
> tail(baseavocat)
CASENUM ATTORNEY CLMSEX MARITAL CLMINSUR SEATBELT CLMAGE  LOSS
1335   34204        2      2       2        2        1     26 0.161
1336   34210        2      1       2        2        1     NA 0.576
1337   34220        1      2       1        2        1     46 3.705
1338   34223        2      2       1        2        1     39 0.099
1339   34245        1      2       2        1        1     18 3.277
1340   34253        2      2       2        2        1     30 0.688

On dispose d’une variable dichotomique indiquant si un assuré – suite à un accident de la route – a été représenté par un avocat (1 si oui, 2 si non). On connaît le sexe de l’assuré (1 pour les hommes et 2 pour les femmes), le statut marital (1 s’il est marié, 2 s’il est célibataire, 3 pour un veuf, et 4 pour un assuré divorcé). On sait aussi si l’assuré portait ou non une ceinture de sécurité lorsque l’accident s’est produit (1 si oui, 2 si non et 3 si l’information n’est pas connue). Enfin, une information pour savoir si le conducteur du véhicule était ou non assuré (1 si oui, 2 si non et 3 si l’information n’est pas connue). On va recoder un peu les données afin de les rendre plus claires à lire.

Les transparents du cours sont en ligne ici,

Sur les arbres de régression, je mettrais en ligne un billet, afin d’illustrer la méthode. En attendant des compléments théoriques peuvent se trouver en ligne http://genome.jouy.inra.fr/…, http://ensmp.fr/…, ou http://ujf-grenoble.fr/… (pour information, nous ne verrons que la méthode CART). Je peux renvoyer au livre (et au blog) de Stéphane Tuffery, ou (en anglais) au livre de Richard Berk, dont un résumé se trouve en ligne sur http://crim.upenn.edu/….

La semaine suivante, nous aborderons la régression de Poisson, les méthodes de biais minimal, et introduirons les modèles linéaires généralisés. Je renvoie au chapitre sur la tarification a priori du Denuit & Charpentier (2005), aux chapitres 12 et 13 de Frees (2010) ou encore les chapitres 5 et 6 du De Jong  & Heller (2008). Pour les plus curieux qui veulent comprendre les liens entre les modèles linéaires généralisés et la tarification par crédibilité, je renvoie à l’article de Klinker (2010)

# Segmentation en tarification, compléments

Dans le premier cours d’actuariat IARD, nous avons vu l’importance de la ségmentation, et son implication sur le calcul des primes (passer d’une espérance mathématique à une espérance conditionnelle). Pour aller un peu plus loin, quelques compléments,

pour une lecture plus économique de la problématique de la segmentation en assurance

ou pour une lecture plus légale

Sinon, plusieurs articles de vulgarisation peuvent être lu sur internet,

La première démo aura lieu lundi, en salle informatique. Karim sera une introduction au langage R, à la manipulation des variables (qualitatives et quantitatives). Je mettrais en ligne les transparents en fin de semaine, et les codes seront mis en ligne dans le courant de la semaine prochaine.

Comme annoncé hier, il n’y aura pas cours mercredi prochain. Le mercredi suivant, nous verrons la modélisation des variables indicatrices, i.e. la régression logistique, et les arbres de régression. On supposera que le modèle linéaire aura été vu, je mets un lien vers les transparents du cours ACT6420 de la session passée, notes de cours transparents1 et transparents2. Il est aussi possible de relire Frees (2010), chapitres 3, 4, 5 et 6.

Pour commencer à pratiquer la régression logistique, on utilisera la petite fonction suivante

logit = function(formula, lien="logit", data=NULL) {
}

Sinon, la Casualty Actuarial Society a mis en ligne plusieurs documents en ligne sur les arbres de régression (qui sont peu abordés dans les livres mentionnés auparavant),

pour une comparaison de toutes les méthodes

Les transparents seront mis en ligne en fin de semaine prochaine. A suivre donc…

# Incréments négatifs dans les triangles de paiements

Considérons le triangle d’incréments de payements suivants,

> PAID
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 3209 4372 4411 4428 4435 4456
[2,] 3367 4659 4696 4720 4730   NA
[3,] 3871 5345 5338 5420   NA   NA
[4,] 4239 5917 6020   NA   NA   NA
[5,] 4929 6794   NA   NA   NA   NA
[6,] 5217   NA   NA   NA   NA   NA
> INC=PAID
> INC[,2:6]=PAID[,2:6]-PAID[,1:5]
> INC
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 3209 1163   39   17    7   21
[2,] 3367 1292   37   24   10   NA
[3,] 3871 1474   -7   82   NA   NA
[4,] 4239 1678  103   NA   NA   NA
[5,] 4929 1865   NA   NA   NA   NA
[6,] 5217   NA   NA   NA   NA   NA

Comme on peut le voir, sur la troisième ligne, on a un incrément de paiementnégatif. A priori, ce n’est pas forcément gênant. En particulier, on peut faire tourner la méthode Chain Ladder,

> lambda=rep(NA,5)
> for(k in 1:5){
+ lambda[k]=(sum(PAID[1:(6-k),k+1])/
+             sum(PAID[1:(6-k),k]))}
> lambda
[1] 1.380933 1.008476 1.008515 1.001858 1.004735
> PROJECTION=PAID
> for(k in 1:5){
+ PROJECTION[((7-k):6),k+1]=
+ PROJECTION[((7-k):6),k]*lambda[k]}
> sum(PROJECTION[,6]-
+ diag(PROJECTION[,6:1]))
[1] 2469.703

et la sortie coïncide avec la fonction R,

> MackChainLadder(PAID)

Latest Dev.To.Date Ultimate    IBNR Mack.S.E CV(IBNR)
1  4,456       1.000    4,456     0.0    0.000      NaN
2  4,730       0.995    4,752    22.4    0.146  0.00652
3  5,420       0.993    5,456    35.8    2.405  0.06721
4  6,020       0.985    6,111    91.3   41.679  0.45629
5  6,794       0.977    6,956   161.5   71.620  0.44334
6  5,217       0.707    7,376 2,158.6   95.750  0.04436

Totals
Latest:            32,637.00
Dev:                    0.93
Ultimate:          35,106.70
IBNR:               2,469.70
Mack S.E.:            146.62
CV(IBNR):  0.059366227164502
Message d'avis :
In Mack.S.E(CL[["Models"]], FullTriangle, est.sigma = est.sigma,
'loglinear' model to estimate sigma_n doesn't appear appropriate
p-value > 5.
est.sigma will be overwritten to 'Mack'.
Mack's estimation method will be used instead.

On notera le message d’avis, laissant entendre qu’il peut y avoir un petit soucis. En fait, le soucis apparaît clairement si on souhaite faire une régression de Poisson. Car autant les valeurs non-entières ne sont pas trop gênantes, autant les valeurs négatives le sont !

> Y=as.vector(INC)
> D=rep(1:6,each=6)
> A=rep(2001:2006,6)
> base=data.frame(Y,D,A)
> reg=glm(Y~as.factor(D)+as.factor(A),
Erreur dans eval(expr, envir, enclos) :
les valeurs négatives sont interdites pour la famille poisson

Il faut alors trouver une solution si on a des incréments négatifs. Et aucune solution ne sera intellectuellement satisfaisante, car avec une valeur négative, il n’y a aucune légitimité pour continuer à utiliser un modèle Poissonnien. Donc on va bricoler…

• Prendre de l’argent à droite et à gauche

La première solution consiste à dire que cet incrément n’a pas de raison d’être. On va donc aller chercher de l’argent dans la colonne avant (sur la même ligne, on va supposer qu’il s’agit d’un problème sur les cadences de paiements) ou sur la colonne après. Au lieu d’avoir

[3,] 3871 1474   -7   82   NA   NA

On peut renflouer en prenant à gauche,

[3,] 3871 1467   0    82   NA   NA

ou à droite

[3,] 3871 1474   0    75   NA   NA

En faisant ces opérations, on change seulement localement la courbe de paiements pour la troisième année de survenance. On peut se demander si le fait de prendre à gauche a un impact sur l’estimation du montant de provisions, ou sur l’incertitude associée. Pour ça, on peut utiliser les fonctions suivantes,

library(ChainLadder)
CL=function(T){
sum(
MackChainLadder(T)$FullTriangle[,ncol(T)] -diag(T[,ncol(T):1]))} CL.SE=function(T){ MackChainLadder(T)$Total.Mack.S.E}
CumInc=function(T){
m=T
for(i in 2:nrow(T)){m[,i]=apply(T[,1:i],1,sum)}
return(m)}

On peut alors regarder, si on prend à droite ou à gauche ce qui se passe, ou ce qui se passe si on renfloue non pas à 0 mais à une valeur strictement positive (e.g. 1),

VCL=rep(NA,8)
for(j in 1:8){
INCd=INC
INCd[3,2:4]=INC[3,2:4]+c(1-j,+7,j-1-7)
VCL[j]=CL(CumInc(INCd))}
plot(1:8,VCL,type="b",col="blue",xlim=c(0,9),ylim=c(2462,2474))
VCL=rep(NA,10)
for(j in 1:10){
INCd=INC
INCd[3,2:4]=INC[3,2:4]+c(1-j,+9,j-1-9)
VCL[j]=CL(CumInc(INCd))}
lines(0:9,VCL,type="b",pch=0,col="red")

avec la courbe bleu si on renfloue à 0, et rouge si on renfloue à 1, pour l’estimation du montant total de provisions,

VCL=rep(NA,8)
for(j in 1:8){
INCd=INC
INCd[3,2:4]=INC[3,2:4]+c(1-j,+7,j-1-7)
VCL[j]=CL.SE(CumInc(INCd))}
plot(1:8,VCL,type="b",col="blue",xlim=c(0,9),ylim=c(130,145))
VCL=rep(NA,10)
for(j in 1:10){
INCd=INC
INCd[3,2:4]=INC[3,2:4]+c(1-j,+9,j-1-9)
VCL[j]=CL.SE(CumInc(INCd))}
lines(0:9,VCL,type="b",pch=0,col="red")

avec cette fois l’impact des transferts sur l’écart-type. En abscisse, on a (à peu de choses près) le montant enlevé sur la colonne de droite. Bref, les transferts venant de la gauche ou de la droite sont une solution, mais le choix d’où vient l’argent ne sera pas neutre, ni sur l’estimation, ni sur la variance de l’estimation (et l’intervalle de confiance).

• Jouer à faire des translations…

Une autre piste pourrait être de noter que, dans le modèle linéaire, si on translate nos données (vers le haut), on ne change pas la pente, et la constante est juste augmenté (de la taille de la translation).

> lm(dist~speed,data=cars)

Call:
lm(formula = dist ~ speed, data = cars)

Coefficients:
(Intercept)        speed
-17.579        3.932

> lm((dist+10)~speed,data=cars)

Call:
lm(formula = (dist + 10) ~ speed, data = cars)

Coefficients:
(Intercept)        speed
-7.579        3.932

Autrement dit, la prédiction faite par notre modèle n’est pas modifiée par la translation, à condition d’opérer la translation ensuite sur la prédiction. Sur le dessin ci-dessous, on fait pareil, mais avec une régression de Poisson: la courbenoire est sur les données brutes, la courbe bleu est obtenue sur les points translatés vers le haut, et la rouge si on translate la prédiction sur les points translatés,

Si la translation est trop importante, la prédiction finale (en rouge) se retrouve davantage éloignée de la vraie valeur (en noir).
Pareillement, on peut translater nos paiements vers le haut dans les triangles, de manière à avoir des paiements positifs. On peut commencer par translater tout le triangle, par exemple en translatant juste assez pour que tous les incréments soient positifs ou nuls,

> k=7
> Y=as.vector(INC)+k
> D=rep(1:6,each=6)
> A=rep(2001:2006,6)
> base=data.frame(Y,D,A)
> reg=glm(Y~as.factor(D)+as.factor(A),
> Yp=predict(reg,type="response",
+ newdata=base)-k
> sum(Yp[is.na(Y)==TRUE])
[1] 2508.620

On peut aussi se demander ce qui se passerait si on décalait davantage vers le haut. On l’a vu sur l’animation, plus la translation est importante, plus la prédiction se décale. Un stratégie pourrait être de faire une estimation pour plusieurs valeurs de translations – de telle sorte que tous les incréments soient positifs ou nuls – et ensuite d’extrapoler pour savoir ce qui se passerait si on ne translatait pas,

> K=7:20
> R=rep(NA,length(K))
> for(i in 1:length(K)){
+ k=K[i]
+ Y=as.vector(INC)+k
+ D=rep(1:6,each=6)
+ A=rep(2001:2006,6)
+ base=data.frame(Y,D,A)
+ reg=glm(Y~as.factor(D)+as.factor(A),
+ Yp=predict(reg,type="response",
+ newdata=base)-k
+ R[i]=sum(Yp[is.na(Y)==TRUE])}
> plot(K,R,xlim=c(0,20),ylim=c(2465,max(R)))
> abline(lm(R~K),col="blue")
> (yp=predict(lm(R~K),newdata=(K=0)))
1
2470.199
> points(0,yp,col="red",pch=19)

On a cette fois une estimation plus proche de celle que nous avions en utilisant directement, sur le triangle brut, la méthode chain ladder. Mais on peut aller un peu plus loin. Car dans le triangle, on a translaté toutes les valeurs, alors que seule une était négative. On pourrait par exemple translater seulement les valeurs de la troisième colonne,

> K=7:20
> R=rep(NA,length(K))
> for(i in 1:length(K)){
+ k=K[i]
+ Y=as.vector(INC)
+ D=rep(1:6,each=6)
+ A=rep(2001:2006,6)
+ Y[D==3]=Y[D==3]+k
+ base=data.frame(Y,D,A)
+ reg=glm(Y~as.factor(D)+as.factor(A),
+ Yp=predict(reg,type="response",
+ newdata=base)
+ Yp[D==3]=Yp[D==3]-k
+ R[i]=sum(Yp[is.na(Y)==TRUE])}
> predict(lm(R~K),newdata=(K=0))
1
2469.703

Tiens, on retombe exactement sur l’estimateur de la méthode chain ladder… Étonnant, non ?

Une petite précision. Dans la vraie vie, on ne devrait pas avoir d’incréments négatifs dans les triangles de paiements. Maintenant, il faut reconnaître qu’on en voit apparaître lorsque l’on bootstrappe les résidus, et que l’on génère de pseudo-triangles…
Si des praticiens ont des commentaires sur les incréments négatifs, les commentaires sont ouverts (et peu modérés), donc racontez nous comment vous faites, je suis preneur…

# ACT2040: examen final

L’examen final du cours ACT2040 aura lieu le 13 décembre, durant trois heures. L’examen sera en deux parties: une de tarification, et une de provisionnement. Les deux porteront sur la lecture et l’analyse critique de sorties informatiques. Pour la tarification, les bases de données utilisées sont en ligne ci-dessous

> BASEN=read.table("http://freakonometrics.free.fr/baseN.txt",header=TRUE,sep=";")
ageconducteur agepermis sexeconducteur situationfamiliale  habitation zone
1            57        39              F             Celiba peri-urbain    8
2            54        35              H             Celiba      urbain    3
3            51        32              F             Celiba      urbain    1
4            53        35              H              Marie       rural    4
5            61        43              H              Marie      urbain    8
6            60        29              F              Marie peri-urbain    1
agevehicule proprietaire    payment  marque         poids     usage
1          12    locataire     Annuel  AUTRES     8.>3500kg PROMENADE
2          20     sans mrp Semestriel PEUGEOT 4.3100-3199kg PROMENADE
3           4     sans mrp     Annuel  RAPIDO     1.<2700kg PROMENADE
4           1     sans mrp     Annuel  AUTRES 3.3000-3099kg PROMENADE
5           1 proprietaire     Annuel    FIAT 6.3300-3399kg PROMENADE
6          10     sans mrp    Mensuel    FIAT     8.>3500kg PROMENADE
exposition nombre   voiture
1          1      0 Monospace
2          1      0   Berline
3          1      0  sans avp
4          1      0  sans avp
5          1      1 Monospace
6          1      0  sans avp

Parmi les variables, la description (sommaire) est la suivante,

• ageconducteur: âge du conducteur principal du véhicule
• agepermis: ancienneté du permis de conduire du conducteur principal du véhicule
• sexeconducteur: sexe du conducteur principal (H ou F)
• situationfamiliale: situation familiale du conducteur principal (“Celiba”, “Marie” ou “Veuf/Div”)
• habitation: zone d’habitation du conducteur principal (“peri-urbain”, “rural” ou “urbain” )
• zone: zone d’habitation (allant de 1 à 8)
• agevehicule: age du véhicule
• proprietaire: si le conducteur principal possède un contrat Habitation, son statut (“locataire” ou “proprietaire”)  Sinon “sans mrp”
• payment:type de fractionnement de la prime d’assurance automobile (“Annuel”, “Mensuel” ou “Semestriel”)
• marque: marque du véhicule
> levels(BASEN[,10])
[1] "ADRIA"       "AUTOSTAR"    "AUTRES"      "BURSTNER MOBIL"
[5] "CHALLENGER"  "CHAUSSON"    "CITROEN"     "FIAT"
[9] "FORD"        "HYMERMOBIL"  "MERCEDES"    "PEUGEOT"
[13] "PILOTE"     "RAPIDO"      "RENAULT"     "VOLKSWAGEN"
• poids: classe de poids du véhicule
> levels(BASEN[,11])
[1] "1.<2700kg"    "2.2700-2999kg""3.3000-3099kg""4.3100-3199kg"
[5] "5.3200-3299kg""6.3300-3399kg""7.3400-3499kg""8.>3500kg"
• usage: utilisation du véhicule principal (“PROMENADE” ou “TOUS_DEPLACEMENTS”)
• exposition: exposition, en années
• nombre: nombre d’accident responsabilité civile du conducteur principal, pendant l’année passée
• cout: cout du sinistre
• voiture: type de véhicule
> levels(BASEN[,15])
[1] "Berline"            "Break"              "Buggy"
[4] "Cabriolet"          "Combispace"         "Coup\xe9"
[7] "Coup\xe9 Cabriolet" "Jeep"               "Minibus"
[10] "Minispace"          "Monospace"          "sans avp"

Concernant la partie sur le provisionnement, des sorties associées au triangle suivant seront proposées,

> TRIANGLE=read.table("http://freakonometrics.blog.free.fr/public/code/triangleACT2040.txt")
> TRIANGLE
X1      X2      X3      X4      X5      X6      X7      X8      X9     X10
1  357848 1124788 1735330 2218270 2745596 3319994 3466336 3606286 3833515 3901463
2  352118 1236139 2170033 3353322 3799067 4120063 4647867 4914039 5339085      NA
3  290507 1292306 2218525 3235179 3985995 4132918 4628910 4909315      NA      NA
4  310608 1418858 2195047 3757447 4029929 4381982 4588268      NA      NA      NA
5  443160 1136350 2128333 2897821 3402672 3873311      NA      NA      NA      NA
6  396132 1333217 2180715 2985752 3691712      NA      NA      NA      NA      NA
7  440832 1288463 2419861 3483130      NA      NA      NA      NA      NA      NA
8  359480 1421128 2864498      NA      NA      NA      NA      NA      NA      NA
9  376686 1363294      NA      NA      NA      NA      NA      NA      NA      NA
10 344014      NA      NA      NA      NA      NA      NA      NA      NA      NA