> mean(BASEN$nombre)/mean(BASEN$exposition)
[1] 0.05241692
> exp(glm(nombre~1,offset=log(exposition),data=BASEN,
+ family=poisson(link="log"))$coefficients)
(Intercept)
0.05241692

(je vais essayer de donner à chaque question deux méthodes pour obtenir le résultat). Si on suppose que le nombre de sinistres par an suit une loi de Poisson, alors la probabilité d’avoir au moins un accident est

> lambda=mean(BASEN$nombre)/mean(BASEN$exposition)
> 1-ppois(0,lambda)
[1] 0.05106684
> 1-exp(-lambda)
[1] 0.05106684

Si on cherche à distinguer par sexe du conducteur principal

> mean(BASEN$nombre[BASEN$sexeconducteur=="H"])/
+ mean(BASEN$exposition[BASEN$sexeconducteur=="H"])
[1] 0.05218143
> mean(BASEN$nombre[BASEN$sexeconducteur=="F"])/
+ mean(BASEN$exposition[BASEN$sexeconducteur=="F"])
[1] 0.05457013
> exp(glm(nombre~0+sexeconducteur,offset=log(exposition),
+ data=BASEN,family=poisson(link="log"))$coefficients)
sexeconducteurF sexeconducteurH
0.05457013      0.05218143

Ces deux valeurs ne sont pas significativement différentes. Si on regarde la sortie de la régression

Coefficients:
                Estimate Std. Error z value Pr(>|z|)    
sexeconducteurF -2.90827    0.05464  -53.23   <2e-16 ***
sexeconducteurH -2.95303    0.01848 -159.82   <2e-16 ***

les deux coefficients sont quasiment identiques (en enlevant la constante). Si on rajoute la constante, la modalité de référence devientle sexe féminin

Coefficients:
                Estimate Std. Error z value Pr(>|z|)    
(Intercept)     -2.90827    0.05464 -53.230   <2e-16 ***
sexeconducteurH -0.04476    0.05768  -0.776    0.438

et la modalité homme n’est alors plus significative. Bref, le sexe du conducteur n’est pas une variable explicative.

Ensuite, parmi les variables qui semble pouvoir etre utilisée pour modéliser le nombre de sinistre, on pourra retenir l’age du conducteur (décroissant jusqu’à 30 ans, puis croissant), l’age du véhicule (décroissant), l’ancienneté du permis de conduire (décroissant les 10 premières années puis croissant – ou indépendant), éventuellement certains type de véhicules (le coupé cabriolet semble avoir – en moyenne – plus d’accident que la Jeep), la situation familiale, le lieu d’habitation (en particulier rural semble avoir moins d’accident qu’urbain), le statut d’habitation, le type de paiement de la prime (mensuel semblant plus risqué qu’annuel), le poids du véhicule, ou la marque (ou plutot certaines marques, comme Peugeot versus Burstner Mobil). Bref, toutes, sauf le sexe du conducteur.

Pour un conducteur de 45 ans qui a son permis de conduire depuis seulement 5 ans, les trois modèles donnent la fréquence annuelle suivante

> nbase=data.frame(ageconducteur=45,agepermis=5,exposition=1)
> exp(t(c(1,5))%*%reg6$coefficients)
[,1]
[1,] 0.03395972
> predict(reg6,newdata=nbase,type="response")
1
0.03395972
> exp(t(c(1,45))%*%reg7$coefficients)
[,1]
[1,] 0.04100356
> predict(reg7,newdata=nbase,type="response")
1
0.04100356
> exp(t(c(1,5,45))%*%reg8$coefficients)
[,1]
[1,] 0.05395398
> predict(reg8,newdata=nbase,type="response")
1
0.05395398

Les trois modèles sont respectivement un modèle de Poison avec un lien log, un modèle binomiale négative avec un lien log, et un modèle binomial, avec un lien logit. Les deux premiers sont très proches: ils ont la meme fonction de lien, et comme nous l’avions vu à l’examen intra, le modèle binomiale négative est un modèle de Poisson avec une variable latente non observée qui suivrait une loi Gamma. Donc pour ces deux modèles, si les classes sont homogènes, et que par classe, la loi de Poisson est adaptée, il est normal d’avoir des coefficients très proches. Pour le modèle binomial avec un lien logistique, il faut se souvenir que la fonction lien est

car ici la probabilité d’avoir un accident est faible. Donc les deux fonctions liens sont proches. De plus, pour une loi de Poisson

Bref, pour notre loi de Poisson, on a presque une loi binomiale de paramètre si est petit.

Pour la prédiction suivant les trois modèles, on utilise classiquement

> exp(t(c(1,40,0,0,1,0,1,0))%*%reg9$coefficients)
           [,1]
