#####################
#FUNBCTIONS FOR HUYBERS REPLY 

#
#1. MANN's SHORT SEGMENT STANDARDIZATION PRIOR TO TREE RING PRINCIPAL COMPONENT CALCULATIONS
	#this procedure can be observed in FORTRAN code at Professor Mann's FTP site ftp://holocene.evsc.virginia.edu/pub/MBH98/TREE/ITRDB/NOAMER/pca-noamer

  sd.detrend<-function(x) {
   	t<-c(1:length(x))
       	fm<-lm(x~t)
      	sd.detrend<-sd(fm$residuals)
	sd.detrend
      	}


  mannomatic<-function(x, M=78) {N<-length(x);  xstd<- (x- mean( x[(N-M):N]))/sd(x[(N-M):N]);
	sdprox<-sd.detrend(xstd[(N-M):N]); mannomatic<- xstd/sdprox; mannomatic }

	adjust<-function(x) {
		adjust<-normalize(x)
		adjust<- (adjust*sd.obs) +mean.obs 
		adjust<-ts(adjust,start=tsp(x)[1],end=tsp(x)[2])
		adjust}

#3. FUNCTION TO EXTEND SERIES BY PERSISTENCE
		#=extension by persistence is common practice in MBH98
  extend.persist<-function(tree) {
	extend.persist<-tree
	for (j in 1:ncol(tree) )  {
		test<-is.na(tree[,j])
		end1<-max ( c(1:nrow(tree)) [!test])
		test2<-( c(1:nrow(tree))>end1) & test
		extend.persist[test2,j]<-tree[end1,j]
		}
	extend.persist
	}	
		#JONES LOOKUP
   jones<-function(lat,long) {
	i<-18-floor(lat/5);
	j<- 36+ceiling(long/5);
	jones<-72*(i-1) + j
	jones
	}

#4. FUNCTION TO RETURN VERIFICATION PERIOD STATISTICS
	#calibration and verificaiton are both vectors of length 2 with start and end period
	#estimator and observed are time series which may have different start and end periods

  verification.stats<-function(estimator,observed,calibration,verification)  {
        combine<-ts.union(estimator,observed)
	index.cal<- (calibration[1]-tsp(combine)[1]+1):(calibration[2]-tsp(combine)[1]+1)       
	index.ver<-(verification[1]-tsp(combine)[1]+1):(verification[2]-tsp(combine)[1]+1)
	xmean<-mean ( combine[index.cal,"observed"],na.rm=TRUE  )
       	mx<-mean(  combine[index.ver,"observed"],na.rm=TRUE  ) # 
	test.ver<-cov ( combine[index.ver,],use=use0)
	
	test<-cov(combine[index.cal,],use=use0)
	R2.cal<-(test[1,2]*test[2,1])/(test[1,1]*test[2,2])
	R2.ver<-(test.ver[1,2]*test.ver[2,1])/(test.ver[1,1]*test.ver[2,2])
	RE.cal<-1-sum( (combine[index.cal,"estimator"]- combine[index.cal,"observed"]) ^2,na.rm=TRUE)/sum( (xmean- combine[index.cal,"observed"]) ^2,na.rm=TRUE)
	
	RE.ver<-1-sum( (combine[index.ver,"estimator"]- combine[index.ver,"observed"]) ^2,na.rm=TRUE)/sum( (xmean- combine[index.ver,"observed"]) ^2,na.rm=TRUE)
	CE<- 1-sum( (combine[index.ver,"estimator"]- combine[index.ver,"observed"]) ^2,na.rm=TRUE)/sum( (mx- combine[index.ver,"observed"]) ^2,na.rm=TRUE)
	test<-sign(diff(combine[index.ver,],lag=1))
	temp<-	( test[,1]*test[,2] ==1)
	sign.test<-sum(temp)
	test1<-scale(combine[index.ver,],scale=FALSE)
	temp<-sign(test1)
	temp<-(temp[,1]*temp[,2]==1)
	test2<-test1[,1]*test1[,2]
	prod.test<-  ( mean(test2[temp])+mean(test2[!temp]))/  sqrt( var(test2[temp])/sum(temp)  + var(test2[!temp])/sum(!temp) )
	verification.stats<-data.frame(RE.cal,RE.ver,R2.cal,R2.ver,CE,sign.test,prod.test)
	names(verification.stats)<-c( "RE.cal","RE.ver","R2.cal","R2.ver","CE","sign.test","prod.test")
	verification.stats
	}


#5.	FUNCTION TO NORMALIZE TIME SERIES ON SPECIFIED PERIODS 
	#usually 1902-1980 here

  standardize<-function(x,M1,M2) {
	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)
	standardize<- ts(  (x-m1)/sd1, start=tsp(x)[1],end=tsp(x)[2])
	standardize
	}
  normalize=standardize	

