Parmi les problèmes qui ont connu leur heure de gloire il y a quelques
années, il y a eu le sujet des triangles multivariés, ou plutôt de la
dépendance qui pourrait exister entre des triangles de paiements.
Considérons les deux triangles suivants, obtenus auprès d'un même
assureur, pour des sinistres en auto matériel entre 1994 et 2003, i.e.

et pareil pour les sinistres en auto-corporel,

On
peut alors envisager deux techniques: travailler sur le triangle
agrégé, ou travailler sur les triangles séparer, puis les agréger les
prédictions ensuites,
- Comparaison des prédictions par la méthode de Mack (i.e. chain-ladder)
On
peut commencer par utiliser la formule de Mack, en supposant que le
montant des provisions suit une loi lognormale. Rappelons que si

alors les deux premiers moments vérifient
et
ce qui s'inverse simplement, et permet d'écrire
et
Aussi, à partir du montant prédit et de l'écart-type, on peut en déduire les paramètres de la loi sous-jacente.
> P.corp=read.table("http://perso.univ-rennes1.fr/arthur.charpentier/auto-corporel.csv",
+ header=FALSE,sep=";",na.strings = "NA",dec=",")
> P.corp=as.matrix(P.corp)
> n=nrow(P.corp)
> P.mat =read.table("http://perso.univ-rennes1.fr/arthur.charpentier/auto-materiel.csv",
+ header=FALSE,sep=";",na.strings = "NA",dec=",")
> P.mat=as.matrix(P.mat)
> P.mat=P.mat[1:n,1:n]
> P.mat = P.mat[2:10,1:9]
> P.corp= P.corp[2:10,1:9]
> n=9
> P.tot = P.mat + P.corp
> library(ChainLadder)
> V=MackChainLadder(P.tot)$Total.Mack.S.E
> E=sum(MackChainLadder(P.tot)$FullTriangle[,n]-
+ diag(MackChainLadder(P.tot)$FullTriangle[n:1,]))
> mu = log(E) - .5*log(1+V^2/E^2)
> sigma2 = log(1+V^2/E^2)
Les trois distributions sont alors les suivantes (sous hypothèse de lognormalité)

