# 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))
+ 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))
+ 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…

# Bootstrap et régression

Lors du dernier cours, on a évoqué l’utilisation du bootstrap pour obtenir des intervalles de confiance sur des prévisions. Je mets en ligne les codes tapés en cours (très sommairement commentés, je peux renvoyer vers des vieux billets du cours ACT6420 pour des compléments). On va travailler sur ma base préférée pour évoquer la régression linéaire (avant de parler triangles de provisionnement, revenons cinq minutes sur des choses simples).

> plot(cars)
> reg=lm(dist~speed,data=cars)
> abline(reg,col="red")
> n=nrow(cars)
> x=21
> points(x,predict(reg,newdata= data.frame(speed=x)),pch=19,col="red")

On cherche ici à faire une prédiction en un point. Comme rappelé en cours (mais aussi dans le cours de modèles de prévision), quand on veut donner un intervalle de confiance pour la prévision, il convient de distinger l’intervalle de confiance pour le prédicteur (qui va dépendre de l’erreur d’estimation des paramétres) et l’intervalle de confiance pour une potentielle valeur (on peut parler de génération de scénarios, qui va dépendre en plus de l’erreur de modèle, c’est à dire de la dispersion des résidus). Commençons par l’intervalle de confiance sur la prédiction, sur le best estimate comme on dit en provisionnement

> Yx=rep(NA,500)
> B=matrix(NA,500,2)
> for(s in 1:500){
+ indice=sample(1:n,size=n,
+ replace=TRUE)
+ base=cars[indice,]
+ #points(base,pch=3)
+ reg=lm(dist~speed,data=base)
+ abline(reg,col="light blue")
+ points(x,predict(reg,newdata=data.frame(speed=x)),pch=19,col="blue")
+ Yx[s]=predict(reg,newdata=data.frame(speed=x))
+ B[s,]=coefficients(reg)
+ }

Les valeurs bleues sont ici des prévisions possibles, obtenues en rééchantillonnant dans notre base d’observations. Pour rappel, l’intervalle de confiance (à 90%), sous hypothèse de normalité des résidus (et donc des estimateurs de la pente et de la constante de la droite de régression) s’obtient de la manière suivante

> reg=lm(dist~speed,data=cars)
> U=predict(reg,interval ="confidence",
+ level=.9,newdata=
+ data.frame(speed=0:30))
> lines(0:30,U[,2],col="red",lwd=2)
> lines(0:30,U[,3],col="red",lwd=2)

On peut comparer ici la distribution des valeurs obtenues sur nos 500 jeux de données générées, et comparer les quantiles empiriques, avec les quantiles sous hypothèse de normalité,

