Tag Archives: France

The U.S. Has Been At War 222 Out of 239 Years

This morning, I discovered an interesting statistic, “America Has Been At War 93% of the Time – 222 Out of 239 Years – Since 1776“,  i.e. the U.S. has only been at peace for less than 20 years total since its birth. I wanted to check, get a better understanding and look at other countries in the world.

As always, we can try to extract information from wikipedia, since there are pages dedicated to that information

url="https://en.wikipedia.org/wiki/List_of_wars_involving_the_United_States"
download.file(url,destfile = "warUS.html")
url="https://en.wikipedia.org/wiki/List_of_wars_involving_France"
download.file(url,destfile = "warFR.html")
url="https://fr.wikipedia.org/wiki/Liste_des_guerres_de_la_France#Premi.C3.A8re_R.C3.A9publique"
download.file(url,destfile = "guerre.html")
url="https://en.wikipedia.org/wiki/List_of_wars_involving_Canada"
download.file(url,destfile = "warCAN.html")

If we look at the US page, there are tables, so it should be easy to extract it. For instance,

Even if the war did last 1 day, we will say that the US were at war in 1811. The information we want to confirm can be “there were 21 full years – from Jan 1st till Dec 31st – where the US were not at war, once, during those years“. From the row above, we can claim that the US were at war in 1811. Most of the time, we have

I.e. there is a beginning (here 1775) and an end (1783). So here, the US are said to be at war in 1775, 1776, 1777, 1778, 1779, 1780, 1781, 1782, 1783. To extract the information, we look for regular expressions in the first column, with number, on 4 digits.

https://freakonometrics.hypotheses.org/files/2017/03/guerre-us1.png

Well, sometimes it can be a bit tricky, since we have 3 dates, 1941, 1945 and (in the legend) 1944. But if we consider the minimal and the maximal dates, we have our range of dates.

Now that we we how to extract information, let’s do it. The code will be

library(stringr)
ext_date=function(x){
dates12="[0-9]{4}"
#grep(pattern = dates2, x = col1[1])
L=str_extract_all(as.character(x),dates12)
return_L=list()
if(length(L)>0){
for(j in 1:length(L))
if(length(L[[j]])==1) return_L[[j]]=as.numeric(L[[j]])
if(length(L[[j]])>=2) return_L[[j]]=seq(min(as.numeric(L[[j]])),max((as.numeric(L[[j]]))))
}
return(return_L)}

For the US, we get the following years

library(XML)
tables=readHTMLTable("warUS.html")
list_dates=list()
for(i in 1:length(tables)){
if(!is.null(dim(tables[[i]]))){
if(ncol(tables[[i]])>1){
col1=tables[[i]][,1]
list_dates[[i]]=lapply(col1,ext_date)
}
}}
d=unique(unlist(list_dates))

(red means at war, while green means no-war) and indeed,

> length(d)
[1] 222

there were 222 years with war.  Now, what about another country. Like France. Here I use the French wiki page, since information is not in tables in the English one.

tables=readHTMLTable("guerre.html")
list_dates=list()
for(i in 1:length(tables)){
if(!is.null(dim(tables[[i]]))){
if(ncol(tables[[i]])>1){
col1=tables[[i]][,1]
col2=tables[[i]][,2]
col12=paste(col1,col2)
list_dates[[i]]=lapply(col12,ext_date)
}
}}
d=unique(unlist(list_dates))

On the same period of time (starting in 1775), France was also on war most of the time.

Less than the US, but still: 185 years with war,

> length(d[d>=1775])
[1] 185

And on a longer period of time? Why not start, say, around the Hundred Years’s War,

meaning that since 1337, there were (only) 174 years without a single war where France was involved.

Let’s try another one. Like Canada,

tables=readHTMLTable("warCAN.html")
list_dates=list()
for(i in 1:length(tables)){
if(!is.null(dim(tables[[i]]))){
if(ncol(tables[[i]])>1){
col1=tables[[i]][,1]
list_dates[[i]]=lapply(col1,ext_date)
}
}}
d=unique(unlist(list_dates))

Guess what… there’s a lot of green on that graph. Surprised?

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

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

Spatial and Temporal Viz of Gas Price, in France

A great think in France, is that we can play with a great database with gas price, in all gas stations, almost eveyday. The file is rather big, so let’s make sure we have enough memory to run our codes,

> rm(list=ls())

To extract the data, first, we should extract the xml file, and then convert it in a more common R object (say a list)

> year=2014
> loc=paste("http://donnees.roulez-eco.fr/opendata/annee/",year,sep="")
> download.file(loc,destfile="oil.zip")

Content type 'application/zip' length 15248088 bytes (14.5 MB)

> unzip("oil.zip", exdir="./")
> fichier=paste("PrixCarburants_annuel_",year,
".xml",sep="")
> library(plyr)
> library(XML)
> library(lubridate)
> l=xmlToList(fichier)

We have a large dataset, with prices, for various types of gaz, for almost any gas station in France, almost every day, in 2014. It is a 1.4Gb list, with 11,064 elements (each of them being a gas station)

> length(l)
[1] 11064

There are two ways to look at the data. A first idea is to consider a gas station, and to extract the time series.

