Freakonometrics

To content | To menu | To search

Teaching › ACT2040-2011

Entries feed - Comments feed

Tuesday, December 13 2011

ACT2040: examen final (suite et fin)

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).

On suppose importée les vases (sinon, le code est ici). Pour la première partie - sur la tarification (a priori) - la fréquence annuelle de sinistres est tout simplement

Continue reading...

Thursday, December 8 2011

ACT2040: quels résidus et quelle loi simuler

Suite à une question sur les résidus, je vais essayer de prendre deux minutes pour reprendre un point que je n'avais pas abordé en cours, faute de temps. Les résidus les plus intéressants en provisionnement (et dans les GLM) sont les résidus de Pearson,

http://freakonometrics.blog.free.fr/public/perso4/res-boot-01.gif

Mais il existe également des résidus dits ajustés, permettant de prendre en compte le fait que le nombre de paramètres est ici relativement grand, ramené au nombre d'observations. On pose alors

http://freakonometrics.blog.free.fr/public/perso4/res-boot-2.gif

A première vue, utiliser l'un ou l'autre devrait donner la même chose si l'on génère des pseudo-triangles, puisque dans un cas, on utiliserait

http://freakonometrics.blog.free.fr/public/perso4/res-boot.gif

et dans l'autre

http://freakonometrics.blog.free.fr/public/perso4/res-boot-3.gif

Bref, multiplier par une constante pour ensuite diviser par la même constante, ça revient au même. Sauf qu'il peut sembler légitime de poser

http://freakonometrics.blog.free.fr/public/perso4/res-boot-05.gif
(on reviendra là dessus tout à l'heure). En attendant, voilà nos deux résidus
> (E=residuals(regp,"pearson"))
1             2             3             4
9.488238e-01  2.404895e-02  1.168421e-01 -1.082940e+00
5             6             7             8
1.302749e-01 -1.007348e-13 -1.128013e+00  2.773332e-01
9            10            11            13
5.669707e-02  8.919633e-01 -2.110748e-01 -1.533031e+00
14            15            16            19
-2.213449e+00 -1.024162e+00  4.237393e+00 -4.899687e-01
20            21            25            26
7.929194e-01 -2.972380e-01 -4.275912e-01  4.140426e-01
31
-6.202125e-15
> n=sum(is.na(Y)==FALSE)
> k=ncol(PAID)+nrow(PAID)-1
> (R=residuals(regp,"pearson")*sqrt(n/(n-k)))
1             2             3             4
1.374976e+00  3.485024e-02  1.693203e-01 -1.569329e+00
5             6             7             8
1.887862e-01 -1.459787e-13 -1.634646e+00  4.018940e-01
9            10            11            13
8.216186e-02  1.292578e+00 -3.058764e-01 -2.221573e+00
14            15            16            19
-3.207593e+00 -1.484151e+00  6.140566e+00 -7.100321e-01
20            21            25            26
1.149049e+00 -4.307387e-01 -6.196386e-01  6.000048e-01
31
-8.987734e-15 

L'autre point sur lequel je voulais revenir sur les simulations, est que la loi de Poisson n'est peut-être pas adaptée pour simuler des scénarios de paiements futurs. En effet, si on fait une régression quasipoisson, on voit que le paramètre de surdispersion n'est pas négligeable,

> regqp=glm(Y~as.factor(D)+as.factor(A),
+     data=base,family=quasipoisson(link="log"))
> summary(regqp)
 
Call:
glm(formula = Y ~ as.factor(D) + as.factor(A), 
family = quasipoisson(link = "log"), data = base)   Deviance Residuals: Min 1Q Median 3Q Max -2.3426 -0.4996 0.0000 0.2770 3.9355   Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 8.05697 0.02769 290.995 < 2e-16 *** as.factor(D)2 -0.96513 0.02427 -39.772 2.41e-12 *** as.factor(D)3 -4.14853 0.11805 -35.142 8.26e-12 *** as.factor(D)4 -5.10499 0.22548 -22.641 6.36e-10 *** as.factor(D)5 -5.94962 0.43338 -13.728 8.17e-08 *** as.factor(D)6 -5.01244 0.39050 -12.836 1.55e-07 *** as.factor(A)2002 0.06440 0.03731 1.726 0.115054 as.factor(A)2003 0.20242 0.03615 5.599 0.000228 *** as.factor(A)2004 0.31175 0.03535 8.820 4.96e-06 *** as.factor(A)2005 0.44407 0.03451 12.869 1.51e-07 *** as.factor(A)2006 0.50271 0.03711 13.546 9.28e-08 *** --- 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

On peut alors envisager de simuler non pas des lois de Poisson, mais des lois quasipoissons, en reprenant un vieux code,

> 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)
+ }

A partir de ces deux remarques, on peut reprendre le code qui permettait de générer des scénarios de paiements pour les années futures,

> Yp=predict(regp,type="response",newdata=base)
>
Rs=rep(NA,20000) > for(s in 1:20000){ + serreur=sample(erreur, + size=36,replace=TRUE) + E=matrix(serreur,6,6) + sY=matrix(Yp,6,6)+E*sqrt(matrix(Yp,6,6)) + sbase=data.frame(sY=as.vector(sY),D,A) + sbase$sY[is.na(Y)==TRUE]=NA + sreg=glm(sY~as.factor(D)+as.factor(A), + data=sbase,family=poisson(link="log")) + sYp=predict(sreg,type="response", + newdata=sbase) + sYpscenario=rqpois(36,sYp,phi=3.18623) + Rs[s]=sum(sYpscenario[is.na(Y)==TRUE]) + }
si on regarde la densité de la distribution des paiements futurs, on obtient,

avec en bleu la distribution obtenue sur les résidus de Pearson brut, et en simulant une loi de Poisson, et en rouge la distribution des provisions avec les deux modifications. On observe clairement que les quantiles (par exemple) ont fortement changé. On peut le voir encore plus précisément sur un box-plot

