# Testing for a causal effect (with 2 time series)

A few days ago, I came back on a sentence I found (in a French newspaper), where someone was claiming that

“… an old variable explains 85% of the change in a new variable. So we can talk about causality”

and I tried to explain that it was just stupid : if we consider the regression of the temperature on day $t+1$ against the number of cyclist on day $t$, the $R^2$ exceeds 80%… but it is hard to claim that the number of cyclists on specific day will actually cause the temperature on the next day… Nevertheless, that was frustrating, and I was wondering if there was a clever way to test for causality in that case. A popular one is Granger causality (I can mention a paper we published a few years ago where we use such a test, Tents, Tweets, and Events: The Interplay Between Ongoing Protests and Social Media). To explain that test, consider a bivariate time series (just like the one we have here), $\boldsymbol{z}_t=(x_t,y_t)$, and consider some bivariate autoregressive model
$${\begin{bmatrix}x_{t}\\y_{t}\end{bmatrix}}={\begin{bmatrix}c_{1}\\c_{2}\end{bmatrix}}+{\begin{bmatrix}a_{1,1}&\textcolor{red}{a_{1,2}}\\\textcolor{blue}{a_{2,1}}&a_{2,2}\end{bmatrix}}{\begin{bmatrix}x_{t-1}\\y_{t-1}\end{bmatrix}}+{\begin{bmatrix}u_{t}\\v_{t}\end{bmatrix}}$$where $\boldsymbol{\varepsilon}_t=(u_t,v_t)$ is some bivariate white noise, in the sense that (i) $\mathbb{E} (\boldsymbol{\varepsilon}_{t})=\boldsymbol{0}$ (the noise is centered) (ii) $\mathbb{E} (\boldsymbol{\varepsilon}_{t}\boldsymbol{\varepsilon}_{t}^\top)=\Omega }$, so the variance matrix is constant, but possibly non-diagonal (iii) $\mathbb{E} (\boldsymbol{\varepsilon}_{t}\boldsymbol{\varepsilon}_{t-h}^\top)=\boldsymbol{0} }$ for all $h\neq 0$. Note that we can use the simplified expression$${\boldsymbol{z}_t=\boldsymbol{c}+\boldsymbol{A}\boldsymbol{z}_{t-1}+\boldsymbol{\varepsilon}_t}$$Now, Granger test is based on several quantities. With off-diagonal terms of matrix $\Omega$, we have a so-called instantaneous causality, and since $\Omega$ is symmetry, we will write $x\leftrightarrow y$. With off-diagonal terms of matrix $\boldsymbol{A}$, we have a so-called lagged causality, with either $\textcolor{blue}{x\rightarrow y}$ or $\textcolor{red}{x\leftarrow y}$ (and possibly both, if both terms are significant).

So I wanted to try on my two-variable problem.

df = read.csv("cyclistsTempHKI.csv") dfts = cbind(C=ts(df$cyclists,start = c(2014, 1,2), frequency = 365), T=ts(df$meanTemp,start = c(2014, 1,2), frequency = 365)) library(vars)

I now have “time series” objects, and we can fit a VAR model,

var2 = VAR(dfts, p = 1, type = "const") coefficients(var2) $C Estimate Std. Error t value Pr(&gt;|t|) C.l1 0.8684009 0.02889424 30.054460 8.080226e-107 T.l1 70.3042012 20.07247411 3.502518 5.102094e-04 const 807.6394001 187.75472482 4.301566 2.110412e-05$T Estimate Std. Error t value Pr(&gt;|t|) C.l1 0.0003865391 6.257596e-05 6.177118 1.540467e-09 T.l1 0.6611135594 4.347074e-02 15.208241 6.086394e-42 const -1.6413074565 4.066184e-01 -4.036481 6.446018e-05

For instant, we can run a causality, to test if the number of cyclists can cause the temperature (on the next day)

causality(var2, cause = "C") $Granger Granger causality H0: C do not Granger-cause T data: VAR object var2 F-Test = 38.157, df1 = 1, df2 = 842, p-value = 1.015e-09 Here, we should clearly reject $H_0$, which is that there is no causal effect. Which is the way statisticians say that there should be some causal effect between the number of cyclist and the temperature… So clearly, something is wrong here. Either it is some sort of superpower that cyclists are not aware of. Or this test that was used for forty years (Clive Granger even got a Nobel price for it) is not working. Or we missed something. Actually… I think we missed something here. Possibly because the series are not stationary. We can almost see it with Phi = matrix(c(coefficients(var2)$C[1:2,1],coefficients(var2)$T[1:2,1]),2,2) eigen(Phi) eigen() decomposition$values  0.9594810 0.5700335

where the highest eigenvalue is very close to one. But actually, we look here at the temperature…

