##REPLICATING JUCKES SI FIGURE 1

##FUNCTIONS TO SMOOTH AND NORMALIZE
  filter.combine.pad<-function(x,a,M=NA) {
        temp<-!is.na(x)  
        year<-c(tsp(x)[1]:tsp(x)[2])[temp]
        w<-x[temp]
        N<-length(w)
        if(is.na(M)) M<-trunc (length(a)/2)+1
        w<-c (rep(mean(w[1:M]),M),w,rep(mean(w[(N-M+1):N]),M))
        y<-filter(w,a,method="convolution",sides=2)
        y<-y[(M+1):(length(y)-M)]
        z<-ts (rep(NA,length(x)),start=tsp(x)[1],end=tsp(x)[2])
        z[temp]<-y
        filter.combine.pad<-ts.union(x,z)
        #filter.combine<-y[,2]
        filter.combine.pad
        }
 f<-function(x){ f<-filter.combine.pad(x,rep(1/21,21))[,2];f}

 normalize<- function(x,M1=1902,M2=1980) {
       	 m1<-mean(x[(M1-tsp(x)[1]+1):(M2-tsp(x)[1]+1)],na.rm=TRUE)
       	 sd1<-sd(x[(M1-tsp(x)[1]+1):(M2-tsp(x)[1]+1)],na.rm=TRUE)
       	 normalize<- ts(  (x-m1)/sd1, start=tsp(x)[1],end=tsp(x)[2])
       	 normalize
       	 }

 sm <- function(x,y){
            l <- length( x )
            z = c( rep( x[1], y ), x, rep(x[l], y) )
            a = filter( z, rep( 1./(2*y+1), (2*y+1) ), method="convolution", sides=2, circular = TRUE )
            a[(1+y):(l+y)]
            }

##LOAD PC EXTRACT
	#first 5 are extracted from mitrie_new_proxy_pcs_1400_v01.nc; 6th is from archived MBH98 BACKTO_1400/pc01.out
	#first 5 are series 16,18,20,22,24 #labeled mbh,mbhl,mbhx,std,cen

	mbhpc<-read.table("http://www.climateaudit.org/data/mitrie/mitrie.pcs.txt",sep="\t")

##EMULATE SI FIGURE 1
	mbhpc<-ts(mbhpc,start=1400)
	mbhpc.scale<-ts( scale(mbhpc),start=1400)

	#try with normalizing on full period

	plot(1400:1980, -sm(mbhpc.scale[,"mbh"],10),lwd=2,type="l",col="red",ylim=c(-2,2.9)) #mbh98
	lines(1400:1980,-sm(mbhpc.scale[,"mbhx"],10),lwd=2,lty=2) #mbhx
  	lines(1400:1980,-sm(mbhpc.scale[,"std"],10),lwd=2,col="blue") #std
   	lines(1400:1980,-sm(mbhpc.scale[,"cen"],10),lwd=2,col="maroon") #cen
  	lines(1400:1980,sm(mbhpc.scale[,"archived"],10),lwd=2,col="green",lty=2) #cen
  	legend(1400,2.7,legend=c("mbh","mbhx","std","cen","archived"),fill=c("red","black","blue","maroon","green"),lty=c(1,2,1,1,2))

	#try with scaling on 1856-1980	
	#m0<-apply(mbhpc[(1856:1980)-1399,],2,mean);sd0<-apply(mbhpc[(1856:1980)-1399,],2,sd)
	##mbhpc.scale<-ts( scale(mbhpc,center=m0,scale=sd0),start=1400)
	##this doesn't work either.