Tag Archives: Poisson

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 ), on a un processus de Poisson de paramètre  et sur la seconde période (de durée ), alors

et

Un test de  s’écrit  Notons que la loi
conditionnelle de  sachant  suit une loi binomiale, donc le paramètre de probabilité est fonction de ce ratio. Plus précisément, conditionnellement à,  suit un loi binomiale  avec

A partir de là, il est facile de mettre en place une stratégie de test. Une autre stratégie est d’utiliser

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

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.

Continue reading Computing AIC on a Validation Sample

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…

Continue reading I Fought the (distribution) Law (and the Law did not win)

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

Continue reading The odds of a cluster of airplane accidents

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 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 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  grid, and draw  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

And  is a common probability. To be more specific, we can write it

which is  when  is a Poisson distribution with parameter  (or with mean ). In our case,  is the number of points in a given cell, and since there are   points for   cells (it is a  grid), then . And the Poisson distribution arise here as a limit of a binomial distribution. More precisely, here . 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_1X_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\vertX=x\sim\mathcal{N}(\beta_0+\beta_1x,\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)
reggig=glm(Y~X,data=df,family=gaussian(link="identity"))
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_1x 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_1x}, and the variance is also increasing with the explanatory variable \text{Var}(Y\vert X=x)=e^{\beta_0+\beta_1x}, since with a Poisson regression,
Y\vert X=x\sim\mathcal{P}(e^{\beta_0+\beta_1x})

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.

Régression de Poisson

Mercredi, on finira les arbres de classification, et on commencera la modélisation de la fréquence de sinistre. Les transparents sont en ligne.

Comme annoncé lors du premier cours, je suggère de commencer la lecture du Practicionner’s Guide to Generalized Linear Models. Le document correspond au minimum attendu dans ce cours.

Loi de Poisson

Suite du cours ACT2121, de préparation pour l’examen P de la SOA (probability). Un nouveaux jeu d’exercices, sur le thème 6 (tel que classifié dans le livre de Jacques Labelle, qui servira de référence pour ce cours)

Pour rappels, l’intra du 27 septembre portera sur les thèmes 1-6, c’est à dire sur ce qui a été abordé dans les feuilles mises en ligne.

Exposure as a possible explanatory variable

Iin insurance pricing, the exposure is usually used as an offset variable to model claims frequency. As explained many times on this blog (e.g. here), and in my notes, if we have to identical drivers, but one with an exposure of 6 months, and the other one of one year, it should be natural to assume that, on average, the second driver will have two times more accidents. This is the motivation to use a standard (homogeneous) Poisson process to model claim frequency. One can also see here legal issue, since, in case of a (partial) reinbursement of a premium, it would be done prorata temporis. The risk is proportional to the exposure. Thus, if https://latex.codecogs.com/gif.latex?Y_i denote the number of claims of insured https://latex.codecogs.com/gif.latex?i, with characteristics https://latex.codecogs.com/gif.latex?\boldsymbol{X}_{i}=(X_{i,1},\cdots,X_{i,k}) and exposure https://latex.codecogs.com/gif.latex?E_i, with a Poisson regression, we would write

https://latex.codecogs.com/gif.latex?Y_i\sim\mathcal{P}(E_i\cdot%20\exp(\boldsymbol{X}_i%27\boldsymbol{\beta}))

or equivalently

https://latex.codecogs.com/gif.latex?Y_i\sim\mathcal{P}(\exp(\log(E_i)+\boldsymbol{X}_i%27\boldsymbol{\beta}))

From this expression, the logarithm of the exposure is an explanatory variable, but there should be no coefficient (the coefficient here is taken to be one). Can’t we use the exposure as an explanatory variable ? Will we get a unit parameter ?

Of course, in the context of ratemaking, it is probably not a relevant question, since actuaries are required to predict annual claim frequency (since insurance contract are supposed to provide a one year coverage). But it might be interesting to get a better understanding of why people might be leaving our portfolio (i.e. are cancelling their insurance policy before term, or not renew someday).

To be more specific and get a better understanding, consider the following model: consider a Poisson process to model claims arrival, and people dedicated to their insurance company (they never leave). Let us generate scenarios over twenty years

> n=983
> D1=as.Date("01/01/1993",'%d/%m/%Y')
> D2=as.Date("31/12/2013",'%d/%m/%Y')
> L=D1+0:(D2-D1)
> set.seed(1)
> arrival=sample(L,size=n,replace=TRUE)
> exposure=N=rep(NA,n)
> departure=rep(D2,n)
> set.seed(2)
> for(i in 1:n){
+   expo=D2-arrival[i]
+   w=0
+   while(max(w)<expo) w=c(w,max(w)+1+trunc(rexp(1,1/1000)))
+   exposure[i]=departure[i]-arrival[i]
+   N[i]=max(0,length(w)-2)}
> df=data.frame(N=N,E=exposure/365)