plot(dfts) so, at least, we should expect some seasonal unit root here. So let us use two techniques. The first one is a classical one-year difference, $\Delta_{365}\boldsymbol{z}_t=\boldsymbol{z}_t-\boldsymbol{z}_{t-365}$

var2 = VAR(diff(dfts,365), p = 1, type = "const") coefficients(var2) $C Estimate Std. Error t value Pr(&gt;|t|) C.l1 0.8376424 0.07259969 11.537823 1.993355e-16 T.l1 42.2638410 28.58783276 1.478386 1.449076e-01 const -507.5514795 219.40240747 -2.313336 2.440042e-02$T Estimate Std. Error t value Pr(&gt;|t|) C.l1 0.000518209 0.0003277295 1.5812096 1.194623e-01 T.l1 0.598425288 0.1290511945 4.6371154 2.162476e-05 const 0.547828079 0.9904263469 0.5531235 5.823804e-01

The test on the fited VAR model yields

causality(var2, cause = "C") $Granger Granger causality H0: C do not Granger-cause T data: VAR object var2 F-Test = 2.5002, df1 = 1, df2 = 112, p-value = 0.1167 i.e., with a 11% $p$-value, we should reject the assumption that the number of cyclists cause the temperature (on the next day), and actually, we should also reject the other way causality(var2, cause = "T")$Granger   Granger causality H0: T do not Granger-cause C   data: VAR object var2 F-Test = 2.1856, df1 = 1, df2 = 112, p-value = 0.1421

Nevertheless, if we look at the instantaneous causality, this one makes more sense

$Instant H0: No instantaneous causality between: T and C data: VAR object var2 Chi-squared = 13.081, df = 1, p-value = 0.0002982 The second idea would be to use a one day difference, $\Delta_{1}\boldsymbol{z}_t=\boldsymbol{z}_t-\boldsymbol{z}_{t-1}$ and to fit a VAR model on that one VARselect(diff(dfts,1), lag.max = 4, type="const")$selection AIC(n) HQ(n) SC(n) FPE(n) 3 3 2 3

but on that one, a VAR(1) model – with only one lag – might not be sufficient. It might be better to consider a VAR(3)

var2 = VAR(diff(dfts,1), p = 3, type = "const")

and on that one, one more time, we should reject the causal effect of the number of cyclists on the temperature (on the next day)