> hist(Yx,proba=TRUE,col="light blue",border="white")
> abline(v=U[x+1,2:3],col="red",lwd=2)
> D=density(Yx)
> lines(D)
> I=which(D$x<=quantile(Yx,.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(Yx,.95)) > polygon(c(D$x[I],rev(D$x[I])),c(D$y[I],rep(0,length(I))),col="blue",border=NA)

On peut noter que les ordres de grandeur sont comparables.

> reg=lm(dist~speed,data=cars)
> quantile(Yx,c(.05,.95))
5%      95%
58.63689 70.31281
> predict(reg,interval ="confidence",
+ level=.9,newdata=data.frame(speed=x))
fit      lwr      upr
1 65.00149 59.65934 70.34364

Regardons maintenant l’autre type d’intervalle de confiance, sur la valeur possible de la variable d’intérêt. Cette fois, en plus de tirer des nouveaux échantillons et calculer des prédictions, on va en plus rajouter un bruit à chaque tirage, qui permettra d’obtenir une valeur possible.

> Yx=rep(NA,500)
> for(s in 1:500){
+ indice=sample(1:n,size=n,
+ replace=TRUE)
+ base=cars[indice,]
+ #points(base,pch=3)
+ reg=lm(dist~speed,data=base)
+ erreur=residuals(reg)
+ #abline(reg,lty=2)
+ E=sample(erreur,size=1)
+ Yx[s]=predict(reg,newdata=data.frame(speed=x))+E
+ points(x,Yx[s],pch=19,col="red")
+ }

Là encore, on peut comparer (graphiquement pour commencer) les valeurs obtenues par rééchantillonnage, et celle obtenues sous hypothèse de normalité,

> hist(Yx,proba=TRUE,col="light blue",border="white")
> abline(v=U[2:3],col="red",lwd=2)
> D=density(Yx)
> lines(D)
> I=which(D$x<=quantile(Yx,.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(Yx,.95)) > polygon(c(D$x[I],rev(D$x[I])),c(D$y[I],rep(0,length(I))),col="blue",border=NA)

Ce qui donne, numériquement, les comparaisons suivantes

> quantile(Yx,c(.05,.95))
5%      95%
44.43468 96.01357
> (U=predict(reg,interval ="prediction",level=.9,newdata=data.frame(speed=x)))
fit      lwr      upr
1 67.63136 45.16967 90.09305

On observe cette fois une légère asymméytrie vers la droite. Manifestement, on ne peut pas supposer les résidus Gaussien, car il y a plus de grandes valeurs positives, que négatives. Ce qui fait du sens compte tenu de la nature des données (une distance ne peut être négative).

On avait ensuite commencé à discuter de l’utilisation des modèles de régression en provisionnement. Afin d’avoir des données présentant de l’indépendance, on avait rappelé qu’il fallait travailler avec les incréments de paiments, et non pas les paiements cumulés.

> T
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 3209 4372 4411 4428 4435 4456
[2,] 3367 4659 4696 4720 4730   NA
[3,] 3871 5345 5398 5420   NA   NA
[4,] 4239 5917 6020   NA   NA   NA
[5,] 4929 6794   NA   NA   NA   NA
[6,] 5217   NA   NA   NA   NA   NA
> n=ncol(T)
> Y=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

On peut alors constituer une base de données, avec comme variables explicatives la ligne et la colonne.

> y=as.vector(as.matrix(Y))
> base=data.frame(
+ y,
+ ai=rep(2000:2005,n),
+ bj=rep(0:(n-1),each=n))
>
y   ai bj
1  3209 2000  0
2  3367 2001  0
3  3871 2002  0
4  4239 2003  0
5  4929 2004  0
6  5217 2005  0
7  1163 2000  1
8  1292 2001  1
9  1474 2002  1
10 1678 2003  1
11 1865 2004  1
12   NA 2005  1
> tail(base,12)
y   ai bj
25  7 2000  4
26 10 2001  4
27 NA 2002  4
28 NA 2003  4
29 NA 2004  4
30 NA 2005  4
31 21 2000  5
32 NA 2001  5
33 NA 2002  5
34 NA 2003  5
35 NA 2004  5
36 NA 2005  5

On peut alors commencer par utiliser le modèle Regression models based on log-incremental payments de Stavros Christofides, basé sur une modélisation lognormale, introduite initialement par Etienne de Vylder en 1978 (Markus en parle, en trois parties, sur son blog http://lamages.blogspot.ca/Barnett%20Zehnwirth)

> reg1=lm(log(y)~
+ as.factor(ai)+
+ as.factor(bj),data=base)
> summary(reg1)

Call:
lm(formula = log(y) ~ as.factor(ai) + as.factor(bj), data = base)

Residuals:
Min       1Q   Median       3Q      Max
-0.26374 -0.05681  0.00000  0.04419  0.33014

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)         7.9471     0.1101  72.188 6.35e-15 ***
as.factor(ai)2001   0.1604     0.1109   1.447  0.17849
as.factor(ai)2002   0.2718     0.1208   2.250  0.04819 *
as.factor(ai)2003   0.5904     0.1342   4.399  0.00134 **
as.factor(ai)2004   0.5535     0.1562   3.543  0.00533 **
as.factor(ai)2005   0.6126     0.2070   2.959  0.01431 *
as.factor(bj)1     -0.9674     0.1109  -8.726 5.46e-06 ***
as.factor(bj)2     -4.2329     0.1208 -35.038 8.50e-12 ***
as.factor(bj)3     -5.0571     0.1342 -37.684 4.13e-12 ***
as.factor(bj)4     -5.9031     0.1562 -37.783 4.02e-12 ***
as.factor(bj)5     -4.9026     0.2070 -23.685 4.08e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1753 on 10 degrees of freedom
(15 observations deleted due to missingness)
Multiple R-squared: 0.9975,	Adjusted R-squared: 0.9949
F-statistic: 391.7 on 10 and 10 DF,  p-value: 1.338e-11

> base$py=exp(predict(reg1, + newdata=base)+summary(reg1)$sigma^2/2)
> round(matrix(base$py,n,n),1) [,1] [,2] [,3] [,4] [,5] [,6] [1,] 2871.2 1091.3 41.7 18.3 7.8 21.3 [2,] 3370.8 1281.2 48.9 21.5 9.2 25.0 [3,] 3768.0 1432.1 54.7 24.0 10.3 28.0 [4,] 5181.5 1969.4 75.2 33.0 14.2 38.5 [5,] 4994.1 1898.1 72.5 31.8 13.6 37.1 [6,] 5297.8 2013.6 76.9 33.7 14.5 39.3 > sum(base$py[is.na(base$y)]) [1] 2481.857 On obtient un montant un peu différent de celui obtenu par la méthode Chain Ladder, mais néanmoins comparable. On peut aussi tenter une régression de Poisson (avec un lien logarithmique), comme suggéré en 1975 par Hachemeister et Stanard, > 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) Deviance Residuals: Min 1Q Median 3Q Max -2.3426 -0.4996 0.0000 0.2770 3.9355 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 La prédiction coïncide avec l’estimateur obtenu par la méthode Chain Ladder. Le lien avec les méthodes de biais minimal a été établi par Klaus Schmidt et Angela Wünsche en 1998, dans Chain ladder, marginal sum and maximum likelihood estimation. La semaine prochaine, on parlera des méthodes de bootstrap pour obtenir des intervalles de confiance, ou des quantiles, sur les montants de réserve. Je ne sais pas si j’aurais le temps de taper des transparents, je préfère, sur cette partie du cours taper au fur et à mesure, et écrire au tableau. Je renvoie au Chapitre 3 du livre avec Christophe Dutang – en ligne sur http://cran.r-project.org/doc/contrib/ – pour le détail. C’est le code que je tape en cours, tout en essayant de répondre aussi aux questions. # 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. # Examen intra, régression logistique et de Poisson L’examen intra du cours ACT2040 aura lieu mercredi matin, de 9:00 à 12:00. Aucun document autorisé, sauf les calculatrices (modèle standard, cf plan de cours), et les téléphones seront formellement interdits. Il y aura 34 questions portant sur la première partie du cours (jusqu’à la fin des modèles de comptage, sections 1 à 5 des transparents). 15 questions porteront sur la base décrite dans un précédant billet, sur le nombre de relations extra-conjugales. Il s’agira de décrire les sorties en ligne ici. Je laisse 36 heures pour prendre connaissance de ces sorties. Une version sera donnée lors de l’examen (imprimée 2 pages par feuille, comme dans la version en ligne: si quelqu’un a besoin d’un exemplaire imprimé plus gros, merci de me le faire savoir avant l’examen). # 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.

# De l’interprétation d’un effet nonlinéaire en régression

Je vais poursuivre un peu le TD de statistique de l’actuariat 2, où nous avions parlé tarification et régression Poissonnienne.

• Analyser un effet nonlinéaire

Dans un premier temps, nous avions regardé un modèle simple, où l’on voulait voir si l’âge du conducteur expliquait la fréquence de sinistres. Bref, on peut utiliser une régression log-Poisson, qui donne la prédiction suivante, en fonction de l’âge,

Mais l’âge est une variable un peu particulière car c’est une fausse variable continue: seule les valeurs entières apparaissent. Comme on chercher à calculer

notons qu’il est possible de considérer l’estimateur empirique naturel de cette grandeur, i.e.

Graphiquement, on obtient la courbe suivante,

Le côté ératique de la série à droite et à à gauche vient de la faible représentation dans la base, comme en atteste l’histogramme de l’âge du conducteur,

Je peux en profiter pour revenir à un commentaire posté ici sur les modèles GAM. En effet, au lieu de regarder la fréquence en fonction de l’âge, on peut s’intéresser à la la différence par rapport à la fréquence globale.

> m=mean(base\$nbsin)
> Y=predict(REG,newdata=data.frame(age=X,expo=E),type = "response")
> lines(X,Y/m-1,col="red")

On obtient alors le graphique suivant (en superposant à l’ajustement obtenu par GAM)

Notons que si l’on ajuste un modèle log-Poisson uniquement sur la tranche d’âge [50-70], on est très proche de ce que propose le modèle GAM, d’où l’interprétation de l’ajustement GAM comme un ajustement local.

• Changer de modèle (et les comparer)

On a choisi un modèle log-Poisson, mais on se doute que d’autres modèles seraient possible. En particulier dans notre base, les assurés ont eu 0 ou 1 sinistre l’année passée. Autrement dit, des lois binomiales peuvent être ajustée,

> REG0 = glm(nbsin~age,base,
+  offset=expo)
> REG1a = glm(nbsin~age,base,
+  offset=expo)
> REG1b = glm(nbsin~age,base,
+  offset=expo)
> REG1c = glm(nbsin~age,base,
+  offset=expo)

Le critère AIC est souvent invoqué pour juger de l’adéquation d’un modèle. Il vaut ici

> AIC(REG0)
[1] 36470.7
> AIC(REG1a)
[1] 36044.45
> AIC(REG1b)
[1] 37304.02
> AIC(REG1c)
[1] 32549.83

Mais on peut également regarder le critère BIC,

> AIC(REG0, k = log(nrow(base)))
[1] 36488.81
> AIC(REG1a, k = log(nrow(base)))
[1] 36062.57
> AIC(REG1b, k = log(nrow(base)))
[1] 37322.14
> AIC(REG1c, k = log(nrow(base)))
[1] 32577

Autrement dit, en utilisant le règle “the smaller, the better” on retiendrait le dernier modèle qui est un modèle Gaussien, manifestement assez mauvais: sur le graphique ci-dessous, la régression log-Poisson est en mauve, presque confondu avec la régression logistique1, en bleu, la régression probit en vert et le modèle linéaire (Gaussien) en rouge,

Bref, ce dernier modèle n’est pas convainquant… Par contre, notons que si l’on compare le modèle log-Poisson et logit, ce dernier conduit à tarifer moins cher les âges extrêmes,

Autrement dit, au lieu de comparer des modèles actuariels sur la base d’outils statistiques, et uniquement de ces outils, il peut être instructif d’essayer de comprendre qui verra sa prime baisser avec tel ou tel modèle… Mais je reviendrais plus longuement, une autre fois, sur les histoires de choix de modèles…..