Si on revient cinq minutes sur ce qu'on vient de faire. Dans le premier rappelons que l'on a besoin d'utiliser, pour prédire les paiements futurs, que
http://freakonometrics.blog.free.fr/public/perso4/siiiim-01.gif
(avec la filtration correspondant à la partie supérieure du triangle disponible), i.e.
http://freakonometrics.blog.free.fr/public/perso4/siiiiim-03.gif
Comme on utilise un modèle de Poisson, on suppose aussi que
http://freakonometrics.blog.free.fr/public/perso4/siiiim-5.gif
i.e.
http://freakonometrics.blog.free.fr/public/perso4/siiiiim-08.gif
c'est à dire que
http://freakonometrics.blog.free.fr/public/perso4/siiiim-04.gif
Il faut avoir des résidus centrés et de variance unitaire. A première vue, c'est le cas pour nos résidus de Pearson (à peu près),
> mean(residuals(regp,"pearson"))
[1] -0.02462518
> sd(residuals(regp,"pearson"))
[1] 1.261934
sauf que l'estimateur de la variance est biaisé. Certes, classiquement, l'estimateur est basé sur une normalisation par un facteur http://freakonometrics.blog.free.fr/public/perso4/siiiim-11.gif (classique en statistique quand on possède http://freakonometrics.blog.free.fr/public/perso4/siiim-10.gif observations)
> sqrt(sum((residuals(regp,"pearson")-
+ mean(residuals(regp,"pearson")))^2)/
+ (length(residuals(regp,"pearson"))-1))
[1] 1.261934
Sauf qu'ici on ne perd pas un degré de liberté (remplacer l'espérance des observations par leur moyenne empirique): en régression, il faut corriger par le nombre de variables explicatives. Et donc, pour avoir un estimateur sans biais de la variance de nos résidus, il faut corriger par un facteur http://freakonometrics.blog.free.fr/public/perso4/siiiim-12.gif. D'où l'utilisation des résidus dits ajustés afin d'avoir des résidus de variance (vraiment) unitaire.
Pour le second point, on utilise le fait qu'en pratique, la dispersion des paiements est plus grande que ce que donnerait un modèle de Poisson. Et il convient d'en tenir compte dans nos simulations. En l'occurrence, on veut que  http://perso.univ-rennes1.fr/arthur.charpentier/latex/qp001.png, mais aussi quehttp://perso.univ-rennes1.fr/arthur.charpentier/latex/qp002.png. 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
http://perso.univ-rennes1.fr/arthur.charpentier/latex/qp003.png et http://perso.univ-rennes1.fr/arthur.charpentier/latex/qp004.png
ou bien simuler une loi binomiale négative, de moyenne http://perso.univ-rennes1.fr/arthur.charpentier/latex/qp008.png, et de paramètre de surdispersion http://perso.univ-rennes1.fr/arthur.charpentier/latex/qp009.png telle que la variance s'écrive
http://perso.univ-rennes1.fr/arthur.charpentier/latex/qp010.png
Aussi, on prend http://perso.univ-rennes1.fr/arthur.charpentier/latex/qp011.png alors que
http://perso.univ-rennes1.fr/arthur.charpentier/latex/qp012.png
Avec ces méthodes, on obtient quelque chose... qui est programmé sous R,
> BootChainLadder(PAID,10000,"od.pois")
BootChainLadder(Triangle = PAID, R = 10000, 
process.distr = "od.pois")   Latest Mean Ultimate Mean IBNR SD IBNR IBNR 75% IBNR 95% 1 4,456 4,456 0.0 0.0 0 0 2 4,730 4,752 22.2 12.2 28 44 3 5,420 5,456 35.5 15.2 43 64 4 6,020 6,086 65.7 19.8 78 102 5 6,794 6,946 151.9 28.7 170 202 6 5,217 7,364 2,147.4 111.0 2,218 2,339   Totals Latest: 32,637 Mean Ultimate: 35,060 Mean IBNR: 2,423 SD IBNR: 132 Total IBNR 75%: 2,506 Total IBNR 95%: 2,653   > quantile(Rs,c(.75,.95)) 75% 95% 2509 2653
et qui donne les mêmes résultats que ce qu'on vient de reprogrammer...

Wednesday, December 7 2011

ACT2040: tail factor et cloture de la première ligne des triangles de paiements

La grosse hypothèse que l'on a fait quand on travaillait sur les triangles, était que la première année était close. Et donc, que pour la première ligne de paiements cumulés, le montant final était la charge totale des sinistres survenus cette année là. En fait, en 1999, Thomas Mack avait proposé un modèle relativement simple permettant l'inclusion d'un tail factor. L'idée est que l'on peut supposer que les facteurs de développement dans la méthode Chain Ladder décroissent tranquillement vers 1. Toujours sur le même triangle (mais en utilisant une régression pondérée sans constante - ce qui permet d'éviter les problèmes des valeurs manquantes à n'importe quel endroit dans le triangle), nous avions

> LAMBDA <- rep(NA,nc-1)
> for(k in 1:(nc-1)){
+ LAMBDA[k]=lm(PAID[,k+1]~0+PAID[,k],
+ weights=1/PAID[,k])$coefficients}
> LAMBDA
[1] 1.380933 1.011433 1.004343 1.001858 1.004735
Par "tranquillement", on entend une décroissance exponentielle des coefficients avec le temps. On peut estimer un modèle linéaire sur le logarithme des facteurs de développement
> logL <- log(LAMBDA-1)
> tps <- 1:(nc-1)
> modele <- lm(logL~tps)
> plot(tps,logL,xlim=c(1,20),ylim=c(-30,0),col="red")
> abline(modele)
> tpsP <- seq(6,1000)
> logP <- predict(modele,newdata=data.frame(tps=tpsP))
> points(tpsP,logP ,pch=0,col="blue")

Si on regarde les coefficients de transition, on a

Pour passer de la dernière colonne observé à la charge ultime, il suffit de faire le produit des facteurs de transitions pour les années non encore observées, c'est à dire les points bleus
> (facteur <- prod(exp(logP)+1))
[1] 1.000707
Autrement dit, il faudrait rajouter 0.07% à la charge ultime. Sans ce facteur ultime, le montant de provisions était
> DIAG <- diag(PAID[,6:1])
> PRODUIT <- c(1,rev(LAMBDA))
> sum((cumprod(PRODUIT)-1)*DIAG)
[1] 2426.985
mais si on rajoute le facteur que l'on vient de calculer, on passe à
> sum((cumprod(PRODUIT)*facteur-1)*DIAG)
[1] 2451.764
On retrouve ici la même chose que ce qui est obtenu sous R,
> MackChainLadder(PAID,tail=TRUE)
MackChainLadder(Triangle = PAID, tail = TRUE)
 
Latest Dev.To.Date Ultimate     IBNR Mack.S.E CV(IBNR)
1  4,456       0.999    4,459     3.15    0.299   0.0948
2  4,730       0.995    4,756    25.76    0.712   0.0277
3  5,420       0.993    5,460    39.64    2.528   0.0638
4  6,020       0.988    6,090    70.37    5.064   0.0720
5  6,794       0.977    6,952   157.99   31.357   0.1985
6  5,217       0.708    7,372 2,154.86   68.499   0.0318
 
Totals
Latest:             32,637.00
Dev:                     0.93
Ultimate:           35,088.76
IBNR:                2,451.76
Mack S.E.:              79.37
CV(IBNR):  0.0323730538389919

ACT2040: triangles de paiements, ou charge dossier par dossier

Comme on l'avait vu en début de cours, on a en fait souvent plus d'information que les paiements effectués. On a aussi les estimations de charge effectuées par les gestionnaires de sinistres, pour les sinistres présents dans la base (i.e. les sinistres déclarés). Par exemple, sur le jeu de données utilisé en cours, on a les paiements, et l'estimation de la charge dossier par dossier,

> PAID
     [,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
> INCURRED
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,] 4975 4629 4497 4470 4456 4456
[2,] 5135 4949 4783 4760 4750   NA
[3,] 5681 5631 5492 5470   NA   NA
[4,] 6272 6198 6131   NA   NA   NA
[5,] 7326 7087   NA   NA   NA   NA
[6,] 7353   NA   NA   NA   NA   NA

On peut alors utiliser la méthode Chain Ladder non seulement sur les paiements, mais aussi sur les estimation de charge dossier par dossier. Visuellement, on retrouve des courbes assez proches de celles vues en cours (certes, avec des développements mensuels),

