# Base de données de tarification

Pour compléter le cours de ce matin, un mot rapide sur les bases, et plus particulièrement la base de contrats. Au sujet des variables,

• densite est la densité de population dans la commune où habite le conducteur principal,
• zone : zone A B C D E ou F, selon la densité en nombre d’habitants par km2 de la commune de résidence (A =”1-50″, B=”50-100″, C=”100-500″, D=”500-2,000″, E=”2,000-10,000″, F=”10,000+”.

A titre d’information, la répartition de la population en France se fait de la manière suivante • marque : marque du véhicule selon la table suivante (1 Renault Nissan; 2 Peugeot Citroën ; 3 Volkswagen Audi Skoda Seat ; 4 Opel GM; 5 Ford ; 6 Fiat ; 10 Mercedes Chrysler ; 11 BMW Mini ;12 Autres japonaises et coréennes ; 13 Autres européennes ; 14 Autres marques et marques inconnues). Cette variable n’est pas une variable numérique
• region : code à 2 chiffres (ce qui n’est pas une valeur numérique) donnant les 22 régions françaises (code INSEE), soit géographiquement • ageconducteur : âge du conducteur principal en début de la couverture,
• agevehicule : âge du véhicule en début de période.

Je demande de ne pas utiliser la variable de bonus, qui fait intervenir une information utilisée en tarification a posteriori (qui ne fait pas l’objet de ce cours).

# Regression on categorical variables

This morning, Stéphane asked me tricky question about extracting coefficients from a regression with categorical explanatory variates. More precisely, he asked me if it was possible to store the coefficients in a nice table, with information on the variable and the modality (those two information being in two different columns). Here is some code I did to produce the table he was looking for, but I guess that some (much) smarter techniques can be used (comments – see below – are open). Consider the following dataset

> base
x sex   hair
1  1   H  Black
2  4   F  Brown
3  6   F  Black
4  6   H  Black
5 10   H  Brown
6  5   H Blonde

with two factors,

> levels(base$hair)  "Black" "Blonde" "Brown" > levels(base$sex)
 "F" "H"

Let us run a (standard linear) regression,

> reg=lm(x~hair+sex,data=base)

which is here

> summary(reg)

Call:
lm(formula = x ~ hair + sex, data = base)

Residuals:
1          2          3          4          5          6
-3.714e+00 -2.429e+00  2.429e+00  1.286e+00  2.429e+00 -2.220e-16

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)   3.5714     3.4405   1.038    0.408
hairBlonde    0.2857     4.8655   0.059    0.959
hairBrown     2.8571     3.7688   0.758    0.528
sexH          1.1429     3.7688   0.303    0.790

Residual standard error: 4.071 on 2 degrees of freedom
Multiple R-squared: 0.2352,	Adjusted R-squared: -0.9121
F-statistic: 0.205 on 3 and 2 DF,  p-value: 0.886

If we want to extract the names of the factors (assuming here that there are no numbers in the name of the factor), and the values of the associated modality, one can use

> VARIABLE=c("",gsub("[-^0-9]", "", names(unlist(reg$xlevels)))) > MODALITY=c("",as.character(unlist(reg$xlevels)))
> names=data.frame(VARIABLE,MODALITY,NOMVAR=c(
+ "(Intercept)",paste(VARIABLE,MODALITY,sep="")[-1]))
> regression=data.frame(NOMVAR=names(coefficients(reg)),
+ COEF=as.numeric(coefficients(reg)))
> merge(names,regression,all.x=TRUE)
NOMVAR VARIABLE MODALITE      COEF
1 (Intercept)                   3.5714286
2   hairBlack     hair    Black        NA
3  hairBlonde     hair   Blonde 0.2857143
4   hairBrown     hair    Brown 2.8571429
5        sexF      sex        F        NA
6        sexH      sex        H 1.1428571

or, if we want modalities exluding references,

> merge(names,regression)
NOMVAR VARIABLE MODALITE      COEF
1 (Intercept)                   3.5714286
2  hairBlonde     hair   Blonde 0.2857143
3   hairBrown     hair    Brown 2.8571429
4        sexH      sex        H 1.1428571

In order to reproduce the table Stéphane sent me, let us use the following code to produce an html table,

> library(xtable)
> htlmtable <- xtable(merge(names,regression))
> print(htlmtable,type="html")
NOMVAR VARIABLE MODALITY COEF
1 (Intercept) 3.57
2 hairBlonde hair Blonde 0.29
3 hairBrown hair Brown 2.86
4 sexH sex H 1.14

So yes, it is possible to build a table with the variable, modalities, and coefficients. This function can be interesting on prospective mortality, when we do have a large number of modalities per factor (years, ages and year of birth). Consider the following datasets

