Tag Archives: Montréal

Econometrics Seminars at Université de Montréal

This Thursday, I will present at the CIREQ Séminaire Marcel-Dagenais en Économétrie at Université de Montréal, ou paper Optimal Transport for Counterfactual Estimation: A Method for Causal Inference, written with Emmanuel Flachaire and Ewen Gallic.

Many problems ask a question that can be formulated as a causal question: “what would have happened if…?” For example, “would the person have had surgery if he or she had been Black?” To address this kind of questions, calculating an average treatment effect (ATE) is often uninformative, because one would like to know how much impact a variable (such as skin color) has on a specific individual, characterized by certain covariates. Trying to calculate a conditional ATE (CATE) seems more appropriate. In causal inference, the propensity score approach assumes that the treatment is influenced by x, a collection of covariates. Here, we will have the dual view: doing an intervention, or changing the treatment (even just hypothetically, in a thought experiment, for example by asking what would have happened if a person had been Black) can have an impact on the values of x. We will see here that optimal transport allows us to change certain characteristics that are influenced by the variable we are trying to quantify the effect of. We propose here a mutatis mutandis version of the CATE, which will be done simply in dimension one by saying that the CATE must be computed relative to a level of probability, associated to the proportion of x (a single covariate) in the control population, and by looking for the equivalent quantile in the test population. In higher dimension, it will be necessary to go through transport, and an application will be proposed on the impact of some variables on the probability of having an unnatural birth (the fact that the mother smokes, or that the mother is Black).

Slides are available online. I will try to mention additional papers published this year, such as Fairness in Multi-Task Learning via Wasserstein Barycenters, Mitigating Discrimination in Insurance with Wasserstein Barycenters or more recently A Sequentially Fair Mechanism for Multiple Sensitive Attributes.

Snow in Montréal (Canada)

Winter started a bit more than one month ago… but we have already experienced many snow storms… there is still a lot snow in gardens and in the streets,

I was wondering if it was that unusual, but apparently not. Compared with last year, it is (for the first months of winter, until the end of Januray), it +50%, but it is comparable with previous years

Yes, we a simple loop, we can easily extract data from official wesite https://climat.meteo.gc.ca/ (but not too far away, even 2015 contains a lot of missing observations). For this month, we use

url = "https://climat.meteo.gc.ca/climate_data/daily_data_f.html?StationID=51157&timeframe=2&StartYear=1840&EndYear=2023&Day=30&Year=2023&Month=1#"
download.file(url,destfile = "M.html")
k = which(tables[[1]]$`JOUR `=="Somme")
neige = tables[[1]]$`Neige tot. Definitioncm `[k]
x = as.numeric(sub(",", ".", strsplit(neige, "LegendCarer")[[1]][1], fixed = TRUE))