[1,] 0.04351037

Comme on a accès aux sorties, pour les trois modèles, les trois prédictions sont très proches (mais on s’en doutait d’après la question précédante)

> predict(reg9,newdata=nbase,type="response")
         1 
0.04351037 
> predict(reg10,newdata=nbase,type="response")
         1 
0.04351056 
> predict(reg11,newdata=nbase,type="response")
         1 
0.04352741

Pour la méthode de biais minimal, il s’agit du principe de balancement, qui coïncide avec la régression de Poisson sur les classes. Pour les régressions stepwise, dans un cas, on élimine (étape par étape) des variables, et dans l’autre cas, des modalités de variables (pour les variables qualitatives).

Les deux modèles semblent valides, au sens où toutes les variables sont significatives. Le découpage en classes pour l’age du véhicule semble justifié par l’arbre de régression, qui suggèrait de découper à 8 et 14 ans. Pour les prédictions, on a pour le premier modèle

> nbase15=data.frame(ageconducteur=40,agepermis=20,
+ habitation="urbain",locataire=TRUE,paymensuel=FALSE,
+ marie=TRUEFORD=TRUE,RENAULT=FALSE,PEUGEOT=FALSE,
+ exposition=1)
> predict(reg15,newdata=nbase15,type="response")
         1 
0.03398312

alors que le second donne

> nbase16=data.frame(ageconducteur=40,agevehF="NEUF",rural=FALSE,
+ locataire=TRUE,paymensuel=FALSE,marie=TRUE,
+ jeune=FALSE,ancien=FALSE,exposition=1)
> predict(reg16,newdata=nbase16,type="response")
         1 
0.05208656

Moralité, nos deux modèles (qui semblent autant valides l’un que l’autre) donnent des prédictions relativement différentes. Si l’assuré change de voiture, cela n’a pas d’influence pour le second modèle (le véhicule est toujours dans la catégorie neuf). Pour le premier, il perd le rabais qu’il avait en conduisant un véhicule Ford, et donc la fréquence annuelle de sinistre augmente de 23.5%

> 1/exp(reg15$coefficients["FORDTRUE"])
FORDTRUE 
1.235329

Bref, la prédiction devient

> predict(reg15,newdata=nbase15,type="response")/
+ exp(reg15$coefficients["FORDTRUE"])
         1 
0.04198033

On peut d’ailleurs le vérifier simplement,

> nbase15b=data.frame(ageconducteur=40,agepermis=20,
+ habitation="urbain",locataire=TRUE,paymensuel=FALSE,
+ marie=TRUE,FORD=FALSE,RENAULT=FALSE,PEUGEOT=FALSE,
+ exposition=1)
> predict(reg15,newdata=nbase15b,type="response")
         1 
0.04198033

Si on passe maintenant à la modélisation de la charge de sinistres, le coût moyen par sinistre serait

> mean(BASEY$cout)
[1] 1468.522

Cette fois, contrairement à la modélisation de la fréquence, il est plus difficile de trouver des variables explicatives. A priori, au vu des graphiques, aucune variable ne semblerait convenir…
Si on regarde la prédiction pour le premier modèle obtenu,

> nbase7=data.frame(VOLKSWAGEN=FALSE,agepermis=20,
+ sexeconducteur="H",usage="TOUS_DEPLACEMENTS")
> predict(regr7,newdata=nbase7,type="response")
       1 
1388.622

alors que pour le second,

> exp(predict(regr8,newdata=nbase7))*
+ exp(summary(regr8)$sigma^2/2)
       1 
1359.552

Si c’est la femme de l’assuré qui avait pris l’assurance, cela n’aurait rien changé pour le second modèle. En revanche, pour le premier, l’assuré perdrait le rabais qu’il avait en étant un homme, et donc le cout moyen attendu augmenterait de 23%,

> 1/exp(regr7$coefficients["sexeconducteurH"])
sexeconducteurH 
       1.237832

i.e.

> predict(regr7,newdata=nbase7,type="response")/
+ exp(regr7$coefficients["sexeconducteurH"])
       1 
1718.88

La encore, on retrouve le meme résultat directement

> nbase7b=data.frame(VOLKSWAGEN=FALSE,agepermis=20,
+ sexeconducteur="F",usage="TOUS_DEPLACEMENTS")
> predict(regr7,newdata=nbase7b,type="response")
       1 
1718.88

Passons maintenant aux calculs de la prime pure. Sans segmentation, la prime pure serait

> sum(BASEN$nombre)/sum(BASEN$exposition)*
+ mean(BASEY$cout)
[1] 76.9754

Les prédictions des modèles pour la fréquence donnaient

