All posts by Arthur Charpentier

Arthur Charpentier, professor in Montréal, in Actuarial Science. Former professor-assistant at ENSAE Paristech, associate professor at Ecole Polytechnique and assistant professor in Economics at Université de Rennes 1.  Graduated from ENSAE, Master in Mathematical Economics (Paris Dauphine), PhD in Mathematics (KU Leuven), and Fellow of the French Institute of Actuaries.

Econometrics: Learning from Statistical Learning Techniques

In two weeks, I will be invited as a keynote speaker in London, to give a talk on what can central bank policymakers learn from other disciplines. Which is an interesting question. Initially, I wanted to give a talk on actuarial science, large risks, and connexion with finance, but I will finally give a talk on connexions between Econometrics and Machine Learning, and how we can – as Econometricans – actually learn a lot from people coming from the Statistical Learning community – in the spirit of Varian (2013). I will upload the slides within the next ten days….

Et si je tweetais en anglais pour une mauvaise raison

Sur Twitter, assez naturellement, j’ai tendance à tweeter en anglais. J’ai toujours eu l’impression qu’en anglais, la limite des 140 caractères était moins contraignante. Qu’en anglais, la même chose pouvait être dite de manière plus concise.

Histoire de vérifier cette vague intuition, je suis allé chercher des livres sur Gutenberg. Le plus simple est de récupérer un livre en français, au format texte (txt), de compter le nombre de symboles, et de trouver le même livre en anglais. La fonction, à partir d’un code de livre, ressemble à ça

> cpte=function(no){
+     url0=paste(
+         "http://www.gutenberg.org/cache/epub/",
+         no,"/pg",no,".txt",sep="")
+     download.file(url0,destfile="u.txt")
+     s=scan("u.txt",what="character")
+     sum(nchar(s))
+ }

Par exemple, pour David Copperfield, de Dickens, le  code pour la version originale, en anglais, est

> gb=766

et en français, il faut récupérer 2 fichiers

> fr_a=17868
> fr_b=17869

Aussi,

> g=cpte(gb)
essai de l'URL 'http://www.gutenberg.org/cache/epub/766/pg766.txt'
Content type 'text/plain; charset=utf-8' length 1992535 bytes (1.9 MB)
=================================================
downloaded 1.9 MB

Read 288181 items
> f=cpte(fr_a)+cpte(fr_b)
essai de l'URL 'http://www.gutenberg.org/cache/epub/17868/pg17868.txt'
Content type 'text/plain; charset=utf-8' length 1161492 bytes (1.1 MB)
=================================================
downloaded 1.1 MB

Read 191311 items
essai de l'URL 'http://www.gutenberg.org/cache/epub/17869/pg17869.txt'
Content type 'text/plain; charset=utf-8' length 1135661 bytes (1.1 MB)
=================================================
downloaded 1.1 MB

Read 186515 items

le ratio du nombre de caractère en anglais, et en français, est de l’ordre de 90%. Autrement dit, la version anglaise contient 10% de caractères en moins que la version française,

> g/f
[1] 0.9090498

Si on regarde une pièce de théâtre, comme Othello, la version anglaise est 25% plus courte que la version française.

> fr=18179
> gb=1793
> g=cpte(gb)
essai de l'URL 'http://www.gutenberg.org/cache/epub/1793/pg1793.txt'
Content type 'text/plain; charset=utf-8' length 188634 bytes (184 KB)
=================================================
downloaded 184 KB

Read 24951 items
> f=cpte(fr)
essai de l'URL 'http://www.gutenberg.org/cache/epub/18179/pg18179.txt'
Content type 'text/plain; charset=utf-8' length 239158 bytes (233 KB)
=================================================
downloaded 233 KB

Read 36157 items
> g/f
[1] 0.746805

Plus généralement pour les pièces de théâtre de Shakespeare, on valide mon intuition : les versions en anglais sont toujours plus courtes que les versions en français,

Pour faire ce graphique, on peut utiliser la trentaine de livres que j’ai regardé, à l’aide du fichier

> download.file(
+ "http://freakonometrics.free.fr/livre.RData",
+ "livre.RData")
essai de l'URL 'http://freakonometrics.free.fr/livre.RData'
Content type 'text/plain' length 1095 bytes
=================================================
downloaded 1095 bytes

> load("livre.RData")

La base ressemble à ça

> base
           nom        FR      UK language   ratio
1  3 Mousquet.    1110255 1156011       FR 104.12
2  Dorian Grey     398653  404799       UK 101.54
3  Pierre-Jean     248299  228838       FR  92.16
4      Odyssée     691648  651044       GR 118.17
5       Iliade     550924  514953       GR  74.45
6  Jules César     157726  105685       UK  67.00
7      Macbeth     160571   96276       UK  59.95
8      Othello     189984  141881       UK  74.68
9       Hamlet     275629  145380       UK  52.74
10   Com. Mép.     115583   97704       UK  84.53
11  Ant.-Cléo.     185806  130857       UK  70.42
12   Rom.-Jul.     195273  126323       UK  64.69
13    Roi Lear     196185  124483       UK  63.45
14  Sganarelle      50408   63387       FR 125.74
15 OliverTwist     778914  770784       UK  98.95
16  D.Copperf.    1808770 1644262       UK  90.90
17    Assomoir     585237  565444       FR  96.61
18   T. Raquin     351525  353837       FR 100.65
19    20000 L.     745559  518932       FR  69.60
20    80 Jours     367388  329073       FR  89.57
21  5 Semaines     491510  483289       FR  98.32
22 CentreTerre     362617  370584       FR 102.19
23          93     585073  637499       FR 108.96
24  Notre Dame     885048  928256       FR 104.88
25   Jane Eyre     922233  919521       UK  99.70
26    Possédés    1241293 1293319       RU 104.19
27 Guerre-Paix    2409360 2739725       RU 113.71
28     Corbeau      27884   21886       UK  78.48
29  Mme Bovary     585237  565444       FR  96.61

Si on visualise tous les titres, et pas seulement les pièces de Shakespeare, on obtient

Autrement dit, la différence n’est pas aussi significative qu’elle y paraissait. Plus troublant, si on regarde les livres écrit en anglais (en rouge ci-dessous) on observe effectivement que la version française est plus longue. Par contre, pour les livres écrits en français (en bleu), la version anglaise est sensiblement de la même taille…

A suivre donc…

Fibonacci, les lapins, le nombre d’or et les calculs actuariels

Léonard de Pise, souvent plus connu sous le nom de Fibonacci, a laissé sa trace dans les livres de mathématiques pour avoir étudié, ou plutôt modélisé, la reproduction des lapins. De la dynamique du nombre de lapins est né le nombre d’or, le fameux ratio magique, que l’on retrouve dans toute la peinture et l’architecte classique. Mais la contribution la plus importante est probablement le Liber Abaci, datant de 1202, dans lequel Fibonacci a introduit explicitement le principe d’actualisation financière en Europe. Fibonacci a grandi dans la colonie pisane de Bougie, en Algérie, où son père était officier des douanes. Et c’est probablement grâce à la culture arabe que les premiers calculs actuariels ont été formalisés.

