# Actuariat de l’Assurance non-Vie #3

Après plusieurs jours d’interruption, le cours d’actuariat non-vie reprendra ce lundi. Les slides sont en ligne. Il s’agit d’une introduction aux modèles de comptage.

# Econometrics vs. Machine Learning with Temporal Patterns

A few months ago, I did publish a (long) post entitled ‘some thoughts on economics, mathematics, econometrics, machine learning, etc‘. In that post, I was discussing possible differences between foundations of econometrics, and machine learning. I wanted to get back today on an important point, related to training/sampling datasets, when we have temporal data.

I was discussing this morning, with a student of the Data Science for Actuaries program, an interesting point related to claim frequency models, for insurance ratemaking. Since the goal is to predict claims frequency (to assess the level of the insurance premium), he suggested to use old data to train the model, and more recent one to test it. The problem is that the model did not incorporate any temporal pattern, and we got surprising results.

Consider here a simple dataset,

> set.seed(1)
> n=50000
> X1=runif(n)
> T=sample(2000:2015,size=n,replace=TRUE)
> L=exp(-3+X1-(T-2000)/20)
> E=rbeta(n,5,1)
> Y=rpois(n,L*E)
> B=data.frame(Y,X1,L,T,E)

Claims frequency is driven by a Poisson process, with one covariate, X1, and we assume that the intensity decreases (with an exponential rate). Consider here a standard linear regression, without any time effect

> reg=glm(Y~X1+offset(log(E)),data=B,
+ family=poisson)

We can also compute the empirical annualized claims frequency

> u=seq(0,1,by=.01)
> v=predict(reg,newdata=data.frame(X1=u,E=1))
> p=function(x){
+   B=B[abs(B$X1-x)<.1,] + sum(B$Y)/sum(B$E) + } > vp=Vectorize(p)(seq(.05,.95,by=.1)) and plot the two curves on the same graph, > plot(seq(.05,.95,by=.1),vp,type="b") > lines(u,exp(v),lty=2,col="red") This is what we usually do in econometrics. In machine learning, and more specifically to assess the quality of the model, and for model selection, it is common to split the dataset in two parts. A training sample, and a validation sample. Consider some randomized training/validation samples, then fit a model on the training sample, and finally use it to get a prediction, > idx=sample(1:nrow(B),size=nrow(B)*7/8) > B_a=B[idx,] > B_t=B[-idx,] > reg=glm(Y~X1+offset(log(E)),data=B_a, + family=poisson) > u=seq(0,1,by=.01) > v=predict(reg,newdata=data.frame(X1=u,E=1)) > p=function(x){ + B=B_a[abs(B_a$X1-x)<.1,]
+   sum(B$Y)/sum(B$E)
+ }
> vp_a=Vectorize(p)(seq(.05,.95,by=.1))
> plot(seq(.05,.95,by=.1),vp_a,col="blue")
> lines(u,exp(v),lty=2)
> p=function(x){
+   B=B_t[abs(B_t$X1-x)<.1,] + sum(B$Y)/sum(B$E) + } > vp_t=Vectorize(p)(seq(.05,.95,by=.1)) > lines(seq(.05,.95,by=.1),vp_t,col="red") The blue curve is the prediction on the training sample (as we usually do in econometrics), but then the red curve is the prediction on the testing sample. Here, volatility probably comes from the small size of the testing sample (1 observation out of 8). Now, what if we use the year as a splitting criteria : we fit a model on old years to fit a model, and we test it on recent years, > B_a=subset(B,T<2014) > B_t=subset(B,T>=2014) > reg=glm(Y~X1+offset(log(E)),data=B_a,family=poisson) > u=seq(0,1,by=.01) > v=predict(reg,newdata=data.frame(X1=u,E=1)) > p=function(x){ + B=B_a[abs(B_a$X1-x)<.1,]
+   sum(B$Y)/sum(B$E)
+ }
> vp_a=Vectorize(p)(seq(.05,.95,by=.1))
> plot(seq(.05,.95,by=.1),vp_a,col="blue")
> lines(u,exp(v),lty=2)
> p=function(x){
+   B=B_t[abs(B_t$X1-x)<.1,] + sum(B$Y)/sum(B$E) + } > vp_t=Vectorize(p)(seq(.05,.95,by=.1)) > lines(seq(.05,.95,by=.1),vp_t,col="red") Clearly, we miss something here… We were looking at such a graph this morning, and it took me some time to understand how training and validation samples were designed, and that there was a possible temporal effect (actually, this morning, it was based on a 3 year training sample, and a 1 year validation sample). Since there is a temporal pattern, let us capture it. As an econometrician, let me use a regression model > reg=glm(Y~X1+T+offset(log(E)),data=B, + family=poisson) > C=coefficients(reg) > u=seq(1999,2016,by=.1) > v=exp(-(u-2000)/20-3) > plot(2000:2015,exp(C[1]+C[3]*(2000:2015))) > lines(u,v,lty=2,col="red") (I focus only on the evolution of the temporal variate on that graph). Here, we use a linear model, but there are usually no reason to assume linearity. So we might consider splines > library(splines) > reg=glm(Y~X1+bs(T)+offset(log(E)), + data=B,family=poisson) > u=seq(1999,2016,by=.1) > v=exp(-(u-2000)/20-3) > v2=predict(reg,newdata=data.frame(X1=0, + T=2000:2015,E=1)) > plot(2000:2015,exp(v2),type="b") > lines(u,v,lty=2,col="red") But here again, why should we assume that there is an underlying smooth function? There might be some ruptures… So let us consider a regression on factors > reg=glm(Y~0+X1+as.factor(T)+offset(log(E)), + data=B,family=poisson) > C=coefficients(reg) > u=seq(1999,2016,by=.1) > v=exp(-(u-2000)/20-3) > plot(2000:2015,exp(C[2:17]),type="b") > lines(u,v,lty=2,col="red") An alternative might be to consider some more general model, like a regression tree > library(rpart) > reg=rpart(Y~X1+T+offset(log(E)),data=B, + method="poisson",cp=1e-4) > p=function(t){ + B=B[B$T==t,]
+   B$E=1 + mean(predict(reg,newdata=B)) + } > y_m=Vectorize(function(t) p(t))(2000:2015) > u=seq(1999,2016,by=.1) > v=exp(-(u-2000)/20-3+.5) > plot(2000:2015,y_m,ylim=c(.02,.085),type="b") > lines(u,v,lty=2,col="red") Here, it seems that something went wrong. I guess it’s coming from the exposure. So consider a simplier model, on the annualized frequency, and with weights that are related to the exposure > reg=rpart(Y/E~X1+T,data=B,weights=B$E,cp=1e-4)
> p=function(t){
+   B=B[B$T==t,] + B$E=1
+   mean(predict(reg,newdata=B))
+ }
> y_m=Vectorize(function(t) p(t))(2000:2015)
> u=seq(1999,2016,by=.1)
> v=exp(-(u-2000)/20-3+.5)
> plot(2000:2015,y_m,ylim=c(.02,.085),type="b")
> lines(u,v,lty=2,col="red")

