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).

http://freakonometrics.blog.free.fr/public/perso4/plate-tekto.gif

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,

> max(which(1-game>.5))
[1] 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,

> plot(agenaissance$ANNEE+.5,agenaissance$EV,
+ type="l",lwd=2,col="purple")

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

http://freakonometrics.blog.free.fr/public/perso4/act2040-1.gif

car ici la probabilité d’avoir un accident est faible. Donc les deux fonctions liens sont proches. De plus, pour une loi de Poisson

http://freakonometrics.blog.free.fr/public/perso4/act2040-2.gif
http://freakonometrics.blog.free.fr/public/perso4/act2040-3.gif

Bref, pour notre loi de Poisson, on a presque une loi binomiale de paramètre http://freakonometrics.blog.free.fr/public/perso4/act2040-4.gif si http://freakonometrics.blog.free.fr/public/perso4/act2040-4.gif est petit.

Pour la prédiction suivant les trois modèles, on utilise classiquement

> exp(t(c(1,40,0,0,1,0,1,0))%*%reg9$coefficients)
           [,1]
[1,] 0.04351037

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

> sum(BASEN$nombre)/sum(BASEN$exposition)*
+ mean(BASEY$cout)
[1] 76.9754

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

> c(n1,n2)%*%t(c(c1,c2))
            1        1
[1,] 47.18972 46.20183
[2,] 72.32855 70.81440

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,

http://freakonometrics.blog.free.fr/public/perso4/res-boot-01.gif

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

http://freakonometrics.blog.free.fr/public/perso4/res-boot-2.gif

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

http://freakonometrics.blog.free.fr/public/perso4/res-boot.gif

et dans l’autre

http://freakonometrics.blog.free.fr/public/perso4/res-boot-3.gif

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

http://freakonometrics.blog.free.fr/public/perso4/res-boot-05.gif

(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,

plot(density(Rs))

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

http://freakonometrics.blog.free.fr/public/perso4/siiiim-01.gif

(avec la filtration correspondant à la partie supérieure du triangle disponible), i.e.

http://freakonometrics.blog.free.fr/public/perso4/siiiiim-03.gif

Comme on utilise un modèle de Poisson, on suppose aussi que

http://freakonometrics.blog.free.fr/public/perso4/siiiim-5.gif

i.e.

http://freakonometrics.blog.free.fr/public/perso4/siiiiim-08.gif

c’est à dire que

http://freakonometrics.blog.free.fr/public/perso4/siiiim-04.gif

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),

> mean(residuals(regp,"pearson"))
[1] -0.02462518
> sd(residuals(regp,"pearson"))
[1] 1.261934

sauf que l’estimateur de la variance est biaisé. Certes, classiquement, l’estimateur est basé sur une normalisation par un facteur http://freakonometrics.blog.free.fr/public/perso4/siiiim-11.gif (classique en statistique quand on possède http://freakonometrics.blog.free.fr/public/perso4/siiim-10.gif 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 http://freakonometrics.blog.free.fr/public/perso4/siiiim-12.gif. 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  https://perso.univ-rennes1.fr/arthur.charpentier/latex/qp001.png, mais aussi quehttps://perso.univ-rennes1.fr/arthur.charpentier/latex/qp002.png. 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

https://perso.univ-rennes1.fr/arthur.charpentier/latex/qp003.png et https://perso.univ-rennes1.fr/arthur.charpentier/latex/qp004.png

ou bien simuler une loi binomiale négative, de moyenne https://perso.univ-rennes1.fr/arthur.charpentier/latex/qp008.png, et de paramètre de surdispersion https://perso.univ-rennes1.fr/arthur.charpentier/latex/qp009.png telle que la variance s’écrive

https://perso.univ-rennes1.fr/arthur.charpentier/latex/qp010.png

Aussi, on prend https://perso.univ-rennes1.fr/arthur.charpentier/latex/qp011.png alors que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/qp012.png

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 à

> sum((cumprod(PRODUIT)*facteur-1)*DIAG)
[1] 2451.764

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

http://freakonometrics.blog.free.fr/public/perso4/pseudo-triangle.gifOn 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,
http://freakonometrics.blog.free.fr/public/perso4/glm-translation.gif
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.