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,
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
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
et dans l’autre
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
(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
(avec la filtration correspondant à la partie supérieure du triangle disponible), i.e.
Comme on utilise un modèle de Poisson, on suppose aussi que
i.e.
c’est à dire que
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),
sauf que l’estimateur de la variance est biaisé. Certes, classiquement, l’estimateur est basé sur une normalisation par un facteur 
(classique en statistique quand on possède 
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 
. 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 
, 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 
ou bien simuler une loi binomiale négative, de moyenne 
, et de paramètre de surdispersion 
telle que la variance s’écrive
Aussi, on prend 
alors que
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…