That was for the econometrician perspective. With a machine learning perspective, consider a training sample (here based on old data) and a validation sample (based on more recent ones)

> B_a=subset(B,T<2014)
> B_t=subset(B,T>=2014)

If we consider a model, it is easy to get a prediction on recent years, even if the model was designed to model older ones,

> reg_a=glm(Y~X1+T+offset(log(E)),
+ data=B_a,family=poisson)
> C=coefficients(reg_a)
> u=seq(1999,2016,by=.1)
> v=exp(-(u-2000)/20-3)
> plot(2000:2015,exp(C[1]+C[3]*c(2000:2013,
+ NA,NA)),type="b")
> lines(u,v,lty=2,col="red")
> points(2014:2015,exp(C[1]+C[3]*2014:2015),
+ pch=19,col="blue")

But if we use years as factors, things are more complicated.

> reg_a=glm(Y~0+X1+as.factor(T)+offset(log(E)),
+ data=B_a,family=poisson)
> C=coefficients(reg_a)
> RMSE=function(A){
+   L=exp(C[1]*B_t$X1+ A[1]*(B_t$T==2014) + A[2]*(B_t$T==2015)) + Y_t=L*B_t$E
+   sum( (Y_t - B_t$Y )^2)} > i=optim(c(.4,.4),RMSE)$par
> plot(2000:2015,c(exp(C[2:15]),NA,NA),)
> u=seq(1999,2016,by=.1)
> v=exp(-(u-2000)/20-3)
> lines(u,v,lty=2,col="red")
> points(2014:2015,exp(i),pch=19,col="blue")

becase we need to get a prediction on levels that were not in our training sample. Here, we minimize the RMSE to quantify factor levels for recent years. And the output is not that bad.

So yes, it is possible to get a training dataset on older data, and test it on recent years. But one should be careful, and take into account, properly, temporal patterns.

# Simple Distributions for Mixtures?

The idea of GLMs is that given some covariates$X$$Y|X$ has a distribution in the exponential family (Gaussian, Poisson, Gamma, etc). But that does not mean that $Y$ has a similar distribution… so there is no reason to test for a Gamma model for $Y$ before running a Gamma regression, for instance. But are there cases where it might work? That the non-conditional distribution is the same (same family at least) than the conditional ones?

For instance, if $(X,Y)$ has a joint Gaussien distribution, then both marginals are Gaussian, but also $Y|X$. So, in that case, if the covariate is normally distributed, it is possible to have a Gaussian distribution also for $Y$. The econometric interpretation is that with a standard Gaussian linear model, if$X$ is normally distributed, not only the conditional distribution $Y|X$ is Gaussian but also the non-conditional distribution of $Y$.

> set.seed(1)
> n=1e3
> X=rnorm(n,10,2)
> Y=1+3*X+rnorm(n)
> plot(X,Y,xlim=c(4,20))

Indeed, here the distribution of $Y$ is also Gaussian

> library(nortest)

Anderson-Darling normality test

data:  Y
A = 0.23155, p-value = 0.802

> shapiro.test(Y)

Shapiro-Wilk normality test

data:  Y
W = 0.99892, p-value = 0.8293

(not only from a statistical point of view, the thoery of Gaussian random vectors confirms that the non-conditional distribution is Gaussian actually)

