Tag Archives: France

Colloque, au Centre Culturel International de Cerisy, en France

Je serai pour les 6 prochains jours au Centre Culturel International de Cerisy, pour un colloque sur “l’assurance face à ses ruptures“. Au programme (13/9) histoire de l’assurance, (14/9) la financiarisation, (15/9) les données massives, (16/9) le changement climatique et (17/9) le rôle de l’État. Je ferais une intervention dimanche matin, sur le thème Certitudes collective et incertitudes individuelles, les données massives changent-elles la donne ?

Dans cette intervention, nous commencerons par un détour historique et sociologique en revenant sur les travaux d’Émile Durkheim et Max Weber. Ces chercheurs ont montré que bien que les actions individuelles puissent être imprévisibles, les comportements collectifs suivent des schémas réguliers. En agrégeant les actions individuelles, nous observons des régularités au sein de groupes plus larges, une idée illustrée par la prédictibilité du suicide dans le contexte des données massives. Ensuite, nous aborderons les réflexions de Herbert Simon et sa théorie de la rationalité limitée, suggérant que malgré les limites cognitives des individus, les rendant individuellement imprévisibles, des modèles prévisibles émergent dans la prise de décision globale. Nous explorerons ensuite des questions épistémologiques en discutant l’interprétation fréquentiste des probabilités, laquelle requiert la répétition pour être quantifiée, et la difficulté d’associer une probabilité à un événement unique. La réponse bayésienne, qui interprète les probabilités comme des croyances ou des scores, complique la relation avec l’équité actuarielle. Nous questionnerons la calibration des modèles et l’interprétation des scores à travers des exemples simples, tels que “avoir 70% de chances qu’une opération militaire réussisse” ou “avoir 70% de chances de pluie entre 14 et 15 heures”. Enfin, nous conclurons en prenant en compte le caractère temporel de la prévision: prévoir un accident 15 minutes avant qu’il ne survienne n’est pas équivalent à le prévoir un an à l’avance.

Les slides sont en ligne. Sinon Le Monde, quotidien français, publiait un portrait des colloques (organisés depuis 1952) cet été, il y a un mois,

Fairness and discrimination, PhD Course, #1 Motivation

This week, we will start our MAT998P course, in Montréal, entitled “équité et discrimination des modèles prédictifs“. It will mainly be based on the forthcoming textbook,

I can also mention the R package

> library(devtools)
> install_github("freakonometrics/InsurFair")

And because it is the first course, this week, I will start with some motivations this week… First of all, let me recall a definition, from Schauer (2006)

To be an actuary is to be a specialist in generalization, and actuaries engage in a form of decision making that is sometimes called actuarial. Actuaries guide insurance companies in making decisions about large categories that have the effect of attributing to the entire category certain characteristics that are probabilistically indicated by membership in the category, but that still may not be possessed by a particular member of the category.

Motivation #1 Redlining

In 1937, the HOLC (Home Owners’ Loan Corporation) produced the following map of Philadelphia, related to “residential security”

These maps were related to concept of “redlining”. According to Merriam Webster dictionary,

to redline is (1) to withhold home-loan funds or insurance from neighborhoods considered poor economic risks; (2) to discriminate against in housing or insurance.

On the (fictitious) maps below, we have three variables, ploted

  • some red and green areas (risky-non risky)
  • some unsanitary index (on a 0-100 scale)
  • the proportion of Black inhabitants per neiborhood

In an insurance context, risky areas (with a higher premium) should be correlated with unsanitarity index (or any risk-related variable), and those variables are legitimate predictive variables. But they can also be related to less-legitimate variable, that could be racial, here. The challenge here is that a lot of variables are correlated…

I could mention here that, for  Glenn (2000), insurer’s risk selection process has two sides:

  • the one presented to regulators and policyholders (numbers, statistics and objectivity),
  •  the other presented to underwriters (stories, character and subjective judgment).

The rhetoric of insurance exclusion – numbers, objectivity and statistics – forms what Brian Glenn calls

the myth of the actuary (…) a powerful rhetorical situation in which decisions appear to be based on objectively determined criteria when they are also largely based on subjective ones (…) or the subjective nature of a seemingly objective process.

