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 #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.
Être – en même temps – moyen et extraordinaire
Juste avant 1850, Adolphe Quételet publiait Sur l’appréciation des documents statistiques, et en particulier sur l’application des moyennes, et introdusait alors le concept d’homme moyen (l’ouvrage est cité dans La théorie des moyennes et son emploi dans les sciences d’observation, de Joseph Lottin, paru en 1909 – je n’ai pas trouvé en ligne cet ouvrage de Quételet). Ce concept sera détaillé encore davantage dans Anthropométrie, ou mesure des différentes facultés de l’homme, concept qui sera abondament illustré
Mais si la moyenne est un concept fondamental en statistique, il convient de garder en mémoire que la moyenne n’existe peut être pas, au sein de la population. C’est ce que rappelle la fameuse ‘blague’ “un mec qui à la tête dans le four et les pieds dans le frigo, en moyenne, il est bien“. La moyenne peut être quelquechose d’extraordinaire (au sens premier du terme).
Continue reading Être – en même temps – moyen et extraordinaire
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.
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.
Visualizing Inequalities in a 3-Person Economy
Yesterday, in the course on inequalities, I mentioned (too) briefly the 3-person Economy. I wanted to spend some time in a short post on visualisations of inequalities in such a context. As mentioned in the slides, it is possible to use a ternary plot representation. In the case where we believe that the scale independence principle makes sense, i.e. . A distribution of incomes can be represented as a barycenter in an equilateral triangle (also called de Finetti diagram). The midpoint is the equal situation: the three agents share the same wealth. Because of the scale independence property, we can look at distribution of wealth on the simplex. A wealth distribution is a vector
where each component is one of the (red) distance below. A is on top of the triangle, and the vertical distance is proportional to the wealth of A. The closer to the bottom line, the poorer A is.
To visualize this distribution of wealth, we can use the trifield package. To add the point and the three segments, the code is
tripoint=function(s){ p=s/sum(s) p1=c(0,s[2]+s[1]/2,s[3]+s[1]/2)/sum(s) p2=c(s[1]+s[2]/2,0,s[3]+s[2]/2)/sum(s) p3=c(s[1]+s[3]/2,s[2]+s[3]/2,0)/sum(s) C <- abc2xy(matrix(p,1,3)) points(C,pch=19,col="red",cex=2) C1 <- abc2xy(matrix(p1,1,3)) C2 <- abc2xy(matrix(p2,1,3)) C3 <- abc2xy(matrix(p3,1,3)) segments(C1[1],C1[2],C[1],C[2],lwd=2,col="red") segments(C2[1],C2[2],C[1],C[2],lwd=2,col="red") segments(C3[1],C3[2],C[1],C[2],lwd=2,col="red") }
For instance, to visualize the equal case (inequality indices are defined as a distance to this situation)
tripoint(c(1,1,1))
For a case where there is inequality, use for instance
tripoint(c(1,2,3))
Continue reading Visualizing Inequalities in a 3-Person Economy
Bootstrap et régression
Lors du dernier cours, on a évoqué l’utilisation du bootstrap pour obtenir des intervalles de confiance sur des prévisions. Je mets en ligne les codes tapés en cours (très sommairement commentés, je peux renvoyer vers des vieux billets du cours ACT6420 pour des compléments). On va travailler sur ma base préférée pour évoquer la régression linéaire (avant de parler triangles de provisionnement, revenons cinq minutes sur des choses simples).
> plot(cars) > reg=lm(dist~speed,data=cars) > abline(reg,col="red") > n=nrow(cars) > x=21 > points(x,predict(reg,newdata= data.frame(speed=x)),pch=19,col="red")
On cherche ici à faire une prédiction en un point. Comme rappelé en cours (mais aussi dans le cours de modèles de prévision), quand on veut donner un intervalle de confiance pour la prévision, il convient de distinger l’intervalle de confiance pour le prédicteur (qui va dépendre de l’erreur d’estimation des paramétres) et l’intervalle de confiance pour une potentielle valeur (on peut parler de génération de scénarios, qui va dépendre en plus de l’erreur de modèle, c’est à dire de la dispersion des résidus). Commençons par l’intervalle de confiance sur la prédiction, sur le best estimate comme on dit en provisionnement
> Yx=rep(NA,500)
> B=matrix(NA,500,2)
> for(s in 1:500){
+ indice=sample(1:n,size=n,
+ replace=TRUE)
+ base=cars[indice,]
+ #points(base,pch=3)
+ reg=lm(dist~speed,data=base)
+ abline(reg,col="light blue")
+ points(x,predict(reg,newdata=data.frame(speed=x)),pch=19,col="blue")
+ Yx[s]=predict(reg,newdata=data.frame(speed=x))
+ B[s,]=coefficients(reg)
+ }
Les valeurs bleues sont ici des prévisions possibles, obtenues en rééchantillonnant dans notre base d’observations. Pour rappel, l’intervalle de confiance (à 90%), sous hypothèse de normalité des résidus (et donc des estimateurs de la pente et de la constante de la droite de régression) s’obtient de la manière suivante
> reg=lm(dist~speed,data=cars) > U=predict(reg,interval ="confidence", + level=.9,newdata= + data.frame(speed=0:30)) > lines(0:30,U[,2],col="red",lwd=2) > lines(0:30,U[,3],col="red",lwd=2)
On peut comparer ici la distribution des valeurs obtenues sur nos 500 jeux de données générées, et comparer les quantiles empiriques, avec les quantiles sous hypothèse de normalité,
> hist(Yx,proba=TRUE,col="light blue",border="white") > boxplot(Yx,horizontal=TRUE,at=.07,boxwex = 0.02,add=TRUE,col="light green") > abline(v=U[x+1,2:3],col="red",lwd=2) > D=density(Yx) > lines(D) > I=which(D$x<=quantile(Yx,.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(Yx,.95)) > polygon(c(D$x[I],rev(D$x[I])),c(D$y[I],rep(0,length(I))),col="blue",border=NA)
On peut noter que les ordres de grandeur sont comparables.
> reg=lm(dist~speed,data=cars)
> quantile(Yx,c(.05,.95))
5% 95%
58.63689 70.31281
> predict(reg,interval ="confidence",
+ level=.9,newdata=data.frame(speed=x))
fit lwr upr
1 65.00149 59.65934 70.34364
Regardons maintenant l’autre type d’intervalle de confiance, sur la valeur possible de la variable d’intérêt. Cette fois, en plus de tirer des nouveaux échantillons et calculer des prédictions, on va en plus rajouter un bruit à chaque tirage, qui permettra d’obtenir une valeur possible.
> Yx=rep(NA,500)
> for(s in 1:500){
+ indice=sample(1:n,size=n,
+ replace=TRUE)
+ base=cars[indice,]
+ #points(base,pch=3)
+ reg=lm(dist~speed,data=base)
+ erreur=residuals(reg)
+ #abline(reg,lty=2)
+ E=sample(erreur,size=1)
+ Yx[s]=predict(reg,newdata=data.frame(speed=x))+E
+ points(x,Yx[s],pch=19,col="red")
+ }
Là encore, on peut comparer (graphiquement pour commencer) les valeurs obtenues par rééchantillonnage, et celle obtenues sous hypothèse de normalité,
> hist(Yx,proba=TRUE,col="light blue",border="white") > boxplot(Yx,horizontal=TRUE,at=.025,boxwex = 0.005,add=TRUE,col="light green") > abline(v=U[2:3],col="red",lwd=2) > D=density(Yx) > lines(D) > I=which(D$x<=quantile(Yx,.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(Yx,.95)) > polygon(c(D$x[I],rev(D$x[I])),c(D$y[I],rep(0,length(I))),col="blue",border=NA)
Ce qui donne, numériquement, les comparaisons suivantes
> quantile(Yx,c(.05,.95))
5% 95%
44.43468 96.01357
> (U=predict(reg,interval ="prediction",level=.9,newdata=data.frame(speed=x)))
fit lwr upr
1 67.63136 45.16967 90.09305
On observe cette fois une légère asymméytrie vers la droite. Manifestement, on ne peut pas supposer les résidus Gaussien, car il y a plus de grandes valeurs positives, que négatives. Ce qui fait du sens compte tenu de la nature des données (une distance ne peut être négative).
On avait ensuite commencé à discuter de l’utilisation des modèles de régression en provisionnement. Afin d’avoir des données présentant de l’indépendance, on avait rappelé qu’il fallait travailler avec les incréments de paiments, et non pas les paiements cumulés.
> T
[,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
> n=ncol(T)
> Y=T
> Y[,2:n]=T[,2:n]-
+ T[,1:(n-1)]
> 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
On peut alors constituer une base de données, avec comme variables explicatives la ligne et la colonne.
> y=as.vector(as.matrix(Y))
> base=data.frame(
+ y,
+ ai=rep(2000:2005,n),
+ bj=rep(0:(n-1),each=n))
>
> head(base,12)
y ai bj
1 3209 2000 0
2 3367 2001 0
3 3871 2002 0
4 4239 2003 0
5 4929 2004 0
6 5217 2005 0
7 1163 2000 1
8 1292 2001 1
9 1474 2002 1
10 1678 2003 1
11 1865 2004 1
12 NA 2005 1
> tail(base,12)
y ai bj
25 7 2000 4
26 10 2001 4
27 NA 2002 4
28 NA 2003 4
29 NA 2004 4
30 NA 2005 4
31 21 2000 5
32 NA 2001 5
33 NA 2002 5
34 NA 2003 5
35 NA 2004 5
36 NA 2005 5
On peut alors commencer par utiliser le modèle Regression models based on log-incremental payments de Stavros Christofides, basé sur une modélisation lognormale, introduite initialement par Etienne de Vylder en 1978 (Markus en parle, en trois parties, sur son blog http://lamages.blogspot.ca/Barnett%20Zehnwirth)
> reg1=lm(log(y)~
+ as.factor(ai)+
+ as.factor(bj),data=base)
> summary(reg1)
Call:
lm(formula = log(y) ~ as.factor(ai) + as.factor(bj), data = base)
Residuals:
Min 1Q Median 3Q Max
-0.26374 -0.05681 0.00000 0.04419 0.33014
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.9471 0.1101 72.188 6.35e-15 ***
as.factor(ai)2001 0.1604 0.1109 1.447 0.17849
as.factor(ai)2002 0.2718 0.1208 2.250 0.04819 *
as.factor(ai)2003 0.5904 0.1342 4.399 0.00134 **
as.factor(ai)2004 0.5535 0.1562 3.543 0.00533 **
as.factor(ai)2005 0.6126 0.2070 2.959 0.01431 *
as.factor(bj)1 -0.9674 0.1109 -8.726 5.46e-06 ***
as.factor(bj)2 -4.2329 0.1208 -35.038 8.50e-12 ***
as.factor(bj)3 -5.0571 0.1342 -37.684 4.13e-12 ***
as.factor(bj)4 -5.9031 0.1562 -37.783 4.02e-12 ***
as.factor(bj)5 -4.9026 0.2070 -23.685 4.08e-10 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1753 on 10 degrees of freedom
(15 observations deleted due to missingness)
Multiple R-squared: 0.9975, Adjusted R-squared: 0.9949
F-statistic: 391.7 on 10 and 10 DF, p-value: 1.338e-11
> base$py=exp(predict(reg1,
+ newdata=base)+summary(reg1)$sigma^2/2)
> round(matrix(base$py,n,n),1)
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 2871.2 1091.3 41.7 18.3 7.8 21.3
[2,] 3370.8 1281.2 48.9 21.5 9.2 25.0
[3,] 3768.0 1432.1 54.7 24.0 10.3 28.0
[4,] 5181.5 1969.4 75.2 33.0 14.2 38.5
[5,] 4994.1 1898.1 72.5 31.8 13.6 37.1
[6,] 5297.8 2013.6 76.9 33.7 14.5 39.3
> sum(base$py[is.na(base$y)])
[1] 2481.857
On obtient un montant un peu différent de celui obtenu par la méthode Chain Ladder, mais néanmoins comparable. On peut aussi tenter une régression de Poisson (avec un lien logarithmique), comme suggéré en 1975 par Hachemeister et Stanard,
> 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)
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) 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
La prédiction coïncide avec l’estimateur obtenu par la méthode Chain Ladder. Le lien avec les méthodes de biais minimal a été établi par Klaus Schmidt et Angela Wünsche en 1998, dans Chain ladder, marginal sum and maximum likelihood estimation. La semaine prochaine, on parlera des méthodes de bootstrap pour obtenir des intervalles de confiance, ou des quantiles, sur les montants de réserve. Je ne sais pas si j’aurais le temps de taper des transparents, je préfère, sur cette partie du cours taper au fur et à mesure, et écrire au tableau. Je renvoie au Chapitre 3 du livre avec Christophe Dutang – en ligne sur http://cran.r-project.org/doc/contrib/ – pour le détail. C’est le code que je tape en cours, tout en essayant de répondre aussi aux questions.
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.
Pigeonholes and triangles
Once again, there was a nice maths puzzle on http://www.futilitycloset.com/ last week (but without further reference). The question was the following, “Five points are located in an equilateral triangle with 10-inch sides (or on its perimeter). What’s the maximum distance between the two closest points?” Actually, this is simply an application of Dirichlet’s pigeonholes theorem, as mentioned in the answer of the puzzle, “Connect the midpoints of the triangle’s sides to make four smaller triangles. Because there are five points, two of them must fall within one of these triangles. The maximum distance between these two is 5 inches.”
Thus, with Dirichlet’s pigeonholes theorem, we know not only the maximal minimum distance, but also where points must be (on corners of inside triangles). Here, there might be two possibilities (with also the different shapes obtained using rotations),
 |
 |
Further, we also observe that this result is valid not only with five points, but with six. And if we go further, e.g. with nine points, we have the following
 |
 |
So actually, it is possible to have a simple conjecture: let denote the number of points, and let
be so that
then the minimal distance is . Which can be visualized below,
Based on that pigeonhole theorem, I have the intuition that this result is valid (here we just count the number of inside-triangles), but can we check if this is correct or not ? One idea might be to draw points randomly, and so see where points might end, or at least the maximal distance obtained over millions of random draws… But standard monte carlo might take a while… so we can use two ideas. One idea is from quasi-monte carlo techniques: since we want points to be to be as separate as possible, we do not need to draw randomly in the triangle, but perhaps we can draw randomly points on some grid. A second idea from the latin hypercube technique (and that pigeonholes theorem): instead to generating points randomly in the triangle, perhaps we can draw them in specific regions. For instance, with five points, we know that 4 points have to be in those sub-triangles, and the additional point in any one of those triangle. And because of symmetry, with five points, we can claim that this additional point has to be specifically in one sub-triangle. A random sample with five points within the same sub-triangle will be useless (and a waste of computational time).
With the following code, we define a grid, for a triangle, either upward, or downward, starting from some point (on the left), with a given length, and a given number of subdivision.
TRIANGLES=function(xinf,yinf,l,n,updown="up"){ X=NA;Y=NA for(i in n:1){ u=xinf+seq(0+(n-i)/2/(n-1),1-(n-i)/2/(n-1),length=i)*l if(updown=="up") v=rep(yinf+(n-i)*sqrt(3)/2*l/(n-1),i) if(updown=="down") v=rep(yinf-(n-i)*sqrt(3)/2*l/(n-1),i) X=c(X,u);Y=c(Y,v)} return(cbind(X[-1],Y[-1]))}
 |
 |
Here are grid with respectively 20 and 50 points on the lower side. It is then possible to define 4 grids, corresponding to the four sub-triangles,
k=3;st=6 firstgrid=TRIANGLES(0,0,(k-2)/(k-1),k-1) secondgrid1=TRIANGLES( firstgrid[1,1],firstgrid[1,2],1*(k-2)/(k-1),st) secondgrid2=TRIANGLES( firstgrid[2,1],firstgrid[1,2],1*(k-2)/(k-1),st) secondgrid3=TRIANGLES( firstgrid[3,1],firstgrid[3,2],1*(k-2)/(k-1),st) secondgrid4=TRIANGLES( firstgrid[3,1],firstgrid[3,2],1*(k-2)/(k-1),st,updown="down")
Then, we just draw randomly five points on that grid, in the four sub-triangles,
N=5 Dmax=0 setpointmax=matrix(0,4,2) indice=c(1:4,sample(1:4,size=N-4,replace=FALSE)) tindice=table(indice) indice1=sample(1:nrow(secondgrid1),size= tindice[1],replace=FALSE) indice2=sample(1:nrow(secondgrid2),size= tindice[2],replace=FALSE) indice3=sample(1:nrow(secondgrid3),size= tindice[3],replace=FALSE) indice4=sample(1:nrow(secondgrid4),size= tindice[4],replace=FALSE) setpoint=rbind(secondgrid1[indice1,],secondgrid2[indice2,], secondgrid3[indice3,],secondgrid4[indice4,])
No, we can run a code, where we keep in mind locations of the five points each time we beak a record,
D=min(dist(setpoint,"euclidean")) if(D>Dmax){Dmax=D setpointmax=setpoint}
Here are some locations obtained after running the algorithm a few times (with five points)
On the graph below, we can visualize the time it takes before having a record, and the convergence towards 1/2 (which is the true value of the maximal distance)
The convergence is slow… extremely slow… However, we can run the same code for more than five points, e.g. seven points (actually, here sub-triangles are not used here, and it looks like we have been lucky here, since the convergence was rather fast),
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,
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…
ACT2040: examen final
L’examen final du cours ACT2040 aura lieu le 13 décembre, durant trois heures. L’examen sera en deux parties: une de tarification, et une de provisionnement. Les deux porteront sur la lecture et l’analyse critique de sorties informatiques. Pour la tarification, les bases de données utilisées sont en ligne ci-dessous
> BASEN=read.table("http://freakonometrics.free.fr/baseN.txt",header=TRUE,sep=";") > BASEY=read.table("http://freakonometrics.free.fr/baseY.txt",header=TRUE,sep=";") > head(BASEN) ageconducteur agepermis sexeconducteur situationfamiliale habitation zone 1 57 39 F Celiba peri-urbain 8 2 54 35 H Celiba urbain 3 3 51 32 F Celiba urbain 1 4 53 35 H Marie rural 4 5 61 43 H Marie urbain 8 6 60 29 F Marie peri-urbain 1 agevehicule proprietaire payment marque poids usage 1 12 locataire Annuel AUTRES 8.>3500kg PROMENADE 2 20 sans mrp Semestriel PEUGEOT 4.3100-3199kg PROMENADE 3 4 sans mrp Annuel RAPIDO 1.<2700kg PROMENADE 4 1 sans mrp Annuel AUTRES 3.3000-3099kg PROMENADE 5 1 proprietaire Annuel FIAT 6.3300-3399kg PROMENADE 6 10 sans mrp Mensuel FIAT 8.>3500kg PROMENADE exposition nombre voiture 1 1 0 Monospace 2 1 0 Berline 3 1 0 sans avp 4 1 0 sans avp 5 1 1 Monospace 6 1 0 sans avp
Parmi les variables, la description (sommaire) est la suivante,
- ageconducteur: âge du conducteur principal du véhicule
- agepermis: ancienneté du permis de conduire du conducteur principal du véhicule
- sexeconducteur: sexe du conducteur principal (H ou F)
- situationfamiliale: situation familiale du conducteur principal (“Celiba”, “Marie” ou “Veuf/Div”)
- habitation: zone d’habitation du conducteur principal (“peri-urbain”, “rural” ou “urbain” )
- zone: zone d’habitation (allant de 1 à 8)
- agevehicule: age du véhicule
- proprietaire: si le conducteur principal possède un contrat Habitation, son statut (“locataire” ou “proprietaire”) Sinon “sans mrp”
- payment:type de fractionnement de la prime d’assurance automobile (“Annuel”, “Mensuel” ou “Semestriel”)
- marque: marque du véhicule
> levels(BASEN[,10]) [1] "ADRIA" "AUTOSTAR" "AUTRES" "BURSTNER MOBIL" [5] "CHALLENGER" "CHAUSSON" "CITROEN" "FIAT" [9] "FORD" "HYMERMOBIL" "MERCEDES" "PEUGEOT" [13] "PILOTE" "RAPIDO" "RENAULT" "VOLKSWAGEN"
- poids: classe de poids du véhicule
> levels(BASEN[,11]) [1] "1.<2700kg" "2.2700-2999kg""3.3000-3099kg""4.3100-3199kg" [5] "5.3200-3299kg""6.3300-3399kg""7.3400-3499kg""8.>3500kg"
- usage: utilisation du véhicule principal (“PROMENADE” ou “TOUS_DEPLACEMENTS”)
- exposition: exposition, en années
- nombre: nombre d’accident responsabilité civile du conducteur principal, pendant l’année passée
- cout: cout du sinistre
- voiture: type de véhicule
> levels(BASEN[,15]) [1] "Berline" "Break" "Buggy" [4] "Cabriolet" "Combispace" "Coup\xe9" [7] "Coup\xe9 Cabriolet" "Jeep" "Minibus" [10] "Minispace" "Monospace" "sans avp"
Concernant la partie sur le provisionnement, des sorties associées au triangle suivant seront proposées,
> TRIANGLE=read.table("http://freakonometrics.blog.free.fr/public/code/triangleACT2040.txt") > TRIANGLE X1 X2 X3 X4 X5 X6 X7 X8 X9 X10 1 357848 1124788 1735330 2218270 2745596 3319994 3466336 3606286 3833515 3901463 2 352118 1236139 2170033 3353322 3799067 4120063 4647867 4914039 5339085 NA 3 290507 1292306 2218525 3235179 3985995 4132918 4628910 4909315 NA NA 4 310608 1418858 2195047 3757447 4029929 4381982 4588268 NA NA NA 5 443160 1136350 2128333 2897821 3402672 3873311 NA NA NA NA 6 396132 1333217 2180715 2985752 3691712 NA NA NA NA NA 7 440832 1288463 2419861 3483130 NA NA NA NA NA NA 8 359480 1421128 2864498 NA NA NA NA NA NA NA 9 376686 1363294 NA NA NA NA NA NA NA NA 10 344014 NA NA NA NA NA NA NA NA NA
Uncertainty in claims reserving (WIM)
I will be talking on Friday, at the 1st Québec-Ontario Workshop on Insurance Mathematics (so called WIN, already mentioned here). The program is now online here. My talk will be on Solvency II requirements in claims reserving (slides can be found here). Even if Canada will not adapt the Solvency II European capital test, it looks like Solvency II will matter to Canadian insurance companies (as mentioned here).
Discussions géométriques du soir
J’avais raconté l’autre jour l’aventure de mon fils qui souhaitait faire des triangles avec des bouts de bois (ici). Il a récidivé ce soir. Sauf qu’il n’avait qu’un grand baton, et que j’avais fermé le portillon qui l’empêchait d’aller dans la forêt en chercher (la nuance entre romancer et devenir mythomane doit être mince). Je lui ai suggéré de le casser en trois portions, puis de constituer un triangle avec les trois morceaux restant. Manque de bol, l’expérience rate (comme l’autre fois). Et me souvenant de mes pénibles calculs de la dernière fois, je me suis rappelé qu’il avait une chance sur deux de pouvoir effectivement construire un triangle… sauf que la manière d’obtenir les trois bouts de bois me semble différente… donc peut être que ce 1/2 n’est plus valide ?
En fait, je crois qu’il faut définir clairement comment on va couper les bouts de bois.
- méthode 1: le coupage simultané
Je suppose que je peux couper d’un coup le bout de bois en deux. Formellement, je tire deux points sur l’intevalle correspondant à mon bout de bois, les tirages étant uniformes, et indépendants. On arrive alors à trois morceaux, 
, 
et 
, en sachant que
la longueur du bout de bois dont je diposais initialement. On supposera que 
=1 pour simplifier (ce qui ne changera rien au résultat final, on s’en doute). L’ensemble des points que je peux obtenir peuvent alors se visualiser comme des points du simplexe,
autrement dit, on peut les représenter dans le triangle équilatéral ci-dessous.
Or on ne peut constituer un triangle que si la plus grande des longeurs n’excède pas la somme des deux autres. Supposons que soit la plus grande longueur (et qu’elle soit visualisée graphiquement par rapport à la distance au côté en bas), le point rouge ci-dessous (correpondant à un triplet de longueurs possibles) ne convient pas. En fait, dès lors que l’on entre dans l’aire rosée, on ne peut plus constituer un triangle…
Moralité, à moins de tomber dans le triangle jaune, on ne peut pas constituer de triangle…. Ce qui veut dire qu’on n’a plus qu’une chance sur quatre !
- méthode 2: le coupage séquentiel
Une autre méthodde consister à couper d’abord en deux (en tirant un point au hasard sur mon intervalle), puis de recouper en deux le plus grand morceau… Damned, mon fils a de bonne idée, mais ça sent les calculs à n’en plus finir ça…Ou alors va falloir continuer à essayer de faire des dessins, car mine de rien, c’est plus simple que les calculs de la dernière fois. Bon, la longueur de mon bout de bois est la hauteur de mon triangle. Je commence par le couper en deux, la distance au côté en bas devant être la plus petite des deux.
- couper d’un coup comme découpage séquentiel aléatoire ?
Mais attends deux minutes mon bonhomme (oui, je parle à mon fils). Si j’avais choisi de découper le plus petit morceau, on sait qu’il est impossible de constituer un triangle car par définition 
Supposons que l’on casse le baton en deux, et que l’on tire ensuite au hasard (avec probabilité 1/2 et 1/2) quelle partie sera alors coupée en deux, on doit se ramener au cas précédant non ? La probabilité de constituer un triangle doit alors être
Damned, mais ça ne fait pas un quart ça ! Bon, je crois que je vais prendre une aspirine et aller demander à mon fils d’aller ranger ses bâtons !