Freakonometrics

To content | To menu | To search

Thursday, March 22 2012

Do we appreciate sunbathing in Spring ?

We are currently experiencing an extremely hot month in Montréal (and more generally in North America). Looking at people having a beer, and starting the first barbecue of the year, I was wondering: if we asked people if global warming was a good or a bad thing, what do you think they will answer ? Wearing a T-shirt in Montréal in March is nice, trust me ! So how can we study, from a quantitative point of view, depending on the time of the year, what people think of global warming ?

A few month ago, I went quickly through

score.sentiment = function(sentences, pos.words,
neg.words, .progress='none')
{
require(plyr)
require(stringr)
scores = laply(sentences, 
function(sentence, pos.words, neg.words) { sentence = gsub('[[:punct:]]', '', sentence) sentence = gsub('[[:cntrl:]]', '', sentence) sentence = gsub('\\d+', '', sentence) sentence = tolower(sentence) word.list = strsplit(sentence, '\\s+') words = unlist(word.list) pos.matches = match(words, pos.words) neg.matches = match(words, neg.words) pos.matches = !is.na(pos.matches) neg.matches = !is.na(neg.matches) score = sum(pos.matches) - sum(neg.matches) return(score) }, pos.words, neg.words, .progress=.progress ) scores.df = data.frame(score=scores, text=sentences) return(scores.df) }

hu.liu.pos = scan("positive-words.txt", what="character",
comment.char=';')
hu.liu.neg = scan('negative-words.txt', what='character',
comment.char=';')

> score.sentiment("It's awesome I am so happy,
thank you all",
+ hu.liu.pos,hu.liu.neg)$score
[1] 3

> score.sentiment("I'm desperate, life is a nightmare,
I want to die",
+ hu.liu.pos,hu.liu.neg)$score
[1] -3

But one can easy see a big problem with this methodology. What if the sentence included negations ? E.g.

> score.sentiment("I'm no longer desperate, life is
not a nightmare anymore I don't want to die",
+ hu.liu.pos,hu.liu.neg)$score
[1] -3

Here the sentence is negative, extremely negative, if we look only at the score. But it should be the opposite. I simple idea is to change (slightly) the function, so that once a negation is found in the sentence, we take the opposite of the score. Hence, we just add at the end of the function

if("not"%in%words){score=-score}

Here we obtain

> score.sentiment.neg("I'm no longer desperate,
life is not a nightmare anymore I don't want to die",
+ hu.liu.pos,hu.liu.neg)$score
[1] 3

But does it really work ? Let us focus on Tweets,

library(twitteR)

Consider the following tweet-extractions, based on two words, a negative word, and the negation of a positive word,

> tweets=searchTwitter('not happy',n=1000)
> NH.text= lapply(tweets, function(t) t$getText() )
> NH.scores = score.sentiment(NH.text,
+ hu.liu.pos,hu.liu.neg)
 
> tweets=searchTwitter('unhappy',n=1000)
> UH.text= lapply(tweets, function(t) t$getText() )
> UH.scores = score.sentiment(UH.text,
+ hu.liu.pos,hu.liu.neg)

> plot(density(NH.scores$score,bw=.8),col="red")
> lines(density(UH.scores$score,bw=.8),col="blue")

> UH.scores = score.sentiment.neg(UH.text,
+ hu.liu.pos,hu.liu.neg)
> NH.scores = score.sentiment.neg(NH.text,
+ hu.liu.pos,hu.liu.neg)
> plot(density(NH.scores$score,bw=.8),col="red") > lines(density(UH.scores$score,bw=.8),col="blue")

> w.tweets=searchTwitter("snow",since= LISTEDATE[k],
+ until= LISTEDATE[k+1],geocode="40,-100,2000mi")
> W.text= lapply(w.tweets, function(t) t$getText() )
> W.scores = score.sentiment.neg(W.text,
+ hu.liu.pos,hu.liu.neg, .progress='text')
> M[k]=mean(W.scores$score)
We obtain here the following score function, over three years, on Twitter,

Well, we have to admit that the pattern is not that obvious. There might me small (local) bump in the winter, but it is not obvious...

Let us get back to the point used to introduce this post. If we study what people "feel" when they mention global warming, let us run the same code, again in North America
> w.tweets=searchTwitter("global warming",since= LISTEDATE[k],
+ until= LISTEDATE[k+1],geocode="40,-100,2000mi")
Actually, I was expecting a nice cycle, with positive scores in Spring, and perhaps negative scores during heat waves, in the middle of the Summer...

