# Forecasting Natural Catastrophes (is rather difficult)

Following my previous post, I wanted to spend more time, on the time series with “global weather-related disaster losses as a proportion of global GDP” over the time period 1990-2016 that Roger Pilke sent me last night.

db=data.frame(year=1990:2016, ratio=c(.23,.27,.32,.37,.22,.26,.29,.15,.40,.28,.14,.09,.24,.18,.29,.51,.13,.17,.25,.13,.21,.29,.25,.2,.15,.12,.12))

In my previous post, I spend some time explaining that we should provide some sort of ‘confidence interval’ when we try to predict a pattern. That was what we call ‘model uncertainty’. But there are two (important) issues that I did not mention. (1) it is a time series, so why not use techniques dedicated to time series objects ? (2) we do not really care actually about ‘model uncertainty’ (unless we want to assess if a decreasing trend is significant, or not), and we care more about real prediction uncertainty: in the next ten years, what could be the range for the this ratio, with some given probability (say 95%)? Could we say that with 95% chance the global weather-related disaster losses as a proportion of global GDP should be (each year) within 0 and 0.35 or 0 and 0.7?

A first idea might be to use exponential smoothing techniques (without a seasonal component here).

ratio=ts(db$ratio,start=1990,frequency=1) plot(ratio,xlim=c(1990,2030)) hw=HoltWinters(ratio,gamma=FALSE) phw=predict(hw,n.ahead=15,prediction.interval = TRUE) plot(hw,phw,xlim=c(1990,2030)) polygon(c(2017:2031,rev(2017:2031)), c(phw[,2],rev(phw[,3])),border=NA,col=rgb(0,0,1,.2)) The decreasing trend is coming from the fact that exponential smoothing is here a linear regression, with weight exponentially decaying with time (the older, the smaller the weight). But we cannot use that prediction, since the ratio cannot (obviously) be negative. So why not consider, here, the logarithm of the ratio plot(db$year,db$ratio,type="l",xlim=c(1990,2030),ylim=c(-.2,.7),xlab="year",ylab="ratio") hw=HoltWinters(log(ratio),gamma=FALSE) phw=predict(hw,n.ahead=15,prediction.interval = TRUE) abline(v=2016,lty=2,col="grey") lines(2017:2031,exp(phw[,2]),col="blue") lines(2017:2031,exp(phw[,3]),col="blue") lines(c(1992:2016,2017:2031),c(exp(hw$fitted[,1]),exp(phw[,1])),col="red") polygon(c(2017:2031,rev(2017:2031)),exp(c(phw[,2],rev(phw[,3]))),border=NA,col=rgb(0,0,1,.2)) abline(h=0,lty=2)

The confidence band is huge, here. What if we consider some ARIMA model here?

fit=auto.arima(ratio) farma=forecast(fit,15) farma=cbind(as.numeric(farma$fitted)[1:15],as.numeric(farma$lower[,1]),as.numeric(farma$upper[,1]),as.numeric(farma$lower[,2]),as.numeric(farma$upper[,2])) plot(db$year,db$ratio,type="l",xlim=c(1990,2030),ylim=c(-.2,.7),xlab="year",ylab="ratio") abline(v=2016,lty=2,col="grey") lines(2017:2031,farma[,4],col="blue") lines(2017:2031,farma[,5],col="blue") lines(2017:2031,farma[,1],col="red") polygon(c(2017:2031,rev(2017:2031)),c(farma[,4],rev(farma[,5])),border=NA,col=rgb(0,0,1,.2)) abline(h=0,lty=2) Here, there is an intercept, but no dynamics for the time series (which is considered, here, as a pure white noise). We get exactly the same if we consider the average value of the series fit=lm(ratio~1,data=db) s=summary(fit)$sigma plot(db$year,db$ratio,type="l",xlim=c(1990,2030),ylim=c(-.2,.7),xlab="year",ylab="ratio") abline(v=2016,lty=2,col="grey") ndb=data.frame(year=2017:2031) pf=predict(fit,newdata=ndb) farma=cbind(pf,pf-1.96*s,pf+1.96*s) lines(2017:2031,farma[,2],col="blue") lines(2017:2031,farma[,3],col="blue") lines(1990:2031,c(predict(fit),farma[,1]),col="red") polygon(c(2017:2031,rev(2017:2031)),c(farma[,2],rev(farma[,3])),border=NA,col=rgb(0,0,1,.2)) abline(h=0,lty=2)

Here, we get back to my previous post, if we want to consider a possible trend (and not only an intercept)

