Tag Archives: IBNR

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.

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.

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,

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

https://latex.codecogs.com/gif.latex?\varepsilon_i=\frac{Y_i-\widehat{Y}_i}{\sqrt{\text{Var}(Y_i)}}

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

https://latex.codecogs.com/gif.latex?\varepsilon_i=\frac{Y_i-\widehat{Y}_i}{\sqrt{\widehat{Y}_i}}

> 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 https://latex.codecogs.com/gif.latex?\widehat{Y}_i est – asymptotiquement – un bon estimateur pour https://latex.codecogs.com/gif.latex?\text{Var}(Y_i), 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

https://latex.codecogs.com/gif.latex?\varepsilon_i=\sqrt{\frac{n}{n-k}}\cdot\frac{Y_i-\widehat{Y}_i}{\sqrt{\widehat{Y}_i}}

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,

  1. 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
  2. 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)
  3. 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

Chain Ladder, avec R

Un billet rapide pour mettre en ligne des parties du code tapé en cours, mercredi dernier. On avait commencé par convertir la feuille du classeur excel en un fichier texte, pour faciliter la lecture,

> setwd("C:\\Users\\savsalledecours\\Desktop")
> triangle=read.table("exACT2040.csv",header=TRUE,sep=";")	
> triangle
  ANNEE   X0   X1   X2   X3   X4   X5
1  2000 3209 4372 4411 4428 4435 4456
2  2001 3367 4659 4696 4720 4730   NA
3  2002 3871 5345 5398 5420   NA   NA
4  2003 4239 5917 6020   NA   NA   NA
5  2004 4929 6794   NA   NA   NA   NA
6  2005 5217   NA   NA   NA   NA   NA

L’idée – quand on importe un triangle – est de récupérer une base sous la forme précédente, avec des valeurs manquantes dans la partie inférieure du triangle (on verra l’intérêt quand on fait une régression). On avait ensuite calculé les facteurs de transition, et en même temps complété le triangle,

> T=triangle[,2:7]
> rownames(T)=triangle$ANNEE
> T2=T
> n=ncol(T)
> L=rep(NA,n-1)
> for(j in 1:(n-1)){
+ L[j]=sum(T[1:(n-j),j+1])/sum(T[1:(n-j),j])
+ T2[(n-j+1):n,j+1]=L[j]*T2[(n-j+1):n,j]
+ }

Les facteurs de transition sont ici,

> L
[1] 1.380933 1.011433 1.004343 1.001858 1.004735

et le triangle complété

> T2
       X0       X1       X2       X3       X4       X5
2000 3209 4372.000 4411.000 4428.000 4435.000 4456.000
2001 3367 4659.000 4696.000 4720.000 4730.000 4752.397
2002 3871 5345.000 5398.000 5420.000 5430.072 5455.784
2003 4239 5917.000 6020.000 6046.147 6057.383 6086.065
2004 4929 6794.000 6871.672 6901.518 6914.344 6947.084
2005 5217 7204.327 7286.691 7318.339 7331.939 7366.656

Le montant de provision est ici en faisant la différence entre la charge ultime (dans la dernière colonne) et les derniers paiements observés (sur la seconde diagonale)

> CU=T2[,n]
> Pat=diag(as.matrix(T2[,n:1]))
> Ri=CU-Pat
> R=sum(Ri)

soit, numériquement

> R
[1] 2426.985

On avait alors vu que l’on pouvait calculer un tail factor, en supposant une décroissance exponentielle des facteurs de transition, et on rajoutait alors une colonne correspondant au montant ultime, par année d’accident,

> logL=log(L-1)
> t=1:5
> b=data.frame(logL,t)
> reg=lm(logL~t,data=b)
> logLp=predict(reg,newdata=data.frame(t=6:100))
> Lp=exp(logLp)+1
> Linf=prod(Lp)
> T3=T2
> T3$Xinf=T3$X5*Linf

On a ici

> T3
       X0       X1       X2       X3       X4       X5     Xinf
2000 3209 4372.000 4411.000 4428.000 4435.000 4456.000 4459.149
2001 3367 4659.000 4696.000 4720.000 4730.000 4752.397 4755.755
2002 3871 5345.000 5398.000 5420.000 5430.072 5455.784 5459.639
2003 4239 5917.000 6020.000 6046.147 6057.383 6086.065 6090.366
2004 4929 6794.000 6871.672 6901.518 6914.344 6947.084 6951.993
2005 5217 7204.327 7286.691 7318.339 7331.939 7366.656 7371.862

(je laisse reprendre le code pour calculer le montant de provisions). Enfin, on avait montré comment utiliser une régression pondérée, pour calculer les facteurs de transition,

> T4=as.matrix(T$X0,n,1)
> for(j in 1:(n-1)){
+ Y=T[,j+1]
+ X=T[,j]
+ base=data.frame(X,Y)
+ reg=lm(Y~0+X,weights=1/X)
+ T4=cbind(T4,
+ predict(reg,
+ newdata=data.frame(X=T4[,j]
+ )))
+ }

Ce qui donnait la même projection que la méthode Chain Ladder

> T4
  [,1]     [,2]     [,3]     [,4]     [,5]     [,6]
1 3209 4431.414 4482.076 4501.543 4509.909 4531.263
2 3367 4649.601 4702.758 4723.184 4731.961 4754.367
3 3871 5345.591 5406.705 5430.188 5440.279 5466.039
4 4239 5853.775 5920.698 5946.414 5957.464 5985.673
5 4929 6806.619 6884.435 6914.337 6927.186 6959.986
6 5217 7204.327 7286.691 7318.339 7331.939 7366.656

La suite mercredi prochain, même si on risque d’aller très vite sur la méthode de Mack (et les calculs d’erreur quadratique moyenne pour arriver à la régression de Poisson). A suivre donc…

Solvabilité et provisionnement

Mercredi, nous allons aborder en cours les aspects de solvabilité 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.

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

Triangles et provisionnement

La première partie du cours sur le provisionnement (calcul des provisions pour sinistres à payer) aura lieu dans 10 jours. Les transparents sont en ligne ici, et portent sur la construction des triangles de paiements. La méthode chain ladder (et la formalisation proposée par Thomas Mack) ainsi que les extensions seront présentées. La seconde partie portera sur les méthodes basées sur la régression de Poisson.