> PCL$FullTriangle
1        2        3        4        5        6
1 3209 4372.000 4411.000 4428.000 4435.000 4456.000
2 3367 4659.000 4696.000 4720.000 4730.000 4752.397
3 3871 5345.000 5398.000 5420.000 5430.072 5455.784
4 4239 5917.000 6020.000 6046.147 6057.383 6086.065
5 4929 6794.000 6871.672 6901.518 6914.344 6947.084
6 5217 7204.327 7286.691 7318.339 7331.939 7366.656
> ICL$FullTriangle
1        2        3        4        5        6
1 4975 4629.000 4497.000 4470.000 4456.000 4456.000
2 5135 4949.000 4783.000 4760.000 4750.000 4750.000
3 5681 5631.000 5492.000 5470.000 5455.777 5455.777
4 6272 6198.000 6131.000 6101.117 6085.253 6085.253
5 7326 7087.000 6920.146 6886.416 6868.510 6868.510
6 7353 7129.075 6961.230 6927.300 6909.288 6909.288
> library(ChainLadder)
> PCL <- MackChainLadder(PAID)
> ICL <- MackChainLadder(INCURRED)
> k <- 5
> plot(0:(nc+1-k),c(0,PCL$FullTriangle[k,1:(nc+1-k)]),
+ pch=19,type="b",
+ ylim=c(0,max(c(PCL$FullTriangle[k,],ICL$FullTriangle[k,]))),
+ xlim=c(0,nc), ylab="",xlab="",col="blue")
> lines(0:(nc+1-k),c(0,ICL$FullTriangle[k,1:(nc+1-k)]),
+ pch=19,type="b",col="red")
> lines((nc+1-k):nc,PCL$FullTriangle[k,(nc+1-k):nc],
+ pch=1,type="b",col="blue")
> lines((nc+1-k):nc,ICL$FullTriangle[k,(nc+1-k):nc],
+ pch=1,type="b",col="red")

Il existe une méthode combinant les deux triangles, appelées Munich Chain Ladder, dont la théorie a été développée par Gerhard Quarg et Thomas Mack. Numériquement, on a

> (MNCL <- MunichChainLadder(Paid=PAID,Incurred=INCURRED))
MunichChainLadder(Paid = PAID, Incurred = INCURRED)
 
Latest Paid Latest Incurred Latest P/I Ratio Ult. Paid Ult.
1  4,456      4,456     1.000    4,456      4,456
2  4,730      4,750     0.996    4,753      4,750
3  5,420      5,470     0.991    5,455      5,454
4  6,020      6,131     0.982    6,086      6,085
5  6,794      7,087     0.959    6,983      6,980
6  5,217      7,353     0.710    7,538      7,533
 
Totals
Paid Incurred P/I Ratio
Latest:   32,637   35,247      0.93
Ultimate: 35,271   35,259      1.00
> k <- 5
> plot(0:(nc+1-k),c(0,MNCL$MCLPaid[k,1:(nc+1-k)]),
+ pch=19,type="b", ylim=c(0,max(c(MNCL$MCLPaid[k,],
+ MNCL$MCLIncurred[k,]))),xlim=c(0,nc),
+ ylab="",xlab="",col="blue")
> lines(0:(nc+1-k),c(0,MNCL$MCLIncurred[k,1:(nc+1-k)]),
+ pch=19,type="b",col="red")
> lines((nc+1-k):nc,MNCL$MCLPaid[k,(nc+1-k):nc],
+ pch=1,type="b",col="blue")
> lines((nc+1-k):nc,MNCL$MCLIncurred[k,(nc+1-k):nc],
+ pch=1,type="b",col="red")

ACT2040: simulations et incréments négatifs

Pour faire suite (rapidement) à mon dernier billet, le cas le plus fréquent pour observer des incréments négatifs est lorsque l'on fait des simulations pour obtenir la distribution du montant de paiements futurs. Rappelons que pour générer un pseudo-triangle, on utilise les résidus de Pearson, et on pose

http://freakonometrics.blog.free.fr/public/perso4/pseudo-triangle.gif

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

> source("http://perso.univ-rennes1.fr/
arthur.charpentier/bases.R")
> INC=PAID
> INC[,2:6]=PAID[,2:6]-PAID[,1:5]
> 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),
+     data=base,family=poisson(link="log"))
> Yp=predict(reg,type="response",
+ newdata=base)
> erreurs=residuals(reg,"pearson")
> min(sqrt(Yp[is.na(Y)==FALSE]))
[1] 2.868171
> min(erreurs)
[1] -2.213449

autrement dit, on n'aura jamais d'incrément négatif lors de nos boucle, si l'on génère les résidus par bootstrap. Sinon, avec des lois paramétriques non-bornées inférieurement (loi normale, loi de Student), il est tout a fait possible d'obtenir des incréments négatifs. La méthode pour éviter le problème des incréments négatifs dans les boucles... est probablement de ne pas prendre en compte les scénarios où des incréments négatifs ont été obtenus, i.e.

R=rep(NA,10000)
for(s in 1:10000){
serreur=sample(erreurs,
size=36,replace=TRUE)
E=matrix(serreur,6,6)
sY=matrix(Yp,6,6)+E*sqrt(matrix(Yp,6,6))
if(min(sY[is.na(Y)==FALSE])>=0){
sbase=data.frame(sY=as.vector(sY),D,A)
sbase$sY[is.na(Y)==TRUE]=NA
sreg=glm(sY~as.factor(D)+as.factor(A),
data=sbase,family=poisson(link="log"))
sYp=predict(sreg,type="response",
newdata=sbase)
R[s]=sum(sYp[is.na(Y)==TRUE])}
}

Une stratégie un peu plus propre pourrait être de mettre des 0 dès qu'on obtient un incrément négatif... Mais encore une fois, c'est de la cuisine.

ACT2040: incréments négatifs dans les triangles

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 paiement né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)
MackChainLadder(Triangle = 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),
+     data=base,family=poisson(link="log"))
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 courbe noire 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,
http://freakonometrics.blog.free.fr/public/perso4/glm-translation.gif
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),
+ data=base,family=poisson(link="log"))
> 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),
+ data=base,family=poisson(link="log"))
+ 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),
+ data=base,family=poisson(link="log"))
+ 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...

Tuesday, December 6 2011

ACT2040: examen final (suite)

L'examen final pour le cours d'actuariat IARD aura lieu mardi prochain, et il sera basé sur des données mises en ligne sur le blog depuis une dizaine de jours. Je mets aujourd'hui en ligne les sorties qui seront fournies le jour de l'examen, et sur lesquelles des questions seront posées. J'encourage les étudiants à survoler ces sorties avant l'examen. En revanche, je ne répondrai à aucune question sur ce document ! Amenez vos calculatrices à l'examen.

ACT2040: de la robustesse des arbres

Pour répondre à une question rapide, la robustesse des arbres dépend de deux choses: du nombre d'observations (plus il y a d'observations, plus les nœuds seront stables), et de l'hétérogénéité entre les classes (plus les classes sont différentes, plus il sera facile de les distinguer). Pour faire simple, considérons un petit modèle binomial à trois classes (avec une unique variable explicative), par exemple si on simule 1000 valeurs,

m=c(.2,.5,.3)
moyenne=function(x){
return(m[1]*(x<3)+m[2]*((x>=3)&(x<=7))+m[3]*(x>7))
}
MY=Vectorize(moyenne)
set.seed(1)
n=1000
U=runif(n)*10
Y=rbinom(n,size=1,prob=MY(U))
plot(U,Y,ylim=c(-.1,1.1))
library(splines)
base=data.frame(Y,U)
reg=glm(Y~bs(U,10),data=base,family=binomial)
u=seq(0,10,by=.01)
y=predict(reg,newdata=data.frame(U=u),type="response")
lines(u,y,lwd=2,col="blue")
lines(u,MY(u),lwd=2,col="red",type="s")

Si on simule quelques milliers d'échantillons, on peut regarder la distribution de mes nœuds (si le découpage par arbre converge), en gardant - au mieux - que deux nœuds,

