Uncertainty in claims reserving (WIM)

I will be talking on Friday, at the 1st Québec-Ontario Workshop on Insurance Mathematics (so called WIN, already mentioned here). The program is now online here. My talk will be on Solvency II requirements in claims reserving (slides can be found here). Even if Canada will not adapt the Solvency II European capital test, it looks like Solvency II will matter to Canadian insurance companies (as mentioned here).

You find it cold in Montréal ? trust me, it is even worse

As people say in Montréal, “aujourd’hui, il fait frette”. And I have been surprised recently when some people told my that we would reach -35°C Sunday evening… I checked around, and I found -25°C on all weather forecast websites. But nowhere -35°C. I asked some friends, and they told me that those people were not really looking at the air temperature (as we observe on the thermometer), but they were looking at the wind chill, also called “felt air temperature on exposed skin due to the wind” (température ressentie).
And indeed, such a quantity does exist, and can be found on theclimate.weatheroffice.gc.ca website. There is also a physical background for that quantity. Hence, the windchill is http://freakonometrics.blog.free.fr/public/maths/windchill2.png defined as

http://freakonometrics.blog.free.fr/public/maths/windchill1.png

where http://freakonometrics.blog.free.fr/public/maths/windchill3.png is the air temperature (in °C), and http://freakonometrics.blog.free.fr/public/maths/windchill4.png the wind speed (in km/h). Please don’t ask me how to interpret this power 0.16 (I already find difficult to explain a square root in an econometric equation). If we look at the past previous days we observe the following observations,
where points on top are temperature, while below we have felt temperature.So, basically, winters are even colder than what you might think..

And the story is not over, yet. The same thing holds for summer: if you take into account humidity, summer are even hotter than what you think… There is thehumidexhttp://freakonometrics.blog.free.fr/public/maths/humidex2.png, defined here as

http://freakonometrics.blog.free.fr/public/maths/humidex.png

where http://freakonometrics.blog.free.fr/public/maths/humidex3.png denotes a dewpoint (see here for more details).

That index appeared in the 70’s, with a work of Masterson and Richardson entitled a method of quantifying human discomfort due to excessive heat and humidity (published in 1979).By that time, in Canada, on average, 22 people died, per year, because of those excessive heat and humidity. For those interested by the origin of that index, you can have a look here.
Recently, @Annmaria (here) told me that one might expect variance to increase, i.e. maximas should be increasing faster than minimas. I just wonder if this intuition can be related to the fact that more and more people (including some medias) now talk more about felt temperatures than measured temperatures. And if we compare past temperatures to felt temperature we have today, it looks like the difference between extremes is increasing….

Estimer une variance… ou un écart-type ?

Lors du cours de demain, nous allons parler (longuement) de l’estimation de la variance (comme cela avait commencé en TD vendredi dernier). Mais en pratique, c’est souvent l’écart-type qui nous intéresse. On va passer par mal de temps à chercher un estimateur sans biais de la variance, et montrer que

devrait convenir, mais cela ne nous donne pas vraiment d’information sur le biais de l’estimateur naturel de l’écart-type, à savoir

car http://freakonometrics.blog.free.fr/public/maths/variance-est03.png.
Bref, il y a sûrement un biais, et ça pourrait être intéressant de le quantifier. Classiquement, on suppose disposer d’un échantillon http://freakonometrics.free.fr/blog/sigma01.png i.i.d. On va faire une hypothèse forte, à savoir que leur loi est une loi normale http://freakonometrics.free.fr/blog/sigma03.png, où http://freakonometrics.free.fr/blog/sigma02.png est inconnu.
Naturellement, on voudrait considérer

http://freakonometrics.free.fr/blog/sigma04.png

On va poser, afin de simplifier un peu,

http://freakonometrics.free.fr/blog/sigma07.png
et on va chercher à calculer

http://freakonometrics.free.fr/blog/sigma06.png

où http://freakonometrics.free.fr/blog/sigma08.png est la densité de cette variable. Or on sait que http://freakonometrics.free.fr/blog/sigma20.png suit une loi du chi-deux à http://freakonometrics.free.fr/blog/sigma11.png degrés de liberté, dont la densité est la loi gamma suivante

http://freakonometrics.free.fr/blog/sigma13.png

En notant que l’on peut réécrire

http://freakonometrics.free.fr/blog/sigma22.png

on a que cette espérance s’écrit

http://freakonometrics.free.fr/blog/sigma24.png

soit en faisant un peu de calculs,

http://freakonometrics.free.fr/blog/sigma26.png

ou encore

http://freakonometrics.free.fr/blog/sigma27.png

Ce qui va se simplifier en utilisant la fonction gamma

http://freakonometrics.free.fr/blog/sigma28.png

On a ainsi

http://freakonometrics.blog.free.fr/public/perso/stdvariance.png

http://freakonometrics.free.fr/blog/sigma35.png

Quelques ordres de grandeurs sur les accidentés de la route

Hier soir, Benoît m’a posé une colle par courriel, et je dois avouer que j’ai horreur de ça… Il me demandait un ordre de grandeur descoûts moyens en France (du point de vue de l’assureur) d’un blessé et d’un mort pour les sinistres automobiles“. La colle m’a embarrassé parce que déteste ne pas avoir la moindre idée des ordres de grandeurs sur des sujets qui pourtant m’intéressent. En fait, je connais les coûts extrêmes, les plus gros accidents, mais pas les coûts moyens… Moralité je me suis mis au travail, et comme toujours sur ces sujets, Audrey m’a beaucoup aidé. Alors avant de commencer, il faut savoir que quand on parle de gravité d’un accident corporel (dans les précédant billets, ici ou , je ne parlais que des nombres) il faut parler de taux d’AIPP, atteinte permanente à l’intégrité physique et/ou psychique (anciennement appelées taux d’IPP, i.e. taux d’incapacité permanente partielle). Alors c’est très technique, c’est fixé par des docteurs, des caisses d’assurance maladie, des juges, bref, il y a tellement de paramètres que c’est dur de s’y retrouver. Sauf qu’il y a des espèces de grilles… Elles sont évoquées ici par exemple,

  • Perte, paralysie d’un doigt ou perte de l’odorat ou d’une partie de la voix =  1 à 5 %
  • Perte de l’ ouie d’une oreille ou d’un orteil genou remplacé par une prothèse =  5 à 10 %
  • Perte de tous les orteils, perte de toutes les dents, perte de la voix =  10 à 15 %