Here $X$ is continuous. What if we consider a finite mixture here, i.e.$X$ takes only a finite number of values? Actually, Teicher (1963) proved that it is not possible to have a non-conditional Gaussian distribution for $Y$. But in practice, would we really reject the Gaussian assumption, for $Y$? If the number of classes is to small, yes. But with a large number of classes (a sufficiently large number of mixture components), it is possible,

> pv=function(k=2){
+ n=1e4
+ X=rnorm(n,10,2)
+ Q=quantile(X,(0:k)/k)
+ Q[1]=0
+ Xc=cut(X,Q,labels=1:k)
+ XcN=tapply(X,Xc,mean)
+ Xn=XcN[as.numeric(Xc)]
+ Y=1+3*Xn+rnorm(n)
+ ad.test(Y)$p.value} > plot(2:100,Vectorize(pv)(2:100),type="l") > abline(h=.05,col="red") So here, it could be possible to have also a Gaussian distribution, for $Y$. As least to accept that assumption, statistically. In the context of a Poisson regression, it is well know that it’s not possible to have at the same time $Y|X$ that is Poisson distributed (that’s a Poisson regression) and also $Y$ that is Poisson distributed. That simply comes from the fact that $\mathbb{E}[Y]=\mathbb{E}[\mathbb{E}[Y\vert X]]$ while $\text{Var}[Y]=\mathbb{E}[\text{Var}[Y\vert X]]+\text{Var}[\mathbb{E}[Y\vert X]]$ and because of the conditional Poisson distribution, then $\text{Var}[Y\vert X]=\mathbb{E}[Y\vert X]$ Thus, $\text{Var}[Y]=\mathbb{E}[Y]+\underbrace{\text{Var}[\mathbb{E}[Y\vert X]]}_{>0}$ So $Y$ cannot be Poisson distribution. But again, it could be possible, if heterogeneity is not too large, to accept the null assumption of a Poisson distribution for $Y$. More generally, it is very difficult to have a distribution family for $Y|X$ that is also the distribution of the non-conditional variable $Y$. In the context of a finite mixture ($X$ takes a finite number of values),Teicher (1963) proved that it was not not possible, neither for the Gaussian distribution nor the Gamma distribution. An to go further, check Monfrini (2002) (thanks Romuald for point out the reference). Hence, as a keep saying, before running a regression model on$Y|X$ with some given family, it is never a good idea to check if the non-conditional distribution $Y$ has the same distribution. Because there is no reason, usually, to remain in the same family. # De la fréquence des ouragans aux Etats-Unis Pour la fin du cours de statistique, on joue enfin avec des données. Après les taux de réussite au brevet des collèges, on a travaillé un peu sur la survenance des ouragans aux Etats-Unis. La base utilisée est la suivante > base=read.table( "http://freakonometrics.free.fr/ouragan.csv", sep=";",header=TRUE,dec=",") > base=base[1:207,1:7] Le nombre annuel d’ouragan est la série suivante > TB <- table(base$Year)
>  years <- as.numeric(names(TB))
>  counts <- as.numeric(TB)
>  years0=(1900:2005)[which(!(1900:2005)%in%years)]
>  db <- data.frame(years=c(years,years0),
counts=c(counts,rep(0,length(years0))))
>  X=db$counts > plot(db,type="h") On peut comparer les fréquences observées avec ce que donnerait une loi de Poisson (paramètre obtenu par maximum de vraisemblance) > lambda=mean(X) > n=length(X) > cbind(X,dpois(0:6,lambda)*n) [,1] [,2] 0 16 15.038430 1 34 29.367500 2 27 28.674870 3 13 18.665717 4 5 9.112744 5 6 3.559128 6 5 1.158396 On peut alors faire un test d’ajustement, > library(vcd) > fit=goodfit(X,"poisson") > plot(fit) > summary(fit) Goodness-of-fit test for poisson distribution X^2 df P(> X^2) Likelihood Ratio 14.17686 5 0.01452415 On voit que le test d’ajustement d’une loi Poisson n’est pas très concluant ici. Notons que ni la loi binomiale > fit=goodfit(X,"binomial") > plot(fit) > summary(fit) Goodness-of-fit test for binomial distribution X^2 df P(> X^2) Likelihood Ratio 46.21196 5 8.222959e-09 ni la loi binomiale négative ne semble marcher. > fit=goodfit(X,"nbinomial") > plot(fit) > summary(fit) Goodness-of-fit test for nbinomial distribution X^2 df P(> X^2) Likelihood Ratio 10.91507 4 0.02753523 On peut faire un test de stabilité temporelle du paramètre d’une loi de Poisson, à l’aide du test de Przyborowski & Wilenski (1940). Un ancien package permettait de mettre en place ce test. Mais on peut heureusement utiliser une autre fonction, avec une version différente du test (on pourra aussi consulter Krishnamoorthy & Thomson (2004) pour aller plus loin). L’idée du premier test est assez simple: notons que si sur la première période (de durée $T_1$), on a un processus de Poisson de paramètre $\lambda_1$ et$\lambda_2$ sur la seconde période (de durée $T_2$), alors $N_1=\sum_{t\in T_1} N_t \sim \mathcal{P}(\lambda_1 T_1)$ et $N_2=\sum_{t\in T_2} N_t \sim \mathcal{P}(\lambda_2 T_2)$ Un test de $H_0:\lambda_1=\lambda_2$ s’écrit $H_0:\lambda_1/\lambda_2=1$ Notons que la loi conditionnelle de $X_1$ sachant $X_1+X_2$ suit une loi binomiale, donc le paramètre de probabilité est fonction de ce ratio. Plus précisément, conditionnellement à$N_1+N_2=k$,$X_1$ suit un loi binomiale $\mathcal{B}(k,p_{1\vert 2})$ avec $p_{1\vert 2}=\frac{\pi}{1+\pi} \text{ avec }\pi=\frac{T_1}{T_2}\frac{\lambda_1}{\lambda_2}$ A partir de là, il est facile de mettre en place une stratégie de test. Une autre stratégie est d’utiliser $T=\frac{1}{\sqrt{N_1+N_2}}\left(\frac{N_1}{T_1}-\frac{N_2}{T_2}\right)$ > seuil=1960 > X1=db[db$years<seuil,"counts"]
>  X2=db[db$years>=seuil,"counts"] > N1=length(db$years<seuil)
>  N2=length(db$years>=seuil) > poisson.test(c(sum(X1),sum(X2)),c(N1,N2)) Comparison of Poisson rates data: c(sum(X1), sum(X2)) time base: c(N1, N2) count1 = 85, expected count1 = 103.5, p-value = 0.01216 alternative hypothesis: true rate ratio is not equal to 1 95 percent confidence interval: 0.5218662 0.9265666 sample estimates: rate ratio 0.6967213 Bref, on peut accepter que le nombre annuel d’ouragan n’est pas un processus de Poisson homogène dans le temps… # Autres Données pour le cours de Statistiques Encore quelques données pour faire un TP de statistique. Ce sont des données de coûts d’ouragans, aux Etats-Unis, > library(gdata) > db=read.xls( + "http://sciencepolicy.colorado.edu/publications/special/public_data_may_2007.xls", + sheet=1) > stupidcomma = function(x){ + x=as.character(x) + for(i in 1:10){x=sub(",","",as.character(x))} + return(as.numeric(x))} > base=db[,1:4] > base$Base.Economic.Damage=
Vectorize(stupidcomma)(db$Base.Economic.Damage) > base$Normalized.PL05=Vectorize(stupidcomma)(db$Normalized.PL05) > base$Normalized.CL05=Vectorize(stupidcomma)(db$Normalized.CL05) Lire la base est un peu pénible, à cause de virgules qui traînent comme séparateur de milliers. Mais peu importe, le petit code permet de nettoyer la base. • Fréquence d’ouragans La première étape pourra être de travailler sur la fréquence. On pourra ajuster des lois de Poisson, ou binomiales, ou binomiales négatives, et faire des tests de stabilité du paramètre, avec le temps (on pourra couper en 4 ou 5 blocs). > TB <- table(base$Year)
> years <- as.numeric(names(TB))
> counts <- as.numeric(TB)
> years0=(1900:2005)[which(!(1900:2005)%in%years)]
> db <- data.frame(years=c(years,years0),
+       counts=c(counts,rep(0,length(years0))))
> db[88:93,]
years counts
88  2003      3
89  2004      6
90  2005      6
91  1902      0
92  1905      0
93  1907      0
> plot(years,counts,type='h')

