Category Archives: Actuarial science

Il y a des jours avec, et des jours sans

Je ne sais pas si ça vous est déjà arrivé, mais quand on fait des crêpes, la probabilité d’en réussir une est proportionnelle à la probabilité que l’on vient de réussir les dernières. Autrement dit si j’ai réussi les premières, je suis dans un bon jour, et je réussi davantage celles qui suivent. Alors que si je commence par rater les premières, je suis dans un mauvais jour, et je les rate toutes (pour les sportifs, il est aussi possible de noter que, par exemple, au tennis, la probabilité de réussir son service est d’autant plus forte que j’ai réussi les derniers: c’est pareil).

  • Une première formalisation du problème

Supposons que sur les http://freakonometrics.blog.free.fr/public/perso3/crepes-01.gif premiers essais (qui servent d’étalonnage), http://freakonometrics.blog.free.fr/public/perso3/crepes-02.gif fois j’ai réussi ma crêpe, et http://freakonometrics.blog.free.fr/public/perso3/crepes-03.gif fois, je l’ai ratée (trop de trous, trop irrégulière, trop cuite, que sais-je encore). Mon a priori sur la probabilité de ne pas avoir de chance (et de rater ma crêpe) est http://freakonometrics.blog.free.fr/public/perso3/crepes-04.gif. Supposons que la probabilité de réussir à l’instant http://freakonometrics.blog.free.fr/public/perso3/crepes-06.gif (on commence à compter après les http://freakonometrics.blog.free.fr/public/perso3/crepes-01.gif premiers essais) soit la fréquence empirique observée sur le passé. Autrement dit

http://freakonometrics.blog.free.fr/public/perso3/crepes-05.gif

http://freakonometrics.blog.free.fr/public/perso3/crepees-07.gif est le nombre de succès. En particulier http://freakonometrics.blog.free.fr/public/perso3/crepes-09.gif. On retrouve ici le fait que plus j’en ai raté, plus j’ai de chances de rater les prochaines… c’est la mauvaise main en quelque sorte.
Afin de mieux comprendre, supposons http://freakonometrics.blog.free.fr/public/perso3/crepes-10.gif. Autrement dit, sur mes deux premiers essais (qui ne servent qu’à étalonner), j’ai raté et réussi une fois. Juste avant l’instant 1, j’ai deux possibilité: réussir ou rater, avec la probabilité http://freakonometrics.blog.free.fr/public/perso3/crepes-11.gif. I.e.

http://freakonometrics.blog.free.fr/public/perso3/crepes-12.gif

Si on se place une crêpe plus loin, et que l’on regarde les probabilités conditionnelles

http://freakonometrics.blog.free.fr/public/perso3/crepes-13b.gif

avec quelque-chose de similaire si on conditionne par l’autre événement,

http://freakonometrics.blog.free.fr/public/perso3/crepes-14.gif

Et donc la probabilité non conditionnelle s’écrit

http://freakonometrics.blog.free.fr/public/perso3/crepes-15.gif

i.e.

http://freakonometrics.blog.free.fr/public/perso3/crepes-16.gif
http://freakonometrics.blog.free.fr/public/perso3/crepes-18.gif

et

http://freakonometrics.blog.free.fr/public/perso3/crepes-19.gif

Moralité

http://freakonometrics.blog.free.fr/public/perso3/crepes-20.gif

Autrement dit, on poursuivant la récurrence, on a que http://freakonometrics.blog.free.fr/public/perso3/crepees-07.gif est uniformément distribué sur http://freakonometrics.blog.free.fr/public/perso3/crepes-21.gif
Étonnant, non ? Autrement dit, si j’ai commencé par un échec et une réussite, la probabilité d’avoir 50 réussites sur 100 essais sera très exactement http://freakonometrics.blog.free.fr/public/perso3/crepes-23.gif. Qui est la même probabilité que d’avoir 1 réussite, ou 99 !
On a aussi que

http://freakonometrics.blog.free.fr/public/perso3/crepes-22.gif

