library(strucchange)
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)
	
	}
	

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