La dynamique de reproduction des lapins et le nombre d’or

Fibonacci est connu par tous les étudiants du monde pour avoir posé le problème suivant:

partant d’un couple, combien de couples de lapins obtiendrons-nous après un nombre donné de mois sachant que chaque couple produit chaque mois un nouveau couple, lequel ne devient productif qu’après deux mois ?

Formellement, on cherche à décrire le comportement de la suite définie par la relation http://latex.codecogs.com/gif.latex?u_n=u_{n-1}+u_{n-2} (avec deux valeurs initiales). Comme cette suite est linéairement récurrente d’ordre deux, on peut écrire son équation caractéristique. On obtient une équation du second degré, http://latex.codecogs.com/gif.latex?x^2-x-1, dont une des racines est connue sous le nom de « nombre d’or »,

http://latex.codecogs.com/gif.latex?\varphi=\frac{1+\sqrt{5}}{2}

qui vaut environ 1,618. On peut alors montrer que l’on obtient la formule dite de Binet,

http://latex.codecogs.com/gif.latex?u_n=\frac{1}{\sqrt{5}}\left(\varphi^n-\left(\frac{-1}{\varphi}\right)^n\right)

Comme l’avait déjà remarqué Johannes Kepler, vers 1611, le taux de croissance des nombres de Fibonacci, c’est-à-dire , converge vers le nombre d’or, http://latex.codecogs.com/gif.latex?\varphi. Ce nombre est appelé en anglais « golden ratio », car il s’agit aussi d’une proportion, d’un ratio, défini comme l’unique rapport http://latex.codecogs.com/gif.latex?a/b entre deux longueurs http://latex.codecogs.com/gif.latex?a et http://latex.codecogs.com/gif.latex?b, telles que la somme http://latex.codecogs.com/gif.latex?a+b sur la plus grande des deux (http://latex.codecogs.com/gif.latex?a sur la Figure 1) soit identique, autrement dit

http://latex.codecogs.com/gif.latex?\frac{a+b}{a}=\frac{a}{b}

Alors forcément, la plus grande longueur est 61,8% plus grande que la petite. Euclide évoquait ce ratio très particulier dans ses Eléments. Comme le montre la Figure 1, ce ratio est utilisé pour construire une spirale. Luca Pacioli, un moine franciscain italien du Moyen Age, l’appelle la “divine proportion“. En effet, comme l’avait énoncé Saint Thomas d’Aquin, « les choses qui sont dotées de proportions correctes réjouissent les sens », ce qui semble conférer un caractère quasi-divin à ce nombre.

Figure 1 : interprétations géométriques du nombre d’or

Comme le montre la Figure 1, à gauche, à partir de rectangles dont le rapport des longueurs des deux côtés est égal au nombre d’or, on peut construire itérativement une succession de carrés emboîtés, sous forme de spirale. On retrouve ce découpage dans la composition d’un grand nombre de tableaux de la Renaissance, mais aussi dans plusieurs ouvrages architecturaux. Et si l’on colle deux rectangles dont le ratio des longueurs des côtés est égal au nombre d’or, un verticalement et un horizontalement, le prolongement de la diagonale du rectangle horizontal passera par le sommet du rectangle verticale. Cette figure se retrouve dans de nombreux croquis de Le Corbusier.

La notion de récursivité

Au-delà du fait que c’est à Fibonacci que l’on doit le fameux nombre d’or, cette suite récursive est surtout très populaire en informatique pour introduire la notion de récursivité. En effet, un algorithme pour calculer les termes de cette suite est

F = function(n) {

if(n>1) Fn = F(n-1) + F(n-2)

if(n<=1) Fn = n

return(Fn) }

Cet algorithme est simple, mais le temps d’exécution est long, très long. Si certains notent qu’il croit exponentiellement vite, on peut montrer assez simplement que le temps d’exécution est donné par une relation de type de celle de Fibonacci.

La comptabilité, en chiffres romains

Mais souvenons-nous que Fibonacci est né à la fin du XIIème siècle. Si la géométrie n’a pas beaucoup changé, l’arithmétique, elle a beaucoup changé. Pour rappel, on comptait alors encore bien souvent à l’aide de la numérotation romaine. Cette dernière est basée sur 7 symboles, I (pour 1), V (pour 5), X (10), L (50), C (100), D (pour 500) et enfin M (1000). A partir de ces symboles, certains nombres se construisent par addition (par exemple 11 va s’écrire XI) et d’autres par soustraction (4 va s’écrire IV): un symbole plus petit placé avant correspond à une soustraction, et un symbole plus petit placé après, à une addition. On s’attend ainsi à ce que tous les symboles soient placés par ordre décroissant (et si c’est le cas, on les ajoute), mais s’il y en a un qui est placé avant, on va le soustraire. Aussi, DCXIV contient les symboles D (500), C (100), X (10), I (1) et V (5). Comme 1 est avant 5, c’est qu’il faut lire -1, et pas +1. Et donc DCXIV vaut 500+100+10-1+5 soit 614. De même XLV correspond à -10+50+5, soit 45. Il existait aussi une forme d’écriture multiplicative. Par exemple 300 pouvait s’écrire CCC ou IIIC. Et 80 IVXX (d’où une explication possible de la forme française de “quatre-vingt”).

L’origine de ces symboles se trouve dans les entailles pratiques par les étrusques, sur des bouts de bois, ou des os. Mais chez les étrusques, seule l’addition de symboles était autorisée. Par la suite, les égyptiens ont adopté un système décimal, et introduit l’addition (uniquement), ce qui a permis aux scribes de faire des observations statistiques, et des prévisions. Lorsqu’ils adoptèrent la langue et les pratiques grecques, vers 300 avant Jésus Christ, les égyptiens ont commencé à tenir de véritables documents comptables, avec des stocks de données agricoles, et des flux monétaires. A la même époque en Grèce, les lettres de change existaient, et le prêt à intérêt était fréquent. A Athènes, la cour des logistes avait le rôle d’expert-comptable. Les romains ont poursuivi et amélioré cette tradition, avec des livres de recettes et de dépenses. L’arithmétique était complexe, comme en témoigne Schärlig (2014), et pourtant, ce système d’encoches a été utilisé probablement jusqu’au XVIIème siècle, en Europe.

Le livre des calculs, Liber Abaci, et la règle de trois