Here the expected time between claims is considered to be 1000 days. The (annual) intensity of the Poisson process is here

> 365/1000
[1] 0.365

so if we run a Poisson regression on the logarithm of the exposure (please feel free to had other covariates if you want, the example here is just to see what could happen when exposure is considered as a standard covariate), we should get a parameter close to

> log(365/1000)
[1] -1.007858

Here, the regression on a constant, with the offset variable is

> reg=glm(N~1+offset(log(E)),data=df,family=poisson)
> summary(reg)

Call:
glm(formula = N ~ 1 + offset(log(E)), family = poisson, data = df)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.4145  -0.4673   0.2367   0.8770   3.6828  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.04233    0.02532  -41.17   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 1116.9  on 982  degrees of freedom
Residual deviance: 1116.9  on 982  degrees of freedom
AIC: 3282.9

Number of Fisher Scoring iterations: 5

which is consistent with what we just said. If we run the regression with the logarithm of the exposure as a possible explanatory variable, we would expect to have a coefficient close to 1. And indeed…

> reg=glm(N~log(E),data=df,family=poisson)
> summary(reg)

Call:
glm(formula = N ~ log(E), family = poisson, data = df)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.0810  -0.8373  -0.1493   0.5676   3.9001  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.03350    0.08546  -12.09   <2e-16 ***
log(E)       1.00920    0.03292   30.66   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 2553.6  on 982  degrees of freedom
Residual deviance: 1064.2  on 981  degrees of freedom
AIC: 3762.7

Number of Fisher Scoring iterations: 5

If we keep the offset, and add the variable, we can see that it become useless (which is a test of a unit parameter, somehow)

> reg=glm(N~log(E)+offset(log(E)),data=df,family=poisson)
> summary(reg)

Call:
glm(formula = N ~ log(E) + offset(log(E)), family = poisson, 
    data = df)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.0810  -0.8373  -0.1493   0.5676   3.9001  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.033503   0.085460 -12.093   <2e-16 ***
log(E)       0.009201   0.032920   0.279     0.78    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 1064.3  on 982  degrees of freedom
Residual deviance: 1064.2  on 981  degrees of freedom
AIC: 3762.7

Number of Fisher Scoring iterations: 5

Here, we do have pure Poisson processes, so exposure is crucial, since the parameter of the Poisson distribution is proportional to the exposure. But we cannot learn anything else from the exposure.

Consider some real data.

> head(baseFREQ)
  nocontrat exposition zone puissance agevehicule
1        27       0.87    C         7           0
2       115       0.72    D         5           0
3       121       0.05    C         6           0
4       142       0.90    C        10          10
5       155       0.12    C         7           0
6       186       0.83    C         5           0
  ageconducteur bonus marque carburant densite region nbre
1            56    50     12         D      93     13    0
2            45    50     12         E      54     13    0
3            37    55     12         D      11     13    0
4            42    50     12         D      93     13    0
5            59    50     12         E      73     13    0
6            75    50     12         E      42     13    0

What do we get if we consider a Poisson regression on the logarithm of the exposure ?

> reg=glm(nbre~log(exposition),data=baseFREQ,family=poisson)
> summary(reg)

Call:
glm(formula = nbre ~ log(exposition), family = poisson, data = baseFREQ)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.3988  -0.3388  -0.2786  -0.1981  12.9036  

Coefficients:
                Estimate Std. Error z value Pr(>|z|)    
(Intercept)     -2.83045    0.02822 -100.31   <2e-16 ***
log(exposition)  0.53950    0.02905   18.57   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 12931  on 49999  degrees of freedom
Residual deviance: 12475  on 49998  degrees of freedom
AIC: 16150

Number of Fisher Scoring iterations: 6

If we add the exposure to the offset, what’s happening ? (let us use a nonparametric transformation, so visualize what’s going on)

> library(gam)
> reg=gam(nbre~offset(log(exposition))+s(exposition),data=baseFREQ,family=poisson)
> plot(reg,se=TRUE)

There is a clear and significant effect. The more insured stay, the less likely they get a claim. Actually, it can be observed without running a regression.

> i1=which(baseFREQ$nbre>0)
> i0=which(baseFREQ$nbre==0)
> h1=hist(baseFREQ$exposition[i1],probability=TRUE)
> h0=hist(baseFREQ$exposition[i0],probability=TRUE)
> plot(h1$mids,h1$density,type='s',lwd=2,col="red")
> lines(h0$mids,h0$density,type='s',col='blue',lwd=2)

In blue, we have the density of the exposure for those who did not have claims, and in red, the density of those who did have one claim (or more)

So here, we cannot assume a unit value for the parameter. What does that mean ? Can we reproduce such a behavior ?

