To reject random walk in climate

December 6, 2012
By

(This article was first published on Wiekvoet, and kindly contributed to R-bloggers)

I read the post The surprisingly weak case for global warming and the rejection; Climate: Misspecified. Based on the first, I wanted to make a post, just to write I agree with the second.

The post features a number of plots like this

For me, one of the noticeable features of this figure is how much the red line does not deviate from middle. Hence, in this post I want to examine how extreme it is in the middle.
The red line, is
# cumsum(changes)
The blue lines are (repeated):

# jumps = sample(changes, 130, replace=T)
# lines(cumsum(c(0,jumps)), col=rgb(0, 0, 1, alpha = .1))

For how close the line is to the centre I use
sd(cumsum(changes))
[1] 12.81393
And the simulated distributions:
simdes <- sapply(1:10000,
   function(x) sd(c(0,cumsum(sample(changes, 130, replace=TRUE)))))
Plot of distributions:
plot(density(simdes))
abline(v=sd(cumsum(changes)),col='red')
All of the simulations have a bigger sd than the original
sum(stdes>sd(cumsum(changes)))
[1] 10000
And that is why I can only say, this random jump, it is not the correct model for these data. 

R code

The first part of the code is (obviously) copied, final lines are mine 
theData = read.table("C:\\Users\\...\\GLB.Ts+dSSTcleaned.txt", 
    header=TRUE,na.strings='**') 
theData <- theData[complete.cases(theData),]

# There has to be a more elegant way to do this
theData$means = rowMeans(aggregate(theData[,c("DJF","MAM","JJA","SON")], by=list(theData$Year), FUN="mean")[,2:5])

# Get a single vector of Year over Year changes
rawChanges = diff(theData$means, 1)

# SD on yearly changes
sd(rawChanges,na.rm=TRUE)

# Subtract off the mean, so that the distribution now has an expectaion of zero
changes = rawChanges - mean(rawChanges)

# Find the total range, 1881 to 2011
(theData$means[131] - theData$means[1])/100

# Year 1 average, year 131 average, difference between them in hundreths
y1a = theData$means[1]/100 + 14
y131a = theData$means[131]/100 + 14
netChange = (y131a - y1a)*100 

# First simulation, with plotting
plot.ts(cumsum(c(0,rawChanges)), col="red", ylim=c(-300,300), lwd=3, xlab="Year", ylab="Temperature anomaly in hundreths of a degrees Celsius")

trials = 1000
finalResults = rep(0,trials)

for(i in 1:trials) {
  jumps = sample(changes, 130, replace=T)
  # Add lines to plot for this, note the "alpha" term for transparency
  lines(cumsum(c(0,jumps)), col=rgb(0, 0, 1, alpha = .1))
  finalResults[i] = sum(jumps)
}
# Re-plot red line again on top, so it's visible again
lines(cumsum(c(0,rawChanges)), col="red", ylim=c(-300,300), lwd=3) 
#
# new lines 
sd(cumsum(changes))
simdes <- sapply(1:10000,function(x) sd(c(0,cumsum(sample(changes, 130, replace=TRUE)))))
plot(density(simdes))
abline(v=sd(cumsum(changes)),col='red')
sum(stdes>sd(cumsum(changes)))

To leave a comment for the author, please follow the link and comment on his blog: Wiekvoet.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.