set.seed(1)
N=matrix(NA,10000,2)
n=1000
for(s in 1:1000){
U=runif(n)*10
Y=rbinom(n,size=1,prob=MY(U))
T=tree(Y~U)
noeud=T$frame$splits[,1]
noeud=noeud[noeud!=""]
NOEUD=as.numeric(substr(noeud,2,nchar(noeud)))
N[s,]=NOEUD[1:2]
}
vn=as.vector(N)
vn=vn[is.na(vn)==FALSE]
plot(density(vn,bw=.2),col="blue",lwd=2,xlim=c(0,10))
abline(v=c(3,7),col="red")

Rappelons que l'on a ici 1000 observations pour faire notre régression, avec trois classes, dont les valeurs moyennes par classe sont respectivement 20%, 50% et 30%. On note que l'on coupe bien entre les deux premières classes, un peu moins bien entre les deux dernières.

Si on suppose que les classes sont moins différentes, par exemple avec des valeurs moyennes à 30%, 50% et 40%, on obtient visuellement la distribution suivante pour les nœuds,

En l’occurrence, on a beaucoup de mal à distinguer entre les deux dernières classes.

Si on change la taille de notre échantillon, en passant de 1000 à 200, on s'aperçoit (en revenant aux valeurs précédentes des taux par classes) que l'on a beaucoup de mal à estimer correctement les nœuds,

De manière assez surprenante, il va prendre des nœuds proches des bords (proches de 0 et de 10).

Moralité, les classes obtenues par arbre ne sont robustes qu'à condition d'avoir suffisamment d'observations, et d'avoir des classes réellement différentes.

ACT2040: devoir maison 2

Pour le second devoir maison, les bases sont disponibles depuis quelques jours sur un billet introduisant les cadences de paiement. Pour reprendre ce que je disais, moyennant une transformation des noms des colonnes afin de pouvoir importer le fichier xls, on peut utiliser le code suivant