In order to get a better understandung, consider two possible behaviors for the insured. The first one will be the following : if the company does not offer substantial discounts after no several years with no claims, the insured might leave the company. For instance, if the insured has no claim during 5 years, then after 5 years, he will leave the company (to get a better price somewhere else, say). The code will be

> for(i in 1:n){
+   expo=D2-arrival[i]
+   w=c(0,0)
+   while((max(w)<expo) & (max(diff(w))<1500)) w=c(w,max(w)+trunc(rexp(1,1/1000)))
+   if(max(diff(w))>1500) departure[i]=arrival[i]+max(w[-length(w)])+1500
+   exposure[i]=departure[i]-arrival[i]
+   N[i]=max(0,length(w)-3)}
> df=data.frame(N=N,E=exposure/365)

Here, I consider 1500 days, instead of 5 years,, but it is the same idea. So, what do we have here ?

> reg=glm(N~log(E),data=df,family=poisson)
> summary(reg)

Call:
glm(formula = N ~ log(E), family = poisson, data = df)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.5684  -0.9668  -0.2321   0.4244   3.6265  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -2.50844    0.10286  -24.39   <2e-16 ***
log(E)       1.65738    0.04494   36.88   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 2567.31  on 982  degrees of freedom
Residual deviance:  885.71  on 981  degrees of freedom
AIC: 2897.9

Here, the coefficient is (significantly) larger than 1. More precisely,

> reg=glm(N~log(E)+offset(log(E)),data=df,family=poisson)
> summary(reg)

Call:
glm(formula = N ~ log(E) + offset(log(E)), family = poisson, 
    data = df)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.5684  -0.9668  -0.2321   0.4244   3.6265  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -2.50844    0.10286  -24.39   <2e-16 ***
log(E)       0.65738    0.04494   14.63   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 1114.24  on 982  degrees of freedom
Residual deviance:  885.71  on 981  degrees of freedom
AIC: 2897.9

There is clearly a bias here : people staying long are more like likely to have an accident. Which is consistent with our story, since clients with low risks left.

The second behavior will be the following : sometimes, the insured are not satisfied with the way claims are handled, and they might leave after the first claim. Consider the case where, after one claim, it is likely (e.g. with probability 50%) that the insured leaves the company. Instead of assuming that the insured did not like claims management, consider the case were the car is so damaged that he cannot drive it anymore. So it will be useless to pay an insurance premium. The code here will be

> for(i in 1:n){
+   expo=D2-arrival[i]
+   w=0
+   stay=TRUE
+   while((max(w)<expo) & (stay==TRUE)) { w=c(w,max(w)+trunc(rexp(1,1/1000)))
+   stay=sample(c(TRUE,FALSE),prob=c(.5,.5),size=1)}
+   N[i]=length(w)-2
+   if(stay==FALSE) {departure[i]=arrival[i]+max(w)
+   N[i]=length(w)-1}
+   exposure[i]=departure[i]-arrival[i]}
> df=data.frame(N=N,E=exposure/365)

Here, after each claim, the insured toss a coin to see if he cancels the contract, or not.

> reg=glm(N~log(E),data=df,family=poisson)
> summary(reg)

Call:
glm(formula = N ~ log(E), family = poisson, data = df)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.28402  -0.47763  -0.08215   0.33819   2.37628  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.09920    0.04251   2.334   0.0196 *  
log(E)       0.30640    0.02511  12.203   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 666.92  on 982  degrees of freedom
Residual deviance: 498.29  on 981  degrees of freedom
AIC: 2666.3

This time, the parameter is (again significantly) smaller than one.

> reg=glm(N~log(E)+offset(log(E)),data=df,family=poisson)
> summary(reg)

Call:
glm(formula = N ~ log(E) + offset(log(E)), family = poisson, 
    data = df)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.28402  -0.47763  -0.08215   0.33819   2.37628  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.09920    0.04251   2.334   0.0196 *  
log(E)      -0.69360    0.02511 -27.625   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 1116.87  on 982  degrees of freedom
Residual deviance:  498.29  on 981  degrees of freedom
AIC: 2666.3

The story is now rather different, since those who stay long should not have encountered a lot of opportunities to leave. So clearly, they did not have much claims. If someone has a long exposure, the negative sign in the output above means that he should not have much claims, on average.

As we can see, those models produce rather difference outputs. Note that it is possible much more interpretations. For instance, depending on the way data were extracted,

  • all policies observed, over those twenty years,
  • all policies in force at some specific date, until now
  • all policies in force at some specific date, until one year after
  • all policies in force now

So far, we have been using the first method, but the other ones will yield different interpretations, e.g. because of survivor bias. But that’s another story… And one can read Boucher and Denuit (2008) to go further.