library(strucchange)
library(zoo)

readGISS<-function(source="../data/GISS.txt",temps=14,plot=T) {
f<-read.table(source,fill=TRUE,skip=10)
f<-f[f$V1 %in% 1880:2008,]
d<-as.numeric(t(as.matrix(f[seq(1,length(f[,1])),temps])))*0.01
ts(na.omit(d),start=1880,frequency=1)
}

readCRU<-function(source="../data/crutem3vgl.txt",temps=14,plot=T) {
	f<-read.table(source,fill=TRUE)
	d<-as.vector(t(as.matrix(f[seq(1,length(f[,1]),by=2),temps])))
	d<-d[d!= 0]
#d<-c(d,-0,-0)
	ts(d,frequency=length(temps),start=1850)
}




fm<-function(x,l,limit=1) {
	m<-l$coefficients[1]
	m1<-l$coefficients[2]
	m2<-l$coefficients[3]
	s<-summary(l)$coefficients[1,2]
	s1<-summary(l)$coefficients[2,2]
	s2<-summary(l)$coefficients[3,2]
	y1<-(m1+limit*s1)*x+(m2+limit*s2)*(x^2)
	y1
	}
	
fm1<-function(x,l,limit=1) {
	m<-l$coefficients[1]
	m1<-l$coefficients[2]
	s<-summary(l)$coefficients[1,2]
	s1<-summary(l)$coefficients[2,2]
	y1<-m+limit*s+(m1+limit*s1)*x
	y1
	}

TEMP<-readCRU()
GISS<-readGISS()

sea<-read.table("gslGRL2008.txt",skip=14)
SL<-ts(sea$V2,start=1700)
x<-time(SL)
SLt<-window(SL,start=1870)
Tt<-window(TEMP,end=2002)
#SLt<-(SLt-mean(SLt))/sd(SLt)
Tt<-(Tt-mean(Tt))/sd(Tt)
x<-time(SLt)

Gt<-window(GISS,start=1870,end=2002)
#SLt<-(SLt-mean(SLt))/sd(SLt)
Tt<-(Tt-mean(Tt))/sd(Tt)
x<-time(SLt)

sea<-read.table("church_white_grl_gmsl.txt")
cSL<-ts(sea$V2,start=1870,frequency=12)
cSL<-aggregate(cSL, nfrequency = 1, FUN = mean)
#cSL<-(cSL-mean(cSL))/sd(cSL)



trend<-function(y) {
	x<-1:length(y)
	lm(y~x)$coefficients[2]
	}

trendline<-function(T,n=11) {
	l<-c()
	d<-1:(length(T)-n)
	x<-1:n
	for (i in d) {
	l<-c(l,trend(T[i:(i+n-1)]))
	}
	ts(l,start=start(T)+n/2)
	}
	
ssatrend<-function(T,m=11,pad=TRUE) {
	x<-1:m
	y1<-T[(length(T)-m+1):length(T)]
	l1<-lm(y1~x)$coefficients[2]
	y2<-T[1:m]
	l2<-lm(y2~x)$coefficients[2]
	if (pad) {
		Tt<-c(T,T[length(T)]+l1*(1:(1.5*m)))
		}
	else Tt<-T
	s<-c(seq(1:m)/m,1-seq(1:m)/m)/m
	Ttf<-filter(Tt,s)
	Ttf<-na.omit(Ttf)
    ts(Ttf,start=start(T)+m/2)
	}
	
rfilt<-function(T,m=11,filter=TRUE,pad=TRUE,trend=15) {
	l<-c()
	if (filter) Ttf<-ssatrend(T,m,pad)
	else Ttf<-T
	#browser()
	d<-1:(length(Ttf)-trend)
	for (i in d) {
	l<-c(l,trend(Ttf[i:(i+trend)]))
	}
	
	t<-ts(l,start=start(Ttf)+trend/2)

	t
		}
	
run1<-function(cSL,m=15) {
t1<-ssatrend(cSL,m,pad=TRUE)
t2<-ssatrend(cSL,m,pad=FALSE)
plot(t1,xlim=c(1870,2010),ylim=c(-100,150),col=2)
lines(t2-10,col=3)
lines(cSL)
text(1880,150,"Church Pad",col=2)
text(1880,130,"Church NPad",col=3)
text(1880,110,"Church (raw)",col=1)
t1<-ssatrend(SLt,m,pad=TRUE)
t2<-ssatrend(SLt,m,pad=FALSE)
lines(t1,col=2)
lines(t2-10,col=3)
lines(SLt,col=6)
text(1880,90,"Jevrejeva Pad",col=2)
text(1880,70,"Jevrejeva NPad",col=3)
text(1880,50,"Jevrejeva (raw)",col=6)
}

run2<-function(m=15) {
t1<-rfilt(cSL,m,filter=TRUE,pad=TRUE)
t2<-rfilt(cSL,m,filter=TRUE,pad=FALSE)
t3<-rfilt(cSL,m,filter=TRUE,trend=2)
t4<-rfilt(cSL,m,filter=FALSE,trend=30)
plot(t1,col=1,ylim=c(0,4),xlim=c(1880,2010),lwd=3)
lines(t2+0.1,col=2,lwd=3)
lines(t3+0.2,col=3)
lines(t4+0.3,col=4)
text(1885,4,"Filter+Pad",col=1)
text(1885,3.5,"Filter only",col=2)
text(1885,3,"NFilter Trend=2",col=3)
text(1885,2.5,"NFilter Trend=30",col=4)
#browser()
#t<-ssatrend(GISS,m)
#ta<-aggregate(t,nfrequency=1/5,FUN=mean)
#t4a<-aggregate(t4,nfrequency=1/5,FUN=mean)
#plot(ta,t4a)

}

run3<-function(m=15) {
f<-1
t1<-rfilt(cSL,m,filter=TRUE,pad=TRUE)
t2<-rfilt(cSL,m,filter=TRUE,pad=FALSE)
 
t<-ssatrend(Gt,m)
ta<-ts(as.vector(aggregate(t,nfrequency=1/f,FUN=mean)), start=1900, frequency=1/f)

padded<-aggregate(t1,nfrequency=1/f,FUN=mean)
unpadded<-aggregate(t2,nfrequency=1/f,FUN=mean)
par(mfrow=c(2,1),mar=c(2,4,4,3))
plot(ta,padded,col="red",xlim=c(-0.5,0.5),ylim=c(0,3.5))
par(mar=c(4, 4, 2, 3))
plot(ta,unpadded,col="red",xlim=c(-0.5,0.5),ylim=c(0,3.5))
ta0<-window(ta,end=end(padded),start=start(padded))
browser()
l1<-lm(padded~ta0)
ta0<-window(ta,end=end(unpadded),start=start(unpadded))
browser()
l2<-lm(unpadded~ta0)
summary(l1)
summary(l2)
browser()
}