Mais que se passe-t-il si http://freakonometrics.blog.free.fr/public/perso3/crepes-02.gif et http://freakonometrics.blog.free.fr/public/perso3/crepes-03.gif ne sont plus égaux à 1 ?

  • Lien avec l’approche bayésienne et la théorie de la crédibilité

Ce qui est amusant, c’est que la probabilité

http://freakonometrics.blog.free.fr/public/perso3/crepes-05.gif

est bien connue par les personnes qui connaissent l’assurance et la théorie de la crédibilité. Supposons que je note http://freakonometrics.blog.free.fr/public/perso3/crepes-24.gif la variable indicatrice correspondant à une réussite ou un succès au http://freakonometrics.blog.free.fr/public/perso3/crepes-26.gifème essai (avec probabilité http://freakonometrics.blog.free.fr/public/perso3/crepes-25.gif). D’un point de vu bayésien, la loi conjuguée de la loi binomiale (on va compter les réussites sur n essais) est une loi Beta. Rappelons que la loi Beta http://freakonometrics.blog.free.fr/public/perso3/crepes-27.gif a pour densité

http://freakonometrics.blog.free.fr/public/perso3/crepes-28.gif

sur http://freakonometrics.blog.free.fr/public/perso3/crepes-29.gif, et que l’espérance de cette variable est http://freakonometrics.blog.free.fr/public/perso3/crepes-30.gif, i.e. http://freakonometrics.blog.free.fr/public/perso3/crepes-31.gif. Donc cette modélisation devrait parfaitement convenir à notre problème. Ce qui est intéressant, c’est que comme on travaille avec la loi conjuguée, la loi ex post est encore une loi Beta. Autrement dit, ici, si je suppose que http://freakonometrics.blog.free.fr/public/perso2/testic-07.gif suit a priori une loi Beta http://freakonometrics.blog.free.fr/public/perso3/crepes-27.gif, et que la loi de http://freakonometrics.blog.free.fr/public/perso3/crepees-07.gif sachant http://freakonometrics.blog.free.fr/public/perso2/testic-07.gif suit une loi binomiale http://freakonometrics.blog.free.fr/public/perso3/crepes-32.gif, alors la loi a posteriori de http://freakonometrics.blog.free.fr/public/perso2/testic-07.gif sachant http://freakonometrics.blog.free.fr/public/perso3/crepes-33.gif est une loi http://freakonometrics.blog.free.fr/public/perso3/crepes-34.gif (c’est ce que l’on avait vu ici). Et en particulier

http://freakonometrics.blog.free.fr/public/perso3/crepes-35.gif

Autrement dit on obtient une combinaison convexe entre la moyenne de la loi a priori (à gauche) et la fréquence empirique (à droite).
On retrouve ici ce que nous avions noté auparavant, en notant que http://freakonometrics.blog.free.fr/public/perso3/crepes-36.gif, et http://freakonometrics.blog.free.fr/public/perso3/crepes-37.gif.
On peut toutefois se demander si le résultat que l’on avait obtenu sur la loi de S_n reste valide ici…
Notons que

http://freakonometrics.blog.free.fr/public/perso3/crepes-38.gif

i.e. (moyennant quelques calculs)

http://freakonometrics.blog.free.fr/public/perso3/crepes-40.gif

Dans le cas où http://freakonometrics.blog.free.fr/public/perso3/crepes-10.gif, en réécrivant aussi le premier terme, on a

http://freakonometrics.blog.free.fr/public/perso3/crepes-41.gif

i.e.

http://freakonometrics.blog.free.fr/public/perso3/crepes-42.gif

On retrouve bien la formule obtenue auparavant… Ce qui est intéressant aussi, c’est que si on note que

http://freakonometrics.blog.free.fr/public/perso3/crepes-43.gif

or

http://freakonometrics.blog.free.fr/public/perso3/crepes-44.gif

et donc

http://freakonometrics.blog.free.fr/public/perso3/crepes-45.gif

et par récurrence, on obtient que

http://freakonometrics.blog.free.fr/public/perso3/crepes-46.gif

et que de plus

http://freakonometrics.blog.free.fr/public/perso3/crepes-47.gif

