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.

## 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)`

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

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)

Pingback: wypozyczalnia samochodow Gliwice

Pingback: polecam

Pingback: strona firmy

Pingback: kliknij link

Pingback: depilacja sprawdz

Pingback: pity2015program.pl - pit 2015 format

Pingback: strona www

Pingback: wideofilmowanie Lublin