Arthur Charpentier

Aller au contenu | Aller au menu | Aller à la recherche

Keyword - bootstrap

Fil des billets

mercredi 31 mars 2010

Triangles multiples et dépendance

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 http://perso.univ-rennes1.fr/arthur.charpentier/latex/ln-tr-01.png
alors les deux premiers moments vérifient
http://perso.univ-rennes1.fr/arthur.charpentier/latex/ln-tr-02.png
et
http://perso.univ-rennes1.fr/arthur.charpentier/latex/ln-tr-03.png
ce qui s’inverse simplement, et permet d’écrire
http://perso.univ-rennes1.fr/arthur.charpentier/latex/ln-tr-04.png
et
http://perso.univ-rennes1.fr/arthur.charpentier/latex/ln-tr-05.png
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+LC

On 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.9
qu’il faudrait comparer à
> qlnorm(.9,mu,sqrt(sigma2))
[1] 365395.1

Bref, 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).
  • Approche par les GLM
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.3051383

On 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.9

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

lundi 30 mars 2009

Talk for students in actuarial science, Belo Horizonte

Talk for the students of the undergraduate program in Actuarial Science in Belo Horizonte, Tuesday, on "claims reserving (IBNR) and statistical issues". Slides can be found here.

Based on recent studies, such as the report by Swiss Re (here), the CAS Working Group on Quantifying Variability in Reserve Estimates (here), or the ROC/GIRO report on Best Estimates and Reserving Uncertainty (there), I will stress on important issues in statistical estimation of claim reserves. Note that a report by Tillinghast can also be found here.

I will also mention statutory elements, such as MCR and SCR, i.e. insurance company are now requiered not only to provide a best estimate for there IBNR, but also have to quantify uncertainty.

I will also mention practical issues for actuaries (and statisticians), using package ChainLadder from R, and showing how it becomes possible to interface R and MSExcel (see e.g. slides here, based on RExcel-Addin).

jeudi 12 février 2009

Satellite Meeting: New Trends in Actuarial Practice, Rio de Janeiro

Lecture on "Claims Reserving with R" at the Satellite Meetings in April 2009. For motivations on claims reservings in non-life insurance, see the Swiss Re technical document (here).

mardi 20 janvier 2009

Econométrie: monte carlo et bootstrap

Un dernier commentaire sur le cours d’économétrie, toujours pour reprendre des erreurs commises dans les projets. Une partie du cours était dédiée à l’utilisation des méthodes de Monte Carlo (cf également un ancien billet).

Si on cherche à calculer le biais d’un estimateur, au lieu d’un calcul numérique de l’espérance (qui n’est parfois pas possible analytiquement), il est possible de simuler des échantillons, puis d’approcher l’espérance par une moyenne empirique (via la loi des grands nombres). L’idée n’est donc pas de simuler un échantillon tel que les résidus soient gaussiens, mais plusieurs centaines, puis de regarder le comportement moyen sur ces scénarios.

Pour aller un peu plus loin, je renvoie vers des notes de cours sur internet, en particulier celles de Jean-Marie Dufour,

    1. "General considerations on finite-sample inference in econometrics and statistics", 2002. Slides: PS; PDF
    2. "Finite-sample inference and bounds methods in econometrics and statistics", 2002. Slides: PS; PDF
    3. "Finite-sample inference in econometrics and statistics", 2006. Slides: PS; PDF
    4. "Finite-sample inference, weak identification and macroeconometrics", 2006. Slides: PS; PDF

Sinon Bernard Salanié évoque aussi l’utilisation des méthodes de simulations en théorie des tests, en notant que les propriétés à distances finies (petits échantillons) peuvent être sensiblement différentes des théories asymptotiques de la plupart des tests usuels.

mercredi 10 décembre 2008

Formation sur les méthodes de provisionnement

Formation Caritat sur les méthodes de provisionnement en assurance non-vie.

La formation aura lieu sur deux jours. En guise d’introduction, je renvois à un l’article paru en 2004 dans l’Argus de l’Assurance. Sinon les slides que je vais utiliser son maintenant accessibles. Ils reposent sur une feuille excel et des codes R (ces derniers seront donnés lors de la formation).

1. Introduction aux techniques de provisionnement : Vers une classification des modèles
2. L’approche Chain-Ladder et quelques modèles déterministes

L’approche Chain-Ladder
Le problème de l’inflation dans les triangles
3/ Méthodes récursives model-free : l’approche de Mack
Erreur par année de survenance
Erreur par triangle
Impact sur les boni-mali
4/ Modèles factoriels

Quelles variables explicatives ?
Régression lognormale
Régression poissonnienne
Modélisation GLM
Uniformisation par les modèles Tweedie
Introduction aux modèles GAM
5/ Méthodes de simulation
Le bootstrap, ou simulation non paramétrique
Les simulations paramétriques
Calculs de VaR ou de TVaR
6/ Approche bayésienne
Introduire un avis d’expert
L’approche bayésienne du provisionnement
7/ La difficulté du choix de modèle
8/ Agréger les triangles
9/ La problématique du provisionnement dans l’optique Solvabilité 2

Sur le plan pratique, les premiers modèles seront développés sous Excel, directement, et ensuite nous verrons l’utilisation de R. Un fichier source pour R sera bientôt mis en ligne, ainsi qu’un fichier xls. Pour motiver un peu (si besoin était) l’utilisation de R pour les actuaires, je renvois simplement aux slides présentés au Meeting de la CAS mi-novembre. Nous évoquerons également l’interface possible Excel/R, i.e. utiliser les fonctionalités de R depuis Excel, via le RExcel-Addin. Sinon, côté biblio, je renvoie à un document du GIRO 2007 sur une comparaison entre les méthodes de provisionnement.

- page 1 de 2