> (n1=predict(reg15,newdata=nbase15,type="response"))
         1 
0.03398312 
> (n2=predict(reg16,newdata=nbase16,type="response"))
         1 
0.05208656

alors que pour les couts moyens par sinistre, nous avions

> (c1=predict(regr7,newdata=nbase7,type="response"))
       1 
1388.622 
> (c2=exp(predict(regr8,newdata=nbase7))
+ *exp(summary(regr8)$sigma^2/2))
       1 
1359.552

Les quatre primes pures obtenues sont alors

> c(n1,n2)%*%t(c(c1,c2))
            1        1
[1,] 47.18972 46.20183
[2,] 72.32855 70.81440

Passons maintenant à la seconde partie – sur le provisionnement des sinistres à payer. On commence par travailler sur le triangle

TRIANGLE2010
          D1      D2      D3      D4      D5      D6      D7      D8      D9
A2002 357848 1124788 1735330 2218270 2745596 3319994 3466336 3606286 3833515
A2003 352118 1236139 2170033 3353322 3799067 4120063 4647867 4914039 5339085
A2004 290507 1292306 2218525 3235179 3985995 4132918 4628910 4909315      NA
A2005 310608 1418858 2195047 3757447 4029929 4381982 4588268      NA      NA
A2006 443160 1136350 2128333 2897821 3402672 3873311      NA      NA      NA
A2007 396132 1333217 2180715 2985752 3691712      NA      NA      NA      NA
A2008 440832 1288463 2419861 3483130      NA      NA      NA      NA      NA
A2009 359480 1421128 2864498      NA      NA      NA      NA      NA      NA
A2010 376686 1363294      NA      NA      NA      NA      NA      NA      NA

Pour obtenir le montant payé en 2008 pour les sinistres survenus en 2004, on utilise la différence entre

> TRIANGLE2010["A2004",c("D4","D5")]
           D4      D5
A2004 3235179 3985995

i.e.

> TRIANGLE2010["A2004","D5"]-TRIANGLE2010["A2004","D4"]
[1] 750816

Pour le montant payé au total en 2005, on utilise

> CUMULPAIEMENTS2010=TRIANGLE2010
> INCREMENTPAIEMENTS2010=CUMULPAIEMENTS2010
> INCREMENTPAIEMENTS2010[,2:9]=CUMULPAIEMENTS2010[,2:9]-
+ CUMULPAIEMENTS2010[,1:8]
> sum(diag(as.matrix(INCREMENTPAIEMENTS2010[1:4,4:1])))
[1] 2729241

On attaque maintenant la méthode Chain Ladder. Les coefficients de développement sont

> lambda=rep(NA,8)
> for(k in 1:8){
+ lambda[k]=(sum(TRIANGLE2010[1:(9-k),k+1])/
+ sum(TRIANGLE2010[1:(9-k),k]))}
> lambda[1:5]
[1] 3.474193 1.704149 1.460866 1.161765 1.095763

On peut alors obtenu ce qui sera payé en 2011 pour les sinistres survenus en 2010,

> TRIANGLE2010["A2010","D1"]*(lambda[1]-1)
[1] 931993.8

alors que pour le montant payé en 2011 pour les sinistres survenus en 2009,

> TRIANGLE2010["A2009","D2"]*(lambda[2]-1)
[1] 1000686

Pour le montant payé en 2010 et 2011 pour les sinistres survenus en 2010,

> TRIANGLE2010["A2010","D1"]*lambda[1]*lambda[2]
[1] 2230186

Pour les sinistres survenus en 2008, la charge total estimée sera

> TRIANGLE2010["A2008","D3"]*prod(lambda[3:8])
[1] 5531128

ce qui fait un montant de provision

> TRIANGLE2010["A2008","D3"]*(prod(lambda[3:8])-1)
[1] 3111267

On rajoute maintenant une diagonale dans notre triangle, i.e. une nouvelle année de paiements (en 2011).
Le montant payé en 2011 pour les sinistres survenus en 2010 est

> TRIANGLE2011["A2010","D2"]-TRIANGLE2011["A2010","D1"]
[1] 986608

à comparer avec

> TRIANGLE2010["A2010","D1"]*(lambda[1]-1)
[1] 931993.8

prédit par la méthode Chain-Ladder, soit une erreur de l’ordre de 5%,

> (TRIANGLE2011["A2010","D2"]-
+ TRIANGLE2011["A2010","D1"])/
+ (TRIANGLE2010["A2010","D1"]*(lambda[1]-1))-1
[1] 0.05859927

Si on compare maintenant les coefficients de transition de la méthode Chain Ladder,

