Big Data : la recherche s’expose

Vendredi matin, j’interviendrai à la journée “Big Data : la recherche s’expose” organisée par l’Institut Louis Bachelier au petit auditorium du palais Brongniart. Au programme

  • 9h00 : Accueil par Bertrand Villeneuve, Université Paris-Dauphine
  • 9h15-9h45 : René Aïd, IdR « Laboratoire Finance des Marchés de l’Energie », Fondation Institut Europlace de Finance.
  • 9h45 – 10h15 : Florence Jusot, Chaire « Santé », Fondation du Risque.
  • 10h15 – 10h45 : Pause-café
  • 10h45 – 11h15 : Emmanuel Bacry, Chaire « Marchés en Mutation », Fondation Institut Europlace de Finance.
  • 11h15 – 11h45 : Arthur Charpentier, IdR « Valorisation et Nouveaux Usages Actuariels de l’Information », Fondation du Risque.
  • 11h45 – 12h15 : Serge Darolles, IdR « Private Equity et Venture Capital », Fondation Institut Europlace de Finance.
  • 12h15 : Conclusion par Nizar Touzi, Ecole Polytechnique.
  • 12h30 : Buffet déjeuner
  • 14h00 : Fin de la Journée des Chaires 2016.

Je présenterais les résultats obtenus lors du pricing game obtenu cet automne (qui sera relancé, sous une forme actualisé, cet été.. à suivre donc).

Ce que la courbe ROC (et l’AUC) ne raconte pas

En préparant une intervention pour mardi prochain, j’épluchais les résultats renvoyés pour un exercice, et j’ai eu un résultat assez étrange avec un modèle de classification. J’avais donné la même base cet automne à l’ensae, et j’avais donc près d’une trentaine d’autres modèles, pour comparer (disons plutôt que sur la même base de test, j’ai une trentaine de prévisions). Les observations noires sont celles obtenues cet automne (le trait correspond aux meilleurs AUC sur la base de test), et les observations rouges sont celles que j’ai obtenu pour l’intervention de mardi (là encore, le trait vertical correspond aux meilleurs modèles, au sens du AUC), sur une observation de la base de test,

Ce sont les probabilités prédites (mais j’ai enlevé l’échelle).

Pour presque toutes mes observations, les poids rouges sont bien au dessus des autres… Mais cela n’enlève en rien le fait que l’AUC obtenu (pour les deux modèles rouges) est très bon. C’est effectivement un résultat important (et connu) : le critère AUC (mais plus généralement la courbe ROC en entier) n’indique en rien si la valeur prédite est bonne, ou pas. Il nous dit juste si l’ordre obtenu est correct. Si les valeurs les plus importantes sont effectivement les valeurs pour lesquelles on a un 1, l’AUC sera très bon.

C’est ce qu’on peut observer sur le petit exemple ci-dessous. Considérons un modèle logistique simulé assez simple,

> n=1e3
> set.seed(1)
> x1=rnorm(n)
> x2=runif(n)
> u=-3+x2+x1
> p=exp(u)/(1+exp(u))
> y=rbinom(n,prob=p,size=1)
> library(ROCR)
> df=data.frame(y,x1,x2)
> mean(df$y)
[1] 0.116
> reg=glm(y~.,data=df,family=binomial)
> p=predict(reg,type="response")
> mean(p)
[1] 0.116
> pred1=prediction(p, df$y)
> L=performance(pred1, "tpr", "fpr")

L’AUC est ici

> auc=performance(pred1, "auc")@y.values[[1]]
> auc
[1] 0.7681191

et la courbe ROC est la suivante,

> plot(unlist(L@x.values),unlist(L@y.values),
+ type="s",col="blue")

Supposons que l’on change la constante du modèle logistique,

> reg$coefficients[1]=0

Dans ce cas, notre prévision est assez mauvaise, car la probabilité moyenne prédite est ici

> u=reg$coefficients[1]+reg$coefficients[2]*
+ df$x1+reg$coefficients[3]*df$x2
> p=exp(u)/(1+exp(u))
> mean(p)
[1] 0.6060676

(on est loin des 11,6% de 1 dans la base). Pourtant le AUC est bon

