##SCRIPT TO CALCULATE AR1 COEFFICIENTS FOR MBH98 NOAMER NETWORK

   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
	}	


   ##LOAD TREE RING NETWORK
	url<-"ftp://ftp.agu.org/apend/gl/2004GL021750"
	tree<-read.table(file.path(url,"2004GL021750-NOAMER.1400.txt"),header=TRUE,sep="\t") #collated AD1400 network
	tree<-ts(tree[1:581,2:ncol(tree)],start=1400)
	tree<-extend.persist(tree)  
	temp<-!is.na(tree[1,])
	tree<-tree[,temp]
	dim(tree) #581 70
	  
   ##AR1
   	stat<-array(NA,dim=c(70,2))
	for (i in 1:70) { stat[i,]<-arima(tree[,i],order=c(1,0,0))$coef}
	row.names(stat)<-dimnames(tree)[[2]]
	stat[order(stat[,1],decreasing=TRUE),] #bristlecones are all high AR1
	range(stat[,1]) #0.02679753 0.79111163
   	hist(stat[,1],breaks=seq(0,1,.1),xlab="AR1", main=c("NOAMER AR1 Coefficients","ARMA(1,0) Model"),col="grey80")
 	mean(stat[,1])#0.4145068
	sum(stat[,1]<.15) #3

  ##ARMA (1,1)	
   	stat<-array(NA,dim=c(70,3))
	for (i in 1:70) {stat[i,]<-arima(tree[,i],order=c(1,0,1))$coef}
	row.names(stat)<-dimnames(tree)[[2]]
	stat[order(stat[,1],decreasing=TRUE),] #bristlecones are all high AR1
	range(stat[,1]) # 0.3301175 0.9866811
   	hist(stat[,1],breaks=seq(0,1,.1),xlab="AR1", main=c("NOAMER AR1 Coefficients","ARMA(1,1) Model"),col="grey80")
	mean(stat[,1])# 0.7850041
	sum(stat[,1]<.15) #0