Fibonacci est l’auteur d’un Liber Abaci, correspondant à un livre des calculs (Abaci donnera par la suite le terme abaque). En fait, alors que la norme à l’époque était d’utiliser l’arithmétique romaine, Leonardo s’est énormément inspiré des mathématiques arabes. Les nombres arabes présentaient l’avantage d’utiliser le système décimal, ce qui a fait gagner énormément de temps de calcul (on peut aussi remonter à l’arithmétique indienne, basé sur le même système). Il y avait ainsi neuf symboles correspondant à 1, 2, 3, 4, 5, 6, 7, 8, 9 (et le 0 a été introduit au Vème siècle). Afin d’illustrer la simplicité des calculs, Fibonacci évoque intensément l’utilisation de la « règle de trois » (même s’il ne l’appelle pas ainsi, il faudra attendre Le Gendre (1684) pour voir évoqué cette terminologie). Cette règle, relative à la théorie des proportions, se retrouve aussi beaucoup utilisée chez Euclide (au IIIème siècle avant Jésus Christ): « si quatre nombres sont en proportions, le produit des extrêmes [le plus grand et le plus petit] est égal au produit des moyens ». En fait, la règle de trois permet de trouver ce qui s’appelle parfois la « quatrième proportionnelle », qui est le quatrième nombre à mettre dans un tableau de proportionnalité dont 3 cases sur 4 sont déjà remplies. Certains ont appelé la règle de trois la « règle d’or ».

En actuariat, le calcul des provisions pour sinistres à payer par conservation des cadences de paiements, repose historiquement sur cette règle. Le Tableau 1 illustre ainsi l’utilisation de la règle de trois pour estimer les paiements futurs : fin 2014, 57 millions d’euros avaient été payé au titre des sinistres de 2014, et fin 2015, 83 millions ; si fin 2015, 63 millions avaient été payés au titre des sinistres de 2015, on peut estimer qu’un peu plus de 90 millions seront payés fin 2016. Ce montant http://latex.codecogs.com/gif.latex?\color{red}{x} est obtenu, dans la méthode dite Chain Ladder (décrite dans Charpentier et al. (2010)), par la règle de trois.

paiements

année 1

paiements

année 2

2014

http://latex.codecogs.com/gif.latex?57

http://latex.codecogs.com/gif.latex?83

2015

http://latex.codecogs.com/gif.latex?63

http://latex.codecogs.com/gif.latex?{\color{red}{x}}=\frac{63\times%2083}{57}

Table1 : illustration de la règle de trois en provisionnement

On retrouve de nombreuses utilisations de cette règle dans les calculs comptables (on pourra consulter Taulane (1689)), mais c’est surtout à Fibonacci qu’on doit la formalisation de la règle, dans un contexte financier.

Les bases de l’actualisation

Mais on doit surtout à Fibonacci d’avoir formalisé le calcul actuariel. Rappelons que la valeur actuelle d’un ensemble de flux financiershttp://latex.codecogs.com/gif.latex?\lbrace%20F_1,\cdots,F_n\rbrace, versés dans le futur, à différentes dates s’écrit

http://latex.codecogs.com/gif.latex?V=\sum_{t=1}^n\frac{F_t}{(1+i)^t}

où http://latex.codecogs.com/gif.latex?i est le taux d’actualisation. Cette à cette formule que va aboutir Fibonacci, dans le Liber Abaci, lorsqu’il va tenter de résoudre le problème suivant

« Un homme se rend à Lucca pour affaires, où il double son argent. Il dépense là-bas 12 dinars. Ensuite, il va à Florence où il peut là encore doubler son argent, et il dépense encore 12 dinars. Puis il va ensuite à Pise, double encore son investissement, et dépense à nouveau 12 dinars. Quelle somme avait-il au début? »

Ce calcul de la somme initiale correspond fondamentalement au calcul d’une valeur actuelle. Pour reprendre les notations de Fibonacci pour investir 1 à Pise (la troisième ville), il faut disposer de ½½½, qu’il interprète comme un facteur d’escompte. Pour investir 1 à Florence (la seconde ville), il fallait disposer de ½½ au début. Etc. Pour connaître la somme initiale, on notera que le facteur d’actualisation est ½+½½+½½½, soit 7/8. Le capital initial est alors 7/8×12, soit 10 ½ . Avec des notations plus contemporaines, avec http://latex.codecogs.com/gif.latex?i = 100%, on veut calculer,

http://latex.codecogs.com/gif.latex?\frac{12}{(1+i)}+\frac{12}{(1+i)^2}+\frac{12}{(1+i)^3}

On retrouve ici la définition classique de la valeur présente d’un ensemble de flux financier. Mais plus sophistiqué, Fibonacci propos ensuite le problème suivant

« Un soldat se voit offrir une annuité par un roi, de 300 par an, payée à une cadence trimestrielle, à hauteur de 75 par trimestre. Soudain, le roi change pour 300, payé en fin d’année. Le soldat est capable de gagner 2 pour 100 placés, par mois (à la fin de chaque trimestre) à partir de son investissement »

On retrouve ici un problème qui pourrait être posé à un examen sur les calculs actuariels.

Conclusion

Depuis Goetzmann (2005) et la redécouverte des travaux de Léonard de Pise, beaucoup considèrent que Fibonacci a introduit le calcul actuariel en Europe, ou plutôt son formalise, avec des notations très contemporaines. Mais il convient de rappeler que l’actualisation existait depuis bien longtemps. Des traces de formalisation se retrouvent dans toutes les cultures, à commencer probablement en Chine, où Les Neuf Chapitres sur l’Art des Mathématiques, datant de 263, évoquent un problème proche de celui présenté par Fibonacci (avec là encore un voyageur qui part dépenser son argent). On retrouvera ensuite ce problème (et sa résolution) dans des textes arméniens, dans des papyrus égyptiens, mais aussi sur le Lilivati de Bhāskarācārya, écrit en Inde vers 1150, juste avant la formalisation par Léonard de Pise. Mais la lecture de Fibonacci est facilité par l’utilisation de chiffres arabes, qui rendent les calculs plus proches de ceux que nous avons l’habitude de voir, tous les jours.

Références

Charpentier, A., Devineau, L. & Nessi, J.M. (2010) Mesurer le risque lors du calcul des provisions pour sinistres à payer. Risques 83, 93-100.

Fibonacci (1202). Fibonacci’s Liber abaci, traduit par Laurence E. Sigler A translation into modern English of Leonardo Pisano’s Book of calculation, Springer-Verlag,‎ (2002).

Goetzmann, W. N. (2005) Fibonacci and the Financial Revolution. In Goetzmann and Rouwenhorst, The Origins of Value: The financial innovations that created modern capital markets, Oxford University Press.

Le Gendre, F. (1684). L’arithmétique en sa perfection.

Le Goff, J.P. (2000). Une mise en perspective historique de la quatrième proportionnelle à la fonction linéaire, in Cahiers de la Maison de la Recherche en Sciences Humaines, n° 24, L’épistémologie de la règle de trois. Caen, Presses Universitaires de Caen, déc. 2000.

Schärlig, A. (2014). Les moyens de calcul du comptable antique. Revue d’Histoire des Comptabilités, 6.

Taulane, J. (1689) Trésor d’arithmétique des marchands et des géomètres.