> pred1=prediction(p, df$y)
> L=performance(pred1, "tpr", "fpr")
> auc=performance(pred1, "auc")@y.values[[1]]
> auc
[1] 0.7681191

(c’est en fait la même valeur qu’auparavant, ce qui a du sens puisque la courbe ROC est identique)

> lines(unlist(L@x.values),unlist(L@y.values),
+ type="s",col="red")

Autrement dit, ces outils, classiquement utilisés pour juger la qualité d’un classifieur, ne permettent en aucun cas de dire que la probabilité prédite a du sens. Ces critères permettent juste de dire qu’on identifie assez bien les personnes qui ont le plus de chance d’avoir la réponse 1. Ce qui n’est pas si mal… mais c’est un autre problème que celui d’avoir une probabilité qui soit pertinente.

Popularité, ou reconnaissance ?

Assez souvent, quand on est enseignant/chercheur, on remplit des demandes de financements, ou de promotions, on nous demande de parler de l’impact de nos travaux et de notre reconnaissance. La première fois que j’ai fait face à cette question, c’était il y a 5 ans, à Montréal, pour une demande de financement auprès du NSERC, au niveau fédéral. Sur le moment, la question m’a déstabilisé. Non pas qu’elle n’ait pas de sens, mais qu’on me la pose à moi. Qui suis-je pour savoir si ce que je fais a le moindre impact ?

Classiquement, l’impact, on le mesure avec les citations par exemple. On peut alors utiliser l’indice synthétique de Hirsch. Ces citations, elles commencent à figurer sur nos CV, ou tout du moins sur les CV en ligne que sont google scholarorcidresearch gateacademia,  etc. Un autre endroit où on retrouve la reconnaissance, et c’est probablement lié aux citations, ce sont les invitations en séminaire ou en conférence, les demandes pour faire partie de jurys, etc. Je devrais préciser que vu de l’extérieur, ce qu’on voit, ce sont les invitations “acceptées”. Les invitations en conférence sont probablement une information importante, mais là par exemple – pour des raisons personnelles – je me retrouve à décliner des invitations. Mon activité va baisser, mais comment l’interpréter ? Dans la liste, on devrait peut être rajouter les invitations pour arbitrer des publications d’articles dans des journaux et revues. On est alors très clairement identifié comme un pair sur un thème précis.

Mais pour comprendre en quoi nous avons un quelconque impact, il est toujours agréable quand on nous informe sur la raison qui a motivé une invitation. Malgré tout, c’est assez rare. Je pourrais toutefois noter quelques cas, où on nous invite pour présenter un papier particulier dans une conférence, ou quand l’éditeur explique brièvement la raison de la sollicitation, lors d’une demande d’arbitrage. Un exemple récent que je pourrais mentionner (histoire d’illustrer un peu mon billet) est l’invitation que j’avais reçu, en début d’année, pour aller parler à la Bank of England,

Plus précisément, comme le dit l’invitation (je reprends)

The Chief Economists Workshop is an annual forum that we host to facilitate high-level discussion between the Chief Economists from around thirty-five of the major central banks in the world, and academic and professional experts. The format is a mixture of presentations, discussion and debate. As such it offers a great opportunity for you to interact with and influence leading policymakers. There is no media presence. I attach a copy of the programmes from the previous two years to give you an idea of the format.

Many of us in the Bank of England are avid followers of your Twitter feed, and your ability to offer multi-disciplinary insights to current policy problems. I would like to invite you to give a talk along these lines – for example on actuarial and statistical insights into Big Data. It is not necessary to have a formal paper. The idea is more to discuss the issues and generate debate. We would be looking for a presentation of about 45 minutes followed by 45 minutes of Q and A and can schedule your talk to fit in with your commitments.

J’ai eu le plaisir d’être invité avec Andrew Gelman (NYU) et Andrew Lo (MIT). J’ai l’impression qu’un universitaire devrait être humilié de recevoir une telle lettre d’invitation (même si mon compte twitter est l’occasion de partager des informations en lien avec les thèmes de recherche qui m’intéressent – la statistique, l’économie, l’assurance, les risques, le climat, etc). Mais si on y pense deux minutes, la raison invoquée a l’avantage d’être honnête, et finalement flatteuse (on m’invite pour mon côté polymathe, ce qui me fait infiniment plaisir). Je doute que beaucoup de monde lise mes papiers (vu le prix qu’il faut payer pour y accéder, ça se comprendrait presque), et j’ai l’impression que quand on me demande de venir parler dans une conférence, c’est pour partager des connaissances, donner mon point de vue sur tel ou tel point, sur lequel j’ai pu réfléchir au cours de mes travaux de recherche.