##HAC VARIANCE
	hac<-function(x) {
		fm<-lm(x~1)
		test<-NeweyWest(fm,prewhite=FALSE)
		hac<-test*ssq(x)
		hac<-sqrt(hac)
		hac
		}


#6.  	FUNCTION FOR ARIMA COEFFICIENTS

  arima.data<-function (x)  {
	temp<- !is.na(x)
    	Y<-arima(x[temp ],order=c(1,0,0))
    	arima.data<-Y$coef[1]
  	arima.data
	} #function


#7 	FUNCTION TO Create a directory and subdirectories if they do not exist.
   createsubdirs<-function(dirpath1) {
	if( ! file.exists(dirpath1) ) {
		if( ! file.exists(dirname(dirpath1)) ) {
			# If parent directory is missing, request its creation.
			createsubdirs( dirname(dirpath1) )
		}
		# the parent directory exists but not this dir
		dir.create( dirpath1 )
	   }
	}

#8. 	FUNCTION TO CALCULATE HOCKEYSTAT
hockeystat<-function(test,N=581,M=79) { 
	hockeystat<-  ( mean(test[(N-M+1):N],na.rm=TRUE)- mean(test[1:N],na.rm=TRUE) )/sd(test[1:N],na.rm=TRUE); 
	hockeystat 
	}

#9.	FUNCTION TO CALCULATE CORRELATION TO TEMPERATURE PC1

corfunction<-function(x,N=581,M=79) { 
	corfunction<-cor(x[(N-M+1):N],pc[1:M,1]); 
	corfunction
	}

#FUNCTIONS

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

	mannomatic.extended<-function(b,M=78,sdmethod="detrended") {
        N<-nrow(b)
        m2<-apply(b[(N-M):N,],2,mean)
        s2<-apply(b[(N-M):N,],2,sd)
        if (sdmethod=="unscaled") d<-scale(b,center= m2,scale=FALSE)
        if (sdmethod=="undetrended")  d<-scale(b,center= m2, scale= apply(b[(N-M):N,],2,sd) )   
        if (sdmethod=="detrended")  { d<-scale(b,center= m2, scale= apply(b[(N-M):N,],2,sd) );
                sdprox<-apply(d[(N-M):N,],2,sd.detrend)
                d<-scale(d,center= FALSE,scale= sdprox)}
        mannomatic.extended<-list(d,sdprox)
        mannomatic.extended
        }
	test1<-function(x) {test1<-c(mean(x[(1902-tsp(x)[1]+1):(1980-tsp(x)[1]+1)]),sd(x[(1902-tsp(x)[1]+1):(1980-tsp(x)[1]+1)])); test1}


	covj<-function(x,j) {
		N<-length(x)
		index<-1:(N-j)
		covj<- cov( x[index],x[index+j])
		covj}

	hac.variance<-function(x,method=method.bw) {
		temp<-!is.na(x)
		x<-x[temp]-mean(x[temp])
		if(method=="onethird") L<-trunc(length(x)^1/3 )
		if (method=="sqrt") L<-trunc(length(x)/2)
		if(method=="long") L<-trunc(length(x)/2)
		acf.cov<-acf(x,lag.max=(L-1),type="covariance")$acf	
		index<-c(1:L)
		hac.variance<-  sqrt( var(x)+ 2*sum(  (1 - index/(L+1))*acf.cov)) 	
		hac.variance}
	
	kernel.sun<-function(t,p=8) {kernel.sun<- 1-  abs(t)^p; kernel.sun}

	bartlett.variance<-function(x) {
		x<-x[temp]
		L<-length(x)
		acf.cov<-acf(x,lag.max=(L-1),type="covariance")$acf	
		index<-c(1:L)/L
		bartlett.variance<-  sqrt ( var(x)+ 2*sum (  kernel.sun(index/L)*acf.cov  )  ) 	
		bartlett.variance}

		

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

	mannomatic.extended<-function(b,M=78,sdmethod="detrended") {
        N<-nrow(b)
        m2<-apply(b[(N-M):N,],2,mean)
        s2<-apply(b[(N-M):N,],2,sd)
        if (sdmethod=="unscaled") d<-scale(b,center= m2,scale=FALSE)
        if (sdmethod=="undetrended")  d<-scale(b,center= m2, scale= apply(b[(N-M):N,],2,sd) )   
        if (sdmethod=="detrended")  { d<-scale(b,center= m2, scale= apply(b[(N-M):N,],2,sd) );
                sdprox<-apply(d[(N-M):N,],2,sd.detrend)
                d<-scale(d,center= FALSE,scale= sdprox)}
        mannomatic.extended<-list(d,sdprox)
        mannomatic.extended
        }