Table of contents for SSA

  1. Another Copenhagen Synthesis Report Error?
  2. R code for SSA example
  3. SSA decomposition of HadCRU
  4. Proof of AGW

The code for plotting the non-linear temperature trend, using SSA (singular spectrum analysis) in the figure below is here – ssa-demo. I have made it as turnkey as I can. The steps are:

1. Get and Install package ssa (http://r-forge.r-project.org/projects/ssa/). I had to hand-compile and move the C shared library around so it would find it, not sure why.

2. Run the script below with source(”filename”). Uncomment line indicated after the first time to speed it up. You should get the following replication of the Rahmstorf figure with 11 and 14 embedding period.

rahm11-141

Other things to do

SSA (Singular Spectrum Analysis) is like principle components analysis, and fourier analysis, for time series. The output is a signal decomposition into EOF’s that can be trending or periodic.

The plots produced by Rahmsdorf only use the first EOF, like using the first eigenvector. It seems be a filter that throws out the periodic and white noise components.

So we can look at these EOF’s on the original series as follows:

plot(s$y)

ssa1

We only get two EOF’s out of the original CRU series from 1970-2008. But if we run it on the series with the ‘minimum roughness criterion’ data appended to the end of the series, we get more than 10 EOF’s, of which 6 are shown here. I don’t have a clue why this happens!

plot(s2$y[,1:6])

ssa2

Appending that data to the end (the diagonal line in the first figure) appears to add a whole lot of information to the series, according to SSA.

Here is the code:

library(ssa)
library(zoo)

readCRU<-function(file="http://www.cru.uea.ac.uk/cru/data/temperature/hadcrut3gl.txt" , temps=2:13,plot=T) {
f<-read.table(file,fill=TRUE)
d<-as.vector(t(as.matrix(f[seq(1,length(f[,1]),by=2),temps])))
d<-d[d!= 0]
ts(d,frequency=length(temps),start=1850)
}

mssa<-function(t,em=11,lty=1,lwd=3) {
# Do first ssa then work out the slope in the window
s<-ssa(t,L=em)
f<-s$y[(length(t)-em+1):length(t),1]
x<-1:em
l<-lm(f~x)
a<-((1:em)*l$coefficients[2]+t[length(t)])

# Append points to end 'minimum roughness criterion
t1<-ts(c(t,a),start=start(t))
s2<-ssa(t1,L=em)
s3<-window(s2$y,end=end(t))

# Work out the displacement
x<-window(s2$y[,1],end=end(t))
intercept<-lm(t~x)$coefficients[1]

# Plot the line
lines(t1,lty=lty)
#lines(s$y[,1],col=2,lty=lty,lwd=lwd)
lines(s2$y[,1]+intercept,col=3,lty=lty,lwd=lwd)
lines(s3[,1]+intercept,col=4,lty=lty,lwd=lwd)
browser()
}

#Uncomment this line after first source
TEMP<-readCRU(temps=14)

t<-window(TEMP,start=1970)
plot(t,xlim=c(1970,2020),ylab="Anomaly C")
# Plot 11 window ssa
mssa(t,lty=2)
#Plot 14 window ssa
mssa(t,em=14)

Previous in series Next in series