Autrement dit, effectivement, si j’en ai raté beaucoup au début, ma probabilité d’en rater beaucoup est grande (et elle est d’ailleurs constante). Certains parleraient de no-free lunch.

Some stylized facts about large risk covers

A couple of weeks ago, David Cummins (here) was giving a talk in Laval University. And we’ve seen a series of extremely interesting graphs and figures about catastrophe reinsurance market, as well as Cat Bonds prices. The first one was the rate one line index for catastrophe reinsurance (the rate on line is the excess of loss premium expressed as a percentage of the reinsurance cover), from Guy Carpenter (2010, page 10 here).

Following hurricane Andrew in 1992, prices went up quite high. But following hurricane Katrina (which is, so far, the most costly insured disaster following the second World War, with a cost exceeding 70 billions US$ – 2008 $ – while Andrew was only 24 billions – again 2008 $), the bump is much smaller. I though cycles where much larger in the reinsurance industry.

Then there was a discussion about cat bond pricing, with a graph from Lane Financial (2010, page 13, here) with the ratio premium over expected loss

This is extremely interesting, even if it is only about cat bond, and not about reinsurance covers. Usually, when we introduce premium principles in actuarial courses, we start with the pure premium, i.e.

http://freakonometrics.blog.free.fr/public/perso2/chargement-PP-02.gifThen we explain that with such a price, ruin probability is certain (with an infinite time horizon), so we need to add a safety margin, and a standard idea (but that can be criticized since the expected value has – usually – nothing to do with the variability) is to add a loading proportional to the pure premium. Then the premium is

http://freakonometrics.blog.free.fr/public/perso2/chargement-PP-01.gifFor small risks, like motor insurance, the loading is not huge. Actually, if risks have finite variance, it can be obtained simply using the central limit theorem (but I’ll get back on that point in a couple of weeks). Here, we see that loading http://freakonometrics.blog.free.fr/public/perso2/thetaloading.gif can be large (up to 400% in 2009).

An finally an updated graph with a comparison between BB corporate bondscoupon, and BB catastrophe bonds coupon,

(I guess the source is again Morton Lane). I found surprising the recent gap (following Katrina) between the two spreads. I guess financial market started to be scared and understood that catastrophes are not that rare…. I wonder what 2008 and 2009 prices looked like.

Variable annuities is not a systemic risk ?

The Geneva Association just published on its website an interesting report on variable annuities and systemic risk (online here). Based on a definition of potentially systemically risky activities, on interconnectedness or substitutability, the report claims that since “none of the criteria is triggered”, variables annuities is “not a potentially systemically risk activity”. Even if “short-term effects are conceivable”. I guess it is a diplomatic way to say it…

Note that a series of slides can also be downloaded (there) on insurance and systemic risk. But that deserves a more detailed post.

 

Pour une tarification de l’assurance automobile à l’aide du tour de poitrine !

Plusieurs sites spécialisés en assurance commencent à évoquer un arrêté probable de la cour européenne sur la discrimination en assurance (par exemple ici ou ). Une des bases (économiques) de l’assurance est le principe d’Akerlof qui pousse les assureurs à segmenter par classe de risque. Afin de segmenter, et de révéler les classes de risques, on utilise l’historique de sinistralité (information dite a posteriori), ou bien des informations exogènes (dites a priori) sur le conducteur, le véhicule, son usage, etc. Par exemple on peut utiliser l’ancienneté du véhicule, et le nombre de kilomètre effectués (en moyenne) par le conducteur, comme sur le graphique ci-dessous (retrouvé dans les transparents que l’on utilisait avec François Bucchini quand on donnait le cours d’assurance dommage à l’ENSAE, les probas étant “normalisées” dans une espère de base 100)

ou encore le type de carburant utilisé (diesel ou essence)

On retrouve que plus on conduit, plus la probabilité d’avoir un accident augmente, mais le carburant et l’âge du véhicule semblent être aussi des variables discriminantes. Et parmi les variables qui semblent significatives (pour expliquer la probabilité d’avoir un accident), il y a le sexe (croisé ici avec le kilométrage, comme auparavant),

