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
sd(cumsum(changes))
 12.81393
And the simulated distributions:
simdes <- sapply(1:10000,
function(x) sd(c(0,cumsum(sample(changes, 130, replace=TRUE)))))
Plot of distributions:
sum(stdes>sd(cumsum(changes)))
 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 <- 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 – theData\$means)/100
# Year 1 average, year 131 average, difference between them in hundreths
y1a = theData\$means/100 + 14
y131a = theData\$means/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)))