> time_series=function(no,type_gas="Gazole"){
+   prix=list()
+   date=list()
+   nom=list()
+   j=0
+   for(i in 1:length(l[[no]])){
+     v=names(l[[no]])
+     if(!is.null(v[i])){
+       if(v[i]=="prix"){
+         j=j+1
+  date[[j]]=as.character(l[[no]][[i]]["maj"])
+  prix[[j]]=as.character(l[[no]][[i]]["valeur"])
+  nom[[j]]=as.character(l[[no]][[i]]["nom"])
+       }}
+   }
+   id=which(unlist(nom)==type_gas)
+   n=length(id)
+   jour=function(j) as.Date(substr(date[[id[j]]],1,10),"%Y-%m-%d")
+   jour_heure=function(j) as.POSIXct(substr(date[[id[j]]],1,19), format = "%Y-%m-%d %H:%M:%S", tz = "UTC")
+   ext_y=function(j) substr(date[[id[j]]],1,4)
+   ext_m=function(j) substr(date[[id[j]]],6,7)
+   ext_d=function(j) substr(date[[id[j]]],9,10)
+   ext_h=function(j) substr(date[[id[j]]],12,13)
+   ext_mn=function(j) substr(date[[id[j]]],15,16)
+   prix_essence=function(i) as.numeric(prix[[id[i]]])/1000
+   base1=data.frame(indice=no,
+            id=l[[no]]$.attrs["id"],
+            adresse=l[[no]]$adresse,
+            ville=l[[no]]$ville,
+  lat=as.numeric(l[[no]]$.attrs["latitude"])
/100000,
+  lon=as.numeric(l[[no]]$.attrs["longitude"])
/100000,
+       cp=l[[no]]$.attrs["cp"],
+       saufjour=l[[no]]$ouverture["saufjour"], 
+       Y=unlist(lapply(1:n,ext_y)),
+       M=unlist(lapply(1:n,ext_m)),
+       D=unlist(lapply(1:n,ext_d)),
+       H=unlist(lapply(1:n,ext_h)),
+       MN=unlist(lapply(1:n,ext_mn)),
+    prix=unlist(lapply(1:n,prix_essence)))
+   
+   base1=base1[!is.na(base1$prix),]
+   
+   date_d=paste(year,"-01-01 12:00:00",sep="")
+   date_f=paste(year,"-12-31 12:00:00",sep="")
+   vecteur_date=seq(as.POSIXct(date_d, format =
+                 "%Y-%m-%d %H:%M:%S"),
+                    as.POSIXct(date_f, format = 
+                 "%Y-%m-%d %H:%M:%S"),by="days")
+   date=paste(base1$Y,"-",base1$M,"-",base1$D,
+   " ",base1$H,":",base1$MN,":00",sep="")
+   date_base=as.POSIXct(date, format = 
+                "%Y-%m-%d %H:%M:%S", tz = "UTC")
+   idx=function(t) sum(vecteur_date[t]>=date_base)
+   vect_idx=Vectorize(idx)(1:length(vecteur_date))
+   P=c(NA,base1$prix)
+   base2=ts(P[1+vect_idx],
+         start=year,frequency=365)
+   list(base=base1,
+        ts=base2)
+ }

To get the time series, extrapolation is necessary, since we have here observation at irregular dates. Here, for instance, for the second gas station, we get

> plot(time_series(2)$ts,ylim=c(1,1.6),col="red")
> lines(time_series(2,"SP98")$ts,col="blue")

An alternative is to study gas price from a spatial perspective. Given a date, we want the price in all stations. As previously, we keep the last price observed, in each station,

> spatial=function(dt){
+   base=NULL
+   for(no in 1:length(l)){  
+     prix=list()
+     date=list()
+     j=0
+     for(i in 1:length(l[[no]])){
+     v=names(l[[no]])
+     if(!is.null(v[i])){
+       if(v[i]=="prix"){
+   j=j+1
+   date[[j]]=as.character(l[[no]][[i]]["maj"])
+       }}
+   }
+   n=j
+   D=as.Date(substr(unlist(date),1,10),"%Y-%m-%d")
+   k=which(D==D[which.max(D[D<=dt])])
+ if(length(k)>0){
+   B=Vectorize(function(i) l[[no]][[k[i]]])(1:length(k))
+ if("nom" %in%  rownames(B)){  
+   k=which(B["nom",]=="Gazole")
+   prix=as.numeric(B["valeur",k])/1000
+   if(length(prix)==0) prix=NA
+   base1=data.frame(indice=no,
+   lat=as.numeric(l[[no]]$.attrs["latitude"])
/100000,
+   lon=as.numeric(l[[no]]$.attrs["longitude"])
/100000,
+   gaz=prix)
+   base=rbind(base,base1)
+ }}}
+ return(base)}

For instance, for the 5th of May, 2014, we get the following dataset

> B=spatial(as.Date("2014-05-05"))

To visualize prices, consider only mainland France (excluding islands in the Pacific, or close to the Caribeans)

> idx=which((B$lon>(-10))&(B$lon<20)&
+ (B$lat>35)&(B$lat<55))
> B=B[idx,]
> Q=quantile(B$gaz,seq(0,1,by=.01),na.rm=TRUE)
> Q[1]=0
> x=as.numeric(cut(B$gaz,breaks=unique(Q)))
> CL=c(rgb(0,0,1,seq(1,0,by=-.025)),
+ rgb(1,0,0,seq(0,1,by=.025)))
> plot(B$lon,B$lat,pch=19,col=CL[x])

Red dots are the most expensive gas stations, that particular day.

If we add contours of the French regions, we get

> library(maps)
> map("france")
> points(B$lon,B$lat,pch=19,col=CL[x])

 

We can also focus on some specific region, say the South of Brittany.

> library(OpenStreetMap)
> map <- openmap(c(lat= 48,   lon= -3),
+                c(lat= 47,   lon= -2))
> map <- openproj(map) 
> plot(map)
> points(B$lon,B$lat,pch=19,col=CL[x])

