###TRY TO REPLICATE MANN AND JONES 2003 PC1

###
##COLLATE INFORMATION
###

  #1. LOAD "TARGET" PC1 VERSION FROM JONES AND MANN 2004
    #this proxy is identified as series 1 in Table S1 in the O&B SI
    #download version archived at WDCP in connection with Jones and Mann [2004]
 	url="ftp://ftp.ncdc.noaa.gov/pub/data/paleo/contributions_by_author/jones2004/jonesmannrogfig4a.txt"
	g=read.table(url,skip=9,fill=TRUE,nrow=2000)
	dimnames(g)[[2]]=c("year","wusa.mann","jacoby","wusa.briffa","jasper")
	g=ts(g,start=g[1,1],end=g[nrow(g),1])
	temp=(g== -99.99)
	g[temp]=NA
	pc.mj04=ts(g[200:1980,2],start=200) #recovers PC1 as time series

  #2. LOAD MANNOMATIC PC1 for AD200 NETWORK
    #this was previously calculated using script http://www.climateaudit.org/scripts/mann.jones.2003/collation.pc1.txt 
	pc1.mannomatic=read.table("http://www.climateaudit.org/data/mann.jones.2003/pc1.mannomatic.dat",sep="\t",header=TRUE)
	pc1.mannomatic=ts(pc1.mannomatic[,2],start=pc1.mannomatic[1,1])

  #3. LOAD RELEVANT PCs from BACKTO_1000 NETWORK
     #these appear to contribute to splicing - determined by much experimentation
    #these were originally at  Mann's old FTP site (now deleted) holocene.evsc.virginia.edu/MBH98/TREE/ITRDB/NOAMER	
	mbhdata="http://www.climateaudit.org/data"
	pc1f=read.table(file.path(mbhdata,"MBH99/BACKTO_1000-FIXED/pc1-fixed.dat"))
	pc1f=ts(pc1f[,2:4],start=1000);dimnames(pc1f)[[2]]=c("adj","raw","delta")	
		#the three columns are the "adjusted", "raw" and adjustment delta from old MBH archive

   #COLLATE
	pcversions=ts.union(pc.mj04,pc1.mannomatic,pc1f)
	dimnames(pcversions)[[2]]=c("mj04","PC1","pc1f","pc1e","delta")
	use0="pairwise.complete.obs"
	  #these are the Jones and Mann 2004 Archive; Mannomatic Emulation; MBH AD1000 "Fixed"

##MATCHING EFFORTS 
#1. MJ04 Version is rescaling of mannomatic PC1 up to 1700

   cor(pcversions[(200:1700)-199,1:4],use=use0)
	#        mj04    PC1   pc1f   pc1e
	#mj04 1.0000 0.9998 0.8234 0.8234
	#PC1  0.9998 1.0000 0.8217 0.8217
	#pc1f 0.8234 0.8217 1.0000 1.0000
	#pc1e 0.8234 0.8217 1.0000 1.0000

#2. MJ04 Version is rescaling of AD1000 "pc1-fixed"
   cor(pcversions[(1700:1980)-199,],use=use0)
	#        mj04    PC1   pc1f   pc1e  delta
	#mj04  1.0000 0.8087 0.9999 0.9158 0.1837
	#PC1   0.8087 1.0000 0.8092 0.9057 0.5518
	#pc1f  0.9999 0.8092 1.0000 0.9160 0.1839
	#pc1e  0.9158 0.9057 0.9160 1.0000 0.5628
	#delta 0.1837 0.5518 0.1839 0.5628 1.0000

#3. PC1 version archived for MJ04 is standardized over 1751-1950 as stated in Jones and Mann 2004 Figure 4 caption
	#Jones and Mann 2004 Figure 4 Caption:Each series has been normalized over the period 1751–1950 and then smoothed with a 50-year Gaussian filter.

  	m0=apply(pcversions[(1751:1950)-199,],2,mean,na.rm=T);m0 # -0.00035 -0.01489 -0.03702  0.16017  0.19720 
  	sd0=apply(pcversions[(1751:1950)-199,],2,sd,na.rm=T);sd0 # 1.00060 0.01195 0.22818 0.25640 0.10389 

#4. AD1000 pc1-fixed rescaled on 1751-1950 is close match to JM04 archived PC1;
	#but similarly rescaled AD200 PC1 doesn't match
	#ad hoc match by scaling on 1600-1950 and centering on 1000-1650

    X=scale(pcversions[,2:4],center=m0[2:4],scale=sd0[2:4])	
    temp= (time(X)>=1700)	
		#apply(X[(1751:1950)-199,],2,mean,na.rm=T) # 7.367e-17  8.973e-18 -2.524e-17 
		#apply(X[(1751:1950)-199,],2,sd,na.rm=T) #   1    1    1 
 
   nf=layout(array(1:3,dim=c(3,1)),heights=c(1.1,1,1.3))
   par(mar=c(0,3,2,1))	
   	plot(200:1980,pcversions[,1]- X[,2],type="l",xlab="",ylab="",axes=FALSE,ylim=c(-1.5,1.5)) # -0.02721  0.02523  #very close
	axis(side=1,labels=FALSE);axis(side=2,las=1);box();abline(h=0,lty=2)
	text(200,1.4,font=2,pos=4,"To Adjusted AD1000 PC1") 
	title(main="Mann and Jones 2003 PC1")
   	range(pcversions[temp,1]- X[temp,2]) #  -0.02644  0.02489   #very close
   par(mar=c(0,3,0,1))	
   	plot(200:1980,pcversions[,1]- X[,1],type="l",xlab="",ylab="",axes=FALSE,ylim=c(-1.5,1.5)) # -0.02721  0.02523  #very close
	axis(side=1,labels=FALSE);axis(side=2,las=1);box();abline(h=0,lty=2)
	text(200,1.4,font=2,pos=4,"To Rescaled AD200 PC1") 
   par(mar=c(3,3,0,1))	
   	#centered on 1000-1650 or so  #scaled on 1750-1950
	sd1=sd(X[(1600:1950)-199,1]);sd1 # ] 1.048
	m1= mean(X[(1000:1650)-199,1]);m1  -0.6128
	x=scale(X[,1],center=m1,scale=sd1)
	plot(200:1980,pcversions[,1]- x,type="l",xlab="",ylab="",axes=FALSE,ylim=c(-1.5,1.5)) # -0.02721  0.02523  #very close
	axis(side=1);axis(side=2,las=1);box();abline(h=0,lty=2)
	text(200,1.4,font=2,pos=4,"Scale: 1600-1950; Center: 1000-1650") 
	#abline(h= mean(pcversions[index-199,1]- Y[index-199,1]) ,col="grey80")
	c(sd(pcversions[index-199,1]),sd(Y[index-199,1])) #  0.9277 1.0000