fit=lm(ratio~year,data=db) s=summary(fit)$sigma plot(db$year,db$ratio,type="l",xlim=c(1990,2030),ylim=c(-.2,.7),xlab="year",ylab="ratio") abline(v=2016,lty=2,col="grey") ndb=data.frame(year=2017:2031) pf=predict(fit,newdata=ndb) farma=cbind(pf,pf-1.96*s,pf+1.96*s) lines(2017:2031,farma[,2],col="blue") lines(2017:2031,farma[,3],col="blue") lines(1990:2031,c(predict(fit),farma[,1]),col="red") polygon(c(2017:2031,rev(2017:2031)),c(farma[,2],rev(farma[,3])),border=NA,col=rgb(0,0,1,.2)) abline(h=0,lty=2) Again, the confidence region is not based on inference related error, but on model uncertainty: we try to visualize where future observations might be with (say) 95% chance. Note we can also consider (why not?) a quadratic regression fit=lm(ratio~poly(year,2),data=db) s=summary(fit)$sigma plot(db$year,db$ratio,type="l",xlim=c(1990,2030),ylim=c(-.2,.7),xlab="year",ylab="ratio") abline(v=2016,lty=2,col="grey") ndb=data.frame(year=2017:2031) pf=predict(fit,newdata=ndb) farma=cbind(pf,pf-1.96*s,pf+1.96*s) lines(2017:2031,farma[,2],col="blue") lines(2017:2031,farma[,3],col="blue") lines(1990:2031,c(predict(fit),farma[,1]),col="red") polygon(c(2017:2031,rev(2017:2031)),c(farma[,2],rev(farma[,3])),border=NA,col=rgb(0,0,1,.2)) abline(h=0,lty=2)

I am usually not a huge fan of those polynomial regression, but recently, I’ve seen that a lot in economic papers (like “if it’s not linear, add a squared version of the explanatory variable”, which is a rather odd strategy, I’ll publish some posts on that issue this year).

Here again, it might be more clever to consider a logarithmic transformation of the ratio, to insure that the ratio remains positive

fit=lm(log(ratio)~year,data=db) s=summary(fit)$sigma plot(db$year,db$ratio,type="l",xlim=c(1990,2030),ylim=c(-.2,.7),xlab="year",ylab="ratio") abline(v=2016,lty=2,col="grey") ndb=data.frame(year=2017:2031) pf=predict(fit,newdata=ndb) farma=cbind(exp(pf+s^2/2),exp(pf-1.96*s),exp(pf+1.96*s)) lines(2017:2031,farma[,2],col="blue") lines(2017:2031,farma[,3],col="blue") lines(1990:2031,c(exp(predict(fit)+s^2/2),farma[,1]),col="red") polygon(c(2017:2031,rev(2017:2031)),c(farma[,2],rev(farma[,3])),border=NA,col=rgb(0,0,1,.2)) abline(h=0,lty=2) Observe that future trend is mainly driven by the three latest observations, that were rather low (compared with older observations). What if we remove them? dbna=db db$ratio[25:27]=NA fit=lm(ratio~1,data=dbna) s=summary(fit)$sigma plot(db$year,db$ratio,type="l",xlim=c(1990,2030),ylim=c(-.2,.7),xlab="year",ylab="ratio") abline(v=2016-3,lty=2,col="grey") ndb=data.frame(year=2014:2031) pf=predict(fit,newdata=ndb) farma=cbind(pf,pf-1.96*s,pf+1.96*s) lines(2014:2031,farma[,2],col="blue") lines(2014:2031,farma[,3],col="blue") lines(1990:2031,c(predict(fit)[1:24],farma[,1]),col="red") polygon(c(2014:2031,rev(2014:2031)),c(farma[,2],rev(farma[,3])),border=NA,col=rgb(0,0,1,.2)) abline(h=0,lty=2) More funny, if we consider a quadratic regression, we obtain an increasing trend for the future fit=lm(ratio~poly(year,2),data=dbna) s=summary(fit)$sigma plot(db$year,db$ratio,type="l",xlim=c(1990,2030),ylim=c(-.2,.7),xlab="year",ylab="ratio") abline(v=2016-3,lty=2,col="grey") ndb=data.frame(year=2014:2031) pf=predict(fit,newdata=ndb) farma=cbind(pf,pf-1.96*s,pf+1.96*s) lines(2014:2031,farma[,2],col="blue") lines(2014:2031,farma[,3],col="blue") lines(1990:2031,c(predict(fit)[1:24],farma[,1]),col="red") polygon(c(2014:2031,rev(2014:2031)),c(farma[,2],rev(farma[,3])),border=NA,col=rgb(0,0,1,.2)) abline(h=0,lty=2)