What we simply observe is that global warming was related to "negative words" on Twitter a few years ago, but we have reached a neutral position nowadays.

And to be honest, I do not really know how to interpret that: is there a problem with the technique I use (obviously, I use here a very simple scoring function, even after integrating a minor correction to take into consideration negations) or is there something going one that can be interpreted ?

Saturday, June 4 2011

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 http://freakonometrics.blog.free.fr/public/perso3/correlation-pluie03.gif la variable indicatrice indiquant s'il pleut, ou pas (1 s'il pleut, et 0 sinon) à l'instant http://freakonometrics.blog.free.fr/public/perso3/correlation-pluie04.gif. Supposons qu'il existe un processus latent sous-jacent, http://freakonometrics.blog.free.fr/public/perso3/correlation-pluie05.gif tel que

http://freakonometrics.blog.free.fr/public/perso3/correlation-pluie06.gif

A une date http://freakonometrics.blog.free.fr/public/perso3/correlation-pluie04.gif, on se demande ici s'il va pleuvoir dans la matinée, i.e. on veut

http://freakonometrics.blog.free.fr/public/perso3/correlation-pluie07.gif

qui peut encore s'écrire
http://freakonometrics.blog.free.fr/public/perso3/correlation-pluie08.gif

ou, ce qui sera pleut être plus utile ensuite
http://freakonometrics.blog.free.fr/public/perso3/correlation-pluie09.gif

Avec notre modèle latent, cette probabilité peut s'écrire
http://freakonometrics.blog.free.fr/public/perso3/correlation-pluie10.gif

c'est à dire en notant http://freakonometrics.blog.free.fr/public/perso3/correlation-pluie13.gif la fonction de répartition (jointe) du vecteur aléatoire http://freakonometrics.blog.free.fr/public/perso3/correlation-pluie11.gif cette probabilité s'écrit
http://freakonometrics.blog.free.fr/public/perso3/correlation-pluie12.gif
Or Wassily Hoeffding a montré qu'il était possible d'encadrer une fonction de répartition, avec les bornes suivantes: pour tout http://freakonometrics.blog.free.fr/public/perso3/correlation-pluie02.gif,
http://freakonometrics.blog.free.fr/public/perso3/correlation-pluie14.gif


http://freakonometrics.blog.free.fr/public/perso3/correlation-pluie15.gif

et
http://freakonometrics.blog.free.fr/public/perso3/correlation-pluie16.gif

(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 http://freakonometrics.blog.free.fr/public/perso3/correlation-pluie05.gif. 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.
http://freakonometrics.blog.free.fr/public/perso3/correlation-pluie18.gif
http://freakonometrics.blog.free.fr/public/perso3/correlation-pluie19.gif 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 http://freakonometrics.blog.free.fr/public/perso3/correlation-pluie20.gif est alors
http://freakonometrics.blog.free.fr/public/perso3/correlation-pluie01.gif

On peut d'ailleurs visualiser facilement la forme de cette matrice en fonction de http://freakonometrics.blog.free.fr/public/perso3/correlation-pluie21.gif,
http://freakonometrics.blog.free.fr/public/perso3/animationcorrelationpluie.gif
Avec ce formalisme, on peut alors en déduire http://freakonometrics.blog.free.fr/public/perso3/correlation-pluie12.gif car on connaît très bien la fonction de répartition d'un vecteur gaussien multivarié via

On peut alors étudier l'impact de http://freakonometrics.blog.free.fr/public/perso3/slatex.gif (qui est associée à la probabilité d'avoir de la pluie pendant une heure en moyenne) et de la corrélation http://freakonometrics.blog.free.fr/public/perso3/correlation-pluie21.gif (qui est associée à la dynamique du processus). Pour simplifier, si on prend http://freakonometrics.blog.free.fr/public/perso3/slatex.gif 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 à
http://freakonometrics.blog.free.fr/public/perso3/pluiecorrel-31.gif

qui est toujours un vecteur Gaussien, mais avec cette fois une structure de corrélation un peu différente,
http://freakonometrics.blog.free.fr/public/perso3/ppluie-coreel.gif

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.blog.free.fr/public/perso3/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ù http://freakonometrics.blog.free.fr/public/perso3/rhofleche.gif1, on a une corrélation proche de 1 partout dans la matrice. En revanche, si http://freakonometrics.blog.free.fr/public/perso3/rhofleche.gif -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.

