Tail factor et clôture 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

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

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.gifOn 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("https://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.

Incréments négatifs dans les triangles de paiements

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 paiementné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 courbenoire 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…