To reject random walk in climate

[This article was first published on Wiekvoet, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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
[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:
All of the simulations have a bigger sd than the original
[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”, 
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
# 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 
simdes <- sapply(1:10000,function(x) sd(c(0,cumsum(sample(changes, 130, replace=TRUE)))))

To leave a comment for the author, please follow the link and comment on their blog: Wiekvoet. offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)