Wednesday, February 23 2011

More and more natural catastrophes...

A few weeks ago (here, or there, among many others) I discussed shortly the increase of the frequency of natural catastrophes. An example is perhaps Australia, with major droughts over the past 10 years; a summer heat wave in Victoria, Australia, caused the massive bushfire in 2009; the flooding in Queensland this winter, as well as a huge cyclone. So it looks like Australia is experiencing that increase of the frequency of natural catastrophes.

Nature published last week a series of interesting papers on natural catastrophes (and its relation with a human factor). Increased flood risk linked to global warming by Quirin Schiermeier (likelihood of extreme rainfall may have been doubled by rising greenhouse-gas levels);Climate change: Human influence on rainfall by Richard P. Allan including some Letter by Min and another Letter by Pall; Human contribution to more-intense precipitation extremes by Seung-Ki Min,et al; and Anthropogenic greenhouse gas contribution to flood risk in England and Wales in autumn 2000 by Pardeep Pall et al.


Tuesday, January 25 2011

You find it cold in Montréal ? trust me, it is even worse than what you can imagine...

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 the climate.weatheroffice.gc.ca website. There is also a physical background for that quantity. Hence, the windchill is http://freakonometrics.blog.free.fr/public/maths/windchill2.png defined as
http://freakonometrics.blog.free.fr/public/maths/windchill1.png
where http://freakonometrics.blog.free.fr/public/maths/windchill3.png is the air temperature (in °C), and http://freakonometrics.blog.free.fr/public/maths/windchill4.png 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 the humidex, http://freakonometrics.blog.free.fr/public/maths/humidex2.png, defined here as
http://freakonometrics.blog.free.fr/public/maths/humidex.png
where http://freakonometrics.blog.free.fr/public/maths/humidex3.png 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....

Tuesday, January 18 2011

Catastrophes naturelles et population mondiale

Pour reprendre une question posée au téléphone par François, et par Eric dans un commentaire (ici), il est possible de comparer l'évolution du nombre de catastrophes et de l'évolution de la population (mondiale). Pour les données de population on peut en trouver sur wikipedia (ici), sur le site des nations unies () après 1950, ou encore sur le site du census.gov (). Ensuite on peut moyenner les doublons.

POP=read.table("http://freakonometrics.blog.free.fr/public/data/pop-mondial.txt",
sep=";",header=TRUE)
population=tapply(POP$POP,POP$YEAR,mean)
years=as.numeric(names(population))
reg=lm(population~bs(years,3))
lines(years,predict(reg),col="red")
plot(years,population)
baseline=predict(reg,newdata=data.frame(years=1900:2010))

Revenons à nos catastrophes "naturelles". Par exemple les tremblements de terre. La première piste (évoquée ) était de lisser cette série pour l'utiliser en baseline (multiplicatif). On retrouve ici une série qui oscille autour de 1.
earthquake=read.table("http://freakonometrics.blog.free.fr/public/data/DataSet.csv",
sep=",",header=TRUE)
flood=read.table("http://freakonometrics.blog.free.fr/public/data/DataSet_2_.csv",
sep=",",header=TRUE)
storm=read.table("http://freakonometrics.blog.free.fr/public/data/DataSet_3_.csv",
sep=",",header=TRUE)
library(splines)
n=nrow(earthquake)
earthquake=earthquake[-((n-3):n),]
DATE=1900:2010
EARTHQUAKE=FLOOD=STORM=rep(0,length(DATE))
EARTHQUAKE[as.numeric(as.character(earthquake$Year))-1899]=
as.numeric(as.character(earthquake$Number))
plot(DATE,EARTHQUAKE,type="b")
reg=lm(EARTHQUAKE~bs(DATE,6))
lines(DATE,predict(reg),col="red")
plot(DATE,EARTHQUAKE/predict(reg),type="b")
abline(h=1,lty=2,col="blue")
Si on utilise cette fois la population mondiale comme baseline, on obtient
POP=read.table("http://freakonometrics.blog.free.fr/public/data/pop-mondial.txt",
sep=";",header=TRUE)
population=tapply(POP$POP,POP$YEAR,mean)
years=as.numeric(names(population))
reg=lm(population~bs(years,3))
lines(years,predict(reg),col="red")
plot(years,population)
baseline=predict(reg,newdata=
data.frame(years=1900:2010))/1000000
plot(DATE,EARTHQUAKE/baseline,type="b")