(bref, un taux de 10% c’est pas très drôle.. mais ce n’est pas ce qui est retenu pour définir un accident grave)

  • Perte du pouce droit (pour un droitier) ou perte de la vision d’un œil =  20 à 25 %
  • Perte de la main droite (pour un droitier) ou perte du bras gauche =  40 à 50 %

Là on va franchir le cap des 50% qui correspondent à une classification particulière, à savoir les graves,

  • Perte du bras droit (pour un droitier), perte de la jambe ou perte totale de la vue ou de l’ouïe =  50 à 60 %
  • Para ou tétraplégie et  traumatisme crâniens =  60 à 100 %

Voilà pour les ordres de grandeurs sur la gravité d’un accidenté.
En 2008, il y a eu de l’ordre de 7 millions de sinistres (avec une indemnisation de l’assureur). Près de la moitié concernent des bris de glaces, et environ 250 000 correspondaient à un accident corporel. Si on creuse un peu, parmi les accidentés, dans 75% des accidents, il s’agit de blessés légers, sans AIPP. 24% sont accidentés avec un AIPP, avec la répartition suivante: 19% ont un taux d’AIPP inférieur à 5%, et 2% dépassent 20%. Enfin, 1% décèdent. Voilà pour la répartition des nombres, regardons maintenant les coûts.
Ces 3% des sinistres (en nombres) représentent 30% du coût total des sinistres, soit 4,5 milliards d’euros (vu à la fin de la première année, selon la FFSA). Et ces milliards d’euros sont répartis de la manière suivante

  • 11% pour les blessés léger
  • 36% pour les AIPP compris entre 1% et 20%
  • 42% pour les AIPP supérieurs à 20%
  • 11% pour les décès

Si on regarde l’indemnité moyenne en 2008 (pour en venir à l’objet de mon billet)

  • les blessés avec AIPP ont été indemnisé à hauteur de 45,680 euros (dont 14,165 euros au titre des dépenses de santé, actuelles et projetées)
  • les ayants droits, en cas de décès, ont été indemnisé à hauteur de 89,800 euros,

Si l’on regarde les blessés légers (sans AIPP), ils ont touché en moyenne 1,500 euros (la médiane se situant plutôt vers 750 euros).
Pour les blessés avec AIPP

  • taux entre 1% et 5%, 9,195 euros
  • taux entre 6% et 9%, 30,828 euros
  • taux entre 10% et 14%, 53,295 euros,
  • taux entre 15% et 19%, 86,950 euros,
  • taux entre 20% et 29%, 130,585 euros,
  • taux entre 20% et 29%, 130,585 euros,
  • taux supérieur à 30%, 644,100 euros,
  • taux supérieur à 50%, 988,000 euros,
  • taux supérieur à 80%, 2,018,000 euros,

Si l’on regarde plus en détail les blessés graves (taux d’AIPP supérieur à 50%) l’indemnité total était proche du million d’euros, répartis de la manière suivante

  • 21% au titre du déficit fonctionnel permanent
  • 6% au titre de préjudices personnels,
  • 25% au titres des dépenses de santé (actuelles et projetées)
  • 5% au titre de pertes de gains professionnels
  • 35% au titre d’assistance par une tierce personne

(le reste étant d’autres frais). De plus, l’indemnité décroît avec l’âge

  • moins de 30 ans, 1,521,000 euros
  • entre 30 et 60 ans, 1,287,000 euros,
  • plus de 60 ans, 580,000 euros,

Comme le note la FFSA, l’évolution de l’indemnité moyenne croît régulièrement depuis 10 ans, avec une croissance annuelle de l’ordre de 9%.

Will I ever be a bayesian statistician ? (part 1)

Last week, during the workshop on Statistical Methods for Meteorology and Climate Change (here), I discovered how powerful bayesian techniques could be, and that there were more and more bayesian statisticians. So, if I was to fully understand applied statisticians in conferences and workshops, I really have to understand basics of bayesian statistics. I have published some time ago some posts on bayesian statistics applied to actuarial problems (here or there), but so far, I always thought that bayesian was a synonym for magician. To be honest, I am a Muggle, and I have not been trained as a bayesian. But I can be an opportunist…

So I decided to publish some posts on bayesian techniques, in order to prove that it is actually not that difficult to implement.

As far as I understand it, in bayesian statistics, the parameter is considered as a random variable (which is also the case, in classical mathematical statistics). But here, here assume that this parameter does have a parametric distribution….
Consider a classical statistical problem: assume we have a sample http://freakonometrics.free.fr/blog/bayy1.png i.i.d. with distribution http://freakonometrics.free.fr/blog/bayy2.png. Here we note

http://freakonometrics.free.fr/blog/bayy3.png

since parameter http://freakonometrics.free.fr/blog/bayyyyy001.png is a random variable. The idea is to assume that http://freakonometrics.free.fr/blog/bayyyyy001.png has a (so called a priori) distribution, e.g.

http://freakonometrics.free.fr/blog/bayy4.png

