Category Archives: Ratemaking

2011, ACT2040, Computer, GLM, Ratemaking, Regression

Confidence interval for predictions with GLMs

04/11/2011 Arthur Charpentier 5 Comments

Consider a (simple) Poisson regression . Given a sample where , the goal is to derive a 95% confidence interval for given , where is the prediction. Hence, we want to derive a confidence interval for the prediction, not the potential observation, i.e. the dot on the graph below

> r=glm(dist~speed,data=cars,family=poisson)
> P=predict(r,type="response",
+ newdata=data.frame(speed=seq(-1,35,by=.2)))
> plot(cars,xlim=c(0,31),ylim=c(0,170))
> abline(v=30,lty=2)
> lines(seq(-1,35,by=.2),P,lwd=2,col="red")
> P0=predict(r,type="response",se.fit=TRUE,
+ newdata=data.frame(speed=30))
> points(30,P1$fit,pch=4,lwd=3)

i.e.

Let denote the maximum likelihood estimator of . Then

where is Fisher information of (from standard maximum likelihood theory). Recall that

where computation of those values is based on the following calculations

In the case of the log-Poisson regression

Let us get back to our initial problem.

confidence interval for the linear combination

A first idea to get a confidence interval for is to get a confidence interval for (by taking exponential values of bounds, since the exponential is a monotone function). Asymptotically, we know that

thus, an approximation for the variance matrix of will be based on , obtained by plugging estimators of the parameters.
Then, since as an asymptotic multivariate distribution, any linear combination of the parameters will also be normal, i.e.
has a normal distribution, centered on , with variance where is the variance of . All those quantities can be easily computed. First, we can get the variance of the estimators

> i1=sum(predict(reg,type="response"))
> i2=sum(cars$speed*predict(reg,type="response"))
> i3=sum(cars$speed^2*predict(reg,type="response"))
> I=matrix(c(i1,i2,i2,i3),2,2)
> V=solve(I)

Hence, if we compare with the output of the regression,

> summary(reg)$cov.unscaled
(Intercept)         speed
(Intercept)  0.0066870446 -3.474479e-04
speed       -0.0003474479  1.940302e-05
> V
[,1]          [,2]
[1,]  0.0066871228 -3.474515e-04
[2,] -0.0003474515  1.940318e-05

Based on those values, it is easy to derive the standard deviation for the linear combination,

> x=30
> P2=predict(r,type="link",se.fit=TRUE,
+ newdata=data.frame(speed=x))
> P2
$fit
1
5.046034

$se.fit
[1] 0.05747075

$residual.scale
[1] 1

> sqrt(V[1,1]+2*x*V[2,1]+x^2*V[2,2])
[1] 0.05747084
> sqrt(t(c(1,x))%*%V%*%c(1,x))
[,1]
[1,] 0.05747084

And once we have the standard deviation, and normality (at least asymptotically), confidence intervals are derived, and then, taking the exponential of the bounds, we get confidence interval

> segments(30,exp(P2$fit-1.96*P2$se.fit),
+ 30,exp(P2$fit+1.96*P2$se.fit),col="blue",lwd=3)

Based on that technique, confidence intervals are no longer centered on the prediction. But who cares ?

delta method

Actually, those who like to use “more or less” expressions for confidence intervals will not like non centered intervals. So, an alternative is to use the delta method. Instead of writing (again) something on the theory, we can use a package which computes that method,

> estmean=t(c(1,x))%*%coef(reg)
> var=t(c(1,x))%*%summary(reg)$cov.unscaled%*%c(1,x)
> library(msm)
> deltamethod (~ exp(x1), estmean, var)
[1] 8.931232
> P1=predict(r,type="response",se.fit=TRUE,
+ newdata=data.frame(speed=30))
> P1
$fit
1
155.4048

$se.fit
1
8.931232

$residual.scale
[1] 1

The delta method gives us (asymptotic) normality, so once we have a standard deviation, we get the confidence interval.

> segments(30,P1$fit-1.96*P1$se.fit,30,
+ P1$fit+1.96*P1$se.fit,col="blue",lwd=3)

Note that those quantities – obtained with two different approaches – are rather close here

> exp(P2$fit-1.96*P2$se.fit)
1
138.8495
> P1$fit-1.96*P1$se.fit
1
137.8996
> exp(P2$fit+1.96*P2$se.fit)
1
173.9341
> P1$fit+1.96*P1$se.fit
1
172.9101

bootstrap techniques

And a third method (but far from what I expect to teach on that course) is to use bootstrap techniques to about those results based on asymptotic normality (we have only 50 observations). The idea is to sample from out dataset, and to run a log-Poisson regression on those new samples, and to repeat a lot of time,

ACT2040, Computer, GLM, Ratemaking, Regression, Spatial

Régression sur des variables catégorielles

26/10/2011 Arthur Charpentier 1 Comment

Petit complément par rapport au cours de mardi. On avait évoqué tout d’abord la lecture des sorties lorsque l’on régresse sur des variables catégorielles (des facteurs). Commençons par supprimer la constante de la régression

> reg0=glm(nbre~0+zone,offset=log(exposition),data=base, 
+ family=poisson(link="log"))
> summary(reg0)

Call:
glm(formula = nbre ~ 0 + zone, family = poisson(link = "log"), 
    data = base, offset = log(exposition))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.5717  -0.3968  -0.2996  -0.1547  12.6722  

Coefficients:
      Estimate Std. Error z value Pr(>|z|)    
zoneB -2.54187    0.06287  -40.43   <2e-16 ***
zoneA -2.54912    0.05285  -48.23   <2e-16 ***
zoneC -2.38525    0.03753  -63.56   <2e-16 ***
zoneD -2.13454    0.03878  -55.05   <2e-16 ***
zoneE -2.00204    0.03965  -50.49   <2e-16 ***
zoneF -2.06932    0.11547  -17.92   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 50966  on 50000  degrees of freedom
Residual deviance: 15692  on 49994  degrees of freedom
AIC: 20800

Number of Fisher Scoring iterations: 6

> predict(reg0,newdata=data.frame(
+ zone=c("A","B","C","D","E"),exposition=rep(1,5)))
        1         2         3         4         5 
-2.549120 -2.541870 -2.385253 -2.134543 -2.002044

On voit que toutes les modalités sont présentes, et toutes sont significatives. Si on régresse sur la constante, il faudra supprimer une modalité pour rendre le modèle identifiable. On peut forcer pour que la modalité de référence soit la seconde,

> base$zone=relevel(base$zone,"B")
> regB=glm(nbre~zone,offset=log(exposition),data=base,
+ family=poisson(link="log"))
> summary(regB)