Alors l’effet peut sembler marginal sur ce graphique… mais c’est loin d’être le cas. Par exemple, sans utiliser des techniques très poussées en économétrie, on peut regarder le nombre moyen de sinistres, et le coût moyen de sinistres, par sexe, et par tranche d’âge (voire aussi par CSP et par puissance du véhicule). Dans une étude faite par un assureur, j’avais trouvé les chiffres suivants

En haut à droite (beaucoup d’accidents, et coût – en moyenne – élevé) on retrouve les jeunes hommes. Donc oui, les jeunes hommes sont significativement beaucoup plus risqués que les autres conducteurs. Et le soucis est que, si on ne segmente pas, Georges Akerlof nous explique que le marché de l’assurance disparait, les “bons” risques ne voulant plus payer pour les “mauvais” risques. Sans pour autant rentrer dans une spirale infernale de la segmentation, il est bon que les primes restent corrélées au risque sous-jacent.

Les assureurs prétendent qu’ils ne «ne font pas de la discrimination, ils font de la différenciation ». Je ne rentrerais pas sur les débats de terminologie (pas aujourd’hui en tous les cas), mais le but n’est pas de trouver des variables “explicatives” de la sinistralité au sens causal (malgré la terminologie usuelle des économètres) mais de trouver des variables “corrélées” avec une forte sinistralité, et de les utiliser pour segmenter. Les assureurs européen avaient, jusqu’alors, bénéficié d’un sursis dans le calcul des primes qui leur permet de pratiquer des tarifs différents « lorsque le sexe est un facteur déterminant dans l’évaluation des risques».

Dans un cours d’analyse de données, j’avais montré (ici) qu’à partir des notes de étudiant(e)s à différents examens, je pouvais prédire le sexe des étudiants. Bon, l’étude avait été faite rapidement, avec un petit jeu de données (et donc sans population d’apprentissage et de test), mais il est facile de trouver des variables permettant de deviner le sexe d’un conducteur. D’aucuns pourraient être tentés d’utiliser la pointure des chaussures, mais personnellement je préférerais le tour de poitrine, ou un tour de poitrine ramené à un tour de hanche. Je suis presque sûr qu’avec de telles observations, on peut avoir des variables fortement corrélées avec la survenance d’accident ! En tous les cas ça promet un peu d’animation chez les agents d’assurance ! voire chez les chirurgiens esthétiques (retirer les implants mammaires pour faire baisser sa prime d’assurance auto, voilà qui est original) !

More and more natural catastrophes…

A few weeks ago (here, or there, among many others) I discussed shortly the increase of the frequency of natural catastrophes. An example is perhaps Australia, with major droughts over the past 10 years; a summer heat wave in Victoria, Australia, caused the massive bushfire in 2009; the flooding in Queensland this winter, as well as a huge cyclone. So it looks like Australia is experiencing that increase of the frequency of natural catastrophes.

Nature published last week a series of interesting papers on natural catastrophes (and its relation with a human factor). Increased flood risk linked to global warming by Quirin Schiermeier (likelihood of extreme rainfall may have been doubled by rising greenhouse-gas levels);Climate change: Human influence on rainfall by Richard P. Allan including some Letter by Min and another Letter by Pall; Human contribution to more-intense precipitation extremes by Seung-Ki Min,et al; and Anthropogenic greenhouse gas contribution to flood risk in England and Wales in autumn 2000 by Pardeep Pall et al.

 

Talk at Desjardins General Insurance

This afternoon, I will give a talk at the seminar of the R&D department at Desjardins General Insurance, on correlation in claims reserving. A lot of interesting papers have been published recently on that topic. On multivariate Chain Ladder, some interesting articles have been published, e.g. the one by Carsten Prohl and Klaus Schmidt (here) or the one by Michael Merz and Mario Wuthrich (there).

But I think another interesting perspective (so far, not in claims reserving, but one should find some time to look at it) should be about multivariate regression (multivariate GLM’s), e.g.

All that will be mentioned in the talk. Slides can be downloaded here,

The dataset used in the example can be obtained with the code below

