##UTILITY SCRIPT REGARDING VON STORCH AND ZORTIA 2005

################
###FUNCTIONS
################
  sd.detrend<-function(x) {
   	t<-c(1:length(x))
       	fm<-lm(x~t)
      	sd.detrend<-sd(fm$residuals)
	sd.detrend
      	}
	#this is Mann's "detrended" standard deviation

  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
	}	
	#this extends a series by persistence as in MBH98


  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 }
	#this is the MBH98 "standardization", subtracting the short-segment mean and dividing by the short segment
	# detrended standard deviation

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

  hockeystat<-function(test,N=502) { hockeystat<-  ( mean(test[N:(N+79)],na.rm=TRUE)- mean(test[1:(N+79)],na.rm=TRUE) )/sd(test[1:(N+79)],na.rm=TRUE); hockeystat }
	# the hockey stick index as defined in MM05a. 

#################
#LOAD THE NOAMER MATRIX
################
	#this is the NOAMER AD1400 matrix. The results are true for any matrix, but this is a useful and easily accessible example.
	url.source<-"ftp://ftp.agu.org/apend/gl/2004GL021750"
	tree<-read.table(file.path(url.source,"2004GL021750-NOAMER.1400.txt"),header=TRUE,sep="\t") #collated AD1400 network
	tree<-ts(tree[,2:ncol(tree)],start=tree[1,1],end=tree[nrow(tree),1])
	dates<-tsp(tree)# [1] 1400 1992 1
	end0<-min(1980,dates[2])
	tree<-ts( tree[(1400-tsp(tree)[1]+1):(end0-dates[1]+1),] ,start=1400,end=end0) 
	tree<-extend.persist(tree) 
	temp<-!is.na(tree[1,])
	X<-tree[,temp]; rm(tree); dim(X) #581 70
	
#####################
#SHOW THAT THERE IS NO DIFFERENCE BETWEEN THE SVD AND PC IN A CENTERED DATA SET AND THEN THAT THERE IS A 
#DIFFERENCE IN A DECENTERED DATA SET 
######################
	
     	X.centered<-scale(X,center=TRUE,scale=FALSE) #centers X without scaling
	m1<-apply(X[503:581,],2,mean) #calculate short-segment mean
	sd1<-apply(X[503:581,],2,sd) #calculate short-segment sd
    	X.shortcentered<-scale(X,center=m1,scale=FALSE) #short centers X without re-scaling
	
	svd.centered<-svd(X.centered)
	pca.centered<-prcomp(X.centered)
	cor(pca.centered$x[,1],svd.centered$u[,1]) #1
	hockeystat(pca.centered$x[,1]) # 0.2901240
	plot.ts(cbind(pca.centered$x[,1],svd.centered$u[,1]))
		#the PC1s are identical as are other elements

	svd.shortcentered<-svd(X.shortcentered)
	pca.shortcentered<-prcomp(X.shortcentered)
	cor(pca.shortcentered$x[,1],svd.shortcentered$u[,1]) # 0.7285055
	plot.ts(cbind(pca.shortcentered$x[,1],svd.shortcentered$u[,1]))
	apply(cbind(pca.shortcentered$x[,1],svd.shortcentered$u[,1]),2,hockeystat) # 0.2901240 1.3295007
		#the PC1s are not identical
		#the hockey-stickness of the shortcentered version is much greater (but less than with additional use of detrended std dev)
		
####################
#SHOW THAT THERE IS NO DIFFERENCE IN COVARIANCE MATRICES OR PCS MERELY FROM RE-CENTERING
########################

	S<-cov(X.centered)
	S1<-cov(X.shortcentered)
	max(S-S1) # 1.387779e-17
		#the covariance matrices are identical

	pca1<-prcomp(X.centered)
	pca2<-prcomp(X.shortcentered)
	max(pca1$x[,1]-pca2$x[,1]) # 1.332268e-15
		#the PC1s are identical

####################
#SHOW THAT THERE IS A DIFFERENCE BETWEEN THE PC1S FROM STANDARDIZING ON 1400-1980 AND 1902-1980
########################

	X.scaled<-scale(X,center=TRUE,scale=TRUE) #centers and scales on 1400-1980
	X.shortscaled<-scale(X,center=m1,scale=sd1)
	sd2<-apply(X.shortscaled,2,sd)

	pca3<-prcomp(X.scaled)
	pca4<-prcomp(X.shortscaled)
	max(abs(pca3$x[,1])-abs(pca4$x[,1])) # 1.363279
	cor(pca3$x[,1],pca4$x[,1])	#[1] 0.9907951
		#the PC1s are highly correlated, but are not identical
	plot.ts(cbind(pca3$x[,1],pca4$x[,1]))
	apply(cbind(pca3$x[,1],pca4$x[,1]),2,hockeystat) # 0.826877 0.943660
		#there is a very slight additional hockeystickness to the shortscaled version, but not much


####################
#SHOW THAT THERE IS EXTRA HOCKEY STICKNESS WITH "DETRENDED" STD DEV
####################
	X.mannomatic<-apply(X,2,mannomatic)
	svd.mannomatic<-svd(X.mannomatic)
	apply(cbind(svd.shortcentered$u[,1],svd.mannomatic$u[,1]),2,hockeystat)
		# 1.329501 1.618661

		#the hockeystat of 1.62 is what is observed in the archived MBH98 PC1 at Mann's FTP site