Mercredi dernier avait lieu l’examen final d’ACT2040. La sortie informatique utilisée (avec quelques valeurs manquantes qu’il fallait retrouver, et qui comprenait des indications numérotées pour des questions spécifiques) et le questionnaire sont en ligne, ainsi que des éléments de réponse. Comme évoqué dans la page de garde de la ‘correction’, le barème sera proche de celui annoncé, mais un coefficient multiplicatif devrait être appliqué car je me suis rendu compte, en tapant la correction que le devoir était long. Je vais profiter des jours à venir pour corriger l’examen, avant d’attaquer la correction des devoirs, pendant les vacances. Toutes les remarques sur la ‘correction’ sont les bienvenues. Sinon, bonnes vacances…
Tag Archives: ACT2040
Examen ACT2040
L’examen pour le cours ACT2040 aura lieu mercredi matin. Comme promis, les sorties informatiques qui seront utilisées sont maintenant en ligne [une coquille s’était glissée dans les dernières pages, une version corrigée a été mise ne ligne mardi a 15:10, a la place de la première version]. Comme je l’avais dit mercredi dernier, la régression binomiale sera utilisée, et la dernière partie portera sur de la construction de tables de mortalité prospective. Des détails sont évoqués dans Actuariat avec R. Sinon, pour ceux qui le souhaitent, les examens des années passées sont en ligne,
- final, décembre 2011, un énoncé basé sur quelques sorties informatiques, avec des éléments de correction
- intra, février 2013, énoncé de l’examen intra en pdf a partir des sorties ici, avec des éléments de correction en pdf (et un petit complément)
- final, avril 2013, un énoncé est en ligne ainsi que des éléments de correction.
Les calculatrices sont fortement recommandées (un téléphone, aussi intelligent soit-il, n’est pas une calculatrice).
Données, examen final
L’examen final du cours ACT2040 aura lieu, comme confirmé par courriel il y a plusieurs jours, le mercredi 11 décembre, pendant 3 heures. Il y aura surtout des sorties à commenter. Les sorties seront obtenues à partir des trois jeux de données suivants
> CORPOREL=read.table("http://freakonometrics.free.fr/corporel-2040.csv",header=TRUE,sep=";") > tail(CORPOREL) degre age cat.age sexe vehicule anciennete alcool cat.alc 76336 indemne 45 40-49 M voiture 6 0 0-20 76337 corporel 59 50-59 F voiture 2 0 0-20 76338 indemne 34 30-39 F voiture 2 0 0-20 76339 indemne 29 26-29 F voiture 5 9 0-20 76340 indemne 64 60+ M voiture 0 0 0-20 76341 indemne 57 50-59 F voiture 1 0 0-20
Il s’agit d’observations d’accident automobiles, en Australie, la variable d’intérêt étant ici la gravité (le degré) de l’accident. Parmi les variables explicatives, l’âge du conducteur (en variable ‘continue‘ et en classes arbitraires) et son degré d’alcoolémie (en variable ‘continue‘ et en classes liées à des critères légaux, en g/10L), l’âge de la voiture, et le sexe du conducteur (M pour les hommes et F pour les femmes).
Le second jeu de données est un triangle de paiements, cumulés, avec la prime acquise, par année calendaire,
> load("http://freakonometrics.free.fr/triangle-intra.txt") > intra $triangle 0 1 2 3 4 5 6 7 8 9 1988 5244 9228 10823 11352 11791 12082 12120 12199 12215 12215 1989 5984 9939 11725 12346 12746 12909 13034 13109 13113 NA 1990 7452 12421 14171 14752 15066 15354 15637 15720 NA NA 1991 7115 11117 12488 13274 13662 13859 13872 NA NA NA 1992 5753 8969 9917 10697 11135 11282 NA NA NA NA 1993 3937 6524 7989 8543 8757 NA NA NA NA NA 1994 5127 8212 8976 9325 NA NA NA NA NA NA 1995 5046 8006 8984 NA NA NA NA NA NA NA 1996 5129 8202 NA NA NA NA NA NA NA NA 1997 3689 NA NA NA NA NA NA NA NA NA $primes [1] 15883 16689 18029 17858 16709 14212 15083 15131 15465 11217
en cas de soucis pour l’importation, il sera possible d’utiliser la commande
> source("http://freakonometrics.free.fr/triangle-intra2.R")
Enfin, le dernier jeux de données sera un peu plus insolite pour ce cours. Il s’agit de données de tables de mortalité, au Canada
> DECES=read.table("http://freakonometrics.free.fr/DECES-CAN.csv",header=TRUE,sep=";") > tail(DECES) D E A Y 772 84 147 105 2010 773 39 76 106 2010 774 23 40 107 2010 775 15 20 108 2010 776 7 8 109 2010 777 5 5 110 2010
Par année Y, on observe un certain nombre de personnes en vie E, d’âge A au 1er janvier, et un certain nombre va décéder dans l’année, D. Par exemple, en 2010, il y avait 147 personnes de 105 ans en début d’année, 84 vont décéder dans l’année. L’idée sera d’utiliser les modèles de provisionnement pour modéliser les probabilités de décéder.
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/…
Pricing reinsurance contracts, another case study
A reinsurance case study for tomorrow’s class. The goal will be to price some nonproportional reinsurance contract, for business interruption claims. Consider the following dataset,
> library(gdata) > db=read.xls( + "https://perso.univ-rennes1.fr/arthur.charpentier/SIN_1985_2000-PE.xls", + sheet=1) Content type 'application/vnd.ms-excel' length 183808 bytes (179 Kb) open URL ================================================== downloaded 179 Kb
As for any (standard) insurance contract, there are two parts in the pricing
- the expected number of claims
- the average cost of individual claims
Here, we do not have covariates (but it might be possible to use some, like the kind of industry, the location, etc).
Let us start with the expected number of claims, per year. Here is the daily frequency,
The data are rather old… but somehow, it is a good thing since after ten years, we can expect that most of the claims have been settled (we’ll discuss claims dynamic starting next week). To plot the graph above, we use
> date=db$DSUR > D=as.Date(as.character(date),format="%Y%m%d") > vD=seq(min(D),max(D),by=1) > sD=table(D) > d1=as.Date(names(sD)) > d2=vD[-which(vD%in%d1)] > vecteur.date=c(d1,d2) > vecteur.cpte=c(as.numeric(sD),rep(0,length(d2))) > base=data.frame(date=vecteur.date,cpte=vecteur.cpte) > plot(vecteur.date,vecteur.cpte,type="h",xlim=as.Date(as.character( + c(19850101,20111231)),format="%Y%m%d"))
Then, we can get a prediction of the daily number of business interruption claims, e.g. for any day in 2010 (assume that we had to price a reinsurance contract a few years ago), using a (standard) Poisson regression
> regdate=glm(cpte~date,data=base,family=poisson(link="log")) > nd2010=data.frame(date=seq(as.Date(as.character(20100101),format="%Y%m%d"), + as.Date(as.character(20101231),format="%Y%m%d"),by=1)) > pred2010 =predict(regdate,newdata=nd2010,type="response") > sum(pred2010) [1] 159.4757
Observe that using old data has drawbacks, since we got much more uncertainty if we use a regression on time (to include some possible trend)
Say we have something like 160 claims over a given year, on average.
> plot(D,db$COUTSIN,type="h")
Let us now focus on the cost of those claims. We have 2,400 claims in our dataset, to fit a model (or at least estimate how much a reinsurance contract might cost us). Assume that we would like to purchase a reinsurance contract for our very large claims. Like the two largest per year. Over 16 years, the decutible should be close to the cost of the 32nd largest claim, which was close to 15 million.
> quantile(db$COUTSIN,1-32/2400)/1e6 98.66667% 15.34579 > abline(h=quantile(db$COUTSIN,1-32/2400),col="blue")
So consider some reinsurance contract with a deductible of 15 million. Unfortunately, we cannot find unlimited covers. So let us assume that a reinsurance company agrees for such a deductible, but with a limited cover of 35 million. The average cost (for the reinsurance company) is where
A first idea is to look at the first cost, i.e. the empirical average of that indemnity, on our portfolio. The indemnity function is
> indemn=function(x) pmin((x-15)*(x>15),50-15)
we can check on a few losses that it is actually what we wish to compute
> indemn(5) [1] 0 > indemn(20) [1] 5 > indemn(50) [1] 35
Now, if the compute the average repayment by the reinsurance company, over 16 years, we get
> mean(indemn(db$COUTSIN/1e6)) [1] 0.1624292
So, per claim, the reinsurance company will pay, on average 162,430. With 160 claims per year, the pure premium should be close to 26 million
> mean(indemn(db$COUTSIN/1e6))*160 [1] 25.98867
(again, for a 35 million cover, for some claims that should occur, on average, twice a year). As we will see, a standard model in reinsurance is the Pareto distribution (or to be more specific, a Generalized Pareto one),
There are three parameters here
- the threshold
(that we will consider as fixed, but we will see its impact on reinsurance pricing)
- the scale parameter
(called
in R)
- the tail index
The strategy is to consider a threshold below our deductible, e.g. 12 million. Then, given that the loss exceed 12 million, we can fit a Generalized Pareto distribution,
> gpd.PL <- gpd(db$COUTSIN,12e6)$par.ests > gpd.PL xi beta 7.004147e-01 4.400115e+06
and compute
> E <- function(yinf,ysup,xi,beta,threshold){ + as.numeric(integrate(function(x) (x-yinf)*dgpd(x,xi,mu=threshold,beta), + lower=yinf,upper=ysup)$value+ + (1-pgpd(ysup,xi,mu=threshold,beta))*(ysup-yinf)) + }
Here, given that a claim exceeds 12 million, the average repayment is close to 6 million
> E(15e6,50e6,gpd.PL[1],gpd.PL[2],12e6) [1] 6058125
Now, we have to take into account the probability to reach 12 million, which is here
> mean(db$COUTSIN>12e6) [1] 0.02639296
So, if we summarize, we have on average 160 claims per year,
> p [1] 159.4757
Only 2.6% will exceed 12 million
> mean(db$COUTSIN>12e6) [1] 0.02639296
So, the yearly frequency of claism larger than 12 million is 4.2 claims
> p*mean(db$COUTSIN>12e6) [1] 4.209036
And for a claim that exceed 12 million, the average repayment is
> E(15e6,50e6,gpd.PL[1],gpd.PL[2],12e6) [1] 6058125
So, the pure premium should be close to
> p*mean(db$COUTSIN>12e6)*E(15e6,50e6,gpd.PL[1],gpd.PL[2],12e6) [1] 25498867
which (hopefully) is close to the empirical value we got. Actually, it is also possible to look at the impact of the threshold parameter, since it is clearly and intermediate value that could be changed. I mean, why 12 and not 10? Consider
> esp=function(threshold=12e6,p=sum(pred2010)){ + (gpd.PL <- gpd(db$COUTSIN,threshold)$par.ests) + return(p*mean(db$COUTSIN>threshold)*E(15e6,50e6,gpd.PL[1],gpd.PL[2],threshold)) + }
We can plot the pure premium as a function of that threshold,
> seuils=seq(1e6,15e6,by=1e6) > plot(seuils,Vectorize(esp)(seuils),type="b",col="red")
which is between 24 and 26 for large thresholds. Again, that is only the first step, and we can price a higher reinsurance layer, like a reinsurance contract with a deductible of 50 million (we have our previous reinsurance contract for claims below that threshold), and a cover of 50 million, for instance. For those high layers, it become interesting to have a parametric model, which should be more robust than the empirical average.
Modélisation des coûts individuels
Cette semaine, même si le réseau de l’UQAM est down, on va continuer le cours et finir la section sur la modélisation de la surdispersion pour la fréquence de sinistres. On devrait ensuite commencer la modélisation des coûts individuels. En particulier, on passera du temps autour de deux points,
- la distinction lognormale et gamma
- l’écrêtement des gros sinistres
Les transparents sont en ligne. Et la base des coûts est celle évoquée au second cours.
Please, never use my codes without checking twice (at least)!
I wanted to get back on some interesting experience, following a discussion I had with Carlos after my class, this morning. Let me simplify the problem, and change also the dataset. Consider the following dataset
> db = read.table("http://freakonometrics.free.fr/db2.txt",header=TRUE,sep=";")
Let me change also one little thing (in the course, we use the age of people as explanatory variables, so let us consider rounded figures to),
> db$X1=round(db$X1*10) > db$X2=round(db$X2*10)
Assume that you want to work with factors, because you don’t see why there should be some linear model (and you did not look at the awesome posts on smoothing techniques)
> db$X1F=cut(db$X1,c(-12,45,75,120)) > db$X2F=cut(db$X2,c(100,200,300))
Then you run your regression,
> reg = glm(Y~X1F+X2F+X3,family=binomial,data=db)
So far, nothing wrong, you can try, no error, no warning. Then Carlos wanted to use a ROC curve to see how the model was performing… so he did use some code I uploaded on the blog, something like
> reg = glm(Y~X1+X2+X3,family=binomial,data=db) > S = predict(reg,type="response") > Y = db$Y
> plot(0:1,0:1,xlab="False Positive Rate",ylab="True Positive Rate",cex=.5) > for(s in seq(0,1,by=.01)){ + Ps=(S>s)*1 + FP=sum((Ps==1)*(nombre==0))/sum(nombre==0) + TP=sum((Ps==1)*(nombre==1))/sum(nombre==1) + points(FP,TP,cex=.5,col="red") + }
To make it nicer, let us use the following code
> ROCcurve=function(s){ + Ps=(S>s)*1 + FP=sum((Ps==1)*(Y==0))/sum(Y==0) + TP=sum((Ps==1)*(Y==1))/sum(Y==1) + return(c(FP,TP))} > u=seq(0,1,by=.001) > vectROC=Vectorize(ROCcurve)(u)
Here, I should mention that I got a warning, but to be honest, when you see the graph, you’re so puzzled that you forget about it…
There were 50 or more warnings (use warnings() to see the first 50)
Carlos ran the code, show me the graph, and asked me “can my model be that bad?”,
My first answer is “no, your model cannot be that bad! you use (almost) all your explanatory variables, you cannot have such a bad model…“. So, where the problem comes from? My first guess is that this is what you get using some random sort of a variable. For instance, if you use the code above on
> reg2 = glm(Y~X1+X2+X3,family=binomial,data=db) > S = sort(predict(reg2,type="response")) > Y = db$Y
Here, one variable is sorted, and we get
This confirms the idea that using a random classifier, the ROC curve is on the diagonal. But here, when you look at the code, you do not see any sorting operation.
So, again, what went wrong? The problem is actually very simple. The dataset looks like
> head(db) Y X1 X2 X3 X1F X2F 1 1 33 163 B (-12,45] (100,200] 2 1 64 185 D (45,75] (100,200] 3 1 53 166 B (45,75] (100,200] 4 1 55 197 C (45,75] (100,200] 5 1 41 184 C (-12,45] (100,200] 6 1 78 196 C (75,120] (100,200]
If we look at the factor variable, the lowest part is not included, in the interval. More precisely,
> levels(db$X1F) [1] "(-12,45]" "(45,75]" "(75,120]"
and if we look more precisely and the range of the two variables, we get
> range(db$X1) [1] -12 120
Wait… the minimum is -12, but it does not appear in the factor (and we do not get any warning)? Yes,
> db[which.min(db$X1)+(-2):2,] Y X1 X2 X3 X1F X2F 429 1 76 227 A (75,120] (200,300] 430 0 35 186 D (-12,45] (100,200] 431 0 -12 109 E <NA> (100,200] 432 1 76 225 B (75,120] (200,300] 433 1 61 206 A (45,75] (200,300]
Now we start so see more clearly what went wrong… there is a missing value in the dataset. And when we get our prediction, there is no missing value: the missing value has be droped.
> length(predict(reg)) [1] 999 > nrow(db) [1] 1000
So when we compare the prediction and the observed value…. there is a problem, since the vectors don’t match. Actually, it was mentioned in the output of the regression (but we did not look at it)
> summary(reg) Call: glm(formula = Y ~ X1F + X2F + X3, family = binomial, data = db) Deviance Residuals: Min 1Q Median 3Q Max -2.9432 0.1627 0.2900 0.4823 1.1377 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.22696 0.23206 0.978 0.328067 X1F(45,75] 1.86569 0.23397 7.974 1.53e-15 *** X1F(75,120] 2.97463 0.65071 4.571 4.85e-06 *** X2F(200,300] 1.11643 0.32695 3.415 0.000639 *** X3B -0.06131 0.31076 -0.197 0.843609 X3C 0.75013 0.35825 2.094 0.036268 * X3D 0.13846 0.31399 0.441 0.659226 X3E -0.13277 0.31859 -0.417 0.676853 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 802.34 on 998 degrees of freedom Residual deviance: 607.61 on 991 degrees of freedom (1 observation deleted due to missingness) AIC: 623.61 Number of Fisher Scoring iterations: 7
Yes, there is a tiny little sentence,
(1 observation deleted due to missingness)
So that was it? Yes! If you drop that missing value, you get something more realistic,
> S= predict(reg,type="response") > Y=db$Y[-which.min(db$X1)]
or (probably better), change the left point of the interval
> db$X1F=cut(db$X1,c(-13,45,75,120)) > reg = glm(Y~X1F+X2F+X3,family=binomial,data=db) > S= predict(reg,type="response") > Y= db$Y
(the output would have been almost the same).
Observe that if we had used a dedicated package, we would not have encountered this problem. For instance (starting with the initial values of the vectors)
> library(ROCR) > prediction(S,Y) Error in prediction(S, Y) : Number of predictions in each run must be equal to the number of labels for each run.
We do have the answer here: both vectors do not have the same length…
So, what is my point? R is great, because of all the packages, but as a teacher, I do not feel comfortable asking my student to use those functions, as black boxes. So I try to write my own codes, to get the same output. So yes, I do write codes to explain what the black box is doing, to simplify the algorithm, and show what’s going on. When working on a (forthcoming) book as the Editor, we had a discussion with Rob Hyndman about that issue. I wanted the contributors to explain the core of the code, with a simplified algorithm, when using a dedicated package. I do truly believe that using simplified codes might help to understand better. Until you start to have problems. Because I write a code to deal with one specific problem, there is no check for possible errors. And once the code is understood, please, please do not use it! use R function that can handle errors…
Introduction aux GLM
Cette semaine, on finit la régression de Poisson (temporairement) avant de présenter la théorie des GLM. Les transparents sont en ligne. On en aura besoin pour aller plus loin sur les modèles avec surdispersion, pour modéliser la fréquence de sinistre, mais aussi pour modéliser les coûts.
Régression de Poisson
Mercredi, on finira les arbres de classification, et on commencera la modélisation de la fréquence de sinistre. Les transparents sont en ligne.
Comme annoncé lors du premier cours, je suggère de commencer la lecture du Practicionner’s Guide to Generalized Linear Models. Le document correspond au minimum attendu dans ce cours.
Introduction à la régression logistique et aux arbres
Pour le second cours ACT2040, on va finir l’introduction (et les rappels de statistique inférentiel) puis attaque la première grosse section, sur la régression logistique et aux arbres de classification. base tirée du livre de Jed Frees, http://instruction.bus.wisc.edu/jfrees/…
> baseavocat=read.table("http://freakonometrics.free.fr/AutoBI.csv", + header=TRUE,sep=",") > tail(baseavocat) CASENUM ATTORNEY CLMSEX MARITAL CLMINSUR SEATBELT CLMAGE LOSS 1335 34204 2 2 2 2 1 26 0.161 1336 34210 2 1 2 2 1 NA 0.576 1337 34220 1 2 1 2 1 46 3.705 1338 34223 2 2 1 2 1 39 0.099 1339 34245 1 2 2 1 1 18 3.277 1340 34253 2 2 2 2 1 30 0.688
On dispose d’une variable dichotomique indiquant si un assuré – suite à un accident de la route – a été représenté par un avocat (1 si oui, 2 si non). On connaît le sexe de l’assuré (1 pour les hommes et 2 pour les femmes), le statut marital (1 s’il est marié, 2 s’il est célibataire, 3 pour un veuf, et 4 pour un assuré divorcé). On sait aussi si l’assuré portait ou non une ceinture de sécurité lorsque l’accident s’est produit (1 si oui, 2 si non et 3 si l’information n’est pas connue). Enfin, une information pour savoir si le conducteur du véhicule était ou non assuré (1 si oui, 2 si non et 3 si l’information n’est pas connue). On va recoder un peu les données afin de les rendre plus claires à lire.
Les transparents sont en ligne sur le blog,
Des compléments théoriques sur les arbres peuvent se trouver en ligne http://genome.jouy.inra.fr/…, http://ensmp.fr/…, ou http://ujf-grenoble.fr/… (pour information, nous ne verrons que la méthode CART). Je peux renvoyer au livre (et au blog) de Stéphane Tuffery, ou (en anglais) au livre de Richard Berk, dont un résumé se trouve en ligne sur http://crim.upenn.edu/….
Assurance IARD
Mercredi aura lieu le premier cours d’assurance IARD de la session (ACT2040). Le plan de cours est en ligne, ainsi que les transparents de la première session. Nous verrons une introduction a la tarification, avec hétérogénéité, et quelques rappels de statistiques.
Des informations suivront, en particulier pour les démonstrations.
R tutorials
My course on non-life insurance (ACT2040) will start in a few weeks. I will use R to illustrate predictive modeling. A nice introduction for those who do not know R can be found online.
Provisionnement et tarification, examen final
L’examen final du cours ACT2040 avait lieu ce matin. L’énoncé est en ligne ainsi que des éléments de correction. En cas d’erreurs, merci de me le faire savoir rapidement, avant que je ne saisisse les notes.
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