La série de points présente une jolie rupture. Afin de trouver la date optimale de rupture, utilisons un test basé sur le critère d'Akaike.
Y=EARTHQUAKE/baseline
X=DATE
DUMMY=1910:2000
VAIC=rep(NA,length(DUMMY))
for(x in DUMMY){
IX=(X>x)
reg=lm(Y~1+IX)
VAIC[x-1909]=AIC(reg)
}
plot(DUMMY,VAIC)
k=which.min(VAIC)
DUMMY[k]
Manifestement, il y a une rupture en 1975 (que je ne sais trop comment interpréter).
Si on regarde une autre série, comme celle des inondations.
n=nrow(flood)
flood=flood[-((n-3):n),]
DATE=1900:2010
FLOOD[as.numeric(as.character(flood$Year))-1899]=
as.numeric(as.character(flood$Number))
plot(DATE,FLOOD,type="b")
reg=lm(FLOOD~bs(DATE,6))
lines(DATE,predict(reg),col="red")
population=tapply(POP$POP,POP$YEAR,mean)
years=as.numeric(names(population))
reg=lm(population~bs(years,3))
baseline=predict(reg,newdata=data.frame(years=1900:2010))
plot(DATE,FLOOD/baseline,type="b")
On retrouve ici clairement le fait que la croissance de la population ne peut expliquer la croissance du nombre d'inondations.

Saturday, January 15 2011

Increasing number of catastrophes and human factor

Peter asked me (here) to discuss the occurrence of natural catastrophes. I have looked at the EM-DAT dataset (here), and more specifically three kinds of risks: earthquakes, floods and storms (in any country over the past century). For earthquakes, we have the following figures,

the blue points are raw data (I mean, public ones, those that can be downloaded on the website), and I added two smoothed patterns, in red (based on splines regressions). This is no reason to assume that there could be an increase of the number of earthquakes. The trend is due to the fact that we look at earthquakes considered as catastrophic, i.e. that caused deaths and severe damages. And as pointed out in Peter's comment, it means that "the increase in 'disasters' might have a lot to do with the development of the population".
If we look now at storms, we have

i.e. also an increasing trend. But here, we can consider that the trend can be explained by two factors: an inflation due to the value at risk (e.g. increase of concentration of population in some areas as for earthquakes) but maybe also a change that can be related to climate change.  If we remove the pattern du to the increase of value at risk on the raw data (i.e. we divide raw data by the trend observed on earthquakes, the red line), we have

Here it looks like the number of catastrophic storms is increasing (removing the human factor). And it might be possible to relate that trend to climate change. If we look finally at major flood events, we also observe an increase of the pattern,

and we can divide by the trend observed on earthquakes, to obtain

Looking at that graph, it is possible to propose another interpretation about the increase: people might have a shorter memory about rivers than the one they have on earthquakes; and that they think they might influence climate risk, and flood risk, (but not physical ones). I mean that people know where earthquake risk is, and did not really concentrate in those areas during the past century. On the other hand, it is possible to assume that people imagine they can change river design, e.g. building dykes, and though they can control the risk. But obviously they couldn't.

Friday, January 14 2011

Warming in Paris: minimas versus maximas ?

Recently, I received comments (here and on Twitter) about my previous graphs on the temperature in Paris. I mentioned in a comment (there) that studying extremas (and more generally quantiles or interquantile evolution) is not the same as studying the variance. Since I am not a big fan of the variance, let us talk a little bit about extrema behaviour.
In order to study the average temperature it is natural to look at the linear (assuming that it is linear, but I proved that it could reasonably be assumed as linear in the paper) regression, i.e. least square regression, which gives the expected value. But if we care about extremes, or almost extremes, it is natural to look at quantile regression.
For instance, below, the green line is the least square regression, the red one is 97.5% quantile, and the blue on the 2.5% quantile regression.

