R code for SSA example

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

<img src="http://landshape.org/enm/wp-content/uploads/2009/06/ssa2-300×249.png" alt="ssa2" title="ssa2" width="300" height="249"

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&quot; , 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)

Advertisements

0 thoughts on “R code for SSA example

  1. Pingback: wypozyczalnia samochodow Gliwice

  2. Pingback: polecam

  3. Pingback: strona firmy

  4. Pingback: kliknij link

  5. Pingback: depilacja sprawdz

  6. Pingback: pity2015program.pl - pit 2015 format

  7. Pingback: strona www

  8. Pingback: wideofilmowanie Lublin

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s