Comme annoncé, l’examen terminal a lieu ce matin, et il est basé sur des données mises en ligne sur le blog depuis une dizaine de jours. L’énoncé est un peu longue, avec une bonne trentaine de pages de sorties informatiques, et un énoncécomprenant une trentaine de questions. Des éléments de correction sont en ligne ci-dessous (plutôt que de taper une correction en pdf, autant mettre le code en ligne).
> 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
(classique en statistique quand on possède 
observations)
. D’où l’utilisation des résidus dits ajustés afin d’avoir des résidus de variance (vraiment) unitaire.
, mais aussi que
. On a vu que le paramètre de surdispersion ne pouvait être supposé unitaire. On a alors le choix, entre simuler une loi Gamma de paramètres
et 
, et de paramètre de surdispersion 
telle que la variance s’écrive
alors que
On voit que l’on risque d’avoir un incrément négatif si le résidu est négatif, et plus grand (en valeur absolue) que la racine carrée de la prédiction (par notre modèle Poissonnien, en l’occurrence). Or dans le triangle vu en cours, on a
L’examen final pour le cours d’actuariat IARD aura lieu mardi prochain, et il sera basé sur des données mises 
sept mots, le journaliste a réussi a utiliser trois clichés qui m’agacent.
En relisant les devoirs maisons, je me suis rendu compte que certains avaient tenté de regrouper les régions (géographiques) par régions homogènes. Sauf que les régions étaient codées par un numéro (selon la 
et 
et 
on en déduit
