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,

plot(density(Rs))

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  https://perso.univ-rennes1.fr/arthur.charpentier/latex/qp001.png, mais aussi quehttps://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

https://perso.univ-rennes1.fr/arthur.charpentier/latex/qp003.png et https://perso.univ-rennes1.fr/arthur.charpentier/latex/qp004.png

ou bien simuler une loi binomiale négative, de moyenne https://perso.univ-rennes1.fr/arthur.charpentier/latex/qp008.png, et de paramètre de surdispersion https://perso.univ-rennes1.fr/arthur.charpentier/latex/qp009.png telle que la variance s’écrive

https://perso.univ-rennes1.fr/arthur.charpentier/latex/qp010.png

Aussi, on prend https://perso.univ-rennes1.fr/arthur.charpentier/latex/qp011.png alors que

https://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…


2 thoughts on “Quels résidus et quelle loi simuler ?”

  1. Normalement elle ne disparait pas… c’est pour ça que je dis qu’utiliser un résidu ou l’autre revient exactement au même lorsque l’on génère des pseudo-triangles. Sauf que plusieurs personnes (je pense à l’Appendix 3 de l’article d’England ou à son slide 18 qui évoque ce point sans le justifier) suggèrent d’utiliser la formule des résidus non-ajustés (sans la racine) pour générer des pseudo-triangles, mais au lieu de mettre les résidus non-ajustés, mettre au contraire les résidus ajustés (les r). Le “peut sembler légitime” signifie que je n’ai pas encore vraiment compris ce qui autorise ou valide ce bricolage (hormis la petite histoire que je raconte à la fin du billet) mais que beaucoup de monde trouve ça normal….

  2. “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”

    Je n’ai pas compris comment la racine [(n-k)/n] disparaît?

    Merci

Leave a Reply

Your email address will not be published. Required fields are marked *