Category Archives: Climate

Using “home made” statistics

21/09/2022 Arthur Charpentier Leave a comment

Since I am still at the Fields Institute in Toronto, enjoying a workshop on Impacts of Climate Change on Economics, Finance, and Insurance, I wanted to share some experience, from this summer. After three years of lockdown because of the covid situation, the family has been able to travel, and we went to France, so that our kids could see their grand-parents (the last visit was a long time ago). And it was hot, very hot. While I was chating with my dad, about the weather, and told me that he had a lot of connected devices in the house, including measures of the temperature. One of the device was in a place where nothing did not really change over time. So I thought it could be sufficient to get robust data. My goal was to see how the popular IPCC graph was on real data

When I got the data, I did plot them, and did compare the distribution back in 2012, and in 2022 (or to be honest, half 2021-half 2022). As for the IPCC graph, I assume a Gaussian distribution.

As expected, there is a clear shift to the right (that is “climate change”). But the most scary part, was actually the linear trend,

Coefficients:
              Estimate Std. Error t value Pr(&gt;|t|)    
(Intercept) -637.30455   80.44650  -7.922 3.01e-15 ***
x              0.32273    0.03988   8.092 7.72e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

with a slope of 0.322, meaning that the average temperature is increasing by 0.322 degrees per year ! That is more than 3°C over the past ten years ! Let me write it again : in a house, +3°C on average over the past ten years.

I thought there were some issues with the data. So I tried to collect some official data, and since there were no official records in their village, I did use the data from Lyon (which is 80 kilometers from their house).

The shift on the right is confirmed here, but unfortuntely, I could not get data after 2020. Now

Coefficients:
              Estimate Std. Error t value Pr(&gt;|t|)    
(Intercept) -567.27953   96.26577  -5.893 4.17e-09 ***
x              0.28803    0.04776   6.031 1.81e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

And here again, I have a slope close to 0.3. So again, mainland, about +3°C over the past 10 years was observed. You might not find that scary, but I do think that it is scary !

Climate, Statistics

Forecasting Natural Catastrophes (is rather difficult)

02/01/2017 Arthur Charpentier 8 Comments

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…

Climate

Radial Graphs for Time Series

24/03/2016 Arthur Charpentier 2 Comments

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

Climate

Variation de Température

22/07/2015 Arthur Charpentier 1 Comment

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

Climate

Winter came

27/12/2013 Arthur Charpentier Leave a comment

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

Climate, Computer, Time Series

A random walk ? What else ?

02/02/2013 Arthur Charpentier 7 Comments

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, Computer, Conferences, Graphics, Time Series

Climate change and insurance