• Coût des ouragans

Dans un second temps, on pourra travailler sur le coût des ouragans

> plot(base$Normalized.PL05/1e9,type="h") On pourra tenter d’ajuster des lois, log normale, gamma, Pareto, etc. Faire des tests d’ajustement. Et tester la stabilité des paramètres, dans le temps (en coupant en 2 blocs, par exemple, ou davantage). # Actuariat de l’Assurance Non-Vie #5 Pour le quatrième cours d’actuariat de l’assurance non-vie à l’ENSAE, la semaine prochaine, on abordera la notion de sur-dispersion dans les modèles de comptage. Les slides sont en ligne (la version pdf téléchargeable est plus complète que celle sur slideshare) # Actuariat de l’Assurance Non-Vie #3 Pour le troisème cours d’actuariat de l’assurance non-vie à l’ENSAE (oui, j’ai pris un peu d’avance), on abordera la modélisation de la fréquence de sinistres. Les slides sont en ligne, # Computing AIC on a Validation Sample This afternoon, we’ve seen in the training on data science that it was possible to use AIC criteria for model selection. > library(splines) > AIC(glm(dist ~ speed, data=train_cars, family=poisson(link="log"))) [1] 438.6314 > AIC(glm(dist ~ speed, data=train_cars, family=poisson(link="identity"))) [1] 436.3997 > AIC(glm(dist ~ bs(speed), data=train_cars, family=poisson(link="log"))) [1] 425.6434 > AIC(glm(dist ~ bs(speed), data=train_cars, family=poisson(link="identity"))) [1] 428.7195 And I’ve been asked why we don’t use a training sample to fit a model, and then use a validation sample to compare predictive properties of those models, penalizing by the complexity of the model. But it turns out that it is difficult to compute the AIC of those models on a different dataset. I mean, it is possible to write down the likelihood (since we have a Poisson model) but I want a code that could work for any model, any distribution…. Hopefully, Heather suggested a very clever idea, using her package And actually, it works well. # I Fought the (distribution) Law (and the Law did not win) A few days ago, I was asked if we should spend a lot of time to choose the distribution we use, in GLMs, for (actuarial) ratemaking. On that topic, I usually claim that the family is not the most important parameter in the regression model. Consider the following dataset > db <- data.frame(x=c(1,2,3,4,5),y=c(1,2,4,2,6)) > plot(db,xlim=c(0,6),ylim=c(-1,8),pch=19) To visualize a regression model, use the following code > nd=data.frame(x=seq(0,6,by=.1)) > add_predict = function(reg){ + prd1=predict(reg,newdata=nd,se.fit = TRUE,type="response") + y1=prd1$fit
+ y1_upp=prd1$fit+prd1$residual.scale*1.96*
prd1$se.fit + y1_low=prd1$fit-prd1$residual.scale*1.96* prd1$se.fit
+ polygon(c(nd$x,rev(nd$x)),c(y1_upp,
rev(y1_low)),col="light green",angle=90,
density=40,border=NA)
+ lines(nd$x,y1,col="red",lwd=2) + } For instance, with a Poisson regression (with a log link function) we get > plot(db) > reg1=glm(y~x,family=poisson(link="log"), + data=db) > add_predict(reg1)  while, with a Gaussian regresion (but still with a log link function), we get > plot(db) > reg2=glm(y~x,family=gaussian(link="log"), + data=db) > add_predict(reg2) If we just care about the expected value of our prediction, the output is more or less the same > plot(db) > lines(nd$x,predict(reg1,newdata=nd,
+ type="response"),col="red",lwd=1.5)
> lines(nd$x,predict(reg2,newdata=nd, + type="response"),col="blue",lwd=1.5) So, indeed, forget about the (distribution) law when running a GLM. Not convinced? Consider – on the same dataset – a Poisson regression (with an identity link function this time) > plot(db) > reg1=glm(y~x,family=poisson(link="identity"), + data=db) > add_predict(reg1)  while, with a Gaussian regresion (but still with an identity link function), we get > plot(db) > reg2=glm(y~x,family=gaussian(link="identity"), + data=db) > add_predict(reg2) Again, if we just plot the expected value of our prediction, the output is more or less the same > plot(db) > lines(nd$x,predict(reg1,newdata=nd,
+ type="response"),col="red",lwd=1.5)
> lines(nd$x,predict(reg2,newdata=nd, + type="response"),col="blue",lwd=1.5) So clearly, the simplistic message you should not care too much about the (distribution) law seems to be valid… # Reinterpreting Lee-Carter Mortality Model Last week, while I was giving my crash course on R for insurance, we’ve been discussing possible extensions of Lee & Carter (1992) model. If we look at the seminal paper, the model is defined as follows # The odds of a cluster of airplane accidents Recently, there have been a lot of airplane accidents. • July, 17th 2014, Hrabove, Ukraine, Malaysia Airlines, Boeing 777, fatalities 298 (/298) • July, 23rd 2014, Magong, Taiwan, TransAsia Airways, ATR 72-500, fatalities 47 (/58) • July, 24th 2014, Aguelhok, Mali, Air Algerie, Mc Donnell Douglas MD-83, fatalities 116 (/116) It is simple to find a lot of datasets about airplane crashes. For instance on http://ntsb.gov/aviationquery. The dataset is nice, with a lot of information, > planes=read.table( + "cbad3ca6-6b8f-4c98-9ee0-601faAviationData.txt", + sep="|",header=TRUE) for instance the exact location of the crashes, > library(maps) > map("world", interior = FALSE) > points(planes$Longitude,planes$Latitude, + pch=19,cex=planes$Total.Fatal.Injuries/50,
+ col="red")