> DEATH=read.table(
+ "http://freakonometrics.free.fr/DeathsSwitzerland.txt",
+ "http://freakonometrics.free.fr/ExposuresSwitzerland.txt",
> DEATH$Age=as.numeric(as.character(DEATH$Age))
> DEATH=DEATH[-which(is.na(DEATH$Age)),] > EXPOSURE$Age=as.numeric(as.character(EXPOSURE$Age)) > EXPOSURE=EXPOSURE[-which(is.na(EXPOSURE$Age)),]
> base=data.frame(y=as.factor(DEATH$Year),a=as.factor(DEATH$Age),
+ c=as.factor(DEATH$Year-DEATH$Age),D=DEATH$Total,E= EXPOSURE$Total)
> base=base[base$E>0,] and the following nonlinear model, based on Lee-Carter model (including a cohort effect), can be estimated using > library(gnm) > reg=gnm(D~a+Mult(a,y)+Mult(a,c),offset=log(E),family=poisson,data=base) In order to extract the 671 coefficients from the regresssion, > length(coefficients(reg))  671 (as properly as possible) we have to be careful: names of coefficients are not that simple to handle. For instance, we can see things like > coefficients(reg) Mult(., year).age98 0.04203519 In order to extract them, define > na=length((reg$xlevels)$age) > ny=length((reg$xlevels)$year) > nc=length((reg$xlevels)$cohort) > VARIABLElong=c("",rep("age",na),rep("Mult(., year).age",na), + rep("Mult(a, .).y",ny), + rep("Mult(., cohort).age",na),rep("Mult(age, .).cohort",nc)) > VARIABLEshort=c("",rep("age",na),rep("age",na),rep("year",ny), + rep("age",na),rep("cohort",nc)) > MODALITY=c("",(reg$xlevels)$age,(reg$xlevels)$age, + (reg$xlevels)$year,(reg$xlevels)$age,(reg$xlevels)$cohort) > names=data.frame(VARIABLElong,VARIABLEshort, + MODALITY,NOMVAR=c("(Intercept)",paste(VARIABLElong,MODALITY,sep="")[-1])) > regression=data.frame(NOMVAR=names(coefficients(reg)), + COEF=as.numeric(coefficients(reg))) Here we go, now we have the coefficients from the regression in a nice table, > outputreg=merge(names,regression) > outputreg[1:10,] NOMVAR VARIABLElong VARIABLEshort MODALITY COEF 1 (Intercept) -8.22225458 2 age1 age age 1 -0.87495451 3 age10 age age 10 -1.67145704 4 age100 age age 100 4.91041650 5 age11 age age 11 -1.00186990 6 age12 age age 12 -1.05953497 7 age13 age age 13 -0.90952859 8 age14 age age 14 0.02880668 9 age15 age age 15 0.42830738 10 age16 age age 16 1.35961403 It is now possible to plot all the coefficients, as functions of the age, the year of observation, or the year of birth. For instance, for the standard average age effect (namely as a function of ), we can use > typevariable=as.character(unique(outputreg$VARIABLElong))
> basegraph=outputreg[outputreg$VARIABLElong==typevariable,] > x=as.numeric(as.character(basegraph$MODALITY))
> y=basegraph$COEF > plot(x,y,type="p",col="blue",xlab="Age") while the cohort effect ( as a function of ) is obtained using > basegraph=outputreg[outputreg$VARIABLElong==typevariable,]
> x=as.numeric(as.character(basegraph$MODALITY)) > y=basegraph$COEF
> plot(x,y,type="p",col="blue",xlab="Cohort (year of birth)",ylim=c(0,10)) # Introduction aux modèles linéaires généralisés

J’ai un peu d’avance dans le cours. Je vais mettre en ligne les transparents pour la semaine prochaine (normalement), où nous aborderons la classe des modèles linéaires généralisés. Les transparents sont en ligne ici.

Je n’ai pas mis de section sur lesGeneralized Additive Models, on se contentera de la section sur le lissage évoquée à la fin des transparents sur la modélisation de la fréquence. Afin de légitimer les méthodes de lissage (sur l’âge de l’assuré en particulier), je renvoie vers un graphique produit il y a plusieurs années par un cabinet de conseil, qui notait que la forme de la fonction de lissage, liant l’âge à la fréquence est identique, dans tous les pays, Mais je pense que je ferais un billet dédié au lissage, dans la problématique de la tarification en assurance IARD.

# The law of small numbers

In insurance, the law of large numbers (named loi des grands nombres initially by Siméon Poisson, see e.g. http://en.wikipedia.org/…) is usually mentioned to legitimate large portfolios, because of pooling and diversification: the larger the pool, the more ‘predictable’ the losses will be (in a given period). Of course, under standard statistical assumption, namely finite expected value, and independence (see http://freakonometrics.blog.free.fr/…. for a discussion, in French). Since in insurance, catastrophes are usually rare – and extremely costly – and actuaries might be interested to model occurrence of that small number of events (see e.g. Aldous’ book on that specific topic, that can be downloaded from http://stat.berkeley.edu/…). The theorem behind is sometimes called the law of small numbers (from the book published by Ladislaus Bortkiewicz, but we’ll get back to that story later on, see also Whitaker (1914) http://biomet.oxfordjournals.org/… or the book recently published by Michael Falk, Jürg Hüsler and Rolf-Dieter Reiss).

• The Poisson distribution

The so-called Poisson distribution (see http://en.wikipedia.org/…) was introduced by Siméon Poisson in 1837 (in Recherches sur la Probabilité des Jugements en Matière Criminelle et en Matière Civile, Précédées des Règles Générales du Calcul des Probabilités, see http://gallica.bnf.fr/…). But it had been defined more than a century before, by Abraham De Moivre, in 1711, in De Mensura Sortis seu; de Probabilitate Eventuum in Ludis a Casu Fortuito Pendentibus (see e.g. the review in http://www.jstor.org/…). Let  denote a counting random variable, then it said to be Poisson distributed if there is  such that

De Moivre obtained that distribution from an approximation of the binomial distribution. Recall that the binomial distribution is a standard distribution in actuarial science, for instance to model the number of deaths among  insured. If individual death probabilities are identical, say , and if deaths are independent events, then

And if  and  , then

Again, this is an asymptotic theorem, which is valid when we have a lot of observations (), but also, the probability of occurrence should be extremely small (since ), which is why to use the term small numbers. Siméon Poisson was not interested by mathematical approximations: his main point was to get a distribution with nice goodness of fit properties for the data he was working on. He wanted to get a better understanding of cours d’assises (jury panel, might be a valid translation of the French term). A jury consists of 12 jurors who voted to determine whether a defendant was guilty. When guilt was predominant, with at least 8 votes against 4, the defendant was convicted (which was 47% of criminal cases). 5 with 7 votes against, the opinion of professional judges was requested (11% of criminal trials again). Using these statistics we can demonstrate that a defendant brought before an assize court is guilty of the order of 68%, and the probability that a juror is not wrong by voting (condemning an innocent or releasing a culprit) was about 54%. He sought to calculate the probability that a defendant is wrongfully convicted, and gets 2%. And 28% of exonerated defendants are in fact guilty. Siméon Poisson introduced this law to get probabilities easily. But the law he considered is central in probability….

• The law of small numbers

The heuristic of the main theorem, related to the Poisson distribution is the following: let  denote i.i.d random variables taking values in (in a general setting, one component can be the time, the other one an upper region of interest, where some stochastic process might be). Let . If  as (or to be a little bit more specific about the assumptions), let  denote the (random variable characterizing) count of events , then  can be approximated by a Poisson distribution with parameter .
The heuristic is that if we consider a large number of observations, and if we count how many are in a given (small) region, then the number of such observations is Poisson distributed.

n=1000
X=runif(n)*10-1.5
Y=runif(n)*10-1.5
plot(X,Y,axis=FALSE,cex=.6)
u=seq(-1,1,by=.01)
v=sqrt(1-u^2)
polygon(c(u,rev(u)),c(v,rev(-v)),col="yellow",border=NA)
I=(X^2+Y^2)<1
points(X[I],Y[I],cex=.6,pch=19,col="red") If we run some simulations,

>  n=1000
>  ns=100000
>  N=rep(NA,ns)
> for(s in 1:ns){
+ X=runif(n)*10-1.5
+ Y=runif(n)*10-1.5
+ I=(X^2+Y^2)<1
+ N[s]=sum(I)
+ }
> hist(N,breaks=0:60,probability=TRUE,col="yellow")
> mean(N)
 31.41257

The parameter of the Poisson distribution is the area of the yellow disk, over the area of the square, i.e.

> (lambda=10*pi)
 31.41593
> lines(0:60-.5,dpois(0:60,lambda),type="b",col="red") To get an interpretation related to insurance modeling, let  denote an upper layer in a reinsurance contract, i.e. for some deductible . Let ‘s denote individual losses. Then the number of claims that hit this upper layer can be modeled with a Poisson distribution. More precisely, if deductible  becomes extremely large (and ), we obtain the point-over-threshold model in extreme value theory (see e.g. http://brale.math.hr/~iugrina/… or  http://fire.nist.gov/bfrlpubs/…): if  has a Poisson distribution and, conditionally on are independent identically distributed generalized Pareto random variables, then  has the generalized extreme value distribution. Thus, exceedances models (for rare events) are closely related to Poisson processes.

• The Poisson process

As mentioned above, the Poisson distribution appears when events occur somehow randomly and independently, over time. It is then natural to study the time between two occurences (or two claims, in an insurance context).

• Poisson distribution, and claims occurrence

It is neither Siméon Poisson nor De Moivre, but Ladislaus Von Bortkiewicz who first mentioned the Poisson distribution as the law of small numbers. In 1898 (see https://archive.org/…), he studied the number number of soldiers killed by being kicked by a horse, from 1875 till 1894, in 200 corps (more precisely 10 corps over 20 ans).

He did obtain the following distribution (here, the parameter of the Poisson distribution is 0.61, i.e. the average number of death per year)

 number of death per year Empirical counts Poisson distribution 0 109 108.67 1 65 66.21 2 22 20.22 3 3 4.11 4 1 0.63 5 and more 0 0.08

It is possible to find a lot of cases where the Poisson distribution fits extremely well. For instance, if we consider the number of hurricanes, that landed in Florida after 1850,

 number of hurricanes per year empirical frequency Poisson frequency 0 30 27.16 1 48 47.99 2 37 42.41 3 29 24.98 4 8 11.03 5 3 3.90 6 3 1.15 7 1 0.29 8 and more 0 0.08
• Poisson distribution, and return period

The return period was introduced by Emil Gumbel, in hydrology, to link probabilities and durations (see e.g. http://freakonometrics.blog.free.fr/…). A decennial event has an occurence probability of 1/10. 10 is then the average waiting time before occurence. This does not mean that the event will not occur before 10 years, or has to occur before 10 years. Consider a return period  (in years), then the yearly probability of non-occurrence is .

And the probability of non-occurence over  years is then . It is standard to summarize this property with the following table,

 return period Number of years () without catastrophes 10 20 50 100 200 10 65.1% 40.1% 18.3% 9.6% 4.9% 20 87.8% 64.2% 33.2% 18.2% 9.5% 50 99.5% 92.3% 63.6% 39.5% 22.5% 100 99.9% 99.4% 86.7% 63.4% 39.5% 200 99.9% 99.9% 98.2% 86.6% 63.3%

The diagonal in the table above is extremely interesting. It looks like there is some kind of convergence towards a limiting value (here  63.2%). Indeed, the number of events observed over n years have a Binomial distribution, with probability , which will converge towards the Poisson distribution with parameter 1. The probability of not having a catastrophe is then , which is equal to 0.632.

• Rare probabilities and the Poisson distribution

The Poisson distribution keeps appearing when computing probabilies of rare events. For instance, the probability to have at least one incident in a nuclear plant in France, over a 50 year period. Assume that the annual probability of an incident in a reactor  is small, e.g. 0.05%. Assume further that reactors are independent among them, and in time. The probability to have an incident over 80 reactors in 50 years is (exactly)

Of course, a linear approximation is not correct (even if it was mentioned in some French newspaper, as explained in an old post http://freakonometrics.blog.free.fr/…)

On the other hand

> p=0.00005
> 1-(1-p)^(50*80)
 0.1812733
> 1-exp(-50*80*p)
 0.1812692

which is the probability that  is null when  has a Poisson distribution with parameter . We clearly see here an application of De Moivre’s approximation in risk management.

Another way of looking at this problem is based on the following idea: given the fact that in 45 years of observations on 450 reactors worldwide (roughly), three major accidents were observed including Three Mile Island (1979) and Fukushima (2011), i.e. the average time between accidents can be estimated at 16 years. For a single reactor, we can assume that the average time to wait before an incident is 450 times 16 years, i.e 7200 years. Or the probability to have one incident, over one year, for one reactor is 1 over 7200 (this is the idea behind the return period concept). If we assume that the arrival of accidents occurs randomly and independently of each other (as defined above) then the number of major accidents observed over a period of 50 years in France follows a Poisson distribution with parameter 50 / (7200/80). Also, the probability of having no major accident over 50 years, with 80 reactors can be estimated by

i.e.

> 1-exp(-50*80/7200)
 0.4262466

(keeping in mind all the uncertainty around the estimated waiting time before a major accident to a single reactor!).

# Regression tree using Gini’s index

In order to illustrate the construction of regression tree (using the CART methodology), consider the following simulated dataset,

> set.seed(1)
> n=200
> X1=runif(n)
> X2=runif(n)
> P=.8*(X1<.3)*(X2<.5)+
+   .2*(X1<.3)*(X2>.5)+
+   .8*(X1>.3)*(X1<.85)*(X2<.3)+
+   .2*(X1>.3)*(X1<.85)*(X2>.3)+
+   .8*(X1>.85)*(X2<.7)+
+   .2*(X1>.85)*(X2>.7)
> Y=rbinom(n,size=1,P)
> B=data.frame(Y,X1,X2)

with one dichotomos varible (the variable of interest, ), and two continuous ones (the explanatory ones  and ).

> tail(B)
Y        X1        X2
195 0 0.2832325 0.1548510
196 0 0.5905732 0.3483021
197 0 0.1103606 0.6598210
198 0 0.8405070 0.3117724
199 0 0.3179637 0.3515734
200 1 0.7828513 0.1478457

The theoretical partition is the following Here, the sample can be plotted below (be careful, the first variate is on the y-axis above, and the x-axis below) with blue dots when  equals one, and red dots when  is null,

> plot(X1,X2,col="white")
> points(X1[Y=="1"],X2[Y=="1"],col="blue",pch=19)
> points(X1[Y=="0"],X2[Y=="0"],col="red",pch=19) In order to construct the tree, we need a partition critera. The most standard one is probably Gini’s index, which can be writen, when ‘s are splited in two classes, denoted here or when ‘s are splited in three classes, denoted etc. Here,  are just counts of observations that belong to partition  such that  takes value . But it is possible to consider other criteria, such as the chi-square distance, where, classically when we consider two classes (one knot) or, in the case of three classes (two knots) Here again, the idea is to maximize that distance: the idea is to discriminate, so we want samples as not independent as possible. To compute Gini’s index consider

> GINI=function(y,i){
+ T=table(y,i)
+ nx=apply(T,2,sum)
+ pxy=T/matrix(rep(nx,each=2),2,ncol(T))
+ vxy=pxy*(1-pxy)
+ zx=apply(vxy,2,sum)
+ n=sum(T)
+ -sum(nx/n*zx)
+ }

We simply construct the contingency table, and then, compute the quantity given above. Assume, first, that there is only one explanatory variable. We split the sample in two, with all possible spliting values , i.e.

Then, we compute Gini’s index, for all those values. The knot is the value that maximizes Gini’s index. Once we have our first knot, we keep it (call it, from now on ). And we reiterate, by seeking the best second choice: given one knot, consider the value that splits the sample in three, and give the highest Gini’s index, Thus, we consider either the following partition

or this one

I.e. we cut either below, or above the previous knot. And we iterate. The code can be something like that,

> X=X2
> u=(sort(X)[2:n]+sort(X)[1:(n-1)])/2
> knot=NULL
> for(s in 1:4){
+ vgini=rep(NA,length(u))
+ for(i in 1:length(u)){
+ kn=c(knot,u[i])
+ F=function(x){sum(x<=kn)}
+ I=Vectorize(F)(X)
+ vgini[i]=GINI(Y,I)
+ }
+ plot(u,vgini)
+ k=which.max(vgini)
+ cat("knot",k,u[k],"\n")
+ knot=c(knot,u[k])
+ u=u[-k]
+ }
knot 69 0.3025479
knot 133 0.5846202
knot 72 0.3148172
knot 111 0.4811517

At the first step, the value of Gini’s index was the following, which was maximal around 0.3. Then, this value is considered as fixed. And we try to construct a partition in three parts (spliting either below or above 0.3). We get the following plot for Gini’s index (as a function of this second knot) which is maximum when the split the sample around 0.6 (which becomes our second knot). Etc. Now, let us compare our code with the standard R function,

> tree(Y~X2,method="gini")
node), split, n, deviance, yval
* denotes terminal node

1) root 200 49.8800 0.4750
2) X2 < 0.302548 69 12.8100 0.7536 *
3) X2 > 0.302548 131 28.8900 0.3282
6) X2 < 0.58462 65 16.1500 0.4615
12) X2 < 0.324591 7  0.8571 0.1429 *
13) X2 > 0.324591 58 14.5000 0.5000 *
7) X2 > 0.58462 66 10.4400 0.1970 *