> library(gdata)
> read.xls("http://freakonometrics.blog.free.fr/
public/data/act2040triangle.xls", sheet = 10)
essai de l'URL 'http://freakonometrics.blog.free.fr/
public/data/act2040triangle.xls'
Content type 'application/octet-stream' length 454656 bytes (444 Kb)
URL ouverte
==================================================
downloaded 444 Kb
 
Year   X12   X24   X36   X48   X60   X72   X84   X96  X108  X120  X132
1  1995  1340  3188  5072  6973  8677 10008 11802 12606 13174 13596 14033
2  1996  1857  4297  6864  9438 11820 13594 14783 15710 16439 16972    NA
3  1997  2024  4891  7790 10733 13792 16071 17695 18886 19735    NA    NA
4  1998  2781  6655 10671 14738 18022 20795 23179 24597    NA    NA    NA
5  1999  3439  8272 13325 18551 23386 26861 29409    NA    NA    NA    NA
6  2000  3714  9039 14638 20326 26117 30643    NA    NA    NA    NA    NA
7  2001  4652 11236 18109 25239 31250    NA    NA    NA    NA    NA    NA
8  2002  5292 12974 21106 29611    NA    NA    NA    NA    NA    NA    NA
9  2003  6818 16984 27677    NA    NA    NA    NA    NA    NA    NA    NA
10 2004  9337 23263    NA    NA    NA    NA    NA    NA    NA    NA    NA
11 2005 15073    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
Encore une fois, une alternative plus simple est de sauver la feuille excel dans un fichier de type csv (avec séparateur de colonne). On précise alors que le séparateur de colonne est ";". Mais il faut éventuellement préciser aussi que le séparateur de décimale est un ",". Il vaut aussi mieux changer les noms des colonnes dans le fichier source (éviter en particulier les +). Bref, le code ressemble à ça,
> base=read.table("triangle.csv",sep=";",header=TRUE,dec=",")
> base[1:5,1:5]
P       A        B      A.1      B.1
1 1994      NA       NA  9347868 10995499
2 1995      NA  9765812 13234375 15257457
3 1996 2554368 10671566 15275057 17326685
4 1997 4189931 13364857 20313611 23964839
5 1998 3574458 13099576 18393003 23361381
Ensuite, l'énoncé est simple: je veux une prédiction du montant de provisions pour toutes les années, en supposant la première année de survenance comme étant close (i.e. la somme des paiements correspond à la charge ultime), et un intervalle de confiance à 95%. Tout simplement. Pour la première partie, merci de renvoyer le triangle complété, avec le montant de provisions par années. Pour la seconde partie, un graphique avec la distribution du montant de provision sera le bienvenu, mais s'il est expliqué, un intervalle de confiance suffira. Comme annoncé en cour (et contrairement à ce qui était annoncé dans le plan de cours) un rapport de 2 ou 3 pages doit être envoyé avant le 19 décembre, minuit, à charpentier.arthur@uqam.ca.

Friday, December 2 2011

ACT2040: lire une sortie de régression

Histoire que les choses soient bien claires, prenons 2 minutes pour revenir sur la lecture d'une sortie de régression,

>  reg=glm(nbre~carburant+zone+ageconducteur+
+
offset(log(exposition)), + data=baseFREQ,family=poisson(link="log")) > summary(reg)   Call: glm(formula = nbre ~ carburant + zone + ageconducteur + offset(log(exposition)), family = poisson(link = "log"), data = baseFREQ)   Deviance Residuals: Min 1Q Median 3Q Max -0.5372 -0.3643 -0.2731 -0.1503 4.5994   Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.360655 0.105711 -22.331 < 2e-16 *** carburantE -0.283956 0.043117 -6.586 4.53e-11 *** zoneB 0.122891 0.108713 1.130 0.258 zoneC 0.224469 0.088310 2.542 0.011 * zoneD 0.411243 0.086804 4.738 2.16e-06 *** zoneE 0.532070 0.088226 6.031 1.63e-09 *** ageconducteur -0.007031 0.001550 -4.537 5.70e-06 *** --- Signif. codes: 0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1   (Dispersion parameter for poisson family taken to be 1)   Null deviance: 13659 on 49999 degrees of freedom Residual deviance: 13538 on 49993 degrees of freedom AIC: 17788   Number of Fisher Scoring iterations: 6

L'écriture formelle de ce modèle serait quelque chose de la forme

http://freakonometrics.blog.free.fr/public/perso4/sortie11.gif

http://freakonometrics.blog.free.fr/public/perso4/sortie02.gif correspond au nombre de sinistres observés pour l'assuré , http://freakonometrics.blog.free.fr/public/perso4/sortie03.gif correspond à l'exposition (i.e. au temps pendant lequel l'assuré a été observé dans la base), http://freakonometrics.blog.free.fr/public/perso4/sortie04.gif est la première variable, ici le carburant (qui est une variable qualitative prenant 2 modalités), http://freakonometrics.blog.free.fr/public/perso4/sortie05.gif est la seconde variable, ici la zone géographique (qui est encore qualitative, avec 5 zones), et enfin http://freakonometrics.blog.free.fr/public/perso4/sortie07.gif est la troisième et dernière variable explicative, ici l'âge du conducteur principal (qui est pris ici comme une variable continue). Le modèle s'écrit, si l'on disjoncte les variables qualitatives http://freakonometrics.blog.free.fr/public/perso4/sortie08.gifhttp://freakonometrics.blog.free.fr/public/perso4/sortie12.gif s'écrit

http://freakonometrics.blog.free.fr/public/perso4/sortie10.gif

Il y alors 6 paramètres à estimer. C'est ce qui est renvoyé dans la sortie de la régression,

               Estimate Std. Error z value Pr(>|z|)
(Intercept)   -2.360655   0.105711 -22.331  < 2e-16 ***
carburantE    -0.283956   0.043117  -6.586 4.53e-11 ***
zoneB          0.122891   0.108713   1.130    0.258
zoneC          0.224469   0.088310   2.542    0.011 *
zoneD          0.411243   0.086804   4.738 2.16e-06 ***
zoneE          0.532070   0.088226   6.031 1.63e-09 ***
ageconducteur -0.007031   0.001550  -4.537 5.70e-06 ***
Dans ce tableau, on a pour chaque paramètre, une estimation, par exemple http://freakonometrics.blog.free.fr/public/perso4/sortie20.gif vaut -2.36, alors que http://freakonometrics.blog.free.fr/public/perso4/sortie21.gif vaut 0.2244. On a ensuite une estimation de l'écart-type de ces estimateurs, par exemple http://freakonometrics.blog.free.fr/public/perso4/sortie23.gif vaut 0.086.  On peut alors utiliser un test de significative basé sur une hypothèse de normalité de ces estimateurs. On note que l'âge du conducteur est significatif dans ce modèle. Pour les facteurs, la significativité est jugée par rapport à la modalité de référence, qui est celle qui n'apparaît pas dans la sortie de la régression. Aussi, le carburant "Essence" est significatif par rapport à la modalité de référence qui est le Diesel. En revanche, la zone "B" n'est pas significative par rapport à la zone de référence qui est la zone "A", ce qui suggèrerait de regrouper les zone "A" et "B" (comme discuté dans un ancien billet). Notons que l'on peut d'ailleurs aller un peu plus loin, en utilisant des estimateurs plus robustes de l'écart-type de ces estimateurs (proposés par Cameron & Trivedi (1998)),
> library(sandwich)
> cov.reg<-vcovHC (reg, type="HC0")
> std.err<-sqrt(diag(cov.reg))
> estimation<-cbind(estimate=reg$coefficients, std.err,
+        pvalue=round(2*(1-pnorm(abs(reg$coefficients/std.err))),5),
+        lower=reg$coefficients-1.96*std.err,
+        upper=reg$coefficients+1.96*std.err)
>
> estimation
estimate     std.err  pvalue       lower        upper
(Intercept)   -2.360654805 0.110541249 0.00000 -2.57731565 -2.143993957
carburantE    -0.283955585 0.044492586 0.00000 -0.37116105 -0.196750117
zoneB          0.122890716 0.111122177 0.26877 -0.09490875  0.340690183
zoneC          0.224468748 0.091032422 0.01367  0.04604520  0.402892295
zoneD          0.411242820 0.089449063 0.00000  0.23592266  0.586562984
zoneE          0.532070373 0.091675282 0.00000  0.35238682  0.711753926
ageconducteur -0.007030691 0.001664597 0.00002 -0.01029330 -0.003768081
On retrouve ici des grandeurs du même ordre que tout à l'heure, avec un intervalle de confiance à 95% pour les estimateurs. Si 0 appartient à l'intervalle de confiance, on dira que le paramètre n'est pas significatif. A partir de ces estimations, il est facile de faire une prédiction, par exemple pour un assuré de 40 ans, résidant dans la zone D et conduisant un véhicule diesel, sa espérance annuelle d'accident est de 10.7%
>  predict(reg,newdata=data.frame(carburant="D",
+  zone="D",ageconducteur=40,exposition=1),type="response")
1
0.1074597
>  cbind(c(1,0,0,0,1,0,40),reg$coefficients)
[,1]         [,2]
(Intercept)      1 -2.360654805
carburantE       0 -0.283955585
zoneB            0  0.122890716
zoneC            0  0.224468748
zoneD            1  0.411242820
zoneE            0  0.532070373
ageconducteur   40 -0.007030691
>  exp(t(c(1,0,0,0,1,0,40))%*%reg$coefficients)
[,1]
[1,] 0.1074597
Notons que dans ce cas, la probabilité de ne pas avoir d'accident est
>  lambda=exp(t(c(1,0,0,0,1,0,40))%*%reg$coefficients)
>  dpois(0,lambda)
[1] 0.8981127
>  exp(-lambda)
[,1]
[1,] 0.8981127
Parmi les autres interprétations possibles de ces valeurs, notons qu'à zone et à âge identique, un assuré conduisant un véhicule essence à une fréquence de sinistres 25% plus faible qu'un assuré conduisant un véhicule diesel,
>  exp(reg$coefficients[2])
carburantE
0.7528001
En bas, on a la déviance du modèle. Dans le cas d'une régression de Poisson, on sait que
http://freakonometrics.blog.free.fr/public/perso4/sortie25.gif
soit ici
>  baseFREQ$pred=predict(reg,type="response")
>  Z=log(baseFREQ$pred^baseFREQ$nbre)
>  2*(sum(log(baseFREQ$nbre^baseFREQ$nbre)-baseFREQ$nbre)-
+  sum(Z- baseFREQ$pred))
[1] 13538.09
> deviance(reg)
[1] 13538.09
Mais la déviance n'est pas forcément pertinente pour comparer des modèles ayant des nombres différentes de paramètres, et donc la sortie inclue le critère d'Akaike.

Thursday, December 1 2011

ACT2040: zonage versus densité de population

Hier, j'avais évoqué le fait que les régions géographiques n'étaient pas de variables continues. Dit comme ça, c'est partiellement faux: il existe une certaine forme de continuité spatiale. Mais pas basée sur des codes arbitraires de régions (tels qu'ils sont utilisés dans la base: la Basse Normandie (25) est plus proche de la Bretagne (53) que de la Bourgogne (26)). En particulier la continuité spatiale est reflétée par la densité de population. Les conducteurs dans les grosses métropoles ont des conduites proches, pareil pour les personnes habitant dans des régions peu peuplées. Si l'on regarde la fréquence de sinistres en fonction de la densité on observe qu'il y a une relative continuité, avec beaucoup plus d'accidents dans les zones à forte concentration,

