#PREISENDORFER ON TREE RING NETWORKS
##THIS SHOWS FUNCTIONS 

###########################
##FUNCTIONS
  #this function calculates eigenvalues for tree ring network plus Rule N 95% eigenvalues for random matrix 
  #default arima method is AR1; ARMA(1,1) is also an option 
  #default method is mannomatic; correlation and covariance are also available
  #default is 100 simulations

 #the AR1 coefficient is calculated for each site and used to generate one column in the network
 #the eigenvalues from 100 simulations are collated and the 95th percentile for each eigenvalue chosen 

  preis_emulation=function(network,step,method="mannomatic",arima_method="ar1",mbhdata="http://www.climateaudit.org/data/mbh98/UVA",K=100) {
	  loc=file.path(mbhdata,paste(network,"dat",sep="."))
	  tree=read.table( loc ,sep="\t",header=TRUE)
	  tree<-ts(tree[,2:ncol(tree)],start=tree[1,1],end=tree[nrow(tree),1])
	  tree=extend.persist(tree)
	  end0=min(1980,tsp(tree)[2])
	  tree=ts(tree[(step:end0)-tsp(tree)[1]+1,],start=step)
	  tree=tree[,!is.na(tree[1,])]
	  if (arima_method=="ar1") {
		g=function(x) arima(x,order=c(1,0,0))$coef[1]
	  	ar0=apply(tree,2,g) }
	  if (arima_method=="arma1") {
		g=function(x) arima(x,order=c(1,0,1))$coef[1:2]
		ar0=apply(tree,2,g) } 
	  lambda=preisf(tree,method)  # normalized eigenvalues
      #do K simulations (base =100)
	N=nrow(tree);M= ncol(tree)
	stat=array(NA,dim=c(K,M))
	for (k in 1:K) {
	  X=array(NA,dim=c(N,M))
	  if(arima_method=="arma1") for (j in 1:M) X[,j]=arima.sim(N,model=list(ar=ar0[1,j],ma=ar0[2,j]) )
	  if(arima_method=="ar1")   for (j in 1:M) X[,j]=arima.sim(N,model=list(ar=ar0[j])  )
	  stat[k,]=preisf(X,method)
	  }
	h=function(x) quantile(x,.95)
	preis=apply(stat,2,h)
	preis_emulation=data.frame(lambda,preis)
	preis_emulation
	}

 #used in preis_emulation
  preisf=function (X,method="mannomatic") {
		N=nrow(X) 
	  	if(method=="semi-mannomatic") {
			m1=apply(X[(N-79):N,],2,mean);sd1=apply(tree[(N-79):N,],2,sd);
	  		X=scale(X,center=m1,scale=sd1)}
 		if(method=="mannomatic") X=mannomatic(X)
		if(method=="covariance") X=scale(X,scale=FALSE)
	  	if (method=="correlation") X=scale(X)
	  	test=svd(X)$d
		preisf=test^2/sum(test^2)
		preisf}	

  #mannomatic standardization
   mannomatic<-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<-d
	mannomatic
	}

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

  #TO EXTEND SERIES BY PERSISTENCE
   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
	}	
	
  #simple plot function showing only mannomatic version
	plotf=function(i,method="mannomatic") {
	 x=preis_emulation(network=target$network[i],step=target$period[i],method) 
	 ylim0=1.2*max(x[1,]);M=nrow(x)
	 plot(1:M,x$preis,type="l",lwd=2,col=4,ylim=c(0,ylim0),las=1,ylab="",xlim=c(0,10.2),xlab="",xaxs="i")
	 points(1:M,x$lambda,col=4)
	 abline(h=2/M,lty=2)
	 title(main=paste(target$network[i],": AD",target$period[i],sep="") )
	 temp=( (x$lambda-x$preis)>0); target$preis[i]=sum(temp)
	 L=target$series[i];points(1:L,x$lambda[1:L],cex=2,font=2,col=2)
	}