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

sm2<- function(x){ sm2<- filter( x, filter=rep(1/21,21) ,method="convolution", sides=2);sm2}

ssq<-function(x){ ssq<-sum(x^2);ssq}

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

	url.source<-"c:/climate/data/new"
	pc.archived<-read.table(file.path(url.source,"TREE/ITRDB/NOAMER/BACKTO_1400/pc01.out"))
	pc.archived<-ts(pc.archived[,2],start=1400)

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

	plot(1400:1980, -sm2(mbhpc.scale[,"mbh"]),lwd=2,type="l",col="red")#,ylim=ylim0,ylab="scale 1856-1980") #mbh98
	lines(1400:1980,-sm2(mbhpc.scale[,"mbhx"]),lwd=2,lty=2) #mbhx
  	lines(1400:1980,-sm2(mbhpc.scale[,"std"]),lwd=2,col="blue") #std
   	lines(1400:1980,-sm2(mbhpc.scale[,"cen"]),lwd=2,col="maroon") #cen
  	lines(1400:1980,sm2(pc.archived.scale),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))



##Re #13: No problem, sorry I missed out two details: the plot in the supplement is made by normalising the uncentred PCs (this only affects “MBH”, “MBHX”, and “Archived”) whereas your code centres first and then normalises. I’ve checked doing it this way round in my code and I’m fairly sure that it explains the difference in amplitude of variation in the plots. As in the manuscript, plots are centred on the calibration period — because the proxies should be centred in this way before input into the regression algorithm.
##http://www.climateaudit.org/?p=893#comment-62708

	mbhpc.scale<-cbind(mbhpc[,c("mbh","mbhx","std","cen")],pc.archived)
	sd0<-apply(mbhpc.scale,2,sd)
	#mbhpc.scale[,c(1,2,5)] <- scale(mbhpc.scale[,c(1,2,5)],center=FALSE,scale=sd0[c(1,2,5)])
	mbhpc.scale<- scale(mbhpc.scale,center=FALSE,scale=sd0)
	m0<-apply(mbhpc.scale[(1856:1980)-1399,],2,mean)
	mbhpc.scale<- scale(mbhpc.scale,center=m0,scale=FALSE)
	mbhpc.scale<-ts(mbhpc.scale,start=1400)
	dimnames(mbhpc.scale)[[2]]<-c("mbh","mbhx","std","cen","archived")

	par(mar=c(3,3,1,1))
	plot(1400:1980, -sm2(mbhpc.scale[,"mbh"]),lwd=2,type="l",col="red",ylim=c(-3,2),xlab="",ylab="") #mbh98
	lines(1400:1980,-sm2(mbhpc.scale[,"mbhx"]),lwd=2,lty=2) #mbhx
  	lines(1400:1980,-sm2(mbhpc.scale[,"std"]),lwd=2,col="blue") #std
   	lines(1400:1980,-sm2(mbhpc.scale[,"cen"]),lwd=2,col="maroon") #cen
  	lines(1400:1980,sm2(mbhpc.scale[,"archived"]),lwd=2,col="green",lty=2) #cen
  	legend(1400,2.1,legend=c("mbh","mbhx","std","cen","archived"),fill=c("red","black","blue","maroon","green"),lty=c(1,2,1,1,2))



#RMS normalization
	mbhpc.scale<-cbind(mbhpc[,c("mbh","mbhx","std","cen")],pc.archived)
	s0<-(apply(mbhpc.scale,2,ssq)/581)^.5
	sd0<-apply(mbhpc.scale,2,sd)
	#mbhpc.scale[,c(1,2,5)] <- scale(mbhpc.scale[,c(1,2,5)],center=FALSE,scale=sd0[c(1,2,5)])
	mbhpc.scale<- scale(mbhpc.scale,center=FALSE,scale=s0)
	m0<-apply(mbhpc.scale[(1856:1980)-1399,],2,mean)
	mbhpc.scale<- scale(mbhpc.scale,center=m0,scale=FALSE)
	mbhpc.scale<-ts(mbhpc.scale,start=1400)
	dimnames(mbhpc.scale)[[2]]<-c("mbh","mbhx","std","cen","archived")

	par(mar=c(3,3,1,1))
	plot(1400:1980, -sm2(mbhpc.scale[,"mbh"]),lwd=2,type="l",col="red",ylim=c(-3,2),xlab="",ylab="") #mbh98
	lines(1400:1980,-sm2(mbhpc.scale[,"mbhx"]),lwd=2,lty=2) #mbhx
  	lines(1400:1980,-sm2(mbhpc.scale[,"std"]),lwd=2,col="blue") #std
   	lines(1400:1980,-sm2(mbhpc.scale[,"cen"]),lwd=2,col="maroon") #cen
  	lines(1400:1980,sm2(mbhpc.scale[,"archived"]),lwd=2,col="green",lty=2) #cen
  	legend(1400,2.1,legend=c("mbh","mbhx","std","cen","archived"),fill=c("red","black","blue","maroon","green"),lty=c(1,2,1,1,2))