Glenn  (2003) claimed that there are many ways to rate accurately. Insurers can rate risks in many different ways depending on the stories they tell on which characteristics are important and which are not.

The fact that the selection of risk factors is subjective and contingent upon narratives of risk and responsibility has in the past played a far larger role than whether or not someone with a wood stove is charged higher premiums (…) virtually every aspect of the insurance industry is predicated on stories first and then numbers

Motivation #2. “Gender directive”, 2004/113/EC

From the Treaty on European Union (26.10.2012)

Art. 2 The Union is founded on the values of respect for human dignity, freedom, democracy, equality, the rule of law and respect for human rights, including the rights of persons belonging to minorities. These values are common to the Member States in a society in which pluralism, non-discrimination, tolerance, justice, solidarity and equality between women and men prevail.

Art. 3 (…) It shall combat social exclusion and discrimination, and shall promote social justice and protection, equality between women and men, solidarity between generations and protection of the rights of the child.

from the Charter of Fundamental Rights of the European Union (18.12.2000)

Art. 21 (Non discrimination): Any discrimination based on any ground such as sex, race, colour, ethnic or social origin, genetic features, language, religion or belief, political or any other opinion, membership of a national minority, property, birth, disability, age or sexual orientation shall be prohibited.

Art. 23 (Equality between men and women) Equality between men and women must be ensured in all areas, including employment, work and pay. The principle of equality shall not prevent the maintenance or adoption of measures providing for specific advantages in favour of the under-represented sex.

and from the EU Directive 2004/113/EC (2004 version)

Art. 5 (Actuarial factors)

1. Member States shall ensure that in all new contracts concluded after 21 December 2007 at the latest, the use of sex as a factor in the calculation of premiums and benefits for the purposes of insurance and related financial services shall not result in differences in individuals’ premiums and benefits.

2. Notwithstanding paragraph 1, Member States may decide before 21 December 2007 to permit proportionate differences in individuals’ premiums and benefits where the use of sex is a determining factor in the assessment of risk based on relevant and accurate actuarial and statistical data. The Member States concerned shall inform the Commission and ensure that accurate data relevant to the use of sex as a determining actuarial factor are compiled, published and regularly updated.

There was initially (2004) an opt-out clause (Article 5(2)), since, where gender is a determining factor in the assessment of
risk based on relevant and accurate actuarial and statistical
data then proportionate differences in individual premiums or
benefits are allowed.

In March 2011, the European Court of Justice issued its judgement into the “Test-Achats case”. The ECJ ruled Article 5(2) was invalid. Thus, insurers were no longer able to use gender as a risk factor when pricing policies.

Other legal documents in Europe can be mentioned such that the “Ten Oever” judgement (Gerardus Cornelis Ten Oever v Stichting Bedrijfspensioenfonds voor het Glazenwassers — en Schoonmaakbedrijf). In April 1993, the Advocate General Vangerven argued that (see De Baere (2012))

the fact that women generally live longer than men has no significance at all for the life expectancy of a specific individual and it is not acceptable for an individual to be penalized on account of assumptions which are not certain to be true in his specific case,

which could be related to the concept of “injustice by generalization”.

Motivation #3. Colorado (September 27, 2023)

In September 27, 2023, the Colorado Division of Insurance exposed a new proposed regulation entitled Concerning Quantitative Testing of External Consumer Data and Information Sources, Algorithms, and Predictive Models Used for Life Insurance Underwriting for Unfairly Discriminatory Outcomes.

Section 4 (Definitions) Bayesian Improved First Name Surname Geocoding, or “BIFSG” means, for the purposes of this regulation, the statistical methodology developed by the RAND corporation for estimating race and ethnicity.

External Consumer Data and Information Source, or “ECDIS” means, for the purposes of this regulation, a data source or an information source that is used by a life insurer to supplement or supplant traditional underwriting factors. This term includes credit scores, credit history, social media habits, purchasing habits, home ownership, educational attainment, licensures, civil judgments, court records, occupation that does not have a direct relationship to mortality, morbidity or longevity risk, consumer-generated Internet of Things data, biometric data, and any insurance risk scores derived by the insurer or third-party from the above listed or similar data and/or information source.