Estimer une Proportion sur une Question Sensible, par Enquête

Comme le rapportait Slate hier (dans un article qui, malheureusement ne mérite pas d’être lu), un joli sujet a été proposé au baccalauréat au lycée français de Pondichéry, en mathématiques.

Ce sujet est d’autant plus intéressant qu’il revient sur une méthode assez classique pour questionner les gens sur des questions sensibles (ici sur le téléchargement pirate, mais on peut poser la même question sur la fraude, ou des pratiques sexuelles).

Comme le notait , cette technique avait été introduite initialement par Stanley Warner, voilà un demi-siècle,

Cette technique est particulièrement intéressante, tant d’un point de vu intellectuel que pratique. Plusieurs sur Twitter me faisaient remarquer que, contrairement à ce qui était raconté dans Slate, cette technique est effectivement utilisée.

Quand j’enseignais cette technique, je prenais la version simple, avec une pièce de monnaie. Supposons qu’une proportion  de personnes fraudent le fisc. Et qu’une proportion  aime Prince, va à la piscine une fois par mois, a un diplôme de l’enseignement supérieur, etc. Peut importe, la seconde proportion doit pouvoir être estimée sans problème avec une autre enquête (voire par le recensement). Bref, on a un estimateur de . Ici, les personnes interrogées tirent une pièce et répondent soit à la question “avez-vous déjà fraudé au Fisc?” s’ils tombent sur face, et à la question “aimez-vous Prince?” s’ils tombent sur pile. On pourrait dire qu’ils tirent une pièce dans leur tête, mais il faut s’assurer que la probabilité de tomber sur face est effectivement 1/2.

La magie est de noter que la proportion de personnes qui fraudent le fisc peut s’écrire

On peut avoir, par une autre enquête un estimateur du second terme, . Et en regardant des personnes qui ont répondu “oui” à une question inconnue, on a un estimateur pour

Et par différence, on en déduit un estimateur de . Les maths sont dans l’article de Stanley Warner, mais il faut reconnaître que la méthode est vraiment jolie…. Dans le cas du sujet du bac, la mise en oeuvre est un peu différente, mais l’idée est un peu la même….

(si certains ont des exemples pratiques d’application, les commentaires sont ouverts, ou sinon par mail pour joindre des documents en pièce jointe)

Partager la donnée : une question de confiance ?

Demain, les Rendez-vous Parlementaires de Bretagne sur l’Economie Numérique vont se poursuivre l’après midi, à partir de 14H dans les locaux de la French Tech Rennes-Saint Malo. Un groupe de travail autour de Partager la donnée : une question de confiance ? sera organisé, animé par Simon, avec une trentaine de participants venus du monde de l’entreprise mais aussi de la recherche, du milieu associatif et bien sur des élus de Rennes Métropole. A suivre donc….

Questions autour du mécanisme Cat Nat, en France

Il y a quelques semaines, Thierry Gouby de la Tribune de l’Assurance me contactait pour me poser quelques questions sur le mécanisme cat nat. J’ai commencé à regrouper mes idées, et comme j’avais un texte 6 fois trop long, par rapport à ce qui m’était demandé, je me suis dit que je pourrais mettre en ligne les éléments qui ne resteront pas, dans la restitution finale…

Est-on en mesure aujourd’hui en France d’évaluer techniquement et avec précision les zones inondables ou les zones de sécheresse ?

Les inondations désignent un ensemble assez vaste de périls, allant d’un ruissellement à un débordement d’un cours d’eau, à une remontée de nappe phréatiques , voire une rupture de barrage (si cette dernière est causée par un phénomène naturel). Depuis quelques années, l’évaluation préliminaire des risques d’inondation (EPRI) a permis d’identifier 122 territoires à risques d’inondation important (TRI). Une cartographie fine des risques a été effectuée sur chacun des territoires, afin de prendre des actions – localement – permettant d’anticiper et de réduire l’impact des inondations, dans l’esprit de la stratégie nationale de gestion des risques d’inondation (SNGRI). Entre 1998 et 2013, près d’un million et demi de sinistres ont été indemnisés au titre du risque inondation, pour un coût total de 17 milliards d’euros. Comme le montre la carte ci-dessous, les TRI sont répartis de la manière suivante

Le terme “sécheresse” désigne un risque un peu particulier, car il concerne le phénomène de subsidence. Dans le cadre du régime dit cat nat, il s’agit d’un risque de tassements différentiels des argiles (les sols argileux peuvent gonfler ou se tasser selon leur teneur en eau, et ces variations peuvent engendrer des dégâts importants aux constructions). On parlera de sécheresse géotechnique. Il ne fait généralement pas de victimes. Entre 1998 et 2013, 600 mille sinistres ont été indemnisé, pour un coût proche de 8 milliards d’euros. Mais la sécheresse peut être associée à des épisodes de canicules, qui en plus d’engendrer davantage de sécheresses géotechniques (pouvant occasionner des déplacements du sol de plusieurs centimètres) peut provoquer une surmortalité, des pertes agricoles conséquentes, mais aussi des incendies de forêt. Des mesures peuvent être prises en compte de manière préventive, afin de réduire les effets des sécheresse géotechnique, comme assurer un ancrage plus profond des fondations, désolidariser autant que possible des bâtiments, ou s’assurer que les drains d’évacuation passe loin des habitations. La répartition des régions risquées peut se visualiser sur la figure suivante

A-t-on assez de données ou d’historiques chiffrés pour prédire avec justesse le risque de tempêtes par exemple ?

Enfin, le risque tempête est un risque qui ne relève pas du régime cat nat, en France métropolitaine, mais depuis une quinzaine d’années, le législateur a rendu la garantie tempête obligatoire pour l’ensemble des contrats de dommage aux biens. Le risque tempête est souvent considéré comme non spatialisé, au sens où toute la France métropolitaine peut être uniformément touché (même si les régions littorales ont un risque souvent plus important, comme la pointe Finistère). Les tempêtes de 1999 ont toutefois soulevé beaucoup de questions (69 départements ont été déclarés en état de catastrophe naturelle, mais au titre des glissements de terrains, des inondations).

La quantification du risque de type tempête (mais la méthodologie est la même pour les autres types de risques) repose sur deux éléments importants: la modélisation du phénomène naturel sous-jacent, et un modèle de vulnérabilité, permettant de passer de l’événement à une estimation du coût des dommages.

En matière des catastrophes naturelles, plusieurs types de données sont souvent utilisées : des données fines, mais sur une courte période de temps (observations journalières de niveau de fleuves par exemple), des traces historiques de catastrophes passées (rares, irrégulièrement espacées dans le temps, et probablement non exhaustives), ainsi que des catalogues d’événements associés aux scénarios de conditions climatiques pour un horizon temporel plus ou moins lointain (tels que proposés par le GIEC).