It looks like the slope is the same, i.e. extremas are increasing as fast as the average...
tmaxparis=read.table("temperature/TG_SOUID100845.txt",
skip=20,sep=",",header=TRUE)
head(tmaxparis)
Dparis=as.Date(as.character(tmaxparis$DATE),"%Y%m%d")
Tparis=as.numeric(tmaxparis$TG)/10
Tparis[Tparis==-999.9]=NA
I=sample(1:length(Tparis),size=5000,replace=FALSE)
plot(Dparis[I],Tparis[I],col="grey")
abline(lm(Tparis~Dparis),col="green")
library(quantreg)
abline(rq(Tparis~Dparis,tau=.025),col="blue")
abline(rq(Tparis~Dparis,tau=.975),col="red")
(here I plot randomly some points to avoid a too heavy figure, since I have too many observations, but I keep all the observations in the regression !).
Now, if we look at the slope for different quantile level (Fig 6 in the paper, here, but on minimum daily temperature, here I look at average daily temperature), the interpretation is different.
s=0
COEF=SD=rep(NA,199)
for(i in seq(.005,.995,by=.005)){
s=s+1
REG=rq(Tparis~Dparis,tau=i)
COEF[s]=REG$coefficients[2]
SD[s]=summary(REG)$coefficients[2,2]
}
with the following graph below,
s=0
plot(seq(.005,.995,by=.005),COEF,type="l",ylim=c(0.00002,.00008))
for(i in seq(.005,.995,by=.005)){
s=s+1
segments(i,COEF[s]-2*SD[s],i,COEF[s]+2*SD[s],col="grey")
}
REG=lm(Tparis~Dparis)
COEFlm=REG$coefficients[2]
SDlm=summary(REG)$coefficients[2,2]
abline(h=COEFlm,col="red")
abline(h=COEFlm-2*SDlm,lty=2,lw=.6,col="red")
abline(h=COEFlm+2*SDlm,lty=2,lw=.6,col="red")

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

Wednesday, January 12 2011

More climate extremes, or simply global warming ?

In the paper on the heat wave in Paris (mentioned here) I discussed changes in the distribution of temperature (and autocorrelation of the time series).
During the workshop on Statistical Methods for Meteorology and Climate Change today (here) I observed that it was still an important question: is climate change affecting only averages, or does it have an impact on extremes ? And since I've seen nice slides to illustrate that question, I decided to play again with my dataset to see what could be said about temperature in Paris.
Recall that data can be downloaded here (daily temperature of the XXth century).
tmaxparis=read.table("/temperature/TX_SOUID100124.txt",
skip=20,sep=",",header=TRUE)
Dmaxparis=as.Date(as.character(tmaxparis$DATE),"%Y%m%d")
Tmaxparis=as.numeric(tmaxparis$TX)/10
tminparis=read.table("/temperature/TN_SOUID100123.txt",
skip=20,sep=",",header=TRUE)
Dminparis=as.Date(as.character(tminparis$DATE),"%Y%m%d")
Tminparis=as.numeric(tminparis$TN)/10
Tminparis[Tminparis==-999.9]=NA
Tmaxparis[Tmaxparis==-999.9]=NA
annee=trunc(tminparis$DATE/10000)
MIN=tapply(Tminparis,annee,min)
plot(unique(annee),MIN,col="blue",ylim=c(-15,40),xlim=c(1900,2000))
abline(lm(MIN~unique(annee)),col="blue")
abline(lm(Tminparis~unique(Dminparis)),col="blue",lty=2)
annee=trunc(tmaxparis$DATE/10000)
MAX=tapply(Tmaxparis,annee,max)
points(unique(annee),MAX,col="red")
abline(lm(MAX~unique(annee)),col="red")
abline(lm(Tmaxparis~unique(Dmaxparis)),col="red",lty=2)
On the plot below, the dots in red are the annual maximum temperatures, while the dots in blue are the annual minimum temperature. The plain line is the regression line (based on the annual max/min), and the dotted lines represent the average maximum/minimum daily temperature (to illustrate the global tendency),

It is also possible to look at annual boxplot, and to focus either on minimas, or on maximas.
annee=trunc(tminparis$DATE/10000)
boxplot(Tminparis~as.factor(annee),ylim=c(-15,10),
xlab="Year",ylab="Temperature",col="blue")
x=boxplot(Tminparis~as.factor(annee),plot=FALSE)
xx=1:length(unique(annee))
points(xx,x$stats[1,],pch=19,col="blue")
abline(lm(x$stats[1,]~xx),col="blue")
annee=trunc(tmaxparis$DATE/10000)
boxplot(Tmaxparis~as.factor(annee),ylim=c(15,40),
xlab="Year",ylab="Temperature",col="red")
x=boxplot(Tmaxparis~as.factor(annee),plot=FALSE)
xx=1:length(unique(annee))
points(xx,x$stats[5,],pch=19,col="red")
abline(lm(x$stats[5,]~xx),col="red")
Plain dots are average temperature below the 5% quantile for minima, or over the 95% quantile for maxima (again with the regression line),