We do obtain similar knots: the first one is 0.302 and the second one 0.584. So, constructing tree is not that difficult…

Now, what if we consider our two explanatory variables? The story remains the same, except that the partition is now a bit more complex to write. To find the first knot, we consider all values on the two components, and again, keep the one that maximizes Gini’s index,

> n=nrow(B)
> u1=(sort(X1)[2:n]+sort(X1)[1:(n-1)])/2
> u2=(sort(X2)[2:n]+sort(X2)[1:(n-1)])/2
> gini=matrix(NA,nrow(B)-1,2)
> for(i in 1:length(u1)){
+ I=(X1<u1[i])
+ gini[i,1]=GINI(Y,I)
+ I=(X2<u2[i])
+ gini[i,2]=GINI(Y,I)
+ }
> mg=max(gini)
> i=1+sum(mg==max(gini[,2]))
> par(mfrow = c(1, 2))
> plot(u1,gini[,1],ylim=range(gini),col="green",type="b",xlab="X1",ylab="Gini index")
> abline(h=mg,lty=2,col="red")
> if(i==1){points(u1[which.max(gini[,1])],mg,pch=19,col="red")
+          segments(u1[which.max(gini[,1])],mg,u1[which.max(gini[,1])],-100000)}
> plot(u2,gini[,2],ylim=range(gini),col="green",type="b",xlab="X2",ylab="Gini index")
> abline(h=mg,lty=2,col="red")
> if(i==2){points(u2[which.max(gini[,2])],mg,pch=19,col="red")
+          segments(u2[which.max(gini[,2])],mg,u2[which.max(gini[,2])],-100000)}
> u2[which.max(gini[,2])]
 0.3025479