So far it was simple. The idea is then to consider the posterior distribution of http://freakonometrics.free.fr/blog/bayyyyy001.png, given the observations http://freakonometrics.free.fr/blog/bayyyyyy02.png. Thus, we need to compute the distribution of http://freakonometrics.free.fr/blog/bayyyyyy03.png which is here extremely simple (due to properties of the Gaussian distribution), i.e.

http://freakonometrics.free.fr/blog/bayyyyyy04.png

where

http://freakonometrics.free.fr/blog/bayyyyyy05.png

And them, it becomes extremely natural to consider http://freakonometrics.free.fr/blog/bayy20.png as an estimator of given our sample data (and thus, we also have a confidence interval since we know the distribution of http://freakonometrics.free.fr/blog/bayyyyy001.png given the observations http://freakonometrics.free.fr/blog/bayyyyyy02.png).
In order to be sure that we understood, consider now a heads and tails problem, i.e. http://freakonometrics.free.fr/blog/bayy5.png. Note, first, that \theta has support http://freakonometrics.free.fr/blog/bayy6.png. So we need a distribution on that support. Why not a beta distribution ? E.g.

http://freakonometrics.free.fr/blog/bayy7.png

Thus,

http://freakonometrics.free.fr/blog/bayy8.png

and

http://freakonometrics.free.fr/blog/bayy9.png

From Bayes formula,

http://freakonometrics.free.fr/blog/bayy10.png

and we get easily

http://freakonometrics.free.fr/blog/bayy11.png

which is the density of a Beta distribution, i.e.

http://freakonometrics.free.fr/blog/bayy12.png
prior=dbeta(u,a,b)
posterior=dbeta(u,a+y,n-y+b)

The estimator proposed is then the expected value of that conditional distribution,

http://freakonometrics.blog.free.fr/public/perso/bayyyyyyyyyyy.png

Note that

http://freakonometrics.free.fr/blog/bayy13.png

Further, it is possible to derive confidence intervals using quantiles of the posterior distribution.
On the graphs below, we consider the following heads/tails sample

A first idea is to consider a uniform prior distribution.

http://freakonometrics.free.fr/blog/bayes-cv-1.gif

A second idea is to consider an asymmetric beta distribution. First, with an asymmetry on the left,

http://freakonometrics.free.fr/blog/bayes-cv-3.gif

or on the right
http://freakonometrics.free.fr/blog/bayes-cv-2.gif

Finally a third idea is simply to get back to the standard Gaussian approximation,

http://freakonometrics.free.fr/blog/bayes-cv-gauss.gif

If we compare the four models, we obtain (the plain black line is the Gaussian approximated distribution for the empirical mean), and red lines are obtained from prior beta distributions

http://freakonometrics.free.fr/blog/bayes-cv-all.gif

The code to generate those graphs is the following
a1=1; b1=1
D1[1,]=dbeta(u,a,b)
a2=4; b2=2
D2[1,]=dbeta(u,a,b)
a3=2; b3=4
D3[1,]=dbeta(u,a,b)
setseed(1)
S=sample(0:1,size=100,replace=TRUE)
COULEUR=rev(rainbow(120))
D1=D2=D3=D4=matrix(NA,101,length(u))
for(s in 1:100){
y=sum(S[1:s])
D1[s+1,]=dbeta(u,a1+y,s-y+b1)
D2[s+1,]=dbeta(u,a2+y,s-y+b2)
D3[s+1,]=dbeta(u,a3+y,s-y+b3)
D4[s+1,]=dnorm(u,y/s,sqrt(y/s*(1-y/s)/s))
plot(u,D1[1,],col="black",type="l",ylim=c(0,8),
xlab="",ylab="")
for(i in 1:s){lines(u,D1[1+i,],col=COULEUR[i])}
points(y/s,0,pch=3,cex=2)
plot(u,D2[1,],col="black",type="l",ylim=c(0,8),
xlab="",ylab="")
for(i in 1:s){lines(u,D2[1+i,],col=COULEUR[i])}
points(y/s,0,pch=3,cex=2)
plot(u,D3[1,],col="black",type="l",ylim=c(0,8),
xlab="",ylab="")
for(i in 1:s){lines(u,D3[1+i,],col=COULEUR[i])}
points(y/s,0,pch=3,cex=2)
plot(u,D4[1,],col="white",type="l",ylim=c(0,8),
xlab="",ylab="")
for(i in 1:s){lines(u,D4[1+i,],col=COULEUR[i])}
points(y/s,0,pch=3,cex=2)
plot(u,D4[s+1,],col="black",lwd=2,type="l",
ylim=c(0,8),xlab="",ylab="")
lines(u,D1[1+i,],col="blue")
lines(u,D2[1+i,],col="red")
lines(u,D3[1+i,],col="purple")
points(y/s,0,pch=3,cex=2)
}

Here, we can see that computations are simple if the prior distribution has a distribution which is the conjugate of the observations’ distribution (see here for the list of prior and posterior standard distributions).
So far, I have two questions that naturally show up

  • is it possible to start with a neutral prior distribution, non informative ?
  • what if we are no longer working with conjugate distributions ?

Well, I guess I have to work a bit more to answer those questions…. to be continued

Convergences de variables aléatoires

Le schéma – ou diagramme – des relations entre les modes convergences ressemble à ce qui est présenté sur la droite. C’est grosso modo ce que j’ai dit en cours (sans le représenter explicitement).

Je ne vais pas démontrer ces implications car ce sont des théorèmes (qui peuvent donc se trouver dans n’importe quel textbook de probabilité).

En revanche, un étudiant m’avait demandé si les réciproques étaient vraies. Comme toujours en maths, si les réciproques étaient toujours vraies, on le mettrait. Donc ici les réciproques sont généralement fausses. On peut toutefois se demander quand ces réciproques peuvent être vraies…et quand elles sont fausses.

  • de la convergence en proba à la convergence en moyenne