Then we have different sections, where insurers are asked to “estimate” the race or ethnicity of policyholders

Section 5 (Estimating Race and Ethnicity) : Insurers shall estimate the race or ethnicity of all proposed insureds that have applied for coverage on or after the insurer’s initial adoption of the use of ECDIS, or algorithms and predictive models that use ECDIS, including a third party acting on behalf of the insurer that used ECDIS, or algorithms and predictive models that used ECDIS, in the underwriting decision-making process, by utilizing:

1. BIFSG and the insureds’ or proposed insureds’ name and geolocation information included in the applications) for life insurance shall be used to estimate the race and ethnicity of each insured or proposed insured.

2. For the purposes of BIFSG, the following racial and ethnic categories shall be used: Hispanic, Black, Asian Pacific Islander (API), and White.

Section 6 (Application Approval Decision Testing Requirements) : Using the BIFSG estimated race and ethnicity of proposed insureds and the following methodology, insurers shall calculate whether Hispanic, Black, and API proposed insureds are disapproved at a statistically significant different rate relative to White applicants for whom the insurer, or a third party acting on behalf of the insurer, used ECDIS, or an algorithm or predictive model that used ECDIS, in the underwriting decision-making process.

1. Logistic regression shall be used to model the binary underwriting outcome of either approved or denied.

2. The following factors may be accounted for as control variables in the regression model: policy type, face amount, age, gender, and tobacco use.

3. The estimated race or ethnicity of the proposed insureds shall be accounted for by including Hispanic, Black, and Asian Pacific Islander (API) as separate dummy variables in the regression model.

4. Determine if there is a statistically significant difference in approval rates for each BIFSG estimated race or ethnicity variable as indicated by a p-value of less than .05.

a. If there is not a statistically significant difference in approval rates, no further testing is required.

b. If there is a statistically significant difference in approval rates, the insurer shall determine whether the difference in approval rates is five (5) percentage points or greater as indicated by the marginal effects value of each BIFSG estimated race or ethnicity variable. (…)

or

Section 7 (Premium Rate Testing Requirements) : Using the insureds’ BIFSG estimated race and ethnicity, insurers shall determine if there is a statistically significant difference in the premium rate per $1,000 of face amount for policies issued to Hispanic, Black, and API insureds relative to White insureds for whom the insurer, or a third party acting on behalf of the insurer, used ECDIS, or an algorithm or predictive model that used ECDIS, in the underwriting decision-making process.

1. Linear regression shall be used to model the continuous numerical outcome of premium rate per $1,000 of face amount.

2. The following factors may be accounted for as control variables in the regression model: policy type, face amount, age, gender, and tobacco use.

3. The estimated race or ethnicity of the proposed insureds shall be accounted for by including Hispanic, Black, and Asian Pacific  Islander (API) as separate dummy variables in the regression model.

4. Determine if there is a statistically significant difference in the premium rate per $1,000 of face amount for each BIFSG estimated race or ethnicity variable as indicated by a p-value of less than .05.

a. If there is not a statistically significant difference in premium rate per $1,000 of face amount, no further testing is required.

b. If there is a statistically significant difference in premium rate per $1,000 of face amount, determine whether the premium rate per $1,000 of face amount is at least 5% more than the average premium rate per $1,000 for all policies.

i. If the difference in premium rate per $1,000 of face amount is less than 5%, no further testing is required.

ii. If the difference in premium rate per $1,000 of face amount is 5% or greater, further testing is required as described in Section 8.

(etc). In order to illustrate, we can use some data, in the region of Atlanta

 

We can change the first and last name of people (and keep other relevant information, including the ZIP code) and compare “predictions” of race (white, black, hispanic, asian, etc)

Motivation #4. Motor Insurance in the U.S.

In the context of motor insurance in the U.S., recall that legal restrictions are per states, and we can observe some diversity about what “sensitive” could mean (via thezebra)

(etc). We will also discuss Avraham et al. (2013) that provides a long discussion accross US states.