The graphs are the following: either we split on the first component (and we obtain the partition on the right, below), or we split on the second one (and we get the following partition), Here, it is optimal to split on the second variate, first. And actually, we get back to the one-dimensional case discussed previously: as expected, it is optimal to split around 0.3. This is confirmed with the code below,

> library(tree)
> arbre=tree(Y~X1+X2,data=B,method="gini")
> arbre$frame[1:4,] var n dev yval splits.cutleft splits.cutright 1 X2 200 49.875000 0.4750000 <0.302548 >0.302548 2 X1 69 12.811594 0.7536232 <0.800113 >0.800113 4 <leaf> 57 8.877193 0.8070175 5 <leaf> 12 3.000000 0.5000000 For the second knot, four cases should be considered: spliting on the second variable (again), either above, or below the previous knot (see below on the left) or spliting on the first one. Then whe have wither a partition below or above the previous knot (see below on the right), Etc. To visualize the tree, the code is the following > plot(arbre) > text(arbre) > partition.tree(arbre) Note that we can also visualize the partition. Nice, isn’t it? To go further, the book Classification and Regression Trees by Leo Breiman (and co-authors) is awesome. Note that there are also interesting sections in the bible Elements of Statistical Learning: Data Mining, Inference, and Prediction by Trevor Hastie, Robert Tibshirani and Jerome Friedman (which can be downloaded from http://www.stanford.edu/~hastie/…) # Régression de Poisson, et biais minimal Lors du prochain cours d’actuariat, nous allons finir les arbres de régression, et introduire la régression de Poisson. Les transparents sont en ligne ici, Je vais présenter la régression en Poisson, en faisant un parallèle avec la régression logistique, la session suivante portera sur la généralisation obtenue avec les modèles linéaires généralisés. Sur la régression de Poisson, je suggère de lire Frees (2010) chapitre 12 (p 343-361), Greene (2012), section 18.3 (p 802-828) ou encore de Jong Heller (2008) chapitre 6. Sur les méthodes de biais minimal, de Jong Heller (2008), section 1.3 et l’article de Sholom Feldblum, http://www.casact.org/…. Sur le passage de ces dernières méthodes (introduites par Robert Bailey dans les années 60, http://www.casact.org/… et http://www.casact.org/…), je recommande la lecture de l’article de Ben Zehnwirth, Ratemaking From Bailey and Simon (1960) to Generalized Linear Regression Models, en ligne sur http://www.casact.org/… Comme annoncé au premier cours, j’essaye de mettre en ligne les transparents au fur et à mesure, mais j’avais pris l’habitude d’écrire au tableau ces dernières années. Il faut donc que je tape tout. Pour le devoir un courriel sera envoyé d’ici la fin de semaine à tous les groupes qui se sont inscrits. # Words are important (when dealing with statistics). But still. In statistics, it might be difficult to know what a symbol stands for. For instance, can either be a real value, i.e. the value taken by a statistics from a given sample. But it can also be a random variable, assuming that the sample is now a collection of i.i.d. random variables. We can usually distinguish ‘s (values from a given sample) and ‘s (the underlying random variables, i.e. for some ). But notations might confusing, and it is hard to distinguish random variables, and values taken by random variables (or realizations). But usually, if we look at the context, one can figure out what symbols stand for. But sometimes, it is difficult to get a proper definition, not for some symbols, but for words. And most of the time common words. Recently, I wrote a short paper, claiming that it was difficult to model the number of bodily injuries related to car accident, since it is difficult to define death. Actually, the definition of dead did change a few years ago (as weird as it might sound), which did cause a rupture of some series. I recently had a similar story, discussing with a pharmacist in Montréal who said to me “you French are known to be the world’s champion in terms of drug consumption“, see e.g. • The French are Europe’s champion medicine-takers” in http://economist.com/…, mentioning “heavy drug-consumption culture • The data show that drug consumption in France remains one of the largest in Europe” in http://bizcovering.com/… • France has one of the largest drug markets in the world and the drug consumption per capitahttp://ispor.org/… (among so many articles) I do not think I am a drug addict (I might be – like most of my colleagues – a coffee addict, but as Paul Erdős – or more probably Alfréd Rényi – said once, “a mathematician is a device for turning coffee into theorems“). The main problem here is the notion of “consumption“. The economics interpretation is simply that someone buys a product or a service (see http://dictionary.reference.com/…). There is also the food-related interpretation, where consuming means ingesting, i.e. eating or drinking (see http://dictionary.cambridge.org/…). So pill and drug “consumption” is ambiguous: is it the number of pills purchased, or ingested (actually consumed), or prescribed? The first thing one should remember is that the Social Security in France refunds (almost) all medications prescribed by a doctor. So it is uncommon to leave the office of a doctor without a prescription, at least of aspirin: a visit to the doctor is usually, in France, the opportunity to stock some over-the-counter drugs. The second thing is that there is a major difference between France and North America when we go to the pharmacy. In Montréal for instance, if I have a prescription for 12 pills, then the pharmacist does give me exactly 12 pills (from a big pot). In France, pills are sold in prepackaged boxes, so if the box contains 10 pills, I will get 2 boxes, just to be sure I’ll get my 12 pills. From a medical point of view, I will consume my 12 pills, but from an economic perspective, I will consume 20. So comparing statistics is extremely difficult, not because the the maths, but because it is difficult to define (even simple) concepts. # Inglés ? Französisch ? When I started to blog, I found it “natural” to blog in French. Not “natural” in the way Dominique Bouhours sees it: according to Bouhours, only the French language exactly reflects the natural way of thinking La langue française est peut-être la seule qui suive exactement l’ordre naturel, et qui exprime les pensées en la manière qu’elles naissent dans l’esprit as he wrote it). “Natural” in the sense that French is my mother language. And perhaps also because, as a French, I did not really see the interest of using another language… Yes, somehow, I still feel that we (French people) believe that Antoine de Rivarol was right when he said, in his Discours sur l’universalité de la langue française that French is innately superior to other languages: French has an element of probity attached to its genius, it is not the language of the French, it is the language of humanity Dégagée de tous les protocoles que la bassesse invente pour la vanité et le pouvoir, elle en est plus faite pour la conversation, lien des hommes et charme de tous les âges, et puisqu’il faut le dire, elle est de toutes les langues la seule qui ait une probité attachée à son génie. Sûre, sociale, raisonnable, ce n’est plus la langue française, c’est la langue humaine But I might be exaggerating here: if I can express more subtle things in French, it simply means that my English is too poor. I mean, in my work as a professor and a researcher, I do use English on a daily basis (reading textbooks, discussing with colleagues, or co-autors, writing research papers, etc). But blogging is just for fun… so at first, I started to blog in French. Then, while I was giving a short course in Brazil, I understood that some students there actually did go to my blog. They were using online translation – which translated quite poorly by that time – and some asked me if I could write more posts in English. I did have also a request from Tal, the editor of http://r-bloggers.com/ who asked me to share my posts that contain R codes with the R community (i.e. non French-speaking people). So I decided to write more and more posts in English. Except perhaps posts related to my undergraduate courses. I thought about all that recently because, while I was in Europe last week, some colleagues in Amsterdam mentioned my blog, claiming that there were (still) too many posts in French. While on the other hand, some French colleagues still ask me why I do not write only in French. Somehow, I still find it more comfortable to write in French, especially on French related issues (e.g. politics and pools, or university related debates). But I still like the idea that my post may be interesting for students in Brazil, actuaries in the U.K. or researchers in Norway… So I guess I will try to write more posts in English. Which is not that complicated actually, since, as a French academician (or George Bernard Shaw, the source is not clear) once put it “English is a language that is relatively easy to speak poorly“. So, let’s do it…. # Régression logistique et arbres Pour le cours de mercredi prochain, la base utilisée sera une base tirée du livre de Jed Frees, http://instruction.bus.wisc.edu/jfrees/… > baseavocat=read.table("http://freakonometrics.free.fr/AutoBI.csv",header=TRUE,sep=",") > tail(baseavocat) CASENUM ATTORNEY CLMSEX MARITAL CLMINSUR SEATBELT CLMAGE LOSS 1335 34204 2 2 2 2 1 26 0.161 1336 34210 2 1 2 2 1 NA 0.576 1337 34220 1 2 1 2 1 46 3.705 1338 34223 2 2 1 2 1 39 0.099 1339 34245 1 2 2 1 1 18 3.277 1340 34253 2 2 2 2 1 30 0.688 On dispose d’une variable dichotomique indiquant si un assuré – suite à un accident de la route – a été représenté par un avocat (1 si oui, 2 si non). On connaît le sexe de l’assuré (1 pour les hommes et 2 pour les femmes), le statut marital (1 s’il est marié, 2 s’il est célibataire, 3 pour un veuf, et 4 pour un assuré divorcé). On sait aussi si l’assuré portait ou non une ceinture de sécurité lorsque l’accident s’est produit (1 si oui, 2 si non et 3 si l’information n’est pas connue). Enfin, une information pour savoir si le conducteur du véhicule était ou non assuré (1 si oui, 2 si non et 3 si l’information n’est pas connue). On va recoder un peu les données afin de les rendre plus claires à lire. Les transparents du cours sont en ligne ici, Sur les arbres de régression, je mettrais en ligne un billet, afin d’illustrer la méthode. En attendant des compléments théoriques peuvent se trouver en ligne http://genome.jouy.inra.fr/…, http://ensmp.fr/…, ou http://ujf-grenoble.fr/… (pour information, nous ne verrons que la méthode CART). Je peux renvoyer au livre (et au blog) de Stéphane Tuffery, ou (en anglais) au livre de Richard Berk, dont un résumé se trouve en ligne sur http://crim.upenn.edu/…. La semaine suivante, nous aborderons la régression de Poisson, les méthodes de biais minimal, et introduirons les modèles linéaires généralisés. Je renvoie au chapitre sur la tarification a priori du Denuit & Charpentier (2005), aux chapitres 12 et 13 de Frees (2010) ou encore les chapitres 5 et 6 du De Jong & Heller (2008). Pour les plus curieux qui veulent comprendre les liens entre les modèles linéaires généralisés et la tarification par crédibilité, je renvoie à l’article de Klinker (2010) # Amsterdam, PhD defense Last week, I was involved in the PhD defense of Julien Tomas, with Rob Kaas (promotor), Frédéric Planchet (co-promotor), Katrien AntonioMarc Goovaerts, Ann De Schepper and Michel Vellekoop. The PhD thesis – untitled quantifying biometric life insurance risks with non-parametric smoothing methods – can be dowloaded on http://dare.uva.nl/… and on http://tel.archives-ouvertes.fr/. The R codes will be available soon on my blog (and on Julien’s new website http://www.likelihood.me/).   # Amsterdam I will be in Amsterdam for the end of this week. I will be in the jury of the PhD defense of Julien Tomas, entitled “Quantifying Biometric Life Insurance Risks With Non-Parametric Smoothing Methods” (the thesis will probably be online soon). But before, I will give a talk at the actuarial seminar at UvA. My visit last time was a real pleasure, and it should be the same this time too. I will give a talk this Thursday on “R for actuarial science“. The slides can be downloaded from here. # Lausanne I will be back in Lausanne (I was already there last summer) to spend a few days, visiting Florian, at HEC Lausanne. I will also give a talk on old and new results on (standard) families of copulas. There will be also a discussion on tail dependence. Slides can be downloaded from here, # R for actuarial science As mentioned in the Appendix of Modern Actuarial Risk Theory, “R (and S) is the ‘lingua franca’ of data analysis and statistical computing, used in academia, climate research, computer science, bioinformatics, pharmaceutical industry, customer analytics, data mining, finance and by some insurers. Apart from being stable, fast, always up-to-date and very versatile, the chief advantage of R is that it is available to everyone free of charge. It has extensive and powerful graphics abilities, and is developing rapidly, being the statistical tool of choice in many academic environments. R is based on the S statistical programming language developed by Joe Chambers at Bell labs in the 80’s. To be more specific, R is an open-source implementation of the S language, developed by Robert Gentlemn and Ross Ihaka. It is a vector based language, which makes it extremely interesting for actuarial computations. For instance, consider some Life Tables, > TD[39:52,] > TV[39:52,] Age Lx Age Lx 39 38 95237 38 97753 40 39 94997 39 97648 41 40 94746 40 97534 42 41 94476 41 97413 43 42 94182 42 97282 44 43 93868 43 97138 45 44 93515 44 96981 46 45 93133 45 96810 47 46 92727 46 96622 48 47 92295 47 96424 49 48 91833 48 96218 50 49 91332 49 95995 51 50 90778 50 95752 52 51 90171 51 95488 Those (French) Life Tables can be found here > TD <- read.table( + "https://perso.univ-rennes1.fr/arthur.charpentier/TD8890.csv",sep=";",header=TRUE) > TV <- read.table( + "https://perso.univ-rennes1.fr/arthur.charpentier/TV8890.csv",sep=";",header=TRUE) From those vectors, it is possible to construct the matrix of death probabilities, , using for instance > Lx <- TD$Lx
>  m <- length(Lx)
>  p <- matrix(0,m,m); d <- p
>  for(i in 1:(m-1)){
+  p[1:(m-i),i] <- Lx[1+(i+1):m]/Lx[i+1]
+  d[1:(m-i),i] <- (Lx[(1+i):(m)]-Lx[(1+i):(m)+1])/Lx[i+1]}
>  diag(d[(m-1):1,]) <- 0
>  diag(p[(m-1):1,]) <- 0
>  q <- 1-p

One can compute easily, e.g., the (curtate) expectation of life defined as

and one can compute the vector of life expectancy, at various ages , as

> life.exp = function(x){sum(p[1:nrow(p),x])}
> e = Vectorize(life.exp)(1:m)

An actually, any kind of actuarial quantity can be derived from those matrices. The expected present value (or actuarial value) of a temporary life annuity-due is, for instance,

The code to compute those functions is here

> for(j in 1:(m-1)){ adots[,j]<-cumsum(1/(1+i)^(0:(m-1))*c(1,p[1:(m-1),j])) }

or consider the expected present value of a term insurance

with the following code

> for(j in 1:(m-1)){ A[,j]<-cumsum(1/(1+i)^(1:m)*d[,j]) }

Some more details can be found in the first part of the notes of the crash courses of last summer, in Meielisalp. Vector – or matrices – are extremely convenient to work with, when dealing with life contingencies. It is also possible to model prospective mortality. Here, the mortality is not only function of the age , but also time ,

> t(DTF)[1:10,1:10]
1899  1900  1901  1902  1903  1904  1905  1906  1907  1908
0  64039 61635 56421 53321 52573 54947 50720 53734 47255 46997
1  12119 11293 10293 10616 10251 10514  9340 10262 10104  9517
2   6983  6091  5853  5734  5673  5494  5028  5232  4477  4094
3   4329  3953  3748  3654  3382  3283  3294  3262  2912  2721
4   3220  3063  2936  2710  2500  2360  2381  2505  2213  2078
5   2284  2149  2172  2020  1932  1770  1788  1782  1789  1751
6   1834  1836  1761  1651  1664  1433  1448  1517  1428  1328
7   1475  1534  1493  1420  1353  1228  1259  1250  1204  1108
8   1353  1358  1255  1229  1251  1169  1132  1134  1083   961
9   1175  1225  1154  1008  1089   981  1027  1025   957   885

Thus, we now have a force of mortality matrix , or surface It is also possible to use R packages to estimate a Lee-Carter model of the mortality rate,

> library(demography)
> MUH =matrix(DEATH$Male/EXPOSURE$Male,nL,nC)
> POPH=matrix(EXPOSURE$Male,nL,nC) > BASEH <- demogdata(data=MUH, pop=POPH, ages=AGE, years=YEAR, type="mortality", + label="France", name="Hommes", lambda=1) > RES=residuals(LCH,"pearson") One can easily study residuals, for instance as a function of the age, or a function of the year, Some more details can be found in the second part of the notes of the crash courses of last summer, in Meielisalp. R is also interesting because of its huge number of libraries, that can be used for predictive modeling. One can easily use smoothing functions in regression, or regression trees, > TREE = tree((nbr>0)~ageconducteur,data=sinistres,split="gini",mincut = 1) > age = data.frame(ageconducteur=18:90) > y1 = predict(TREE,age) > reg = glm((nbr>0)~bs(ageconducteur),data=sinistres,family="binomial") > y = predict(reg,age,type="response") Some practitioners might be scared because the legend claims that R is not as good as SAS to handle large databases. Actually, a lot of functions can be used to import datasets. The most convenient one is probably > baseCOUT = read.table("http://freakonometrics.free.fr/baseCOUT.csv", + sep=";",header=TRUE,encoding="latin1") > tail(baseCOUT,4) numeropol debut_pol fin_pol freq_paiement langue type_prof alimentation type_territoire 6512 87291 2002-10-16 2003-01-22 mensuel A Professeur Vegetarien Urbain 6513 87301 2002-10-01 2003-09-30 mensuel A Technicien Vegetarien Urbain 6514 87417 2002-10-24 2003-10-21 mensuel F Technicien Vegetalien Semi-urbain 6515 88128 2003-01-17 2004-01-16 mensuel F Avocat Vegetarien Semi-urbain utilisation presence_alarme marque_voiture sexe exposition age duree_permis age_vehicule i coutsin 6512 Travail-occasionnel oui FORD M 0.2684932 47 29 28 1 1274.5901 6513 Loisir oui HONDA M 0.9972603 44 24 25 1 278.0745 6514 Travail-occasionnel non VOLKSWAGEN F 0.9917808 23 3 11 1 403.1242 6515 Loisir non FIAT F 0.9972603 23 4 11 1 230.9565 But if the dataset is too large, it is also possible to specify which variables might be interesting, using > mycols = rep("NULL", 18) > mycols[c(1,4,5,12,13,14,18)] <- NA > baseCOUTsubC = read.table("http://freakonometrics.free.fr/baseCOUT.csv", + colClasses = mycols,sep=";",header=TRUE,encoding="latin1") > head(baseCOUTsubC,4) numeropol freq_paiement langue sexe exposition age coutsin 1 6 annuel A M 0.9945205 42 279.5839 2 27 mensuel F M 0.2438356 51 814.1677 3 27 mensuel F M 1.0000000 53 136.8634 4 76 mensuel F F 1.0000000 42 608.7267 It is also possible (before running a code on the entire dataset) to import only the first lines of the dataset. > baseCOUTsubCR = read.table("http://freakonometrics.free.fr/baseCOUT.csv", + colClasses = mycols,sep=";",header=TRUE,encoding="latin1",nrows=100) > tail(baseCOUTsubCR,4) numeropol freq_paiement langue sexe exposition age coutsin 97 1193 mensuel F F 0.9972603 55 265.0621 98 1204 mensuel F F 0.9972603 38 9547.7267 99 1231 mensuel F M 1.0000000 40 442.7267 100 1245 annuel F F 0.6767123 48 179.1925 It is also possible to import a zipped file. The file itself has a smaller size, and it can usually be imported faster. > import.zip = function(file){ + temp = tempfile() + download.file(file,temp); + read.table(unz(temp, "baseFREQ.csv"),sep=";",header=TRUE,encoding="latin1")} > system.time(import.zip("http://freakonometrics.free.fr/baseFREQ.csv.zip")) trying URL 'http://freakonometrics.free.fr/baseFREQ.csv.zip' Content type 'application/zip' length 692655 bytes (676 Kb) opened URL ================================================== downloaded 676 Kb user system elapsed 0.762 0.029 4.578 > system.time(read.table("http://freakonometrics.free.fr/baseFREQ.csv", + sep=";",header=TRUE,encoding="latin1")) user system elapsed 0.591 0.072 9.277 Finally, note that it is possible to import any kind of dataset, not only a text file. Even a Microsoft Excel folder. On a Windows computer, one can use SQL queries > sheet = "c:\\Documents and Settings\\user\\excelsheet.xls" > connection = odbcConnectExcel(sheet) > spreadsheet = sqlTables(connection) > query = paste("SELECT * FROM",spreadsheet$TABLE_NAME,sep=" ")
> result = sqlQuery(connection,query)

Then, once the dataset is imported, several functions can be used,

> cost = aggregate(coutsin~ AgeSex,mean, data=baseCOUT)
> frequency = merge(aggregate(nbsin~ AgeSex,sum, data=baseFREQ),
+ aggregate(exposition~ AgeSex,sum, data=baseFREQ))
> frequency$freq = frequency$nbsin/frequency$exposition > base.freq.cost = merge(frequency, cost) Finally, R is interesting for its graphical interface. “If you can picture it in your head, chances are good that you can make it work in R. R makes it easy to read data, generate lines and points, and place them where you want them. Its very flexible and super quick. When youve only got two or three hours until deadline, R can be brilliant” as said Amanda Cox, a graphics editor at the New York Times. “R is particularly valuable in deadline situations when data is scant and time is precious.”. Several cases were considered on the blog http ://chartsnthings.tumblr.com/…. First, we start with a simple graph, here State Government control in the US Then try to find a nice visual representation, e.g. And finally, you can just print it in your favorite newspaper, And you can get any kind of graphs, And not only about politics, Graphs are important. “Its not just about producing graphics for publication. Its about playing around and making a bunch of graphics that help you explore your data. This kind of graphical analysis is a really useful way to help you understand what you’re dealing with, because if you cant see it, you cant really understand it. But when you start graphing it out, you can really see what you’ve got” as said Peter Aldhous, San Francisco bureau chief of New Scientist magazine. Even for actuaries. “The commercial insurance underwriting process was rigorous but also quite subjective and based on intuition. R enables us to communicate our analytic results in appealing and innovative ways to non-technical audiences through rapid development lifecycles. R helps us show our clients how they can improve their processes and effectiveness by enabling our consultants to conduct analyses efficiently”, as explained by John Lucker, team of advanced analytics professionals at Deloitte Consulting Principal, in http://blog.revolutionanalytics.com/r-is-hot/. See also Andrew Gelman’s view, on graphs, http://www.stat.columbia.edu/… So yes, actuaries might be interested to use R for actuarial communication, as mentioned in http ://www.londonr.org/… The Actuarial Toolkit (see http ://www.actuaries.org.uk/…) stresses the interest of R, “The power of the language R lies with its functions for statistical modelling, data analysis and graphics ; its ability to read and write data from various data sources; as well as the opportunity to embed R in excel or other languages like VBA. In the way SAS is good for data manipulations, R is superior for modelling and graphical output“. From 2011, Asia Capital Reinsurance Group (ACR) uses R to Solve Big Data Challenges (see http ://www.reuters.com/…). And Lloyd’s uses motion charts created with R to provide analysis to investors (as discussed on http ://blog.revolutionanalytics.com/…) A lot of information can be found on http ://jeffreybreen.wordpress.com/… Markus Gesmann mentioned on his blog a lot of interesting graphs used for actuarial reporting, http ://lamages.blogspot.ca/… Further, R is free. Which can be compared with SAS,$6,000 per PC, or \$28,000 per processor on a server (as mentioned on http ://en.wikipedia.org/…)

It is also becoming more and more popular, as a programming language. As mentioned on this month Transparent Language Popularity (see http ://lang-index.sourceforge.net/), R is ranked 12. Far away after C or Java, but before Matlab (22) or SAS (27). On StackOverFlow (see http ://stackoverflow.com/) is also far being C++ (399,232 occurrences) or Java (348,418), but with 21,818 occurrences, it appears before Matlab (14,580) and SAS (899). As mentioned on http ://r4stats.com/articles/popularity/ R is becoming more and more popular, on listserv discussion traffic It is clearly the most popular software in data analysis, as mentioned by the Rexer Analytics survey, in 2009 What about actuaries ? In a survey (see http ://palisade.com/…), R was not extremely popular. If we consider only statistical softwares, SAS is still far ahead, among UK and CAS actuaries But, as mentioned by Mike King, Quantitative Analyst, Bank of America, “I cant think of any programming language that has such an incredible community of users. If you have a question, you can get it answered quickly by leaders in the field. That means very little downtime.” This was also mentioned by Glenn Meyers, in the Actuarial Review “The most powerful reason for using R is the community” (in http ://nytimes.com/…). For instance, http ://r-bloggers.com/ has contributions from more than 425 R users.

As said by Bo Cowgill, from Google “The best thing about R is that it was developed by statisticians. The worst thing about R is that it was developed by statisticians.

# Segmentation en tarification, compléments

Dans le premier cours d’actuariat IARD, nous avons vu l’importance de la ségmentation, et son implication sur le calcul des primes (passer d’une espérance mathématique à une espérance conditionnelle). Pour aller un peu plus loin, quelques compléments,

pour une lecture plus économique de la problématique de la segmentation en assurance

ou pour une lecture plus légale

Sinon, plusieurs articles de vulgarisation peuvent être lu sur internet,

La première démo aura lieu lundi, en salle informatique. Karim sera une introduction au langage R, à la manipulation des variables (qualitatives et quantitatives). Je mettrais en ligne les transparents en fin de semaine, et les codes seront mis en ligne dans le courant de la semaine prochaine.

Comme annoncé hier, il n’y aura pas cours mercredi prochain. Le mercredi suivant, nous verrons la modélisation des variables indicatrices, i.e. la régression logistique, et les arbres de régression. On supposera que le modèle linéaire aura été vu, je mets un lien vers les transparents du cours ACT6420 de la session passée, notes de cours transparents1 et transparents2. Il est aussi possible de relire Frees (2010), chapitres 3, 4, 5 et 6.

Pour commencer à pratiquer la régression logistique, on utilisera la petite fonction suivante

logit = function(formula, lien="logit", data=NULL) {
}