##################### ###THIS IS SCRIPT TO PRODUCE SIMULATIONS, FIGURES AND STATISTICS FOR # 2004GL021750 # Hockey Sticks, Principal Components and Spurious Significance # Stephen McIntyre and Ross McKitrick # ################################################# #PRE REQUISITES: see README # PACKAGE waveslim #PART 1 of this script is functions used in the "narrative" portion of the script #PART 2 of the script is the "narrative" generative all figures and statistics used in the script. ############################################### ############################################### ### ## PART 1 # ##FUNCTIONS REQUIRED FOR NARRATIVE SCRIPT TO PRODUCE SIMULATIONS, FIGURES AND STATISTICS FOR # GRL021750 # Hockey Sticks, Principal components and Spurious Significance # Stephen McIntyre and Ross McKitrick # ################################################ #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 } ##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 localfunction<-function(tree,method2,niter,method.princomp="prcomp",nyear=581) { N<-nrow(tree); #CALCULATE ATTRIBUTES OF NETWORK FOR USE IN RED NOISE if (method2=="arima") {Data1<-apply( tree,2,arima.data)} 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*nyear), dim=c(nyear,niter) ) PC1.princomp<-array (rep (NA,niter*nyear), dim=c(nyear,niter) ) eof1.mannomatic<-array (rep (NA,niter*n), dim=c(n,niter) ) for (i in 1:niter) { #SIMULATE RED/WHITE NOISE #arima version (not used here) if (method2=="arima") {b<-array (rep(NA,nyear*n), dim=c(nyear,n) ) for (k in 1:n) {b[,k]<-arima.sim(list(order = c(1,0,0), ar = Data1[k]), n = nyear)} } #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 #white noise version if(method2=="white.noise") { b<-array( rnorm(nyear*n),dim=c(nyear,n) ) } #NOW DO PRINCIPAL COMPONENT CALCULATIONS: BOTH VERSIONS - WITH AND WITHOUT MBH98 TRANSFORMATION d<-apply(b,2,mannomatic) #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] if(method.princomp=="princomp") {z<-princomp(b); eigen.princomp[,i]<-z$sdev; PC1.princomp[,i]<-z$scores[,1]} else {z<-prcomp(b) eigen.princomp[,i]<-z$sdev; PC1.princomp[,i]<-z$x[,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 } #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 } ###################################################### ### ## PART 2 ### NARRATIVE PORTION TO GENERATE FIGURES AND STATISTICS ### ### ### ######################################################## ##### #REQUIRED DIRECTORY IDENTIFICATIONS ###### url.source<-"ftp://ftp.agu.org/apend/gl/2004GL021750" temp.directory<-"c:/climate/data/ee/sim3" pcs.collation.directory<-file.path(temp.directory,"pcs") idtag<-"2004GL021750-" createsubdirs(temp.directory) createsubdirs(pcs.collation.directory) library(waveslim) library(graphics) ps.options() ######################################## ##MONTE CARLO SIMULATIONS OF HOCKEY STICKS ON TRENDLESS PERSISTENT SERIES ######################## #LOAD NOAMER NETWORK, EXTEND BY PERSISTENCE AND TRUNCATE TO 1980 tree<-read.table(file.path(url.source,"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,]) tree<-tree[,temp] dim(tree) #581 70 #NOW DO SIMULATION method2<-"arfima" NN<-10 NM<-1000 for (nn in 1:NN) { Eigen0<-localfunction(tree,method2,NM) save(Eigen0,file=file.path(temp.directory,paste(method2,"sim",nn,"tab",sep="."))) }#nn-iteration #CALCULATE HOCKEY-STICKNESS ON MANNOMATIC AND ON PRINCOMP stat2<-rep(NA,NN*NM) stat3<-rep(NA,NN*NM) stat4<-rep(NA,NN*NM) index<-1:NM for (nn in 1:NN) { load(file.path(temp.directory,paste(method2,"sim",nn,"tab",sep="."))) stat2[NM*(nn-1)+index]<- apply(Eigen0[[3]],2,hockeystat) stat4[NM*(nn-1)+ index]<- apply(Eigen0[[4]],2,hockeystat) } #nn-iteration stats<-list(stat2,stat4) save(stats,file=file.path(temp.directory,"stats.tab")) ########################## #MAKE FIGURES FROM SIMULATION ########################### #FIGURE 2 load(file.path(temp.directory,"stats.tab")) stat4<-stats[[2]] stat2<-stats[[1]] hist4<-hist(stat4,main="",xlab="",breaks=seq(-2.25,2.25,.25),col="grey50",bg="white",freq=FALSE,ylim=c(0,1),axes=FALSE,font.lab=2) #princomp this has no scouring of center axis(side=1, at = seq(-2,2,1),labels=TRUE,font=2,las=1,cex.axis=.75) fred<-seq(0,1,.25) axis(side=2, at=fred,las=2,cex.axis=.75,font=2) hist(stat2,main="",xlab="",breaks=seq(-2.25,2.25,.25),col="grey50",bg="white",freq=FALSE,axes=FALSE,font.lab=2,ylim=c(0,1.2)) #princomp this has no scouring of center axis(side=1, at = seq(-2,2,1),labels=TRUE,font=2,las=1,cex.axis=.75) fred<-seq(0,1,.25) axis(side=2, at=fred,las=2,cex.axis=.75,font=2) postscript(file.path(url.source,"figure2.eps"), width = 3.3, height = 3.2, horizontal = FALSE, onefile = FALSE, paper = "special", family = "Helvetica",bg="white",pointsize=8) nf <- layout(matrix(c(1,2),2,1,byrow=TRUE)) layout.show(nf) par(mar=c(0,5,5,2)) hist(stat4,main="",xlab="",breaks=seq(-2.25,2.25,.25),col="grey50",bg="white",freq=FALSE,ylim=c(0,1.2),axes=FALSE,font.lab=2) #princomp this has no scouring of center #axis(side=1, at = seq(-2,2,1),labels=TRUE,font=2,las=1,cex.axis=.75) fred<-seq(0,1,.25) axis(side=2, at=fred,las=2,cex.axis=.75,font=2) text(-1.9,1.15,"Centered",font=2) box() par(mar=c(5,5,0,2)) hist(stat2,main="",xlab="'Hockey Stick Index'",breaks=seq(-2.25,2.25,.25),col="grey50",bg="white",freq=FALSE,ylim=c(0,1.2),axes=FALSE,font.lab=2) #princomp this has no scouring of center axis(side=1, at = seq(-2,2,1),labels=TRUE,font=2,las=1,cex.axis=.75) fred<-seq(0,1,.25) axis(side=2, at=fred,las=2,cex.axis=.75,font=2) text(-1.9,1.15,"MBH98",font=2) box() dev.off() #HOCKEY STICK FIGURES CITED IN PAPER #MANNOMATIC #simulations are not hard-wired and results will differ slightly in each run #the values shown here are for a run of 10,000 and results for run of 100 need to be multiplied temp<-(stat2>1)|(stat2< -1); sum(temp) #[1] 9934 temp<-(stat2>1.5)|(stat2< -1.5) ; sum(temp) #[1] 7351 temp<-(stat2>1.75)|(stat2< -1.75); sum(temp) #[1] 2053 temp<-(stat2>2)|(stat2< -2); sum(temp) #[1] 25 #HOCKEY STICK FIGURES CITED IN PAPER #PRCOMP #again these results are not hard-wired temp<-(stat4>1)|(stat4< -1); sum(temp) # #1518 temp<-(stat4>1.5)|(stat4< -1.5); sum(temp) # #120 temp<-(stat4>1.75)|(stat4< -1.75); sum(temp) ##6 temp<-(stat4>2)|(stat4< -2); sum(temp) #[1] 0 ##FIRST EIGENVALUE LOADING localf<-function(x) {localf<-sum(x^2);localf} stat5<-rep(NA,NN*NM) stat6<-rep(NA,NN*NM) stat3A<-rep(NA,NN*NM) stat4A<-rep(NA,NN*NM) index<-1:NM for (nn in 1:NN) { load(file.path(temp.directory,paste(method2,"sim",nn,"tab",sep="."))) stat3A[NM*(nn-1)+ index]<-Eigen0[[1]][1,] stat4A[NM*(nn-1)+ index]<-Eigen0[[2]][1,] stat5[NM*(nn-1)+ index]<-apply(Eigen0[[1]],2,localf) stat6[NM*(nn-1)+ index]<-apply(Eigen0[[2]],2,localf) } #nn iteration stat7<-stat3A^2/stat5 quantile(stat7,c(0.5,0.95,0.99)) # these may differ slightly since simulations not hard-wired # 50% 95% 99% #0.1313634 0.1959358 0.2290063 stat8<-stat4A^2/stat6 quantile(stat8,c(.5,0.95,0.99)) # simulations are not hard-wired #50% 95% 99% #0.04116025 0.04969987 0.05450443 #FIGURE 3 - PC CALCULATIONS #CENTERED CALCULATION z<-princomp(tree) MM.PC1<-ts(-z$scores[,1],start=1400,end=1980) #SIGN-ADJUSTED MM.PC1<-norm(MM.PC1,1902,1980) #MBH98 ARCHIVED PC1 FOR AD1400 STEP url2<-"ftp://holocene.evsc.virginia.edu/pub/MBH98" url2<-"c:/climate/data/new" MBH.PC1<-read.table(file.path(url2,"TREE/ITRDB/NOAMER/BACKTO_1400/pc01.out")) MBH.PC1<- ts( MBH.PC1[,2], start=MBH.PC1[1,1], end=MBH.PC1[nrow(MBH.PC1),1]) MBH.PC1<-norm(MBH.PC1,1902,1980) par(mfrow=c(1,1),mar=c(3,3,2,2)) combine<-ts.union(MBH.PC1,MM.PC1) plot.ts(combine,xlab="",main="") postscript(file.path(url.source,"figure3.eps"), width = 3.3, height = 3.2, horizontal = FALSE, onefile = FALSE, paper = "special", family = "Helvetica",bg="white",pointsize=8) nf <- layout(matrix(c(1,2),2,1,byrow=TRUE)) layout.show(nf) par(mar=c(0,5,5,2)) plot(MBH.PC1,axes=FALSE,type="l",ylab="SD Units",font.lab=2,cex.lab=.75) axis(side=2,at=seq(-6,2,2),las=2,cex.axis=0.75,font=2) box() par(mar=c(5,5,0,2)) plot(MM.PC1,axes=FALSE,type="l",ylab="SD Units",xlab="",font.lab=2,cex.lab=0.75) axis(side=2,at=seq(-3,2,1),las=2,cex.axis=0.75,font=2) axis(side=1, at = seq(1400,2000,100),labels=TRUE,font=2,las=1,cex.axis=.75) box() dev.off() #CALCULATE EXPLAINED VARIANCE MBH.eigen<-read.table("c:/climate/data/new/TREE/ITRDB/NOAMER/BACKTO_1400/eigen.out") MBH.eigen[1,] # 1 0.3818 0.3818 #shows more than 38% z$sdev[1]^2/sum(z$sdev^2) # 0.1911100 #z is carried forward from princomp calculation #this calculation has been benchmarked and replicates the MBH98 stored data ##FIGURE 1 #this graphs a high-RE hockey stick and MBH98 reconstruction #MBH98 RECONSTRUCTION mbhdata<-"ftp://ftp.ngdc.noaa.gov/paleo/paleocean/by_contributor/mann1998/mannnhem.dat" g<-read.table(mbhdata,skip=1) NH.MBH<-ts(g[1:581,2],start=g[1,1],end=g[581,1]) #this plots a number of hockeystick series hockeysticks<-read.table(file.path(url.source,"hockeysticks.txt"),skip=1) #combine<-ts(cbind(hockeysticks[,71],NH.MBH),start=1400,end=1980) postscript(file.path(url.source,"figure1.eps"), width = 3.3, height = 3.2, horizontal = FALSE, onefile = FALSE, paper = "special", family = "Helvetica",bg="white",pointsize=8) nf <- layout(matrix(c(1,2),2,1,byrow=TRUE)) layout.show(nf) par(mar=c(0,5,5,2)) plot(hockeysticks[,71],axes=FALSE,type="l",ylab="",font.lab=2) axis(side=2,at=seq(-.08,.02,.02),las=2,cex.axis=0.75,font=2) box() par(mar=c(5,5,0,2)) plot(NH.MBH,axes=FALSE,type="l",ylab="deg C.",xlab="",font.lab=2,cex.lab=0.75) axis(side=2,at=seq(-.5,.2,.1),las=2,cex.axis=0.75,font=2) axis(side=1, at = seq(1400,2000,100),labels=TRUE,font=2,las=1,cex.axis=.75) box() dev.off() combine<-ts(cbind(hockeysticks[,11],MBH.PC1),start=1400,end=1980) #dimnames(combine)[[2]]<-c("a","b","c","d","e","f","g","MBH") #par(mfrow=c(1,1),mar=c(4,4,2,2)) plot.ts(combine,main="",xlab="",ylab="") par(mfrow=c(2,1),mar=c(4,4,2,2)) plot(1400:1980,hockeysticks[,11],type="l",mar=c(0,4,2,2)) write.table(combine,file="c:/climate/data/ee/examples.txt",sep="\t",quote=FALSE,row.names=FALSE) ################################ ################################# #THE PC1 IN THE MBH98 NORTH AMERICAN NETWORK #IDENTIFY CENSORED SERIES IN MBH98 BACKTO_1400-CENSORED CALCULATION url<-file.path(mbhdata,"TREE/ITRDB/NOAMER") id<-list.files(url) loc<-file.path(url,id[175])#175 is consored list roster<-read.table(loc)#192 loc<-file.path(url,id[176]) #this is main list roster2<-read.table(loc)#192 series censored.list<-roster2[,2][is.na(match(roster2[,2],roster[,2]))] censored.list censored.list<-unlist(strsplit(as.character(censored.list), "\\.")) censored.list<-t(array(censored.list,dim=c(2,length(censored.list)/2)))[,1] #20 series #LIST OF WEIGHTS IN MBH98 PC1 eof01<-read.table(file.path(url,"BACKTO_1400/eof01.out")) name0<-dimnames(tree)[[2]] name0<-unlist(strsplit(name0, "\\.")) name0<-t(array(name0,dim=c(2,length(name0)/2)))[,1] test<-!is.na(match(name0,censored.list)) fred<-cbind(name0[test],eof01[test,2]^2/sum(eof01[,2]^2) ) #20x2 load("c:/climate/data/tree/northamerica.details.tab") #this is collated from ITRDB temp<-!is.na(match(details$id,fred[,1])) #20 selections test2<-details$type temp1<-(test2=="PIAR")|(test2=="PILO")|(test2=="PIBA") #45 selected temp2<-temp&temp1 #15 selected test1<-details$id[temp2] #15 series test3<-!is.na(match(fred[,1],test1)) table1<-fred[test3,]#15 series selected; this can be done manually; I've collected data to selected PIAR, PILO and PIBA automatically test<-!is.na(match(name0,table1[,1])) #15 selected #calculate pct of explained variance from the 15 series F<-as.matrix(eof01[test,2]) m0<-apply(tree,2,mean) A<-scale(tree[,test],center=m0[test],scale=FALSE) %*% F cor(A,MBH.PC1)# 0.9700325 B<-scale(tree,center=m0,scale=FALSE) %*% eof01[,2] z<-data.frame(A,B) fm<-lm(A~B,data=z) summary(fm) #R2 0.9293 #LIST FROM GRAYBILL AND IDSO 1993 #these sites were identified manually graybill.list<-c("co524","co522","co525","co523","ca533","nv510","nv511","nv513","nv512","az510","ca534","ca530","ca528","ca529") graybill.match<-!is.na(match(table1[,1],graybill.list)) table1[,1][graybill.match] #13 picks ########################################## #BENCHMARKING THE REDUCTION OF ERROR STATISTIC FOR THE MBH98 ALGORITHM #use "sparse" NH temperature dataset for verification following MBH98 url<-"ftp://ftp.ngdc.noaa.gov/paleo/paleocean/by_contributor/mann1998/nhem-sparse.dat" use0<-"pairwise.complete.obs" sparse<-read.table(url,skip=1) sparse<-ts(sparse[,2:3],start=sparse[1,1],end=sparse[nrow(sparse),1]) dimnames(sparse)[[2]]<-c("instr","recon") T<-sparse[,"instr"] T<-ts(T[1:(1980-1853)],start=1854,end=1980) calibration<-c(1902,1980) verification<-c(1854,1901) stat<-array (rep (NA,NN*NM*7),dim=c(NN*NM,7)) dimnames(stat)[[2]]<-c( "RE.cal","RE.ver","R2.cal","R2.ver","CE","sign.test","prod.test") for (nn in 1:NN) { load(file.path(temp.directory,paste("arfima.sim",nn,"tab",sep="."))) for (i in 1:NM) { test<-ts(Eigen0[[3]][,i],start=1400,end=1980) z<-ts.union(T,test) N<-tsp(z)[1]-1 z<-data.frame(z) fm<-lm(T~test,data=z[(1902-N):(1980-N),]) predict0<-predict(fm,newdata=z[(1854-N):(1901-N),]) estimator<-ts( c(predict0,fm$fitted.values),start=1854,end=1980) stat[NM*(nn-1)+i,]<- unlist( verification.stats( estimator, T, calibration,verification)) } #i iteration } #nn-iteration quantile (stat[,"RE.ver"],c(.5,.9,.95,.99)) #0.24 0.5102363 0.5440526 0.5910144 quantile(stat[,"R2.ver"],c(.9,.95,.99)) #0.06524158 0.09042600 0.15120589 quantile(stat[,"CE"],c(.9,.95,.99)) #-0.14806575 -0.06879626 0.04128799 quantile(stat[,"sign.test"],c(.9,.95,.99)) # 28 30 32 quantile(stat[,"prod.test"],c(.9,.95,.99)) #1.726455 2.103637 2.807060 #RE.ver 0.46 about 80 percentile ############################################ #SAVE A SELECTION OF HOCKEY STICK SERIES IN ASCII FORMAT order.stat<-order(stat2,decreasing=TRUE)[1:100] order.stat<-sort(order.stat) hockeysticks<-NULL for (nn in 1:NN) { load(file.path(temp.directory,paste("arfima.sim",nn,"tab",sep="."))) index<-order.stat[!is.na(match(order.stat,(1:1000)+(nn-1)*1000))] index<-index-(nn-1)*1000 hockeysticks<-cbind(hockeysticks,Eigen0[[3]][,index]) } #nn-iteration dimnames(hockeysticks)[[2]]<-paste("X",order.stat,sep="") write.table(hockeysticks,file=file.path(url.source,"hockeysticks.txt"),sep="\t",quote=FALSE,row.names=FALSE) ######################################### hockeysticks<-read.table(file.path(url.source,"2004GL021750-hockeysticks.txt"),sep="\t",skip=1) postscript(file.path(url.source,"hockeysticks.eps"), width = 8, height = 11, horizontal = FALSE, onefile = FALSE, paper = "special", family = "Helvetica",bg="white",pointsize=8) nf <- layout(array(1:12,dim=c(6,2)),heights=c(1.1,rep(1,4),1.1)) layout.show(nf) par(mar=c(0,5,1,0)) index<-sample(100,12) plot(hockeysticks[,index[1]],axes=FALSE,type="l",ylab="",font.lab=2,ylim=c(-.1,.03)) axis(side=2,at=seq(-.08,.02,.02),las=2,cex.axis=0.75,font=2) box() for (i in 2:6) { par(mar=c(0,5,0,0)) plot(hockeysticks[,index[i]],axes=FALSE,type="l",ylab="",font.lab=2,ylim=c(-.1,.03)) axis(side=2,at=seq(-.08,.02,.02),las=2,cex.axis=0.75,font=2) box() } par(mar=c(0,0,1,5)) plot(hockeysticks[,index[7]],axes=FALSE,type="l",ylab="",font.lab=2,ylim=c(-.1,.03)) #axis(side=2,at=seq(-.08,.02,.02),las=2,cex.axis=0.75,font=2) box() fm<-lm(MBH.PC1~NH.MBH) hockeysticks[,index[8]]<-predict(fm,newdata=NH.MBH) par(mar=c(0,0,0,5)) plot(hockeysticks[,index[8]],axes=FALSE,type="l",ylab="",font.lab=2,ylim=c(-.1,.03)) box() for (i in 9:12) { par(mar=c(0,0,0,5)) plot(hockeysticks[,index[i]],axes=FALSE,type="l",ylab="",font.lab=2,ylim=c(-.1,.03)) box() } dev.off() index#[1] 98 23 65 97 20 85 46 29 68 89 76 86