Cette semaine, nous allons finir la modélisation du passif d’une compagnie d’assurance non-vie. Les slides sont en ligne.
Actuariat de l’Assurance Non-Vie #10
La semaine dernière, nous avions commencé la modélisation du passif d’une compagnie d’assurance non-vie. Les slides introductifs sont en ligne.
Données sinistres, devoir Actuariat (ENSAE)
Pour le dernier devoir du cours d’actuariat non-vie, à l’ENSAE, quelques données de sinistres. Pour plusieurs milliers de sinistres (garantie dommages tous accidents issu d’un portefeuille automobile d’une société française) on a les données suivantes, entre 1995 et 2014,
- id_sin: identifiant du sinistres (les 4 premiers caractères représentent l’exercice de survenance)
- an_surv: exercice de survenance
- an_gest: exercice de gestion
- etat_sin: code état du sinistre à la fin de l’exercice de gestion (1=mis en jeu, 2= cloture partielle, 3 = cloture totale, 4 = réouvert et 5 = cloture sans suite)
- pmt: paiements effectués sur le sinistre à la fin de l’exercice de gestion
- rec: recours encaissés sur le sinistre à la fin de l’exercice de gestion
- eval_pmt: évaluation des paiements (réglés+à venir) vue à la fin de l’exercice de gestion
- eval_rec: évaluation des recours (réglés+à venir) vue à la fin de l’exercice de gestion
Pour information complémentaire, on dispose des éléments d’information suivants, sur le portefeuille
Exercice | Volume portefeuille (garanties années) |
Chiffre d’affaires (€) | Nb sinistres survenus dans l’année |
2014 | 303 900 | 45 900 000 | 41 678 |
2013 | 296 500 | 45 700 000 | 41 639 |
2012 | 295 600 | 44 000 000 | 42 815 |
2011 | 292 000 | 41 200 000 | 43 351 |
2010 | 282 900 | 40 700 000 | 45 276 |
2009 | 274 300 | 39 800 000 | 42 779 |
2008 | 263 800 | 41 900 000 | 38 941 |
2007 | 256 200 | 41 200 000 | 38 048 |
2006 | 250 800 | 42 400 000 | 38 603 |
2005 | 239 100 | 41 100 000 | 37 539 |
2004 | 226 800 | 37 900 000 | 34 728 |
2003 | 220 300 | 39 200 000 | 34 491 |
2002 | 211 400 | 35 500 000 | 35 984 |
2001 | 202 400 | 34 800 000 | 36 322 |
2000 | 193 400 | 34 000 000 | 36 186 |
1999 | 184 500 | 32 300 000 | 34 123 |
1998 | 175 500 | 29 800 000 | 30 517 |
1997 | 166 500 | 28 000 000 | 28 432 |
1996 | 157 500 | 26 500 000 | 26 852 |
1995 | 148 600 | 25 700 000 | 26 775 |
La base complète est dans le fichier ensae2017.zip. Nous reviendrons sur les attendus lors du prochain cours.
ENSAE, Devoir 3
Pour le dernier devoir du cours d’actuariat de l’assurance dommage, comme annoncé, je demande de faire un calcul de provisions, sur les données du fichier triangle_ensae.ods.
Je veux un rapport, de 5 pages, proposant un montant de provisions pour les 14 années passées, mais aussi pour la dernière année en cours (pour laquelle 103 733 € ont déjà été payés), avec deux méthodes différentes (montant dit best estimate). Dans un second temps, il faudra utiliser une méthode de simulations pour estimer le quantile à 99.5% des paiements à venir, toutes années confondues. Autrement dit, on veut savoir combien de provisions il convient de constituer pour pouvoir indemniser intégralement les assurés dans 99.5% des cas.
Actuariat de l’Assurance non-Vie #11
Dernier jeu de slides, pour compléter la section sur le provisionnement, avec des méthodes plus économétriques. Les slides sont enfin en ligne.
Je complèterais éventuellement s’il y a des questions spécifiques, sur un point, ou un autre.
Actuariat de l’Assurance non-Vie #10
Dernière ligne droite pour le cours d’actuariat de l’assurance non-vie à l’ensae, avec la modélisation de la vie des sinistres, et du calcul des provisions pour sinistres à payer. Les slides sont en ligne.
Macro vs. Micro Methods in Non-Life Claims Reserving (an Econometric Perspective)
With Mathieu Pigeon, we finalized a short paper on ‘macro’ vs. ‘micro’ methods in claims reserving. The paper is now online on arxiv
Traditionally, actuaries have used run-off triangles to estimate reserve (“macro” models, on agregated data). But it is possible to model payments related to individual claims. If those models provide similar estimations, we investigate uncertainty related to reserves, with “macro” and “micro” models. We study theoretical properties of econometric models (Gaussian, Poisson and quasi-Poisson) on individual data, and clustered data. Finally, application on claims reserving are considered.
La vie d’un sinistre…
Pour illustrer un peu le cours d’actuariat de l’assurance dommage, un petit exemple de ce que peut être une base de sinistres (je serais ravi d’en avoir que je puisse partager avec les étudiants, dans le cours).
Regardons par exemple un sinistre,
> X
On a les dates des opérations, avec des paiements et/ou des réestimations de la charge dossier/dossier. Pour visualiser la vie de ce sinistre, on peut utiliser
> D=as.Date(X[,2],"%d/%m/%Y") > plot(D,X[,7]+cumsum((X[,3])),type="l") > lines(D,cumsum((X[,3])),col="red")
Le sinistre est survenu en 2007. Quand on va constituer le triangle, les montants suivants vont être affecté à la ligne correspondant à l’année de survenance 2007,
> X$year=as.POSIXlt(D)$year+1900 > aggregate(X[,3],by=list(X$year), + function(x) sum(x,na.rm=TRUE)) Group.1 x 1 1997 77.597 2 1998 9568.327 3 1999 1398.074 4 2000 16266.760 5 2001 21955.175 6 2002 106.800
correspondant aux incréments observés sur le graphique ci-dessus.
Actuariat de l’Assurance Non-Vie #11
Pour le onzième (et dernier) chapitre du cours d’actuariat de l’assurance non-vie à l’ENSAE, on va finir d’étudier les modèles de provisionnement (provisions pour sinistres à payer) avec les modèles de régression et les modèles bayésien. Je n’aurais pas le temps d’aborder cette année le provisionnement sur données individuelles. Les slides sont en ligne (la version pdf téléchargeable est comme souvent plus complète que celle sur slideshare)
ENSAE, Devoir 3
Pour le troisième (et dernier) devoir du cours d’actuariat de l’assurance dommage, je vais demander de faire un calcul de provisions à l’aide de 2 méthodes vues en cours, à partir des données du fichier triangle_ensae.ods.
Je veux un rapport, de 5 pages, proposant un montant de provisions pour les 14 années passées, avec deux méthodes différentes. Dans un second temps, il faudra utiliser une méthode de simulations pour estimer le quantile à 99.5% des paiements à venir. Le rapport donnera à la fin le montant de provision complémentaire à constituer pour que la compagnie d’assurance puisse indemniser intégralement ses assurés dans 99.5% des cas.
La dead-line est le dimanche 31 janvier.
Actuariat de l’Assurance Non-Vie #10
Pour le dixième chapitre du cours d’actuariat de l’assurance non-vie à l’ENSAE, on commence les modèles de provisionnement (provisions pour sinistres à payer). Les slides sont en ligne (la version pdf téléchargeable est comme souvent plus complète que celle sur slideshare)
Triangles et provisionnement
Pour le quatrième devoir d’actuariat IARD 2, on va travailler sur des (vrais) triangles de paiements de compagnies d’assurance. Les données ont été collectées par Glenn Meyers et Peng Shi, et sont en ligne sur http://casact.org/research/…. A partir du ‘code’ de groupe, pour récupérer ses données, il faut utiliser la petite fonction écrite dans le fichier
> source("http://freakonometrics.free.fr/codeACT2040-4.txt")
Par exemple, pour le groupe 10,
> extract.triangle(10) $triangle.increments 0 1 2 3 4 5 6 7 8 9 1988 1249 2843 4801 6623 8290 10264 1 431 586 714 1989 946 1983 7024 7415 7771 4489 2431 913 638 NA 1990 1765 2978 5111 10617 6285 6898 9628 1 NA NA 1991 1408 3818 4542 11869 3702 6235 6164 NA NA NA 1992 1647 5981 5220 5475 3620 1196 NA NA NA NA 1993 7566 7425 14475 19624 8500 NA NA NA NA NA 1994 2299 4557 17098 37542 NA NA NA NA NA NA 1995 4959 8582 9462 NA NA NA NA NA NA NA 1996 6063 3644 NA NA NA NA NA NA NA NA 1997 6507 NA NA NA NA NA NA NA NA NA $name [1] "Federal Ins Co Grp" $source [1] "prodliab_pos"
Pour le rapport, je veux avoir le nom de la compagnie, ainsi que la ‘source’ (PP Auto, Workers Compensation,Commercial Auto,Medical Malpractice,Product Liability, ouOther Liability). Pour jouer avec les données, il suffit d’utiliser
> T=extract.triangle(10)$triangle.increments > T 0 1 2 3 4 5 6 7 8 9 1988 1249 2843 4801 6623 8290 10264 1 431 586 714 1989 946 1983 7024 7415 7771 4489 2431 913 638 NA 1990 1765 2978 5111 10617 6285 6898 9628 1 NA NA 1991 1408 3818 4542 11869 3702 6235 6164 NA NA NA 1992 1647 5981 5220 5475 3620 1196 NA NA NA NA 1993 7566 7425 14475 19624 8500 NA NA NA NA NA 1994 2299 4557 17098 37542 NA NA NA NA NA NA 1995 4959 8582 9462 NA NA NA NA NA NA NA 1996 6063 3644 NA NA NA NA NA NA NA NA 1997 6507 NA NA NA NA NA NA NA NA NA
Il s’agit de triangle d’incréments, la ‘colonne’ de gauche étant l’année de survenance. Je n’ai pas testé toutes les bases, merci de me dire rapidement s’il y a des soucis.
Claims reserving (introduction)
Mercredi, on commence la modélisation du passif des compagnies d’assurance IARD. Plus particulièrement, nous parlerons des provisions pour sinistres à payer, ou “provision for claims outstanding (PCO)” selon la terminologie anglaise, i.e. “the estimated total cost of ultimate settlement of all claims incurred before the date of record, whether reported or not, less any amounts already paid out in respect thereof.” Je renvoie à la lecture de Le contrôle de la solvabilité des compagnies d’assurance en ligne sur le site de l’OCDE, pour une vision globale des approches de ces provisions. La SOA avait publié un rapport en 2009, Comparison of Incurred But Not Reported IBNR Methods que j’encourage à lire.
Le livre de référence sur le sujet est celui de Mario Wüthrich et Michael Merz. Les premiers chapitres (correspondant a ce qui sera vu en cours) peuvent être téléchargés sur http://actuaries.ch/…
Nous aborderons mercredi les triangles. Parmi les triangles que nous manipulerons
> source("https://perso.univ-rennes1.fr/arthur.charpentier/bases.R")
qui contient plusieurs fichiers, dont
> 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
ainsi que le triangle évoqué sur http://rworkingparty.wikidot.com/
> OthLiabData = read.csv("http://www.casact.org/research/reserve_data/othliab_pos.csv",header=TRUE, sep=",") > library(ChainLadder) > OL = SumData=ddply(OthLiabData,.(AccidentYear,DevelopmentYear,DevelopmentLag),summarise,IncurLoss=sum(IncurLoss_h1-BulkLoss_h1), + CumPaidLoss=sum(CumPaidLoss_h1), EarnedPremDIR=sum(EarnedPremDIR_h1)) > LossTri = as.triangle(OL, origin="AccidentYear", + dev = "DevelopmentLag", value="IncurLoss") > Year = as.triangle(OL, origin="AccidentYear", + dev = "DevelopmentLag", value="DevelopmentYear") > TRIANGLE=LossTri > TRIANGLE[Year>1997]=NA > TRIANGLE dev origin 1 2 3 4 5 6 7 8 9 10 1988 128747 195938 241180 283447 297402 308815 314126 317027 319135 319559 1989 135147 208767 270979 304488 330066 339871 344742 347800 353245 NA 1990 152400 238665 297495 348826 359413 364865 372436 372163 NA NA 1991 151812 266245 357430 400405 423172 442329 460713 NA NA NA 1992 163737 269170 347469 381251 424810 451221 NA NA NA NA 1993 187756 358573 431410 476674 504667 NA NA NA NA NA 1994 210590 351270 486947 581599 NA NA NA NA NA NA 1995 213141 351363 444272 NA NA NA NA NA NA NA 1996 237162 378987 NA NA NA NA NA NA NA NA 1997 220509 NA NA NA NA NA NA NA NA NA
Les transparents sont en ligne,
Pour les lectures complémentaires, Best Estimates for Reserves de Glen Barnett et Ben Zehnwirth est téléchargeable sur http://casact.org/pubs/…. En 2004 , Ben Zehnwirth, Julie Sims et Mark Shapland ont publié Will Your Next Reserve Increase Be Your Last (en ligne sur http://contingencies.org/janfeb04/…) Lors des prochaines séances, nous aborderons les méthodes par simulations. J’encourage la lecture de about the bootstrap, http://insureware.com/Library/… Sinon, des compléments plus spécifiques sont téléchargeables,
- The Chain Ladder and Tweedie Distributed Claims Data par Greg Taylor
- Chain-Ladder Bias: Its Reason and Meaning par Leigh Joseph Halliwell
- Estimating Predictive Distributions for Loss Reserve Models par Glenn Meyers
- Using a Bayesian Approach for Claims Reserving par Mario Wuthrich
- Obtaining Predictive Distributions for Reserves Which Incorporate Expert Opinion par Richard Verrall
- Loss Reserve Estimates: A Statistical Approach for Determining “Reasonableness” par Mark Shapland
- Munich Chain Ladder: A Reserving Method that Reduces the Gap between IBNR Projections Based on Paid Losses and IBNR Projections Based on Incurred Losses par Gerhard Quarg et Thomas Mack
- The Bornhuetter-Ferguson Principle par Klaus Schmidt et Mathias Zocher
- On the Importance of Dispersion Modeling for Claims Reserving: An Application with the Tweedie Distribution par Jean-Philippe Boucher et Danail Davidov
- Quantifying Uncertainty in Reserve Estimates par Zia Rehman et Stuart Klugman
- Robustifying Reserving par Gary Venter et Dumaria Rulina Tampubolon, voir aussi Robustifying Reserving
- Bootstrap Estimation of the Predictive Distributions of Reserves Using Paid and Incurred Claims par Huijuan Liu et Richard Verrall
- Predictive Distributions for Reserves Which Separate True IBNR and IBNER Claims par Huijuan Liu et Richard Verrall
- The Retrospective Testing of Stochastic Loss Reserve Models, par Glenn Meyers
- A GLM-Based Approach to Adjusting for Changes in Case Reserve Adequacy par Larry Decker
- A Method for Modelling Varying Run-off Evolutions in Claims Reserving par Richard Verrall
- A Nonlinear Regression Model of Incurred But Not Reported Losses par Scott Stelljes
- Back-Testing the ODP Bootstrap of the Paid Chain-Ladder Model with Actual Historical Claims Data par Jessica (Weng Kah) Leong, Shaun Wang et Han Chen
- A Practical Way to Estimate One-Year Reserve Risk par Ira Robbin
- Closed-Form Distribution of Prediction Uncertainty in Chain Ladder Reserving by Bayesian Approach par Ji Yao
- The Prediction Error of Bornhuetter-Ferguson par Thomas Mack
- Stochastic Loss Reserving with the Collective Risk Model par Glenn Meyers
- On the Accuracy of Loss Reserving Methodology par Tapio Boles et Andy Staudt
- Bootstrapping Generalized Linear Models for Development Triangles Using Deviance Residuals par Thomas Hartl
- Fitting a GLM to Incomplete Development Triangles par Thomas Hartl
- Gauss—Markov Loss Prediction in a Linear Model par Alexander Ludwig et Klaus Schmidt
- Anatomy of Actuarial Methods of Loss Reserving par Prakash Narayan
- Bootstrap Modeling: Beyond the Basics par Mark R. Shapland et Jessica (Weng Kah) Leong
- GLM Invariants par Fred Klinker
- The Retrospective Testing of Stochastic Loss Reserve Models par Glenn Meyers et Peng Shi
- A GLM-Based Approach to Adjusting for Changes in Case Reserve Adequacy par Larry Decker
- Testing the Assumptions of Age-to-age Factors, par Gary Venter, voir aussi http://rationalargumentator.com/actuaryguide/…
Poisson regression on non-integers
In the course on claims reserving techniques, I did mention the use of Poisson regression, even if incremental payments were not integers. For instance, we did consider incremental triangles
> source("https://perso.univ-rennes1.fr/arthur.charpentier/bases.R") > 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 53 22 NA NA [4,] 4239 1678 103 NA NA NA [5,] 4929 1865 NA NA NA NA [6,] 5217 NA NA NA NA NA
On those payments, it is natural to use a Poisson regression, to predict future payments
> 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) > matrix(Yp,6,6) [,1] [,2] [,3] [,4] [,5] [,6] [1,] 3155.6 1202.1 49.8 19.1 8.2 21.0 [2,] 3365.6 1282.0 53.1 20.4 8.7 22.3 [3,] 3863.7 1471.8 60.9 23.4 10.0 25.7 [4,] 4310.0 1641.8 68.0 26.1 11.2 28.6 [5,] 4919.8 1874.1 77.6 29.8 12.8 32.7 [6,] 5217.0 1987.3 82.3 31.6 13.5 34.7
and the total amount of reserves would be
> sum(Yp[is.na(Y)==TRUE]) [1] 2426.985
Here, payments were in ‘000 euros. What if they were in ‘000’000 euros ?
> a=1000 > INC/a [,1] [,2] [,3] [,4] [,5] [,6] [1,] 3.209 1.163 0.039 0.017 0.007 0.021 [2,] 3.367 1.292 0.037 0.024 0.010 NA [3,] 3.871 1.474 0.053 0.022 NA NA [4,] 4.239 1.678 0.103 NA NA NA [5,] 4.929 1.865 NA NA NA NA [6,] 5.217 NA NA NA NA NA
We can still run a regression here
> reg=glm((Y/a)~as.factor(D)+as.factor(A),data=base,family=poisson(link="log")) > Yp=predict(reg,type="response",newdata=base) > sum(Yp[is.na(Y)==TRUE])*a [1] 2426.985
and the prediction is exactly the same. Actually, it is possible to change currency, and multiply by any kind of constant, the Poisson regression will return always the same prediction, if we use a log link function,
> homogeneity=function(a=1){ + reg=glm((Y/a)~as.factor(D)+as.factor(A), data=base,family=poisson(link="log")) + Yp=predict(reg,type="response",newdata=base) + return(sum(Yp[is.na(Y)==TRUE])*a) + } > Vectorize(homogeneity)(10^(seq(-3,5))) [1] 2426.985 2426.985 2426.985 2426.985 2426.985 2426.985 2426.985 2426.985 2426.985
The trick here come from the fact that we do like the Poisson interpretation. But GLMs simply mean that we do want to solve a first order condition. It is possible to solve explicitly the first order condition, which was obtained without any condition such that values should be integers. To run a simple code, the intercept should be related to the last value of the matrix, not the first one.
> base$D=relevel(as.factor(base$D),"6") > base$A=relevel(as.factor(base$A),"2006") > reg=glm(Y~as.factor(D)+as.factor(A), data=base,family=poisson(link="log")) > summary(reg) Call: glm(formula = Y ~ as.factor(D) + as.factor(A), family = poisson(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 z value Pr(>|z|) (Intercept) 3.54723 0.21921 16.182 < 2e-16 *** as.factor(D)1 5.01244 0.21877 22.912 < 2e-16 *** as.factor(D)2 4.04731 0.21896 18.484 < 2e-16 *** as.factor(D)3 0.86391 0.22827 3.785 0.000154 *** as.factor(D)4 -0.09254 0.25229 -0.367 0.713754 as.factor(D)5 -0.93717 0.32643 -2.871 0.004092 ** as.factor(A)2001 -0.50271 0.02079 -24.179 < 2e-16 *** as.factor(A)2002 -0.43831 0.02045 -21.433 < 2e-16 *** as.factor(A)2003 -0.30029 0.01978 -15.184 < 2e-16 *** as.factor(A)2004 -0.19096 0.01930 -9.895 < 2e-16 *** as.factor(A)2005 -0.05864 0.01879 -3.121 0.001799 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) 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: 209.52
The first idea is to run a gradient descent, as follows (the starting point will be coefficients from a linear regression on the log of the observations),
> YNA <- Y > XNA=matrix(0,length(Y),1+5+5) > XNA[,1]=rep(1,length(Y)) > for(k in 1:5) XNA[(k-1)*6+1:6,k+1]=k > u=(1:(length(Y))%%6); u[u==0]=6 > for(k in 1:5) XNA[u==k,k+6]=k > YnoNA=YNA[is.na(YNA)==FALSE] > XnoNA=XNA[is.na(YNA)==FALSE,] > beta=lm(log(YnoNA)~0+XnoNA)$coefficients > for(s in 1:50){ + Ypred=exp(XnoNA%*%beta) + gradient=t(XnoNA)%*%(YnoNA-Ypred) + omega=matrix(0,nrow(XnoNA),nrow(XnoNA));diag(omega)=exp(XnoNA%*%beta) + hessienne=-t(XnoNA)%*%omega%*%XnoNA + beta=beta-solve(hessienne)%*%gradient} > beta [,1] [1,] 3.54723486 [2,] 5.01244294 [3,] 2.02365553 [4,] 0.28796945 [5,] -0.02313601 [6,] -0.18743467 [7,] -0.50271242 [8,] -0.21915742 [9,] -0.10009587 [10,] -0.04774056 [11,] -0.01172840
We are not too far away from the values given by R. Actually, it is just fine if we focus on the predictions
> matrix(exp(XNA%*%beta),6,6)) [,1] [,2] [,3] [,4] [,5] [,6] [1,] 3155.6 1202.1 49.8 19.1 8.2 21.0 [2,] 3365.6 1282.0 53.1 20.4 8.7 22.3 [3,] 3863.7 1471.8 60.9 23.4 10.0 25.7 [4,] 4310.0 1641.8 68.0 26.1 11.2 28.6 [5,] 4919.8 1874.1 77.6 29.8 12.8 32.7 [6,] 5217.0 1987.3 82.3 31.6 13.5 34.7
which are exactly the one obtained above. And here, we clearly see that there is no assumption such as “explained variate should be an integer“. It is also possible to remember that the first order condition is the same as the one we had with a weighted least square model. The problem is that the weights are function of the prediction. But using an iterative algorithm, we should converge,
> beta=lm(log(YnoNA)~0+XnoNA)$coefficients > for(i in 1:50){ + Ypred=exp(XnoNA%*%beta) + z=XnoNA%*%beta+(YnoNA-Ypred)/Ypred + REG=lm(z~0+XnoNA,weights=Ypred) + beta=REG$coefficients + } > > beta XnoNA1 XnoNA2 XnoNA3 XnoNA4 XnoNA5 XnoNA6 3.54723486 5.01244294 2.02365553 0.28796945 -0.02313601 -0.18743467
XnoNA7 XnoNA8 XnoNA9 XnoNA10 XnoNA11 -0.50271242 -0.21915742 -0.10009587 -0.04774056 -0.01172840
which are the same values as the one we got previously. Here again, the prediction is the same as the one we got from this so-called Poisson regression,
> matrix(exp(XNA%*%beta),6,6) [,1] [,2] [,3] [,4] [,5] [,6] [1,] 3155.6 1202.1 49.8 19.1 8.2 20.9 [2,] 3365.6 1282.0 53.1 20.4 8.7 22.3 [3,] 3863.7 1471.8 60.9 23.4 10.0 25.7 [4,] 4310.0 1641.8 68.0 26.1 11.2 28.6 [5,] 4919.8 1874.1 77.6 29.8 12.8 32.7 [6,] 5217.0 1987.3 82.3 31.6 13.5 34.7
Again, it works just fine because GLMs are mainly conditions on the first two moments, and numerical computations are based on the first order condition, which has less constraints than the interpretation in terms of a Poisson model.
Overdispersed Poisson et bootstrap
Pour le dernier cours sur les méthodes de provisionnement, on s’est arrête aux méthodes par simulation. Reprenons là où on en était resté au dernier billet où on avait vu qu’en faisant une régression de Poisson sur les incréments, on obtenait exactement le même montant que la méthode Chain Ladder,
> Y [,1] [,2] [,3] [,4] [,5] [,6] [1,] 3209 1163 39 17 7 21 [2,] 3367 1292 37 24 10 NA [3,] 3871 1474 53 22 NA NA [4,] 4239 1678 103 NA NA NA [5,] 4929 1865 NA NA NA NA [6,] 5217 NA NA NA NA NA > y=as.vector(as.matrix(Y)) > base=data.frame(y,ai=rep(2000:2005,n),bj=rep(0:(n-1),each=n)) > reg2=glm(y~as.factor(ai)+as.factor(bj),data=base,family=poisson) > summary(reg2) Call: glm(formula = y ~ as.factor(ai) + as.factor(bj), family = poisson, data = base) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 8.05697 0.01551 519.426 < 2e-16 *** as.factor(ai)2001 0.06440 0.02090 3.081 0.00206 ** as.factor(ai)2002 0.20242 0.02025 9.995 < 2e-16 *** as.factor(ai)2003 0.31175 0.01980 15.744 < 2e-16 *** as.factor(ai)2004 0.44407 0.01933 22.971 < 2e-16 *** as.factor(ai)2005 0.50271 0.02079 24.179 < 2e-16 *** as.factor(bj)1 -0.96513 0.01359 -70.994 < 2e-16 *** as.factor(bj)2 -4.14853 0.06613 -62.729 < 2e-16 *** as.factor(bj)3 -5.10499 0.12632 -40.413 < 2e-16 *** as.factor(bj)4 -5.94962 0.24279 -24.505 < 2e-16 *** as.factor(bj)5 -5.01244 0.21877 -22.912 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) 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: 209.52 Number of Fisher Scoring iterations: 4 > base$py2=predict(reg2,newdata=base,type="response") > round(matrix(base$py2,n,n),1) [,1] [,2] [,3] [,4] [,5] [,6] [1,] 3155.7 1202.1 49.8 19.1 8.2 21.0 [2,] 3365.6 1282.1 53.1 20.4 8.8 22.4 [3,] 3863.7 1471.8 61.0 23.4 10.1 25.7 [4,] 4310.1 1641.9 68.0 26.1 11.2 28.7 [5,] 4919.9 1874.1 77.7 29.8 12.8 32.7 [6,] 5217.0 1987.3 82.4 31.6 13.6 34.7 > sum(base$py2[is.na(base$y)]) [1] 2426.985
Le plus intéressant est peut être de noter que la loi de Poisson présente probablement trop peu de variance…
> reg2b=glm(y~as.factor(ai)+as.factor(bj),data=base,family=quasipoisson) > summary(reg2) Call: glm(formula = y ~ as.factor(ai) + as.factor(bj), family = quasipoisson, data = base) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 8.05697 0.02769 290.995 < 2e-16 *** as.factor(ai)2001 0.06440 0.03731 1.726 0.115054 as.factor(ai)2002 0.20242 0.03615 5.599 0.000228 *** as.factor(ai)2003 0.31175 0.03535 8.820 4.96e-06 *** as.factor(ai)2004 0.44407 0.03451 12.869 1.51e-07 *** as.factor(ai)2005 0.50271 0.03711 13.546 9.28e-08 *** as.factor(bj)1 -0.96513 0.02427 -39.772 2.41e-12 *** as.factor(bj)2 -4.14853 0.11805 -35.142 8.26e-12 *** as.factor(bj)3 -5.10499 0.22548 -22.641 6.36e-10 *** as.factor(bj)4 -5.94962 0.43338 -13.728 8.17e-08 *** as.factor(bj)5 -5.01244 0.39050 -12.836 1.55e-07 *** --- 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
Mais on en reparlera dans un instant. Ensuite, on avait commencé à regarder erreurs commises, sur la partie supérieure du triangle. Classiquement, par construction, les résidus de Pearson sont de la forme
On avait vu dans le cours de tarification que la variance au dénominateur pouvait être remplacé par le prévision, puisque dans un modèle de Poisson, l’espérance et la variance sont identiques. Donc on considérait
> base$erreur=(base$y-base$py2)/sqrt(base$py2) > round(matrix(base$erreur,n,n),1) [,1] [,2] [,3] [,4] [,5] [,6] [1,] 0.9 -1.1 -1.5 -0.5 -0.4 0 [2,] 0.0 0.3 -2.2 0.8 0.4 NA [3,] 0.1 0.1 -1.0 -0.3 NA NA [4,] -1.1 0.9 4.2 NA NA NA [5,] 0.1 -0.2 NA NA NA NA [6,] 0.0 NA NA NA NA NA
Le soucis est que si est – asymptotiquement – un bon estimateur pour , ce n’est pas le cas à distance finie, car on a alors un estimateur biaisé pour la variance, et donc la variance des résidus n’a que peu de chances d’être de variance unitaire. Aussi, il convient de corriger l’estimateur de la variance, et on pose alors
qui sont alors les résidus de Pearson tel qu’on doit les utiliser.
> E=base$erreur[is.na(base$y)==FALSE]*sqrt(21/(21-11)) > E [1] 1.374976e+00 3.485024e-02 1.693203e-01 -1.569329e+00 1.887862e-01 [6] -1.459787e-13 -1.634646e+00 4.018940e-01 8.216186e-02 1.292578e+00 [11] -3.058764e-01 -2.221573e+00 -3.207593e+00 -1.484151e+00 6.140566e+00 [16] -7.100321e-01 1.149049e+00 -4.307387e-01 -6.196386e-01 6.000048e-01 [21] -8.987734e-15 > boxplot(E,horizontal=TRUE)
En rééchantillonnant dans ces résidus, on peut alors générer un pseudo triangle. Pour des raisons de simplicités, on va générer un peu rectangle, et se restreindre à la partie supérieure,
> Eb=sample(E,size=36,replace=TRUE) > Yb=base$py2+Eb*sqrt(base$py2) > Ybna=Yb > Ybna[is.na(base$y)]=NA > Tb=matrix(Ybna,n,n) > round(matrix(Tb,n,n),1) [,1] [,2] [,3] [,4] [,5] [,6] [1,] 3115.8 1145.4 58.9 46.0 6.4 26.9 [2,] 3179.5 1323.2 54.5 21.3 12.2 NA [3,] 4245.4 1448.1 61.0 7.9 NA NA [4,] 4312.4 1581.7 68.7 NA NA NA [5,] 4948.1 1923.9 NA NA NA NA [6,] 4985.3 NA NA NA NA NA
Cette fois, on a un nouveau triangle ! on va alors pouvoir faire plusieurs choses,
- compléter le triangle para la méthode Chain Ladder, c’est à dire calculer les montants moyens que l’on pense payer dans les années futures
- générer des scénarios de paiements pour les années futurs, en générant des paiements suivant des lois de Poisson (centrées sur les montants moyens que l’on vient de calculer)
- générer des scénarios de paiements avec des lois présentant plus de variance que la loi de Poisson. Idéalement, on voudrait simuler des lois qusi-Poisson, mais ce ne sont pas de vraies lois. Par contre, on peut se rappeler que dans ce cas, la loi Gamma devrait donner une bonne approximation.
Pour ce dernier point, on va utiliser le code suivant, pour générer des quasi lois,
> 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) + }
Je renvois aux diverses notes de cours pour plus de détails sur la justification, ou à un vieux billet. On va alors faire une petite fonction, qui va soit somme les paiements moyens futurs, soit sommer des générations de scénarios de paiements, à partir d’un triangle,
> CL=function(Tri){ + y=as.vector(as.matrix(Tri)) + base=data.frame(y,ai=rep(2000:2005,n),bj=rep(0:(n-1),each=n)) + reg=glm(y~as.factor(ai)+as.factor(bj),data=base,family=quasipoisson) + py2=predict(reg,newdata=base,type="response") + pys=rpois(36,py2) + pysq=rqpois(36,py2,phi=3.18623) + return(list( + cl=sum(py2[is.na(base$y)]), + sc=sum(pys[is.na(base$y)]), + scq=sum(pysq[is.na(base$y)]))) + }
Reste alors à générer des paquets de triangles. Toutefois, il est possible de générer des triangles avec des incréments négatifs. Pour faire simple, on mettra des valeurs nulles quand on a un paiement négatif. L’impact sur les quantiles sera alors (a priori) négligeable.
> for(s in 1:1000){ + Eb=sample(E,size=36,replace=TRUE)*sqrt(21/(21-11)) + Yb=base$py2+Eb*sqrt(base$py2) + Yb=pmax(Yb,0) + scY=rpois(36,Yb) + Ybna=Yb + Ybna[is.na(base$y)]=NA + Tb=matrix(Ybna,6,6) + C=CL(Tb) + VCL[s]=C$cl + VR[s]=C$sc + VRq[s]=C$scq + }
Si on regarde la distribution du best estimate, on obtient
> hist(VCL,proba=TRUE,col="light blue",border="white",ylim=c(0,0.003)) > boxplot(VCL,horizontal=TRUE,at=.0023,boxwex = 0.0006,add=TRUE,col="light green") > D=density(VCL) > lines(D) > I=which(D$x<=quantile(VCL,.05)) > polygon(c(D$x[I],rev(D$x[I])),c(D$y[I],rep(0,length(I))),col="blue",border=NA) > I=which(D$x>=quantile(VCL,.95)) > polygon(c(D$x[I],rev(D$x[I])),c(D$y[I],rep(0,length(I))),col="blue",border=NA)
Mais on peut aussi visualiser des scénarios basés sur des lois de Poisson (équidispersé) ou des scénarios de lois quasiPoisson (surdispersées), ci-dessous
Dans ce dernier cas, on peut en déduire le quantile à 99% des paiements à venir.
> quantile(VRq,.99) 99% 2855.01
Il faut donc augmenter le montant de provisions de l’ordre 15% pour s’assurer que la compagnie pourra satisfaire ses engagements dans 99% des cas,
> quantile(VRq,.99)-2426.985 99% 428.025