Call:
glm(formula = nbre ~ zone, family = poisson(link = "log"), 
data = base,
offset = log(exposition))

Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.5717  -0.3968  -0.2996  -0.1547  12.6722

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.54187    0.06287 -40.431  < 2e-16 ***
zoneA       -0.00725    0.08213  -0.088 0.929661
zoneC        0.15662    0.07322   2.139 0.032432 *
zoneD        0.40733    0.07387   5.514 3.50e-08 ***
zoneE        0.53983    0.07433   7.263 3.80e-13 ***
zoneF        0.47255    0.13148   3.594 0.000325 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 15809  on 49999  degrees of freedom
Residual deviance: 15692  on 49994  degrees of freedom
AIC: 20800

Number of Fisher Scoring iterations: 6

> predict(regB,newdata=data.frame(
+ zone=c("A","B","C","D","E"),exposition=rep(1,5)))
1         2         3         4         5
-2.549120 -2.541870 -2.385253 -2.134543 -2.002044

On notera que les prédictions ne changent pas. On peut aussi choisir la première comme modalité de référence,

> base$zone=relevel(base$zone,"A")
> reg=glm(nbre~zone,offset=log(exposition),
> data=base,family=poisson(link="log"))
> summary(reg)

Call:
glm(formula = nbre ~ zone, family = poisson(link = "log"), 
data = base,
offset = log(exposition))

Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.5717  -0.3968  -0.2996  -0.1547  12.6722

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.54912    0.05285 -48.232  < 2e-16 ***
zoneB        0.00725    0.08213   0.088 0.929661
zoneC        0.16387    0.06482   2.528 0.011471 *
zoneD        0.41458    0.06555   6.324 2.54e-10 ***
zoneE        0.54708    0.06607   8.280  < 2e-16 ***
zoneF        0.47980    0.12699   3.778 0.000158 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 15809  on 49999  degrees of freedom
Residual deviance: 15692  on 49994  degrees of freedom
AIC: 20800

Number of Fisher Scoring iterations: 6

Le fait que la seconde modalité ne soit pas significative se lit par rapport à la modalité de référence (en l’occurrence la première): non significatif signifie alors non significativement différente. Autrement dit, on peut regrouper les modalités en une seule.

> base$zonesimple=base$zone
> base$zonesimple[base$zone%in%c("A","B")]="A"
> reg=glm(nbre~zonesimple,offset=log(exposition),
+ data=base,family=poisson(link="log"))
> summary(reg)

Call:
glm(formula = nbre ~ zonesimple, family = poisson(link = "log"),
data = base, offset = log(exposition))

Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.5717  -0.3959  -0.2989  -0.1547  12.6722

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.54612    0.04046 -62.937  < 2e-16 ***
zonesimpleC  0.16087    0.05518   2.915  0.00355 **
zonesimpleD  0.41158    0.05604   7.345 2.06e-13 ***
zonesimpleE  0.54408    0.05665   9.605  < 2e-16 ***
zonesimpleF  0.47681    0.12235   3.897 9.74e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 15809  on 49999  degrees of freedom
Residual deviance: 15692  on 49995  degrees of freedom
AIC: 20798

Number of Fisher Scoring iterations: 6

On note qu’avec ce regroupement, les autres modalités sont sensiblement différentes. On peut aussi retenir la troisième comme modalité de référence

> base$zonesimple=relevel(base$zonesimple,"C")
> reg=glm(nbre~zonesimple,offset=log(exposition),
+ data=base,family=poisson(link="log"))
> summary(reg)

Call:
glm(formula = nbre ~ zonesimple, family = poisson(link = "log"),
data = base, offset = log(exposition))

Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.5717  -0.3959  -0.2989  -0.1547  12.6722

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.38525    0.03753 -63.557  < 2e-16 ***
zonesimpleA -0.16087    0.05518  -2.915  0.00355 **
zonesimpleD  0.25071    0.05396   4.646 3.39e-06 ***
zonesimpleE  0.38321    0.05460   7.019 2.24e-12 ***
zonesimpleF  0.31593    0.12142   2.602  0.00927 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 15809  on 49999  degrees of freedom
Residual deviance: 15692  on 49995  degrees of freedom
AIC: 20798

Number of Fisher Scoring iterations: 6

Comme toutes les modalités semblent significatives, on peut tenter de prendre comme modalité de référence une des dernières (dont les estimations des coefficients donnent des résultats très proches)

> base$zonesimple=relevel(base$zonesimple,"F")
> reg=glm(nbre~zonesimple,offset=log(exposition),
+ data=base,family=poisson(link="log"))
> summary(reg)

Call:
glm(formula = nbre ~ zonesimple, family = poisson(link = "log"),
data = base, offset = log(exposition))

Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.5717  -0.3959  -0.2989  -0.1547  12.6722

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.06932    0.11547 -17.921  < 2e-16 ***
zonesimpleC -0.31593    0.12142  -2.602  0.00927 **
zonesimpleA -0.47681    0.12235  -3.897 9.74e-05 ***
zonesimpleD -0.06522    0.12181  -0.535  0.59232
zonesimpleE  0.06727    0.12209   0.551  0.58161
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 15809  on 49999  degrees of freedom
Residual deviance: 15692  on 49995  degrees of freedom
AIC: 20798

Number of Fisher Scoring iterations: 6

Au vue de cette dernière sortie, on peut tenter de fusionner toutes les dernières classes ensembles

> base$zonesimple[base$zone%in%c("D","E","F")]="F"
> reg=glm(nbre~zonesimple,offset=log(exposition),
+ data=base,family=poisson(link="log"))
> summary(reg)

Call:
glm(formula = nbre ~ zonesimple, family = poisson(link = "log"),
data = base, offset = log(exposition))

Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.5660  -0.3959  -0.3004  -0.1547  12.5929

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.07182    0.02696 -76.853  < 2e-16 ***
zonesimpleC -0.31344    0.04621  -6.783 1.18e-11 ***
zonesimpleA -0.47431    0.04861  -9.757  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 15809  on 49999  degrees of freedom
Residual deviance: 15698  on 49997  degrees of freedom
AIC: 20800

Number of Fisher Scoring iterations: 6

Bon, formellement, regrouper deux modalités (i.e. décréter que deux variables sont non significatives simultanément) demande un peu plus qu’un test de Student, ou que deux tests de Student…. Si on remonte un peu en arrière, on aurait pu faire un test multiple avant de fusionner les trois modalités (un test de type Fisher, ou une ANOVA)

> base$zonesimple=relevel(base$zonesimple,"F")
> reg=glm(nbre~zonesimple,offset=log(exposition),
+ data=base,family=poisson(link="log"))
> library(car)
> linearHypothesis(reg,c("zonesimpleD=0","zonesimpleE=0"))
Linear hypothesis test

Hypothesis:
zonesimpleD = 0
zonesimpleE = 0

Model 1: restricted model
Model 2: nbre ~ zonesimple

Res.Df Df  Chisq Pr(>Chisq)
1  49997
2  49995  2 5.7073    0.05763 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Manifestement, on peut accepter l’hypothèse que ces trois catégories n’en font qu’une. La zone géographique peut alors se découper en trois grandes zones, et pas en six. On notera que cela correspond à ce que propose un arbre de régression