Dans le cas du risque tempête la circulation atmosphérique en Atlantique Nord est de mieux en mieux modélisée, et des données très fines peuvent être obtenues. Si la combinaison de modèle numérique de circulation combiné avec des modèles statistiques de vitesse de vent (sur un historique d’une cinquantaine d’années) permettent d’avoir des générations d’événements avec des vents moyens sur des périodes de 10 minutes, la grande difficulté reste d’avoir des modèles pour les vitesses de rafales. Cusack (2013) avait ainsi analysé les liens entre 101 années d’observations de vitesse de vent aux Pays Bas, et les dommages associés au risque tempête.

Concernant les tempêtes historiques, il a été possible d’extrapoler les trajectoires de tempêtes les plus destructrices, y compris des tempêtes qui se sont produites voilà plus d’une centaine d’années. En Mars 1876, une tempête avait par exemple traversé l’Europe, avec une trajectoire très proche de celle de Lothar, en décembre 1999,

ou celle de janvier 1884 qui n’est pas sans rappeler Daria, qui est passé au nord de l’Europe en janvier 1990,

semblent aujourd’hui bien connu par les météorologues, et les assureurs, ce qui semble effectivement confirmer l’assurabilité de ce risque, même si à l’échelle des tempêtes, la France reste un petit territoire pour assurer une parfaite mutualisation du risque. Mais beaucoup ont souligné les caractéristiques imprévisibles des conséquences des tempêtes, car, par exemple, si les tempêtes de décembre 1999 avaient eu lieu lors d’une journée de travail, le nombre de victimes aurait probablement été beaucoup plus important. On notera que depuis, les mécanismes d’alerte ont été amélioré, mais les expériences ont montré que le plus souvent, les bâtiments qui n’ont pas résisté aux tempêtes présentaient, le plus souvent, des défauts de conception ou d’entretien.

La hausse de la fréquence et de la gravité des CatNat, dans l’hexagone comme dans le monde, remet-elle aujourd’hui en cause l’assurabilité des risques climatiques ?

Les catastrophes naturelles sont définis (légalement) comme des “dommages matériels directs non assurables ayant eu pour cause déterminante l’intensité anormale d’un agent naturel”. Autrement dit, ces risques étaient, dans la loi de 1982, considéré comme “non assurables”, d’où le mécanisme un peu particulier qui permet, malgré tout, sa couverture. C’est probablement la présence de l’état, comme assureur en dernier ressort, qui permet effectivement qu’un mécanisme de type assurantiel fonctionne.

Les tempêtes, en revanche, ont été exclues du système cat nat car elles sont – a priori – non spatialisée, et peuvent donc être mutualisées.

Le changement climatique pourrait engendrer des versements, au titre des catastrophes naturelles 2 fois plus élevés qu’aujourd’hui, selon l’estimation faite par la FFSA en 2009. Les catastrophes naturelles seront normalement plus fréquentes. Une sécheresse de l’ampleur de celle de 2003 devrait avoir une période de retour d’une dizaine d’année, et non plus d’une vingtaine d’année. Mais pour la plupart des autres risques (inondation mais aussi tempêtes) la fréquence des événements majeurs ne devrait pas augmenter significativement dans les années à venir. Pour la plupart des risques climatiques, cette étude (confirmée par la CCR) a montré que le quart de la hausse était du à la aléa seul (c’est à dire la fréquence et la gravité de l’événement climatique), et trois quart à un effet de la vulnérabilité (correspond au coût assurantiel).

Un mécanisme d’indemnisation partagé entre assureurs et intervention de l’Etat est-il la bonne solution ? Les assureurs et réassureurs sont-ils finalement à l’abri avec ce système ? D’autres mécanismes sont-ils plus performants ?

Le système français est équilibré, en absence de catastrophe majeure. Mais la crainte d’une inondation majeure à Paris (du type de celle de 1910) ou d’une tempête hivernale (de la puissance de celles de 1990 ou de 1999, mal anticipée) relance régulièrement le débat sur la viabilité du mécanisme.

Le mécanisme français est assez simple, avec trois interactions principales

  1. assurés-assureurs: un contrat d’assurance (obligatoire) assez standard, à ceci près que les compagnies d’assurance ne peuvent fixer librement le montant des primes et des cotisations
  2. assureurs-réassureur: une compagnie de réassurance publique qui propose des traités de réassurance à tous les assureurs (la Caisse Centrale de Réassurance, CCR)
  3. réassureur-état: l’état intervient en tant qu’assureur en dernier ressort, en apportant sa garantie illimitée à la CCR.

Ce système hybride mêlant assurance privée et publique est basé sur une offre d’assurance privée, encadrée par les pouvoirs publics, avec un mécanisme de réassurance public, ou semi-public. On retrouve ce mécanisme général au Japon, et dans certains états américains (Californie, Floride et Texas). La France a toutefois un système plus complet car plusieurs risques sont couverts (sécheresse, mouvement de terrain, lave torrentielle, submersion marine, avalanche, séisme).

Les autres mécanismes sont assez simples,

  • la couverture est facultative, et repose intégralement sur des compagnies d’assurances privées (on peut penser à la Grande Bretagne, l’Irlande ou la Pologne)
  • un mécanisme hybride, combinant un marché concurentiel d’assurance, et des aides publiques, souvent via un fond (comme le Fond Inondation Tempête au Danemark, ou le National Flood Insurance Program aux Etats-Unis) dans le cas où l’intervention publique est programmée ex-ante, ou avec une intervention ex-post (comme en Allemagne, ou en Australie)
  • un dispositif public, obligatoire, et monopolistique d’assurance des catastrophes naturelles (comme en Espagne)
  • aucune prévention ex-ante du risque de catastrophe naturelle, mais le plus souvent, des interventions publiques (ponctuelles) ex-post (comme en Corée, en Italie, ou en Suède)

Le fait que le mode de calcul des primes soit imposé repose sur un principe de solidarité entre les assurés, quel que soit leur niveau d’exposition aux événements naturels. Plus spécifiquement, il s’agit d’un taux uniforme pour la surprime associée à la couverture du risque cat nat, par exemple 12% pour un contrat multirisque habitation et 6% pour un contrat d’assurance d’un véhicule terrestre. Afin d’inciter à la prévention, une modulation des franchises a été introduite il y a une quizaine d’années: les franchises sont plus élevées dans les communes qui ne se sont pas doté d’un plan de prévention des risques naturels (PPRN) alors que la commune a fait l’objet, par le passé, de plusieurs arrêtés cat nat. Aussi, pour une commune non dotée d’un PPRN contre l’inondation, la franchise est doublé au troisième arrêté cat nat (triplé au quatrième, etc). Ces PPRN complètent les plans locaux d’urbanisme (PLU) qui permettent de refuser (ou d’accepter sous condition) un permis de construire en zone inondable, par exemple. Notons que la prime uniforme est une caractéristique qui se retrouve, pas exemple pour le National Flood Insurance Program aux Etats-Unis. Mais contrairement à la France qui utilise les franchises pour tenir compte de plans de préventions, aux Etat-Unis, si un commune (traduction rapide de local communities) entre dans le Community Rating System, elle obtient une note tenant compte des plans de préventions mis en place, et une réduction de prime peut alors être obtenue.