14/11/2011 Arthur Charpentier 1 Comment

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[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 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

http://freakonometrics.hypotheses.org/files/2017/07/Llatex2png.2.php_.png

for lower left tail and

for upper right tail, where

http://freakonometrics.hypotheses.org/files/2017/07/toclatex2png-12.2.php_.png

is the survival copula associated to

http://freakonometrics.hypotheses.org/files/2017/07/toclatex2png-13.2.php_.png

, i.e.

http://freakonometrics.hypotheses.org/files/2017/01/toclatex2png-14.2.php_.png

and

http://freakonometrics.hypotheses.org/files/2017/01/toclatex2png-15.2.php_.png

It looks like there is no tail dependence (in the uper tail). But it is also possible to study weaker tail dependence, through

http://freakonometrics.hypotheses.org/files/2017/01/toc2latex2png.3.php_.png

and

http://freakonometrics.hypotheses.org/files/2017/01/toc2latex2png.4.php_.png

Slides can be visualized below, I will upload them soon,

Climate, Publications

Erratum to the paper on heatwave modelling

06/08/2011 Arthur Charpentier Leave a comment

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.

I regret these inappropriate attributions and inadequate referencing.“

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

Climate, Statistics

Trop de probabilités tue les probabilités

04/06/2011 Arthur Charpentier Leave a comment

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 $https://latex.codecogs.com/gif.latex?\boldsymbol{x}\in\mathbb{R}^d$ ,

où

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ù $https://latex.codecogs.com/gif.latex?(\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 $https://latex.codecogs.com/gif.latex?\boldsymbol{X}^\star$ est alors

On peut d’ailleurs visualiser facilement la forme de cette matrice en fonction de $https://latex.codecogs.com/gif.latex?\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

http://freakonometrics.hypotheses.org/files/2015/12/prooooobacommentaire.gif

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.

Climate, Statistics

ma que calor !

01/06/2011 Arthur Charpentier Leave a comment

J’avais commencé plusieurs analyse (ici ou là) 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…

Climate

You find it cold in Montréal ? trust me, it is even worse

25/01/2011 Arthur Charpentier 2 Comments

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

http://freakonometrics.blog.free.fr/public/maths/windchill1.png

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

http://freakonometrics.blog.free.fr/public/maths/humidex.png

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

Climate, Computer

Warming in Paris: minimas versus maximas ?

14/01/2011 Arthur Charpentier 1 Comment

Recently, I received comments (here and on Twitter) about my previous graphs on the temperature in Paris. I mentioned in a comment (there) that studying extremas (and more generally quantiles or interquantile evolution) is not the same as studying the variance. Since I am not a big fan of the variance, let us talk a little bit about extrema behaviour.

In order to study the average temperature it is natural to look at the linear (assuming that it is linear, but I proved that it could reasonably be assumed as linear in the paper) regression, i.e. least square regression, which gives the expected value. But if we care about extremes, or almost extremes, it is natural to look at quantile regression.

For instance, below, the green line is the least square regression, the red one is 97.5% quantile, and the blue on the 2.5% quantile regression.

It looks like the slope is the same, i.e. extremas are increasing as fast as the average…

tmaxparis=read.table("temperature/TG_SOUID100845.txt",
skip=20,sep=",",header=TRUE)
head(tmaxparis)
Dparis=as.Date(as.character(tmaxparis$DATE),"%Y%m%d")
Tparis=as.numeric(tmaxparis$TG)/10
Tparis[Tparis==-999.9]=NA
I=sample(1:length(Tparis),size=5000,replace=FALSE)
plot(Dparis[I],Tparis[I],col="grey")
abline(lm(Tparis~Dparis),col="green")
library(quantreg)
abline(rq(Tparis~Dparis,tau=.025),col="blue")
abline(rq(Tparis~Dparis,tau=.975),col="red")

(here I plot randomly some points to avoid a too heavy figure, since I have too many observations, but I keep all the observations in the regression !).

Now, if we look at the slope for different quantile level (Fig 6 in the paper, here, but on minimum daily temperature, here I look at average daily temperature), the interpretation is different.

s=0
COEF=SD=rep(NA,199)
for(i in seq(.005,.995,by=.005)){
s=s+1
REG=rq(Tparis~Dparis,tau=i)
COEF[s]=REG$coefficients[2]
SD[s]=summary(REG)$coefficients[2,2]
}

with the following graph below,

s=0
plot(seq(.005,.995,by=.005),COEF,type="l",ylim=c(0.00002,.00008))
for(i in seq(.005,.995,by=.005)){
s=s+1
segments(i,COEF[s]-2*SD[s],i,COEF[s]+2*SD[s],col="grey")
}
REG=lm(Tparis~Dparis)
COEFlm=REG$coefficients[2]
SDlm=summary(REG)$coefficients[2,2]
abline(h=COEFlm,col="red")
abline(h=COEFlm-2*SDlm,lty=2,lw=.6,col="red")
abline(h=COEFlm+2*SDlm,lty=2,lw=.6,col="red")

Here, for minimas (quantiles associated to low probabilities, on the left), the trend has a higher slope than the average, so in some sense, 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…

Climate, Computer, Time Series

More climate extremes, or simply global warming ?

12/01/2011 Arthur Charpentier 3 Comments

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

http://freakonometrics.blog.free.fr/public/perso/gev1.png

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…

Climate, Publications

The return period of the 2003 heat wave

02/01/2011 Arthur Charpentier 1 Comment

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

Climate

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

28/09/2010 Arthur Charpentier Leave a comment

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 !

Freakonometrics

Category Archives: Climate

Using “home made” statistics

Forecasting Natural Catastrophes (is rather difficult)

Radial Graphs for Time Series

Variation de Température

Winter came

A random walk ? What else ?

Climate change and insurance

Erratum to the paper on heatwave modelling

Trop de probabilités tue les probabilités

ma que calor !

You find it cold in Montréal ? trust me, it is even worse

Warming in Paris: minimas versus maximas ?

More climate extremes, or simply global warming ?

The return period of the 2003 heat wave

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

"sendo l'intento mio scrivere cosa utile a chi la intende…"