Motivation #5. Graduate Admission (UC Berkeley)

Another motivation is the popular article, Bickel,  Hammel, and O’Connell (1975)

The dataset mentioned in the article is the following

the bias in the aggregated data stems not from any pattern of discrimination on the part of admissions committees, which seems quite fair on the whole, but apparently from prior screening at earlier levels of the educational system. Women are shunted by their socialization and education toward fields of graduate study that are generally more crowded, less productive of completed degrees, and less well funded, and that frequently offer poorer professional employment prospects

As we can see, if we formalize, we have (almost)

This is Simpson’s paradox. Another simple example is related to mortality : the (overall) mortality rate for women (picked at random in the entiere population) was 0.812% in Costa Rica, lower than 0.929% in Sweden. But as we can see on the left, below, at any age, mortality rates are lower in Sweden than in Costa Rica.

The paradox can easily be explained if we look at age structures in both countries. Long story short, in Costa Rica, picking someone randomly means that the person is very likely to be (very) young, with a low mortality rate; in Sweden, the person is more likely to be older, with a higher mortality rate.

Motivation #6. Propublica, Actuarial Justice

We will also mention actuarial justice, and et al (2016)

Hence, looking at the same data, with difference perspective, could lead to different conclusions. More robust conclusions can be obtained when look at distributions of scores (instead of simple binary predictions)

and we can also consider temporal process (again, instead of simply binary variables, with temporal censoring)

Motivation #7. Insurance in Québec

Two final motivations, in French this time. In Québec, there is the Charte des droits et libertés de la personne (C-12) with some very clear definition of what “discrimination” means,

Art. 10  Toute personne a droit à la reconnaissance et à l’exercice, en pleine égalité, des droits et libertés de la personne, sans distinction, exclusion ou préférence fondée sur la race, la couleur, le sexe, l’identité ou l’expression de genre, la grossesse, l’orientation sexuelle, l’état civil, l’âge sauf dans la mesure prévue par la loi, la religion, les convictions politiques, la langue, l’origine ethnique ou nationale, la condition sociale, le handicap ou l’utilisation d’un moyen pour pallier ce handicap.

Il y a discrimination lorsqu’une telle distinction, exclusion ou préférence a pour effet de détruire ou de compromettre ce droit.

But, interestingly, insurers can almost do anything they want,

Art 20.1 Dans un contrat d’assurance ou de rente, un régime d’avantages sociaux, de retraite, de rentes ou d’assurance ou un régime universel de rentes ou d’assurance, une distinction, exclusion ou préférence fondée sur l’âge, le sexe ou l’état civil est réputée non discriminatoire lorsque son utilisation est légitime et que le motif qui la fonde constitue un facteur de détermination de risque, basé sur des données actuarielles.

Motivation #8. Intention

And finally, I can mention that in many countries (such as France), “indirect discrimination” is considered as discriminatory, so “intention” has nothing to do with the problem… The Loi no 2008-496 du 27 mai 2008 states that

Art. 1 Constitue une discrimination indirecte une disposition, un critère ou une pratique neutre en apparence, mais susceptible d’entraîner, pour l’un des motifs mentionnés au premier alinéa, un désavantage particulier pour des personnes par rapport à d’autres personnes, à moins que cette disposition, ce critère ou cette pratique ne soit objectivement justifié par un but légitime et que les moyens pour réaliser ce but ne soient nécessaires et appropriés.

This law is an extension of Loi no. 72-546 du 1er juillet 1972, which abolished the requirement for specific intent.

Again, following Avraham (2017), keep in mind that insurance is very specific, regarding discrimination

What is unique about insurance is that even statistical discrimination which by definition is absent of any malicious intentions, poses significant moral and legal challenges. Why? Because on the one hand, policy makers would like insurers to treat their insureds equally, without discriminating based on race, gender, age, or other characteristics, even if it makes statistical sense to discriminate (…) On the other hand, at the core of insurance business lies discrimination between risky and non-risky insureds. But riskiness often statistically correlates with the same characteristics policy makers would like to prohibit insurers from taking into account.

That will be the topic of the course…

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,