[This article was first published on

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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

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.