causality(var2, cause = "C") $Granger Granger causality H0: C do not Granger-cause T data: VAR object var2 F-Test = 0.67644, df1 = 3, df2 = 828, p-value = 0.5666 and this time, there could be a (lagged) causal effect of the temperature on the number of cyclists causality(var2, cause = "T")$Granger   Granger causality H0: T do not Granger-cause C   data: VAR object var2 F-Test = 7.7981, df1 = 3, df2 = 828, p-value = 3.879e-05   $Instant H0: No instantaneous causality between: T and C data: VAR object var2 Chi-squared = 55.83, df = 1, p-value = 7.905e-14 but nothing instantaneously… So it looks like Granger causality performs well on that one ! # Devoir, Méthodes de Prévisions (consommation électrique) Pour le projet de l’UE de prévision, je demande une prévision de la consommation d’électricité. Pour construire les modèles, je donne accès à un gros fichier qui donne, pour toute la France, la consommation d’électricité, entre 2012 et 2014, heure par heure. • comme en cours, la première étape sera de proposer une série de 52 prévisions pour chaque semaine de 2015, • la seconde étape sera de proposer deux séries de 31 prévisions, pour chaque jour d’août 2015 et de décembre 2015 • enfin, la dernière étape sera une prévision de 168 prévisions, pour chaque heure de la semaine du 19 au 25 octobre 2015 Pour aider (?), je mets un lien vers la température en France dans quelques villes, paris, toulouse, grenoble, nantes, dunkerque (via ftp.ncdc.noaa.gov). Il sera possible d’utiliser la vraie température de 2015 pour faire les prévisions de consommation électrique. Idéalement, je veux un fichier html (markdown) et 4 fichiers csv avec les 4 séries de prévision. # 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") # On NCDF Climate Datasets Mid november, a nice workshop on big data and environment will be organized, in Argentina, We will talk a lot about climate models, and I wanted to play a little bit with those data, stored on http://dods.ipsl.jussieu.fr/mc2ipsl/. # Variation de Température Hier, je suis tombé (via limportant.fr/) sur un documentaire intéressant, en ligne sur francetvinfo.fr/monde/environnement/. Mais le passage du début (retranscrit sur le site) m’a laissé une impression très étrange, Au Groenland, la glace fond à vue d’œil. Cette année, le thermomètre est passé à 25 degrés au-dessus de 0. Il y a huit ans, pour la même période, le blizzard soufflait et les scientifiques devaient affronter des températures de – 35 degrés. Une amplitude de 65 degrés inquiétante pour les chercheurs qui observent la banquise depuis plus de 25 ans. Ça veut dire quoi “inquiétant” ? Est-ce un phénomène nouveau ? Il y presque 5 ans, quand je suis arrivé à Montréal, j’avais mis en ligne un rapide billet comparant les températures à Rennes et à Montréal. En particulier, il y avait ces deux figures, avec les températures annuelles à Rennes, avec en rouge un quantile supérieur observé, et un bleu un quantile inférieur (je ne prends pas les maximum et minimum pour lisser un peu). Disons que des écarts de 20 degrés, quelle que soit la période de l’année, ne sont pas impossible, loin de là. On peut parfaitement avoir 30 degrés un été, et 10 l’été suivant (ou disons, pour être plus précis, un 14 juillet il peut faire 30 degré, et le 14 juillet suivant 10). A Montréal, on avait L’écart est ici un peu plus important. Contrairement à Rennes, il est plus important l’hiver que l’été. Et encore, à cette époque, le printemps érable n’avait pas eu lieu, avec presque 25 degrés en mars, alors qu’il est possible d’atteindre encore les -20 (ce qui fait un écart de 45 degrés). Je peux d’ailleurs remettre en ligne un billet que j’avais écrit, en demandant si la courbe des températures n’était pas une marche aléatoire. Des écarts de température de plus de 40 degrés ne sont pas rares. Entre un maximum pour une journée, et le minimum une autre années. Le plus troublant (ce n’est que mon expérience) c’était plutôt de gagner 40 degrés en une semaine, passer de -30 à +10 (puis replonger deux jours plus tard à -20) Qu’en est-il du Groenland ? Sur eca.knmi.nl/dailydata, on peut récupérer les données journalières de température, de vent, etc, pour plus de 60 stations météo au Groenland. En faisant un petit code, on peut visualiser toutes les séries (en rouge et en bleu, on a les années les plus chaude et plus froide, en moyenne, moyennant que l’on ne tienne pas compte des valeurs manquantes). Le code est assez simple setwd('/home/Documents/temperature-greenland/') fichiers=list.files() for(i in 1:length(fichiers)){ sc=scan(fichiers[i],what="char")[50:150] i1=which(sc=="[DENMARK],")-1 i2=which(sc=="(STAID:")-1 station=paste(sc[i1:i2],collapse=" ") temp=read.table(fichiers[i],skip=20,sep=',', header=TRUE) date=as.Date(as.character(temp$DATE),format=
"%Y%m%d")
m=format(date, "%m")
d=format(date, "%d")
y=format(date, "%Y")
date2=as.Date(paste("2000",m,d),format="%Y%m%d")
temperature=temp$TG/10 temperature[temperature<(-200)]=NA B=aggregate(temperature,by=list(y),FUN= function(x) length(x[!is.na(x)])) if(length(yr)>2){ yr=B[which(B[,2]>250),1] A=aggregate(temperature[y%in%yr],by=list( y[y%in%yr]),FUN=function(x) mean(x,na.rm=TRUE)) A=A[!is.nan(A$x),]
ymin=A[which.min(A[,2]),1]
ymax=A[which.max(A[,2]),1]
plot(date2,temperature,ylim=c(-60,20))
title(paste(fichiers[i],", ",station,sep=""))
lines(date2[y==ymin],temperature[y==ymin],
col="blue",lwd=3)
lines(date2[y==ymax],temperature[y==ymax],
col="red",lwd=3)
legend(date2,20,c(ymax,ymin),
col=c("red","blue"),lwd=3)
}}

Par exemple, si je prends quelques sorties au hasard (toutes celles avec assez d’observations sont du même ordre)      Maintenant, j’ai mis en ligne suffisamment de billets sur mon blog pour croire que personne ne pensera que je suis sceptique quand au réchauffement climatique, et la diminution du réservoir de glace m’inquiète vraiment (on avait travaillé un temps sur les données des régions arctiques, dans un vieux billet we are winter) Mais les écarts de températures mentionnés dans l’introduction m’étonnent. Même si je n’ai pas 2015 dans mes données, je n’observe pas de tels écarts mentionnés (même si des écarts importants seraient normaux, comme les 30 ou 40 degrés d’écart qu’on peut observer à Montréal, mais plutôt l’hiver). Et j’étais d’autant plus surpris que justement, pour l’instant, l’atlantique nord était justement la région qui ne subissait pas de vague de chaleur pour mai, et pour juin (j’attends sur ncdc.noaa.gov/ ou nsstc.uah.edu/climate/ la carte pour juillet 2015).