Plus récemment, l’autre jour, j’étais particulièrement fier d’être dans les 10 scientifiques (médecins et scientifiques au sens large) les plus populaires sur Twitter, selon Science & Avenir,

D’autant plus que je n’avais rien demandé. Il y a plusieurs mois, j’avais dit – peut être un peu maladroitement – que je ne participerais plus à ces concours, à ces classements (en l’occurrence on nous invitait à participer à un répertoire des meilleurs blogs de sciences). Mais c’est toujours valorisant d’apprendre qu’on est “populaire”. Comme tous les chercheurs, on marche à la carotte, et avoir une médaille est une forme de reconnaissance qui fait forcément plaisir.

Ce n’est certes pas une reconnaissance par mes pairs, mais c’est malgré tout une forme de reconnaissance. Qui permet d’être interrogé (et cité) dans une enquête qui revient sur le financement de la couverture des catastrophes naturelles en France, une chronique sur les cat bonds, ou sur les investissements en un super-calculateur sur le climat. Cette reconnaissance en dehors du cercle classique (et très restreint) d’experts de sa propre discipline a du bon. Pour élargir son champ de vision tout d’abord et avoir un point de vue très différent. Mais aussi parfois rencontrer des gens en dehors du monde académique. Qui ne s’embarrassent pas à payer des éditeurs scientifiques pour avoir accès aux articles, mais qui sont d’avides lecteurs des blogs, qui ont souvent des données, ou des problèmes intéressants à nous soumettre. En revenant sur le financement du régime cat’nat en France pour les besoins de l’entrevue, on a commencé avec des collègues à retravailler sur la prévention. Personnellement, ce sont toujours ces discussions qui me permettent de trouver des thèmes de recherche (en plus de la lecture d’articles académiques finissant par un “tiens, il devrait être possible de regarder…” ou un “personnellement, j’aurais plus fait…”). Et j’adore ça ! peut être parce que c’est toujours l’occasion d’apprendre de nouvelles choses.

Optimal Claiming Strategies in Bonus Malus Systems and Implied Markov Chains

With Arthur David and Romuald Elie, we just wrote a short paper on bonus malus, and optimal strategies to claim a loss (or not)

In this paper, we investigate the impact of the claim reporting strategy of drivers, within a bonus malus system. We exhibit the induced modification of the corresponding class level transition matrix and derive the optimal reporting strategy for rational drivers. The hunger for bonuses induces optimal thresholds under which, drivers do not claim their losses. A numerical algorithm is provided for computing such thresholds and realistic numerical applications are discussed.

The paper is now online on http://papers.ssrn.com/id=2790583 and https://hal.archives-ouvertes.fr/hal-01326798.

Note that we do not discuss here legal issues here (in some contracts, it is compulsory to claim all losses, even small ones), but economic incentives and mathematical issues. Some popular journals in France did mention that issue, of non claims small losses (see http://leparticulier.fr/) but in those standard computations (see below), it is based on some naive model that we improve in our paper,

Où ont été observées les fortes hausses du prix de l’essence ?

Vendredi, dans mon billet “le prix de l’essence en cas de pénurie“, j’essayais de visualiser l’explosion du prix de l’essence dans certaines stations services, ces dernières semaines. Mais plus généralement, on peut se demander où ont été observées les fortes hausses.

A partir de la base téléchargeable sur http://www.prix-carburants.economie.gouv.fr/, il est possible de récupérer les informations les prix, par jour, dans toutes les stations essence de France. Pour ça, on utilise la fonction suivante,

> 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,
+ cp=l[[no]]$.attrs["cp"])
+ base=rbind(base,base1)
+ }}}
+ return(base)}

Par exemple, on peut récupérer les prix fin mai

> B1=spatial(as.Date("2016-05-31"))

et début mai,