Commençons par le plus simple (à mon avis). Il existe un cas particulier où cette réciproque est vraie: le cas de variables aléatoires bornées. Autrement dit, si http://freakonometrics.free.fr/blog/cvvvv20.png et que l’on a http://freakonometrics.free.fr/blog/cvvvvv22.png pour tout http://freakonometrics.free.fr/blog/lgn10.png, alors http://freakonometrics.free.fr/blog/cvvvvvv23.png pour tout http://freakonometrics.free.fr/blog/cvvvvv24.png.
Bon, en fait la condition est un peu trop forte, et il suffit que la suite http://freakonometrics.free.fr/blog/cvvv1.png soit uniformément intégrable.
Maintenant, il existe aussi des cas où cette réciproque sera fausse (et on se doute qu’en construisant une suite non bornée, on devrait y arriver). Par exemple, supposons que pour tout http://freakonometrics.free.fr/blog/lgn10.png,
http://freakonometrics.free.fr/blog/cvvv5.png
On peut montrer que http://freakonometrics.free.fr/blog/cvvvvv7.png, car en effet
http://freakonometrics.free.fr/blog/cvvvv6.png
Toutefois, certes la suite converge dans L1, mais vers une limite différente, étant donné que http://freakonometrics.free.fr/blog/cvvvvv9.png pour tout http://freakonometrics.free.fr/blog/lgn10.png. De plus, http://freakonometrics.free.fr/blog/cvvvv10.png et donc la suite ne converge pas dans L2.

  • de la convergence en loi à la convergence en proba

Un exemple simple ici: si http://freakonometrics.free.fr/blog/cvvv1.png converge en loi vers une constante, alors on aura convergence en probabilité (vers la même constante).
Maintenant la réciproque est relativement simple, l’idée étant que http://freakonometrics.blog.free.fr/public/perso/cvvv555.png n’implique en rien http://freakonometrics.blog.free.fr/public/perso/cvvvvv5666.png, au sens où http://freakonometrics.blog.free.fr/public/perso/cvvvvvv56666.png, pour tout http://freakonometrics.blog.free.fr/public/perso/cvvvvv56565.png. Supposons maintenant que l’on ait à la fois http://freakonometrics.blog.free.fr/public/perso/cvvvv4646.png et http://freakonometrics.free.fr/blog/cvvvv20.png, où on va supposer que http://freakonometrics.free.fr/wesp30.png est une loi symétrique. La symétrie signifie que http://freakonometrics.blog.free.fr/public/perso/cvvv5656655.png. Notons que l’on a alors forcément http://freakonometrics.blog.free.fr/public/perso/cvvvcv_56.png. Supposons que la réciproque soit toujours vraie, alors on devrait avoir http://freakonometrics.blog.free.fr/public/perso/cv65464565645vv.png, et donc on aurait que http://freakonometrics.blog.free.fr/public/perso/cvwgfrwegr656.png, ce qui n’arrive que si http://freakonometrics.free.fr/wesp30.png est la constante http://freakonometrics.free.fr/blog/cvvvv4.png.

  • de la convergence en proba vers la convergence presque sûre

Je ne vais pas m’étendre sur la construction d’exemples pour lesquels la réciproque est fausse (formellement, il suffit de changer d’espace de probabilité sous-jacent, je peux renvoyer vers les notes de cours de Philippe Briand (ici) et Hélène m’a suggéré également les exercices en ligne ). Par contre, il existe des cas importants où la réciproque est vraie.
Considérons une suite de variables indépendantes et identiquement distribuées http://freakonometrics.free.fr/blog/cvvv1.png. Posons alors http://freakonometrics.free.fr/blog/cvvvv30.png On peut alors montrer que http://freakonometrics.free.fr/blog/cvvvvvv26.png si et seulement si http://freakonometrics.free.fr/blog/cvvvvv35.png (ça doit être du Borel-Cantelli).

  • de la convergence en loi à la convergence presque sûre

Je ne vais pas non plus m’étendre davantage, mais on a un résultat qui ressemble à ce que je disais auparavant. En fait, il existe un théorème (dit de Skorohod) qui permet de dire que si on a une convergence en loi, alors on peut construite une suite de variables aléatoires qui
convergera presque sûrement. Formellement on a que si http://freakonometrics.blog.free.fr/public/perso/cvvvv4646.png, alors il existe une suite http://freakonometrics.free.fr/blog/cvvv11.png et une variable http://freakonometrics.free.fr/blog/cvvv13.png, telles que http://freakonometrics.free.fr/blog/cvvv16.png et http://freakonometrics.free.fr/blog/cvvvvv15.png telles que http://freakonometrics.free.fr/blog/cvvvv17.png. Mais de manière générale, si la suite n’est pas construite comme il faut (http://freakonometrics.blog.free.fr/public/perso/cvvvvv56565.png par http://freakonometrics.blog.free.fr/public/perso/cvvvvv56565.png en quelques sortes) alors il est assez rare d’avoir la réciproque.

Un peu de terminologie statistique

Suite à quelques questions (anodines) posées au fil du cours de statistique, je prends 5 minutes pour refaire le point sur les abréviations qui me semblent classiques, mais qui le ne sont peut-être pas (car je n’ai pas l’habitude d’enseigner en deuxième année, pour ceux qui en doutaient encore).

Tout d’abord des locutions latines. Par exemple i.e. n’est pas Internet Explorer mais vient du latin id est et signifie donc c’est-à-dire. Parmi les petites locutions proches, e.g. vient du latin exempli gratia, et signifie littéralement par la grâce de l’exemple mais on retiendra que cela signifie par exemple.