C'est ce qu'on observe quand on regarde les fréquences par code postal, en Belgique (je n'ai pas le code postal dans ma base donc je ne peux pas le faire sur les données à notre disposition), avec des fortes concentrations dans les régions très urbaines, ce qui donne une forme de continuité en lissant spatialement

Moralité, non, on ne peut pas utiliser les codes des départements, par contre, les regroupements basés sur la densité de population sont pertinents.

ACT2040: poids, variable offset et régression binomiale négative

  • Variable offset et processus de Poisson
Comme on l'a vu en cours, le modèle fondamental pour modéliser l'arrivée des sinistres est le processus de Poisson (homogène). Rappelons que si http://freakonometrics.blog.free.fr/public/perso4/ooofset-01.gif est un processus de Poisson (désignant le nombre de sinistres survenant entre http://freakonometrics.blog.free.fr/public/perso4/ooofset-03.gif et http://freakonometrics.blog.free.fr/public/perso4/ooofset-02.gif), alors
http://freakonometrics.blog.free.fr/public/perso4/ooofset-04.gif
i.e. http://freakonometrics.blog.free.fr/public/perso4/ooofset-05.gif suit une loi de Poisson de paramètre http://freakonometrics.blog.free.fr/public/perso4/ooofset-06.gif.
Mais en plus les incréments sont indépendants, au sens où http://freakonometrics.blog.free.fr/public/perso4/ooofset-05.gif et http://freakonometrics.blog.free.fr/public/perso4/ooofset-07.gif sont deux variables de Poisson indépendantes.
Une conséquence est que si le nombre de sinistres pour un assuré présent une année suit une loi de Poisson http://freakonometrics.blog.free.fr/public/perso4/ooofset-08.gif, pour un assuré présent pendant une durée http://freakonometrics.blog.free.fr/public/perso4/ooofset-10.gif, le nombre de sinistre suit une loi de Poisson http://freakonometrics.blog.free.fr/public/perso4/ooofset-09.gif.
Si on revient à la régression de Poisson, on suppose que
http://freakonometrics.blog.free.fr/public/perso4/ooofset-11.gif
mais si on rajoute l'exposition http://freakonometrics.blog.free.fr/public/perso4/ooofset-10.gif, alors
http://freakonometrics.blog.free.fr/public/perso4/ooofset-12.gif
Autrement dit la variable http://freakonometrics.blog.free.fr/public/perso4/ooofset-10.gif intervient dans la régression (ou plutôt son logarithme), mais c'est une variable particulière car le coefficient est ici 1 (et n'est pas à estimer). L'exposition est alors une variable en dehors de l'ensemble des autres variables (offset). Pour intégrer une variable offset, le code est tout simplement
> reg1=glm(nombre~ageconducteur+offset(log(exposition)),
+ data=baseFREQ,family=poisson(link="log"))
  • Variable offset et régression binomiale négative
Dans le cas de la régression binomiale négative, de loi
http://freakonometrics.blog.free.fr/public/perso4/ooofset-13.gif
avec http://freakonometrics.blog.free.fr/public/perso4/ooofset-14.gif et http://freakonometrics.blog.free.fr/public/perso4/ooofset-15.gif.
Le modèle de régression construit à partir de cette loi suppose
http://freakonometrics.blog.free.fr/public/perso4/ooofset-16.gif
Compte tenu du parallèle qui existe entre la loi de Poisson et la loi binomiale négative (on rajoute une variable non-observable suivant une loi Gamma, cf. exercice 3 du partiel), on peut rajouter une variable offset en écrivant
http://freakonometrics.blog.free.fr/public/perso4/ooofset-17.gif
Le code R pour faire une régression binomiale négative est tout simplement
> library(MASS)
> reg2=glm.nb(nombre~ageconducteur+offset(log(exposition)),
+ data=baseFREQ)
On peut d'ailleurs comparer les deux prédictions, sur le portefeuille
>  age=seq(18,85)
>  nb1=predict(reg1,newdata=data.frame(ageconducteur=
+  age,exposition=1),type="response")
>  nb2=predict(reg2,newdata=data.frame(ageconducteur=
+  age,exposition=1), type="response")
> X=seq(10,110,by=5)
> V=tapply(baseFREQ$exposition,cut(baseFREQ$ageconducteur,X),sum)
> plot((X[2:length(X)]+X[1:(length(X)-1)])/2,V,axes=
+ FALSE,col="white",xlab="",ylab="")
> axis(1)
> axis(2)
> for(i in 1:(length(X)-1)){
+ polygon(c(X[i],X[i],X[i+1],X[i+1]),
+         c(0,V[i],V[i],0),col="light green",border=NA)
+ }
> par(new=TRUE)
>  plot(age,nb1,type="l",col="red",axes=FALSE,ylab="",xlab="")
> lines(age,nb2,col="blue")
>  reg3=glm(nombre~as.factor(ageconducteur)+
+ offset(log(exposition)), data=baseFREQ,family=
+ poisson(link="log"))
>  nb3=predict(reg3,newdata=data.frame(ageconducteur=
+ age,exposition=1), type="response")
>  points(age,nb3,pch=3)
> axis(4)
avec en vert l'exposition en fonction de l'âge, en rouge l'espérance du nombre de sinistres par an avec un modèle de Poisson (l'échelle est à droite sur le graphique), en bleu le modèle binomial négatif, et les + correspondent aux estimateurs non-paramétriques de la moyenne du nombre annuel d'accidents par âge. On notera que les sorties de régression sont très proches (les deux modèles étant comparables car on a la même fonction lien),
> summary(reg1)
 
Call:
glm(formula = nombre ~ ageconducteur + offset(log(exposition)),
family = poisson(link = "log"), data = baseFREQ)
 
Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.5068  -0.3759  -0.2776  -0.1521   4.6391
 
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)   -2.135616   0.072877 -29.304  < 2e-16 ***
ageconducteur -0.007913   0.001533  -5.163 2.43e-07 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
(Dispersion parameter for poisson family taken to be 1)
 
Null deviance: 13665  on 49999  degrees of freedom
Residual deviance: 13638  on 49998  degrees of freedom
AIC: 17882
 
Number of Fisher Scoring iterations: 6
 
> summary(reg2)
 
Call:
glm.nb(formula = nombre ~ ageconducteur + offset(log(exposition)),
data = baseFREQ, init.theta = 0.8502027145, link = log)
 
Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.4904  -0.3699  -0.2757  -0.1522   4.0075
 
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)   -2.118145   0.075283 -28.136  < 2e-16 ***
ageconducteur -0.008138   0.001583  -5.141 2.74e-07 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
(Dispersion parameter for Negative Binomial(0.8502) family taken to be 1)
 
Null deviance: 11727  on 49999  degrees of freedom
Residual deviance: 11700  on 49998  degrees of freedom
AIC: 17828
 
Number of Fisher Scoring iterations: 1
 
 
Theta:  0.850
Std. Err.:  0.148
 
2 x log-likelihood:  -17821.593 
On notera que l'utilisation d'un modèle binomial négatif est légitimé par la valeur du paramètre de dispersion dans le modèle quasi-Poisson
> reg1b=glm(nombre~ageconducteur+offset(log(exposition)),
+ data=baseFREQ,family=quasipoisson(link="log"))
> summary(reg1b)
 
Call:
glm(formula = nombre ~ ageconducteur + offset(log(exposition)),
family = quasipoisson(link = "log"), data = baseFREQ)
 
Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.5068  -0.3759  -0.2776  -0.1521   4.6391
 
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)   -2.135616   0.090143 -23.691   <2e-16 ***
ageconducteur -0.007913   0.001896  -4.174    3e-05 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
(Dispersion parameter for quasipoisson family taken to be 1.52996)
 
Null deviance: 13665  on 49999  degrees of freedom
Residual deviance: 13638  on 49998  degrees of freedom
AIC: NA
 
Number of Fisher Scoring iterations: 6

  • Poids versus variable offset dans une régression de Poisson