> library(tree)
> arbre=tree(nbre~zone+offset(log(exposition)),
+ data=base,split="gini")
> plot(arbre)
> text(arbre)

2011, ACT2040, Computer, GLM, Ratemaking, Regression

ACT2040, introduction aux modèles linéaires généralisés

03/10/2011 Arthur Charpentier Leave a comment

On commencera ce mardi les GLM, après avoir introduit les lois exponentielles (qui ont du être revues en démonstration vendredi dernier). La notation utilisée sera que la loi (densité ou fonction de probabilité) de sera de la forme

http://freakonometrics.blog.free.fr/public/perso4/loi-exponentielle.gif

Pour un complément plus exhaustif, je renvoie au chapitre en ligne.

Le modèle linéaire (Gaussien)

Le modèle de base est le modèle Gaussien que l’on avait revu au dernier cours,

> X=c(1,2,3,4)
> Y=c(1,2,5,6)
> base=data.frame(X,Y)
> reg1=lm(Y~1+X,data=base)
> nbase=data.frame(X=seq(0,5,by=.1))
> Y1=predict(reg1,newdata=nbase)

Pour une prédiction (unique), on obtient la prédiction suivante

Le code pour une telle représentation est le suivant

> plot(X,Y,pch=3,cex=1.5,lwd=2,xlab="",ylab="")
> lines(nbase$X,Y1,col="red",lwd=2)
> u=2
> mu=predict(reg1)[2]
> sigma=summary(reg1)$sigma
> y=seq(0,7,.05)
> loi=dnorm(y,mu,sigma)
> segments(u,y,loi+u,y,col="light green")
> lines(loi+u,y)
> abline(v=u,lty=2)
> points(X[2],Y[2],pch=3,cex=1.5,lwd=2)
> points(X[2],predict(reg1)[2],pch=19,col="red")
> arrows(u-.2,qnorm(.05,mu,sigma),
+ u-.2,qnorm(.95,mu,sigma),length=0.1,code=3,col="blue")

On peut multiplier les prédictions, en se basant sur l’hypothèse d’homoscédasticité (la variance sera alors constante)

Mais on peut aller plus loin

Le modèle linéaire généralisé

Plusieurs modèles peuvent etre estimés, en changeant la loi de la variable à expliquer, et la fonction lien,

> reg2=glm(Y~1+X,data=base,family=poisson(link="identity"))
> Y2=predict(reg2,newdata=nbase,type="response")
> reg3=glm(Y~1+X,data=base,family=poisson(link="log"))
> Y3=predict(reg3,newdata=nbase,type="response")
> reg4=glm(Y~1+X,data=base,family=gaussian(link="log"))
> Y4=predict(reg4,newdata=nbase,type="response")
> sigma=sqrt(summary(reg4)$dispersion)

Pour le modèle Poissonnien avec un lien identité, on obtient

On obtient ainsi une variance qui augmente avec la prédiction,

Pour une régression de Poisson avec un lien logarithmique,

i.e. pour nos quatre prédictions

On peut comparer avec une prédiction d’un modèle Gaussien avec un lien logarithmique,

i.e. pour les quatre prédictions

2011, ACT2040, Computer, Graphics, Premium principle, Ratemaking, Statistics

Visualisation en tarification, avec R

21/09/2011 Arthur Charpentier Leave a comment

Chose promise, chose due: mardi prochain, le 27 septembre, le cours d’Actuariat IARD se déroulera en salle informatique (PK-S1525) à l’horaire habituel. Sinon quelques notes de cours sur la tarification a priori sont en ligne ici.

Quelques références en R

Avant de commencer à programmer en R, quelques références qui peuvent être utiles, en français pour commencer,

“R pour les débutants” d’Emmanuel Paradis, (PDF)
“Introduction à la programmation en S” par Vincent Goulet, (PDF)
“Statistique de l’Assurance” par Arthur Charpentier (PDF) qui sont mes notes de cours de l’an passé, mais qui insistent sur l’utilisation de R, pas sur la programmation en R, qui est supposée un peu connue.

Sinon en anglais, les références sont un peu plus nombreuses,

“R for Beginners” d’Emmanuel Paradis (PDF),
“An Introduction to R” par Longhow Lam (PDF)
“The R language — a short companion” par Marc Vandemeulebroecke (PDF),
“The R Guide” par Jason Owen (PDF),
“Econometrics in R” par Grant Farnsworth (PDF) pour aller plus loin sur les régressions,
“Practical Regression and Anova using R” by Julian Faraway (PDF) sur le meme sujet
“Statistics with R and S-Plus” d’Hugo Quené (PDF)
“Statistical Computing and Graphics Course Notes” par Frank Harrell, (PDF).
“Using R for Data Analysis and Graphics – Introduction, Examples and Commentary” par John Maindonald (PDF).

Le code pour importer les données est le suivant,

> sinistre <- read.table("http://freakonometrics.free.fr/sinistreACT2040.txt",
+ header=TRUE,sep=";")
> sinistres=sinistre[sinistre$garantie=="1RC",]
> contrat <- read.table("http://freakonometrics.free.fr/contractACT2040.txt",
+ header=TRUE,sep=";")
> T=table(sinistres$nocontrat)
> T1=as.numeric(names(T))
> T2=as.numeric(T)
> nombre1 = data.frame(nocontrat=T1,nbre=T2)
> I = contrat$nocontrat%in%T1
> T1= contrat$nocontrat[I==FALSE]
> nombre2 = data.frame(nocontrat=T1,nbre=0)
> nombre=rbind(nombre1,nombre2)
> basenb = merge(contrat,nombre)
> head(basenb)
> basesin=merge(sinistres,contrat)
> basesin=basesin[basesin$cout>0,]

Faire des graphiques en R

Dans A practicioner’s guide to Generalized Linear Models (en ligne ici), on peut voir des graphiques comme celui ci-dessous, avec la fréquence de sinistre en fonction de l’age (ou plutôt de classes d’ages). Sur le graphique ci-dessous, on a également une distinction entre les hommes (en bleu) et les femmes (en rouge),

sur ce graphique, les fréquences sont exprimées en logarithme du multiplicateur (i.e. en variation par rapport à la moyenne du portefeuille).
On peut utiliser le code suivant pour générer (automatiquement) des graphiques similaires,