A l’occasion, il m’arrive d’utiliser ca qui est une abréviation du latin circa, signifiant environ.Et quid signifie en latin quoi. Sinon dans les textbook, q.e.d. vient du latin quod erat demonstrandum, que l’on traduit souvent par ce qu’il fallait démontrer, ou cqfd. Enfin, vs est une contraction du latin versus signifiant comparé à. Bon, je pense avoir fait le tour pour les expressions latines.

Du côté des lettres grecques couramment utilisées en statistique, α i.e. alpha, et β i.e. beta, que l’on verra apparaître dans les erreurs de tests statistiques. Γ est une fonction qu’il pourrait être utile d’utiliser dans les calculs d’intégrales (correspondant à une généralisation de la factorielle). δ i.e. delta, apparaît dans la notation de la masse de Dirac. ε i.e. epsilon que l’on retrouve un peu partout en analyse pour dénoter une constante qui peut être aussi petite qu’on le souhaite. θ i.e. theta qui correspond au paramètre de nos lois, mais parfois on utilisera aussi λ i.e. lambda qui apparaît comme paramètre de la loi de Poisson ou de la loi exponentielle. μ i.e. mu qui désigne généralement la moyenne d’une variable aléatoire, ν i.e. nu qui désigne les degrés de liberté dans les lois de Student ou de Fisher. ρ i.e. rho qui désigne parfois la corrélation. σ i.e. sigma qui désigne l’écart-type d’une variable aléatoire. χ i.e. khi (ou chi) qui apparaît souvent avec une puissance 2, dans la loi du chi-deux. Et enfin ω i.e. omega qui désigne un état de la nature, appartenant à Ω (la lettre majuscule).

Voilà pour le premier point de terminologie. Je me laisse la possibilité de poursuivre, en particulier sur les abréviations plus mathématiques (comme les p.s. ou en anglais a.s., ou autres w.r.t.). Toutes les suggestions sont les bienvenues…

Catastrophes naturelles et population mondiale

Pour reprendre une question posée au téléphone par François, et par Eric dans un commentaire (ici), il est possible de comparer l’évolution du nombre de catastrophes et de l’évolution de la population (mondiale). Pour les données de population on peut en trouver sur wikipedia (ici), sur le site des nations unies () après 1950, ou encore sur le site du census.gov (). Ensuite on peut moyenner les doublons.

POP=read.table("http://freakonometrics.blog.free.fr/public/data/pop-mondial.txt",
sep=";",header=TRUE)
population=tapply(POP$POP,POP$YEAR,mean)
years=as.numeric(names(population))
reg=lm(population~bs(years,3))
lines(years,predict(reg),col="red")
plot(years,population)
baseline=predict(reg,newdata=data.frame(years=1900:2010))

Revenons à nos catastrophes “naturelles“. Par exemple les tremblements de terre. La première piste (évoquée ) était de lisser cette série pour l’utiliser en baseline (multiplicatif). On retrouve ici une série qui oscille autour de 1.

earthquake=read.table("http://freakonometrics.blog.free.fr/public/data/DataSet.csv",
sep=",",header=TRUE)
flood=read.table("http://freakonometrics.blog.free.fr/public/data/DataSet_2_.csv",
sep=",",header=TRUE)
storm=read.table("http://freakonometrics.blog.free.fr/public/data/DataSet_3_.csv",
sep=",",header=TRUE)
library(splines)
n=nrow(earthquake)
earthquake=earthquake[-((n-3):n),]
DATE=1900:2010
EARTHQUAKE=FLOOD=STORM=rep(0,length(DATE))
EARTHQUAKE[as.numeric(as.character(earthquake$Year))-1899]=
as.numeric(as.character(earthquake$Number))
plot(DATE,EARTHQUAKE,type="b")
reg=lm(EARTHQUAKE~bs(DATE,6))
lines(DATE,predict(reg),col="red")
plot(DATE,EARTHQUAKE/predict(reg),type="b")
abline(h=1,lty=2,col="blue")

Si on utilise cette fois la population mondiale comme baseline, on obtient

POP=read.table("http://freakonometrics.blog.free.fr/public/data/pop-mondial.txt",
sep=";",header=TRUE)
population=tapply(POP$POP,POP$YEAR,mean)
years=as.numeric(names(population))
reg=lm(population~bs(years,3))
lines(years,predict(reg),col="red")
plot(years,population)
baseline=predict(reg,newdata=
     data.frame(years=1900:2010))/1000000
plot(DATE,EARTHQUAKE/baseline,type="b")

La série de points présente une jolie rupture. Afin de trouver la date optimale de rupture, utilisons un test basé sur le critère d’Akaike.

Y=EARTHQUAKE/baseline
X=DATE
DUMMY=1910:2000
VAIC=rep(NA,length(DUMMY))
for(x in DUMMY){
IX=(X>x)
reg=lm(Y~1+IX)
VAIC[x-1909]=AIC(reg)
}
plot(DUMMY,VAIC)
k=which.min(VAIC)
DUMMY[k]

Manifestement, il y a une rupture en 1975 (que je ne sais trop comment interpréter).

Si on regarde une autre série, comme celle des inondations.

n=nrow(flood)
flood=flood[-((n-3):n),]
DATE=1900:2010
FLOOD[as.numeric(as.character(flood$Year))-1899]=
as.numeric(as.character(flood$Number))
plot(DATE,FLOOD,type="b")
reg=lm(FLOOD~bs(DATE,6))
lines(DATE,predict(reg),col="red")
population=tapply(POP$POP,POP$YEAR,
mean)
years=as.numeric(names(population))
reg=lm(population~bs(years,3))
baseline=predict(reg,newdata=data.frame(years=1900:2010))
plot(DATE,FLOOD/baseline,type="b")

On retrouve ici clairement le fait que la croissance de la population ne peut expliquer la croissance du nombre d’inondations.