As we can see on that map, there are regions that are rather empty, where the closest gas station might be a bit far away. Actually, it is possible to add Voronoi sets on the map,

> dB=data.frame(lon=B$lon,lat=B$lat)
> idx=which(!duplicated(dB))
> dB=dB[idx,]

 

which could help to get the price of the closest gaz station.

> library(tripack)
> V <- voronoi.mosaic(dB$lon[id],dB$lat[id])
> plot(V,add=TRUE)

It is possible to plot each polygon with the color of the gaz station we add. Actually, it is a bit tricky, and I could not find a R function to to this. So I did it manually,

> plot(map)
> P <- voronoi.polygons(V)
> library(sp)
> point_in_i=function(i,point) point.in.polygon(point[1],point[2],P[[i]][,1],P[[i]][,2])
> which_point=function(i) which(Vectorize(function(j) point_in_i(i,c(dB$lon[id[j]],dB$lat[id[j]])))(1:length(id))>0)
> for(i in 1:length(P)) polygon(P[[i]],col=CL[x[id[which_point(i)]]],border=NA)

With this map, we can see that we have blue areas, i.e. all stations in a given area are cheap (because of competition), but in some places, a very expensive one is next to a very cheap one. I guess we should look closer at the dynamics… [to be continued….]

Clusters of (French) Regions

For the data scienec course of tomorrow, I just wanted to post some functions to illustrate cluster analysis. Consider the dataset of the French 2012 elections

> elections2012=read.table(
"http://freakonometrics.free.fr/elections_2012_T1.csv",sep=";",dec=",",header=TRUE)
> voix=which(substr(names(
+ elections2012),1,11)=="X..Voix.Exp")
> elections2012=elections2012[1:96,]
> X=as.matrix(elections2012[,voix])
> colnames(X)=c("JOLY","LE PEN","SARKOZY","MÉLENCHON","POUTOU","ARTHAUD","CHEMINADE","BAYROU","DUPONT-AIGNAN","HOLLANDE")
> rownames(X)=elections2012[,1]

The hierarchical cluster analysis is obtained using

> cah=hclust(dist(X))
> plot(cah,cex=.6)

To get five groups, we have to prune the tree

> rect.hclust(cah,k=5)
> groups.5 <- cutree(cah,5)

We have to zoom-in to visualize the French regions,

It is also possible to use

> library(dendroextras)
> plot(colour_clusters(cah,k=5))

And again, if we zoom-in, we get

The interpretation of the clusters can be obtained using

> aggregate(X,list(groups.5),mean)
  Group.1     JOLY   LE PEN  SARKOZY
1       1 2.185000 18.00042 28.74042
2       2 1.943824 23.22324 25.78029
3       3 2.240667 15.34267 23.45933
4       4 2.620000 21.90600 34.32200
5       5 3.140000  9.05000 33.80000

It is also possible to visualize those clusters on a map, using

> library(RColorBrewer)
> CL=brewer.pal(8,"Set3")
> carte_classe <- function(groupes){
+ library(stringr)
+ elections2012$dep <- elections2012[,2]
+ elections2012$dep <- tolower(elections2012$dep)
+ elections2012$dep <- str_replace_all(elections2012$dep, pattern = " |-|'|/", replacement = "")
+ library(maps)
+ france<-map(database="france")
+ france$dep <- france$names
+ france$dep <- tolower(france$dep)
+ france$dep <- str_replace_all(france$dep, pattern = " |-|'|/", replacement = "")
+ corresp_noms <- elections2012[, c(1,2, ncol(elections2012))]
+ corresp_noms$dep[which(corresp_noms$dep %in% "corsesud")] <- "corsedusud"
+ col2001<-groupes+1
+ names(col2001) <- corresp_noms$dep[match(names(col2001), corresp_noms[,1])]
+ color <- col2001[match(france$dep, names(col2001))]
+ map(database="france", fill=TRUE, col=CL[color])
+ }
> carte_classe(cutree(cah,5))

or, if we simply want 4 clusters

> carte_classe(cutree(cah,4))

 

Douce France

Hier, sur Twitter, je mentionnais une carte parue sur Times Higher Education à propos des “100 meilleurs universités” (dans le monde entier). Je ne sais plus par quel biais j’ai vu des cartes, par pays, permettant de voir où étaient les 100 meilleures universités (d’ordinaire, on a un classement, mais là, pour la première fois, on avait des cartes),

Le premier tweet, en anglais, qui montrait la carte est passé assez inaperçu. Ce n’est pas le cas du second, en français, qui a été beaucoup repris, dans lequel je posais la question de manière plus directe, “sauras-tu remarquer une particularité géographique de l’enseignement supérieur français?“. Car je trouve la carte française assez choquante. disait de manière assez élégante que la “ville lumière” faisait de l’ombre aux autres villes françaises. Ce n’est pas peu dire ! Dans tous les autres pays comparables, avec 4, 5 ou 6 universités dans ce classement, on a une relative dispersion sur le territoire. Mais en France, on a Paris. Et c’est tout.

C’est amusant car cette semaine, j’ai commencé le livre Unknown Quantity de John Derbyshire qui raconte, de manière très vulgarisée (et donc parfaite pour moi) l’histoire de l’algèbre. Et dans un des premiers chapitres (on est encore en Mésopotamie) il y a ce paragraphe qui m’a beaucoup rappelé cette carte

Effectivement, il faudrait que je creuse un peu, mais j’ai du mal à voir ce qui va ressortir de cette méga-concentration géographique typiquement française de l’enseignement supérieur. Rien de bon, je crains…