# La Loi des Petits Nombres

Comme nous l’avions vu dans Charpentier (2010), la loi des grands nombres est souvent évoquée pour justifier la mutualisation des risques indépendants : plus la mutualité sera grande, plus petite sera la variabilité. A condition que les risques ne soient pas trop grands. Un cas problématique sont les risques catastrophiques : ils sont rares, et parfois tellement (potentiellement) coûteux que l’hypothèse d’existence de la variance doit être remise en cause. La modélisation de ces événements rares repose sur la loi des petits nombres, pour reprendre le nom de l’ouvrage de Ladislaus Bortkiewicz. Comme nous allons le voir, la loi de Poisson est à la loi des petits nombres ce que la loi normale est à la loi des grands nombres. Nous verrons donc en détails l’importance de la loi de Poisson pour modéliser les événements rares. Et nous verrons pourquoi la probabilité qu’un événement ne survienne pas est toujours de 37%. Ou presque.

• La loi de Siméon-Denis Poisson

Siméon-Denis Poisson a travaillé sur les calculs des probabilités pendant presque 20 ans, de 1820 à 1840. Comme pour beaucoup de ses contemporains, ses premiers travaux portèrent sur des problèmes de jeux, avec une communication à l’Académie des Sciences sur l’avantage du banquier au jeu de trente et quarante. Mais son premier travail conséquent a porté sur un problème longuement étudié par Laplace (qui a été un de ses professeurs) sur la proportion des naissances des filles et des garçons (problème classique de statistique) publié en 1830. Il y présente en particulier une démonstration de la loi des grands nombres pour une loi de Bernoulli, qu’il modifiera par la suite. Dans ce mémoire, il y présente une loi, qui sera appelée plus tard la loi de Poisson : « la probabilité qu’une événement dont la chance à chaque épreuve est la fraction très petite x/mu n’arrive pas plus de $n$ fois dans un très grand nombre d’épreuve mu d’épreuve (pour reprendre la terminologie utilisée en 1837, ce qui correspond à la fonction de répartition avec une terminologie plus contemporaine) est »