> P.corp=read.table("http://freakonometrics.blog.free.fr/public/data/auto-corporel.csv", +        header=FALSE,sep=";",na.strings = "NA",dec=",") > P.corp=as.matrix(P.corp) > n=nrow(P.corp) > P.mat =read.table("http://freakonometrics.blog.free.fr/public/data/auto-materiel.csv", +        header=FALSE,sep=";",na.strings = "NA",dec=",") > P.mat=as.matrix(P.mat) > P.mat=P.mat[1:n,1:n] >  P.mat = P.mat[2:10,1:9] >  P.corp= P.corp[2:10,1:9] > n=9 > P.tot = P.mat + P.corp

Mais d’où sort la loi de Pareto ?

Un billet pour Jean François qui me demandait d’où sortait vraiment la loi de Pareto… Je ne reviendrais pas sur les applications de la loi de Pareto ici (je peux renvoyer à un ancien billet sur les applications en réassurance ici par exemple), si ce n’est pour noter que l’on apprend des choses surprenantes sur internet. Par exemple dans l’article sur la théorie du portefeuille de Wikipedia (ici) on apprend que dans le modèle moyenne-variance, les rendements sont supposés suivre une loi de Pareto…

Mais passons. La question de Jean François était sur l’origine de la loi de Pareto, et sur son lien avec la règle du 80-20. Pour les amateurs d’italien, on peut trouver les premières traces de la loi dans

(partiellement en ligne ici). Malheureusement, mes compétences en italien sont (très) limitées. Mais pour faire simple Vilfredo Pareto propose une loi en fonction puissance pour modéliser les revenus, en 1895,

Il fait également de l’estimation de son indice de queue (en calculant la pente dans le Pareto-plot comme on dirait maintenant)

Pour la version en français, on peut trouver la première apparition de la loi dans

(que l’on peut lire ici). Là encore, une illustration sur les revenus est proposée,

Le point de départ est une relation affine entre le logarithme de la fonction de répartition (ou plutôt de survie) et le logarithme du revenu,
et là encore une application est proposée (il me semble que la différence d’un point – environ – entre la version en italien et la version en français vient du fait qu’il travaille ici sur les fonctions de répartition, alors qu’auparavant, il travaillait davantage sur les densités),

Il propose alors une interprétation en terme de quantile (ce qui a été traduit par la loi du 80-20),

Alors il semble toutefois que Pareto n’est jamais explicitement noté que 20% des gens possédaient 80% de la richesse (qui a donné naissance au principe 80-20), ce qui pourrait sous-entendre que Pareto aurait devancé les travaux de Lorenz.Ce nom 80-20 semble avoir été attribué par Joseph Juran en 1954, dans un article “Pareto, Lorenz, Cournot, Bernoulli and others“. Il s’en est repenti (et excusé) par la suite, dans un mea culpa (ici).
Bref, tout cela a des implications évidentes sur les politiques de redistribution de la richesse, et donc sur la façon dont doit être mis en place un système d’impôts sur le revenus… Mais c’est une autre histoire…

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

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

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.

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.

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…

Encore une année record…

Comme l’a noté la blogosphère ces dernières heures, il semble que des records aient (encore) été battus en terme de catastrophes naturelles: “catastrophes naturelles : une année 2010 exceptionnelle” (ici) ou “2010, année record des catastrophes naturelles” (). Maintenant on peut remonter un peu dans le temps aussi: “catastrophes naturelles: 2008 bat des records” (ici) ou “2008 a battu des records de catastrophes naturelles” (), “en 2007, les catastrophes naturelles ont atteint un nombre record” (ici), “catastrophes naturelles : 2005 bat des records en matière d’intensité et de coûts” (ici) ou “les catastrophes ont entraîné une facture record en 2005” (), “déjà record, le bilan 2004 des catastrophes naturelles s’alourdit” (ici), “année 2000 : record de catastrophes naturelles” (ici) et “catastrophes naturelles : encore une année record” (en 2000, ) etc, etc.

Bon, en tant que statisticien, les rares résultats que je connaisse sur les records sont des théorèmes en log (i.e. en variation lente). Si http://freakonometrics.hypotheses.org/files/2015/12/record01.png est une suite de variables i.i.d. et que l’on définie les durées entre records,