As we can see, it is rather difficult to get relevant prediction for the future, based on 25 observations…. If anyone has a suggestion, comments are open…

# 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")

# 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=" ")
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[5],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). # Winter came Generated with pictures from http://nohrsc.noaa.gov/nh_snowcover/ data can be downloaded from ftp://sidads.colorado.edu/DATASETS/NOAA/G02158 See also the Global Sea Ice Area time series, > d=read.table("http://arctic.atmos.uiuc.edu/cryosphere/timeseries.global.anom.1979-2008") > tail(d) V1 V2 V3 V4 12778 2013.984 0.7913005 18.42697 17.63567 12779 2013.986 0.8523080 18.39049 17.53818 12780 2013.989 0.8819072 18.30466 17.42275 12781 2013.992 1.0200537 18.33854 17.31848 12782 2013.995 1.0612829 18.27418 17.21289 12783 2013.997 1.0163171 18.14861 17.13230 > plot(d$V1,d$V3,type="l",ylab="Global Sea Ice Area, million km2") # 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 $\rho(h)=\text{cor}(X_t,X_{t-h})\sim \alpha |h|^{2H-2}$ as $h\rightarrow\infty$ for some $H\in(1/2,1)$ (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 $(1/2,1)$

> 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/
D=as.Date(as.character(TEMP$DATE),"%Y%m%d") T=TEMP$TN/10
day=as.POSIXlt(D)$yday+1 an=trunc(TEMP$DATE/10000)
plot(D,T,col="light blue",xlab="Minimum
daily temperature in Paris",ylab="",cex=.5)
abline(R,lwd=2,col="red")

We can clearly see an increasing linear trend. But we do not care (too much) here about the increase of the average temperature, but more dispersion, and tails. Here are decenal box-plots

or quantile-regressions

library(quantreg)
PENTESTD=PENTE=rep(NA,99)
for(i in 1:99){
R=rq(T~D,tau=i/100)
PENTE[i]=R$coefficients[2] PENTESTD[i]=summary(R)$coefficients[2,2]
}
m=lm(T~D)$coefficients[2] plot((1:99)/100,(PENTE/m-1)*100,type="b") segments((1:99)/100,((PENTE-2*PENTESTD)/m-1)*100, (1:99)/100,((PENTE+2*PENTESTD)/m-1)*100, col="light blue",lwd=3) points((1:99)/100,(PENTE/m-1)*100,type="b") abline(h=0,lty=2,col="red") In order to get a better understanding of the graph above, here are slopes of quantile regressions associated to different probabilities, The annualized maxima (of minimum temperature, i.e. warmest night of the year) i.e. the regression of yearly maximas. tail index of a Generalized Pareto distribution Instead of looking at observation over a century (the trend is obviously linear), we can focus on seaonal behavior, B=data.frame(Y=rep(T,3),X=c(day,day-365,day+365), A=rep(an,3)) library(quantreg) library(splines) Q50=rq(Y~bs(X,10),data=B,tau=.5) Q95=rq(Y~bs(X,10),data=B,tau=.95) Q05=rq(Y~bs(X,10),data=B,tau=.05) YP95=predict(Q95,newdata=data.frame(X=1:366)) YP05=predict(Q05,newdata=data.frame(X=1:366)) I=(T>predict(Q95))|(T<predict(Q05)) YP50=predict(Q50,newdata=data.frame(X=1:366)) plot(day[I],T[I],col="light blue",cex=.5) lines(1:365,YP95[1:365],col="blue") lines(1:365,YP05[1:365],col="blue") lines(1:365,YP50[1:365],col="blue",lwd=3) with on red series from 1900 till 1920, and on purple from 1990 till 2010. If we remove the linear trend, and the seasonal cycle, here are the residuals, assume to be stationary, on during the year Obviously, something has been missed, The graph below is the volatility of the residual series, within the year, Instead of looking at volatility, we can focus on tails, with tail index per month, mois=as.POSIXlt(D)$mon+1
Pmax=Dmax=matrix(NA,12,2)
for(s in 1:12){
X=T3[mois==s]
FIT=gpd(X,5)
Pmax[s,1:2]=FIT$par.ests Dmax[s,1:2]=FIT$par.ses
}
plot(1:12,Pmax[,1],type="b",col="blue",
ylim=c(-.6,0))
segments(1:12,Pmax[,1]+2*Dmax[,1],1:12,Pmax[,1]-
2*Dmax[,1],col="light blue",lwd=2)
points(1:12,Pmax[,1],col="blue")
text(1:12,rep(-.5,12),c("JAN","FEV","MARS",
"AVR","MAI","JUIN","JUIL","AOUT","SEPT",
"OCT","NOV","DEV"),cex=.7)

At the end of the talk, I will also mention multiple city models, e.g. Paris and Marseilles,

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,

# Erratum to the paper on heatwave modelling

An erratum to the paper published on heatwave modelling should appear soon, on http://link.springer.com/article/…

In my article, On the return period of the 2003 heat wave, the sentence “summer of2003 was by far the hottest summer since 1500” is wrongly attributed to ‘Luterbacher et al. (2003)’. My probable intent was to cite Casty et al. (2005) which actually states “the years 1994, 2000, 2002 and particularly 2003 were the warmest since 1500”. Note that a similar statement can also be found in Stott et al. (2004), “The summer of 2003 was probably the hottest in Europe since at latest ad 1500”.

Regarding the statement (in reference to Luterbacher et al., which was actually published in 2004) that “their estimate of the return period of that event is 250 years”, the attribution again is incorrect. Such an estimate can be found in Brown et al. (2005), “best estimate of a 1 in 250 year event” or Lowe et al. (2004) “the same event is likely to occur much more frequently; most likely, once every 250 years”, among others.
• Brown S, Stott P, Clark R, Temperature Extremes (2011) the Past and the Future. Hadley Centre for Climate Prediction and Research Working Paper. [online http://stabilisation.metoffice.com/…]
• Casty C, Wanner H, Luterbacher J, Esper J, Böhm R (2005) Temperature and precipitation variability in the European Alps since 1500. Int J Climatol 25(14):1855–1880, 30 November 2005 CrossRef
• Lowe J, Smith F, Jenkins G, Pope V (2004) Uncertainty, risk and dangerous climate change. Hadley Centre [online http://cedadocs.badc.rl.ac.uk/247/
• Luterbacher J, Dietrich D, Xoplaki E, Grosjean M, Wanner H (2004) European seasonal and annual temperature variability, trends, and extremes since 1500. Science 303:1499–1503 CrossRef
• Stott PA, Stone DA, Allen MR (2004) Human contribution to the European heatwave of 2003. Nature 432:610–614 (2 December 2004) CrossRef

# Trop de probabilités tue les probabilités

Hier matin, au moment de quitter la maison, j’ai jeté un œil à la météo (car à Montréal le temps est très changeant, pas seulement d’un jour sur l’autre, comme je l’avis mentionné en arrivant en septembre, ici). En voyant les gros nuages gris dans le ciel, en me réveillant, je me suis demandé s’il fallait mettre un manteau et des bottes en caoutchouc aux enfants (pour aller l’école). Sur le site de météo, on a en effet la probabilité de pluie, i.e. de P.D.P.

Si on suppose que les probabilités sont indépendantes d’heure en heure, la probabilité qu’il va pleuvoir dans la matinée (dans les 5 prochaines heures pour faire simple, ou plus concret), va s’écrire comme le complémentaire de la probabilité qu’il va faire beau toute la matinée, autrement dit, avec des probabilités de pluie de l’ordre de 25%, par heure, la probabilité qu’il pleuve pendant les prochaines heures est de 75% environ

> 1-((1-.25))^5
[1] 0.7626953

Bref, on est loin de 25%. Essayons d’aller un peu plus loin pour mieux comprendre…
Notons la variable indicatrice indiquant s’il pleut, ou pas (1 s’il pleut, et 0 sinon) à l’instant . Supposons qu’il existe un processus latent sous-jacent, tel que

A une date , on se demande ici s’il va pleuvoir dans la matinée, i.e. on veut

qui peut encore s’écrire

ou, ce qui sera pleut être plus utile ensuite

Avec notre modèle latent, cette probabilité peut s’écrire

c’est à dire en notant la fonction de répartition (jointe) du vecteur aléatoire cette probabilité s’écrit
Or Wassily Hoeffding a montré qu’il était possible d’encadrer une fonction de répartition, avec les bornes suivantes: pour tout $\boldsymbol{x}\in\mathbb{R}^d$,

et

(en notant que ces bornes sont atteignables). Autrement dit, ici, la probabilité qu’il fasse beau 5 heures de suite est encadrée par la probabilité qu’il pleuve sur une heure… et 1.
Bref, on n’a pas grand chose à dire*.

Une possibilité est de supposer un modèle Gaussien pour. Supposer les lois marginales n’est en rien contraignant. En revanche, supposer que la dynamique sous-jacente est Gaussienne est une hypothèse très forte. Et pour faire simple, supposons même que la dynamique soit auto-régressive (à l’ordre 1), i.e.
où $(\eta_t)$ est un bruit blanc, indépendant du processus latent. L’intérêt est que l’on a spécifié ainsi la structure de la matrice de variance covariance. En particulier, la matrice de corrélation du vecteur $\boldsymbol{X}^\star$ est alors

On peut d’ailleurs visualiser facilement la forme de cette matrice en fonction de $\rho$,

Avec ce formalisme, on peut alors en déduire car on connaît très bien la fonction de répartition d’un vecteur gaussien multivarié via

> library(mnormt)

On peut alors étudier l’impact de (qui est associée à la probabilité d’avoir de la pluie pendant une heure en moyenne) et de la corrélation (qui est associée à la dynamique du processus). Pour simplifier, si on prend tel que la probabilité de pluie soit de l’ordre de 25%, on a la probabilité d’avoir de la pluie suivante,

> s = .25
> r = .85
> i = 1:n
> I = matrix(rep(i,each=n),n,n)
> J = matrix(rep(i,n),n,n)
> C = r^(abs(I-J))
> 1-pmnorm(rep(qnorm(1-s),n),mean=rep(0,n),varcov=C)[1]
[1] 0.458074

i.e. en faisant une boucle sur la valeur de la corrélation

Et plus généralement si on croise le seuil (i.e. la probabilité horaire de pleuvoir) et la probabilité, on a l’abaque suivante

On pourrait toutefois objecter que l’important n’est pas forcément qu’il ne pleuve pas de la matinée, mais qu’il ne pleuve pas quand ils sont dehors, i.e. (pour simplifier), à 8 heures, midi et 16 heures (c’est à dire toutes les 4 heures). On s’intéresse alors à

qui est toujours un vecteur Gaussien, mais avec cette fois une structure de corrélation un peu différente,

Dans ce cas, la probabilité de ne pas avoir de plus à ces trois moments de la journée est alors de la forme suivante

La corrélation joue alors un rôle relativement faible car il intervient à la puissance 4 dans la matrice d’auto-corrélation.
Dans ce cas, on notera que l’hypothèse d’indépendance n’est peut être pas absurde. Et avec une P.D.P. de 30%, sur toute la journée, on peut probablement se dire que la probabilité d’avoir de la pluie soit le matin, soit le midi, soit le soir sera de l’ordre de

> 1-(1-.3)^3
[1] 0.657

A tester !
* En fait si, on pourrait dire beaucoup de choses (et écrire une dizaine de billets) en particulier sur la borne supérieure, ici 1, dès lors que la probabilité horaire d’avoir de la pluie dépasse 20% – car la borne supérieure est

En fait, la borne inférieur n’est pas une fonction de répartition de manière générale. Toutefois, la borne est atteignable pour une certaine structure de dépendance. Ce qui explique que par la suite, avec un modèle Gaussien, on sera relativement loin de cette borne supérieure. La borne inférieure (i.e. la borne supérieure sur les fonction de répartition) est elle atteinte dès lors que les risques sont variables sont comonotones, c’est à dire dans le cas Gaussien si la corrélation entre deux instants (et entre tous les instants) vaut 1. La petite animation sur la matrice de corrélation permet de mieux comprendre la difficulté de comprendre le cas où la corrélation rho est négative: dans le cas où 1, on a une corrélation proche de 1 partout dans la matrice. En revanche, si -1, on a des corrélations proches de -1 dans la moitié (environ) des cellules de la matrice, et +1 dans l’autre moitié. La borne supérieure sur la probabilité savoir de la pluie dans matinée correspond donc à une structure plus complexe que la borne inférieure, car on ne peut pas avoir une corrélation négative forte avec tous les lags.

# ma que calor !

J’avais commencé plusieurs analyse (ici ou ) sur les appels passé au 911, à la police, à Montréal, en fonction de l’heure et des jours de la semaine afin de détecter d’éventuels cycles. Mais ces derniers jours, maintenant qu’il commence à faire beau et chaud, je me suis demandé s’il n’était pas possible de croiser avec la température… Car j’aurais l’idée naïve qu’on ne va pas commettre un cambriolage quand il fait -25°C, ou que les conflits sont en hausse s’il fait +35°C…

Si on regarde assez rapidement, en comptant le nombre d’appels en fonction de la température, et en faisant une régression splines, on obtient

autrement dit, entre -20°C et +30°C, il y a environ 50% d’appels en plus. Ce qui est loin d’être négligeable ! Mais on peut aussi croiser la température avec l’heure d’appel, avec l’heure en abscisse et la température en ordonnée

Aussi, la région ou le nombre d’appels est le plus élevé est indiquée ci-dessous: on note que l’on a plus d’appels quand il faut chaud, et plutôt en fin de journée,

Pour les conflits, par exemple, on a davantage de conflits en fin d’après midi comme on l’avait relevé, mais surtout, ces derniers sont très liés à la température: il y a d’autant plus d’appels qu’il faut chaud.

Les appels sont d’autant plus importants dans la région bleue ci-dessous,

Si l’on se concentre uniquement sur la température et les conflits, on a l’évolution suivante

après presque deux fois plus d’appels s’il fait +30°C par rapport à -20°C. De même pour les agressions sexuelles, qui augmentent très fortement avec la température (environ +25%),

En revanche les holds ups sont indépendants de la température,

En revanche les alarmes cambrioleurs semblent générer plus d’appels s’il fait chaud, mais aussi s’il fait froid…

# 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",
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, warming 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",
Dmaxparis=as.Date(as.character(tmaxparis$DATE),"%Y%m%d") Tmaxparis=as.numeric(tmaxparis$TX)/10
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…

# The return period of the 2003 heat wave

The paper on “on the return period of the 2003 heat wave” has been published online here, before Christmas. It should be published soon by Climatic Change, http://link.springer.com/article/...

Extremal events are difficult to model since it is difficult to characterize formally those events. The 2003 heat wave in Europe was not characterized by very high temperatures, but mainly the fact that night temperature were no cool enough for a long period of time. Hence, simulation of several models (either with heavy tailed noise or long range dependence) yield different estimations for the return period of that extremal event.

To go further on the impact on mortality (which was not the aim of the paper), there is a paper in Nature (here).

# 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 !

# 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

 jour beau temps pluie jour beau temps 1955 612 2567 pluie 606 723 1329

ce qui donne les probabilités de transition suivantes,

 jour beau temps pluie jour beau temps 76,15 % 23,85 % pluie 45,60 % 54,40 %

Pour Paris, la probabilité de transition jour après jour a la forme suivante

 jour beau temps pluie jour beau temps 2689 959 3648 pluie 946 1466 2412

ce qui donne les probabilités de transition suivantes,

 jour beau temps pluie jour beau temps 73,71 % 26,29 % pluie 39,22 % 60,78 %

Pour Marseille, la probabilité de transition jour après jour a la forme suivante

 jour beau temps pluie jour beau temps 2527 375 2902 pluie 362 216 578

ce qui donne les probabilités de transition suivantes,

 jour beau temps pluie jour beau temps 87,08 % 12,92 % pluie 62,63 % 37,37 %

Pour Strasbourg, la probabilité de transition jour après jour a la forme suivante

 jour beau temps pluie jour beau temps 31 128 159 pluie 132 1464 1596

ce qui donne les probabilités de transition suivantes,

 jour beau temps pluie jour beau temps 19,50 % 80,50 % pluie 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

 jour beau temps pluie jour beau temps 379 25 404 pluie 26 7 33

ce qui donne les probabilités de transition suivantes,

 jour beau temps pluie jour beau temps 93,81 % 6,19 % pluie 78,79 % 21,21 %

Pour Paris, la probabilité de transition semaine après semaine a la forme suivante

 jour beau temps pluie jour beau temps 576 46 622 pluie 53 4 57

ce qui donne les probabilités de transition suivantes,

 jour beau temps pluie jour beau temps 92,60 % 7,40 % pluie 92,98 % 7,02 %

Pour Marseille, la probabilité de transition semaine après semaine a la forme suivante

 jour beau temps pluie jour beau temps 274 59 333 pluie 47 9 56

ce qui donne les probabilités de transition suivantes,

 jour beau temps pluie jour beau temps 82,28 % 17,72 % pluie 83,93 % 16,07 %

Pour Strasbourg, la probabilité de transition semaine après semaine a la forme suivante

 jour beau temps pluie jour beau temps 1494 614 2018 pluie 613 939 1552

ce qui donne les probabilités de transition suivantes,

 jour beau temps pluie jour beau temps 70,87 % 29,13 % pluie 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….