The blog will not be updated during holidays, so Merry Christmas to all my readers (and all those who end up here by accident). And yes, I will spend some time on the Rivera Maya,
Monthly Archives: December 2011
Quel lancer de dé faire à Serpents et Échelles ?
Hier soir
, j’évoquais l’utilisation des chaînes de Markov au jeu Serpents et Échelles. Comme me le faisait remarquer Jean-Philippe, de manière assez étrange, les enfants sont toujours beaucoup plus content après avoir fait un 6 qu’après avoir fait un 1. Mais est-ce réellement la valeur optimale que le dé doit prendre ? Ça dépend bien sûr de la position. Par exemple, au premier lancer, on peut se demander ce que deviendrait le nombre de tours moyen nécessaires pour finir le jeu, conditionnellement aux 6 lancers possibles. On peut calculer, toujours conditionnellement aux 6 valeurs possibles de dés (et des positions où on va se retrouver), l’espérance du nombre de tours pour finir la partie,
esperance=function(h0){ INITIAL = as.numeric(which(M[h0+1,]>0))-1 ESPERANCE=rep(NA,length(INITIAL)) names(ESPERANCE)=INITIAL for(k in 1:length(INITIAL)){ initial=rep(0,n+1); initial[INITIAL[k]]=1 distrib=initial%*%M game=rep(NA,1000) for(h in 1:length(game)){ game[h]=distrib[n+1] distrib=distrib%*%M} ESPERANCE[k]=sum(1-game)} return(ESPERANCE)}
(à partir du code mis en ligne hier), e.g. pour le premier lancer,
> esperance(0) 1 2 3 5 6 14 32.16499 31.99947 31.82348 31.47954 31.36441 29.83020
où la valeur indiquée est la position où l’on pourrait se retrouver (la dernière correspond à un 4). On note que le meilleur premier lancer possible est celui qui nous amène sur la première échelle, i.e. la 4ème case,
Si on regarde case par case, on a les valeurs des “meilleurs” lancers de dés suivant (sans forcément l’unicité) en fonction de sa position sur le plateau
(avec en rouge les serpents et en bleu les échelles). Amusant non ?
Maps with R, and polygon boundaries
With R, it is extremely easy to draw maps. Let us start with something simple, like French regions. Baptiste mentioned on his blog that shapefiles can be downloaded from http://ign.fr/ website. Hence, if you extract the zip file, it is possible to get claims frequency per region (as done in the course ACT2040),
> library(maptools) > library(maps) > departements<-readShapeSpatial("DEPARTEMENT.SHP") > region<-tapply(baseFREQ[,"nbre"], + as.factor(baseFREQ[,"region"]),sum)/ + tapply(baseFREQ[,"exposition"], + as.factor(baseFREQ[,"region"]),sum) > depFREQ=rep(NA,nrow(departements)) > names(depFREQ)=as.character( + departements$CODE_REG) > for(nom in names(region)){ + depFREQ[names(depFREQ)==nom] = + region[nom]} > plot(departements,col=gray((depFREQ-.05)*20)) > legend(166963,6561753,legend=seq(1,0,by=-.1)/20+.05, + fill=gray(seq(1,0,by=-.1)),cex=1.25, bty="n")
Another application is on earthquakes. It is possible to use shapefiles of tectonic plates contour, and to relate earthquakes to plates. Shapefiles can be found onhttp://www.colorado.edu/ (here).
First, we can extract the shapes of the tectonic plates
> plates = readShapePoly("plates.shp",
+ proj4string=CRS("+proj=longlat"))
> PP=SpatialPolygons2PolySet(plates)
Consider Montreal,
> montreal=c(-73.600,45.500)
Given that specific location, it is possible to use the following code to get the associated plate,
> PLATE.loc=function(pt){ + K=NA + for(k in 1:17){ + c=point.in.polygon(pt[1], pt[2], + PP[PP$PID==k,c("X")],PP[PP$PID==k,c("Y")], + mode.checked=FALSE) + if(c>0){K=k} + } + return(K)} > abline(v=montreal[1],col="red") > abline(h=montreal[2],col="red") > PLATE.loc(montreal) [1] 1
and then to plot the associated tectonic plate very easily
> PLATE=function(k0){ + library(maps) + map("world") + polygon(PP[PP$PID==k0,c("X")],PP[PP$PID==k0,c("Y")], + col="red") + for(k in (1:17)[-k0]){polygon(PP[PP$PID==k,c("X")], + PP[PP$PID==k,c("Y")],col="light blue")} + map("world",add=TRUE)} > PLATE(PLATE.loc(montreal))
Those code were used in the paper written with Mathieu, and that will be presented on January 30th at the Geotop seminar.
Basics on Markov Chain (for parents)
Markov chains is a very interesting and powerful tool. Especially for parents. Because if you think about it quickly, most of the games our kids are playing at are Markovian. For instance, snakes and ladders…
It is extremely easy to write down the transition matrix, one just need to define all snakes and ladders. For the one above, we have,
n=100 M=matrix(0,n+1,n+1+6) rownames(M)=0:n colnames(M)=0:(n+6) for(i in 1:6){diag(M[,(i+1):(i+1+n)])=1/6} M[,n+1]=apply(M[,(n+1):(n+1+6)],1,sum) M=M[,1:(n+1)] starting=c(4,9,17,20,28,40,51,54,62, 64,63,71,93,95,92) ending =c(14,31,7,38,84,59,67,34,19, 60,81,91,73,75,78) for(i in 1:length(starting)){ v=M[,starting[i]+1] ind=which(v>0) M[ind,starting[i]+1]=0 M[ind,ending[i]+1]=M[ind,ending[i]+1]+v[ind]}
So, why is it important to have a Markov Chain ? Because, once you’ve noticed that you had a Markov Chain game, you can derive anything you want. For instance, you can get the distribution after some turns,
powermat=function(P,h){ Ph=P if(h>1){ for(k in 2:h){ Ph=Ph%*%P}} return(Ph)} initial=c(1,rep(0,n)) COLOR=rev(heat.colors(101)) u=1:sqrt(n) boxes=data.frame( index=1:n, ord=rep(u,each=sqrt(n)), abs=rep(c(u,rev(u)),sqrt(n)/2)) position=function(h=1){ D=initial%*%powermat(M,h) plot(0:10,0:10,col="white",axes=FALSE, xlab="",ylab="",main=paste("Position after",h,"turns")) segments(0:10,rep(0,11),0:10,rep(10,11)) segments(rep(0,11),0:10,rep(10,11),0:10) for(i in 1:n){ polygon(boxes$abs[i]-c(0,0,1,1), boxes$ord[i]-c(0,1,1,0), col=COLOR[min(1+trunc(500*D[i+1]),101)], border=NA)} text(boxes$abs-.5,boxes$ord-.5, boxes$index,cex=.7) segments(c(0,10),rep(0,2),c(0,10),rep(10,2)) segments(rep(0,2),c(0,10),rep(10,2),c(0,10))}
Here, we have the following (note that I assume that once 100 is reached, the game is over)
Assume for instance, that after 10 turns, your daughter accidentally drops her pawn out of the game. Here is the theoretical (unconditional) position of her pawn after 10 turns,
so, if she claims she was either on 58, 59 or 60, here are the theoretical probabilities to be in each cell after 10 turns,
> h=10 > (initial%*%powermat(M,h))[59:61]/ + sum((initial%*%powermat(M,h))[59:61]) [1] 0.1597003 0.5168209 0.3234788
i.e. it is more likely she was on 59 (60th cell of the vector since we start in 0). You can also look at the distribution of the number of turns (at first with only one player).
distrib=initial%*%M game=rep(NA,1000) for(h in 1:length(game)){ game[h]=distrib[n+1] distrib=distrib%*%M} plot(1-game[1:200],type="l",lwd=2,col="red", ylab="Probability to be still playing")
Once you have that survival distribution, you have the expected number of turns to finish the game,
> sum(1-game) [1] 32.16499
i.e. in 33 turns, on average, your daughter reaches the 100 cell. But in 50% of the games, it takes less than 29,
But assuming that you are playing with your daughter, and that the game is over once one player reaches the 100 cell, it is possible to get the survival distribution of the first time one of us reaches the 100 cell,
plot((1-game[1:200])^2,type="l",lwd=2,col="blue", ylab="Probability to be still playing (2 players)")
Here, the expected number of turns before ending the game is
> sum((1-game)^2) [1] 23.40439
And if you ask your son to join the game, the survival distribution function is
plot((1-game[1:200])^3,type="l",lwd=2,col="purple", ylab="Probability to be still playing (3 players)")
i.e. the expected number of turns before the end is now
> sum((1-game)^3) [1] 20.02098
Combien de temps profite-t-on de ses grands parents ?
Ce Hier matin, je suis tombé un peu par hasard sur deux graphiques de l’INSEE (en France) avec l’age moyen des mères à l’accouchement, en fonction de rang de naissance de l’enfant, avec tout d’abord 1905-1965,
puis 1960-2000
Ces graphiques sont passionnants en soi – comme en ont témoigné pas mal de followers sur Twitter – mais ils m’ont fait m’interroger. En particulier, sur la croissance observée depuis 30 ans, qui me faisait penser à la tendance croissante observée sur les durées de vie. On n’a – malheureusement – pas accès au données complètes sur le site, mais on peut trouver d’autres donnéesintéressantes (en l’occurrence l’age moyen à la naissance).
> agenaissance=read.table("http://freakonometrics.blog.free.fr/ public/data/agenaissance.csv",header=TRUE,sep=",") > agenaissance$Age=as.character(agenaissance$AGE) > agenaissance$AGE=as.numeric(substr(agenaissance$Age,1,2))+ + as.numeric(substr(agenaissance$Age,4,4))/10 > plot(agenaissance$ANNEE+.5,agenaissance$AGE, + type="l",lwd=2,col="blue")
Visuellement, on retrouve la courbe en bleu foncée sur les graphiques ci-dessus,
On peut alors aller en cran plus loin, en se demandant non pas quel était l’âge moyen de la mère, mais de la grand-mère (au sens la mère de la mère)
> agenaissance$NAIS.MERE=(agenaissance$ANNEE+.5)- + agenaissance$AGE > w=(trunc(agenaissance$NAIS.MERE-.5)) > rownames(agenaissance)=agenaissance$ANNEE > a1=agenaissance[as.character(w),]$NAIS.MERE > a2=agenaissance[as.character(w+1),]$NAIS.MERE > p=agenaissance$NAIS.MERE-(w+.5) > agenaissance$NAIS.GRD.MERE=(1-p)*a1+p*a2 > agenaissance$age.GRD.MERE=agenaissance$ANNEE+.5- + agenaissance$NAIS.GRD.MERE > tail(agenaissance) ANNEE AGE Age NAIS.MERE NAIS.GRD.MERE age.GRD.MERE 2000 2000 30.3 30,3 1970.2 1942.87 57.63 2001 2001 30.4 30,4 1971.1 1943.80 57.70 2002 2002 30.4 30,4 1972.1 1944.92 57.58 2003 2003 30.5 30,5 1973.0 1945.95 57.55 2004 2004 30.5 30,5 1974.0 1947.05 57.45 2005 2005 30.6 30,6 1974.9 1948.04 57.46 > plot(agenaissance$ANNEE+.5,agenaissance$age.GRD.MERE, + type="l",lwd=2,col="red")
Là encore, on peut visualiser l’âge de la grand-mère maternelle à la naissance
A partir de là, on peut se demander combien de temps on profite de ses grands-parents (ou tout du moins ici de sa grand mère maternelle), en se basant sur les calculs d’espérance de vie résiduelle. En utilisant le modèle de Lee-Carter pour modéliser les taux de décès annuel, et en extrapolant sur le siècle en cours, on peut extrapoler les espérances de vie résiduelles.
> Deces <- read.table("http://freakonometrics.free.fr/ Deces-France.txt",header=TRUE) > Expo <- read.table("http://freakonometrics.free.fr/ Exposures-France.txt",header=TRUE,skip=2) > Deces$Age <- as.numeric(as.character(Deces$Age)) > Deces$Age[is.na(Deces$Age)] <- 110 > Expo$Age <- as.numeric(as.character(Expo$Age)) > Expo$Age[is.na(Expo$Age)] <- 110 > library(forecast) > library(demography) > YEAR <- unique(Deces$Year);nC=length(YEAR) > AGE <- unique(Deces$Age);nL=length(AGE) > MUF <- matrix(Deces$Female/Expo$Female,nL,nC) > POPF <- matrix(Expo$Female,nL,nC) > BASEF <- demogdata(data=MUF, pop=POPF,ages=AGE, + years=YEAR, type="mortality", + label="France", name="Femmes", lambda=1) > LCF <- lca(BASEF) > LCFf<-forecast(LCF,h=100) > A <- LCF$ax > B <- LCF$bx > K1 <- LCF$kt > K2 <- K1[length(K1)]+LCFf$kt.f$mean > K <- c(K1,K2) > MU <- matrix(NA,length(A),length(K)) > for(i in 1:length(A)){ + for(j in 1:length(K)){ + MU[i,j] <- exp(A[i]+B[i]*K[j]) }} > esp.vie = function(xentier,T){ + s <- seq(0,99-xentier-1) + MUd <- MU[xentier+1+s,T+s-1898] + Pxt <- cumprod(exp(-diag(MUd))) + ext <- sum(Pxt) + return(ext) } > EVIE = function(x,T){ + x1 <- trunc(x) + x2 <- x1+1 + return((1-(x-x1))*esp.vie(x1,T)+(x-x1)*esp.vie(x2,T)) } > agenaissance$EV=NA > for(i in 1:100){ + t <- 2006-i + agenaissance$EV[agenaissance$ANNEE==t]= + EVIE(x=agenaissance$age.GRD.MERE[ + agenaissance$ANNEE==t],t) } > tail(agenaissance) ANNEE AGE Age NAIS.MERE NAIS.GRD.MERE age.GRD.MERE EV 2000 30.3 30,3 1970.2 1942.87 57.63 29.13876 2001 30.4 30,4 1971.1 1943.80 57.70 29.17047 2002 30.4 30,4 1972.1 1944.92 57.58 29.39027 2003 30.5 30,5 1973.0 1945.95 57.55 29.52041 2004 30.5 30,5 1974.0 1947.05 57.45 29.72511 2005 30.6 30,6 1974.9 1948.04 57.46 29.80398
Autrement dit, sur la dernière ligne, l’espérance de vie (résiduelle) pour une femme de 57.46 ans en 2005 était d’environ 29.80 ans. On peut alors visualiser non seulement l’âge moyen de sa grand-mère à la naissance, mais son espérance de vie résiduelle,
On note que depuis 30 ans, en France, la durée (moyenne) pendant laquelle les petits-enfants vont profiter de leur grands parents s’est stabilisé à une trentaine d’années. On peut aussi continuer, et remonter d’un cran (en refaisant tourner le code avec quelques modifications): on a alors l’âge (moyen) de son arrière grand mère à la naissance,
et la durée de vie (résiduelle) des arrière grand mères
On manque ici de données, mais il semble que l’on profite – en moyenne – environ 5 ans de son arrière grand mère. Maintenant on peut aussi s’interroger sur les limites de cette étude rapide. En particulier, de même qu’il existe une corrélation forte entre les durées de vie de conjoints (e.g. broken heart syndrom de Jagger & Sutton (1991)), on peut se demander si la naissance d’enfants et de petits-enfants a un impact sur la durée de vie résiduelle d’une personne (ou si on peut supposer l’indépendance comme on l’a fait ici).
Le paradoxe de l’examen surprise…
La semaine dernière, Steve Landsburg reprenait (sur http://thebigquestions.com) sous une forme amusante un paradoxe que j’avais évoqué voilà deux ans sur le blog. Je ne résiste pas au plaisir de me replonger dans ce paradoxe…
Un professeur rédigeant son plan de cours veut mentionner un examen surprise, qui compterait pour 10% de la note finale. Il annonce donc à ses étudiants, au premier cours, qu’un examen surprise (et un seul) se tiendrait à une date (surprise) parmi les 12 séances de cours. Mais comme le fait noter un élève, l’examen ne pourrait avoir lieu lors de la dernière séance (la 12ème), car cela signifierait qu’il n’y a pas eu d’examen les 11 premières séances, et donc l’examen surprise devrait forcément se tenir au 12ème cours (mais alors il ne serait plus surprise du tout). Bref, impossible d’avoir un examen surprise ce jour là. Mais dans ce cas, surenchérit une étudiante, il ne peut non plus avoir lieu à l’avant dernière séance (la 11ème), car là encore cela signifierait qu’il n’y a pas eu d’examen les 10 premières séances, et donc l’examen surprise devrait forcément au 11ème cours (car il ne peut pas avoir lieu au 12ème comme on vient de le voir, et là encore, il ne serait plus surprise du tout). Bref, de manière récurrente, on peut montrer que l’examen surprise a alors lieu forcément le premier jour. Mais là encore, il n’y a pas de surprise. Et donc impossible qu’il y ait un examen surprise si ce dernier a été annoncé. Moralité, pour faire un examen surprise, il ne faut surtout pas évoquer qu’il y en aura un. Mais comme on ne peut pas faire un examen non annoncé dans les modalités de contrôle des connaissances, il est impossible de faire des surprises aux étudiants à l’université… C’est triste, non ?
Natural Catastrophe Insurance: How Should the Government Intervene?
An updated version of the joint paper with Benoit Le Maux is online on http://papers.ssrn.com/.
“The present paper develops a new theoretical framework for analyzing the decision to provide or buy insurance against the risk of natural catastrophes. In contrast with conventional models of insurance, the insurer has a non-zero probability of insolvency that depends on the distribution of the risks, the premium rate, and the amount of capital in the company. Among several results, we show that risk-averse policyholders will accept to pay higher rates for a government-provided insurance with unlimited guarantee. However, depending on the correlation between and within the regional risks, a government program can be more attractive to high-correlation than to low correlation areas, which may lead to inefficiencies if the insurance ratings are not appropriately chosen.”
ACT2040: examen final (suite et fin)
Comme annoncé, l’examen terminal a lieu ce matin, et il est basé sur des données mises en ligne sur le blog depuis une dizaine de jours. L’énoncé est un peu longue, avec une bonne trentaine de pages de sorties informatiques, et un énoncécomprenant une trentaine de questions. Des éléments de correction sont en ligne ci-dessous (plutôt que de taper une correction en pdf, autant mettre le code en ligne).
On suppose importée les bases (sinon, le code est ici). Pour la première partie – sur la tarification (a priori) – la fréquence annuelle de sinistres est tout simplement
> mean(BASEN$nombre)/mean(BASEN$exposition) [1] 0.05241692 > exp(glm(nombre~1,offset=log(exposition),data=BASEN, + family=poisson(link="log"))$coefficients) (Intercept) 0.05241692
(je vais essayer de donner à chaque question deux méthodes pour obtenir le résultat). Si on suppose que le nombre de sinistres par an suit une loi de Poisson, alors la probabilité d’avoir au moins un accident est
> lambda=mean(BASEN$nombre)/mean(BASEN$exposition) > 1-ppois(0,lambda) [1] 0.05106684 > 1-exp(-lambda) [1] 0.05106684
Si on cherche à distinguer par sexe du conducteur principal
> mean(BASEN$nombre[BASEN$sexeconducteur=="H"])/ + mean(BASEN$exposition[BASEN$sexeconducteur=="H"]) [1] 0.05218143 > mean(BASEN$nombre[BASEN$sexeconducteur=="F"])/ + mean(BASEN$exposition[BASEN$sexeconducteur=="F"]) [1] 0.05457013 > exp(glm(nombre~0+sexeconducteur,offset=log(exposition), + data=BASEN,family=poisson(link="log"))$coefficients) sexeconducteurF sexeconducteurH 0.05457013 0.05218143
Ces deux valeurs ne sont pas significativement différentes. Si on regarde la sortie de la régression
Coefficients:
Estimate Std. Error z value Pr(>|z|)
sexeconducteurF -2.90827 0.05464 -53.23 <2e-16 ***
sexeconducteurH -2.95303 0.01848 -159.82 <2e-16 ***
les deux coefficients sont quasiment identiques (en enlevant la constante). Si on rajoute la constante, la modalité de référence devientle sexe féminin
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.90827 0.05464 -53.230 <2e-16 ***
sexeconducteurH -0.04476 0.05768 -0.776 0.438
et la modalité homme n’est alors plus significative. Bref, le sexe du conducteur n’est pas une variable explicative.
Ensuite, parmi les variables qui semble pouvoir etre utilisée pour modéliser le nombre de sinistre, on pourra retenir l’age du conducteur (décroissant jusqu’à 30 ans, puis croissant), l’age du véhicule (décroissant), l’ancienneté du permis de conduire (décroissant les 10 premières années puis croissant – ou indépendant), éventuellement certains type de véhicules (le coupé cabriolet semble avoir – en moyenne – plus d’accident que la Jeep), la situation familiale, le lieu d’habitation (en particulier rural semble avoir moins d’accident qu’urbain), le statut d’habitation, le type de paiement de la prime (mensuel semblant plus risqué qu’annuel), le poids du véhicule, ou la marque (ou plutot certaines marques, comme Peugeot versus Burstner Mobil). Bref, toutes, sauf le sexe du conducteur.
Pour un conducteur de 45 ans qui a son permis de conduire depuis seulement 5 ans, les trois modèles donnent la fréquence annuelle suivante
> nbase=data.frame(ageconducteur=45,agepermis=5,exposition=1) > exp(t(c(1,5))%*%reg6$coefficients) [,1] [1,] 0.03395972 > predict(reg6,newdata=nbase,type="response") 1 0.03395972 > exp(t(c(1,45))%*%reg7$coefficients) [,1] [1,] 0.04100356 > predict(reg7,newdata=nbase,type="response") 1 0.04100356 > exp(t(c(1,5,45))%*%reg8$coefficients) [,1] [1,] 0.05395398 > predict(reg8,newdata=nbase,type="response") 1 0.05395398
Les trois modèles sont respectivement un modèle de Poison avec un lien log, un modèle binomiale négative avec un lien log, et un modèle binomial, avec un lien logit. Les deux premiers sont très proches: ils ont la meme fonction de lien, et comme nous l’avions vu à l’examen intra, le modèle binomiale négative est un modèle de Poisson avec une variable latente non observée qui suivrait une loi Gamma. Donc pour ces deux modèles, si les classes sont homogènes, et que par classe, la loi de Poisson est adaptée, il est normal d’avoir des coefficients très proches. Pour le modèle binomial avec un lien logistique, il faut se souvenir que la fonction lien est
car ici la probabilité d’avoir un accident est faible. Donc les deux fonctions liens sont proches. De plus, pour une loi de Poisson
Bref, pour notre loi de Poisson, on a presque une loi binomiale de paramètre 
si est petit.
Pour la prédiction suivant les trois modèles, on utilise classiquement
Comme on a accès aux sorties, pour les trois modèles, les trois prédictions sont très proches (mais on s’en doutait d’après la question précédante)
> predict(reg9,newdata=nbase,type="response") 1 0.04351037 > predict(reg10,newdata=nbase,type="response") 1 0.04351056 > predict(reg11,newdata=nbase,type="response") 1 0.04352741
Pour la méthode de biais minimal, il s’agit du principe de balancement, qui coïncide avec la régression de Poisson sur les classes. Pour les régressions stepwise, dans un cas, on élimine (étape par étape) des variables, et dans l’autre cas, des modalités de variables (pour les variables qualitatives).
Les deux modèles semblent valides, au sens où toutes les variables sont significatives. Le découpage en classes pour l’age du véhicule semble justifié par l’arbre de régression, qui suggèrait de découper à 8 et 14 ans. Pour les prédictions, on a pour le premier modèle
> nbase15=data.frame(ageconducteur=40,agepermis=20, + habitation="urbain",locataire=TRUE,paymensuel=FALSE, + marie=TRUEFORD=TRUE,RENAULT=FALSE,PEUGEOT=FALSE, + exposition=1) > predict(reg15,newdata=nbase15,type="response") 1 0.03398312
alors que le second donne
> nbase16=data.frame(ageconducteur=40,agevehF="NEUF",rural=FALSE, + locataire=TRUE,paymensuel=FALSE,marie=TRUE, + jeune=FALSE,ancien=FALSE,exposition=1) > predict(reg16,newdata=nbase16,type="response") 1 0.05208656
Moralité, nos deux modèles (qui semblent autant valides l’un que l’autre) donnent des prédictions relativement différentes. Si l’assuré change de voiture, cela n’a pas d’influence pour le second modèle (le véhicule est toujours dans la catégorie neuf). Pour le premier, il perd le rabais qu’il avait en conduisant un véhicule Ford, et donc la fréquence annuelle de sinistre augmente de 23.5%
> 1/exp(reg15$coefficients["FORDTRUE"]) FORDTRUE 1.235329
Bref, la prédiction devient
> predict(reg15,newdata=nbase15,type="response")/ + exp(reg15$coefficients["FORDTRUE"]) 1 0.04198033
On peut d’ailleurs le vérifier simplement,
> nbase15b=data.frame(ageconducteur=40,agepermis=20, + habitation="urbain",locataire=TRUE,paymensuel=FALSE, + marie=TRUE,FORD=FALSE,RENAULT=FALSE,PEUGEOT=FALSE, + exposition=1) > predict(reg15,newdata=nbase15b,type="response") 1 0.04198033
Si on passe maintenant à la modélisation de la charge de sinistres, le coût moyen par sinistre serait
> mean(BASEY$cout) [1] 1468.522
Cette fois, contrairement à la modélisation de la fréquence, il est plus difficile de trouver des variables explicatives. A priori, au vu des graphiques, aucune variable ne semblerait convenir…
Si on regarde la prédiction pour le premier modèle obtenu,
> nbase7=data.frame(VOLKSWAGEN=FALSE,agepermis=20, + sexeconducteur="H",usage="TOUS_DEPLACEMENTS") > predict(regr7,newdata=nbase7,type="response") 1 1388.622
alors que pour le second,
> exp(predict(regr8,newdata=nbase7))* + exp(summary(regr8)$sigma^2/2) 1 1359.552
Si c’est la femme de l’assuré qui avait pris l’assurance, cela n’aurait rien changé pour le second modèle. En revanche, pour le premier, l’assuré perdrait le rabais qu’il avait en étant un homme, et donc le cout moyen attendu augmenterait de 23%,
> 1/exp(regr7$coefficients["sexeconducteurH"])
sexeconducteurH
1.237832
i.e.
> predict(regr7,newdata=nbase7,type="response")/ + exp(regr7$coefficients["sexeconducteurH"]) 1 1718.88
La encore, on retrouve le meme résultat directement
> nbase7b=data.frame(VOLKSWAGEN=FALSE,agepermis=20, + sexeconducteur="F",usage="TOUS_DEPLACEMENTS") > predict(regr7,newdata=nbase7b,type="response") 1 1718.88
Passons maintenant aux calculs de la prime pure. Sans segmentation, la prime pure serait
Les prédictions des modèles pour la fréquence donnaient
> (n1=predict(reg15,newdata=nbase15,type="response")) 1 0.03398312 > (n2=predict(reg16,newdata=nbase16,type="response")) 1 0.05208656
alors que pour les couts moyens par sinistre, nous avions
> (c1=predict(regr7,newdata=nbase7,type="response")) 1 1388.622 > (c2=exp(predict(regr8,newdata=nbase7)) + *exp(summary(regr8)$sigma^2/2)) 1 1359.552
Les quatre primes pures obtenues sont alors
Passons maintenant à la seconde partie – sur le provisionnement des sinistres à payer. On commence par travailler sur le triangle
TRIANGLE2010
D1 D2 D3 D4 D5 D6 D7 D8 D9
A2002 357848 1124788 1735330 2218270 2745596 3319994 3466336 3606286 3833515
A2003 352118 1236139 2170033 3353322 3799067 4120063 4647867 4914039 5339085
A2004 290507 1292306 2218525 3235179 3985995 4132918 4628910 4909315 NA
A2005 310608 1418858 2195047 3757447 4029929 4381982 4588268 NA NA
A2006 443160 1136350 2128333 2897821 3402672 3873311 NA NA NA
A2007 396132 1333217 2180715 2985752 3691712 NA NA NA NA
A2008 440832 1288463 2419861 3483130 NA NA NA NA NA
A2009 359480 1421128 2864498 NA NA NA NA NA NA
A2010 376686 1363294 NA NA NA NA NA NA NA
Pour obtenir le montant payé en 2008 pour les sinistres survenus en 2004, on utilise la différence entre
> TRIANGLE2010["A2004",c("D4","D5")] D4 D5 A2004 3235179 3985995
i.e.
> TRIANGLE2010["A2004","D5"]-TRIANGLE2010["A2004","D4"] [1] 750816
Pour le montant payé au total en 2005, on utilise
> CUMULPAIEMENTS2010=TRIANGLE2010 > INCREMENTPAIEMENTS2010=CUMULPAIEMENTS2010 > INCREMENTPAIEMENTS2010[,2:9]=CUMULPAIEMENTS2010[,2:9]- + CUMULPAIEMENTS2010[,1:8] > sum(diag(as.matrix(INCREMENTPAIEMENTS2010[1:4,4:1]))) [1] 2729241
On attaque maintenant la méthode Chain Ladder. Les coefficients de développement sont
> lambda=rep(NA,8) > for(k in 1:8){ + lambda[k]=(sum(TRIANGLE2010[1:(9-k),k+1])/ + sum(TRIANGLE2010[1:(9-k),k]))} > lambda[1:5] [1] 3.474193 1.704149 1.460866 1.161765 1.095763
On peut alors obtenu ce qui sera payé en 2011 pour les sinistres survenus en 2010,
> TRIANGLE2010["A2010","D1"]*(lambda[1]-1) [1] 931993.8
alors que pour le montant payé en 2011 pour les sinistres survenus en 2009,
> TRIANGLE2010["A2009","D2"]*(lambda[2]-1) [1] 1000686
Pour le montant payé en 2010 et 2011 pour les sinistres survenus en 2010,
> TRIANGLE2010["A2010","D1"]*lambda[1]*lambda[2] [1] 2230186
Pour les sinistres survenus en 2008, la charge total estimée sera
> TRIANGLE2010["A2008","D3"]*prod(lambda[3:8]) [1] 5531128
ce qui fait un montant de provision
> TRIANGLE2010["A2008","D3"]*(prod(lambda[3:8])-1) [1] 3111267
On rajoute maintenant une diagonale dans notre triangle, i.e. une nouvelle année de paiements (en 2011).
Le montant payé en 2011 pour les sinistres survenus en 2010 est
> TRIANGLE2011["A2010","D2"]-TRIANGLE2011["A2010","D1"] [1] 986608
à comparer avec
> TRIANGLE2010["A2010","D1"]*(lambda[1]-1) [1] 931993.8
prédit par la méthode Chain-Ladder, soit une erreur de l’ordre de 5%,
> (TRIANGLE2011["A2010","D2"]- + TRIANGLE2011["A2010","D1"])/ + (TRIANGLE2010["A2010","D1"]*(lambda[1]-1))-1 [1] 0.05859927
Si on compare maintenant les coefficients de transition de la méthode Chain Ladder,
> lambda2011=rep(NA,9) > for(k in 2:9){ + lambda2011[k]=(sum(TRIANGLE2011[1:(10-k),k+1])/ + sum(TRIANGLE2011[1:(10-k),k]))} > lambda2011[1:5] [1] NA 1.747333 1.457413 1.173852 1.103824 > lambda[1:5] [1] 3.474193 1.704149 1.460866 1.161765 1.095763
Le coefficient permettant de passer de la colonne 2 à la colonne 3 augmente (légèrement). Si on regarde la prédiction du montant de paiements effectué en 2012 pour les sinistres survenus en 2009,
> TRIANGLE2011["A2009","D3"]*(lambda2011[3]-1) [1] 1310258
qui sont à comparer avec les
> TRIANGLE2010["A2009","D2"]*lambda[2]*(lambda[3]-1) [1] 1116132
estimés un an auparavant. La variation entre les estimations est de l’ordre de 17%,
> TRIANGLE2011["A2009","D3"]*(lambda2011[3]-1)/ + (TRIANGLE2010["A2009","D2"]*lambda[2]*(lambda[3]-1))-1 [1] 0.1739278
Entre 2010 et 2012, on peut estimer que l’on payera
> TRIANGLE2011["A2010","D2"]*lambda2011[2] [1] 2382128
pour les sinistres survenus en 2010.
Si on se focalise sur les sinistres survenus en 2008, la prédiction de la charge finale est ici
> TRIANGLE2011["A2008","D4"]*prod(lambda2011[4:9]) [1] 5660771
à comparer avec la prédiction faire en 2010,
> TRIANGLE2010["A2008","D3"]*prod(lambda[3:8]) [1] 5531128
i.e. notre nouvelle prédiction est 2% plus élevée,
> TRIANGLE2011["A2008","D4"]*prod(lambda2011[4:9])/ + (TRIANGLE2010["A2008","D3"]*prod(lambda[3:8]))-1 [1] 0.02343867
Le provisions pour sinistres à payer pour les sinistres survenus en 2008 est estimé, fin 2011, à
> TRIANGLE2011["A2008","D4"]*(prod(lambda2011[4:9])-1) [1] 2177641
alors que fin 2010, un montant de provision avait été constitué, à hauteur de
> TRIANGLE2010["A2008","D3"]*(prod(lambda[3:8])-1) [1] 3111267
On notera que la différence
> TRIANGLE2011["A2008","D4"]*(prod(lambda2011[4:9])-1)- + TRIANGLE2010["A2008","D3"]*(prod(lambda[3:8])-1) [1] -933626.7
vient des paiements effectués en 2011, et du changement dans l’estimation de charge ultime
> TRIANGLE2011["A2008","D4"]*prod(lambda2011[4:9])- + (TRIANGLE2010["A2008","D3"]*prod(lambda[3:8])) [1] 129642.3 > TRIANGLE2011["A2008","D4"]-TRIANGLE2011["A2008","D3"] [1] 1063269 > TRIANGLE2011["A2008","D4"]*prod(lambda2011[4:9])- + (TRIANGLE2010["A2008","D3"]*prod(lambda[3:8]))- + (TRIANGLE2011["A2008","D4"]-TRIANGLE2011["A2008","D3"]) [1] -933626.7
On attaque enfin la dernière partie, avec la régression de Poisson. La prédiction du montant de paiements faits en 2012 pour les sinistres survenus en 2009 est ici
> Y=as.vector(as.matrix(INCREMENTPAIEMENTS2011)) > ANNEE=rep(2002:2010,10) > DEVEL=rep(1:10,each=9) > baseTriangle=data.frame(Y,A=as.factor(ANNEE), + D=as.factor(DEVEL)) > reg=glm(Y~A+D,data=baseTriangle,famil=poisson(link="log")) > nbase=data.frame(A="2009",D="4") > predict(reg,newdata=nbase,type="response") 1 1310258
qui correspond à ce que donnait la méthode Chain Ladder
> TRIANGLE2011["A2009","D3"]*(lambda2011[3]-1) [1] 1310258
(mais on savait que les deux modèles prédisaient la même chose). Pour les sinistres survenus en 2005, le montant de provisions à constituer est
> nbase=data.frame(A="2005",D=c("8","9","10")) > predict(reg,newdata=nbase,type="response") 1 2 3 247190.0 370179.3 92268.5 > sum(predict(reg,newdata=nbase,type="response")) [1] 709637.8
qui, là encore, correspond à ce que prédisait la méthode Chain Ladder.
> TRIANGLE2011["A2005","D7"]*(prod(lambda2011[7:9])-1) [1] 709637.8
Quels résidus et quelle loi simuler ?
Suite à une question sur les résidus, je vais essayer de prendre deux minutes pour reprendre un point que je n’avais pas abordé en cours, faute de temps. Les résidus les plus intéressants en provisionnement (et dans les GLM) sont les résidus de Pearson,
Mais il existe également des résidus dits ajustés, permettant de prendre en compte le fait que le nombre de paramètres est ici relativement grand, ramené au nombre d’observations. On pose alors
A première vue, utiliser l’un ou l’autre devrait donner la même chose si l’on génère des pseudo-triangles, puisque dans un cas, on utiliserait
et dans l’autre
Bref, multiplier par une constante pour ensuite diviser par la même constante, ça revient au même. Sauf qu’il peut sembler légitime de poser
(on reviendra là dessus tout à l’heure). En attendant, voilà nos deux résidus
> (E=residuals(regp,"pearson")) 1 2 3 4 9.488238e-01 2.404895e-02 1.168421e-01 -1.082940e+00 5 6 7 8 1.302749e-01 -1.007348e-13 -1.128013e+00 2.773332e-01 9 10 11 13 5.669707e-02 8.919633e-01 -2.110748e-01 -1.533031e+00 14 15 16 19 -2.213449e+00 -1.024162e+00 4.237393e+00 -4.899687e-01 20 21 25 26 7.929194e-01 -2.972380e-01 -4.275912e-01 4.140426e-01 31 -6.202125e-15 > n=sum(is.na(Y)==FALSE) > k=ncol(PAID)+nrow(PAID)-1 > (R=residuals(regp,"pearson")*sqrt(n/(n-k))) 1 2 3 4 1.374976e+00 3.485024e-02 1.693203e-01 -1.569329e+00 5 6 7 8 1.887862e-01 -1.459787e-13 -1.634646e+00 4.018940e-01 9 10 11 13 8.216186e-02 1.292578e+00 -3.058764e-01 -2.221573e+00 14 15 16 19 -3.207593e+00 -1.484151e+00 6.140566e+00 -7.100321e-01 20 21 25 26 1.149049e+00 -4.307387e-01 -6.196386e-01 6.000048e-01 31 -8.987734e-15
L’autre point sur lequel je voulais revenir sur les simulations, est que la loi de Poisson n’est peut-être pas adaptée pour simuler des scénarios de paiements futurs. En effet, si on fait une régression quasipoisson, on voit que le paramètre de surdispersion n’est pas négligeable,
> regqp=glm(Y~as.factor(D)+as.factor(A), + data=base,family=quasipoisson(link="log")) > summary(regqp) Call: glm(formula = Y ~ as.factor(D) + as.factor(A), family = quasipoisson(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 t value Pr(>|t|) (Intercept) 8.05697 0.02769 290.995 < 2e-16 *** as.factor(D)2 -0.96513 0.02427 -39.772 2.41e-12 *** as.factor(D)3 -4.14853 0.11805 -35.142 8.26e-12 *** as.factor(D)4 -5.10499 0.22548 -22.641 6.36e-10 *** as.factor(D)5 -5.94962 0.43338 -13.728 8.17e-08 *** as.factor(D)6 -5.01244 0.39050 -12.836 1.55e-07 *** as.factor(A)2002 0.06440 0.03731 1.726 0.115054 as.factor(A)2003 0.20242 0.03615 5.599 0.000228 *** as.factor(A)2004 0.31175 0.03535 8.820 4.96e-06 *** as.factor(A)2005 0.44407 0.03451 12.869 1.51e-07 *** as.factor(A)2006 0.50271 0.03711 13.546 9.28e-08 *** --- 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
On peut alors envisager de simuler non pas des lois de Poisson, mais des lois quasipoissons, en reprenant un vieux code,
> 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) + }
A partir de ces deux remarques, on peut reprendre le code qui permettait de générer des scénarios de paiements pour les années futures,
> Yp=predict(regp,type="response",newdata=base) > Rs=rep(NA,20000) > for(s in 1:20000){ + serreur=sample(erreur, + size=36,replace=TRUE) + E=matrix(serreur,6,6) + sY=matrix(Yp,6,6)+E*sqrt(matrix(Yp,6,6)) + sbase=data.frame(sY=as.vector(sY),D,A) + sbase$sY[is.na(Y)==TRUE]=NA + sreg=glm(sY~as.factor(D)+as.factor(A), + data=sbase,family=poisson(link="log")) + sYp=predict(sreg,type="response", + newdata=sbase) + sYpscenario=rqpois(36,sYp,phi=3.18623) + Rs[s]=sum(sYpscenario[is.na(Y)==TRUE]) + }
si on regarde la densité de la distribution des paiements futurs, on obtient,
avec en bleu la distribution obtenue sur les résidus de Pearson brut, et en simulant une loi de Poisson, et en rouge la distribution des provisions avec les deux modifications. On observe clairement que les quantiles (par exemple) ont fortement changé. On peut le voir encore plus précisément sur un box-plot
Si on revient cinq minutes sur ce qu’on vient de faire. Dans le premier rappelons que l’on a besoin d’utiliser, pour prédire les paiements futurs, que
(avec la filtration correspondant à la partie supérieure du triangle disponible), i.e.
Comme on utilise un modèle de Poisson, on suppose aussi que
i.e.
c’est à dire que
Il faut avoir des résidus centrés et de variance unitaire. A première vue, c’est le cas pour nos résidus de Pearson (à peu près),
sauf que l’estimateur de la variance est biaisé. Certes, classiquement, l’estimateur est basé sur une normalisation par un facteur 
(classique en statistique quand on possède 
observations)
> sqrt(sum((residuals(regp,"pearson")- + mean(residuals(regp,"pearson")))^2)/ + (length(residuals(regp,"pearson"))-1)) [1] 1.261934
Sauf qu’ici on ne perd pas un degré de liberté (remplacer l’espérance des observations par leur moyenne empirique): en régression, il faut corriger par le nombre de variables explicatives. Et donc, pour avoir un estimateur sans biais de la variance de nos résidus, il faut corriger par un facteur 
. D’où l’utilisation des résidus dits ajustés afin d’avoir des résidus de variance (vraiment) unitaire.
Pour le second point, on utilise le fait qu’en pratique, la dispersion des paiements est plus grande que ce que donnerait un modèle de Poisson. Et il convient d’en tenir compte dans nos simulations. En l’occurrence, on veut que 
, mais aussi que
. On a vu que le paramètre de surdispersion ne pouvait être supposé unitaire. On a alors le choix, entre simuler une loi Gamma de paramètres
et 
ou bien simuler une loi binomiale négative, de moyenne 
, et de paramètre de surdispersion 
telle que la variance s’écrive
Aussi, on prend 
alors que
Avec ces méthodes, on obtient quelque chose… qui est programmé sous R,
> BootChainLadder(PAID,10000,"od.pois") BootChainLadder(Triangle = PAID, R = 10000, process.distr = "od.pois") Latest Mean Ultimate Mean IBNR SD IBNR IBNR 75% IBNR 95% 1 4,456 4,456 0.0 0.0 0 0 2 4,730 4,752 22.2 12.2 28 44 3 5,420 5,456 35.5 15.2 43 64 4 6,020 6,086 65.7 19.8 78 102 5 6,794 6,946 151.9 28.7 170 202 6 5,217 7,364 2,147.4 111.0 2,218 2,339 Totals Latest: 32,637 Mean Ultimate: 35,060 Mean IBNR: 2,423 SD IBNR: 132 Total IBNR 75%: 2,506 Total IBNR 95%: 2,653 > quantile(Rs,c(.75,.95)) 75% 95% 2509 2653
et qui donne les mêmes résultats que ce qu’on vient de reprogrammer…
Tail factor et clôture de la première ligne des triangles de paiements
La grosse hypothèse que l’on a fait quand on travaillait sur les triangles, était que la première année était close. Et donc, que pour la première ligne de paiements cumulés, le montant final était la charge totale des sinistres survenus cette année là. En fait, en 1999, Thomas Mack avait proposé un modèle relativement simple permettant l’inclusion d’un tail factor. L’idée est que l’on peut supposer que les facteurs de développement dans la méthode Chain Ladder décroissent tranquillement vers 1. Toujours sur le même triangle (mais en utilisant une régression pondérée sans constante – ce qui permet d’éviter les problèmes des valeurs manquantes à n’importe quel endroit dans le triangle), nous avions
> LAMBDA <- rep(NA,nc-1) > for(k in 1:(nc-1)){ + LAMBDA[k]=lm(PAID[,k+1]~0+PAID[,k], + weights=1/PAID[,k])$coefficients} > LAMBDA [1] 1.380933 1.011433 1.004343 1.001858 1.004735
Par “tranquillement”, on entend une décroissance exponentielle des coefficients avec le temps. On peut estimer un modèle linéaire sur le logarithme des facteurs de développement
> logL <- log(LAMBDA-1) > tps <- 1:(nc-1) > modele <- lm(logL~tps) > plot(tps,logL,xlim=c(1,20),ylim=c(-30,0),col="red") > abline(modele) > tpsP <- seq(6,1000) > logP <- predict(modele,newdata=data.frame(tps=tpsP)) > points(tpsP,logP ,pch=0,col="blue")
Si on regarde les coefficients de transition, on a
Pour passer de la dernière colonne observé à la charge ultime, il suffit de faire le produit des facteurs de transitions pour les années non encore observées, c’est à dire les points bleus
> (facteur <- prod(exp(logP)+1)) [1] 1.000707
Autrement dit, il faudrait rajouter 0.07% à la charge ultime. Sans ce facteur ultime, le montant de provisions était
> DIAG <- diag(PAID[,6:1]) > PRODUIT <- c(1,rev(LAMBDA)) > sum((cumprod(PRODUIT)-1)*DIAG) [1] 2426.985
mais si on rajoute le facteur que l’on vient de calculer, on passe à
On retrouve ici la même chose que ce qui est obtenu sous R,
> MackChainLadder(PAID,tail=TRUE) MackChainLadder(Triangle = PAID, tail = TRUE) Latest Dev.To.Date Ultimate IBNR Mack.S.E CV(IBNR) 1 4,456 0.999 4,459 3.15 0.299 0.0948 2 4,730 0.995 4,756 25.76 0.712 0.0277 3 5,420 0.993 5,460 39.64 2.528 0.0638 4 6,020 0.988 6,090 70.37 5.064 0.0720 5 6,794 0.977 6,952 157.99 31.357 0.1985 6 5,217 0.708 7,372 2,154.86 68.499 0.0318 Totals Latest: 32,637.00 Dev: 0.93 Ultimate: 35,088.76 IBNR: 2,451.76 Mack S.E.: 79.37 CV(IBNR): 0.0323730538389919
Triangles de paiements, ou charge dossier par dossier
Comme on l’avait vu en début de cours, on a en fait souvent plus d’information que les paiements effectués. On a aussi les estimations de charge effectuées par les gestionnaires de sinistres, pour les sinistres présents dans la base (i.e. les sinistres déclarés). Par exemple, sur le jeu de données utilisé en cours, on a les paiements, et l’estimation de la charge dossier par dossier,
> PAID
[,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
> INCURRED
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 4975 4629 4497 4470 4456 4456
[2,] 5135 4949 4783 4760 4750 NA
[3,] 5681 5631 5492 5470 NA NA
[4,] 6272 6198 6131 NA NA NA
[5,] 7326 7087 NA NA NA NA
[6,] 7353 NA NA NA NA NA
On peut alors utiliser la méthode Chain Ladder non seulement sur les paiements, mais aussi sur les estimation de charge dossier par dossier. Visuellement, on retrouve des courbes assez proches de celles vues en cours (certes, avec des développements mensuels),
> PCL$FullTriangle 1 2 3 4 5 6 1 3209 4372.000 4411.000 4428.000 4435.000 4456.000 2 3367 4659.000 4696.000 4720.000 4730.000 4752.397 3 3871 5345.000 5398.000 5420.000 5430.072 5455.784 4 4239 5917.000 6020.000 6046.147 6057.383 6086.065 5 4929 6794.000 6871.672 6901.518 6914.344 6947.084 6 5217 7204.327 7286.691 7318.339 7331.939 7366.656 > ICL$FullTriangle 1 2 3 4 5 6 1 4975 4629.000 4497.000 4470.000 4456.000 4456.000 2 5135 4949.000 4783.000 4760.000 4750.000 4750.000 3 5681 5631.000 5492.000 5470.000 5455.777 5455.777 4 6272 6198.000 6131.000 6101.117 6085.253 6085.253 5 7326 7087.000 6920.146 6886.416 6868.510 6868.510 6 7353 7129.075 6961.230 6927.300 6909.288 6909.288 > library(ChainLadder) > PCL <- MackChainLadder(PAID) > ICL <- MackChainLadder(INCURRED) > k <- 5 > plot(0:(nc+1-k),c(0,PCL$FullTriangle[k,1:(nc+1-k)]), + pch=19,type="b", + ylim=c(0,max(c(PCL$FullTriangle[k,],ICL$FullTriangle[k,]))), + xlim=c(0,nc), ylab="",xlab="",col="blue") > lines(0:(nc+1-k),c(0,ICL$FullTriangle[k,1:(nc+1-k)]), + pch=19,type="b",col="red") > lines((nc+1-k):nc,PCL$FullTriangle[k,(nc+1-k):nc], + pch=1,type="b",col="blue") > lines((nc+1-k):nc,ICL$FullTriangle[k,(nc+1-k):nc], + pch=1,type="b",col="red")
Il existe une méthode combinant les deux triangles, appelées Munich Chain Ladder, dont la théorie a été développée par Gerhard Quarg et Thomas Mack. Numériquement, on a
> (MNCL <- MunichChainLadder(Paid=PAID,Incurred=INCURRED)) MunichChainLadder(Paid = PAID, Incurred = INCURRED) Latest Paid Latest Incurred Latest P/I Ratio Ult. Paid Ult. 1 4,456 4,456 1.000 4,456 4,456 2 4,730 4,750 0.996 4,753 4,750 3 5,420 5,470 0.991 5,455 5,454 4 6,020 6,131 0.982 6,086 6,085 5 6,794 7,087 0.959 6,983 6,980 6 5,217 7,353 0.710 7,538 7,533 Totals Paid Incurred P/I Ratio Latest: 32,637 35,247 0.93 Ultimate: 35,271 35,259 1.00 > k <- 5 > plot(0:(nc+1-k),c(0,MNCL$MCLPaid[k,1:(nc+1-k)]), + pch=19,type="b", ylim=c(0,max(c(MNCL$MCLPaid[k,], + MNCL$MCLIncurred[k,]))),xlim=c(0,nc), + ylab="",xlab="",col="blue") > lines(0:(nc+1-k),c(0,MNCL$MCLIncurred[k,1:(nc+1-k)]), + pch=19,type="b",col="red") > lines((nc+1-k):nc,MNCL$MCLPaid[k,(nc+1-k):nc], + pch=1,type="b",col="blue") > lines((nc+1-k):nc,MNCL$MCLIncurred[k,(nc+1-k):nc], + pch=1,type="b",col="red")
Simulations et incréments négatifs
Pour faire suite (rapidement) à mon dernier billet, le cas le plus fréquent pour observer des incréments négatifs est lorsque l’on fait des simulations pour obtenir la distribution du montant de paiements futurs. Rappelons que pour générer un pseudo-triangle, on utilise les résidus de Pearson, et on pose
On voit que l’on risque d’avoir un incrément négatif si le résidu est négatif, et plus grand (en valeur absolue) que la racine carrée de la prédiction (par notre modèle Poissonnien, en l’occurrence). Or dans le triangle vu en cours, on a
> source("https://perso.univ-rennes1.fr/ arthur.charpentier/bases.R") > INC=PAID > INC[,2:6]=PAID[,2:6]-PAID[,1:5] > 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) > erreurs=residuals(reg,"pearson") > min(sqrt(Yp[is.na(Y)==FALSE])) [1] 2.868171 > min(erreurs) [1] -2.213449
autrement dit, on n’aura jamais d’incrément négatif lors de nos boucle, si l’on génère les résidus par bootstrap. Sinon, avec des lois paramétriques non-bornées inférieurement (loi normale, loi de Student), il est tout a fait possible d’obtenir des incréments négatifs. La méthode pour éviter le problème des incréments négatifs dans les boucles… est probablement de ne pas prendre en compte les scénarios où des incréments négatifs ont été obtenus, i.e.
R=rep(NA,10000) for(s in 1:10000){ serreur=sample(erreurs, size=36,replace=TRUE) E=matrix(serreur,6,6) sY=matrix(Yp,6,6)+E*sqrt(matrix(Yp,6,6)) if(min(sY[is.na(Y)==FALSE])>=0){ sbase=data.frame(sY=as.vector(sY),D,A) sbase$sY[is.na(Y)==TRUE]=NA sreg=glm(sY~as.factor(D)+as.factor(A), data=sbase,family=poisson(link="log")) sYp=predict(sreg,type="response", newdata=sbase) R[s]=sum(sYp[is.na(Y)==TRUE])} }
Une stratégie un peu plus propre pourrait être de mettre des 0 dès qu’on obtient un incrément négatif… Mais encore une fois, c’est de la cuisine.
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 (suite)
L’examen final pour le cours d’actuariat IARD aura lieu mardi prochain, et il sera basé sur des données mises en ligne sur le blog depuis une dizaine de jours. Je mets aujourd’hui en ligne les sorties qui seront fournies le jour de l’examen, et sur lesquelles des questions seront posées. J’encourage les étudiants à survoler ces sorties avant l’examen. En revanche, je ne répondrai à aucune question sur ce document ! Amenez vos calculatrices à l’examen.
« l’homme qui avait prédit la crise »
Je vais ressortir un peu ma casquette de schtroumpf grognon aujourd’hui, au risque de me faire encore taxer d’anti-journaliste primaire. Mais j’ai été agacé en lisant dans un article (journalistique, pas académique) qu’un personnage était présenté comme « l’homme qui avait prédit la crise ». En 
sept mots, le journaliste a réussi a utiliser trois clichés qui m’agacent.
Commençons par la fin, à savoir « la crise ». Sans vouloir jouer le schtroumpf grognon, je n’aime pas ce mot, abondamment utilisé par les journalistes… C’est quoi « la crise » ? Surtout que quand on pose la question, la réponse est toujours « allons, tu sais bien… la crise, quoi ». A la rigueur, «la crise financière »… mais laquelle ? « allons, arrête de faire du mauvais esprit… la crise de la dette, tu sais… ». Mais de quoi on parle là ? Des états européens sur-endettés qui risquent de se faire downgrader par les agences de notation ? Des dettes des particuliers qui explosent ? Des endettements des étudiants qui doivent faire face à la hausse des frais d’inscription ? Bref, je n’aime pas ce mot fourre tout qui ne peut être défini de manière univoque. Mais ce n’est pas ce qui m’agace le plus…
Il y a ensuite le mot « prédit ». Car c’est le mot qui est le plus utilisé par les journalistes, en France. Étrangement, le mot « prévu » n’est pas utilisé. Car «prédire » n’est pas « prévoir ». Certes, les deux se rapportent à l’avenir, de part le préfixe pré- qui exprime l’antériorité dans le temps. Le verbe « prédire » est souvent employé pour des pronostics, qui relèvent de l’intuition, d’un sentiment prémonitoire, voire d’une expérience surnaturelle: les voyantes « prédisent » l’avenir. Par contre, il existe des instituts de « prévision », ainsi que des cours de méthodes de « prévision » (ou encore des ouvrages, sur le sujet). On fait des « prévisions » météo, ou budgétaires. Bref, « prédire » relève de la foi, alors que « prévoir » relève davantage de la science. Il y a deux ans, j’avais déjà parlé des prédicateurs dans un billet suite au tremblement de terre qui avait eu lieu en mars 2009 en Italie (et je me permets de reprendre – plus ou moins – l’image extraite du plus formidable des albums de Tintin – d’où est aussi tiré l’image qui figure sur mon blog, et à laquelle je me suis pleinement identifiée). Et ce qui m’agace c’est que les journalistes préfèrent écouter les prédicateurs aux prévisionnistes… Et je ne parle pas des nombreux articles sur Paul le poulpe…
Enfin, il y a l’utilisation du terme « l’homme », ou plus précisément l’emploi du «l’», qui laisse entendre qu’il y a unicité: une seule personne aurait annoncé la crise. Vu de loin, c’est d’ailleurs amusant de noter que chaque journaliste a « son homme » qui aurait annoncé la crise. C’est Paul Jorion pour Rue89 (et France Culture), Nouriel Roubini pour le New York Times (et d’autres médias en France) mais certains pensent à Robert Shiller, Georges Magnus, Melchior Palyi, Victor Maslov, Raghuram Rajan, etc. Manifestement, il n’y a pas unicité, et pourtant tous les journalistes aiment utiliser ce pronom « l’ ». En fait, souvent les journalistes omettent de dire que cet homme est « l’homme » au sein des personnes faisant parti du réseau de personnes qu’ils lisent, et qui font partie de leur cercle médiatique. Ce qui est agaçant, c’est que les journalistes ne prennent pas la peine de lire les articles écrits par les économistes dans les revues d’économie. Malheureusement, dans un article académique, on ne prédit pas des crises, mais il y a de nombreuses analyses critiques qui permettent d’éclairer, et pas seulement ex-post. Mais lire des articles théoriques demande du temps, et des compétences…