$P=\left(1+x+\frac{x^2}{1\cdot 2}+\frac{x^3}{1\cdot 2\cdot 3}+\cdots+\frac{x^n}{1\cdot 2\cdots n}\right)e^{-x }$

On verra réapparaître cette loi dans son fameux traité, paru en 1837, Recherches sur la probabilité des jugements. En particulier, dans le chapitre 8, il obtient sa loi comme limite de la loi binomiale $B(T,\lambda /T)$ lorsque $T$ devient grand, mais passe assez rapidement à d’autres considérations. Il n’étudie pas cette loi limite, et ne propose pas vraiment de l’utiliser. Il faudra attendre les travaux de Ladislaus Bortkiewicz, presque un siècle plus tard, pour voir des applications, dans Das Gesetz der kleinen Zahlen (la loi des petits nombres). Et comme toujours en mathématiques, lui attribuer tout le crédit est un peu excessif, puisqu’un siècle auparavant, en 1718, de Moivre avait obtenu la loi même loi, toujours comme limite de la loi Binomiale. Si le nom de Poisson est resté à la postérité, c’est essentiellement parce que Boltzmann le cite, en 1868, ainsi que Seidel en 1876, et surtout Tchebychev. Cela dit, Poisson était loin d’être un inconnu, en tant que scientifique. Ses travaux l’ont amené à travailler sur des problèmes d’électrostatique (la fameuse équation de Poisson), à l’équation de la chaleur avec Fourier, il s’est opposé à Fresnel sur des problèmes d’optique, et il a présidé à deux reprises l’Académie des Sciences.

La distribution d’un nombre d’événements obtenu comme somme de variables de Bernoulli, dans une grande population, peut s’approcher par la loi de Poisson. Ou lorsque la probabilité de survenance est faible, par rapport à la taille de l’échantillon. Formellement, si $N$ suit une loi binomiale $B(n,p)$, avec $p\sim\lambda/n$, alors

$P[N=k]=\binom{n}{k}p^k(1-p)^{n-k}\sim e^{-\lambda}\frac{\lambda^k}{k!}$

On parlera de loi des petits nombres car on compte ici des évènements rares (la probabilité de survenance d’un évènement étant inversement proportionnelle à $n$).

Ce résultat se traduit de la manière suivante: si on se donne un échiquier, $10\times 10$ si on lance 100 pièces et si on compte le nombre de pièces par case, la distribution suivra une loi de Poisson de moyenne 1 (le ratio entre le nombre de pièces et de cases). Une illustration est évoquée sur la Figure suivante

 Nombre de pièces par case Fréquence Loi de Poisson 0 36 36,78 1 39 36,78 2 16 18,39 3 7 6,13 4 2 1,53 5 et plus 0 0,37
• L’utilisation de la loi de Poisson, et les mélanges Poissonniens

A la fin du XIXème siècle, à l’université de Göttingen, Wilhelm Lexis (connu en démographie pour les diagrammes qui portent aujourd’hui son nom) avait envie d’étudier les applications de cette loi. Alors qu’il étudiait l’utilisation de la loi de Poisson dans un contexte démographique, un de ces étudiants, Ladislaus Bortkiewicz, proposa de l’utiliser pour modéliser des nombres d’accidents. Dans l’exemple fameux de Bortkiewicz, datant de 1898, il avait étudié le nombre de soldats de cavaliers morts par ruade de cheval, entre 1875 et 1894, dans 10 corps (soit 200 corps annuels). Il avait obtenu la distribution donnée dans la table suivante. La colonne de droite donne le nombre de décès donné par une loi de Poisson de même moyenne. On note que l’ajustement est très bon, et on comprend vite la fascination de Bortkiewicz pour cette loi.

 Nombre de décès, par corps Fréquence Loi de Poisson 0 109 108,67 1 65 66,21 2 22 20,22 3 3 4,11 4 1 0,63 5 et plus 0 0,08

Maintenant, pour être tout à fait honnête, la loi que nous avons énoncée comme une loi des petits nombres est très différente de celle proposée par Bortkiewicz (celle qu’il a énoncé dans son ouvrage lui a valu un cinglant commentaire de Corrado Gini, qui écrivait en 1907, de manière provocatrice que « the law of small numbers does not exist »). Mais l’utilisation de la loi de Poisson pour modéliser les accidents, et plus généralement les petits nombres venait de débuter.

