Cette semaine, nous allons finir la modélisation du passif d’une compagnie d’assurance non-vie. Les slides sont en ligne.
Actuariat de l’Assurance Non-Vie #10
La semaine dernière, nous avions commencé la modélisation du passif d’une compagnie d’assurance non-vie. Les slides introductifs sont en ligne.
Actuariat de l’Assurance Non-Vie #9
Cette semaine, nous avons fini les modèles de tarification, avec une extention sur les modèles globaux (sans passer par une séparation entre fréquence et coût moyen), et sur le choix de variables, et le choix de modèles. Les slides introductifs sont en ligne.
Faire (rapidement) un zonier
Dans le cours d’actuariat de l’assurance non-vie, on avait évoqué rapidement l’idée de faire un zonier. Autrement dit, on souhaite créer une variable polytomique (avec 4 ou 5 classes) représentant le critère spatial du risque. On va tenter de segmenter le territoire, en regroupant l’ensemble des territoires ‘proches’ (en terme de risque) étudiées en un nombre de classes prédéfini, permettant de tenir compte d’un (potentiel) critère spatial dans la tarification. Pour les aspects pratiques, comme on va manipuler des données spatiales, chargeons quelques packages pour commencer
library(maptools) library(rgeos) library(rgdal) library(ggplot2) library(plyr) library(maptools) library(cartography)
On va aussi avoir besoin de données spatialisées (disons en France métropolitaire). Pour faire simple, je vais tirer des codes insee de communes au hasard, un millier, et une variable Y qui va nous intéresser. Et je suppose qu’elle dépend de X que je n’observe pas, mais qui est liée à des caractéristiques spatiales (en gros la latitude, i.e. un positionnement nord/sud)
download.file("http://freakonometrics.free.fr/popfr19752010.csv","popfr.csv")
base = read.csv("popfr.csv",header=TRUE)
base$insee = base$dep*1000+base$com
n=1000
set.seed(123)
id=sample(1:nrow(base),size=n)
simbase=data.frame(insee=base$insee[id])
X=(46-base$lat[id])/2
simbase$Y=X+rnorm(n)
Ici, on a la base suivante (mais on va supposer que je n’observe pas X – qui n’a servi qu’à simuler des données)
> head(simbase) insee Y 1 61499 -1.8573363 2 19181 -0.7059191 3 55307 -0.5923649 4 74030 0.6098773 5 30328 -0.4584795 6 10050 -1.2361240
Classiquement, Y peut être un résidu (normalisé) de régression, par exemple. Pour visualiser, il faut un fond de carte
download.file("http://biogeo.ucdavis.edu/data/gadm2.8/rds/FRA_adm0.rds","FRA_adm0.rds")
FR0=readRDS("FRA_adm0.rds")
Ensuite, on met les points (en fusionnant nos données avec la base insee donnant latitude et longitude des codes insee de communes)
plot(FR0)
simbase = merge(simbase,base[,c("insee","long","lat")])
cols = rev(carto.pal(pal1 = "red.pal",n1 = 10,pal2="green.pal",n2=10))
bk = seq(-5,4.5,length=21)
cuty = cut(simbase$Y,breaks=bk,labels=1:20)
points(simbase$long,simbase$lat, col=cols[cuty],pch=19,cex=.5)
On va voir eux techniques pour faire un zonier. Le premier est de travailler par zone prédéfinie, comme le département,
simbase$dpt=trunc(simbase$insee/1000)
A=aggregate(x = simbase$Y,by=list(simbase$dpt),mean)
names(A)=c("dpt","y")
download.file("http://biogeo.ucdavis.edu/data/gadm2.8/rds/FRA_adm2.rds","FRA_adm2.rds")
FR=readRDS("FRA_adm2.rds")
donnees_carte=data.frame(FR@data)
d=donnees_carte$CCA_2
d[d=="2A"]="201"
d[d=="2B"]="202"
donnees_carte$dpt=as.numeric(as.character(d))
donnees_carte=merge(donnees_carte,A,all.x=TRUE)
donnees_carte=donnees_carte[order(donnees_carte$OBJECTID),]
bk=seq(-2.75,2.75,length=21)
donnees_carte$cuty=cut(donnees_carte$y,breaks=bk,labels=1:20)
cols = rev(carto.pal(pal1 = "red.pal",n1 = 10,pal2="green.pal",n2=10))
plot(FR, col=cols[donnees_carte$cuty],xlim=c(-5.2,12))
autrement dit, on fait une carte choroplèthe
(il faut s’assurer que les couleurs sont mises au bon endroit). Ensuite, on peut définir deux zones,
bk=seq(-2.75,2.75,length=3) donnees_carte$cuty=cut(donnees_carte$y,breaks=bk,labels=1:2) plot(FR, col=cols[c(4,16)][donnees_carte$cuty],xlim=c(-5.2,12))
ou quatre
bk=seq(-2.75,2.75,length=5) donnees_carte$cuty=cut(donnees_carte$y,breaks=bk,labels=1:4) plot(FR, col=cols[c(3,8,12,17)][donnees_carte$cuty],xlim=c(-5.2,12))
On va créer uen variable à 4 modalités, à partir des départements.
Une autre solution consiste à utiliser les données spatiales, obtenues par fusion avec la base insee
simbase=merge(simbase,base[,c("insee","long","lat")])
On va ensuite se donner une grille, en France. C’est un peu pénible. Il nous faut le polynôme de contour de la France
P1=FR0@polygons[[1]]@Polygons[[355]]@coords P2=FR0@polygons[[1]]@Polygons[[27]]@coords plot(FR0,border=NA) polygon(P1) polygon(P2)
(ici on a juste la France métropolitaine, et la Corse – deux polygônes donc)
et on part d’une grille sur un rectangle
grille<-expand.grid(seq(min(simbase$long),max(simbase$long),length=101),seq(min(simbase$lat),max(simbase$lat),length=101)) paslong=(max(simbase$long)-min(simbase$long))/100 paslat=(max(simbase$lat)-min(simbase$lat))/100
On retient alors juste les points qui sont dans les polygônes
f=function(i){ (point.in.polygon (grille[i, 1]+paslong/2 , grille[i, 2]+paslat/2 , P1[,1],P1[,2])>0)+(point.in.polygon (grille[i, 1]+paslong/2 , grille[i, 2]+paslat/2 , P2[,1],P2[,2])>0) }
indic=unlist(lapply(1:nrow(grille),f))
grille=grille[which(indic==1),]
points(grille[,1]+paslong/2,grille[,2]+paslat/2,cex=.4,pch=19,col="blue")
voilà le résultat
Maintenant, on peut utiliser du krigeage mais j’ai plutôt voulu tenter des plus proches voisins. Pour chaque point de la grille, on prend la moyenne des plus proches voisins (à vol d’oiseau, i.e. avec une norme Euclidienne – sur la sphère)
library(geosphere)
knn=function(i,k=20){
d=distHaversine(grille[i,1:2],simbase[,c("long","lat")], r=6378.137)
r=rank(d)
ind=which(r<=k)
mean(simbase[ind,"Y"])
}
grille$y=Vectorize(knn)(1:nrow(grille))
bk=seq(-2.75,2.75,length=21)
grille$cuty=cut(grille$y,breaks=bk,labels=1:20)
cols <- rev(carto.pal(pal1 = "red.pal",n1 = 10,pal2="green.pal",n2=10))
points(grille[,1]+paslong/2,grille[,2]+paslat/2,col=cols[grille$cuty],pch=19)
Ici, ca donne la carte suivante
mais là encore, on peut retenir deux niveaux
bk=seq(-2.75,2.75,length=3) grille$cuty=cut(grille$y,breaks=bk,labels=1:2) plot(FR0,border=NA) polygon(P1) polygon(P2) points(grille[,1]+paslong/2,grille[,2]+paslat/2,col=cols[c(4,16)][grille$cuty],pch=19)
ce qui donne les zones suivantes
bk=seq(-2.75,2.75,length=5) grille$cuty=cut(grille$y,breaks=bk,labels=1:4) plot(FR0,border=NA) polygon(P1) polygon(P2) points(grille[,1]+paslong/2,grille[,2]+paslat/2,col=cols[c(3,8,12,17)][grille$cuty],pch=19)
ou bien, si on retient quatre niveaux
A partir de là, pour créer une variable de zone, le plus simple est à un point donné de chercher le point de la grille le plus proche, et de lui donner la valeur de la couleur associée
pred=function(z){
d=distHaversine(z,grille[,1:2], r=6378.137)
grille[which.min(d),"cuty"]
}
Facile, non ?
Données sinistres, devoir Actuariat (ENSAE)
Pour le dernier devoir du cours d’actuariat non-vie, à l’ENSAE, quelques données de sinistres. Pour plusieurs milliers de sinistres (garantie dommages tous accidents issu d’un portefeuille automobile d’une société française) on a les données suivantes, entre 1995 et 2014,
- id_sin: identifiant du sinistres (les 4 premiers caractères représentent l’exercice de survenance)
- an_surv: exercice de survenance
- an_gest: exercice de gestion
- etat_sin: code état du sinistre à la fin de l’exercice de gestion (1=mis en jeu, 2= cloture partielle, 3 = cloture totale, 4 = réouvert et 5 = cloture sans suite)
- pmt: paiements effectués sur le sinistre à la fin de l’exercice de gestion
- rec: recours encaissés sur le sinistre à la fin de l’exercice de gestion
- eval_pmt: évaluation des paiements (réglés+à venir) vue à la fin de l’exercice de gestion
- eval_rec: évaluation des recours (réglés+à venir) vue à la fin de l’exercice de gestion
Pour information complémentaire, on dispose des éléments d’information suivants, sur le portefeuille
| Exercice | Volume portefeuille (garanties années) |
Chiffre d’affaires (€) | Nb sinistres survenus dans l’année |
| 2014 | 303 900 | 45 900 000 | 41 678 |
| 2013 | 296 500 | 45 700 000 | 41 639 |
| 2012 | 295 600 | 44 000 000 | 42 815 |
| 2011 | 292 000 | 41 200 000 | 43 351 |
| 2010 | 282 900 | 40 700 000 | 45 276 |
| 2009 | 274 300 | 39 800 000 | 42 779 |
| 2008 | 263 800 | 41 900 000 | 38 941 |
| 2007 | 256 200 | 41 200 000 | 38 048 |
| 2006 | 250 800 | 42 400 000 | 38 603 |
| 2005 | 239 100 | 41 100 000 | 37 539 |
| 2004 | 226 800 | 37 900 000 | 34 728 |
| 2003 | 220 300 | 39 200 000 | 34 491 |
| 2002 | 211 400 | 35 500 000 | 35 984 |
| 2001 | 202 400 | 34 800 000 | 36 322 |
| 2000 | 193 400 | 34 000 000 | 36 186 |
| 1999 | 184 500 | 32 300 000 | 34 123 |
| 1998 | 175 500 | 29 800 000 | 30 517 |
| 1997 | 166 500 | 28 000 000 | 28 432 |
| 1996 | 157 500 | 26 500 000 | 26 852 |
| 1995 | 148 600 | 25 700 000 | 26 775 |
La base complète est dans le fichier ensae2017.zip. Nous reviendrons sur les attendus lors du prochain cours.
Actuariat de l’Assurance Non-Vie #8
Pour poursuivre sur la modèlisation des coûts individuels, quelques slides sur les grands risques. Les slides sont en ligne.
Actuariat de l’Assurance Non-Vie #7
Après pas mal de temps passé sur les modèles de comptage, on va commencer la modèlisation des coûts individuels. Les slides sont en ligne.
Actuariat de l’Assurance Non-Vie #6
On reviendra sur les modèles de comptage, après la digression sur les GLM et la sur-dispersion, en parlant un peu de tarification a posteriori. En particulier, nous verrons deux approches assez différentes, avec les modèles de crédibilité et l’approche par bonus-malus. Les slides sont en ligne.
Actuariat de l’Assurance Non-Vie #5
Petit retour sur les modèles de comptage, avec la digression sur les GLM, avec les modèles permettant d’avoir de la surdispersion (théoriquement impossible dans un modèle de Poisson, mais plus réaliste compte tenu des particularités des données d’assurance). Les slides sont en ligne.
Actuariat de l’Assurance Non-Vie #4
Ce mardi, poursuite du cours d’actuariat, avec les modèles linéaires généralisés (qui englobent la régression logistique et la régression de Poisson, sur lesquelles nous avons passé du temps). Les slides sont en ligne.
Multinomial Logit as an Iterated Logit Regression
For the second section of the course at ENSAE, yesterday, we’ve seen how to run a multinomial logistic regression model. It is simply an extension of the binomial logistic regression. But actually, it is also possible to consider iterative binomial regressions.
Consider here a response variable Y with a multinomial distribution (3 factors to have something more general than the binomial), taking values \{A,B,C\}, with respective probabilities \mathbf{p}=(p_A,p_B,p_C). Here is a code to generate some multinomial variables
msample=function(A,B,C){
Y=rep(NA,B)
for(i in 1:B){Y[i]=sample(A,size=1,prob=C[i,])}
return(Y)
}
and here is a code to generate a dataset with n rows,
generate3=function(n,x,pb=c(-2,0)){
set.seed(x)
X1=runif(n)
X2=runif(n)
X3=runif(n)
s1=pb[1]+X1+X2
s2=pb[2]-X1+X2
P1=exp(s1)/(1+exp(s1)+exp(s2))
P2=exp(s2)/(1+exp(s1)+exp(s2))
Y=msample(0:2,n,cbind(1-P1-P2,P1,P2))
df=data.frame(Y=Y,X1=X1,X2=X2,X3=X3)
return(df)
}
Let us generate a training dataset and a validation one
pb=c(.31,.42)
DF1=generate3(1000,1,pb=pb)
DF2=generate3(500,2,pb=pb)
With a multivariate logistic regression
\mathbb{P}[Y=A|\mathbf{x}]=\frac{\exp[\mathbf{x}^{\text{T}}\mathbf{\alpha}]}{1+\exp[\mathbf{x}^{\text{T}}\mathbf{\alpha}]+\exp[\mathbf{x}^{\text{T}}\mathbf{\beta}]}
\mathbb{P}[Y=B|\mathbf{x}]=\frac{\exp[\mathbf{x}^{\text{T}}\mathbf{\beta}]}{1+\exp[\mathbf{x}^{\text{T}}\mathbf{\alpha}]+\exp[\mathbf{x}^{\text{T}}\mathbf{\beta}]}
\mathbb{P}[Y=B|\mathbf{x}]=\frac{1}{1+\exp[\mathbf{x}^{\text{T}}\mathbf{\alpha}]+\exp[\mathbf{x}^{\text{T}}\mathbf{\beta}]}
For convenience, consider the most popular factor in our training dataset
modalite=names(sort(table(DF1$Y),decreasing = TRUE))
Consider a regression model on the simulated dataset (with several covariates), let us estimate it, and let us get predictions.
library(nnet)
reg=multinom(as.factor(Y) ~ ., data = DF1)
mp1=predict (reg, DF1, "probs")
mp2=predict (reg, DF2, "probs")
An alternative can be the following.
consider a first regression model on the Bernoulli variable Y_A=\mathbf{1}(Y=A). Actually, we will consider the most important factor, but for convenience, assume that it is A.
\mathbb{P}[Y_A=A|\mathbf{x}]=\frac{\exp[\mathbf{x}^{\text{T}}\mathbf{a}]}{1+\exp[\mathbf{x}^{\text{T}}\mathbf{a}]}
On our dataset, estimate that model, and get predictions. In the case where Y\neq A, define another Bernoulli variable Y_B=\mathbf{1}(Y=B|Y\neq A). We can estimate that model and derive two probabilities, \mathbb{P}(Y=B|Y\neq A) and \mathbb{P}(Y=C|Y\neq A) (the sum of the two being equal to 1). Based on those two models, it is possible to compute the three probabilities we are looking for. \mathbb{P}[Y=A] is obtained from the first model, and we can derive the other two from \mathbb{P}[Y=B|Y\neq A]\cdot\mathbb{P}[Y\neq A] and \mathbb{P}[Y=C|Y\neq A]\cdot\mathbb{P}[Y\neq A].
reg1=glm((Y==modalite[1])~.,data=DF1,family=binomial)
reg2=glm((Y==modalite[2])~.,data=DF1[-which(DF1$Y==modalite[1]),],family=binomial)
p11=predict (reg1, newdata=DF1, type="response")
p12=predict (reg2, newdata=DF1, type="response")
p21=predict (reg1, newdata=DF2, type="response")
p22=predict (reg2, newdata=DF2, type="response")
mmp1=cbind(p11,(1-p11)*p12,(1-p11)*(1-p12))
mmp2=cbind(p21,(1-p21)*p22,(1-p21)*(1-p22))
colnames(mmp1)=colnames(mmp2)=modalite
Let us compare the predicted probabilites, on the same dataset (here the training dataset)
> mmp1[1:9,c("0","1","2")]
0 1 2
1 0.19728737 0.4991805 0.3035321
2 0.17244580 0.5648537 0.2627005
3 0.19291753 0.5971058 0.2099767
4 0.09087176 0.7787304 0.1303978
5 0.23400225 0.4083022 0.3576955
6 0.18063647 0.6637352 0.1556283
7 0.13188881 0.7402710 0.1278401
8 0.13776970 0.6524959 0.2097344
9 0.12325864 0.6790336 0.1977078
> mp1[1:9,c("0","1","2")]
0 1 2
1 0.19691036 0.5022692 0.3008205
2 0.17123189 0.5680647 0.2607034
3 0.19293066 0.5984402 0.2086291
4 0.08821851 0.7813318 0.1304497
5 0.23470739 0.4109990 0.3542936
6 0.18249687 0.6602168 0.1572863
7 0.13128711 0.7400898 0.1286231
8 0.13525341 0.6553618 0.2093848
9 0.12090016 0.6815915 0.1975084
The two are very close. So yes, it is possible to see the multinomial regression as some sequential binomial regressions.
Variables Catégorielles et Modèle Logistique
Petit complément, suite aux coquilles qu’il y avait dans les slides, sur une propriété de la régression logistique quand on régresse sur des variables catégorielles. On était sur la base des avocats
> avocat <- read.table("http://freakonometrics.free.fr/AutoBI.csv",header=TRUE,sep=",") > avocat$CLMSEX <- factor(avocat$CLMSEX, labels=c("M","F")) > avocat$MARITAL <- factor(avocat$MARITAL, labels=c("M","C","V","D")) > avocat=avocat[!is.na(avocat$CLMSEX),]
> attach(avocat)
> sum((ATTORNEY==2)&(CLMSEX=="F"),na.rm=TRUE)/sum(CLMSEX=="F",na.rm=TRUE)
[1] 0.5256065
> sum((ATTORNEY==2)&(CLMSEX=="M"),na.rm=TRUE)/sum(CLMSEX=="M",na.rm=TRUE)
[1] 0.4453925
Autrement dit, la proportion d’hommes qui se sont fait représentés par un avocat est de 44.539%, et la proportion de femmes 52.56%. On peut visualiser le tableau croisé ci-dessous
> tab=xtabs(~ATTORNEY+CLMSEX,data=avocat)
> require(vcd)
> mosaic(tab, shade=TRUE, legend=TRUE)
Quand on fait une régression linéaire (Gaussienne), on retrouve ces probabilités dans les valeurs de coefficients: pour les femmes, on a 44.539%, et pour les hommes, on a la différence avec les femmes (cette dernière modalité étant la modalité de référence ici)
> reglm = lm((ATTORNEY==2) ~ CLMSEX, data=avocat)
> summary(reglm)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.44539 0.02060 21.620 < 2e-16 ***
CLMSEXF 0.08021 0.02756 2.911 0.00367 **
C’est ce que l’on a ci-dessous
> sum(coefficients(reglm))
[1] 0.5256065
> coefficients(reglm)[1]
(Intercept)
0.4453925
On a un résultat similaire avec une régression logistique
> reglogit = glm((ATTORNEY==2) ~ CLMSEX, data=avocat,family=binomial)
> summary(reglogit)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.21930 0.08312 -2.639 0.00833 **
CLMSEXF 0.32182 0.11097 2.900 0.00373 **
même si les coefficients n’ont pas la même interprétation. En transformant ces coefficients, on retrouve très exactement les mêmes valeurs
> exp(sum(reglogit$coefficients[1:2])) /(1+exp(sum(reglogit$coefficients[1:2])))
[1] 0.5256065
> exp(reglogit$coefficients[1]) /(1+exp(reglogit$coefficients[1]))
(Intercept)
0.4453925
Désolé pour la typo.
Actuariat de l’Assurance Non-Vie #3
Mardi, seconde partie du cours d’actuariat, avec les modèles de classification le matin, mais l’après midi, on devrait commencer les modèles de comptage. Les slides sont en ligne.
Le matin, je dois intervenir 15 minutes à l’IHP vers 9 heures, dans un colloque sur Artificial Intelligence for Fintech and Insurtech: si le RER traîne un peu trop, entre Luxembourg et Lozère, j’aurais peut-être 5 minutes de retard….
Base de Données, Devoir d’Actuariat
Pour les deux premiers devoirs du cours d’actuariat de l’assurance non-vie, on utilisera les mêmes données, dont le descriptif complet est en ligne ici. La première base d’apprentissage (sa version corrigée) est aussi en ligne
> base=read.csv2("http://freakonometrics.free.fr/PG3/PG_2017_SURV_YEAR0.csv",sep=",")
> dim(base)
[1] 100000 32
> table(base$survenance)
FALSE TRUE
85636 14364
On reviendra en cours mardi prochain sur les attendus, mais l’idée sera de modélisation la variable survenance. Pour se second projet, on passera à la modélisation de la prime pure.
Actuariat de l’Assurance Non-Vie #2
Mardi, on commencera le cours d’actuariat, avec un cours d’introduction le matin, mais l’après midi, on attaquera plus sérieusement, avec une mise en bouche, avec les modèles de classification (de type {survenance} vs {non survenance} d’accident). Là encore, les slides sont en ligne.
Je parlerais aussi du premier devoir, qui portera sur les modèles que nous (re)verrons ici.