> B2=spatial(as.Date("2016-05-01"))

En fusionnant ces deux bases, on peut récupérer la variation du prix, sur un mois,

> names(B1)=c("indice","lat","long","fin","cp")
> names(B2)=c("indice","lat","long","debut","cp")
> B=merge(B1,B2)
> B$var=(B$fin-B$debut)/B$debut*100

Plus particulièrement, on va s’intéresser aux stations dans lesquelles le prix a augmenté de plus de 10%,

> idx=which((B$lon>(-10))&(B$lon<20)&
+             (B$lat>35)&(B$lat<55))
> B=B[idx,]
> B=B[!is.na(B$var),]
> B$Y=(B$var>10)*1

Notons déjà que dans plus de 90% des stations, le prix a augmenté entre début mai, et fin mai

> mean(B$var > 0)
[1] 0.9288522

mais les hausses de plus de 10% restent rares (moins de 5% des stations essence)

> mean(B$var > 10)
[1] 0.03305452

La distribution de la variation des prix (en %) est la suivante

> plot(density(B$var),xlim=c(-5,20))

Si on regarde maintenant où ont été observées ces fortes hausses, on observe qu’elles peuvent avoir lieu n’importe où

> library(maps)
> map("france")
> points(B$lon,B$lat,pch=19,col=c(rgb(0,0,1,.25),
+ rgb(1,0,0,.95))[1+(B$var>10*1)],cex=.5)

et c’est encore plus flagrant quand on regarde les hausses de plus de 5%

> points(B$lon,B$lat,pch=19,col=c(rgb(0,0,1,.35),
+ rgb(1,0,0,.95))[1+(B$var>5*1)],cex=.5)

Et si on regarde les endroits où la proportion des stations ayant fortement augmenté les prix est la plus forte, on obtient la carte suivante,

Autrement dit la hausse des prix a pu être observée un peu partout en France. On peut alors comparer avec la carte des pénuries, mise en ligne sur https://mon-essence.fr/

Si certaines régions avec des fortes hausses de prix correspondent à des pénuries d’essence, ce n’est pas le cas partout. Bref, comprendre ce qui s’est passé sur les prix ce mois de mai est loin d’être simple… à suivre donc….

Le prix de l’essence en cas de pénurie

Cet après midi, sur le site de 60 millions de consommateurs, Lionel Maugain revenait sur les stations services qui avaient de très fortes variations de prix pendant les jours de pénurie.

En utilisant la base des prix des carburants à la pompe, téléchargeable sur http://www.prix-carburants.economie.gouv.fr/, il est possible de visualiser ces fortes variations de prix. Pour lire la base, rien de plus simple

> rm(list=ls())
> annee=year=2016
> fichier="PrixCarburants_annuel_2016.xml"
> library(plyr)
> library(XML)
> library(lubridate)
> l=xmlToList(fichier)

Ensuite, on transforme la grosse liste en une série temporelle, à l’aide de la fonction suivante

> 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,"-06-01 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)
+      vp=P[1+vect_idx]
+      vp=vp*100/(vp[!is.na(vp)][1])
+      base2=ts(vp,start=year,frequency=365)
+      list(base=base1,ts=base2)
+    }

Dans l’article, la plus forte variation a été observée à Sevran. Pour récupérer la liste des stations services, on utilise

> ville=list()
> for(no in 1:length(l)) 
      ville[[no]]=l[[no]]$ville
> ville=unlist(ville)
> idx=which(ville=="Sevran")

Pour la première station, l’évolution du prix est relativement modérée

> plot(time_series(idx[1])$ts,col="red",type="l")
> abline(h=100,lty=2)

Même si le prix semble avoir augmenté de 5% pendant les dernières semaines. Mais ce n’est effectivement rien en comparaison de l’autre station,

> plot(time_series(idx[2])$ts,col="red",type="l")
> abline(h=100,lty=2)

avec un prix qui a augmenté de 50% en 2 mois, et presque 30% en quelques jours.

A Villeneuve-d’Asq, on retrouve un comportement similaire, mais moins important

> idx=which(ville=="Villeneuve-d'Ascq")
> plot(time_series(idx[3])$ts,col="red",type="l")
> abline(h=100,lty=2)