Analyse des Données et Cartes

Vendredi, on continuera en cours la classification (non supervisée) et en particulier les méthodes hiérarchiques, et les arbres. On va utiliser la base suivante, avec les résultats des élections présidentielles (premier tour) de 2012,

elections2012 <- read.csv("http://komodo.regardscitoyens.org/public/presidentielles2012/resultats/resultats_departements_final_T1.csv",
sep=";",header=TRUE,dec=",")

La base qu’on garde contient juste les voix exprimées, en pourcentage

voix <- which(substr(names(elections2012),1,12)=="X..Voix.Exp.")
X <- as.matrix(elections2012[,voix])
colnames(X) <- substr(names(elections2012)[voix],13,nchar(names(elections2012)[voix]))
rownames(X) <- elections2012[,3]

Parfois, il convient de normaliser les données. Ici, sur les département, ça n’est pas nécessaire. Par contre, sur les candidats, ça le serait (comme pour l’ACP en fait).  Pour visualiser les distances (entre départements), on peut utiliser

heatmap(X)

On peut faire une CAH, sur notre matrice de distance entre rangs

cah <- hclust(dist(X))
plot(cah,cex=.6)

et si on souhaite garder 5 groupes (par exemple) on utilise

rect.hclust(cah,k=5)
groups.5 <- cutree(cah,5)

On peut aussi visualiser les classes avec la fonction

library(dendroextras)
plot(colour_clusters(cah,k=5))

Pour comprendre qui se trouve dans nos groupes, on peut utiliser

aggregate(X,list(groups.5),mean)

qui va nous renvoyer les “votes  moyens” pour chaque candidat, dans chaque groupe. Maintenant, pour en finir avec la visualisation, on peut utiliser le code suivant (un peu long, certes) qui permettra de visualiser sur un carte les différents groupes

carte_classe <- function(groupes){
library(stringr)
elections2012$dep <- elections2012$Libellé.du.département
elections2012$dep <- tolower(elections2012$dep)
elections2012$dep <- str_replace_all(elections2012$dep, pattern = " |-|'|/", replacement = "")
library(maps)
france<-map(database="france")
france$dep <- france$names
france$dep <- tolower(france$dep)
france$dep <- str_replace_all(france$dep, pattern = " |-|'|/", replacement = "")
corresp_noms <- elections2012[, c("Libellé.du.département", "dep")]
corresp_noms$dep[which(corresp_noms$dep %in% "corsesud")] <- "corsedusud"
col2001<-groupes+1
names(col2001) <- corresp_noms$dep[match(names(col2001), corresp_noms$Libellé.du.département)]
color <- col2001[match(france$dep, names(col2001))]
map(database="france", fill=TRUE, col=color)
}

carte_classe(groups.5)

Avec ces fonctions, on devrait pouvoir tester différentes méthodes pour constituer des groupes.

Men set to live as long as women by 2030?

A few months ago, in Men set to live as long as women, figures show, it was mentioned that (in the U.K.)

the gap between male and female life expectancy is closing and men could catch up by 2030, according to an adviser for the Office for National Statistics.

