#################################
##FUNCTIONS REQUIRED FOR NARRATIVE SCRIPT TO PRODUCE SIMULATIONS, FIGURES AND STATISTICS FOR
#			XXXXXXXXXX
#	Hockey Sticks, Principal components and Spurious Significance
#			Stephen McIntyre and Ross McKitrick
#				XXXXXXXX, 
#			doi: XXXXXXXXXXXXXXXXXXXX, 2004
##################################################################

#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

   mannomatic<-function(b) {
	N<-nrow(b)
	d<-array ( rep(NA,N*ncol(b)),dim=c(N,ncol(b)))
	m2<-apply(b[(N-78):N,],2,mean)
    	s2<-apply(b[(N-78):N,],2,sd)
	d<-scale(b,center=m2,scale=s2)
  		#A2. CALCULATE AND DIVIDE BY DETRENDED VARIANCE
	sdprox<-rep(NA,ncol(tree))
    	for (j in 1:ncol(d) ) {
      		X<-c(1:(1980-1901))
      		Y<-d[(N-78):N,j]
      		z<-cbind(X,Y)
      		z<-data.frame(z)
      		fm<-lm(Y~X,data=z)
      		sdprox[j]<-sd(fm$residuals)
      		}
  	d<-scale(d,center=FALSE,scale=sdprox)
	mannomatic<-d
	mannomatic
	}

##2. SIMULATION FUNCTION
   #two alternative methods provided: arima and arfima; only arfima used in simulations here
	#NITER IS parameTER
	#returns eigen1 - mannomatic eigenvalues; eigen 2- princomp eigenvalues and 3 - PC1s mannomatic (hockeysticks)  4- PC1s princomp
	#this is done in chunks to fit using my machine

  localfunction<-function(tree,method2,niter) {
	N<-nrow(tree);
      #CALCULATE ATTRIBUTES OF NETWORK FOR USE IN RED NOISE
	if (method2=="arima")	{Data1<-apply( tree,2,arima.data)}
	
	if(method2=="arima2") {Data1<-apply( tree,2,arima.data2)}

	if (method2=="arfima") {Data<-array(rep(NA,N*ncol(tree)), dim=c(N,ncol(tree)));
		for (k in 1:ncol(tree)){ 
			Data[,k]<-acf(tree[,k][!is.na(tree[,k])],N)[[1]][1:N]
			}#k
		} #arfima	

      #MAKE RED NOISE WITH ATTRIBUTES OF NOAMER SITES AND DO SIMULATION niter TIMES
	n<-ncol(tree)
  	eigen.mannomatic<-array (rep (NA,niter*n), dim=c(n,niter) )
  	eigen.princomp<-array (rep (NA,niter*n), dim=c(n,niter) )
  	PC1.mannomatic<-array (rep (NA,niter*N), dim=c(N,niter) )
  	PC1.princomp<-array (rep (NA,niter*N), dim=c(N,niter) )
  	eof1.mannomatic<-array (rep (NA,niter*n), dim=c(n,niter) )
 
	for (i in 1:niter) {

	#SIMULATE RED NOISE  
	     #arima version (not used here)	
		if (method2=="arima") {b<-array (rep(NA,581*n), dim=c(581,n) )
		for (k in 1:n) {b[,k]<-arima.sim(list(order = c(1,0,0), ar = Data1[k]), n = 581)}
		} #arima bracket


	#alternate arima version 
		if (method2=="arima2") {b<-array (rep(NA,581*n), dim=c(581,n) )
		for (k in 1:n) {b[,k]<-arima.sim(list(order = c(1,0,0), ar = Data1[k]), n = 581)}
		} #arima bracket


	      #arfima version (used here)
		if (method2=="arfima") {N<-nrow(tree);
		b<-array (rep(NA,N*n), dim=c(N,n) )
		for (k in 1:n) {
			b[,k]<-hosking.sim(N,Data[,k])
     			}#k
   		}#arfima

   	#NOW DO PRINCIPAL COMPONENT CALCULATIONS: BOTH VERSIONS - WITH AND WITHOUT MBH98 TRANSFORMATION
   		d<-mannomatic(b) #this does Mannomatic transformation of data 
		w<-svd(d)  
	   	eigen.mannomatic[,i]<-w$d
		PC1.mannomatic[,i]<-w$u[,1]
		eof1.mannomatic[,i]<-w$v[,1]
		z<-princomp(b)  #do calculation without Mannomatic transformation
		eigen.princomp[,i]<-z$sdev
		PC1.princomp[,i]<-z$scores[,1]
	} #end of i-iteration

	#COLLECT AND SAVE RESULTS	
	eigen<-list(eigen.mannomatic,eigen.princomp,PC1.mannomatic,PC1.princomp,eof1.mannomatic)
	names(eigen)<-c("eigen.mannomatic","eigen.princomp","PC1.mannomatic","PC1.princomp","eof1.mannomatic")
	localfunction<-eigen
	localfunction
	} #function


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

#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

  norm<-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)
	norm<- ts(  (x-m1)/sd1, start=tsp(x)[1],end=tsp(x)[2])
	norm
	}


arima.data2<-function(x,lag.max0=250) {
	test<-acf(x,lag.max=lag.max0)
	acf.coef<-c(test$acf)
	temp<-(acf.coef<0)
	N<-min(c(1:(lag.max0+1))[temp])
	z<-data.frame(1:(N-1), acf.coef[1:(N-1)])
	names(z)<-c("lag","acf")
	fm.nls<-nls(acf~rho^lag,data=z,start=list(rho=test$acf[2]))
	arima.data2<-coef(fm.nls)
	arima.data2
	}



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



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