Thursday, December 6, 2012

To reject random walk in climate

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

No comments:

Post a Comment