> graphique=function(nom="ageconducteur",
+ niveau=c(17,21,24,29,34,44,64,84,110),
+ continu=TRUE,type=1){
+ if(continu==TRUE){X=cut(basenb[,nom],niveau)}
+ if(continu==FALSE){X=as.factor(basenb[,nom])}
+ E=basenb$exposition
+ Y=basenb$nbre
+ FREQ=levels(X)
+ moyenne=variance=n=rep(NA,length(FREQ))
+ for(k in 1:length(FREQ)){
+ moyenne[k] =weighted.mean(Y[X==FREQ[k]]/E[X==FREQ[k]],
+ E[X==FREQ[k]])
+ variance[k]=weighted.mean((Y[X==FREQ[k]]/E[X==FREQ[k]]-
+ moyenne[k])^2,E[X==FREQ[k]])
+ n[k]       =sum(E[X==FREQ[k]])
+}
+ w=barplot(n,names.arg=FREQ,col="light green",axes=FALSE,
+ xlim=c(0,1.2*length(FREQ)+.5))
+ mid=w[,1]
+ axis(2)
+ par(new=TRUE)
+ IC1=moyenne+1.96/sqrt(n)*sqrt(variance)
+ IC2=moyenne-1.96/sqrt(n)*sqrt(variance)
+ moyenneglobale=sum(Y)/sum(E)
+ 
+ if(type==1){
+ plot(mid,moyenne,ylim=range(c(IC1,IC2)),type="b",
+ col="red",axes=FALSE,xlab="",ylab="",
+ xlim=c(0,1.2*length(FREQ)+.5))
+ segments(mid,IC1,mid,IC2,col="red")
+ segments(mid-.1,IC1,mid+.1,IC1,col="red")
+ segments(mid-.1,IC2,mid+.1,IC2,col="red")
+ points(mid,moyenne,pch=19,col="red")
+ axis(4)
+ abline(h=moyenneglobale,lty=2,col="red")}
+
+ if(type==2){
+ plot(mid,log(moyenne/moyenneglobale),ylim=
+ range(c(log(IC1/moyenneglobale),log(IC2/moyenneglobale))),
+ type="b",col="red",axes=FALSE,xlab="",ylab="",
+ xlim=c(0,1.2*length(FREQ)+.5))
+ segments(mid,log(IC1/moyenneglobale),mid,
+ log(IC2/moyenneglobale),col="blue")
+ segments(mid-.1,log(IC1/moyenneglobale),mid+.1,
+ log(IC1/moyenneglobale),col="blue")
+ segments(mid-.1,log(IC2/moyenneglobale),mid+.1,
+ log(IC2/moyenneglobale),col="blue")
+ points(mid,log(moyenne/moyenneglobale),pch=19,col="red")
+ axis(4)
+ abline(h=0,lty=2,col="red")}
+
+ mtext("Exposition", 2, line=2, cex=1.2,col="light green")
+ if(type==1){mtext("Fréquence annualisée", 
+     4, line=-2, cex=1.2,col="red")}
+ if(type==2){mtext("Fréquence annualisée (log multiplicateur)", 
+    4, line=-2, cex=1.2,col="red")}
+ }

Par exemple en utilisant un découpage arbitraire par classe d’age (comme cela est fait par défaut dans la fonction),

> graphique()

Mais on pourrait aussi utiliser un découpage assurant d’avoir des classes plus grandes, par exemple en utilisant les quantiles,

> Q=quantile(basenb[,"ageconducteur"],(0:10)/10)
> Q[1]=Q[1]-1
> graphique(nom="ageconducteur",niveau=Q,continu=TRUE)

Ces deux graphiques permettent de visualiser la fréquence empirique de sinistre, par classe d’age, sans modèle paramétrique. Des intervalles de confiance sont également représentés (basés sur une hypothèse de normalité). Notons que l’on peut aussi faire une représentation relative à la fréquence moyenne (en log des multiplicateurs),

> graphique(type=2)

ou encore, si on cherche à analyser la fréquence en fonction de la zone géographique d’habitation

> graphique(nom="zone",continu=FALSE,type=2)

Une autre utilisation peut être faite sur la sévérité des sinistres (coût moyen) et la fréquence, par exemple un assureur,

Il est possible de modifier un peu la fonction pour ajouter au graphique la sévérité des sinistres, e.g.

> graphiquecout=function(nom="ageconducteur",
+ niveau=c(17,21,24,29,34,44,64,84,110),
+ continu=TRUE,type=1){
+ if(continu==TRUE){X=cut(basenb[,nom],niveau)}
+ if(continu==FALSE){X=basenb[,nom]}
+ E=basenb$exposition
+ Y=basenb$nbre
+ FREQ=levels(X)
+ moyennen=variancen=nn=rep(NA,length(FREQ))
+ for(k in 1:length(FREQ)){
+ moyennen[k] =weighted.mean(Y[X==FREQ[k]]/E[X==FREQ[k]],
+ E[X==FREQ[k]])
+ variancen[k]=weighted.mean((Y[X==FREQ[k]]/E[X==FREQ[k]]-
+ moyennen[k])^2,E[X==FREQ[k]])
+ nn[k]       =sum(E[X==FREQ[k]])
+ }
+ moyenneglobalen=sum(Y)/sum(E)
+ 
+ if(continu==TRUE){X=cut(basesin[,nom],niveau)}
+ if(continu==FALSE){X=basesin[,nom]}
+ Y=basesin$cout
+ FREQ=levels(X)
+ moyennes=variances=ns=rep(NA,length(FREQ))
+ for(k in 1:length(FREQ)){
+ moyennes[k] =mean(Y[X==FREQ[k]])
+ variances[k]=var(Y[X==FREQ[k]])
+ ns[k]=length(Y[X==FREQ[k]])
+ }
+ moyenneglobales=mean(Y)
+ 
+ w=barplot(nn,names.arg=FREQ,col="light green",
+ axes=FALSE,xlim=c(0,1.2*length(FREQ)+.5))
+ mid=w[,1]
+ 
+  par(new=TRUE)
+ IC1=moyennen+1.96/sqrt(nn)*sqrt(variancen)
+ IC2=moyennen-1.96/sqrt(nn)*sqrt(variancen)
+ plot(mid,moyennen,ylim=range(c(IC1,IC2)),type="b",
+ col="red",axes=FALSE,xlab="",ylab="",
+ xlim=c(0,1.2*length(FREQ)+.5))
+ segments(mid,IC1,mid,IC2,col="red")
+ segments(mid-.1,IC1,mid+.1,IC1,col="red")
+ segments(mid-.1,IC2,mid+.1,IC2,col="red")
+ points(mid,moyennen,pch=19,col="red")
+ axis(4)
+ abline(h=moyenneglobalen,lty=2,col="red")
+ 
+ par(new=TRUE)
+ IC1=moyennes+1.96/sqrt(ns)*sqrt(variances)
+ IC2=moyennes-1.96/sqrt(ns)*sqrt(variances)
+ plot(mid,moyennes,ylim=range(c(IC1,IC2)),type="b",
+ col="blue",axes=FALSE,xlab="",ylab="",
+ xlim=c(0,1.2*length(FREQ)+.5))
+ segments(mid,IC1,mid,IC2,col="blue")
+ segments(mid-.1,IC1,mid+.1,IC1,col="blue")
+ segments(mid-.1,IC2,mid+.1,IC2,col="blue")
+ points(mid,moyennes,pch=19,col="blue")
+ axis(2)
+ abline(h=moyenneglobales,lty=2,col="blue")
+ 
+ mtext("Cout moyen", 2, line=2, cex=1.2,col="blue")
+ mtext("Fréquence annualisée", 4, line=-2, cex=1.2,col="red")
+ }