http://freakonometrics.hypotheses.org/files/2015/12/records02.png

http://freakonometrics.hypotheses.org/files/2015/12/record03.png, avec http://freakonometrics.hypotheses.org/files/2015/12/record04.png. Alfred Rényi a montré que

http://freakonometrics.hypotheses.org/files/2015/12/record08.png

alors que Marcel Neuts a montré que

http://freakonometrics.hypotheses.org/files/2015/12/record09.png

En particulier, on note que

http://freakonometrics.hypotheses.org/files/2015/12/record11.png

Autrement dit, plus le temps avance, plus les records doivent devenir rares.
Par exemple, regardons la simulation suivante. Considérons un processus Poissonnien homogène, http://freakonometrics.hypotheses.org/files/2015/12/catnat01.png, avec les dates des catastrophes, et http://freakonometrics.hypotheses.org/files/2015/12/catnat02.png leur coût.

On note

http://freakonometrics.hypotheses.org/files/2015/12/catnat04.png

le nombre annuel de catastrophes (éventuellement, on ne considère que les catastrophes dont le cout dépasse un seuil prédéfini, comme cela semble généralement être le cas).

On peut définir le processus de records,

http://freakonometrics.hypotheses.org/files/2015/12/catnat07.png

qui ressemble à ça

Bon, maintenant, si on a des records tous les ans (ou presque) c’est probablement que l’hypothèse i.i.d. n’était pas valide.
Supposons par exemple qu’il y ait de l’inflation, i.e. les coût sont http://freakonometrics.hypotheses.org/files/2015/12/catnat03.png.

Dans ce cas, le processus de nombre d’événements annuels est

et celui des records

Si on répète l’expérience un grand nombre de fois, le nombre moyen d’événements par année me cesse de croître (en rouge ci-dessous)

et le processus de records moyen (avec le cas i.i.d. en bleu),

Autrement dit, si l’inflation est trop importante, il est possible d’avoir régulièrement de nouveaux records.
Une autre possibilité est de supposer qu’il n’y a pas d’inflation, mais que le nombre de sinistres ne cesse de croître (ici un processus de Poisson à intensité linéaire), c’est à dire une fréquence plus importante des événements climatiques extrêmes

Le nombre annuel de catastrophes est alors

soit, si l’on considère le processus de records

Là encore, si l’on répète l’expérience un grand nombre de fois, le nombre moyen de catastrophes ne cesse de croître

alors que le processus de records a la moyenne suivante

Bref, il est possible d’avoir effectivement des records presque tous les ans. Cela signifie probablement que les sinistres coûtent de plus en plus cher (hausse de la valeur à risque), mais aussi probablement que les catastrophes sont de plus en plus fréquentes. On retrouve finalement ce que je disais il y a quelques mois maintenant (ici).

Is it that stupid to make extremely long term forecast when studying mortality ?

I received recently a comment by FCA (here) who raised an important question, about forecast in dynamic mortality models. (S)he mentioned that from his(her) point of view, the econometric models I considered were “good to predict for the next, say, 3 or 4 years. Not for the next 50 years…”. Which was the message I tried to stress last year in a conference about retirement in France (here). But from a quantitativepoint of view, how inconsistent were forecasts made 35 years ago, or 60 years ago ?

Consider here the Lee Carter model, obtained on the periods 1816-1950 (in black below), 1816-1975 (in red) and 1816-2000 (in blue), unfortunately, it is difficult to compare http://freakonometrics.blog.free.fr/public/maths/viekt.png‘s since we have identifiability problems here. Nevertheless, we if consider affine transformation so that  http://freakonometrics.blog.free.fr/public/maths/viekt.png‘s are equal in 1900 and 1950 (say), we obtain

On that graph, we considered an ETS (AAN) forecast. If we do not consider the entire series for forecasting, but only observations following WWI (1945), we obtain

For sketches of the R code,

