Tomorrow, François will present our joint work at the 17th Financial Risks International Forum “Big Data & Algorithmic Finance”, A Sequentially Fair Mechanism for Multiple Sensitive Attributes.
His slides are now available.
Tag Archives: Paris
Talk at the Groupe de travail ARC (Actuariat et Risques Contemporains), in Paris
This Friday, I will give a talk in Paris, on using optimal transport to mitigate unfair predictions, at the ARC Seminar.
The insurance industry is heavily reliant on predictions of risks based on characteristics of potential customers. Although the use of said models is common, researchers have long pointed out that such practices perpetuate discrimination based on sensitive features such as gender or race. Given that such discrimination can often be attributed to historical data biases, an elimination or at least mitigation is desirable. With the shift from more traditional models to machine-learning based predictions, calls for greater mitigation have grown anew, as simply excluding sensitive variables in the pricing process can be shown to be ineffective. In this talk, we first investigate why predictions are a necessity within the industry and why correcting biases is not as straightforward as simply identifying a sensitive variable. We then propose to ease the biases through the use of Wasserstein barycenters instead of simple scaling. To demonstrate the effects and effectiveness of the approach we employ it on real data and discuss its implications. The talk will be based on recent work with François Hu and Philipp Ratz (2310.20508, 2309.06627, 2306.12912 and 2306.10155).
Slides are now available online.
Equité et discrimination en assurance
Mercredi 18 mai, j’interviendrai au séminaire de la chaire PARI, à Paris, dans les locaux de la CCR, pour parler d’équité et de discrimination en assurance… Les slides sont en ligne.
En poursuivant la discussion abordée dans “L’équité de l’apprentissage machine en assurance” (co-écrit avec Laurence Barry), nous présentons un aperçu des problèmes auxquels les actuaires sont confrontés lorsqu’ils traitent de la discrimination. Depuis le début de leur histoire, les assureurs sont connus pour utiliser des données pour classifier et tarifer les risques. En tant que tels, ils ont été confrontés très tôt au problème de l’équité et de la discrimination associées aux données. Cette question devient de plus en plus importante avec l’accès à des données plus granulaires et comportementales, et évolue pour refléter les technologies et les préoccupations sociétales actuelles. En examinant les débats antérieurs sur la discrimination, nous montrons que certains préjugés algorithmiques sont une version renouvelée de préjugés plus anciens, tandis que d’autres semblent inverser l’ordre précédent. Paradoxalement, alors que la pratique de l’assurance n’a pas profondément changé et que la plupart de ces biais ne sont pas nouveaux, l’ère de l’apprentissage automatique ébranle encore profondément la conception de l’équité en matière d’assurance. (Cette présentation s’appuiera également sur un rapport qui sera bientôt publié par l’Institut Louis Bachelier).
Changement Climatique et Assurance
Mercredi, je participerais à la Journée du Marché organisée par AON, pour la table ronde du matin sur le thème du changement climatique, à Paris.
La première partie de la table-ronde portera sur l’impact du changement climatique sur l’assurance et la réassurance, et la partie sur l’impact de l’assurance sur le changement climatique. J’interviendrais principalement dans la première partie.
J’en profite pour remettre un lien vers un court article que nous avions publié il y a quelques temps, ‘Changement Climatique et Assurance‘, dans Variances.
Antisélection en Assurance Automobile et Impact des Jours de Carence sur les Arrêts Maladies
Cet après midi, je suis à Paris-Jourdan pour la soutenance de thèse d’Alexandre Godzinski , sur Three empirical essays on moral hazard identification in insurance, supervisée par Pierre-Yves Geoffard.
Pour l’anecdote, le premier chapitre (Effects of the one-day waiting period for sick leave on health-related absences in the French central civil service) pour sur le jour de carence pour l’assurance maladie (sujet qui fait couler de l’encre ces jours-ci).
Ce jour de carence a été introduit au 1er janvier 2012 et supprimé
au 1er janvier 2014. Durant ces deux années d’application, le premier jour d’arrêts maladie était non rémunéré. Il y avait un groupe de contrôle (le secteur privé) et un groupe test (le secteur public). Le jour de carence conduit à une baisse plus de 50 % de la prévalence des absences de 2 jours,
mais aussi une hausse de 25 % de la prévalence des absences de 1 semaine à 3 mois
Aucun effet n’est trouvé sur la prévalence totale pour la population globale. Les agents ne sont donc ni plus ni moins absents pour raison de santé…
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. @BrunoFERRY1 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…
Crowded Cities, Paris, Hong Kong and Montréal
Over the past years, I’ve been living in different cities, all of them being completely different, compared with the others. I have been living in Paris, which is a big city in Europe, with a large suburban neighborhood, too (la banlieue).
Then, I’ve been living in Hong Kong, which is a larger city, in Asia.
It was crowded. I mean, it was the feeling I had, while I was living there. And more recently, I’ve been living in Montréal, in North America. Montreal is a large city. Or to be more specific, an island,
The three cities are quite different. Paris, 2.211 million unhabitants, and 105,4 km² (density 21,057 unhabitants per km²). Montréal, 1.621 million unhabitants, and three times wider 365.1 km² (density 4,441 unhabitants per km²). Hong Kong, 7.234 million unhabitants, and again three times wider 1,104 km² (density 6,553 unhabitants per km²). In Hong Kong, there are several hill where it is not possible to build anything: on a large part of the island, the density is null.
Continue reading Crowded Cities, Paris, Hong Kong and Montréal
Short (but Buzy) Trip in Europe
I will flight tomorrow evening to go to Belgium. I will attend a PhD defense this Monday, at Louvain-la-Neuve. On Tuesday, I will give a crash course on risk measures, for PhD Students, in Louvain-la-Neuve. Then, Wednesday morning, I have a meeting with some friends in an insurance company, in Brussels, and I will be around noon in Paris, to give a talk on Big Data, for the Big Data Task Force of the (French) Institute of Actuaries. I will then have some meetings in the afternoon. I should flight back on Thursday, and give my lecture on Predictive Modeling in the evening…
I guess that Saturday will be a mix between two very different activities, enjoying some time with the kids (and watch the game), and trying to finalize my slides…
Le Vote par Procuration en France
La Vie des Idées a mis en ligne, ce matin, un court texte, écrit par Baptiste Coulmont (a.k.a. @coulmont) et Joël Gombin (a.k.a. @joelgombin), auquel j’ai très modestement contribué, intitulé “Un homme, deux voix. Le vote par procuration“.
Alors que sur son blog, Baptiste a rajouté pas mal d’information sur le vote par procuration en France (et le contexte général, en particulier pourquoi autant de partis courtisent certaines personnes en les incitant à voter par procuration), et sur les bases de données, je voulais en profiter pour mettre en ligne quelques codes utilisés dans l’article, et en particulier, mentionner des graphiques non-utilisés car plus difficile à interpréter, mais à mon avis plus juste en terme de modèle (comme les conclusions étaient les mêmes, on a retenu des graphiques plus classiques). Rappelons tout d’abord qu’on analyse non pas le vote à partir de données individuelles (ceci ne peut s’obtenir, le vote étant encore secret en France), mais à partir des résultats des différents bureaux de vote (c’est la notion de corrélation écologique évoquée dans le texte, à cause du problème potentiel d’ecological fallacy). Moyennant toutes ces précautions d’usage, on a essayé d’analyser les données à notre disposition.
- Modèle de régression, et recherche de variables explications
Pour faire une régression, et expliquer le taux de procurations dans un bureau de vote, ma première idée était de dire que où
est le nombre de procurations dans le bureau de vote
, et où
est (au choix) le nombre d’électeurs inscrits ou le nombre d’électeurs ayant voté. On suppose ici que
, la proportion d’électeurs qui a voté par procuration peut être fonction de divers variables explicatives. Les variables, elles sont dans la base suivante (je renvoie vers le blog de Baptiste pour les bases que l’on utilise, en particulier à partir des données de insee.fr et d’opendata.paris.fr)
> bt1=read.table("paris2007-pres-t1.csv",header=TRUE,sep=";")
> bt2=read.table("paris2007-pres-t2.csv",header=TRUE,sep=";")
> bv=read.table("paris-bv-insee-07.csv",header=TRUE,sep=";")
> bv$BV=bv$BVCOM
> baset1=merge(bt1,bv,by="BV")
> baset2=merge(bt2,bv,by="BV")
> baset1$LOGEMENT=baset1$PROPRIO+baset1$LOCNONHLM+baset1$LOCHLM+baset1$GRATUIT
> baset2$LOGEMENT=baset2$PROPRIO+baset2$LOCNONHLM+baset2$LOCHLM+baset2$GRATUIT
Si on suppose que est fonction
le taux de logements occupés par leur propriétaire, dans le quartier (associé à un bureau de vote),
> variable="PROPRIO" > reference="LOGEMENT" > baset1$taux=baset1[,variable]/baset1[,reference] > baset2$taux=baset2[,variable]/baset2[,reference]
il est légitime de tenter une régression logistique,
voire un lissage par splines, si on pense que le lien peut ne pas être linéaire,
Ceci se fait à l’aide du code suivant, pour la version lissée (par splines cubiques)
> b=hist(baset1$taux,plot=FALSE)
> library(splines)
> regt1=glm(PROCURATIONS/INSCRITS~bs(taux,6),family=binomial,weights=INSCRITS,data=baset1)
> regt2=glm(PROCURATIONS/INSCRITS~bs(taux,6),family=binomial,weights=INSCRITS,data=baset2)
> u=seq(min(baset1$taux)+.015,max(baset1$taux)-.015,by=.001)
> ND=data.frame(taux=u)
> ug=seq(0,max(baset1$taux)+.05,by=.001)
> pt1=predict(regt1,newdata=ND,se=TRUE,type="response")
> pt2=predict(regt2,newdata=ND,se=TRUE,type="response")
> library(RColorBrewer)
> CL=brewer.pal(6, "RdBu")
> plot(ug,ug*1,col="white",xlab=nom,ylab="Taux de procuration",
+ ylim=c(0,.1))
> for(i in 1:(length(b$breaks)-1)){
+ polygon(b$breaks[i+c(0,0,1,1)],c(0,b$counts[i],b$counts[i],0)
+ /max(b$counts)*.05,col="light yellow",border=NA)}
> polygon(c(u,rev(u)),c(pt1$fit+2*pt1$se.fit,rev(pt1$fit-2*pt1$se.fit)),
+ border=NA,density=30,col=CL[4])
et, pour la régression logistique standard (linéaire)
> lines(u,pt1$fit,col=CL[6],lwd=2)
> polygon(c(u,rev(u)),c(pt2$fit+2*pt2$se.fit,rev(pt2$fit-2*pt2$se.fit)),
+ border=NA,density=30,col=CL[3])
> lines(u,pt2$fit,col=CL[1],lwd=2)
> regt1l=glm(PROCURATIONS/INSCRITS~taux,family=binomial,weights=INSCRITS,data=baset1)
> regt2l=glm(PROCURATIONS/INSCRITS~taux,family=binomial,weights=INSCRITS,data=baset2)
> ND=data.frame(taux=ug)
> pt1l=predict(regt1l,newdata=ND,se=TRUE,type="response")
> pt2l=predict(regt2l,newdata=ND,se=TRUE,type="response")
> lines(ug,pt1l$fit,col=CL[5],lty=2)
> lines(ug,pt2l$fit,col=CL[2],lty=2)
> legend(0,.1,c("Second Tour","Premier Tour"),col=CL[c(1,6)],
+ lwd=2,lty=1,border=NA)
(en rajoutant une petite légende, avec une visualisation pour les deux tours de l’élection présidentielle, avec un intervalle de confiance sur la prévision de mon taux de procuration).
On peut faire la même chose sur le taux de logement HLM, dans le quartier,
Les dessins sont parlant, mais dans la sortie du modèle de régression, l’interprétation de laisse à désirer (toute suggestion est la bienvenue !).
> summary(regt1l)
Call:
glm(formula = PROCURATIONS/INSCRITS ~ taux, family = binomial,
data = baset1, weights = INSCRITS)
Deviance Residuals:
Min 1Q Median 3Q Max
-12.9549 -1.5722 0.0319 1.6292 13.1303
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.70811 0.01516 -244.6 <2e-16 ***
taux 1.49666 0.04012 37.3 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 12507 on 836 degrees of freedom
Residual deviance: 11065 on 835 degrees of freedom
AIC: 15699
Number of Fisher Scoring iterations: 4
> summary(regt2l)
Call:
glm(formula = PROCURATIONS/INSCRITS ~ taux, family = binomial,
data = baset2, weights = INSCRITS)
Deviance Residuals:
Min 1Q Median 3Q Max
-15.4872 -1.7817 -0.1615 1.6035 12.5596
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.24272 0.01230 -263.61 <2e-16 ***
taux 1.45816 0.03266 44.65 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 9424.7 on 836 degrees of freedom
Residual deviance: 7362.3 on 835 degrees of freedom
AIC: 12531
Number of Fisher Scoring iterations: 4
On a alors voulu comparer avec un modèle qui me semble moins juste, mais qui est plus simple à interpréter, où on suppose que le taux de procuration (par bureau) est expliqué par un modèle linéaire
(que l’on peut aussi lisser, pour vérifier que le lien est effectivement linéaire). Le code est ici
> regt1=lm(PROCURATIONS/INSCRITS~bs(taux,6),weights=INSCRITS,data=baset1)
> regt2=lm(PROCURATIONS/INSCRITS~bs(taux,6),weights=INSCRITS,data=baset2)
> u=seq(min(baset1$taux)+.015,max(baset1$taux)-.015,by=.001)
> ND=data.frame(taux=u)
> ug=seq(0,max(baset1$taux)+.05,by=.001)
> pt1=predict(regt1,newdata=ND,se=TRUE,type="response")
> pt2=predict(regt2,newdata=ND,se=TRUE,type="response")
> library(RColorBrewer)
> CL=brewer.pal(6, "RdBu")
> plot(ug,ug*1,col="white",xlab=nom,ylab="Taux de procuration",
+ ylim=c(0,.1))
> for(i in 1:(length(b$breaks)-1)){
+ polygon(b$breaks[i+c(0,0,1,1)],c(0,b$counts[i],b$counts[i],0)
+ /max(b$counts)*.05,col="light yellow",border=NA)}
> polygon(c(u,rev(u)),c(pt1$fit+2*pt1$se.fit,rev(pt1$fit-2*pt1$se.fit)),
+ border=NA,density=30,col=CL[4])
> lines(u,pt1$fit,col=CL[6],lwd=2)
> polygon(c(u,rev(u)),c(pt2$fit+2*pt2$se.fit,rev(pt2$fit-2*pt2$se.fit)),
+ border=NA,density=30,col=CL[3])
> lines(u,pt2$fit,col=CL[1],lwd=2)
> regt1l=lm(PROCURATIONS/INSCRITS~taux,weights=INSCRITS,data=baset1)
> regt2l=lm(PROCURATIONS/INSCRITS~taux,weights=INSCRITS,data=baset2)
> ND=data.frame(taux=ug)
> pt1l=predict(regt1l,newdata=ND,se=TRUE,type="response")
> pt2l=predict(regt2l,newdata=ND,se=TRUE,type="response")
> lines(ug,pt1l$fit,col=CL[5],lty=2)
> lines(ug,pt2l$fit,col=CL[2],lty=2)
> legend(0,.1,c("Second Tour","Premier Tour"),col=CL[c(1,6)],
+ lwd=2,lty=1,border=NA)
(j’ai tout mis d’un coup cette fois, les modèles lissés et linéaires, l’un à la suite de l’autre)
Cette fois, on a une sortie de régression plus classique,
> summary(regt1l)
Call:
lm(formula = PROCURATIONS/INSCRITS ~ taux, data = baset1, weights = INSCRITS)
Weighted Residuals:
Min 1Q Median 3Q Max
-1.9994 -0.2926 0.0011 0.3173 3.2072
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.021268 0.001739 12.23 <2e-16 ***
taux 0.054371 0.004812 11.30 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.646 on 835 degrees of freedom
Multiple R-squared: 0.1326, Adjusted R-squared: 0.1316
F-statistic: 127.7 on 1 and 835 DF, p-value: < 2.2e-16
> summary(regt2l)
Call:
lm(formula = PROCURATIONS/INSCRITS ~ taux, data = baset2, weights = INSCRITS)
Weighted Residuals:
Min 1Q Median 3Q Max
-2.9029 -0.4148 -0.0338 0.4029 3.4907
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.033909 0.001866 18.17 <2e-16 ***
taux 0.079749 0.005165 15.44 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.6934 on 835 degrees of freedom
Multiple R-squared: 0.2221, Adjusted R-squared: 0.2212
F-statistic: 238.4 on 1 and 835 DF, p-value: < 2.2e-16
On note plusieurs choses de ces graphiques. (i) les deux types de régression donnent des modèles pour les taux de procuration très très proches. Donc autant prendre le plus simple à interpréter. (ii) le lissage n’apporte rien, et le modèle linéaire semble pertinent. On a ainsi regardé plusieurs variables, et on en a retenu un certain nombre, pour tenter un modèle multiple. Avant de parler des résidus de notre modèle, je devrais peut être prendre quelques lignes pour parler un peu de cartographie (une nouvelle fois, Baptiste m’a fait découvrir de belles fonctions).
- Visualiser les bureaux de vote à Paris
Pour récupérer le fond de carte, avec les bureaux de vote, on utilise (là encore, je renvoie au blog de Baptiste, qui explique l’utilisation de données de cartelec.net)
> library(maptools)
> library(rgdal)
> library(classInt)
> paris=readShapeSpatial("paris-cartelec.shp")
Si on veut visualiser, par exemple, le taux de procuration (disons la moyenne entre les deux tours), on utilise les données suivantes
> elec=data.frame()
> elec=cbind(bt1$BV,(bt1$PROCURATIONS+bt2$PROCURATIONS),(bt1$EXPRIMES+bt2$EXPRIMES))
> colnames(elec)=c("BV","PROCURATIONS","EXPRIMES")
> elec=as.data.frame(elec)
> elec$BV=bt1$BV
Ensuite, viennent les fonctions graphiques, où on va passer d’un taux à une classe et d’une classe à une couleur,
> m=match(paris$BUREAU,elec$BV) > plotvar=100*elec$PROCURATIONS/elec$EXPRIMES > nclr=7 > plotclr=brewer.pal(nclr,"RdYlBu")[nclr:1] > class=classIntervals(plotvar[m], nclr, style="fisher",dataPrecision=1) > colcode=findColours(class, plotclr)
Reste à conclure, en faisant une visualisation graphique de nos données
> par(mar=c(1,1,1,1)) > plot(paris,col=colcode,border=colcode) > legend(656274.9, 6867308,legend=names(attr(colcode,"table")), + fill=attr(colcode, "palette"), cex=1, bty="n", + title="Frequence procurations (%)")
Histoire de conclure, on peut regarder un peu nos résidus, obtenus sur un modèle linéaire. Considérons un modèle sur seulement trois variables explicatives,
> regt1=lm(PROCURATIONS/INSCRITS~I(POP65P/POP)+ + I(PROPRIO/LOGEMENT)+I(CS3/POP1564),weights=INSCRITS,data=baset1)
Dans ce cas, la visualisation des résidus donne
> m=match(paris$BUREAU,elec$BV) > plotvar=100*residuals(regt1) > nclr=7 > plotclr=brewer.pal(nclr,"RdYlBu")[nclr:1] > class=classIntervals(plotvar[m], nclr, style="fisher",dataPrecision=1) > colcode=findColours(class, plotclr) > par(mar=c(1,1,1,1)) > plot(paris,col=colcode,border=colcode) > legend(656274.9, 6867308,legend=names(attr(colcode,"table")), + fill=attr(colcode, "palette"), cex=1, bty="n",title="Residus")
Idéalement, il faudrait avoir un beau bruit (spatial), c’est à dire avoir des couleurs réparties aléatoirement sur Paris. Il reste encore pas mal de régions dont les voisins sont de la même couleurs, et on repère quelques quartiers atypiques, avec soit des résidus importants négativement, ou positivement. Comme toujours, en modélisation, on pourrait passer des heures pour essayer de capturer tous les effets, mais aucun modèle avec les variables à notre disposition nous a permis de faire réellement mieux.
A few days in France
Next week, I will spend a few days in France, in Lyon first, and then in Paris. I will try to post the slides and the code that I will use during my talk in Lyon (workshop on climate change) and also some codes on vine copulas (used to get a better understanding of that tool before the PhD defense of Pierre-André).
Playing with robots
My son would be extremely proud if I tell him I can spend hours building robots. Well, my robots are not as fancy as Dr Tenma’s, but they usually do what I ask them to do. For instance, it is extremely simple to build a robot with R, to extract data from websites. I have mentioned it here (one tennis matches), but it failed there (on NY Marathon). To illustrate the use of robots, assume that one wants to build his own dataset to study prices of airline tickets. First, we have to choose a departure city (e.g. Paris) and an arrival city (e.g. Montreal). Then, one wants to look at all possible dates from April first (I ran it last month) till the end of December (so we create a vector with all leaving dates, namely a vector for the day, one for the month, and one for the year). Then, we choose a return date (say 3 days after).
DEP="Paris" ARR="Montreal" DATE1D=rep(c(1:30,1:31,1:30,1:31,1:31,1:30,1:31,1:30, 1:31,1:31,1:29),3) DATE1M=rep(c(rep(4,30),rep(5,31),rep(6,30),rep(7,31), rep(8,31),rep(9,30),rep(10,31),rep(11,30),rep(12,31), rep(1,31),rep(2,29)),3) DATE1Y=rep(c(rep(2011,30+31+30+31+31+30+31+ 30+31+31+28),rep(2012,31+29)),3) k=3 DATE3D=c((1+k):30,1:31,1:30,1:31,1:31,1:30,1:31, 1:30,1:31,1:31,1:29,1:k) DATE3M=c(rep(4,30-k),rep(5,31),rep(6,30),rep(7,31),rep(8,31), rep(9,30),rep(10,31),rep(11,30),rep(12,31),rep(1,31),rep(2,29), rep(3,k)) DATE3Y=c(rep(2011,30+31+30+31+31+30+31+30+31+ 31+28-k),re p(2012,31+29+k))
It is also possible (for a nice robot), to skip all prior dates
skip=max(as.numeric(Sys.Date()-as.Date("2011-04-01")),1)
Then, we need a website where requests can be written nicely (with cities and dates appearing explicitly). Here, I cannot not mention the website that I used since it is stated on the website that it is strictly forbidden to run automatic requests… Anyway, consider a loop create a url address (actually I chose the value of the date randomly, since I had been told that those websites had memory: if you ask too many times for the same thing during a short period of time, prices would go up),
URL=paste("http://www.♦♦♦♦/dest.dll?qscr=fx&flag=q&city1=", DEP,"&citd1=",ARR,"&", "date1=",DATE1D[s],"/",DATE1M[s],"/",DATE1Y[s], "&date2=",DATE3D[s],"/",DATE3M[s],"/",DATE3Y[s], "&cADULT=1",sep="")
then, we just have to scan the webpage, looking for ticket prices (just looking for some specific names)
page=as.character(scan(URL,what="character")) I=which(page%in%c("Price0","Price1","Price2")) if(length(I)>0){ PRIX=substr(page[I+1],2,nchar(page[I+1])) if(PRIX[1]=="1"){PRIX=paste(PRIX,page[I+2],sep="")} if(PRIX[1]=="2"){PRIX=paste(PRIX,page[I+2],sep="")}
Here, we have to be a bit cautious, if prices exceed 1000. Then, it is possible to start a statistical study. For instance, if we compare to destination (from Paris), e.g. Montréal and New York, we obtain the following patterns (with high prices during holidays),
It is also possible to run the code twice (here it was run last month, and a couple of days ago), for the same destination (from Paris to Montréal),
Of course, it would be great if I could run that code say every week, to build up a nice dataset, and to study the dynamic of prices…
The problem is that it is forbidden to do this. In fact, on the website, it is mentioned that if we want to extract data (for an academic purpose), it is possible to ask for an extraction. But if we do tell that we study specific prices, data might be biased. So the good idea would be to use several servers, to make several requests, randomly, and to collect them (changing dates and destination). But here, my computing skills – unfortunately – reach a limit….
Warming in Paris: minimas versus maximas ?
Recently, I received comments (here and on Twitter) about my previous graphs on the temperature in Paris. I mentioned in a comment (there) that studying extremas (and more generally quantiles or interquantile evolution) is not the same as studying the variance. Since I am not a big fan of the variance, let us talk a little bit about extrema behaviour.
In order to study the average temperature it is natural to look at the linear (assuming that it is linear, but I proved that it could reasonably be assumed as linear in the paper) regression, i.e. least square regression, which gives the expected value. But if we care about extremes, or almost extremes, it is natural to look at quantile regression.
For instance, below, the green line is the least square regression, the red one is 97.5% quantile, and the blue on the 2.5% quantile regression.
It looks like the slope is the same, i.e. extremas are increasing as fast as the average…
tmaxparis=read.table("temperature/TG_SOUID100845.txt", skip=20,sep=",",header=TRUE) head(tmaxparis) Dparis=as.Date(as.character(tmaxparis$DATE),"%Y%m%d") Tparis=as.numeric(tmaxparis$TG)/10 Tparis[Tparis==-999.9]=NA I=sample(1:length(Tparis),size=5000,replace=FALSE) plot(Dparis[I],Tparis[I],col="grey") abline(lm(Tparis~Dparis),col="green") library(quantreg) abline(rq(Tparis~Dparis,tau=.025),col="blue") abline(rq(Tparis~Dparis,tau=.975),col="red")
(here I plot randomly some points to avoid a too heavy figure, since I have too many observations, but I keep all the observations in the regression !).
Now, if we look at the slope for different quantile level (Fig 6 in the paper, here, but on minimum daily temperature, here I look at average daily temperature), the interpretation is different.
s=0 COEF=SD=rep(NA,199) for(i in seq(.005,.995,by=.005)){ s=s+1 REG=rq(Tparis~Dparis,tau=i) COEF[s]=REG$coefficients[2] SD[s]=summary(REG)$coefficients[2,2] }
with the following graph below,
s=0 plot(seq(.005,.995,by=.005),COEF,type="l",ylim=c(0.00002,.00008)) for(i in seq(.005,.995,by=.005)){ s=s+1 segments(i,COEF[s]-2*SD[s],i,COEF[s]+2*SD[s],col="grey") } REG=lm(Tparis~Dparis) COEFlm=REG$coefficients[2] SDlm=summary(REG)$coefficients[2,2] abline(h=COEFlm,col="red") abline(h=COEFlm-2*SDlm,lty=2,lw=.6,col="red") abline(h=COEFlm+2*SDlm,lty=2,lw=.6,col="red")
Here, for minimas (quantiles associated to low probabilities, on the left), the trend has a higher slope than the average, so in some sense, w
arming of minimas is stronger than average temperature, and on other hand, for maximas (high probabilities on the right), the slope is smaller – but positive – so summer are warmer, but not as much as winters.
Note also that the story is different for minimal temperature (mentioned in the paper) compared with that study, made here on average daily temperature (see comments)… This is not a major breakthrough in climate research, but this is all I got…
More climate extremes, or simply global warming ?
In the paper on the heat wave in Paris (mentioned here) I discussed changes in the distribution of temperature (and autocorrelation of the time series).
During the workshop on Statistical Methods for Meteorology and Climate Change today (here) I observed that it was still an important question: is climate change affecting only averages, or does it have an impact on extremes ? And since I’ve seen nice slides to illustrate that question, I decided to play again with my dataset to see what could be said about temperature in Paris.
Recall that data can be downloaded here (daily temperature of the XXth century).
tmaxparis=read.table("/temperature/TX_SOUID100124.txt", skip=20,sep=",",header=TRUE) Dmaxparis=as.Date(as.character(tmaxparis$DATE),"%Y%m%d") Tmaxparis=as.numeric(tmaxparis$TX)/10 tminparis=read.table("/temperature/TN_SOUID100123.txt", skip=20,sep=",",header=TRUE) Dminparis=as.Date(as.character(tminparis$DATE),"%Y%m%d") Tminparis=as.numeric(tminparis$TN)/10 Tminparis[Tminparis==-999.9]=NA Tmaxparis[Tmaxparis==-999.9]=NA annee=trunc(tminparis$DATE/10000) MIN=tapply(Tminparis,annee,min) plot(unique(annee),MIN,col="blue",ylim=c(-15,40),xlim=c(1900,2000)) abline(lm(MIN~unique(annee)),col="blue") abline(lm(Tminparis~unique(Dminparis)),col="blue",lty=2) annee=trunc(tmaxparis$DATE/10000) MAX=tapply(Tmaxparis,annee,max) points(unique(annee),MAX,col="red") abline(lm(MAX~unique(annee)),col="red") abline(lm(Tmaxparis~unique(Dmaxparis)),col="red",lty=2)
On the plot below, the dots in red are the annual maximum temperatures, while the dots in blue are the annual minimum temperature. The plain line is the regression line (based on the annual max/min), and the dotted lines represent the average maximum/minimum daily temperature (to illustrate the global tendency),
It is also possible to look at annual boxplot, and to focus either on minimas, or on maximas.
annee=trunc(tminparis$DATE/10000) boxplot(Tminparis~as.factor(annee),ylim=c(-15,10), xlab="Year",ylab="Temperature",col="blue") x=boxplot(Tminparis~as.factor(annee),plot=FALSE) xx=1:length(unique(annee)) points(xx,x$stats[1,],pch=19,col="blue") abline(lm(x$stats[1,]~xx),col="blue") annee=trunc(tmaxparis$DATE/10000) boxplot(Tmaxparis~as.factor(annee),ylim=c(15,40), xlab="Year",ylab="Temperature",col="red") x=boxplot(Tmaxparis~as.factor(annee),plot=FALSE) xx=1:length(unique(annee)) points(xx,x$stats[5,],pch=19,col="red") abline(lm(x$stats[5,]~xx),col="red")
Plain dots are average temperature below the 5% quantile for minima, or over the 95% quantile for maxima (again with the regression line),
We can observe an increasing trend on the minimas, but not on the maximas !
Finally, an alternative is to remember that we focus on annual maximas and minimas. Thus, Fisher and Tippett theory (mentioned here) can be used. Here, we fit a GEV distribution on a blog of 10 consecutive years. Recall that the GEV distribution is
install.packages("evir") library(evir) Pmin=Dmin=Pmax=Dmax=matrix(NA,10,3) for(s in 1:10){ X=MIN[1:10+(s-1)*10] FIT=gev(-X) Pmin[s,]=FIT$par.ests Dmin[s,]=FIT$par.ses X=MAX[1:10+(s-1)*10] FIT=gev(X) Pmax[s,]=FIT$par.ests Dmax[s,]=FIT$par.ses }
The location parameter 
is the following, with on the left the minimas and on the right the maximas,
while the scale parameter 
is
and finally the shape parameter 
is
On those graphs, it is very difficult to say anything regarding changes in temperature extremes… And I guess this is a reason why there is still active research on that area…
Des étés pluvieux en Bretagne ? une réalité statistique…
Pour compléter le précédant billet (ici) on peut se demander en quoi la Bretagne est différente des autres régions françaises… Nous avions vu ici le niveau de précipitation moyen, jour après jours pendant les mois d’été, en Bretagne. A Rennes. En revanche, à Paris on obtient la moyenne suivante,
que l’on peut comparer à Marseille,
ou encore à Strasbourg,
Sur la figure ci-dessous, on voit que la probabilité d’avoir de la pluie à Paris (au sens au moins 0.1 mm d’eau dans la journée, en trait gras bleu, au moins 2 mm d’eau dans la journée, en trait bleu) est supérieure à la probabilité d’avoir de la pluie à Rennes (respectivement en bleu clair gras, et en bleu clair fin)
On est certes très au dessus de Marseille,
mais très en dessous de Strasbourg,
Mais au delà des lois marginales, ces villes sont différentes de la Bretagne si l’on regarde les matrices de transition.
- Transition d’un jour sur l’autre
Pour Rennes, si on regarde jour après jour, on obtient
|
||||
|
|
|||
|
|
1955 | 612 | 2567 |
|
606 | 723 | 1329 | |
ce qui donne les probabilités de transition suivantes,
|
|||
|
|
||
|
|
76,15 % | 23,85 % |
|
45,60 % | 54,40 % | |
Pour Paris, la probabilité de transition jour après jour a la forme suivante
|
||||
|
|
|||
|
|
2689 | 959 | 3648 |
|
946 | 1466 | 2412 | |
ce qui donne les probabilités de transition suivantes,
|
|||
|
|
||
|
|
73,71 % | 26,29 % |
|
39,22 % | 60,78 % | |
Pour Marseille, la probabilité de transition jour après jour a la forme suivante
|
||||
|
|
|||
|
|
2527 | 375 | 2902 |
|
362 | 216 | 578 | |
ce qui donne les probabilités de transition suivantes,
|
|||
|
|
||
|
|
87,08 % | 12,92 % |
|
62,63 % | 37,37 % | |
Pour Strasbourg, la probabilité de transition jour après jour a la forme suivante
|
||||
|
|
|||
|
|
31 | 128 | 159 |
|
132 | 1464 | 1596 | |
ce qui donne les probabilités de transition suivantes,
|
|||
|
|
||
|
|
19,50 % | 80,50 % |
|
8,27 % | 91,73 % | |
- Transition d’une semaine sur l’autre
Si en revanche on regarde les matrices de transition semaine par semaine, on a des résultats assez différents. Une bonne semaine signifie aucun jour avec plus de 2 mm de pluie.
Pour Rennes, si on regarde semaine après semaine
|
||||
|
|
|||
|
|
379 | 25 | 404 |
|
26 | 7 | 33 | |
ce qui donne les probabilités de transition suivantes,
|
|||
|
|
||
|
|
93,81 % | 6,19 % |
|
78,79 % | 21,21 % | |
Pour Paris, la probabilité de transition semaine après semaine a la forme suivante
|
||||
|
|
|||
|
|
576 | 46 | 622 |
|
53 | 4 | 57 | |
ce qui donne les probabilités de transition suivantes,
|
|||
|
|
||
|
|
92,60 % | 7,40 % |
|
92,98 % | 7,02 % | |
Pour Marseille, la probabilité de transition semaine après semaine a la forme suivante
|
||||
|
|
|||
|
|
274 | 59 | 333 |
|
47 | 9 | 56 | |
ce qui donne les probabilités de transition suivantes,
|
|||
|
|
||
|
|
82,28 % | 17,72 % |
|
83,93 % | 16,07 % | |
Pour Strasbourg, la probabilité de transition semaine après semaine a la forme suivante
|
||||
|
|
|||
|
|
1494 | 614 | 2018 |
|
613 | 939 | 1552 | |
ce qui donne les probabilités de transition suivantes,
|
|||
|
|
||
|
|
70,87 % | 29,13 % |
|
39,50 % | 60,50 % | |
Les tests du chi deux, d’indépendance d’une semaine sur l’autre donnent
- à Rennes, une statistique du chi-deux de 8,054, soit une p-value de 0,45%
- à Paris, une statistique du chi-deux de 0,025, soit une p-value de 87,26%
- à Marseille, une statistique du chi-deux de 0,012, soit une p-value de 91,24%
- à Strasbourg, une statistique du chi-deux de 0,7649, soit une p-value de 38,18%
autrement dit l’hypothèse d’indépendance est acceptée partout, sauf à Rennes….
- Moralité ?
De manière assez paradoxale, on prétend que la Bretagne a un temps changeant, et pour reprendre le titre du précédant billet, effectivement, en Bretagne, il peut faire beau plusieurs fois par jour. Mais sur le long terme, d’une semaine sur l’autre, le temps est au contraire très corrélé, contrairement aux autres régions. A Paris, Marseille ou Strasbourg, qu’il ait fait beau ou qu’il ait plus la semaine précédente, cela n’apporte aucune information sur la probabilité d’avoir de la pluie la semaine où l’on vient en vacances…. Mais pas en Bretagne: manifestement, il existe donc des étés pourris, où il pourra pleuvoir toutes les semaines, et des étés superbes où il ne pleut jamais….
Soutenance de thèse, Paris II
Soutenance de thèse de Noureddine Ben Lagha, à Paris II Assas, sur l’assurance automobile,
Son directeur de thèse était Michel Grun-Rehomme, avec dans le jury, Georges Bresson, Alain Trognon, Michel, Marc Fleurbaey et Denis Fougère.