Rappelons que la loi (conditionnelle) dans une régression log-Poisson est
http://freakonometrics.blog.free.fr/public/perso4/ooofset-20.gif
de telle sorte que la log-vraisemblance s'écrit
http://freakonometrics.blog.free.fr/public/perso4/ooofset-21.gif
La fonction que l'on cherche donc à maximiser est alors
Si l'on rajoute des poids, on va chercher à minimiser la fonction
\
Dans le cas où l'on rajoute une variable offset, on suppose que
http://freakonometrics.blog.free.fr/public/perso4/ooofset-25.gif
de telle sorte que
http://freakonometrics.blog.free.fr/public/perso4/ooofset-24.gif
et donc la log-vraisemblance devient
http://freakonometrics.blog.free.fr/public/perso4/ooofset-26.gif
et la fonction que l'on va maximiser est alors
Les deux fonctions sont manifestement différentes, le poids n'intervenant que sur le second terme dans le cas de la variable offset. La régression de Poisson pondérée donne ici
> reg4=glm(nombre~ageconducteur,weight=exposition,
+  data=baseFREQ,family=poisson(link="log"))
Les ordres de grandeur des prédictions ne sont pas valides ici... Néanmoins, il est possible d'utiliser des poids, comme on le voit sur le calcul rapide ci-dessous,
> M=rep(NA,85-18+1)
> for(i in 18:85){
+ M[i-17]=weighted.mean(x=(baseFREQ$nombre/baseFREQ$exposition)
+ [baseFREQ$ageconducteur==i],
+ w=baseFREQ$exposition[baseFREQ$ageconducteur==i])
+ }
> M
[1] 0.25314943 0.20047905 0.20409416 0.23199038 0.17182485
[6] 0.17236472 0.10258972 0.10115476 0.10784129 0.10133346
[11] 0.08609449 0.05025678 0.08433958 0.06400527 0.06116849
[16] 0.07542661 0.09284957 0.06258664 0.08251126 0.07963923
[21] 0.08038180 0.08674745 0.06699507 0.09803112 0.09619517
[26] 0.08866225 0.09627715 0.09710725 0.08207893 0.09162368
[31] 0.07207799 0.08073084 0.09107094 0.09631251 0.06723020
[36] 0.06862383 0.08284661 0.08548334 0.05930442 0.08185568
[41] 0.07097164 0.07058023 0.06625962 0.08141103 0.10569807
[46] 0.04538514 0.05289389 0.09161058 0.04131018 0.05323120
[51] 0.08052976 0.06241604 0.06418933 0.04717143 0.11227675
[56] 0.05325793 0.04484281 0.04115995 0.11487048 0.06827367
[61] 0.07847543 0.08735320 0.07047837 0.05950256 0.10550113
[66] 0.04496403 0.07965526 0.06195787
> nb3
1          2          3          4          5
0.25314943 0.20047905 0.20409416 0.23199038 0.17182485
6          7          8          9         10
0.17236472 0.10258972 0.10115476 0.10784129 0.10133346
11         12         13         14         15
0.08609449 0.05025678 0.08433958 0.06400527 0.06116849
16         17         18         19         20
0.07542661 0.09284957 0.06258664 0.08251126 0.07963923
21         22         23         24         25
0.08038180 0.08674745 0.06699507 0.09803112 0.09619517
26         27         28         29         30
0.08866225 0.09627715 0.09710725 0.08207893 0.09162368
31         32         33         34         35
0.07207799 0.08073084 0.09107094 0.09631251 0.06723020
36         37         38         39         40
0.06862383 0.08284661 0.08548334 0.05930442 0.08185568
41         42         43         44         45
0.07097164 0.07058023 0.06625962 0.08141103 0.10569807
46         47         48         49         50
0.04538514 0.05289389 0.09161058 0.04131018 0.05323120
51         52         53         54         55
0.08052976 0.06241604 0.06418933 0.04717143 0.11227675
56         57         58         59         60
0.05325793 0.04484281 0.04115995 0.11487048 0.06827367
61         62         63         64         65
0.07847543 0.08735320 0.07047837 0.05950256 0.10550113
66         67         68
0.04496403 0.07965526 0.06195787 
L'idée est d'annualiser les nombres de sinistres, en considérant , mais de mettre un poids plus faible si l'assuré n'a pas été observé longtemps, i.e. un poids . Si on regarde maintenant dans la régression, on obtient exactement l'estimation obtenue avec la variable offset,
i.e
Bref, la régression avec une variable offset est une régression pondérée, mais sur le nombre de sinistres annualisé. Numériquement, on a
> reg5=glm((nombre/exposition)~ageconducteur,weight=exposition,
+ data=baseFREQ,family=poisson(link="log"))
> nb5=predict(reg5,newdata=data.frame(ageconducteur=
+ age,exposition=1), type="response")
> nb1
1          2          3          4          5
0.10248395 0.10167620 0.10087482 0.10007975 0.09929095
6          7          8          9         10
0.09850836 0.09773194 0.09696165 0.09619742 0.09543922
11         12         13         14         15
0.09468699 0.09394070 0.09320028 0.09246570 0.09173691
16         17         18         19         20
0.09101387 0.09029652 0.08958483 0.08887874 0.08817823
21         22         23         24         25
0.08748323 0.08679371 0.08610962 0.08543093 0.08475759
26         27         28         29         30
0.08408955 0.08342678 0.08276923 0.08211687 0.08146965
31         32         33         34         35
0.08082752 0.08019046 0.07955842 0.07893136 0.07830925
36         37         38         39         40
0.07769204 0.07707969 0.07647217 0.07586943 0.07527145
41         42         43         44         45
0.07467818 0.07408959 0.07350564 0.07292628 0.07235150
46         47         48         49         50
0.07178124 0.07121548 0.07065418 0.07009730 0.06954482
51         52         53         54         55
0.06899668 0.06845287 0.06791334 0.06737807 0.06684701
56         57         58         59         60
0.06632014 0.06579742 0.06527883 0.06476432 0.06425386
61         62         63         64         65
0.06374743 0.06324499 0.06274651 0.06225196 0.06176131
66         67         68
0.06127452 0.06079157 0.06031243
> nb5
1          2          3          4          5
0.10248395 0.10167620 0.10087482 0.10007975 0.09929095
6          7          8          9         10
0.09850836 0.09773194 0.09696165 0.09619742 0.09543922
11         12         13         14         15
0.09468699 0.09394070 0.09320028 0.09246570 0.09173691
16         17         18         19         20
0.09101387 0.09029652 0.08958483 0.08887874 0.08817823
21         22         23         24         25
0.08748323 0.08679371 0.08610962 0.08543093 0.08475759
26         27         28         29         30
0.08408955 0.08342678 0.08276923 0.08211687 0.08146965
31         32         33         34         35
0.08082752 0.08019046 0.07955842 0.07893136 0.07830925
36         37         38         39         40
0.07769204 0.07707969 0.07647217 0.07586943 0.07527145
41         42         43         44         45
0.07467818 0.07408959 0.07350564 0.07292628 0.07235150
46         47         48         49         50
0.07178124 0.07121548 0.07065418 0.07009730 0.06954482
51         52         53         54         55
0.06899668 0.06845287 0.06791334 0.06737807 0.06684701
56         57         58         59         60
0.06632014 0.06579742 0.06527883 0.06476432 0.06425386
61         62         63         64         65
0.06374743 0.06324499 0.06274651 0.06225196 0.06176131
66         67         68
0.06127452 0.06079157 0.06031243 
La seule difficulté est qu'il faut admettre que l'on puisse avoir des valeurs non-entière dans une régression de Poisson (ce qui est une toute autre histoire....).

Wednesday, November 30 2011

ACT2040: une région géographique n'est pas une variable continue

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 codification officielle). Par exemple, dans une des bases, nous avions des assurés dans 4 zones géographiques, à savoir la région 82 (région Rhône-Alpes en rouge) la région 54 (région Poitou-Charentes en vert) la région 73 (région Midi-Pyrénées en bleu) et enfin la région 41 (région Lorraine en mauve).

> unique(baseFREQ$region)
[1] 82 54 73 41

Une idée intéressante pour regrouper les régions pouvait être d'utiliser les arbres. Les régions étant des couleurs (on le voit bien sur la carte) et pas des variables quantitatives, il est normal de travailler sur des facteurs. D'ailleurs le code pour faire la carte est le suivant,