Échantillonnage, TCL et LGN

Mercredi matin, nous allons aborder l’échantillonnage. Nous allons partir d’un échantillon http://freakonometrics.free.fr/blog/lgn01.png, supposé indépendant et identiquement distribué. Le modèle sous-jacent est la collection de variables aléatoires http://freakonometrics.free.fr/blog/lgn02.pngindépendantes et identiquement distribuées. Dans un premier temps, on s’intéressera à la moyenne de l’échantillon, correspondant à une réalisation de la variable aléatoire http://freakonometrics.free.fr/blog/lgn04.png définie par

http://freakonometrics.free.fr/blog/lgn05.png

Afin d’avoir des résultats sur le comportement de la moyenne, on rappellera deux théorèmes importants: la loi des grands nombres et le théorème central limite. Afin d’illustrer ces théorèmes, considérons le cas de tirages de PILE et FACE. On suppose que http://freakonometrics.free.fr/blog/lgn06.png désigne la variable aléatoire correspondant au fait d’avoir obtenu PILE au http://freakonometrics.free.fr/blog/lgn07.pngème tirage. On remet des pièces à une centaine d’élèves (chaque élève correspond à une couleur différente), et on regarde la moyenne du nombre de PILE obtenus sur les http://freakonometrics.free.fr/blog/lgn10.png premiers lancers. Et on regarde ce qui se passe quand n augmente. La loi de grands nombres permet d’affirmer que

http://freakonometrics.free.fr/blog/lgn11.png

ici, si l’on suppose que la pièce n’est pas biaisée, soit

http://freakonometrics.free.fr/blog/lgn12.png

pour la version faible, ou pour la version forte

http://freakonometrics.free.fr/blog/lgn13.png

soit

http://freakonometrics.free.fr/blog/lgn16.png

C’est ce que l’on observe ci-dessous, chaque trait correspondant à la trajectoire http://freakonometrics.free.fr/blog/lgn21.png obtenue par un élève. Un histogramme est représenté à droite.

http://freakonometrics.free.fr/animation-cv-binomiale.gif