Si on regarde par classe d’age,

> graphiquecout()

ou encore, en fonction de l’age du véhicule,

> Q=quantile(basenb[,"agevehicule"],(0:10)/10)
> Q[1]=Q[1]-1
> graphiquecout(nom="agevehicule",niveau=Q,continu=TRUE)

ou enfin, en fonction de la zone géographique

> graphiquecout(nom="zone",continu=FALSE)

2011, ACT2040, Ratemaking

ACT2040, exercices de tarification

07/09/2011 Arthur Charpentier Leave a comment

Les exercices commenceront ce vendredi, et les énoncés (ainsi que pour quelques uns des pistes de correction) sont en ligne ici. Les exercices (et les corrections) ont été tapés pour ce cours par Jean-Philippe Boucher. La premiere série d’exercices sont des exercices de remise en route sur les calculs d’espérances, et d’espérance conditionnelle. On reviendra dessus mardi prochain.

Actuarial science, Maths, Premium principle, Probability, Publications, Ratemaking

Loi des grands nombres et théorème central limite en assurance

17/06/2011 Arthur Charpentier 4 Comments

Dans le dernier numéro de la revue Risques figure un petit article sur la loi des grands nombres et le théorème central limite en assurance (ici) Le but était de remettre un peu les pendules à l’heure sur la différence entre les implications de ces deux résultats, et surtout d’amorcer une discussion sur l’importance des hypothèses (en particulier celle de risques indépendants).

Ratemaking

Pour une tarification de l’assurance automobile à l’aide du tour de poitrine !

25/02/2011 Arthur Charpentier 4 Comments

Plusieurs sites spécialisés en assurance commencent à évoquer un arrêté probable de la cour européenne sur la discrimination en assurance (par exemple ici ou là). Une des bases (économiques) de l’assurance est le principe d’Akerlof qui pousse les assureurs à segmenter par classe de risque. Afin de segmenter, et de révéler les classes de risques, on utilise l’historique de sinistralité (information dite a posteriori), ou bien des informations exogènes (dites a priori) sur le conducteur, le véhicule, son usage, etc. Par exemple on peut utiliser l’ancienneté du véhicule, et le nombre de kilomètre effectués (en moyenne) par le conducteur, comme sur le graphique ci-dessous (retrouvé dans les transparents que l’on utilisait avec François Bucchini quand on donnait le cours d’assurance dommage à l’ENSAE, les probas étant “normalisées” dans une espère de base 100)

ou encore le type de carburant utilisé (diesel ou essence)

On retrouve que plus on conduit, plus la probabilité d’avoir un accident augmente, mais le carburant et l’âge du véhicule semblent être aussi des variables discriminantes. Et parmi les variables qui semblent significatives (pour expliquer la probabilité d’avoir un accident), il y a le sexe (croisé ici avec le kilométrage, comme auparavant),

Alors l’effet peut sembler marginal sur ce graphique… mais c’est loin d’être le cas. Par exemple, sans utiliser des techniques très poussées en économétrie, on peut regarder le nombre moyen de sinistres, et le coût moyen de sinistres, par sexe, et par tranche d’âge (voire aussi par CSP et par puissance du véhicule). Dans une étude faite par un assureur, j’avais trouvé les chiffres suivants

En haut à droite (beaucoup d’accidents, et coût – en moyenne – élevé) on retrouve les jeunes hommes. Donc oui, les jeunes hommes sont significativement beaucoup plus risqués que les autres conducteurs. Et le soucis est que, si on ne segmente pas, Georges Akerlof nous explique que le marché de l’assurance disparait, les “bons” risques ne voulant plus payer pour les “mauvais” risques. Sans pour autant rentrer dans une spirale infernale de la segmentation, il est bon que les primes restent corrélées au risque sous-jacent.

Les assureurs prétendent qu’ils ne «ne font pas de la discrimination, ils font de la différenciation ». Je ne rentrerais pas sur les débats de terminologie (pas aujourd’hui en tous les cas), mais le but n’est pas de trouver des variables “explicatives” de la sinistralité au sens causal (malgré la terminologie usuelle des économètres) mais de trouver des variables “corrélées” avec une forte sinistralité, et de les utiliser pour segmenter. Les assureurs européen avaient, jusqu’alors, bénéficié d’un sursis dans le calcul des primes qui leur permet de pratiquer des tarifs différents « lorsque le sexe est un facteur déterminant dans l’évaluation des risques».

Dans un cours d’analyse de données, j’avais montré (ici) qu’à partir des notes de étudiant(e)s à différents examens, je pouvais prédire le sexe des étudiants. Bon, l’étude avait été faite rapidement, avec un petit jeu de données (et donc sans population d’apprentissage et de test), mais il est facile de trouver des variables permettant de deviner le sexe d’un conducteur. D’aucuns pourraient être tentés d’utiliser la pointure des chaussures, mais personnellement je préférerais le tour de poitrine, ou un tour de poitrine ramené à un tour de hanche. Je suis presque sûr qu’avec de telles observations, on peut avoir des variables fortement corrélées avec la survenance d’accident ! En tous les cas ça promet un peu d’animation chez les agents d’assurance ! voire chez les chirurgiens esthétiques (retirer les implants mammaires pour faire baisser sa prime d’assurance auto, voilà qui est original) !

ACT6420, Ratemaking, SAS

La tarification avec SAS

22/09/2010 Arthur Charpentier 5 Comments

En tarification, il est possible d’utiliser d’autres logiciels que R, en particulier, il semble que l’on puisse faire deux ou trois choses avec SAS…. J’en parle un peu car il semble que, paradoxalement, les asssureurs préfèrent encore SAS à R (par exemple). Et comme plusieurs étudiants m’avaient demandé “et comment on fait avec SAS ?“. Bon, par contre je ne mets que les choses de base, parce que SAS est assez limité sur ce qu’il peut faire….

Pour suivre un peu le plan du cours, la première étape est de définir une variable d’exposition dans la table,

DATA contrats;
SET lib.contrats;
lnexpo = log(expo);
run;

Pour faire une régression de Poisson, ce n’est pas forcément compliqué,

PROC GENMOD DATA = base;
ODS OUTPUT ParameterEstimates=Genmod1_Param
           Type3=Genmod1_Var
           Modelfit=Genmod1_InfoModele; 
MODEL nbsin = ageconducteur /
                  dist = poisson   
                  link = log   
                  offset = lnexpo 
                  type3;
RUN; QUIT;

