############## #this is version with padding #this needs to be tidied up and I'll try to do so sometime. The end tails off a bit in clarity #MOBERG ET AL COLLATION #WAVELET ANALYSIS OF LOWF library(waveslim) source("c:/climate/scripts/treeconfig.functions.txt") #COLLATE MOBERG SUMMARY #SI at http://www.nature.com/nature/journal/v433/n7026/suppinfo/nature03265.html url<-"http://www.nature.com/nature/journal/v433/n7026/extref/nature03265-s6.doc" ##fred<-readLines(url) recon<-read.table(url,skip=19,nrow=1998-19) #1979 9 #from 1 to 1979 name<-c("year","recon","lowf","lowf.lb","lowf.ub","ci.lb","ci.ub","ci2.lb","ci2.ub") dimnames(recon)[[2]]<-name recon<-ts(recon[,2:ncol(recon)],start=1,end=1979) col0<-c("grey80","green","blue","yellow","red","black","orange","magenta","salmon") temp<-(recon== -9.99) recon[temp]<-NA temp<-(recon< -9) recon[temp]<-NA ts.plot(recon[,1:5],col=col0) #this replicates part of Supplementary Figure 3. #presumably the main figure tacks on the CRU instrumental series #rm(recon) #LOAD MOBERG DATA #previously collated #load("c:/climate/data/moberg/mob.tab") #id<-dimnames(mob)[[2]] #mob<-cbind(1:1996,mob) #dimnames(mob)[[2]]<-c("year",id) #write.table(mob,file="c:/climate/data/moberg/moberg.proxy.txt",sep="\t") mob<-read.table("http://www.climateaudit.org/data/moberg.proxy.txt",fill=TRUE,sep="\t",header=TRUE) dim(mob) #1996 18 #CENTER AND SCALE DATA #Moberg: From the wavelet transform, #which was calculated using standardized (that is, unitless) proxy #series, a dimensionless reconstruction was obtained by calculating #the inverse wavelet transform (11) Mallat. #page 624 right column mob<-scale(mob) #center and scale the data apply(mob,2,mean,na.rm=T) #all zero apply(mob,2,sd,na.rm=T) #all 1 #PAD SERIES #this function pads the end of a series with mean over the last 200 f<-function(x,K=200) { N<-length(x) m1<-mean(x,na.rm=TRUE) y<-x[(N-K):N]-m1 temp<-!is.na(y) z<-y[temp] M<-length(z) m0<-mean(z[(M-49):M]) y[!temp]<-m0 x[(N-K):N]<-y+m1 f<-x f} #this function pads the series with mean over the first 50 f1<-function(x,K=200) { N<-length(x) m1<-mean(x,na.rm=TRUE) y<-x[1:K]-m1 temp<-!is.na(y) z<-y[temp] M<-length(z) m0<-mean(z[1:50]) y[!temp]<-m0 x[1:K]<-y+m1 f1<-x f1} apply(!is.na(mob),2,sum) #check number of readings # Agassiz Melt % GRIP Borehole Conroy Lake Chesapeake Mg/Ca # 1966 1977 1934 1995 # Sargasso Sea dO18 Caribbean dO18 Tsuolbmajavri Sediments Norway Stalagmite # 1970 1818 1985 1862 # Beijing Stalagmite China (Yang) Arabian Sea Upwelling negative Methesulah Walk # 1985 1990 1847 1980 #negative Indian Garden negative White Mt Tornetrask (Grudd RW) Yamal # 1980 1962 1993 1996 # Taimyr Yakutia # 1996 595 #pad the two series for (i in 1:17) { mob[,i]<-f(mob[,i]) mob[,i]<-f1(mob[,i]) } apply(!is.na(mob),2,sum)#all 1996 except for Yakutia which is 595 round(apply(mob,2,mean,na.rm=T),2) ##this is non-zero now round(apply(mob,2,sd,na.rm=T),2) ##this is non-zero now #RECENTER AFTER PADDING WITHOUT ADJUSTING SCALE mob<-scale(mob,scale=FALSE) round(apply(mob,2,mean,na.rm=T),2) ##this 0 now #EXCLUDE YAKUTIA PENDING DATA RECEIPT mob0<-mob mob<-mob[,1:17] ################ #WAVELET DECOMPOSITIONS ############################3 #this is a plot of high-freq and low-freq decompositions by proxies posted up at climateaudit #Nine low-resolution proxy series (Table 1) are available as calibrated local/regional #temperature indicators (numbers 2–10). Five are annual mean (nos 2, 5, 6, 8, 10), one is #a spring (no. 4) and three are summer (nos 3, 7, 9) temperature indicators. #Multicentennial variability (highlighting .340-yr scales) in these series was extracted #by first applying the WT to each series and then the inverse WTonly for wavelet voices 14–22.: mob.lowf<-array( rep(NA,18*1996),dim=c(1996,18)) mob.highf<-array( rep(NA,18*1996),dim=c(1996,18)) for (k in c(1:17)) { temp<-!is.na(mob[,k]) year<-(1:(1996))[temp] #133 1925 x<-mob[temp,k] x<-x-mean(x) N<-length(x) J<-trunc(log(N,2)) model.la8 <- mra(x, "la8", J, "modwt") #1024 is max fred<-NULL for (i in 1:(J+1)) { fred<-cbind(fred,model.la8[[i]])} test1<-apply(fred[,1:5],1,sum) test2<-apply(fred[,6:ncol(fred)],1,sum) ylim0<-range(c(test1,test2)) #this is a plot of high-freq and low-freq decomposition by proxy nf<-layout(array(1:3,dim=c(3,1)),heights=c(1.2,1,1.4)) par(mar=c(0,4,2,2)) plot(year,x,type="l", axes=FALSE, ylab="", main=paste(k,". ",dimnames(mob)[[2]][k],sep=""),xlim=c(0,2000)) axis(side=2) axis(side=1,labels=FALSE) box() par(mar=c(0,4,0,2)) plot(year,test1 , type="l",axes=FALSE, ylab="<=32 Years",xlim=c(0,2000),font=2,ylim=ylim0) axis(side=2) axis(side=1,labels=FALSE) box() par(mar=c(4,4,0,2)) plot(year,test2 ,type="l",axes=FALSE,ylab=">=64 Years",font=2,xlim=c(0,2000),xlab="",ylim=ylim0) abline(h=0) axis(side=2) axis(side=1) box() mob.lowf[,k]<-test2 mob.highf[,k]<-test1 } save(mob.lowf,file="c:/climate/data/moberg/mob.lowf.tab") save(mob.highf,file="c:/climate/data/moberg/mob.highf.tab") #STATISTICS #check hockeystick-ness #most prononuced hockey sticks are Agassiz melt %, Arabian Sea upwelling, Yamal Chinayang f<-function(x) { f<-hockeystat(x,N=(1996-79));f} round(apply(mob,2,f),2) # Agassiz Melt % GRIP Borehole Conroy Lake Chesapeake Mg/Ca # 2.85 -0.48 -1.88 0.38 #Sargasso Sea dO18 Caribbean dO18 Tsuolbmajavri Sediments Norway Stalagmite # -0.01 0.56 -0.17 -0.45 #Beijing Stalagmite China (Yang) Arabian Sea Upwelling negative Methesulah Walk # 1.16 1.92 2.91 -0.38 #negative Indian Garden negative White Mt Tornetrask (Grudd RW) Yamal # 0.34 0.05 0.98 1.25 # Taimyr # 0.61 ################ #NOW DO WAVELET DECOMPOSITION BY VOICES #in this case this is done through discrete wavelet decomposition rather than CWT #Moberg description #The multi-proxy reconstruction in Fig. 2b was constructed as follows. First, the data set of #19 proxy series in Fig. 1 was divided into two half-hemispheric subsets (western and #eastern half). After wavelet transformation of each individual proxy series into 22 voices as #described above (here using standardized series obtained by subtracting the mean and #dividing by the standard deviation, after linear interpolation to annual resolution), a #scale-by-scale averaging of wavelet voices was undertaken. For each wavelet scale from 1 to #9 (Fourier timescales,80 yr), the voices from the individual tree-ring series were averaged #and, similarly, for each wavelet scale from 10 to 22 (Fourier timescales .80 yr) the #corresponding voices from the individual low-resolution proxy series were averaged. #Merging of the tree-ring-data average voices 1–9 with the low-resolution-data average #voices 10–22 gave one complete WTwith voices 1–22 for each half-hemispheric subset. #The average of these twoWTs then gave the whole-hemispheric WTshown in Fig. 2e. The #rationale for dividing into half-hemispheric subsets is similar to that used in gridding; #namely, to ensure that regions with many data series are weighted the same as regions (of #the same size) with few data series. However, for our data set, little difference would be #obtained if all data were placed in one whole-hemispheric set. A dimensionless Northern #Hemisphere temperature reconstruction was obtained by calculating the inverse #transform. Finally, this time series was calibrated by adjusting its variance and mean value #to be the same as in the instrumental data19 in the overlapping period AD 1856–1979 #(Fig. 2b). Separate estimates of uncertainties due to variance among individual lowresolution #proxy data and in the determination of the variance scaling factor and constant #adjustment term are explained in Supplementary Information. ################# ##CALCULATION OF VOICE AVERAGES #his makes 3-D composition by scale proxy and year #this is repetitive of earlier calculation but saves results at each individual scale ratther than low-freq and high-freq aggregates mob.voice<-array( rep(NA,18*1996,11),dim=c(1996,11,18)) for (k in c(1:17)) { temp<-!is.na(mob[,k]) year<-(1:1996)[temp] #133 1925 x<-mob[temp,k] x<-x-mean(x) N<-length(x) J<-trunc(log(N,2)) x<-scale(x) model.la8 <- mra(x, "la8", J, "modwt") #1024 is max fred<-NULL for (i in 1:(J+1)) { fred<-cbind(fred,model.la8[[i]])} #this is the same plot as before par(mfcol=c(3,1), pty="m", mar=c(5-2,4,4-2,2)) plot(year,x,type="l", axes=FALSE, ylab="", main=paste(k,". ",dimnames(mob)[[2]][k],sep=""),xlim=c(0,2000)) axis(side=2) box() plot(year, apply(fred[,1:5],1,sum), type="l",axes=FALSE, ylab="Scale<=32 Years",xlim=c(0,2000)) axis(side=2) box() plot(year, apply(fred[,6:ncol(fred)],1,sum),type="l",axes=FALSE,ylab="Scale >=64 Years",xlim=c(0,2000)) axis(side=2) axis(side=1) box() for (i in 1:11) {mob.voice[1:1996,i,k]<-fred[1:1996,i]} #collects results into data matrix } #NOW DO VOICE AVERAGES mob.voice.average<-array(rep(NA,11*1996),dim=c(1996,11)) for (i in 1:5) {mob.voice.average[,i]<-apply(mob.voice[,i,12:17],1,mean,na.rm=TRUE)} #high frequency are series 12-17 for(i in 6:11) {mob.voice.average[,i]<-apply(mob.voice[,i,1:11],1,mean,na.rm=TRUE)} #low frequency are series 1-11 scale0<-2^(6:11) mob.recon.lowf<-apply(mob.voice.average[,6:11],1,sum,na.rm=TRUE) #take sum of low frequency mob.recon.highf<-apply(mob.voice.average[,1:5],1,sum,nar.rm=TRUE) #take sum of high-frequency mob.final<-mob.recon.lowf+mob.recon.highf #calculate total #PLOT RESULTS ylim0<-range(mob.voice.average[,6:11]) year<-1:1996 nf<-layout(array(1:7,dim=c(7,1)),heights=c(1.4,1,1,1,1,1,1.3)) par(mar=c(0,4,2,2)) x<-mob.final plot(year,x,type="l", axes=FALSE, ylab="low-f",xlim=c(0,2000)) axis(side=2) axis(side=1,labels=FALSE) box() for (j in 6:10) { par(mar=c(0,4,0,2)) x<-mob.voice.average[,j] plot(year,x, type="l",axes=FALSE, ylab=paste(scale0[j-5],"yrs"),xlim=c(0,2000),ylim=ylim0) axis(side=2) axis(side=1,labels=FALSE) box() } par(mar=c(4,4,0,2)) x<-mob.voice.average[,11] plot(year,x, type="l",axes=FALSE, ylab=paste(scale0[6],"yrs"),xlim=c(0,2000),ylim=ylim0) axis(side=2) axis(side=1) box() #NEXT FIGURE nf<-layout(array(1:3,dim=c(3,1)),heights=c(1.2,1,1.3)) ylim0<-range(c(mob.recon.lowf,mob.recon.highf)) par(mar=c(0,4,2,2)) x<-mob.final plot(year,x,type="l", axes=FALSE, ylab="low-f",xlim=c(0,2000)) axis(side=2) axis(side=1,labels=FALSE) box() par(mar=c(0,4,0,2)) x<-mob.recon.lowf plot(year,x, type="l",axes=FALSE, ylab=paste("low-freq"),xlim=c(0,2000),ylim=ylim0) axis(side=2) axis(side=1,labels=FALSE) box() par(mar=c(4,4,0,2)) x<-mob.recon.highf plot(year,x, type="l",axes=FALSE, ylab=paste("high-freq"),xlim=c(0,2000),ylim=ylim0) axis(side=2) axis(side=1) box() ############## ##CALIBRATION TO NH #Moberg: #To calibrate the reconstruction, its mean value and variance #were adjusted to agree with the instrumental record of Northern #Hemisphere annual mean temperatures19 in the overlapping period #AD 1856–1979 (Fig. 2b). ##LOAD CRU url<-"ftp://ftp.cru.uea.ac.uk/data/tavenh2v.dat" #this is different than tavenh2.dat h<-scan(url) N<-length(h) start0<-h[1] h<-array(h,dim=c(27,N/27)) h<-t(h) CRU<-ts(h[,14],start=start0) #1856 - 2005 tsp(CRU) ##[1] 1856 2005 1 ##DO DECOMPOSITION OF CRU INTO HIGH-FREQ and LOW-FREQ temp<-!is.na(CRU[1:(1996-1855)]) year<-1856:1996# x<-CRU[temp] mean.CRU<-mean(CRU[temp]) x<-x-mean(x) # -0.09636 N<-length(x) J<-trunc(log(N,2)) model.la8 <- mra(x, "la8", J, "modwt") #1024 is max fred<-NULL for (i in 1:(J+1)) { fred<-cbind(fred,model.la8[[i]])} nf<-layout(array(1:3,dim=c(3,1)),heights=c(1.2,1,1.3)) par(mar=c(0,4,2,2)) plot(year,x,type="l", axes=FALSE, ylab="", main=paste("CRU"),xlim=c(1850,2000)) axis(side=2) box() par(mar=c(0,4,0,2)) plot(year, apply(fred[,1:5],1,sum), type="l",axes=FALSE, ylab="Scale<=32 Years",xlim=c(1850,2000)) axis(side=2) box() par(mar=c(4,4,0,2)) plot(year, apply(fred[,6:ncol(fred)],1,sum),type="l",axes=FALSE,ylab="Scale >=64 Years",xlim=c(1850,2000)) axis(side=2) axis(side=1) box() sd.CRU.highf<-sd(apply(fred[,1:5],1,sum,na.rm=T));sd.CRU.highf # 0.1900008 sd.CRU.lowf<-sd(apply(fred[,6:ncol(fred)],1,sum,na.rm=T));sd.CRU.lowf # 0.1763978 mean.CRU.highf<-mean(apply(fred[,1:5],1,sum,na.rm=T));mean.CRU.highf # 0 mean.CRU.lowf<-mean(apply(fred[,6:ncol(fred)],1,sum,na.rm=T));mean.CRU.lowf # 0 sd.CRU<-sd(x) # 0.2671070 sd.recon.lowf<-sd(mob.recon.lowf);sd.recon.lowf # 0.3748422 sd.recon.highf<-sd(mob.recon.highf);sd.recon.highf # 0.4646824 #now rescale moberg anomaly to CRU data mob.lowf.anomaly<-(mob.recon.lowf-mean(mob.recon.lowf))* sd.CRU.lowf/sd.recon.lowf+mean.CRU mob.highf.anomaly<-(mob.recon.highf-mean(mob.recon.highf))* sd.CRU.highf/sd.recon.highf mob.final.anomaly<-mob.lowf.anomaly+mob.highf.anomaly+mean.CRU #COMPARE EMULATION TO REPORTED c(mean(mob.final.anomaly),mean(recon[1:1979,"recon"]))# -0.0963600 -0.3537393 mob.final.anomaly<-mob.final.anomaly+mean(recon[1:1979,"recon"])-mean(mob.final.anomaly) test<-min(mob.final.anomaly[1:1979]-recon[,"recon"]);test# -0.93 (1:1979)[(mob.final.anomaly[1:1979]-recon[,"recon"])==test]# 1489 test<-max(mob.final.anomaly[1:1979]-recon[,"recon"]);test# 1.15 (1:1979)[(mob.final.anomaly[1:1979]-recon[,"recon"])==test]# 1580 test<-((mob.final.anomaly[1:1979]-recon[,"recon"]) >1) ; (1:1979)[test]# 1580 test<- mob.final.anomaly[1:1979]-recon[,"recon"] quantile(test,c(.025,.975)) #-0.5535304 0.6309730 #SHOW EMULATION VS REPORTED year<-1:1996 nf<-layout(array(1:3,dim=c(3,1)),heights=c(1.2,1,1.3)) par(mar=c(0,4,2,2)) plot(year,mob.final.anomaly,type="l", axes=FALSE, ylab="Emulation", xlim=c(0,2000)) axis(side=2);axis(side=1,labels=FALSE); box() par(mar=c(0,4,0,2)) plot(1:1979, recon[,"recon"],type="l",axes=FALSE,ylab="Moberg",xlim=c(0,2000)) axis(side=2);axis(side=1,labels=FALSE); box() par(mar=c(4,4,0,2)) plot(1:1979, mob.final.anomaly[1:1979]-recon[1:1979,"recon"],type="h",axes=FALSE,xlab="",ylab="Residual",xlim=c(0,2000)) axis(side=2); axis(side=1) box() #COMPARE TO CRU nf<-layout(array(1:2,dim=c(2,1)),heights=c(1.1,1.3)) par(mar=c(0,4,2,2)) plot(year,mob.final.anomaly,type="l", axes=FALSE, ylab="",xlim=c(0,2000)) axis(side=2);axis(side=1,labels=FALSE); box() par(mar=c(4,4,0,2)) plot(1:1979,recon[,"recon"],type="l", axes=FALSE, ylab="",xlim=c(0,2000),col="red") axis(side=2);axis(side=1); box() ###PLOT RESIDUALS nf<-layout(array(1:2,dim=c(2,1)),heights=c(1.1,1.3)) par(mar=c(0,3,1,1)) plot(1856:1979,recon[1856:1979,"recon"],type="l",xlab="",ylab="",axes=FALSE) lines(1856:1979,CRU[1:(1979-1855)],col="red") legend(1850,0.28,fill=c("black","red"),legend=c("Moberg","CRU")) axis(side=2);axis(side=1,labels=FALSE); box() par(mar=c(3,3,0,1)) plot(1856:1979,resid.moberg,type="h",xlab="",ylab="",axes=FALSE) axis(side=2);axis(side=1); box() abline(h=0.59*.23*.5,col="red") abline(h=-0.59*.23*.5,col="red") #STATISTICS ON MODEL cor(recon[1856:1979,"recon"],CRU[1:(1979-1855)]) # 0.2976473 resid.moberg<-recon[1856:1979,"recon"]-CRU[1:(1979-1855)] fm<-lm(resid.moberg~1) summary(fm) #Coefficients: # Estimate Std. Error t value Pr(>|t|) #(Intercept) 0.0001919 0.0185207 0.01 0.992 #Residual standard error: 0.2062 on 123 degrees of freedom dwtest(fm) #DW = 0.6045, p-value = 2.381e-15 arima(fm$residuals,order=c(1,0,0)) # 3Coefficients: # ar1 intercept # 0.6938 -0.0049 #s.e. 0.0639 0.0425 #sigma^2 estimated as 0.02172: log likelihood = 61.16, aic = -116.32 arima(fm$residuals,order=c(1,0,1)) # ar1 ma1 intercept # 0.7702 -0.1483 -0.0063 #s.e. 0.0779 0.1183 0.0476 #sigma^2 estimated as 0.02143: log likelihood = 61.98, aic = -115.96