Mais ces plans restent locaux, et dans le cas du risque d’inondation, ce sont souvent les actions de prévention en amont, au niveau des bassins versants, qui ont le plus d’impact sur la réduction du risque.

Mais si le mécanisme n’est pas propre à la France, c’est probablement la définition d’un sinistre relevant du mécanisme cat nat qui rend le mécanisme typiquement français. C’est aujourd’hui la publication d’un arrêté interministériel, reconnaissant l’état de catastrophe naturelle qui fait foi. Cette définition va à l’encontre du fonctionnement même de l’assurance, où les périls sont rarement “non-dénommés”. Afin de prévoir les risques (et les quantifier actuariellement) il convient qu’ils soient identifié aussi clairement que possible. La notion d’“intensité anormale d’un agent naturel” semble s’opposer à une intensité “normale”. Les inondations récurrentes n’ont alors pas vocation à être indemnisées au titre du mécanisme cat nat. Il ne s’agit pas d’agir sur la franchise, mais bien d’exclure un risque. Un autre élément propre à la France est probablement la couverture du risque sécheresse, unique en Europe. Il ne serait pas absurde de penser que ce risque ne relève pas du mécanisme cat nat, mais pourrait être associé à l’assurance décennale du constructeur du bâtiment, les premières années. Ce mécanisme pourrait être suffisamment incitatif pour que les constructeurs tiennent compte davantage des spécificités du sol.

 

La parade des nains, ou la visualisation des inégalités

Dans mon précédant billet, sur la répartition de la population sur le territoire français, j’utilisais l’indice de Gini, comme mesure d’inégalité. Et je notais que l’indice était passé de 0.45 au début du XIXème siècle à 0.75 aujourd’hui. Mais force est de constater que cet indice n’est pas très visuel (à moins que l’on ne soit à l’aise avec les courbes de Lorenz). En 1970, Jan Pen proposa introduisit la “parade des nains” afin d’illustrer les inégalités. Pour décrire des inégalités de revenus, on représente des individus dont la taille est proportionnelle au revenu.

Il y a quelques (rares) géants, et beaucoup de nains. On peut tracer la courbe associée à cette parade des nains. Formelle, pour , on trace

où  est la fonction de répartition de la variable que l’on étudie. Rappelons que pour une loi lognormale , l’indice de Gini est

De manière duale, à partir d’un indice de Gini, on peut obtenir les paramètres d’une loi lognormale sous-jacente. Par exemple, pour une loi lognormale, associée à un indice de Gini de l’ordre de 0.25, on obtient

> library(ineq)
> gini=0.25
> sigma=qnorm((1+gini)/2)*sqrt(2)
> x=rlnorm(n,0,sigma)
> Pen(x)

Pour des indices de Gini de 0.50 et de 0.75, on obtient les courbes suivantes

> gini=0.5
> sigma=qnorm((1+gini)/2)*sqrt(2)
> n=1e4
> y1=qlnorm((1:n)/(n+1),0,sigma)/
+ exp(0+sigma^2/2)

> gini=0.75
> sigma=qnorm((1+gini)/2)*sqrt(2)
> y2=qlnorm((1:n)/(n+1),0,sigma)/
+ exp(0+sigma^2/2)
> plot((1:n)/(n+1),y1,ylim=c(-.5,10))
> rect(0,0,1,1,col="light yellow",border=NA)
> polygon(c(0,(1:n)/(n+1),1),c(0,y1,0),
+ col=rgb(1,0,0,.5))
> polygon(c(0,(1:n)/(n+1),1),c(0,y2,0),
+ col=rgb(0,0,1,.5))
> abline(h=1,lty=2,col="grey")

La courbe jaune est une distribution parfaitement égalitaire. La courbe rouge est obtenue avec un indice de Gini de l’ordre de 0.50 alors que celle en bleu est obtenue avec un indice de 0.75.

Pour reprendre notre exemple de population, répartie sur le territoire, on peut aussi imaginer un territoire linéaire de 1 km de long. On suppose que la zone de gauche est la moins dense, et celle de droite est la plus dense.

  • avec une répartition uniforme, supposons que tout le monde puisse loger dans des immeubles de 5 étages.
  • avec un indice de Gini de 0.50,  les premiers 200m sont des immeubles de 1 étage, les 200m suivants, ce sont des immeubles de 2 étages, etc. A 700m, on obtient les derniers immeubles de 5 étages. Sur la fin, on a des immeubles de 15 étages sur les 50 derniers mètres, et avant un immeuble de 12 étages, et avant encore des immeubles de 10 étages.
  • dans le cas d’un indice de Gini de l’ordre de 0.75, il n’y a rien sur 25% du territoire, puis sur les 25% suivant, des immeubles de 1 étage, etc. 80% des logements ne dépassent pas 5 étages.  Par contre, sur la fin, on des tours relativement haute, dépassant 20 étages.
> n=19
> gini=0.5
> sigma=qnorm((1+gini)/2)*sqrt(2)
> y1=qlnorm((1:n)/(n+1),0,sigma)/
+ exp(0+sigma^2/2)

> gini=0.75
> sigma=qnorm((1+gini)/2)*sqrt(2)
> y2=qlnorm((1:n)/(n+1),0,sigma)/
+ exp(0+sigma^2/2)

> plot((1:n)/(n+1),y1*5,ylim=c(-.5,20))
> for(i in 1:n){
+ rect(i/(n+1)-.02,0,i/(n+1)+.02,1*5,
+ col="light yellow",border=NA)
+ rect(i/(n+1)-.02,0,i/(n+1)+.02,round(y1[i]*5),
+ col=rgb(1,0,0,.5),border=NA)
+ rect(i/(n+1)-.02,0,i/(n+1)+.02,round(y2[i]*5),
+ col=rgb(0,0,1,.5),border=NA)
+ }
> abline(h=1*5,lty=2,col="grey")
> axis(2)

ou, si on zoome, on obtient

La courbe rouge, c’est la répartition de la population, en France au XIXème siècle, ou en Allemagne aujourd’hui (tel que décrit dans mon précédant billet). La courbe bleue représente la France aujourd’hui.

La centralisation française, et la répartition de la population sur le territoire

Suite à mon billet d’hier soir sur la distribution en France, j’ai eu plusieurs commentaires – sur Twitter – qui notaient que ce n’était pas surprenant que la France soit aussi concentrée, en terme de population, compte tenu de l’importance de la centralisation en France (par opposition à l’Allemagne, par exemple).  Et comme le note wikipedia, à la page sur le centralisme, “depuis la Révolution française (et même avant), la France est un État très centraliste“.