La sortie SAS a alors l’allure suivante

                                  The GENMOD Procedure
                    Critère pour évaluer la qualité de l'ajustement
              Critère                   DF          Valeur       Valeur/DF
              Deviance                63E3      26872.5334          0.4237
              Scaled Deviance         63E3      26872.5334          0.4237
              Pearson Chi-Square      63E3      73275.5362          1.1553
              Scaled Pearson X2       63E3      73275.5362          1.1553
              Log Likelihood                   -18474.2667

       Algorithm converged.
                       Analyse des résultats estimés de paramètres

                              Erreur      Wald 95Limites
Paramètre    DF   Estimation   standard      de confiance %    Khi 2   Pr > Khi 2
Intercept     1      -3.5164     0.0851    -3.6832  -3.3496   1708.02       <.0001
ageconducteur 1       0.0168     0.0014     0.0141   0.0195    146.73       <.0001
Scale         0       1.0000     0.0000     1.0000   1.0000
NOTE: The scale parameter was held fixed.

                         Statistiques LR pour Analyse de Type 3
                      Source           DF      Khi 2    Pr > Khi 2
                      ageconducteur     1     148.72        <.0001

Il est aussi possible de faire des GAM (i.e. du lissage de la variable explicative – continue – avec des fonctions splines)

PROC GAM DATA = base;
MODEL nbsin = spline(ageconducteur) / dist = Poisson;
OUTPUT OUT=gam PREDICTED; 
RUN; QUIT;
PROC SORT DATA = gam NODUPKEY; BY age_cond; RUN; QUIT;

et on peut faire des prédictions avec ce modèle (la sortie n’apporte pas grand chose, en pratique),

DATA gam;
SET gam;
pred_nbsin_gam = exp(P_nbsin);
KEEP ageconducteur pred_nbsin_gam;
RUN;

Enfin, on peut tenter de faire un joli graphique. Pour cela, on calcule les prédictions de trois modèles, le premier étant des nombres moyens de sinistres par âge

PROC SORT DATA = base; BY ageconducteur; RUN; QUIT;
PROC MEANS DATA = base NOPRINT;
BY ageconducteur;
VAR nbsin;
WEIGHT expo;
OUTPUT OUT = nbsin_age (DROP = _TYPE_ _FREQ_) MEAN=mo
y_uni_nbsin;
RUN; QUIT;

ensuite, on fait un modèle GLM, et un modèle GAM, et on récupère les sorties

PROC SORT DATA = nbsin_age; BY age_cond; RUN; QUIT;
PROC SORT DATA = gam; BY age_cond; RUN; QUIT;
DATA nbsin_age;
MERGE nbsin_age
      gam;
BY age_cond;
RUN;

On essaye de faire le dessin (je passe les lignes de commande, il y en a une vingtaine)

Pour faire une régression quasiPoisson, le code a l’allure suivante,s

PROC GENMOD DATA = base;
ODS OUTPUT ParameterEstimates=Genmod1bis_Param
           Type3=Genmod1bis_Var
           Modelfit=Genmod1bis_InfoModele; 
MODEL nbsin = ageconducteur /
                 dist = poisson  
                 link = log      
                 offset = lnexpo 
                 type3           
                 scale = deviance;
RUN; QUIT;

La sortie donne alors l’estimation du paramètre de surdispersion (ou sur cet exemple de sousdispersion)

                      Analyse des résultats estimés de paramètres

                                   Erreur    Wald 95Limites
Paramètre      DF   Estimation   standard    de confiance %     Khi 2   Pr > Khi 2

Intercept        1     -3.5164     0.0554  -3.6249  -3.4078   4031.29       <.0001
ageconducteur    1      0.0168     0.0009   0.0150   0.0186    346.32       <.0001
Scale            0      0.6509     0.0000   0.6509   0.6509

On notera que pour calculer le critère d’Akaike, ça n’est pas forcément trivial,

%MACRO CALCUL_AIC_BIC(infomodel=, param=);
    DATA _null_;
    SET &infomodel.;
    IF Criterion = "Log Likelihood" THEN CALL SYMPUT("Loglike", Value);
    IF Criterion = "Deviance" THEN CALL SYMPUT("n_etoile", Df);
    RUN;
    DATA _null_;
    SET &param.  end=fin;
    RETAIN nb_df 0;
    nb_df = nb_df + df;
    IF fin THEN CALL SYMPUT("k", nb_df);
    RUN;
    DATA genmod_aic_bic;
    SET &param.;
    FORMAT Loglike 12.2 K 10. N 10. AIC_CALC 12.2 BIC_CALC 12.2;
    Loglike = 0; K = 0; N = 0; AIC_CALC = 0; BIC_CALC = 0;
    IF Parameter = "Intercept";
    KEEP Loglike K N AIC_CALC BIC_CALC;
    RUN;
    DATA genmod_aic_bic;
    SET genmod_aic_bic;
    Loglike = &loglike.;
    K = &k.;
    N = %eval(&n_etoile. + &k.);
    AIC_CALC = 2 * Loglike + 2 * K;
    BIC_CALC = 2 * Loglike + K * log(N);
    RUN;
    PROC PRINT DATA = genmod_aic_bic;
    RUN; QUIT;
%MEND CALCUL_AIC_BIC;
%CALCUL_AIC_BIC(infomodel=Genmod2_InfoModele, param=Genmod2_Param);

Ratemaking

Les méchanismes bonus-malus

28/10/2009 Arthur Charpentier 1 Comment

Les chaînes de Markov

L’outil de base pour l’étude des systèmes bonus-malus est l’utilisation des chaînes de Markov, dès lors que l’on travaille sur des classes de bonus. Pour plus de compléments sur les chaînes de Markov, je renvoie à des vieilles notes de cours (inachevée, et sûrement remplis de coquillettes) ici.
Pour reprendre les grandes idées, on considère une suite de variables aléatoires à valeurs dans $https://perso.univ-rennes1.fr/arthur.charpentier/latex/calE.png$ , $https://perso.univ-rennes1.fr/arthur.charpentier/latex/suiteXn.png$ . Cette suite est sera une chaîne de Markov si pour tout $https://perso.univ-rennes1.fr/arthur.charpentier/latex/ngeq1.png$ , et pour tout

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

telle que

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/markov-proba1.png$

on ait

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/proab1.png$
= $https://perso.univ-rennes1.fr/arthur.charpentier/latex/proab2.png$

On parlera d’absence de mémoire (au sens où seule la dernière valeure intervient dans le calcul de la loi conditionnelle, et pas l’intégralité de la trajectoire). Si $https://perso.univ-rennes1.fr/arthur.charpentier/latex/suiteXn.png$ est une chaîne de Markov de distribution initiale $https://perso.univ-rennes1.fr/arthur.charpentier/latex/lambdaV.png$ et de matrice de transition P, pour tout $https://perso.univ-rennes1.fr/arthur.charpentier/latex/nkentier.png$ ,

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/proba-markov-statioonnaire.png$

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