Si on suppose que les deux branches sont
indépendantes, on peut alors regarder comment se comporte la convolée des deux lois log-normales,
> library(distr)
> V=MackChainLadder(P.mat)$Total.Mack.S.E
> E=sum(MackChainLadder(P.mat)$FullTriangle[,n]-
+-diag(MackChainLadder(P.mat)$FullTriangle[n:1,]))
> mu = log(E) - .5*log(1+V^2/E^2)
> sigma2 = log(1+V^2/E^2)
> LM = Lnorm(meanlog=mu,sdlog=sqrt(sigma2))
> V=MackChainLadder(P.corp)$Total.Mack.S.E
> E=sum(MackChainLadder(P.corp)$FullTriangle[,n]-
+ diag(MackChainLadder(P.corp)$FullTriangle[n:1,]))
> mu = log(E) - .5*log(1+V^2/E^2)
> sigma2 = log(1+V^2/E^2)
> LC = Lnorm(meanlog=mu,sdlog=sqrt(sigma2))
> LT=LM+LCOn peut alors tracer les densités
> u=seq(0,qlnorm(.95,mu,sqrt(sigma2)),length=1000)
> vtotal=dlnorm(u,mu,sqrt(sigma2))
> vconvol=d(LT)(u)
Si
on compare les quantiles (on se contentera du quantile à 90% car on a
des soucis numériques sur le calcul du quantile de la convolée pour un
niveau trop élevé),
> q(LT)(.9)
[1] 271463.7
> qlnorm(.9,mu,sqrt(sigma2))
[1] 365395.1
Autrement dit, on a le choix sur l'interprétation,
- travailler sur les triangles agrégés apporte beaucoup trop d'incertitudes,
- supposer les triangles matériels et corporels indépendants est une hypothèse trop forte.
Un
approche un peu extrême consisterait à supposer les risques
comonotones, autrement dit le quantile de la somme serait la somme des
quantiles,
> qlnorm(.9,mu,sqrt(sigma2))
[1] 64262.87
> qlnorm(.9,mu,sqrt(sigma2))
[1] 230645
dont la somme fait
[1] 294907.9qu'il faudrait comparer à
> qlnorm(.9,mu,sqrt(sigma2))
[1] 365395.1Bref, supposer les risques comontones ne suffit manifestement pas (je renvoie
ici pour un billet expliquant que le quantile d'une somme est un objet assez complexe).
Une
autre idée est de supposer que les résidus obtenus suite à un GLM sont
- éventuellement - corrélés. Commençons par faire de simples
régressions de Poisson (on va mettre les incréments négatifs à 0, comme
on le verra ci-dessous, ça ne cahnge pas grand chose sur le montant
total de provision.... ceux qui le souhaitent rafineront),
> an <- n; ligne = rep(1:an, each=an); colonne = rep(1:an, an)
> passe = (ligne + colonne - 1)<=an; n = sum(passe)
> PAID=P.corp; INC=PAID
> INC[,2:an]=PAID[,2:an]-PAID[,1:(an-1)]
> I.corp = INC
> PAID=P.mat; INC=PAID
> INC[,2:an]=PAID[,2:an]-PAID[,1:(an-1)]
> I.mat = INC
> Ym = as.vector(I.mat)
> Ym[Ym<0]=0
> Yc = as.vector(I.corp)
> lig = as.factor(ligne)
> col = as.factor(colonne)
> base = data.frame(Ym,Yc,col,lig)
> regm=glm(Ym~col+lig,data=base,family="poisson")
Il y a eu 35 avis (utilisez warnings() pour les visionner)
> regc=glm(Yc~col+lig,data=base,family="poisson")
Il y a eu 38 avis (utilisez warnings() pour les visionner)
> MackChainLadder(P.mat)
MackChainLadder(Triangle = P.mat)
Latest Dev.To.Date Ultimate IBNR Mack.S.E CV(IBNR)
1 2,897,694 1.000 2,897,694 0.00 0.00 NaN
2 353,866 1.000 353,873 6.96 4.52 0.649
3 358,922 1.000 358,926 3.97 21.89 5.515
4 342,988 1.000 343,020 32.27 47.55 1.473
5 288,867 1.000 288,920 52.64 52.77 1.002
6 297,708 1.000 297,775 66.67 371.73 5.576
7 270,537 0.999 270,891 354.26 820.88 2.317
8 290,748 0.997 291,570 822.22 1,180.80 1.436
9 276,447 0.993 278,320 1,872.75 2,486.60 1.328
10 224,499 0.896 250,571 26,071.84 39,252.98 1.506
Totals
Latest: 5,602,275.93
Ultimate: 5,631,559.51
IBNR: 29,283.58
Mack S.E.: 39,374.90
CV(IBNR): 1.34
> sum(exp(predict(regm,newdata=base))[passe!=TRUE])
[1] 29669.72
> MackChainLadder(P.corp)
MackChainLadder(Triangle = P.corp)
Latest Dev.To.Date Ultimate IBNR Mack.S.E CV(IBNR)
1 2,231,325 1.000 2,231,325 0 0 NaN
2 182,004 0.991 183,683 1,678 1,060 0.632
3 196,756 0.981 200,589 3,833 2,731 0.712
4 205,591 0.964 213,166 7,575 4,534 0.599
5 176,427 0.943 187,096 10,669 6,128 0.574
6 156,175 0.912 171,227 15,053 9,091 0.604
7 138,483 0.867 159,709 21,226 12,717 0.599
8 109,483 0.803 136,359 26,876 17,557 0.653
9 73,108 0.704 103,835 30,727 29,956 0.975
10 24,322 0.532 45,723 21,401 59,465 2.779
Totals
Latest: 3,493,676.24
Ultimate: 3,632,713.12
IBNR: 139,036.88
Mack S.E.: 72,058.66
CV(IBNR): 0.52
> sum(exp(predict(regc,newdata=base))[passe!=TRUE])
[1] 139036.9
Mais le point important est que l'on peut récupérer les résidus, et étudier leur structure de dépendance,
> plot(residuals(regc,type="pearson"),
+ residuals(regm,type="pearson"),pch=19,col="blue")

et pour les amateurs de copules, on trouvera un beau cas de dépendance très très forte
> plot(rank(residuals(regc,type="pearson"))/56,
+ rank(residuals(regm,type="pearson"))/56,pch=19,col="green")

Notons toutefois que si on supprime la première ligne des triangles (
combined all years), on obtient des représentations sensiblement différentes,

et une dépendance également très différente,

malgré l'effet ligne qu iaurait dû être
présent dès la première ligne. Visiblement, cette
première ligne est trop complexe pour être traitée
comme les autres (agrégation de processus de gestion de
sinistres trop différents). Notons que la
corrélation au sens de Spearman (corrélation de
rang) vaut ici
> u=rank(residuals(regc,type="pearson"))/46
> v=rank(residuals(regm,type="pearson"))/46
> cor(u,v)
[1] 0.3051383On peut alors tenter de faire du boostrap en tirant les
paires de résidus à chaque fois, ce qui permet de prendre en compte cette dépendance qui semble exister entre les triangles.
Dans un second temps, on peut aussi introduire de la dépendance
sur les scenarios à venir. Lors des tirages des lois normales
(évoqués
ici), on peut très bien considérer
un
vecteur Gaussien bivarié, avec une corrélation
arbitrairement fixée. La fonction de simulation dans la
méthode Chain-Ladder devient alors,
> CLboot2=function(triangle1,l1,s1,triangle2,l2,s2,rho){
+ m=nrow(triangle1)
+ for(i in 2:m){
+ for(j in (m-i+2):m){
+ E=rmnorm(1, mean=c(triangle1[j,i-1]*l1[i-1],
+ triangle2[j,i-1]*l2[i-1]),
+ varcov=matrix(c(triangle1[j,i-1]*s1[i-1]^2,
+ rho*sqrt(triangle1[j,i-1]*triangle2[j,i-1])*
+ s1[i-1]*s2[i-1],
+ rho*sqrt(triangle1[j,i-1]*triangle2[j,i-1])*
+ s1[i-1]*s2[i-1],triangle2[j,i-1]*s2[i-1]^2),2,2))
+ triangle1[j,i]=E[1,1]
+ triangle2[j,i]=E[1,2]
+ }}
+ ULT1=triangle1[,m]
+ DIAG1=diag(triangle1[,m:1])
+ ULT2=triangle2[,m]
+ DIAG2=diag(triangle2[,m:1])
+ return(c(sum(ULT1-DIAG1),sum(ULT2-DIAG2)))
+ }
(le code est sûrement loin d'être optimal, mais au moins on
comprend ce qui se passe). En prenant un peu de temps, on peut alors
faire tourner différents scenarios, avec des dépendances
plus ou moins fortes,
> regm=glm(Ym~col+lig,data=base,family="poisson")
> regc=glm(Yc~col+lig,data=base,family="poisson")
> Ec=residuals(regc,type="pearson")
> Em=residuals(regm,type="pearson")
> YPc=predict(regc,newdata=base)
> YPm=predict(regm,newdata=base)
> n=9
> base=data.frame(YPc,YPm,lig,col)
> nsim=50000
> PROVISION=rep(NA,nsim)
> PROVISIONm=rep(NA,nsim)
> PROVISIONc=rep(NA,nsim)
> PROVISIONc2=rep(NA,nsim)
> PROVISION2=rep(NA,nsim)
> for(k in 1:nsim){
+ I=sample(1:45,size=n^2,replace=TRUE)
+ simEm = Em[I]
+ simEc = Ec[I]
+ I=sample(1:45,size=n^2,replace=TRUE)
+ simEc2= Ec[I]
+
+ bruitm=simEm*sqrt(exp(YPm))
+ bruitc=simEc*sqrt(exp(YPc))
+ bruitc2=simEc2*sqrt(exp(YPc))
+ INCsm=exp(YPm)+bruitm
+ INCsc=exp(YPc)+bruitc
+ INCsc2=exp(YPc)+bruitc2
+ INCMm=matrix(INCsm,n,n)
+ INCMc=matrix(INCsc,n,n)
+ INCMc2=matrix(INCsc2,n,n)
+ CUMMm=INCMm
+ CUMMc=INCMc
+ CUMMc2=INCMc2
+ for(j in 2:n){CUMMm[,j]=CUMMm[,j-1]+INCMm[,j]
+ CUMMc[,j]=CUMMc[,j-1]+INCMc[,j]
+ CUMMc2[,j]=CUMMc2[,j-1]+INCMc2[,j]}
+
+ PROVISIONm[k]=CLb(CUMMm,lambdam,sigmam)
+ PROVISIONc[k]=CLb(CUMMc,lambdac,sigmac)
+ PROVISIONc2[k]=CLb(CUMMc2,lambdac,sigmac)
+ PROVISION[k]=sum(CLb2(CUMMm,lambdam,sigmam,CUMMc,lambdac,sigmac,0.8))
+ PROVISION2[k]=CLb(CUMMm,lambdam,sigmam)+CLb(CUMMc2,lambdac,sigmac)}
où
l'on récupère les prédictions marginales (y compris en générant des
scénarios indépendants dans la procédure de bootstrap). Ensuite, on
regarde la somme des triangles, tout d'abord en tirant les paires
d'erreurs obtenus par GLM, puis en tirant un vecteur Gaussien lors des
prédiction, et enfin dans le cas complémentement indépendant. On peut
commencer par comparer aux montants Chain Ladder obtenu intialement,
> sum(MackChainLadder(P.mat)$FullTriangle[,n]
+ -diag(MackChainLadder(P.mat)$FullTriangle[n:1,]))
[1] 63895.45
> sum(MackChainLadder(P.corp)$FullTriangle[,n]
+ -diag(MackChainLadder(P.corp)$FullTriangle[n:1,]))
[1] 324070.2
Bref, marginalement, on est sur des choses comparables. Si on regarde maintenant au niveau agrégé,
> sum(MackChainLadder(P.tot)$FullTriangle[,n]
+ -diag(MackChainLadder(P.tot)$FullTriangle[n:1,]))
[1] 415998.4
> mean(PROVISION)
[1] 387822.1
> mean(PROVISION2)
[1] 387989.9Autrement dit, on retrouve la fait que la dépendance n'influence pas le
best estimate (ce qui est rassurant), mais il devrait influcer l'erreur,
> plot(density(PROVISION),col="red")
> lines(density(PROVISION2),col="orange")
On obtient la distribution en
rouge est le cas corrélé, et en
orange la prédiction des provisions sous hypothèse d'indépendance. Et si l'on regarde les quantiles
> quantile(PROVISION,.995)
99.5%
464792.2
> quantile(PROVISION2,.995)
99.5%
452022.6
Bref, on obtient un quantile plus élevé que ce que nous avions en
supposant l'indépendance (ce qui nous conforte dans la justesse de
notre algorithme). Comme la distribution semble plutôt Gaussienne, et pas trop lognormale,
j'ai rajouté la "
distribution" obtenue par Mack en supposant la normalité des provisions sur le triangle cumulé, en
bleu, en trait plein, et en pointillé le cas log-normal.

Bref,
on retrouve l'idée intuitive que provisionner sur des risques
non-homogène apporte d'avantage d'incertitude. Il convient donc d'avoir
les triangles les plus stables possible, puis de les agréger ensuite...