Dans les années 1930, Filip Lunderg avait noté l’intérêt (théorique) du processus de Poisson pour modéliser l’arrivée des sinistres, ainsi que plusieurs actuaires de l’école scandinave (Esscher en 1932, Sgerdahl en 1939, Lüders en 1934, etc). Depuis, tous les actuaires ont utilisé cette loi pour modéliser toute sorte d’évènements “rares“, comme la survenance annuelle d’ouragans, aux États-Unis (Figure ci-dessous).

 Nombre d’ouragans Fréquence Loi de Poisson 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 et plus 0 0,08

• De la loi de Poisson à la période de retour

Un concept fondamental en gestion des risques extrêmes a été introduit par Emil Gumbel en 1958, liant le temps qui s’écoule (en années) entre deux événements consécutifs et la probabilité (annuelle) de survenance. Pour des événements se produisant avec une probabilité annuelle $1/T$, indépendamment les uns des autres, le temps moyen d’attente entre deux événements est $T$, appelé période de retour, et la probabilité qu’aucun événement ne survienne pendant $n$ années (consécutives) est alors

$\mathbb{P}\left(N>n\right)=\left(1-\frac{1}{T}\right)^n$

On peut résumer ceci dans le tableau ci-dessous, où on étudie probabilité qu’un événement ne survienne pas pendant n années (en ligne) en fonction de la période de retour $T$. Pour un événement centenaire, il y a 36,60% chance pour qu’il ne survienne pas en cent ans.

 Période de retour (en années) Nombre d’années sans catastrophe 10 20 50 100 200 10 34,86% 59,87% 90,43% 90,43% 95,11% 20 12,15% 35,84% 66,76% 81,79% 90,46% 50 0,51% 7,69% 36,41% 60,50% 77,83% 100 0,00% 0,59% 13,26% 36,60% 60,57% 200 0,00% 0,00% 1,75% 13,39% 36,69%
• La probabilité qu’un événement ne survienne pas est… 37%

Comme on le voit dans le tableau ci-dessus, sur la diagonale, la probabilité qu’un événement ne survienne pas pendant $T$ années quand sa période de retour est $T$ est de l’ordre de 37%. D’ailleurs, si on revient un instant sur notre échiquier $10\times 10$ évoqué auparavant, la probabilité de n’avoir aucune pièce sur une case était de… 37%

Reprenons ici un exemple qui avait fait polémique il y a quelques années, sur le risque nucléaire. Dans un article (Laponche et Dessus, 2011), on apprenait que la probabilité d’avoir un incident majeur sur réacteur nucléaire, pour une année, était de l’ordre de 0,0003 (3 chances sur 10,000). Or comme il y a 143 réacteurs nucléaires, sur 30, la probabilité d’avoir un incident majeur est (selon le calcul des auteurs)

$\underbrace{143\times30}_n\times\underbrace{0,0003}_p\sim 129\%$

D’où la conclusion savoureuse, « la probabilité d’occurrence d’un accident majeur sur ces parcs serait donc […] de plus de 100% pour l’Union européenne ». Comme on l’a noté, la probabilité d’avoir un incident majeur sur une période $n$ lorsque la probabilité annuelle est $p$ s’écrit

$P[N\leq n]=1-(1-p)^n \sim np$

en utilisant un développement limité non justifié ici ! Il convient ici d’utiliser le modèle de Poisson. La probabilité d’avoir au moins un incident majeur est

$P[N\leq n]=1-(1-p)^n \sim 1-e^{-np}$

soit ici 72,47%. Si la probabilité était traduite en terme de durée

$p=\frac{1}{T}$

avec

$T=\frac{1}{0,0003\times143}\sim 23,31$

on voit qu’ici, le temps moyen d’attente entre deux incidents majeurs en Europe est de 23 ans. Ce qui fait que sur 23 ans, la probabilité de n’avoir aucun incident majeur est de l’ordre de 37%.

• Conclusion

La loi de Poisson est présente partout en assurance, car c’est la loi centrale pour modéliser le comptage d’événements rares. Le nombre de décès dans un portefeuille d’assurance-vie suit une loi de Poisson, tout comme le nombre d’accidents par unité de temps en assurance automobile. On retrouve même l’estimateur Chain Ladder du montant de provisions pour sinistres a payer quand on utilise une régression de Poisson sur les incréments de paiements.

# 37% chance

I don’t know if you ever realized, before, but it is quite common to have 37% chance that something happened (or actually “not happened” if we want to be more rigorous). For instance, consider a $n\times n$ grid, and draw $n^2$ points randomly (and uniformely). Then, around 37% cells are empty. Or if you consider a cell, on that grid, there is 37% chance, that the cell is empty. You can look, on the animation below,

Actually, it is quite simple to prove this result. This come from the fact that

$\exp(-1)=0.36787944117\sim 0.37$

And $\exp(-1)$ is a common probability. To be more specific, we can write it

$\exp(-1)\frac{1^0}{0!}$

which is $\mathbb{P}(N=0)$ when $N$ is a Poisson distribution with parameter $\lambda=1$ (or with mean $1$). In our case, $N$ is the number of points in a given cell, and since there are $n^2$  points for $n^2$  cells (it is a $n\times n$ grid), then $\mathbb{E}(N)=1$. And the Poisson distribution arise here as a limit of a binomial distribution. More precisely, here $\mathcal{B}(n^2,1/n^2)$. We have here the law of small number. Nice isn’t it?