On peut aussi centrer (ici retrancher http://freakonometrics.free.fr/blog/lgn22.png) et réduire (ici mulplier par http://freakonometrics.free.fr/blog/lgn23.png) la moyenne empirique pour étudier une éventuelle convergence. On pose alors

http://freakonometrics.free.fr/blog/lgn24.png

Le théorème central limite garantie que

http://freakonometrics.free.fr/blog/lgn27.png

On a toujours l’histogramme des valeurs obtenues, au bout de http://freakonometrics.free.fr/blog/lgn10.png tirages, par les élèves, avec en rouge la densité de la loi normale centrée réduire. Et là encore, l’expérience semble valider la théorie…

http://freakonometrics.free.fr/animation-cv-CR.gif

Maintenant, pour ceux qui doutent, la base est en ligne ici, et je mets le code pour commencer à jouer avec…

M=read.table("http://freakonometrics.blog.free.fr/public/data/PILEFACE.txt")
nl=ncol(M)
nj=nrow(M)
M1=matrix(M=="PILE",nj,nl)
MS=M1
for(l in 2:nl){MS[,l]=
(apply(M1[,1:l],1,mean)-.5)/sqrt(.25/l)}
arret=100
COULEUR= rev(rainbow(nj))
HIST=hist(MS[,arret])
plot(1:nl,MS[1,],ylim=(c(-3,3)),col="white",
xlim=c(0,250),xlab="Nombre de lancers de pièces",
ylab="Nombre moyen de PILE, centré réduit")
for(j in 1:nj){
lines(1:arret,MS[j,1:arret],col=COULEUR[j])
}
abline(v=arret)

Increasing number of catastrophes and human factor

Peter asked me (here) to discuss the occurrence of natural catastrophes. I have looked at the EM-DAT dataset (here), and more specifically three kinds of risks: earthquakes, floods and storms (in any country over the past century). For earthquakes, we have the following figures,

the blue points are raw data (I mean, public ones, those that can be downloaded on the website), and I added two smoothed patterns, in red (based on splines regressions). This is no reason to assume that there could be an increase of the number of earthquakes. The trend is due to the fact that we look at earthquakes considered as catastrophic, i.e. that caused deaths and severe damages. And as pointed out in Peter’s comment, it means that “the increase in ‘disasters’ might have a lot to do with the development of the population”.
If we look now at storms, we have

i.e. also an increasing trend. But here, we can consider that the trend can be explained by two factors: an inflation due to the value at risk (e.g. increase of concentration of population in some areas as for earthquakes) but maybe also a change that can be related to climate change.  If we remove the pattern du to the increase of value at risk on the raw data (i.e. we divide raw data by the trend observed on earthquakes, the red line), we have

Here it looks like the number of catastrophic storms is increasing (removing the human factor). And it might be possible to relate that trend to climate change. If we look finally at major flood events, we also observe an increase of the pattern,

and we can divide by the trend observed on earthquakes, to obtain

Looking at that graph, it is possible to propose another interpretation about the increase: people might have a shorter memory about rivers than the one they have on earthquakes; and that they think they might influence climate risk, and flood risk, (but not physical ones). I mean that people know where earthquake risk is, and did not really concentrate in those areas during the past century. On the other hand, it is possible to assume that people imagine they can change river design, e.g. building dykes, and though they can control the risk. But obviously they couldn’t.

Warming in Paris: minimas versus maximas ?

Recently, I received comments (here and on Twitter) about my previous graphs on the temperature in Paris. I mentioned in a comment (there) that studying extremas (and more generally quantiles or interquantile evolution) is not the same as studying the variance. Since I am not a big fan of the variance, let us talk a little bit about extrema behaviour.

In order to study the average temperature it is natural to look at the linear (assuming that it is linear, but I proved that it could reasonably be assumed as linear in the paper) regression, i.e. least square regression, which gives the expected value. But if we care about extremes, or almost extremes, it is natural to look at quantile regression.

For instance, below, the green line is the least square regression, the red one is 97.5% quantile, and the blue on the 2.5% quantile regression.

It looks like the slope is the same, i.e. extremas are increasing as fast as the average…

tmaxparis=read.table("temperature/TG_SOUID100845.txt",
skip=20,sep=",",header=TRUE)
head(tmaxparis)
Dparis=as.Date(as.character(tmaxparis$DATE),"%Y%m%d")
Tparis=as.numeric(tmaxparis$TG)/10
Tparis[Tparis==-999.9]=NA
I=sample(1:length(Tparis),size=5000,replace=FALSE)
plot(Dparis[I],Tparis[I],col="grey")
abline(lm(Tparis~Dparis),col="green")
library(quantreg)
abline(rq(Tparis~Dparis,tau=.025),col="blue")
abline(rq(Tparis~Dparis,tau=.975),col="red")

(here I plot randomly some points to avoid a too heavy figure, since I have too many observations, but I keep all the observations in the regression !).

Now, if we look at the slope for different quantile level (Fig 6 in the paper, here, but on minimum daily temperature, here I look at average daily temperature), the interpretation is different.

s=0
COEF=SD=rep(NA,199)
for(i in seq(.005,.995,by=.005)){
s=s+1
REG=rq(Tparis~Dparis,tau=i)
COEF[s]=REG$coefficients[2]
SD[s]=summary(REG)$coefficients[2,2]
}

with the following graph below,

s=0
plot(seq(.005,.995,by=.005),COEF,type="l",ylim=c(0.00002,.00008))
for(i in seq(.005,.995,by=.005)){
s=s+1
segments(i,COEF[s]-2*SD[s],i,COEF[s]+2*SD[s],col="grey")
}
REG=lm(Tparis~Dparis)
COEFlm=REG$coefficients[2]
SDlm=summary(REG)$coefficients[2,2]
abline(h=COEFlm,col="red")
abline(h=COEFlm-2*SDlm,lty=2,lw=.6,col="red")
abline(h=COEFlm+2*SDlm,lty=2,lw=.6,col="red")

Here, for minimas (quantiles associated to low probabilities, on the left), the trend has a higher slope than the average, so in some sense, warming of minimas is stronger than average temperature, and on other hand, for maximas (high probabilities on the right), the slope is smaller – but positive – so summer are warmer, but not as much as winters.
Note also that the story is different for minimal temperature (mentioned in the paper) compared with that study, made here on average daily temperature (see comments)… This is not a major breakthrough in climate research, but this is all I got…

More climate extremes, or simply global warming ?

In the paper on the heat wave in Paris (mentioned here) I discussed changes in the distribution of temperature (and autocorrelation of the time series).

During the workshop on Statistical Methods for Meteorology and Climate Change today (here) I observed that it was still an important question: is climate change affecting only averages, or does it have an impact on extremes ? And since I’ve seen nice slides to illustrate that question, I decided to play again with my dataset to see what could be said about temperature in Paris.
Recall that data can be downloaded here (daily temperature of the XXth century).

tmaxparis=read.table("/temperature/TX_SOUID100124.txt",
skip=20,sep=",",header=TRUE)
Dmaxparis=as.Date(as.character(tmaxparis$DATE),"%Y%m%d")
Tmaxparis=as.numeric(tmaxparis$TX)/10
tminparis=read.table("/temperature/TN_SOUID100123.txt",
skip=20,sep=",",header=TRUE)
Dminparis=as.Date(as.character(tminparis$DATE),"%Y%m%d")
Tminparis=as.numeric(tminparis$TN)/10
Tminparis[Tminparis==-999.9]=NA
Tmaxparis[Tmaxparis==-999.9]=NA
annee=trunc(tminparis$DATE/10000)
MIN=tapply(Tminparis,annee,min)
plot(unique(annee),MIN,col="blue",ylim=c(-15,40),xlim=c(1900,2000))
abline(lm(MIN~unique(annee)),col="blue")
abline(lm(Tminparis~unique(Dminparis)),col="blue",lty=2)
annee=trunc(tmaxparis$DATE/10000)
MAX=tapply(Tmaxparis,annee,max)
points(unique(annee),MAX,col="red")
abline(lm(MAX~unique(annee)),col="red")
abline(lm(Tmaxparis~unique(Dmaxparis)),col="red",lty=2)

On the plot below, the dots in red are the annual maximum temperatures, while the dots in blue are the annual minimum temperature. The plain line is the regression line (based on the annual max/min), and the dotted lines represent the average maximum/minimum daily temperature (to illustrate the global tendency),

It is also possible to look at annual boxplot, and to focus either on minimas, or on maximas.

annee=trunc(tminparis$DATE/10000)
boxplot(Tminparis~as.factor(annee),ylim=c(-15,10),
xlab="Year",ylab="Temperature",col="blue")
x=boxplot(Tminparis~as.factor(annee),plot=FALSE)
xx=1:length(unique(annee))
points(xx,x$stats[1,],pch=19,col="blue")
abline(lm(x$stats[1,]~xx),col="blue")
annee=trunc(tmaxparis$DATE/10000)
boxplot(Tmaxparis~as.factor(annee),ylim=c(15,40),
xlab="Year",ylab="Temperature",col="red")
x=boxplot(Tmaxparis~as.factor(annee),plot=FALSE)
xx=1:length(unique(annee))
points(xx,x$stats[5,],pch=19,col="red")
abline(lm(x$stats[5,]~xx),col="red")

Plain dots are average temperature below the 5% quantile for minima, or over the 95% quantile for maxima (again with the regression line),

We can observe an increasing trend on the minimas, but not on the maximas !
Finally, an alternative is to remember that we focus on annual maximas and minimas. Thus, Fisher and Tippett theory (mentioned here) can be used. Here, we fit a GEV distribution on a blog of 10 consecutive years. Recall that the GEV distribution is

http://freakonometrics.blog.free.fr/public/perso/gev1.png
install.packages("evir")
library(evir)
Pmin=Dmin=Pmax=Dmax=matrix(NA,10,3)
for(s in 1:10){
X=MIN[1:10+(s-1)*10]
FIT=gev(-X)
Pmin[s,]=FIT$par.ests
Dmin[s,]=FIT$par.ses
X=MAX[1:10+(s-1)*10]
FIT=gev(X)
Pmax[s,]=FIT$par.ests
Dmax[s,]=FIT$par.ses
}

The location parameter http://freakonometrics.blog.free.fr/public/perso/gev4.png is the following, with on the left the minimas and on the right the maximas,

while the scale parameter http://freakonometrics.blog.free.fr/public/perso/gev3.png is

and finally the shape parameter http://freakonometrics.blog.free.fr/public/perso/gev2.png is

On those graphs, it is very difficult to say anything regarding changes in temperature extremes… And I guess this is a reason why there is still active research on that area…

Cursed numbers ?

In Lost, Hugo “Hurley” Reyes played the numbers 4, 8, 15, 16, 23 and 42 at the lottery, and ended up winning the $114-million jackpot. And over the ensuing weeks, everyone around him seems to suffer increasingly bad luck: Hurley’s grandfather dies of a heart attack, his brother’s wife walks out, his mother breaks her ankle while the house Hurley bought her goes up in flames, and Hurley himself is falsely arrested.
Anyway, last week (here) 4 numbers (out of 6) appeared at the lottery in LA. As pointed out by Xian (here), the odds were not that small, i.e. it is a 1‰ chance,

http://freakonometrics.blog.free.fr/public/perso/probaloto.png
Hence, with one lottery per week, the return period is 16 years. Note this percentage is very close to what we did observe on the French lottery (below the statistics in ‰, from here, in a zip file),

> loto=read.table("loto.csv",dec=",",header=TRUE,sep=";")
> ntirage=nrow(loto)
> loto=loto[51:ntirage,]
> ntirage=nrow(loto)
> N=as.matrix(loto[,c("boule_1","boule_2","boule_3",
   "boule_4","boule_5","boule_6")])
> P=rep(NA,nrow(N))
> for(s in 1:nrow(N)){
+ P[s]=sum(N[s,1]%in%c(4, 8, 15, 16, 23, 42)+
+          N[s,2]%in%c(4, 8, 15, 16, 23, 42)+
+          N[s,3]%in%c(4, 8, 15, 16, 23, 42)+
+          N[s,4]%in%c(4, 8, 15, 16, 23, 42)+
+          N[s,5]%in%c(4, 8, 15, 16, 23, 42)+
+          N[s,6]%in%c(4, 8, 15, 16, 23, 42))
+ }
> table(P)/nrow(N)*1000
P
         0          1          2          3          4
435.732113 405.366057 137.271215  19.966722   1.663894

But what about the full sequence…? Imagine that in France, at the official lottery, the exact sequence played by Hugo appears. What a coincidence.The probability that the sequence appears, assuming that there are 48 possible numbers in the lottery, is

http://freakonometrics.blog.free.fr/public/perso/lotto2.png
i.e. the expected number of draws we need before seeing that sequence for the first time is almost a billion.
Now if we look at all official lotteries around the world, say 100 per week, what is the probability to see Hurley’s sequence shows up – at least once – in 25 years (assuming that after 25 years, no one will remember Lost and those cursed numbers) ? It looks like it is a 1% chance…
http://freakonometrics.blog.free.fr/public/perso/proba-lost.png
So let us wait and see…

Evolution of the number of natural catastrophes

Recently (here and there) I mentioned Reinsurers comments on the increase on the number of catastrophes. But actually, there are other sources of information about catastrophes. For instance the http://www.cred.be/, Centre for Research on the Epidemiology of Disasters. As the mention it below, there has been an increase in the number of natural catastrophes. For a disaster to be entered into the so-called EM-database at least one of the following criteria must be fulfilled:

• Ten (10) or more people reported killed.
• Hundred (100) or more people reported affected.
• Declaration of a state of emergency.
• Call for international assistance.

Note that it is possible to download data (here), so I guess I will work further with that database. To be continued…

Lecture de blog(s): quelques précautions

En finissant mon billet sur le million de pages vues en deux ans sur ce blog, je me suis rendu compte de l’énormité de la chose…. et je me suis dit qu’il convenait de faire une mise en garde (ça en fait des lecteurs potentiels pour mes élucubrations). L’écriture d’un blog est un exercice particulier. Contrairement à une page ouebe régulièrement mise à jour, dans un blog, on ne revient pas en arrière.

Ou plutôt pas complètement….Et les billets s’accumulent. Si un billet contient des fautes, je peux le reprendre, et je raye les fautes. Mais généralement, des lecteurs font des remarques qui éclairent le billet, et le font gagner en clarté. Contrairement aux blogs très fréquentés où les commentateurs sont là pour critiquer, souvent de manière malsaine, ou pour discuter entre eux, il y a peu de commentaires sur mon blog, et souvent, je trouve qu’ils corrigent bien mes imprécisions ou mes coquilles (comme ici).

Moralité, je recommande aux lecteurs intéressés de ne pas se contenter des billets, mais d’aller faire un tour sur les commentaires. Voilà, maintenant au boulot, il y a des billets à finir !