Last year, I did mention in a post that unit-root tests are dangerous, because they might lead us to strange models. For instance, in a post, I did obtain that the temperature observed in January 2013, in Montréal, might be considered as a random walk process (or at leat an integrated process). The code to extract the data has changed (since the website has been updated), so here, we use

library(RCurl) library(XML) options(RCurlOptions = list(useragent = "R")) HEURE=0:23 extracttemp=function(Y,M,D){ url=paste( "http://climate.weather.gc.ca/climateData/hourlydata_e.html?timeframe=1&Prov=QC&StationID=5415&Year=",Y,"&Month=", M,"&Day=",D,sep="") wp <- getURLContent(url) doc <- htmlParse(wp, asText = TRUE) docName(doc) <- url tmp <- readHTMLTable(doc) basejour=data.frame(Year=Y,Month=M,Day=D, Hour=HEURE,Temp=as.numeric(as.character(data.frame(tmp[2])[,2]))[2:25]) return(basejour)} B=NULL for(y in 1955:2013){ for(d in 1:31){ B=rbind(B,extracttemp(y,1,d))}}

Here are all the temperatures observed, and 2013,

plot(B$X,B$Temp,cex=.5,col="light blue",xlab="January, in Montreal",ylab="Temperature (Celsius)") I=which(B$Year==2013) lines(B$X[I],B$Temp[I],col="red")

In the previous post, one test only was used, and one year was considered. I was wondering if this behavior was observed only with temperature of 2013 (or not), and how the other tests (mentioned in a previous post too) were performing.

I might need a function, because those tests cannot be used if there is a missing value, even only one. So I did use the value observed one hour before (just to make sure that the tests can be done)

correcty=function(Y){ I=which(is.na(Y)) if(length(I)==0){Yc=Y} if(length(I)>0){Yc=Y;for(i in I) Yc[i]=Yc[i-1]} return(Yc) }

Now, we can compute the *p*-values, for all the years, and the three different three (keeping in mind that two test if the series is non-stationary, and one if the series is stationary)

DF=matrix(NA,2013-1954,3) library(urca) for(y in 1955:2013){ Z=B$Temp[which(B$Year==y)] Zc=correcty(Z) DF[y-1954,2]=as.numeric(pp.test(Zc)$p.value) DF[y-1954,1]=as.numeric(kpss.test(Zc)$p.value) DF[y-1954,3]=as.numeric(adf.test(Zc)$p.value) }

Visually, if red means stationary, and blue means non-stationary, we get

DFP=DF DFP[,1]=DF[,1]<.05 DFP[,2:3]=DF[,2:3]>.05 library(RColorBrewer) CL=brewer.pal(6, "RdBu") plot(0:1,0:1,xlim=c(1950,2015),ylim=c(0,3),axes=FALSE,xlab="",ylab="") axis(1) text(1952,.5,"KPSS") text(1952,1.5,"PP") text(1952,2.5,"ADF") for(y in 1955:2013){ for(i in 1:3){ polygon(y+c(-1,-1,1,1)/2.2,i-.5+c(-1,1,1,-1)/2.2,col=CL[1+(DFP[y-1954,i]==1)*5],border=NA)}}

Quite frequently, we conclude that the temperature is a random walk. Which does not make sense (from a physical point of view). But again, it might come from the fact that temperature are stationary, but with some fractional behavior (as suggested in the previous post).

This is all very interesting stuff and the code snippets are very helpful. Sadly Climate Canada has reorganised their data access and the URL no longer works. I’ve tried finding the same data but have not managed to fix it.

Could you suggest an update for the url so that code will work now?

Merci beaucoup. 😉