# 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( "http://climate.weatheroffice.gc.ca/climateData/hourlydata_e.html? timeframe=1&Prov=QC&StationID=5415&hlyRange=1953-01-01|2013-02- 01&Year=",Y,"&Month=",M,"&Day=",D,sep="") > 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… # 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
PENTESTD[i]=summary(R)$coefficients[2,2] } m=lm(T~D)$coefficients
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 for lower left tail and for upper right tail, where is the survival copula associated to , i.e. and  It looks like there is no tail dependence (in the uper tail). But it is also possible to study weaker tail dependence, through and  Slides can be visualized below, I will upload them soon, # You find it cold in Montréal ? trust me, it is even worse As people say in Montréal, “aujourd’hui, il fait frette”. And I have been surprised recently when some people told my that we would reach -35°C Sunday evening… I checked around, and I found -25°C on all weather forecast websites. But nowhere -35°C. I asked some friends, and they told me that those people were not really looking at the air temperature (as we observe on the thermometer), but they were looking at the wind chill, also called “felt air temperature on exposed skin due to the wind” (température ressentie). And indeed, such a quantity does exist, and can be found on theclimate.weatheroffice.gc.ca website. There is also a physical background for that quantity. Hence, the windchill is defined as where is the air temperature (in °C), and the wind speed (in km/h). Please don’t ask me how to interpret this power 0.16 (I already find difficult to explain a square root in an econometric equation). If we look at the past previous days we observe the following observations, where points on top are temperature, while below we have felt temperature.So, basically, winters are even colder than what you might think.. And the story is not over, yet. The same thing holds for summer: if you take into account humidity, summer are even hotter than what you think… There is thehumidex , defined here as where denotes a dewpoint (see here for more details). That index appeared in the 70’s, with a work of Masterson and Richardson entitled a method of quantifying human discomfort due to excessive heat and humidity (published in 1979).By that time, in Canada, on average, 22 people died, per year, because of those excessive heat and humidity. For those interested by the origin of that index, you can have a look here. Recently, @Annmaria (here) told me that one might expect variance to increase, i.e. maximas should be increasing faster than minimas. I just wonder if this intuition can be related to the fact that more and more people (including some medias) now talk more about felt temperatures than measured temperatures. And if we compare past temperatures to felt temperature we have today, it looks like the difference between extremes is increasing…. # 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
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
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…

# Tiens voilà, la pluie. Ah! quel sale temps. Où est-il l’été ? l’été où est-il ?

Dans Les cafards, de Jo Nesbø, Harry Hole (qui vient d’Oslo pour ceux qui n’ont pas dévoré la série) s’entend dire à un moment qu’à Bangkok, les discussions quotidiennes portent rarement sur “la pluie et le beau temps” pour la bonne raison que le climat est assez prévisible à Bangkok (par contre on parle quotidiennement du trafic routier). J’avais évoqué ici l’idée reçue que nous avions, sur le fait que le temps à Rennes n’est pas si changeant que ça. Et je me suis rendu compte depuis que nous sommes à Montréal que le temps pouvait vraiment changer d’un jour sur l’autre (et encore, paraît-il, je n’ai rien vu).
Alors tout d’abord, pour ceux qui en douteraient, la température à Rennes, et à Montréal, ça n’est pas vraiment la même chose, pour Montréal, alors que pour Rennes on obtient (je précise mais je pense qu’on aurait deviné sans) Autrement dit, on a une dispersion presque deux fois plus grande à Montréal. Si regarde également les variations (min/max) dans la journée, on obtient, pour Montréal les courbes suivantes, avec enrouge, la valeur moyenne (sur 15 ans) du maximum observé dans la journée, et enbleu, la valeur moyenne du minimum. A Rennes, on retrouve le fait que la dispersion est relativement faible, En fait, les choses sont encore pire quand on regarde les quantiles, i.e. le pire minimum observé dans 10% des cas (les plus froids) et le pire maximum observé dans 10% des cas (les plus chauds); à  Montréal, alors qu’à Rennes, on est beaucoup plus resserrés, Mais ces pires de cas ne sont pas obtenus dans la même journée, on ne peut pas en dire grand chose. Si on veut aller un peu plus loin, on peut regarder la dispersion dans la journée (i.e. différence entre le minimum et le maximum). On note qu’il est assez stable à Montréal (toujours d’environ 10 degrés) alors qu’à Rennes, la différence est de 5 degrés l’hiver, contre 10 degré l’été, Mais surtout, si on regarde la variation d’un jour sur l’autre, danns 90% des cas, on varie de +/- 8 degrés, sur la température moyenne journalière, avec pas mal de cas un peu extrêmes, avec des variation de 10 ou 15 degré d’un jour sur l’autre sur la température moyenne, alors qu’à Rennes, on a plutôt +/- 5 degrés (avec très peu d’évènements extrêmes, car j’ai tronqué les ordonnées afin de mieux voir au centre) Bref, je ne suis pas prêt de m’arrêter de parler longtemps de la pluie et du beau temps !