and then we loop, and store the number we look for in a data frame (yes, we have to convert “50,8LegendCarer^” into the appropriate numerical value (that would be here 50.8

D = data.frame(annee = c(2023,rep(2022:2015,each=12),c(12,11,10)), mois= c(1,rep(12:1,8),12,11,10), lab = neige, snow = x)
for(i in 2:nrow(D)){
    y = D$annee[i]
    m = D$mois[i]
    url = paste("https://climat.meteo.gc.ca/climate_data/daily_data_f.html?StationID=51157&timeframe=2&StartYear=1840&EndYear=2023&Day=30&Year=",y,"&Month=",m,"#",sep="")
  download.file(url,destfile = "M.html")
  k = which(tables[[1]]$`JOUR `=="Somme")
  neige = tables[[1]]$`Neige tot. Definitioncm `[k]
  x = as.numeric(sub(",", ".", strsplit(neige, "LegendCarer")[[1]][1], fixed = TRUE))
  D[i,3] = neige
  D[i,4] = x

Here are the most recent months

> head(D)
  annee mois              lab snpw
1  2023    1 50,8LegendCarer^ 50.8
2  2022   12             63,0 63.0
3  2022   11             14,6 14.6
4  2022   10              0,0  0.0
5  2022    9              0,0  0.0
6  2022    8              0,0  0.0

Of course, we need some codes to plot, we here, I mainly wanted to keep tracks of the code used to extract meteorological data…


Montréal AI Symposium 2022

In about ten days (Saturday afternoon), I will be presenting a poster on fairness, discrimination and insurance at the Montréal AI Symposium, based on our joint paper The Fairness of Machine Learning in Insurance: New Rags for an Old Man?, written with Laurence Barry. Since the paper was quite literary, I used material from the document Insurance: Discrimination, Biases & Fairness to get more a visual poster. Additional information will come while discussing…

This is a poster used at the conference

Radial Graphs for Time Series

On How to: Weather Radials, there was a nice visualisation of temperatures. Since I am too old fashioned for ggplot2, I wanted to reproduce a similar graph with the old plot style.

Assume that daily temperature is in a vector X (e.g. temperature in Montréal, QC, in 2009). To get a radial plot, use

> n=length(X)
> theta=seq(0,1-1/n,length=n)*2*pi
> r=30+X
> plot(r*cos(pi/2-theta),r*sin(pi/2-theta),type="l",xlab="",ylab="",axes=FALSE)
> for(t in 1:n){
+   if(X[t]>0) CL=rgb(0,0,1,.4)
+   if(X[t]<0) CL=rgb(1,0,0,.4)
+   if(X[t]==0) CL="white"
+   segments((30+X[t])*cos(pi/2-theta[t]),(30+X[t])*sin(pi/2-theta[t]),30*cos(pi/2-theta[t]),30*sin(pi/2-theta[t]),col=CL)
+ }
> for(r in 10*seq(0,6)) lines(r*cos(pi/2-theta),r*sin(pi/2-theta),type="l",col="light blue")

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

A random walk ? What else ?

Consider the following time series,

What does it look like ? I know, this is a stupid game, but I keep using it in my time series courses. It does look like a random walk, doesn’t it ? If we use Philipps-Perron test, yes, it does,

> PP.test(x)

	Phillips-Perron Unit Root Test

data:  x 
Dickey-Fuller = -2.2421, Truncation lag parameter = 6, p-value = 0.4758

If we look at the autocorrelation function, we do observe some persistence,

> acf(x,100)

Perhaps this persistence can be related to long range dependence, or to some fractional random walk. A natural idea could be estimate Hurst parameter, using for instance Beran (1992) estimator – based on Whittle (1956) – where we assume that the autocorrelation function satisfies

as  for some  (the so called Hurst index). But here, we start to observe unexpected ouputs,

> library(longmemo)
> (d  <- WhittleEst(x))
'WhittleEst' Whittle estimator for  fractional Gaussian noise ('fGn');	 call:
WhittleEst(x = x)
	  time series of length  n = 759.

H = 0.9899335
coefficients 'eta' =
    Estimate Std. Error z value   Pr(>|z|)
H 0.98993350 0.02468323 40.1055 < 2.22e-16
 <==> d := H - 1/2 = 0.49 (0.025)

 $ vcov       : num [1, 1] 0.000609
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr "H"
  .. ..$ : chr "H"
 $ periodogr.x: num [1:379] 1479.3 1077.3 371.7 287.2 51.2 ...
 $ spec       : num [1:379] 62.5 31.7 21.3 16.1 12.9 ...

or more precisely some non-expected values for Hurst parameter, which should be in 

> confint(d)
      2.5 %   97.5 %
H 0.9415553 1.038312

Oops, perhaps, we did miss something, because it looks like there is extremely strong persistence on our time series,

> plot(d)

It is probablty time to ask where I found that series… To be honest, I did borrow  it from a great canadian website http://climate.weatheroffice.gc.ca/climateData/. For instance, it you want the temperature we did experience a few days ago, you can use

> Y=2013
> M=1
> D=25
> url=paste(
> page=scan(url,what="character")

Yes, that series is the temperature we did experience in Montréal last month (hourly time seies). On the graph below, you can actually compare it with temperature experienced in Januarys over the past 60 years,

So it is not that surprising to see long range dependence models appearing (I did write a paper on that topic precisely a few years ago). What I found puzzeling is that persistence is large, extremely large. And the problem is that I do not see how we can explain ‘jumps’ that we do observe on that series. For instance the behavior of the series while I was in Europe, before January 20th: within 3 days, the temperature went down, from 0°C to -20°C, and up from -20°C to 0°C, and then down again, from 0°C to -20°C (a nice И if we use cyrillic letters). Or how can we explain the oscillating behavior observed the week after, where the temperature went up, from -25°C to (almost) +10°C in a few days. Within 10 days, we did observe also two ‘jumps’ (or ‘crashes‘ if we want to use the terminology of financial time series) with a decrease of 25 degrees in less than 24 hours ! Obviously, we need to find other classes of model to replicate that kind of behavior we observe on temperatures…

Econometric Modeling in Finance and Insurance with the R language

On February 15th, IFM2, the Institute of Financial Mathematics in Montréal will organize an (one day) Executive workshop on Econometric Modeling in Finance and Insurance with the R language. The event is not yet mentioned in the calendar, but the syllabus can be downloaded here. Additional details (slides and R code) will be available soon, on this blog. In the morning, it will be an introduction to the R langage, and in the afternoon, we will focus on applications,

  • Principal components analysis and application to yield curves
  • Regression tree, logistic regression and application to credit scoring
  • Poisson regression and applications to claims reserving (IBNR) and projected mortality tables (LifeMetrics)

Comparer des notes entre pays

A Montréal, l’an passé, le principal souvenir que je garderai de mon cours est le mal que j’ai eu à comprendre comment mettre des notes. Je ne ferais pas un billet sur la difficulté à évaluer le niveau d’un étudiant (car j’ai toujours détesté les examens, aussi bien en tant qu’étudiant que prof), mais j’ai une une réelle difficulté à comprendre ce que vaut une note… Par exemple, un C, ça correspond à quoi (avec mon référentiel de prof français) ? C’est bien ou pas ? C’est quoi la moyenne ?
La seule information dont on dispose pour mettre les notes est le tableau ci-contre. Et autant je pense voir ce qu’est un travail “excellent“, ou “faible“, mais “passable” ? Car le tableau me dit qu’un élève qui a 50 sur 100 a une note passable, ce qui est aussi la référence que l’on utilise en France. Sauf que les élèves ont commence à me dire que 50, c’est pas terrible… qu’une note moyenne ici c’était plutôt 60, voire 65. Bref, sous-entendu (voire entendu, tout court), “monsieur, vous devez remonter les notes“. Et comme je l’ai prouve ici, il ne faut jamais croire les élèves…
Et effectivement, dans un vieux rapport datant des années 90 (en maths et informatiques), je suis tombé sur le graphique suivant
qui donne la distribution “typique” des notes. Autrement dit, la médiane est à 63,ce qui signifie que plus de moitie a une note qui est “bien“,”très bien” voire “excellent“. En revanche, une distribution comme celle dessous est considérée comme “dure et discriminatoire
Pour mieux comprendre les notations au Canada et en France, j’ai voulu faire une petite étude (non, je ne reviendrai pas sur les notes). On continue aujourd’hui àtrouver ce genre de distribution à Montréal. Par exemple, si on considère les cours, analyse (en première année), algèbre linéaire et calcul, on a les histogrammes suivants

ce qui donne le boxplot suivant

i.e. pour les cours algèbre linéaire et calcul on est effectivement sur des moyennes autour de 60. Prenons maintenant, a titre de comparaison, des cours de maths, dans une université française, en l’occurrence a Rennes, comme analyse (en première année), arithmétique et algèbre linéaire, on a des histogrammes beaucoup plus plats,

A titre de comparaison, j’ai été voir du cote du département de Sciences Économiques, avec macroéconomie, microéconomie et mathématiques appliquées, on a

i.e. on retrouve de belles gaussiennes pour les cours d’économie, mais une distribution assez différente pour le cours de maths,

On a effectivement des distributions assez différentes, généralement, en première année à l’université, entre le Canada et la France. Ce qui explique que les profs, lorsque les étudiants souhaitent passer d’un pays a l’autre ne regarde plus vraiment les notes, mais plus le rang de l’étudiant au sein de sa promotion. D’ailleurs, si on fait un graphique quantile-quantile, avec les cours de maths à Montréal versus les cours de maths a Rennes, on a le graphique suivant

autrement dit, une note de 50 à Montréal correspond à 33 en France (soit 6.5 sur 20), alors que 50 en France (10 sur 20) correspond a 59 à Montréal. On retrouve la encore le fait qu’une note médiane correspond à 60 ici… c’est a dire un C. Il va falloir que je m’adapte l’année prochaine….


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


It is also possible (for a nice robot), to skip all prior dates


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),


then, we just have to scan the webpage, looking for ticket prices (just looking for some specific names)


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

911, jour après jour

Après deux billets (ici puis ) sur les cycles intrajournaliers des appels au 911, on peut se demander comment les crimes évoluent au cours de la semaine.

Pour l’ensemble des appels passés au 911, on a la distribution suivante

i.e. un pic les vendredi soir et samedi soir, et un creux le dimanche. Si on regarde les appels pour des cambriolages, on a la distribution suivante

avec des pics en matinée, les vendredi après midi, et les fins de semaine. On peut aussi suivre les troubles de la paix,

qui surviennent certes vers minuit, mais essentiellement en fin de semaine. Ce qui contraste assez avec les hold-ups,

Manifestement, il y a des tendances assez claires. La prochaine étape sera de regarder un peu les saisons, ou mieux, l’impact du climat…

à suivre donc…