Il y a quelques années, Mattia Bunel (aka @mattiabunel) avait rédigé un joli mémoire sur l’influence des contraintes environnementales sur la distribution de la population française entre 1793 et 1999. Au delà du travail de rédaction, il y a surtout un gros travail sur les données. Mattia a remis en forme les données de cassini.ehess.fr. Et Mattia a eu la gentille de m’envoyer ses données, avec en ligne les villages en France, leur superficie, et leur population, à plusieurs dates, entre 1793 et 1999. Pour récupérer la superficie (Mattia a passer du temps à mettre à jour, pour les villages qui ont fusionné, en particulier) le code est tout simplement

> base=read.csv("/Cassini/export.csv")
> dim(base)
[1] 41409    40
> S=base$Superficie
> n=nchar(as.character(S))
> S=substr(S,1,n-3)
> Surface=as.numeric(gsub(" ", "", as.character(S), fixed = TRUE))
> idx=which(!is.na(Surface))
> B=base[idx,]
> dim(B)
[1] 36576    40

On a ainsi les 36,000 communes françaises. Ensuite, on peut extraire, par année, la population. Histoire d’avoir un calcul de l’indice de Gini qui soit cohérent avec ce que j’avais fait hier – dans Non-Uniform Population Density in some European Countries – l’idée est la suivante : supposons qu’il y ait 2 villages, un de superficie 4, et de population 4, et un autre de superficie 1, et de population 3. Dans le billet d’hier, je raisonnais pas unité spatiale (en l’occurrence un petit carré)

et on va faire pareil ici. Autrement dit, on va répartir la population uniformément au sein du village de taille 4, ce qui va faire 4 “villages unitaires” avec 1 habitant chacun,

On a ainsi 5 unités, 4 avec une population de 1, et 1 avec une population de 3. De cet échantillon – {1,1,1,1,3} – on peut calculer l’indice de Gini. La fonction pour extraire la population, et calculer l’indice de Gini est

> library(ineq)
> LC=function(annee){
+   nom=paste("X",annee,sep="")
+   P=B[,nom]
+   P=gsub(" ", "", as.character(P), fixed = TRUE)
+   Pop=as.numeric(substr(P,3,nchar(P)))
+   P=Pop[which(!is.na(Pop))]
+   S=Surface[which(!is.na(Pop))]
+   D=P/S
+   S1=round(S/20)
+   X=rep(NA,sum(S1))
+   s=0
+   for(i in 1:length(S1)){
+     X[s+1:S1[i]]=D[i]
+     s=s+S1[i]
+   }
+   Gini(X)
+ }

Si on calcule notre indice de Gini pour toutes les années, on obtient

> Y=names(base)[7:ncol(base)]
> Y=as.numeric(substr(Y,2,5))
> gini=Vectorize(LC)(Y)
> plot(Y,gini,type="b")

On retrouve un indice de Gini qui dépasse 0.7 aujourd’hui (que nous avions obtenu hier, avec une méthodologie assez différente), mais surtout, on voit que l’indice de Gini n’a cessé d’augmenter, depuis la révolution française… L’interprétation rapide serait que le centralisme français ne cesse pas d’augmenter, depuis 200 ans.

La donnée, pierre angulaire de notre économie

Jeudi, je participerais aux Rendez-vous Parlementaires de Bretagne sur l’Economie Numérique,

où on m’a demandé de parler de données. Je dois avoir une quinzaine de minutes, mais comme toujours j’ai prévu pour 2 heures d’exposé. Les slides sont maintenant en ligne, et je me suis dis que je pourrais faire un billet pour détailler le fond de mon propos,

Continue reading La donnée, pierre angulaire de notre économie

Non-Uniform Population Density in some European Countries

A few months ago, I did mention that France was a country with strong inequalities, especially when you look at higher education, and research teams. Paris has almost 50% of the CNRS researchers, while only 3% of the population lives there.

It looks like Paris is the only city, in France. And I wanted to check that, indeed, France is a country with strong inequalities, when we look at population density.

Using data from sedac.ciesin.columbia.edu, it is possible to get population density on a small granularity level,

> rm(list=ls())
> base=read.table(
+      "/home/charpentier/glp00ag.asc",
+      skip=6)
> X=t(as.matrix(base,ncol=8640))
> X=X[,ncol(X):1]

The scales for latitudes and longitudes can be obtained from the text file,

> #ncols         8640
> #nrows         3432
> #xllcorner     -180
> #yllcorner     -58
> #cellsize      0.0416666666667

Hence, we have

> library(maps)
> world=map(database="world")
> vx=seq(-180,180,length=nrow(X)+1)
> vx=(vx[2:length(vx)]+vx[1:(length(vx)-1)])/2
> vy=seq(-58,85,length=ncol(X)+1)
> vy=(vy[2:length(vy)]+vy[1:(length(vy)-1)])/2

If we plot our density, as in a previous post, on Where People Live,

> I=seq(1,nrow(X),by=10)
> J=seq(1,ncol(X),by=10)
> image(vx[I],vy[J],log(1+X[I,J]),
+ col=rev(heat.colors(101)))
> lines(world[[1]],world[[2]])

we can see that we have a match, between the big population matrix, and polygons of countries.

Consider France, for instance. We can download the contour polygon with higher precision,

> library(rgdal)
> fra=download.file(
"http://biogeo.ucdavis.edu/data/gadm2.8/rds/FRA_adm0.rds",
+ "fr.rds")
> Fra=readRDS("fr.rds")
> n=length(Fra@polygons[[1]]@Polygons)
> L=rep(NA,n)
> for(i in 1:n) L[i]=nrow(Fra@polygons[[1]]@Polygons[[i]]@coords)
> idx=which.max(L)
> polygon_Fr=
+       Fra@polygons[[1]]@Polygons[[idx]]@coords
> min_poly=apply(polygon_Fr,2,min)
> max_poly=apply(polygon_Fr,2,max)
> idx_i=which((vx>min_poly[1])&(vx<max_poly[1]))
> idx_j=which((vy>min_poly[2])&(vy<max_poly[2]))
> sub_X=X[idx_i,idx_j]
> image(vx[idx_i],vy[idx_j],
+       log(sub_X+1),col=rev(heat.colors(101)),
+       xlab="",ylab="")
> lines(polygon_Fr)