> library(maps)
>  france<-map(database="france")
>  dpt=c("Ain","Ardeche","Drome","Isere","Loire ","Rhone",
+ "Savoie","Haute-Savoie","Charente","Charente-Maritime",
+ "Deux-Sevres","Vienne","Ariege","Aveyron","Haute-Garonne",
+ "Gers","Lot","Hautes-Pyrenees","Tarn","Tarn-et-Garonne",
+ "Meurthe-et-Moselle","Meuse","Moselle","Vosges")
>  couleur=c(rep(2,8),rep(3,4),rep(4,8),rep(6,4))
>  match=match.map(france,dpt)
>  color=couleur[match]
>  map(database="france", fill=TRUE, col=color)

L'arbre sur les régions en tant que facteurs donne le découpage suivant

>  baseFREQ$fregion=as.factor(baseFREQ$region)
>  ARBRE1=tree(nombre~fregion,data=baseFREQ,split="gini")
>  plot(ARBRE1)
>  text(ARBRE1)

Bon, R a la mauvaise idée de recoder les classes (mais il garde l'ordre, i.e. a correspond à la région 41, b à 54, c à 73 et d à 82). Visuellement, on retient qu'il est possible de considérer deux grandes régions, ac (i.e. 41 et 73) et bd (i.e. 54 et 72). L'intérêt des arbres sur des variables qualitatives, des facteurs, c'est que tous les regroupements sont possibles. En revanche, si on fait un arbre sur la région qui est lue en tant que nombre (quantitatif), on obtient

>  ARBRE2=tree(nombre~region,data=baseFREQ,split="gini")
>  plot(ARBRE2)
>  text(ARBRE2)

Il est alors impossible de regrouper dans une même classe deux régions séparées par un nombre, i.e. on ne peut regrouper 41 et 82 dans la même classe. R suggère de distinguer peut être trois régions, à savoir 82 (à droite), puis 73 (au centre) et enfin de mettre éventuellement 41 et 54 ensemble. Ce qui n'est pas la stratégie optimale quand on regroupe des facteurs.

Tuesday, November 29 2011

ACT2040: chain ladder

Pour la fin du cours, on a commencé à voir le méthode dit Chain Ladder. Pour cela, rappelons que l'on dispose du jeu de données suivants, de paiements cumulés,

> source("http://perso.univ-rennes1.fr/
arthur.charpentier/bases.R")
> PAID
     [,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
> PAID/PREMIUM
       [,1]   [,2]   [,3]   [,4]   [,5]   [,6]
[1,] 0.6989 0.9522 0.9607 0.9644 0.9660 0.9705
[2,] 0.7206 0.9972 1.0051 1.0102 1.0124     NA
[3,] 0.7960 1.0991 1.1100 1.1145     NA     NA
[4,] 0.8191 1.1433 1.1632     NA     NA     NA
[5,] 0.8688 1.1976     NA     NA     NA     NA
[6,] 0.8112     NA     NA     NA     NA     NA

si on visualise l'évolution du ratio sinistre/primes. On peut visualiser les facteurs de transition dans le tableau ci-dessous,

> PAID[,2:6]/PAID[,1:5]
[,1]     [,2]     [,3]     [,4]     [,5]
[1,] 1.362418 1.008920 1.003854 1.001581 1.004735
[2,] 1.383724 1.007942 1.005111 1.002119       NA
[3,] 1.380780 1.009916 1.004076       NA       NA
[4,] 1.395848 1.017407       NA       NA       NA
[5,] 1.378373       NA       NA       NA       NA
[6,]       NA       NA       NA       NA       NA

Toutefois, pour faire le ratio moyen, on ne fait pas la moyenne des ratios, ou alors à condition de pondérer correctement, par exemple pour passer de la première à la seconde ligne du tableau,

> k=1
> weighted.mean(x=PAID[,k+1]/PAID[,k],w=PAID[,k],na.rm=TRUE)
[1] 1.380933
> sum(PAID[1:(6-k),k+1])/sum(PAID[1:(6-k),k])
[1] 1.380933
En itérant, on récupère l'ensemble des coefficients de transition,
> 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.011433 1.004343 1.001858 1.004735
Une fois récupérés tous les coefficients de transition, on peut compléter la partie inférieure du triangle (ou de la matrice),
> PROJECTION=PAID
> for(k in 1:5){
+ PROJECTION[((7-k):6),k+1]=PROJECTION[((7-k):6),k]*lambda[k]
+ print(PROJECTION)}
[,1]     [,2] [,3] [,4] [,5] [,6]
[1,] 3209 4372.000 4411 4428 4435 4456
[2,] 3367 4659.000 4696 4720 4730   NA
[3,] 3871 5345.000 5398 5420   NA   NA
[4,] 4239 5917.000 6020   NA   NA   NA
[5,] 4929 6794.000   NA   NA   NA   NA
[6,] 5217 7204.327   NA   NA   NA   NA
[,1]     [,2]     [,3] [,4] [,5] [,6]
[1,] 3209 4372.000 4411.000 4428 4435 4456
[2,] 3367 4659.000 4696.000 4720 4730   NA
[3,] 3871 5345.000 5398.000 5420   NA   NA
[4,] 4239 5917.000 6020.000   NA   NA   NA
[5,] 4929 6794.000 6871.672   NA   NA   NA
[6,] 5217 7204.327 7286.691   NA   NA   NA
[,1]     [,2]     [,3]     [,4] [,5] [,6]
[1,] 3209 4372.000 4411.000 4428.000 4435 4456
[2,] 3367 4659.000 4696.000 4720.000 4730   NA
[3,] 3871 5345.000 5398.000 5420.000   NA   NA
[4,] 4239 5917.000 6020.000 6046.147   NA   NA
[5,] 4929 6794.000 6871.672 6901.518   NA   NA
[6,] 5217 7204.327 7286.691 7318.339   NA   NA
[,1]     [,2]     [,3]     [,4]     [,5] [,6]
[1,] 3209 4372.000 4411.000 4428.000 4435.000 4456
[2,] 3367 4659.000 4696.000 4720.000 4730.000   NA
[3,] 3871 5345.000 5398.000 5420.000 5430.072   NA
[4,] 4239 5917.000 6020.000 6046.147 6057.383   NA
[5,] 4929 6794.000 6871.672 6901.518 6914.344   NA
[6,] 5217 7204.327 7286.691 7318.339 7331.939   NA
[,1]     [,2]     [,3]     [,4]     [,5]     [,6]
[1,] 3209 4372.000 4411.000 4428.000 4435.000 4456.000
[2,] 3367 4659.000 4696.000 4720.000 4730.000 4752.397
[3,] 3871 5345.000 5398.000 5420.000 5430.072 5455.784
[4,] 4239 5917.000 6020.000 6046.147 6057.383 6086.065
[5,] 4929 6794.000 6871.672 6901.518 6914.344 6947.084
[6,] 5217 7204.327 7286.691 7318.339 7331.939 7366.656
où on visualise, étape par étape, le remplissage de la partie inférieure de la matrice. La dernière colonne contient des prédictions des charges ultimes, par année de survenance, et sur la diagonale, on a toujours le montant de paiements déjà effectués par année de survenance. Le montant de provision est alors la différence entre ce que l'on pense payer, et ce que l'on a déjà payé.
> PROJECTION[,6]-diag(PAID[,6:1])
[1]    0.00   22.3968   35.7838   66.0646  153.0835 2149.6564
> sum(PROJECTION[,6]-diag(PAID[,6:1]))
[1] 2426.985



Monday, November 28 2011

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=";")
> 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"  

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

- page 1 of 3