Pour étudier l’évolution sur le long terme, on peut utiliser l’équation de Chapman-Kolmogorov: si $https://perso.univ-rennes1.fr/arthur.charpentier/latex/suiteXn.png$ est une chaîne de Markov de distribution initiale $https://perso.univ-rennes1.fr/arthur.charpentier/latex/lambdaV.png$ et de matrice de transition P, pour tout
$https://perso.univ-rennes1.fr/arthur.charpentier/latex/nkentier.png$ , pour tout $https://perso.univ-rennes1.fr/arthur.charpentier/latex/zvze54v6.png$ ,

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/chapman-kol.png$

Un corrolaire important est que si

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

alors

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

Bref, on peut facilement étudier l’évolution à long terme de la chaîne de Markov en faisant un peu de calcul matriciel, i.e. des calculs de puissances de matrices.
En fait, formellement, on peut soit étudier des matrices, soit des graphes. Par exemple,

Si on veut étudier l’évolution sur deux périodes, il suffit de calculer le carré de la matrice, i.e.

Les systèmes bonus-malus par classes

Les tableaux ci-dessous sont tirés des annexes de la bible de Jean Lemaire, Bonus-Malus systems in automobile insurance. Le plus simple est probablement celui de Hong Kong, qui est Markovien (à l’ordre 1),

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.BM-hong-kong-B-9_m.jpg

L’exemple ci-dessous est celui qui a existé en Allemagne dans les années 80,

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.BM-germany-B-7_m.jpg

puis celui qui existait en Finlande,

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.BM-finnish-B-6_m.jpg

et enfin le système italien du début des années 90,

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.BM-italy-B-11_m.jpg

Pour tous ces sytèmes, une fois spécifiée la loi du nombre de sinistres par an, il est possible de construire la matrice de transition des classes de bonus sur un an, ou plus.
En particulier, l’étude de la loi invariante de ce genre de sytème (que je mentionnais en introduction) permet d’étudier la viabilité du système à plus ou moins long terme, et la vitesse de la convergence vers l’état stationnnaire.

Le système bonus-malus français

En France, le mécanisme de bonus-malus peut paraître un peu plus compliqué, et il fonctionne de la manière suivante.

Le calcul du bonus-malus s’ef
fectue en fonction des sinistres survenus (ou non) pendant la période de douze mois consécutifs précédant de deux mois l’échéance annuelle du contrat d’assurance.
Seuls les accidents responsables donnent lieu à malus. C’est à dire qu’en cas de foce majeur, ou si l’accident est entièrement imputable à un tiers, alors le conducteur n’aura pas de malus. Notons aussi que Le malus est rattaché au contrat d’assurance. Autrement dit, si quelqu’un cause un accident avec un véhicule assuré par un tiers, il se supportera pas le malus.
La prime de référence est la prime qui sera utilisée pour ensuite actualiser le montant de prime en tenant compte du coefficient de majoration ou minoration. Il s’agit de la prime technique calculée par l’assureur en fonction de charactéristiques a priori, comme la puissance du véhicule véhicule, la zone géographique de circulation, l’usage socioprofessionnel ou encore le kilométrage parcouru….
A partir de cette prime de référence, on calcule une prime qui tient compte de l’expérience passée, via un coefficient multiplicatif.
La valeur de référence est de 1 pour le taux de majoration-minoration. Ce taux subie une réduction – corrrespondant à un bonus – pour une année sans accident, d’un niveau de 5 %. Toutefois, la valeur du coefficient est arrêté à la deuxième décimale et arrondi par défaut. Au final, le taux prend donc un nombre fini de valeurs. Pour résumer schématiquement

1 année sans accident	100 x 0.95 = 0.95 de Bonus
2 années sans accident	0.95 x 0.95 = 0.90 de Bonus
3 années sans accident	0.90 x 0.95 = 0.85 de Bonus
4 années sans accident	0.85 x 0.95 = 0.80 de Bonus
5 années sans accident	0.80 x 0.95 = 0.76 de Bonus
6 années sans accident	0.76 x 0.95 = 0.72 de Bonus
7 années sans accident	0.72 x 0.95 = 0.68 de Bonus
8 années sans accident	0.68 x 0.95 = 0.64 de Bonus
9 années sans accident	0.64 x 0.95 = 0.60 de Bonus
10 années sans accident	0.60 x 0.95 = 0.57 de Bonus
11 années sans accident	0.57 x 0.95 = 0.54 de Bonus
12 années sans accident	0.54 x 0.95 = 0.51 de Bonus
13 années sans accident	0.51 x 0.95 = 0.48 plafonné à 0.50

La réduction maximale est donc de 0.50, soit 50% de réduction sur la prime de référence.
En revanche, chaque accident responsable donne lieu à une majoration – on palera de malus – de 25 % du coefficient par accident. Toutefois, cette majoration est réduite de moitié (12,5 %) si la responsabilité n’est que partiellement engagée.Là encore, le coefficient obtenu est arrêté à la deuxième décimale et arrondi par défaut. Aussi, avec deux accidents responsables dans l’année, le coefficient sera multiplié par un coefficient de 1,56. Afin d’éviter l’absence d’assurance, on notera qu’en aucun cas, le coefficient de réduction-majoration ne peut excéder 3,50.
Pour les jeunes conducteurs, l’article A. 335-9-1 du code des assurances prévoit que l’assureur a la possibilité de majorer la prime de référence (dans des limites fixées par la loi là encore).
Notons enfin qu’il existe des avantages pour les bons conducteurs: si un conducteur est au bonus maximal depuis au moins trois ans, il conserve son bonus après le premier accident responsable. Notons que c’est le genre de disposition qui figure dans la loi et que certains assureurs utilisent à des fins commerciales (comme le montre l’affiche ci-dessus ou ci-dessous, qui a été analysée sur le blog de Gizmo (ici)).

Il existe aussi une règle dite de descente rapide qui peut atténuer les effets des malus: après deux années consécutives sans accident, votre coefficient ne peut être supérieur à 1.
Ce qui rend le système intéressant, c’est qu’il est obligatoire, en France. Aussi, en cas de changement d’assureur, le coefficient de bonus-malus sera recalculé à partir d’un relevé d’information: l’assureur est tenu de fournir ce document à la résiliation du contrat (contenant des informations sur le conducteur et le véhicule, ainsi que le nombre des sinistres survenus au cours des cinq années précédant l’établissement du relevé d’information, ainsi que leur nature et la part de responsabilité retenue)
Davantage d’information se trouve ici, là (ce document propose d’ailleurs une petite simulation à la fin sur le coût d’un accident) ou là encore (pour une approche plus actuarielle), voire là (pour une approche plus économique).

Actuarial science, GLM

De l’interprétation d’un effet nonlinéaire en régression