> lambda2011=rep(NA,9)
> for(k in 2:9){
+ lambda2011[k]=(sum(TRIANGLE2011[1:(10-k),k+1])/
+ sum(TRIANGLE2011[1:(10-k),k]))}
> lambda2011[1:5]
[1]       NA 1.747333 1.457413 1.173852 1.103824
> lambda[1:5]
[1] 3.474193 1.704149 1.460866 1.161765 1.095763

Le coefficient permettant de passer de la colonne 2 à la colonne 3 augmente (légèrement). Si on regarde la prédiction du montant de paiements effectué en 2012 pour les sinistres survenus en 2009,

> TRIANGLE2011["A2009","D3"]*(lambda2011[3]-1)
[1] 1310258

qui sont à comparer avec les

> TRIANGLE2010["A2009","D2"]*lambda[2]*(lambda[3]-1)
[1] 1116132

estimés un an auparavant. La variation entre les estimations est de l’ordre de 17%,

> TRIANGLE2011["A2009","D3"]*(lambda2011[3]-1)/
+ (TRIANGLE2010["A2009","D2"]*lambda[2]*(lambda[3]-1))-1
[1] 0.1739278

Entre 2010 et 2012, on peut estimer que l’on payera

> TRIANGLE2011["A2010","D2"]*lambda2011[2]
[1] 2382128

pour les sinistres survenus en 2010.
Si on se focalise sur les sinistres survenus en 2008, la prédiction de la charge finale est ici

> TRIANGLE2011["A2008","D4"]*prod(lambda2011[4:9])
[1] 5660771

à comparer avec la prédiction faire en 2010,

> TRIANGLE2010["A2008","D3"]*prod(lambda[3:8])
[1] 5531128

i.e. notre nouvelle prédiction est 2% plus élevée,

> TRIANGLE2011["A2008","D4"]*prod(lambda2011[4:9])/
+ (TRIANGLE2010["A2008","D3"]*prod(lambda[3:8]))-1
[1] 0.02343867

Le provisions pour sinistres à payer pour les sinistres survenus en 2008 est estimé, fin 2011, à

> TRIANGLE2011["A2008","D4"]*(prod(lambda2011[4:9])-1)
[1] 2177641

alors que fin 2010, un montant de provision avait été constitué, à hauteur de

> TRIANGLE2010["A2008","D3"]*(prod(lambda[3:8])-1)
[1] 3111267

On notera que la différence

> TRIANGLE2011["A2008","D4"]*(prod(lambda2011[4:9])-1)-
+ TRIANGLE2010["A2008","D3"]*(prod(lambda[3:8])-1)
[1] -933626.7

vient des paiements effectués en 2011, et du changement dans l’estimation de charge ultime

> TRIANGLE2011["A2008","D4"]*prod(lambda2011[4:9])-
+ (TRIANGLE2010["A2008","D3"]*prod(lambda[3:8]))
[1] 129642.3
> TRIANGLE2011["A2008","D4"]-TRIANGLE2011["A2008","D3"]
[1] 1063269
> TRIANGLE2011["A2008","D4"]*prod(lambda2011[4:9])-
+ (TRIANGLE2010["A2008","D3"]*prod(lambda[3:8]))-
+ (TRIANGLE2011["A2008","D4"]-TRIANGLE2011["A2008","D3"])
[1] -933626.7

On attaque enfin la dernière partie, avec la régression de Poisson. La prédiction du montant de paiements faits en 2012 pour les sinistres survenus en 2009 est ici

> Y=as.vector(as.matrix(INCREMENTPAIEMENTS2011))
> ANNEE=rep(2002:2010,10)
> DEVEL=rep(1:10,each=9)
> baseTriangle=data.frame(Y,A=as.factor(ANNEE),
+ D=as.factor(DEVEL))
> reg=glm(Y~A+D,data=baseTriangle,famil=poisson(link="log"))
> nbase=data.frame(A="2009",D="4")
> predict(reg,newdata=nbase,type="response")
      1 
1310258

qui correspond à ce que donnait la méthode Chain Ladder

> TRIANGLE2011["A2009","D3"]*(lambda2011[3]-1)
[1] 1310258

(mais on savait que les deux modèles prédisaient la même chose). Pour les sinistres survenus en 2005, le montant de provisions à constituer est

> nbase=data.frame(A="2005",D=c("8","9","10"))
> predict(reg,newdata=nbase,type="response")
       1        2        3 
247190.0 370179.3  92268.5 
> sum(predict(reg,newdata=nbase,type="response"))
[1] 709637.8

qui, là encore, correspond à ce que prédisait la méthode Chain Ladder.

> TRIANGLE2011["A2005","D7"]*(prod(lambda2011[7:9])-1)
[1] 709637.8