We can observe an increasing trend on the minimas, but not on the maximas !
Finally, an alternative is to remember that we focus on annual maximas and minimas. Thus, Fisher and Tippett theory (mentioned here) can be used. Here, we fit a GEV distribution on a blog of 10 consecutive years. Recall that the GEV distribution is
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 http://freakonometrics.blog.free.fr/public/perso/gev4.png is the following, with on the left the minimas and on the right the maximas,

while the scale parameter http://freakonometrics.blog.free.fr/public/perso/gev3.png is

and finally the shape parameter http://freakonometrics.blog.free.fr/public/perso/gev2.png 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...

Friday, January 7 2011

Evolution of the number of natural catastrophes

Recently (here and there) I mentioned Reinsurers comments on the increase on the number of catastrophes. But actually, there are other sources of information about catastrophes. For instance the http://www.cred.be/, Centre for Research on the Epidemiology of Disasters. As the mention it below, there has been an increase in the number of natural catastrophesFor a disaster to be entered into theso-called EM-database at least one of the following criteria must be fulfilled:

• Ten (10) or more people reported killed.
• Hundred (100) or more people reported affected.
• Declaration of a state of emergency.
• Call for international assistance.

Note that it is possible to download data (here), so I guess I will work further with that database. To be continued...

Friday, December 10 2010

Insurance for a warming planet

I read yesterday an interesting review of the book written by Bjørn Lomborg, which proposes an economic analysis (the subtitle being "comparing costs and benefits") of climate change. Actually, it is a collective book written by the Copenhagen Business School. Martin L. Weitzman (here) wrote a nice review in Nature, a few week ago, precisely untitled 'Insurance for a warming planet'. As Martin Weitzman wrote it, surprisingly, there is no real discussion about uncertain in the book: "key parameters
are approximated by firm values, such as the median or mean, rather t
han a probability distribution. The modelling thus becomes a knobt-widdling exercise in optimizing outcomes". And to go further, "The economics of climate change is mainly about decision-making under extreme uncertainty [...] A striking feature of the economics of climate change is that rare but catastrophic events may have unfathomable costs. Deep uncertainty about the unknown unknowns of what might go wrong is therefore coupled with essentially unlimited liability. The resulting battle between declining probabilities and increasing damages is difficult to resolve. Alas, this uncertainty can figure prominently in evaluations of climate change policies. Its absence in a book dealing with economic comparisons of smart solutions is a serious omission [...] When confronted with the possibility of extreme damages at low probabilities, most people do not look to averages. Instead they think about how much insurance they need, and can afford to buy, to survive those events. Climate policy is better viewed as buying insurance for the planet against extreme outcomes than as the solution to a multivariate problem over which we have control".

Sunday, November 7 2010

Updating meteorological forecasts, part 1