T=1980
base0=data.frame(D,E,A,Y,a=as.factor(A),
y=as.factor(Y))
base=base0[base0$Y<=T,]
LC2=gnm(D~a+Mult(a,y),offset=log(E),family=
poisson,data=base)
A=LC2$coefficients[1]+LC2$coefficients[2:110]
B=LC2$coefficients[111:220]
K0=LC2$coefficients[221:length(LC2$coefficients)]
Y=as.numeric(K0)
K1=c(K0,forecast(ets(Y,model="AAN"),h=240)$mean)
K2=c(K0,forecast(auto.arima(Y,allowdrift=TRUE),h=240)$mean)
MU=matrix(NA,length(A),length(K1))
MU1=MU2=MU
for(i in 1:length(A)){
for(j in 1:length(K1)){
MU1[i,j]=exp(A[i]+B[i]*K1[j])
MU2[i,j]=exp(A[i]+B[i]*K2[j])
}}
x=40
s=seq(0,109-x-1)
t=2000
Pxt1=cumprod(exp(-diag(MU1[x+1+s,t+s-base1$Year[1]-1])))
Pxt2=cumprod(exp(-diag(MU2[x+1+s,t+s-base1$Year[1]-1])))
r=.035
m=70
h=seq(0,39)
V1=1/(1+r)^(m-x+h)*Pxt1[m-x+h]
V2=1/(1+r)^(m-x+h)*Pxt2[m-x+h]
M=cbind(V1,V2)
apply(M,2,sum)

Actually, it is not that bad…. even if it is only a qualitative intuition. Again, I am not a demographer, and my interest is more on actuarial science… so if we look at the estimation of annuities (still the same insurance contract, as here) for some insured of age 40 in 2000, we get the following graph (where forecasts http://freakonometrics.blog.free.fr/public/maths/viekt.png‘s were obtained on the complete series, i.e. from 1816 until the year we consider),

(here it means that in 1900, I had to forecast mortality for someone of age 40 in 2000… so we had to forecast mortality with a 150 year horizon). Obviously, even if we are able to forecast improvement of mortality rates, it is not enough since it looks like, each year, improvement are alway higher than what what expected. Note that if we run it twice (since there might be problem with initial values in the econometric procedure) we obtain something similar,

So, the output is consistent. And if we change the way we predict future values, e.g. on focusing only on the past 50 years, i.e.

K1=c(K0,forecast(ets(Y[(length(Y)-50):length(Y)],
model="AAN"),h=240)$mean)
K2=c(K0,forecast(auto.arima(Y[(length(Y)-50):length(Y)],
allowdrift=TRUE),h=240)$mean)

we obtain the following graph for the annuity associated to an insurance contract sold in 2000,

so that relative changes compared with 1980 are (in %)

Hence, over a bit more than 25 years, we underestimated annuities of 25%. We if start to take into account possible investments, it is not so bad, I think….  don’t you think ?

 

Finding roots of functions in actuarial science

The following simple code can be used to find roots of functions (based on the secant algorithm),

secant=function(fun, x0, x1, tolerence=1e-07, niter=500){
for ( i in 1:niter ) {
	x2 <- x1-fun(x1)*(x1-x0)/(fun(x1)-fun(x0))
	if (abs(fun(x2)) < tolerence)
		return(x2)
	x0 <- x1
	x1 <- x2
}}

It can be interesting in actuarial science, e.g. to find the actuarial rate so that to present values are equal. For instance, consider the following capital, given only if the insured is still alive (this example was initially considered here). We would like to find the rate so that the probable discounted value is 600,

> Lx=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/TV8890.csv",
+ header=TRUE,sep=";")
> capital=c(100,100,125,125,150,150)
> n=length(capital)
> x=0.035
> X=45
> f=function(x){
+ capital.act=capital*(1/(1+x))^(1:n)
+ PROBA=Lx[((Lx[,1]>X)*(Lx[,1]<=(X+n)))==1,2]/Lx[(Lx[,1]==X)==1,2]
+ return(sum(capital.act*PROBA))}
>
> f1=function(x){f(x)-600}
> secant(f1,0,0.1)
[1] 0.06022313
> f(0.06022313)
[1] 600

*