# GLM, non-linearity and heteroscedasticity

Last week in the non-life insurance course, we’ve seen the theory of the Generalized Linear Models, emphasizing the two important components

• the link function (which is actually the key component in predictive modeling)
• the distribution, or the variance function

Just to illustrate, consider my favorite dataset

­lin.mod = lm(dist~speed,data=cars)

A linear model means here

$Y_i=\beta_0+\beta_1 X_i +\varepsilon_i$

where the residuals are assumed to be centered, independent, and with identical variance. If we visualize that linear regression, we usually see something like that

The idea here (in GLMs) is to assume

$Y\vert X=x\sim\mathcal{N}( \beta_0+\beta_1 x,\sigma^2)$

which will produce the same model as the one describe previously, based on some error term. That model can be visualized below,

attach(cars)
n=2
X= cars$speed Y=cars$dist
df=data.frame(X,Y)
vX=seq(min(X)-2,max(X)+2,length=n)
vY=seq(min(Y)-15,max(Y)+15,length=n)
mat=persp(vX,vY,matrix(0,n,n),zlim=c(0,.1),theta=-30,ticktype ="detailed", box = FALSE)
x=seq(min(X),max(X),length=501)
C=trans3d(x,predict(reggig,newdata=data.frame(X=x),type="response"),rep(0,length(x)),mat)
lines(C,lwd=2)
sdgig=sqrt(summary(reggig)$dispersion) x=seq(min(X),max(X),length=501) y1=qnorm(.95,predict(reggig,newdata=data.frame(X=x),type="response"), sdgig) C=trans3d(x,y1,rep(0,length(x)),mat) lines(C,lty=2) y2=qnorm(.05,predict(reggig,newdata=data.frame(X=x),type="response"), sdgig) C=trans3d(x,y2,rep(0,length(x)),mat) lines(C,lty=2) C=trans3d(c(x,rev(x)),c(y1,rev(y2)),rep(0,2*length(x)),mat) polygon(C,border=NA,col="yellow") C=trans3d(X,Y,rep(0,length(X)),mat) points(C,pch=19,col="red") n=8 vX=seq(min(X),max(X),length=n) mgig=predict(reggig,newdata=data.frame(X=vX)) sdgig=sqrt(summary(reggig)$dispersion)
for(j in n:1){
stp=251
x=rep(vX[j],stp)
y=seq(min(min(Y)-15,qnorm(.05,predict(reggig,newdata=data.frame(X=vX[j]),type="response"), sdgig)),max(Y)+15,length=stp)
z0=rep(0,stp)
z=dnorm(y, mgig[j], sdgig)
C=trans3d(c(x,x),c(y,rev(y)),c(z,z0),mat)
polygon(C,border=NA,col="light blue",density=40)
C=trans3d(x,y,z0,mat)
lines(C,lty=2)
C=trans3d(x,y,z,mat)
lines(C,col="blue")}

We do have two parts here: the linear increase of the average, $\mathbb{E}(Y\vert X=x)=\beta_0+\beta_1 x$ and the constant variance of the normal distribution $\text{Var}(Y\vert X=x)=\sigma^2$.

On the other hand, if we assume a Poisson regression,

poisson.reg = glm(dist~speed,data=cars,family=poisson(link="log"))

we have something like

This time, two things have changed simultaneously: our model is no longer linear, it is an exponential one $\mathbb{E}(Y\vert X=x)=e^{\beta_0+\beta_1 x}$, and the variance is also increasing with the explanatory variable $\text{Var}(Y\vert X=x)=e^{\beta_0+\beta_1 x}$, since with a Poisson regression,

$Y\vert X=x\sim\mathcal{P}(e^{ \beta_0+\beta_1 x})$

If we adapt the previous code, we get

The problem is that we changed two things when we introduced the Poisson regression from the linear model. So let us look at what happens when we change the two components independently. First, we can change the link function, with a Gaussian model but this time a multiplicative model (with a logarithm link function)

gaussian.reg = glm(dist~speed,data=cars,family=gaussian(link="log"))

which is still, here, an homoscedasctic model, but this time non-linear. Or we can change the link function in the Poisson regression, to get a linear model, but heteroscedastic

poisson.lin = glm(dist~speed,data=cars,family=poisson(link="identity"))

So this is basically what GLMs are about….

# Surdispersion et comptage

Cette semaine, au cours d’assurance non-vie, on abordera la surdispersion, qui clôturera la partie du cours sur la modélisation de la fréquence de sinistres. Les transparents sont en ligne. Mais avant de parler de surdispersion, on finira la présentation des GLM. Je mets un lien vers le chapitre 15 du livre de John Fox Applied regression analysis and generalized linear models ainsi que le livre de James K. Lindsey Applying Generalized Linear Models. Je voudrais aussi renvoyer vers les notes de cours de Germán Rodríguez, avec des notes sur la régression de Poisson (avec un petit complément sur la notion de overdispersion).

Les instructions pour le second devoir seront envoyées par courriel.