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.Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
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))
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 their blog: Wiekvoet.
R-bloggers.com 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.