01/10/2009 Arthur Charpentier Leave a comment

Je vais poursuivre un peu le TD de statistique de l’actuariat 2, où nous avions parlé tarification et régression Poissonnienne.

Analyser un effet nonlinéaire

Dans un premier temps, nous avions regardé un modèle simple, où l’on voulait voir si l’âge du conducteur expliquait la fréquence de sinistres. Bref, on peut utiliser une régression log-Poisson, qui donne la prédiction suivante, en fonction de l’âge,

Mais l’âge est une variable un peu particulière car c’est une fausse variable continue: seule les valeurs entières apparaissent. Comme on chercher à calculer

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-glm-nb-1.png$

notons qu’il est possible de considérer l’estimateur empirique naturel de cette grandeur, i.e.

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-glm-nb-2.png$

où

$https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-glm-nb-3.png$

Graphiquement, on obtient la courbe suivante,

Le côté ératique de la série à droite et à à gauche vient de la faible représentation dans la base, comme en atteste l’histogramme de l’âge du conducteur,

Je peux en profiter pour revenir à un commentaire posté ici sur les modèles GAM. En effet, au lieu de regarder la fréquence en fonction de l’âge, on peut s’intéresser à la la différence par rapport à la fréquence globale.

> m=mean(base$nbsin)
> Y=predict(REG,newdata=data.frame(age=X,expo=E),type = "response")
> lines(X,Y/m-1,col="red")

On obtient alors le graphique suivant (en superposant à l’ajustement obtenu par GAM)

Notons que si l’on ajuste un modèle log-Poisson uniquement sur la tranche d’âge [50-70], on est très proche de ce que propose le modèle GAM, d’où l’interprétation de l’ajustement GAM comme un ajustement local.

Changer de modèle (et les comparer)

On a choisi un modèle log-Poisson, mais on se doute que d’autres modèles seraient possible. En particulier dans notre base, les assurés ont eu 0 ou 1 sinistre l’année passée. Autrement dit, des lois binomiales peuvent être ajustée,

> REG0 = glm(nbsin~age,base,
+  family=poisson(link="log"),
+  offset=expo)
> REG1a = glm(nbsin~age,base,
+  family=binomial(link = "logit"),
+  offset=expo)
> REG1b = glm(nbsin~age,base,
+  family=binomial(link = "probit"),
+  offset=expo)
> REG1c = glm(nbsin~age,base, 
+  family=gaussian(link = "identity"),
+  offset=expo)

Le critère AIC est souvent invoqué pour juger de l’adéquation d’un modèle. Il vaut ici

> AIC(REG0)
[1] 36470.7
> AIC(REG1a)
[1] 36044.45
> AIC(REG1b)
[1] 37304.02
> AIC(REG1c)
[1] 32549.83

Mais on peut également regarder le critère BIC,

> AIC(REG0, k = log(nrow(base)))
[1] 36488.81
> AIC(REG1a, k = log(nrow(base)))
[1] 36062.57
> AIC(REG1b, k = log(nrow(base)))
[1] 37322.14
> AIC(REG1c, k = log(nrow(base)))
[1] 32577

Autrement dit, en utilisant le règle “the smaller, the better” on retiendrait le dernier modèle qui est un modèle Gaussien, manifestement assez mauvais: sur le graphique ci-dessous, la régression log-Poisson est en mauve, presque confondu avec la régression logistique¹, en bleu, la régression probit en vert et le modèle linéaire (Gaussien) en rouge,

Bref, ce dernier modèle n’est pas convainquant… Par contre, notons que si l’on compare le modèle log-Poisson et logit, ce dernier conduit à tarifer moins cher les âges extrêmes,

Autrement dit, au lieu de comparer des modèles actuariels sur la base d’outils statistiques, et uniquement de ces outils, il peut être instructif d’essayer de comprendre qui verra sa prime baisser avec tel ou tel modèle… Mais je reviendrais plus longuement, une autre fois, sur les histoires de choix de modèles…..

Actuarial science, Climate, Ratemaking

Quand les assureurs souhaiteraient faire la pluie et le beau temps

23/07/2009 Arthur Charpentier Leave a comment

Damned, je me suis fais devancé par un très très beau billet d’Anne (ici) sur l’assurance proposée par AXA et AON. Ca m’apprendra à ne plus lire régulièrement Le Monde

Tout comme Anne, je suis un peu agacé par le nom donné sur le site d’AON: un assureur ne peut pas garantir la pluie et le beau temps, de même qu’il ne peut pas garantir que vous ne tomberez pas malade ou que vous n’aurez pas d’accident de voiture dans l’année ! Ou alors c’est autre chose que de l’assurance…

Je suis également un peu sceptique car ces risques ne sont pas vraiment indépendants, surtout au niveau local (dans un même département, il y a des chances pour que tout le monde soit sinistré la même semaine). Peut être qu’au niveau national, il existe des dipôles climatiques qui font que s’il pleut dans un département, il y a de très fortes chances pour qu’il ne pleuve pas dans un autre. C’était une étude qu’on souhaitait mener avec des collègues de Paris 6, et nous n’avions riens trouvé sur le sujet. Mais les gens d’AXA et de Metnext savent ce qu’ils font, je pense… Car il faut reconnaître que les garanties (ici) excluent beaucoup de choses si on regarde en détails le contrat, en particulier quand on lit qu’il faut “3 heures consécutives” de pluie, “entre 8 heures et 20 heures“…etc.
En tous les cas AON a l’honnêteté d’annoncer clairement la couleur sur son site (ici). Tout d’abord avec de joli petit logo que j’ai reproduit ci dessous, avec cette superbe accroche (évoquée en détails ici) cette garantie est, pour les voyagistes, “une source de revenus complémentaires“…
Sinon il faut reconnaître également que les assureurs offraient déjà ce genre de garanties, mais seulement dans les grands évènements, comme des championnats du monde de football, des tournois de tennis, ou des concerts de grosses stars….pour la première fois, la garantie n’est plus pour Michael Jackson qui souhaiterait remplir des stades, mais pour Robert et Raymonde qui voudraient faire du camping au sec….

Actuarial science, Bayesian, Publications, Ratemaking

Article (de vulgarisation) sur la crédibilité

18/10/2008 Arthur Charpentier 1 Comment

Parution d’un court article sur la Théorie de la Crédibilité dans la revue Risques.

La théorie de la crédibilité est une jolie application de l’approche bayésienne de la statistique, puisque l’on cherche à déterminer la probabilité d’appartenir à une classe de risque conditionnellement à une expérience de sinistralité observée. La théorie développée par Hans Bühlmann propose d’ailleurs de jolies applications, en particulier sur les jeux de fléchettes, où l’on compare la moyenne sur un sous-groupe homogène (ou supposé l’être) avec la moyenne de la population totale. En faisant intervenir des variance intra et des variances inter, on arrive à un résultat relativement simple et élégant.

Freakonometrics