As Mark Twain said "the art of prophecy is very difficult, especially about the future" (well, actually I am not sure Mark Twain was the  first one to say so, but if you're interested by that sentence, you can look here). I have been rather surprised to see how Canadians can be interested in weather, and weather forecasts (see e.g. here for a first study). For instance, on that website (here), we can have forecasts of the temperature (but also rain and wind, which is interesting when you ride a bike) over the next week, almost on an hourly basis (I show here only GFS forecasts (Global Forecast System,here), which is quite standard, see here or there, but other meteorological models are considered on the website)

As a statistician (or maybe more as a bike rider) I was wondering if I should trust those forecasts... and if a mistake today means an error tomorrow, too...
So, to look at this problem, Vincent (alias @Vicnent) helped me (again, since he helped me already to get backups of betting websites during the soccer world cup, in June) to make backups of those tables...
  • A descriptive study
For instance, if we look at forecasts made on the 21st of October (00:00) we have the following for the temperature in °C, for different time horizon,
6 hours latter (21st of October (06:00)), we have the following

On the other hand, consider forecasts made for the 25th of October (11:00) we have

while for 3 hours later (25th of October (14:00)) we have

and again for 3 hours later (25th of October (17:00)) we have

So obviously, 5 days earlier, there has been a bump of optimism, or perhaps simply a typing mistake since we jump from 4 before and after to 14 during 6 hours*. If we can expect to have different predictions for different days (even within the day since it is usually cooler during the night), we can try to see if bumps, or forecast errors can be understood.
Here, it is difficult to observe the scatterplot of temperature,

That graph was based on the following points, with here the scatterplot of prediction made on the 21st of October (00:00)

with a smoothed version below (using standard splines),

or, if we consider forecasts made for the 25th of October (11:00) we have

with again, a smoothed version,

The code is the  following,
donnees=read.table("http://perso.univ-rennes1.fr/
arthur.charpentier/pred-meteo2.txt"
)
x=donnees$DATE1+donnees$HEURE1/24
x[x<15]=x[x<15]+31
y=donnees$DATE2+donnees$HEURE2/24
y[y<15]=y[y<15]+31
z=y-x
h=donnees$TEMPERATURE
D=data.frame(x,y,z,h)
X0=sort(unique(x))
Y0=sort(unique(y))
Z0=sort(unique(z))
matY.X=matrix(NA,length(X0),length(Y0))
library(splines)
for(i in 1:length(X0)){
dd=D[D$x==X0[i],]
ax=dd$y
ay=dd$h
a=data.frame(ax,ay)
ra=lm(ay~bs(ax,df=5),data=ba)
X=Y0[(Y0>min(ax))&(Y0<max(ax))]
iX=which((Y0>min(ax))&(Y0<max(ax)))
Y=predict(ra,newdata=data.frame(ax=X))
matY.X[i,iX]=Y
}
So, let us look at levels of iso-temperature, as a function of the date the prediction is made (the x-axis) and the horizon (the y-axis). First, we we fix the date of the prediction, we get (from the code above)

and then, when we fix the horizon of the prediction

It looks like the small heat wave we experienced at the end of October has been perfectly predicted... And here, it looks like the only pattern we can observe is the standard evolution of temperature: there are no bumps due to forecasting problems..
  • Modeling predicted temperatures
Let us use some additive or multiplicative models, i.e. either
http://perso.univ-rennes1.fr/arthur.charpentier/latex/predictemp01.png
where http://perso.univ-rennes1.fr/arthur.charpentier/latex/predictemp0x.png denotes the date the prediction was made and http://perso.univ-rennes1.fr/arthur.charpentier/latex/predictemp0y.png the date the horizon, or
http://perso.univ-rennes1.fr/arthur.charpentier/latex/predictemp02.png
This was easy with the gam framework... For instance, consider the following code (some regressions mentioned below will be used later on)
library(mgcv)
library(splines)
REG1=gam(h~s(x,y),data=D,family=gaussian)
REG2=lm(h~bs(x)+bs(y),data=D)
REG3=gam(h~s(x)+s(y),data=D,family=gaussian)
REG4=gam(h~s(x)+s(y),data=D,family=gaussian(link="log"))
REG5=gam(h~s(y),data=D,family=gaussian)
D$e=residuals(REG5)
REG6=gam(e~s(x,y),data=D,family=gaussian)
Here is the prediction we had with the additive model
and for the multiplicative model

(colors are different because of the log scaling but the shape is similar).

For the component on http://perso.univ-rennes1.fr/arthur.charpentier/latex/predictemp0x.png (in blue), we have, with the additive model

and for the multiplicative model

On the other hand, for the component on http://perso.univ-rennes1.fr/arthur.charpentier/latex/predictemp0y.png (in red) we have, with the additive model

and for the multiplicative model 

 Note that we observe something which looks like the true temperature (we actually had)

So finally, on average, those predictions are consistent with the true temperature we experienced. But so far, we cannot say anything about prediction errors. If we consider a bivariate spline regression, i.e.
http://perso.univ-rennes1.fr/arthur.charpentier/latex/predictemp04.png
for some smoothed function, we have, again, the same kind of predictions,

It looks like what we've seen before. And because we have those horizontal lines, only the horizon http://perso.univ-rennes1.fr/arthur.charpentier/latex/predictemp0y.png should be considered.... So let us study the following model
http://perso.univ-rennes1.fr/arthur.charpentier/latex/predictemp05.png
where we obtain the following prediction

(which is similar with previous graph) but if we fit a spline on residuals obtained from that regression, i.e.
http://perso.univ-rennes1.fr/arthur.charpentier/latex/predictemp07.png
where
http://perso.univ-rennes1.fr/arthur.charpentier/latex/predictemp06.png
we have the following distribution for http://perso.univ-rennes1.fr/arthur.charpentier/latex/predictemp0h.png

This is where we can see prediction errors. We can observe that meteorologists have been a bit surprised by the heat-wave we observed just after the 25th. But let us wait a few more days to get more data, and perhaps we'll have a look at something that I do care about: prediction of rainfall....
* perhaps I should mention that we do not have backups of tables, but backups of html pages... then I use simple linguistic tools available on R to extract the information where it should be....

Tuesday, September 28 2010

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 en rouge, la valeur moyenne (sur 15 ans) du maximum observé dans la journée, et en bleu, 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 !

Saturday, August 21 2010

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 cm 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 plu 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....

En Bretagne, il fait beau plusieurs fois par jour

bon, et à l'occasion il peut pleuvoir un peu.... Il peut donc être intéressant pour planifier un peu ses vacances de calculer la probabilité d'avoir de la pluie.

Les calculs ont été fait sur les données de pluviométrie à Rennes, en ligne sur le site eca&d, ici (données de qualité gratuites).
  • Probabilité d'avoir de la pluie, pendant les vacances d'été
Considérons ici la série du niveau de précipitation par jour, en par 0,1 mm,

et la probabilité d'avoir de la pluie dans la journée (régression logistique),

Bref, il peut pleuvoir à Rennes l'été. Mais ça n'aide pas vraiment pour planifier ses vacances. Car si ça se trouve, il y a des étés sans pluie, et des étés où il pleut sans cesse.
  • Dynamique et matrice de transition (par jour et par semaine)
Autrement dit, au lieu de regarder les lois marginales (comme la probabilité d'avoir de la pluie dans la journée), on peut s'intéresser à la dynamique de la série, modélisée sous la forme d'une chaîne de Markov. Si on regarde jour après jour, avec les 30 mois de juillet et août entre 1980 et 2009,


jour 


beau
temps
pluie

jour 
beau temps
8712981169
pluie
292335627
ce qui donne les probabilités de transition suivantes,


jour 


beau
temps
pluie
jour 
beau temps
75,51 %25,49 %
pluie
46,57 %53, 43 %
autrement, s'il fait beau aujourd'hui, on a 3 chances sur 4 d'avoir du beau temps demain.
Si on regarde semaine après semaine, où l'intérêt sont les semaines sans pluie,


semaine 


beau
temps
pluie

semaine 
beau temps
132639
pluie
23140163
avec là aussi les probabilités de transition suivantes,


semaine 


beau
temps
pluie
semaine 
beau temps
33,33 %66,67 %
pluie
14,11 %85,89 %
Si on regarde semaine après semaine, où l'intérêt sont les semaines avec six jours sans pluie (on s'autorise une journée de pluie),


semaine 


beau
temps
pluie

semaine 
beau temps
363672
pluie
30100130
ce qui donne les probabilités de transition suivantes,


semaine 


beau
temps
pluie
semaine 
beau temps
50,00 %50,00 %
pluie
23,08 %76,92 %
Bon, maintenant, si on pense qu'avoir 2 mm de pluie dans la journée, c'est juste un peu d'humidité dans l'air, les matrices de transitions sont sensiblement différentes, dans le cas où on s'autorise une journée dans la semaine avec un peu d'humidité pour parler de beau temps,


semaine 


beau
temps
pluie

semaine 
beau temps
106 36 142
pluie
31 29 60
ce qui donne les probabilités de transition suivantes,


semaine 


beau
temps
pluie
semaine 
beau temps
75,65 % 25,35 %
pluie
51,67 % 48,33 %
Autrement dit, on retrouve une matrice proche de celle que nous avions sur les données journalière: s'il a fait beau la semaine avant de venir en Bretagne, on a 3 chances sur 4 d'avoir du beau temps. En revanche, s'il a fait mauvais, on a une chance sur deux d'avoir beau temps la semaine suivante.
Si notre définition de pluie est encore plus laxiste (il faut qu'il y ait eu un déluge une journée dans la semaine, à savoir plus de 2 cm d'eau dans la journée), alors cette fois, on obtient,


semaine 


beau
temps
pluie

semaine 
beau temps
171 12 183
pluie
14 5 130
ce qui donne les probabilités de transition suivantes,


semaine 


beau
temps
pluie
semaine 
beau temps
93,44 % 6,56 %
pluie
73,68 % 26,32 %
Autrement dit, s'il y a eu au moins une très mauvaise journée la semaine passée, on a 1 chance sur 4 d'en avoir également une la semaine suivante. En revanche, s'il a fait beau tout le temps, on a presque 95% de chances d'avoir du beau temps la semaine suivante.
Moralité, quelle que soit la définition retenue pour définir le beau temps, le temps à Rennes n'est pas indépendant d'une semaine sur l'autre: il y a manifestement des étés pluvieux, et d'autre non.