##COMBINED VERSION

	mbhpc.scale<-cbind(mbhpc[,c("mbh","mbhx","std","cen")],pc.archived)
	sd0<-apply(mbhpc.scale,2,sd)
	#mbhpc.scale[,c(1,2,5)] <- scale(mbhpc.scale[,c(1,2,5)],center=FALSE,scale=sd0[c(1,2,5)])
	mbhpc.scale<- scale(mbhpc.scale,center=FALSE,scale=sd0)
	m0<-apply(mbhpc.scale[(1856:1980)-1399,],2,mean)
	mbhpc.scale<- scale(mbhpc.scale,center=m0,scale=FALSE)
	mbhpc.scale<-ts(mbhpc.scale,start=1400)
	dimnames(mbhpc.scale)[[2]]<-c("mbh","mbhx","std","cen","archived")

	par(mar=c(3,3,1,1))
	plot(1400:1980, -sm2(mbhpc.scale[,"mbh"]),lwd=2,type="l",col="blue",ylim=c(-3,2),xlab="",ylab="") #mbh98
	#lines(1400:1980,-sm2(mbhpc.scale[,"mbhx"]),lwd=2,lty=2) #mbhx
  	lines(1400:1980,-sm2(mbhpc.scale[,"std"]),lwd=2,col="maroon") #std
   	#lines(1400:1980,-sm2(mbhpc.scale[,"cen"]),lwd=2,col="maroon") #cen
  	#lines(1400:1980,sm2(mbhpc.scale[,"archived"]),lwd=2,col="green",lty=2) #cen
  	#legend(1400,2.1,legend=c("mbh","mbhx","std","cen","archived"),fill=c("red","black","blue","maroon","green"),lty=c(1,2,1,1,2))

	mbhpc.scale<-cbind(mbhpc[,c("mbh","mbhx","std","cen")],pc.archived)
	s0<-(apply(mbhpc.scale,2,ssq)/581)^.5
	sd0<-apply(mbhpc.scale,2,sd)
	#mbhpc.scale[,c(1,2,5)] <- scale(mbhpc.scale[,c(1,2,5)],center=FALSE,scale=sd0[c(1,2,5)])
	mbhpc.scale<- scale(mbhpc.scale,center=FALSE,scale=s0)
	m0<-apply(mbhpc.scale[(1856:1980)-1399,],2,mean)
	mbhpc.scale<- scale(mbhpc.scale,center=m0,scale=FALSE)
	mbhpc.scale<-ts(mbhpc.scale,start=1400)
	dimnames(mbhpc.scale)[[2]]<-c("mbh","mbhx","std","cen","archived")

	lines(1400:1980, -sm2(mbhpc.scale[,"mbh"]),lwd=2,type="l",col="red",ylim=c(-3,2),xlab="",ylab="") #mbh98
	#lines(1400:1980,-sm2(mbhpc.scale[,"mbhx"]),lwd=2,lty=2) #mbhx
  	#lines(1400:1980,-sm2(mbhpc.scale[,"std"]),lwd=2,col="blue") #std
   	#lines(1400:1980,-sm2(mbhpc.scale[,"cen"]),lwd=2,col="maroon") #cen
  	lines(1400:1980,sm2(mbhpc.scale[,"archived"]),lwd=2,col="green",lty=2) #cen
  	legend(1400,2.1,legend=c("mbh","mbhx","std","cen","archived"),fill=c("red","black","blue","maroon","green"),lty=c(1,2,1,1,2))