(the slides are available online http://cass.city.ac.uk/…).

Continue reading Men set to live as long as women by 2030?

French dataset: population and GPS coordinates

A short post today based on recent work by @3wen (Ewen Gallic, graduate Student in Rennes, spending a year in Montreal). Since we were working on a detailed French dataset (per commune), we needed a dataset containing a list allcommunes, with population and location. GPS coordinates were extracted from Google, using the following php file, inspired by http://www.andrew-kirkpatrick.com/ on Google geocoding api with php webpage. Population was interpolated from INSEE’s datasets, i.e. http://www.insee.fr/ (since data are over a 35 year period, from 1975 to 2010, changes have been taken into account as carefully are possible – e.g. merges and splits of cities – based on thatdescription). A spline model has been used for all cities (with three degrees of freedom, and null and negative interpolation became one, since we’ll be using loglinear models afterwards). Names are from that dataset, still on INSEE’s website, http://www.insee.fr/.

A zipped file can be downloaded here popfr19752010.zip, but it is also possible to use the code below (it is a 24Mo dataset). Since it was hard to find such a dataset online (different files can be found, but we found none with population and location), we have decided to upload that dataset. Please let us know if there are problems with those data…

> base=read.csv(
+ "http://freakonometrics.free.fr/popfr19752010.csv",
+ header=TRUE)

Using that code, it is possible to locate all the communes in France (metropolitan), for instance

> library(maps)
> map("france")
> points(base$long,base$lat,cex=.1,col="red",pch=19)
> points(base$long,base$lat,cex=2*base$pop_2010/
+ max(base$pop_2010),col="blue",pch=19)

Several additional lines of code on that dataset (and also others) will be uploaded, soon.

Cette oeuvre est mise à disposition sous licence Paternité – Partage à l’Identique 3.0 non transposé. Pour voir une copie de cette licence, visitez http://creativecommons.org/. Date : 24 mai 2012, par Ewen GALLIC. Sources : INSEE, API Google Maps v3 et GeoHack (coordonnées GPS), propres calculs (estimation de population à partir des données INSEE).

  • reg : code region INSEE (character)
  • dep : code departement INSEE (character, corse 201 et 202 au lieu de 2A et 2B)
  • com : code commune INSEE (character)
  • article : article du nom de la commune (character)
  • com_nom : nom de la commune (character)
  • long : longitude (numeric)
  • lat : latitude (numeric)
  • pop_i : estimation de la population à la date i (ramenée à 1 si <=0), i=1975,…,2010 (numeric)

Births and week-ends, in France

This week, I have seen on the internet (sorry, I cannot find proper references) the graph produced here on the right: which birthday is most likely ? The fact that I have no further information is important, since I do not know in which country such a graph was obtained. At least, I know it should not be France…

In France, I have already mentioned that there is a strong week-end effect: nowadays, there is 25% less deliveries during week-ends than during the week. Calot (1981) observed already that there were less deliveries on Sundays. This has been confirmed more recently, e.g. in http://www.lepoint.fr/ or http://www.prepabl.fr/, with a significant difference between week days, and week-ends. Here  is the number of birth per day, over 40 years, with in blue the average trend during the week, and in red, during week-ends,

naissance=read.table(
"http://freakonometrics.free.fr/naissanceFR2.txt")
attach(naissance)
date=as.Date(date)
plot(date, nbre,cex=.5)
t2=as.POSIXlt(date)
jour=t2$wday
X=naissance$date
Y=naissance$nbre
J=jour
df=data.frame(X,Y,J)
library(splines)
regs=lm(Y~bs(X,df=20),data=df[jour%in%c(0,6),])
Yp=predict(regs,newdata=df)
lines(X,Yp,col="red",lwd=3)
regs=lm(Y~bs(X,df=20),data=df[jour%in%1:5,])
Yp=predict(regs,newdata=df)
lines(X,Yp,col="blue",lwd=3)

If we look at the evolution of the ratio week-ends over weeks days, we have the following graph

t2=as.POSIXlt(date)
jour=t2$wday
jour=jour[1:(1982*7)]
nbre2=jour
for(i in 1:1982){
taux=sum(nbre[6:7+7*(i-1)])/
sum(nbre[1:5+7*(i-1)])/2*5
nbre2[1:5+7*(i-1)]=nbre[1:5+7*(i-1)]*taux
nbre2[6:7+7*(i-1)]=nbre[6:7+7*(i-1)]
nbre2[1:7+7*(i-1)]=
mean(nbre[1:7+7*(i-1)])/mean(nbre2[1:7+7*(i-1)])*
nbre2[1:7+7*(i-1)]
}
nbretaux=jour
for(i in 1:1982){
taux=sum(nbre[6:7+7*(i-1)])/
sum(nbre[1:5+7*(i-1)])/2*5
nbretaux[1:7+7*(i-1)]=taux
}
plot(date[1:length(nbre2)],nbretaux)
X= date[1:length(nbre2)]
Y=nbretaux
library(splines)
reg=lm(Y~bs(X,df=20))
Yp=predict(reg)
lines(X,Yp,col="red",lwd=3)

In the beginning of the 70’s, during week-ends, there were 5% less deliveries, but 25% less around 2000. It is then possible to produce the same kind of graphs as the one above, per year of birth. And here, we clearly observe the importance of the week end effect (maybe also because of color choice)

naissance=read.csv(
"http://freakonometrics.free.fr/naissanceFR.csv",
sep=";")
M=as.matrix(naissance[,3:ncol(naissance)])
BIRTH=as.vector(t(M))
YEAR=rep(1968:2005,each=12*31)
MONTH=rep(rep(1:12,each=31),38)
DAY=rep(1:31,12*38)
X=NA
for(y in 1968:2005){
sbase=base[YEAR==y,]
X=c(X,sbase$BIRTH/sum(sbase$BIRTH,
na.rm=TRUE))
}
base=data.frame(YEAR,MONTH,DAY,
BIRTH,BIRTHDAYPROB=X[-1])

m1=min(base$BIRTHDAYPROB,na.rm=TRUE)
m2=max(base$BIRTHDAYPROB,na.rm=TRUE)
y=1980
colr=rev(heat.colors(100))
sbase=base[YEAR==y,]
plot(0:1,0:1,col="white",xlim=c(-1,12),
ylim=c(-31,1),axes=FALSE,xlab=
paste("Naissance en",y,sep=" "),ylab="")
for(x in 1:nrow(sbase)){
a=sbase$MONTH[x];b=sbase$DAY[x]
polygon(c(a-.9,a-.9,a-.1,a-.1),-c(b-.9,b-.1,
b-.1,b-.9),col=colr[(sbase$BIRTHDAYPROB[x]-m1)/
(m2-m1)*100],border=NA)
}
text((1:12)-.5,.5,c("J","F","M","A","M","J","J",
"A","S","O","N","D"),cex=.7)
text(-.5,-(1:31)+.5,1:31,cex=.7)

Une région géographique n’est pas une variable continue

En relisant les devoirs maisons, je me suis rendu compte que certains avaient tenté de regrouper les régions (géographiques) par régions homogènes. Sauf que les régions étaient codées par un numéro (selon la codification officielle). Par exemple, dans une des bases, nous avions des assurés dans 4 zones géographiques, à savoir la région 82 (région Rhône-Alpes en rouge) la région 54 (région Poitou-Charentes en vert) la région 73 (région Midi-Pyrénées en bleu) et enfin la région 41 (région Lorraine en mauve).

> unique(baseFREQ$region) 
[1] 82 54 73 41

Une idée intéressante pour regrouper les régions pouvait être d’utiliser les arbres. Les régions étant des couleurs (on le voit bien sur la carte) et pas des variables quantitatives, il est normal de travailler sur des facteurs. D’ailleurs le code pour faire la carte est le suivant,

> library(maps) 
>  france<-map(database="france") 
>  dpt=c("Ain","Ardeche","Drome","Isere","Loire ","Rhone",  
+ "Savoie","Haute-Savoie","Charente","Charente-Maritime", 
+ "Deux-Sevres","Vienne","Ariege","Aveyron","Haute-Garonne",  
+ "Gers","Lot","Hautes-Pyrenees","Tarn","Tarn-et-Garonne", 
+ "Meurthe-et-Moselle","Meuse","Moselle","Vosges") 
>  couleur=c(rep(2,8),rep(3,4),rep(4,8),rep(6,4))  
>  match=match.map(france,dpt) 
>  color=couleur[match] 
>  map(database="france", fill=TRUE, col=color)

L’arbre sur les régions en tant que facteurs donne le découpage suivant

>  baseFREQ$fregion=as.factor(baseFREQ$region) 
>  ARBRE1=tree(nombre~fregion,data=baseFREQ,split="gini")  
>  plot(ARBRE1) 
>  text(ARBRE1)

Bon, R a la mauvaise idée de recoder les classes (mais il garde l’ordre, i.e. a correspond à la région 41, b à 54, c à 73 et d à 82). Visuellement, on retient qu’il est possible de considérer deux grandes régions, AC (i.e. 41 et 73) et BD (i.e. 54 et 72). L’intérêt des arbres sur des variables qualitatives, des facteurs, c’est que tous les regroupements sont possibles. En revanche, si on fait un arbre sur la région qui est lue en tant que nombre (quantitatif), on obtient

>  ARBRE2=tree(nombre~region,data=baseFREQ,split="gini") 
>  plot(ARBRE2) 
>  text(ARBRE2)

Il est alors impossible de regrouper dans une même classe deux régions séparées par un nombre, i.e. on ne peut regrouper 41 et 82 dans la même classe. R suggère de distinguer peut être trois régions, à savoir 82 (à droite), puis 73 (au centre) et enfin de mettre éventuellement 41 et 54 ensemble. Ce qui n’est pas la stratégie optimale quand on regroupe des facteurs.

Climate change and insurance

I will be in Lyon next Monday to give a talk on “Modeling heat-waves: return period for non-stationary extremes” in a workshop entitled “Changement climatique et gestion des risques“. An interesting reference might be some pages from Le Monde (2010). The talk will be more a discussion about modeling series of temperatures (daily temperatures). A starting point might be the IPCC Third Assessment graph which illustrates the effect on extreme temperatures when (a) the mean temperature increases, (b) the variance increases, and (c) when both the mean and variance increase for a normal distribution of temperature.

I will add here some code used to generate some graphs I will comment. The graph below it the daily minimum temperature,

TEMP=read.table("http://freakonometrics.blog.free.fr/
public/data/TN_STAID000038.txt",header=TRUE,sep=",")
D=as.Date(as.character(TEMP$DATE),"%Y%m%d")
T=TEMP$TN/10
day=as.POSIXlt(D)$yday+1
an=trunc(TEMP$DATE/10000)
plot(D,T,col="light blue",xlab="Minimum
daily temperature in Paris",ylab="",cex=.5)
abline(R,lwd=2,col="red")

We can clearly see an increasing linear trend. But we do not care (too much) here about the increase of the average temperature, but more dispersion, and tails. Here are decenal box-plots

or quantile-regressions

library(quantreg)
PENTESTD=PENTE=rep(NA,99)
for(i in 1:99){
R=rq(T~D,tau=i/100)
PENTE[i]=R$coefficients[2]
PENTESTD[i]=summary(R)$coefficients[2,2]
}
m=lm(T~D)$coefficients[2]
plot((1:99)/100,(PENTE/m-1)*100,type="b")
segments((1:99)/100,((PENTE-2*PENTESTD)/m-1)*100,
(1:99)/100,((PENTE+2*PENTESTD)/m-1)*100,
col="light blue",lwd=3)
points((1:99)/100,(PENTE/m-1)*100,type="b")
abline(h=0,lty=2,col="red")

In order to get a better understanding of the graph above, here are slopes of quantile regressions associated to different probabilities,

The annualized maxima (of minimum temperature, i.e. warmest night of the year)

i.e. the regression of yearly maximas.

tail index of a Generalized Pareto distribution

Instead of looking at observation over a century (the trend is obviously linear), we can focus on seaonal behavior,

B=data.frame(Y=rep(T,3),X=c(day,day-365,day+365),
A=rep(an,3))
library(quantreg)
library(splines)
Q50=rq(Y~bs(X,10),data=B,tau=.5)
Q95=rq(Y~bs(X,10),data=B,tau=.95)
Q05=rq(Y~bs(X,10),data=B,tau=.05)
YP95=predict(Q95,newdata=data.frame(X=1:366))
YP05=predict(Q05,newdata=data.frame(X=1:366))
I=(T>predict(Q95))|(T<predict(Q05))
YP50=predict(Q50,newdata=data.frame(X=1:366))
plot(day[I],T[I],col="light blue",cex=.5)
lines(1:365,YP95[1:365],col="blue")
lines(1:365,YP05[1:365],col="blue")
lines(1:365,YP50[1:365],col="blue",lwd=3)

with on red series from 1900 till 1920, and on purple from 1990 till 2010. If we remove the linear trend, and the seasonal cycle, here are the residuals, assume to be stationary,

on during the year

Obviously, something has been missed,

The graph below is the volatility of the residual series, within the year,

Instead of looking at volatility, we can focus on tails, with tail index per month,

mois=as.POSIXlt(D)$mon+1
Pmax=Dmax=matrix(NA,12,2)
for(s in 1:12){
X=T3[mois==s]
FIT=gpd(X,5)
Pmax[s,1:2]=FIT$par.ests
Dmax[s,1:2]=FIT$par.ses
}
plot(1:12,Pmax[,1],type="b",col="blue",
ylim=c(-.6,0))
segments(1:12,Pmax[,1]+2*Dmax[,1],1:12,Pmax[,1]-
2*Dmax[,1],col="light blue",lwd=2)
points(1:12,Pmax[,1],col="blue")
text(1:12,rep(-.5,12),c("JAN","FEV","MARS",
"AVR","MAI","JUIN","JUIL","AOUT","SEPT",
"OCT","NOV","DEV"),cex=.7)

At the end of the talk, I will also mention multiple city models, e.g. Paris and Marseille,

If we look at residuals (once we have removed the linear trend and the seasonal cycle) we observe some positive dependence

In order to study (strong) tail dependence, define

http://freakonometrics.hypotheses.org/files/2017/07/Llatex2png.2.php_.png

for lower left tail and

http://freakonometrics.hypotheses.org/files/2017/07/Clatex2png.2.php_.png

for upper right tail, where http://freakonometrics.hypotheses.org/files/2017/07/toclatex2png-12.2.php_.png is the survival copula associated to http://freakonometrics.hypotheses.org/files/2017/07/toclatex2png-13.2.php_.png, i.e.
http://freakonometrics.hypotheses.org/files/2017/01/toclatex2png-14.2.php_.png

and

http://freakonometrics.hypotheses.org/files/2017/01/toclatex2png-15.2.php_.png

It looks like there is no tail dependence (in the uper tail). But it is also possible to study weaker tail dependence, through

http://freakonometrics.hypotheses.org/files/2017/01/toc2latex2png.3.php_.png

and

http://freakonometrics.hypotheses.org/files/2017/01/toc2latex2png.4.php_.png


Slides can be visualized below, I will upload them soon,

Maudits français

Tous les malouins le savent, Jacques Cartier est parti de Saint Malo pour découvrir le Canada en 1536 (si l’on considère son second voyage, c’est à dire celui où il a vraiment pénétré le royaume de Kanata en voguant sur le Saint Laurent). Il y a d’ailleurs un musée du Québec dans les remparts. Et comme s’en vente le musée, les québécois viennent de Saint-Malo.

Or, pour ceux qui ne connaissent pas Saint-Malo, c’est à deux pas du Mont Saint Michel, et du Couesnon qui sépare la Bretagne et la Normandie. Étant né en Normandie, avant de venir vivre à Rennes,  j’étais convaincu que j’avais des racines communes avec les Québécois*. Mais j’ai été surpris en voyant les noms de famille des copines de classe de ma fille: les noms ne m’étaient pas du tout familiers… Pourtant, si nous avons des racines communes, il ne serait pas surprenant que les noms de familles se retrouvent (c’est le principe de base de la généalogie, si j’ai bien tout suivi). Aussi j’ai voulu creuser davantage….

  • Que disent les experts en généalogie ?
http://freakonometrics.blog.free.fr/public/perso3/aunis.gifhttp://freakonometrics.blog.free.fr/public/perso3/bretagne.gifhttp://freakonometrics.blog.free.fr/public/perso3/Saintonge.gifhttp://freakonometrics.blog.free.fr/public/perso3/Normandie.gifhttp://freakonometrics.blog.free.fr/public/perso3/anjou.gif

On peut trouver des éléments de réponse ici ou là. Parmi les 3800 Français qui ont immigré (ou émigrer, ça dépend du point de vue) entre 1608 et 1700 :

  • La Normandie a fourni 1 pionnier sur 5.
  • Le Poitou, l’Aunis et la Saintonge ont, ensemble, donné 1 pionnier sur 4.
  • Paris et l’Ile-de-France, 1 pionnier sur 7.
  • La Bretagne, un peu plus de 3 pionniers sur 100.
  • L’Anjou, 3 pionniers sur 100.
  • La Champagne, presque 3 pionniers sur 100.
  • La Picardie, un peu plus de 2 pionniers sur 100.

 

Entre 1700 et 1765, près de 5000 autres hommes, femmes et enfants originaires de France viennent s’établir au Canada, mais “très peu proviennent de la Bretagne et des régions baignées par la mer“. Ce sont d’ailleurs les ordres de grandeur que j’ai pu avoir dans l’ouvrage de Michel Lambert (qui n’est pas un livre de généalogie à proprement parler mais qui est incroyablement intéressant, même s’il ne donne pas de source – ou de détails sur le sens – pour ces chiffres),

  • Normandie, 22.6%
  • Aunis, 16.4%
  • Perche, 11.4%
  • Ile de France (région parisienne) 10.5%
  • Poitou, 7.5%
  • Maine, 5.2%
  • Saintonge 5.2%
  • Bretagne 3.5%

Moralité, même si Jacques Cartier est parti de Bretagne, rares semblent être les bretons qui ont émigrés au Québec.

  • Calculs de probabilités conditionnelles à partir des noms de famille

Il existe des sites qui donnent les “classements” des noms de familles, au Québec par exemple (ici, avec les 1000 noms les plus portés, avec les proportions respectives), ou dans les régions françaises (, où j’ai pris à chaque fois les 2000 noms les plus portés avec les proportions respectives, en changeant de région). Les données peuvent être récupérées avec le code suivant

quebec=read.table("http://freakonometrics.free.fr/nom-quebec2.txt",   
                   header=TRUE,dec=",")
bretagne=read.table("http://freakonometrics.free.fr/nom-bretagne.txt",
                    header=TRUE,dec=",",sep="\t")
poitou=read.table("http://freakonometrics.free.fr/nom-poitou.txt",
                  header=TRUE,dec=",",sep="\t")
normandie=read.table("http://freakonometrics.free.fr/nom-basse-normandie.txt",
                     header=TRUE,dec=",",sep="\t")
aquitaine=read.table("http://freakonometrics.free.fr/nom-aquitaine.txt",
                     header=TRUE,dec=",",sep="\t")
alsace=read.table("http://freakonometrics.free.fr/nom-alsace.txt",  
                  header=TRUE,dec=",",sep="\t")
loire=read.table("http://freakonometrics.free.fr/nom-loire.txt",
                 header=TRUE,dec=",",sep="\t")

J’ai retenu 6 régions françaises, et j’ai calculé la proportions de québécois dont le nom de famille était dans le top 2000 d’une région donnée, i.e.

> head(quebec)
  Rang Nomdefamille Pourcentage
1    1     Tremblay       1.076
2    2       Gagnon       0.790
3    3          Roy       0.753
4    4         Cote       0.692
5    5     Bouchard       0.530
6    6     Gauthier       0.522

> minquebec=tolower(quebec$Nomdefamille)
> rangquebec=quebec$Rang
> pctquebec=quebec$Pourcentage/sum(quebec$Pourcentage)
> head(bretagne)
  Rang Patronyme Naissances
1    1   LE GALL      18126
2    2   LE GOFF      16093
3    3   LE ROUX      15905
4    4    THOMAS      15705
5    5    MARTIN      12094
6    6    TANGUY      12045

> minbretagne=tolower(bretagne$Patronyme)
> rangbretagne=bretagne$Rang
> pctbretagne=bretagne$Naissances/sum(bretagne$Naissances)
> I=minquebec%in%minbretagne
> sum(pctquebec[I])/sum(pctquebec)
[1] 0.2579641
Région Proportion
Bretagne 23,39%
Normandie 31,07%
Loire 32,83%
Aquitaine 27,65%
Poitou 32,50%
Alsace 9,71%

Si on retrouve peu de noms alsacien au Québec, on peut être un peu “surpris” que les régions qui ressortent le plus sont la Loire et le Poitou et pas la Bretagne et la Normandie…. Bon, j’en conviens, je travaille sur les noms de familles bruts, c’est à dire que je ne tiens pas compte du fait que certains noms se sont déformés. Il s’agit de la version simple… on peut bien entendu aller plus loin, par exemple en utilisant des modèles de mélange par exemple… à suivre donc…

*en fait, mes origines sont davantage bourguignonnes, de part mes quatre grands parents, comme le savent tous ceux qui m’ont déjà payé un verre de vin… mais ayant passé toute mon enfance en Normandie, je peux me considérer un peu comme normand. Et adorant la galette saucisse, je suis définitivement breton….

Du sex-ratio en France

En novembre dernier, Baptiste @Coulmont m’avait envoyé un courriel correspondant à ce qu’il a mis en ligne sur son blog (ici) sur l’utilisation du fichier des prénoms français pour analyser le sex ratio à la naissance en France. J’attendais qu’il publie ses commentaires avant de mettre les miens (car le graphique qu’il a mis en ligne ce matin m’avait fait m’interroger).
Pour faire simple, la base de prénoms inclue un sexe. Mais si on regarde le rapport du nombre de garçon sur le nombre de filles, à la naissance, on a le graphique ci-dessous,

dat=read.table("nat2004.csv",sep=";",header=TRUE)
naissancesm=rep(NA,105)
for (i in 1900:2004) {
naissancesm[i-1899]=sum(dat[dat$annais==i&dat$sexe==1,"nombre"],
na.rm=TRUE)
}
naissancesf=rep(NA,105)
for (i in 1900:2004) {
naissancesf[i-1899]=sum(dat[dat$annais==i&dat$sexe==2,"nombre"],
na.rm=TRUE)
}
plot(1900:2004,naissancesm/naissancesf,col="red")

La tendance du début est très surprenante. Par exemple,  si on reprend les chiffres donnés par Pierre Simon Laplace sur les naissances à Paris entre 1750 et 1800 (mentionné hier, ici) on est déjà sur ratio de l’ordre de 1.05 (que l’on retrouve sur la fin de notre graphique, mais pas le début).
> 393386/377555
[1] 1.041930
 
> 251527/241945
[1] 1.039604

Donc il n’y a pas de raison d’avoir cette diff1rence. La conclusion semble être qu’il y a un soucis sur la base des prénoms. En effet, dans un rapport de Anouch Chahnazarian (ici), on retrouve l’évolution suivante,

que l’on peut rapprocher des données d’Éric Brian & Marie Jaisson (ici pour les données et quelques pages) utilisées dans leur ouvrage Le sexisme de la première heure, hasard et sociologie, qui mesure ici la fréquence de garçons à la naissance

Si l’on compare ces dernières données (via le fichier ici), calculé sur les données de l’INSEE et de l’INED, on retrouve un niveau très proche de celui que l’on a sur le fichier des prénoms, avec toutefois un biais constamment négatif.

b=read.table("http://freakonometrics.blog.free.fr/public/data/sex-ratio.txt")
X=b$V1
Y=b$V2/(100-b$V2)
lines(X,Y,col="blue")

Attention donc à la variable de sexe dans la base surtout avant guerre (en espérant que ce soit le seul soucis). Et je renvoie au blog de Baptiste qui publie toujours des choses très amusantes sur les prénoms, ici (je n’ai plus trop eu le temps de travailler dessus depuis les 5 billets en ligne ).