################ ##functions shown here3 including some utility smoothing functions ## later on, Jones et al 1998 data is loaded and functions applied ################### #VARIANCE ADJUSTMENT ACCORDING TO OSBORN ET AL 1999 #page 461: variance reduced according to Osborn et al. (1998) Dendrochronologia #see also NAture 1998 # "Each regional mean thus obtained tended to have greater variance during years when few chronologies were available to contribute to the average; #this effect was corrected for by scaling by the square root of the ffective number of independent samples available in each year # n'= n/(1+(n-1)*rbar) #this is a function to implement rbar-adjustment #rbar0 is a vector of length of the total series, calculated in various ways externally bo.adjust<-function(js.mean,rbar0){ count.eff<- count/(1+(count-1)*rbar0);#plot.ts(count.eff) #goes from 2.3 at the beginning to maximum of 5.35 NN<-max(count,na.rm=T) count.eff.max<- NN/(1+(NN-1)*rbar0) ;#count.eff.max# 5.35 var.adj<-sqrt ( count.eff/count.eff.max) ; #equation 7 of Osborn et al. #factor goes from about 0.7 at beginning to 1 bo.adjust<-js.mean*var.adj bo.adjust } ############### ## FUNCTIONS ############## #CALCULATES MEAN CORRELATION FOR A TIME SERIES MATRIX WITH OVERLAP>M TO BE INCLUDED IN CALCULATION #uses functions overlap and template rbar<-function (recon,M=10,method="short") { if(method =="long") recon<-collate(recon,1) N<-ncol(recon) corry<-cor(recon,use="pairwise.complete.obs") temp<-overlap(recon,M) & template (N) rbar<-corry[temp] rbar<-c(rbar) rbar<-rbar[!is.na(rbar)] rbar<-mean(rbar) rbar } #MAKES TEMPLATE IN NXN MATRIX TO ASSIST IN MEAN CORRELATION CALCULATION #this returns a logical matrix with TRUE in left-diagonal to simplify mean correlation calculations template<-function(N) { template<-NULL for (k in 1:N) { template<- c(template, rep(FALSE,k), rep(TRUE,N-k) ) } template<-array(template, dim =c(N,N)) template } #CALCULATES OVERLAPS FOR SPARSE TIME SERIES recon larger than M #this returns a logical matrix of same shape as correlation matrix showing whether the two series have >M overlap overlap<-function (recon,M) { test<-!is.na(recon) temp<-t(test) %*% test overlap<-(temp>M) & !is.na(temp) overlap } ########## ##SMOOTHING UTILTY FUNCTIONS ############## ##gaussian.filter.weights #see http://www.ltrr.arizona.edu/~dmeko/notes_8.pdf gaussian.filter.weights<-function(N,year) { #N number of points; year - is number of years sigma<-year/6 i<-( (-(N-1)/2):((N-1)/2) ) /sigma w<- ((2*pi)^-0.5) * exp(-0.5* i^2) gaussian.filter.weights<-w/sum(w) gaussian.filter.weights } ##truncated.gauss.weights truncated.gauss.weights<-function(year) { a<-gaussian.filter.weights(2*year+1,year) temp<-(a>0.05*max(a)) a<-a[temp] a<-a/sum(a) truncated.gauss.weights<-a truncated.gauss.weights } ##filter with truncated gauss weights with end-period mean padding filter.combine.pad<-function(x,a,M=NA) { temp<-!is.na(x) w<-x[temp] N<-length(w) if(is.na(M)) M<-trunc (length(a)/2)+1 w<-c (rep(mean(w[1:M]),M),w,rep(mean(w[(N-M+1):N]),M)) y<-filter(w,a,method="convolution",sides=2) y<-y[(M+1):(length(y)-M)] z<- rep(NA,length(x)) z[temp]<-y if(class(x)=="ts") z<-ts(z,start=tsp(x)[1]) filter.combine.pad<-cbind(x,z) #filter.combine<-y[,2] filter.combine.pad } ########## ###LOAD COLLATED JONES ET AL 1998 DATA ######### jser<-read.table("http://www.climateaudit.org/data/jser.txt",header=TRUE) jser<-ts(jser[(1000:1991)-999,2:ncol(jser)],start=1000) #end matches archived #NORMALIZE TO 1901-1950 per Jones 1998 Figure 3 #this is described on page 461 col 2 js<-jser m1<-apply(jser[902:951,],2,mean,na.rm=TRUE) sd1<-apply(jser[902:951,],2,sd,na.rm=TRUE) for (k in 1:ncol(js)) { js[,k]<-(jser[,k]-m1[k])/sd1[k] } #COUNT OF SERIES count<-apply(!is.na(js[,1:10]),1,sum); #count of available series ##MEAN OF NH SERIES js.mean<-ts(apply(js[,1:10],1,mean,na.rm=TRUE),start=1000,end=1991) #average of available series #end matches archived N<-length(js.mean) #992 #CONSTANT RBAR0 #this yields quite close matching of variances as shown below rbar0<-rbar(js[,1:10],M=10,method="short");rbar0# 0.0967 #this is the average mean correlation of the 10 series as available js.adj.constant<-ts(bo.adjust(js.mean,rep(rbar0,N)),start=1000) plot.ts(js.adj.constant) #TIME_VARYING RBAR #calculate rbar by individual series for NH #also done for GLB corry<-cor(js[,1:10],use="pairwise.complete.obs") diag(corry)<-NA rbar.jones<-apply(corry,2,mean,na.rm=TRUE);rbar.jones #rbar.jones # Fenno Urals Jasper Svalbard CEng # 0.07537571 0.08835517 0.09195917 0.05443478 0.12338519 # CEur Kameda.melt Jacoby.treeline Briffa.WUSA Crete.O18 # 0.10209179 0.12758733 0.18236210 0.08980311 0.03158346 mean(rbar.jones) # [1] 0.09669378 NH test<-t(array( rep(rbar.jones,N),dim=c(ncol(js),N))) test[is.na(js[,1:ncol(js)])]<-NA rbar.eff<-apply(test,1,mean,na.rm=TRUE) js.adj.timevarying<-ts(bo.adjust(js.mean,rbar.eff),start=1000) #COMPARE TO ARCHIVED #LOAD JONES RECONSTRUCTION url<-"ftp://ftp.ncdc.noaa.gov/pub/data/paleo//contributions_by_author/jones1998/jonesdata.txt" g<-read.table(url,skip=48,fill=TRUE) dimnames(g)[[2]]<-c("year","NH.norm","NH.anom","NH.count","SH.norm","SH.anom","SH.count") archived<-ts(g,start=1000,end=1991) archived<-archived[,2:ncol(archived)] f<-function(x){f<-filter.combine.pad(x,truncated.gauss.weights(51))[,2];f} year<-tsp(archived)[1]:tsp(archived)[2] png("d:/climate/images/proxy/jones1998.png",width=600,height=480) nf<-layout (array( 1:2,dim=c(2,1)),heights=c(1,1.3)) par(mar=c(0,3,1,1)) plot(year,f(archived[,1]),type="l",xlab="",ylab="",axes=FALSE,col="blue")#ylim=c(-1,.5)) axis(side=1,labels=FALSE);axis(side=2,las=1);box(); abline(h=0) lines(1000:1991,f(ts(js.adj.timevarying,start=1000)),lty=2) lines(1000:1991,f(ts(js.adj.constant,start=1000)),lty=3,col="red") par(mar=c(3,3,0,1)) plot(1000:1988,js.adj.timevarying[1:989]-archived[1:989,1],type="l",col="grey60",xlab="",ylab="",axes=FALSE,ylim=c(-1,1)) axis(side=1);axis(side=2,las=1);box(); abline(h=0) #this shows that smoothed version matches pretty well, but that discrepancies in ordinary average dev.off() test<-max(abs( js.adj.timevarying[1:989]-archived[1:989,1])); c( (1000:1988)[ ( js.adj.timevarying[1:989]-archived[1:989,1])==test],test) # 1176.0000000 0.4718136 #this shows that discrepancies up to nearly 0.5 just in emulation #RBAR VALUES ARE TIME-SENSITIVE rbar(js,M=10,method="short");rbar1# 0.05121467 rbar(js[1:901,1:10],M=10,method="short");# 0.0595238 rbar(js[1:801,1:10],M=10,method="short");# 0.07437782 rbar(js[1:601,1:10],M=10,method="short");# 0.01486390 rbar(js[1:901,c(1:4,7:10)],M=10,method="short");# 0.05080707 rbar(js[1:801,c(1:4,7:10)],M=10,method="short");# 0.06646925 rbar(js[1:701,c(1:4,7:10)],M=10,method="short");# 0.04636852 rbar(js[1:601,c(1:4,7:10)],M=10,method="short");# 0.009523411 rbar(js[1:601,1:10],M=10,method="short");# 0.01486390 rbar(js[1:601,],M=10,method="short");# 0.02750064 rbar(js[601:991,],M=10,method="short");# 0.05001516 rbar(js[901:991,c(1:4,7:10)],M=10,method="short");# 0.01855876 rbar(js[601:900,c(1:4,7:10)],M=10,method="short");# 0.05691586 rbar(js[801:991,c(1:4,7:10)],M=10,method="short");# 0.1306393 rbar(js[701:991,c(1:4,7:10)],M=10,method="short");# 0.1016828 rbar(js[751:991,c(1:4,7:10)],M=10,method="short");# 0.1146699