We are now able to extract information about population for France, only (actually, it is only mainland France, islands are not considered here… to avoid complicated computations

> library(sp)
> xy=expand.grid(x = vx[idx_i], y = vy[idx_j])
> dim(xy)
[1] 65730     2

Here, we have 65,730 small squares, in France.

> pip=point.in.polygon(xy[,1],xy[,2],
+     polygon_Fr[,1],polygon_Fr[,2])>0
> dim(pip)=dim(sub_X)
> Fr=sub_X[pip]
> sum(Fr)
[1] 58105272

Observe that the total population within the French polygon is close to 60 million people, which is consistent with actual figures. Now, if we look more carefully at repartition over the French territory

> library(ineq)
> Gini(Fr)
[1] 0.7296936

Gini coefficient is rather high (over 70%), but it is also possible to visualize Lorenz curve,

> plot(Lc(Fr))

Observe that in 5% of the territory, we can find almost 54% of the population

> 1-min(LcF$L[LcF$p>.95])
[1] 0.5462632

In order to compare with other countries, consider the

> LC=function(rds="fr.rds"){
+ Fra=readRDS(rds)
+ n=length(Fra@polygons[[1]]@Polygons)
+ L=rep(NA,n)
+ for(i in 1:n) 
L[i]=nrow(Fra@polygons[[1]]@Polygons[[i]]@coords)
+ idx=which.max(L)
+ polygon_Fr=
+      Fra@polygons[[1]]@Polygons[[idx]]@coords
+ min_poly=apply(polygon_Fr,2,min)
+ max_poly=apply(polygon_Fr,2,max)
+ idx_i=which((vx>min_poly[1])&(vx<max_poly[1]))
+ idx_j=which((vy>min_poly[2])&(vy<max_poly[2]))
+ sub_X=X[idx_i,idx_j]
+ xy=expand.grid(x = vx[idx_i], y = vy[idx_j])
+ dim(xy)
+ pip=point.in.polygon(xy[,1],xy[,2],
+     polygon_Fr[,1],polygon_Fr[,2])>0
+ dim(pip)=dim(sub_X)
+ Fr=sub_X[pip]
+ return(list(gini=Gini(Fr),LC=Lc(Fr))
+ }
> FRA=LC()

For instance, consider Germany, or Italy

> deu=download.file(
"http://biogeo.ucdavis.edu/data/gadm2.8/rds/DEU_adm0.rds","deu.rds")
> DEU=LC("deu.rds")
> ita=download.file(
"http://biogeo.ucdavis.edu/data/gadm2.8/rds/ITA_adm0.rds","ita.rds")
> ITA=LC("ita.rds")

It is possible to plot Lorenz curve, together,

> plot(FRA$LC,col="blue")
> lines(DEU$LC,col="black")
> lines(ITA$LC,col="red")

Observe that France is clearly below the other ones. Compared with Germany, there is a significant difference

> FRA$gini
[1] 0.7296936
> DEU$gini
[1] 0.5088853

More precisely, if 54% of French people live in 5% of the territory, only 40% of Italians, and 32% of the Germans,

> 1-min(FRA$LC$L[FRA$LC$p>.95])
[1] 0.5462632
> 1-min(ITA$LC$L[ITA$LC$p>.95])
[1] 0.3933227
> 1-min(DEU$LC$L[DEU$LC$p>.95])
[1] 0.3261124

Exercises in R

Following the course we had on Monday, I upload here some exercises – that I did send a few weeks ago by email – an some correction. I can also mention a nice article, to scrape or not to scrape: technical and ethical challenges …

The final part of the course, on graphs and regression models will be in May. I will try to upload slides on the regression part as soon as possible.

How long could it take to run a regression

This afternoon, while I was discussing with Montserrat (aka @mguillen_estany) we were wondering how long it might take to run a regression model. More specifically, how long it might take if we use a Bayesian approach. My guess was that the time should probably be linear in , the number of observations. But I thought I would be good to check.

Let us generate a big dataset, with one million rows,

> n=1e6
> X=runif(n)
> Y=2+5*X+rnorm(n)
> B=data.frame(X,Y)

Consider as a benchmark the standard linear regression,

> lm_freq = function(n){
+   idx = sample(1:1e6,size=n)
+   reg = lm(Y~X,data=B[idx,])
+   summary(reg)
+ }

Here the regression is a subset of smaller size. We can do the same with a Bayesian approach, using stan,

> stan_lm ="
+ data {
+ int N;
+ vector[N] x;
+ vector[N] y;
+ }
+ parameters {
+ real alpha;
+ real beta;
+ real tau;
+ }
+ transformed parameters {
+ real sigma;
+ sigma <- 1 / sqrt(tau);
+ }
+ model{
+ y ~ normal(alpha + beta * x, sigma);
+ alpha ~ normal(0, 10);
+ beta ~ normal(0, 10);
+ tau ~ gamma(0.001, 0.001);
+ }
+ "

Define then the model

> library(rstan)
> system.time( 
  stanmodel <<- stan_model(model_code = stan_lm))
utilisateur     système      écoulé 
      0.043       0.000       0.043

We want to see how long it might take to run a regression,

> lm_bayes = function(n){
+   idx = sample(1:1e6,size=n)
+   fit = sampling(stanmodel,
+       data = list(N=n,
+                   x=X[idx],
+                   y=Y[idx]),
+       iter = 1000, warmup=200)
+   summary(fit)
+ }

We use the following package to see how long it takes

> library(microbenchmark)
> time_lm = function(n){
+  M = microbenchmark(lm_freq(n),
+      lm_bayes(n),times=50)
+  return(apply( matrix(M$time,nrow=2),1,mean))
+ }

We can now compare the time it took with ten, one hundred, on thousand, and ten thousand observations,

> vN = c(10,100,1000,10000)
> T = Vectorize(time_lm)(vN)

we can then plot it

> plot(vN,T[2,]/1e6,log="xy",col="red",type="b",
+      xlab="Number of Observations",ylab="Time")
> lines(vN,T[1,]/1e6,col="blue",type="b")

It looks like (if we forget about the very small sample) that the time it takes to run a regression is linear, with the two techniques (the frequentist and the Bayesian ones).

And actually, the same story olds for logistic regressions. Consider the following dataset

> n=1e6
> X=runif(n)
> S=-3+2*X+rnorm(n)
> Y=rbinom(n,size=1,prob=exp(S)/(1+exp(S)))
> B=data.frame(X,Y)

The frequentist version of the logistic regression is

> glm_freq = function(n){
+   idx = sample(1:1e6,size=n)
+   reg = glm(Y~X,data=B[idx,],family=binomial)
+   summary(reg)
+ }

and the Bayesian one, using stan,

> stan_glm = "
+ data {
+ int N;
+ vector[N] x;
+ int<lower=0,upper=1> y[N];
+ }
+ parameters {
+ real alpha;
+ real beta;
+ }
+ model {
+ alpha ~ normal(0, 10);
+ beta ~ normal(0, 10);
+ y ~ bernoulli_logit(alpha + beta * x);
+ }
+ "
> stanmodel = stan_model(model_code = stan_glm) )
> glm_bayes = function(n){
+   idx = sample(1:1e6,size=n)
+   fit = sampling(stanmodel,
+        data = list(N=n,
+        x = X[idx],
+        y = Y[idx]),
+        iter = 1000, warmup=200)
+   summary(fit)
+ }

Again, we can see how long it takes to run those regression models

> time_gl m= function(n){
+   M = microbenchmark(glm_freq(n),
+   glm_bayes(n),times=50)
+   return(apply( matrix(M$time,nrow=2),1,mean))
+ }