###MAKE INTO A FUNCTION #this is a utility function to calculate the items in Santer Table 1 for a given time series #this is for monthly data f=function(x,id,end0=2050,start0=1979,M=19, trend.model=trend.model.T2LT, sd.model=sd.model.T2LT ) { f=rep(NA,7) temp= (time(x)< (end0+1))&(time(x)>=start0) x=ts(x[temp],start=c(start0,1),freq=12) year=c(time(x))/10; N=length(x) f[3]=sd(x) fm= lm (x~year) trend.obs=f[1]= fm$coef[2] r=f[4]= arima (fm$residuals,order=c(1,0,0))$coef[1]; # f[5]= N * (1-r)/(1+r) ; se.obs=f[2] = sqrt((N-2)/(f[5]-2))* summary(fm)$coef[2,"Std. Error"] f[7]= (trend.model-trend.obs)/ (sd.model/sqrt(M-1)) f[6]= (trend.model-trend.obs)/ sqrt( (sd.model/sqrt(M-1))^2 + se.obs^2) names(f)=c("trend","se","sd","r1","neff","d1star","d1") info=rbind( round(f,3),santer[paste(id),]) row.names(info)=c("emulation","santer") C= (sd.model/sqrt(M-1))^2 + se.obs^2;C # 0.0195 D=(sd.model/sqrt(M-1))^2 /(M-1) +se.obs^2/(santer[id,"neff"]-2) ;D # 0.00155 f=list(x=x,fm=fm,info=info,adj_df=C/D) f } ###LOAD DATA #1. Santer 2008 Table 1 loc="http://www.climateaudit.org/data/models/santer_2008_table1.dat" santer=read.table(loc,skip=1) names(santer)=c("item","layer","trend","se","sd","r1","neff") row.names(santer)=paste(santer[,1],santer[,2],sep="_") santer=santer[,3:ncol(santer)] #insert Table # info santer$d1=NA;santer$d1star=NA santer["UAH_T2LT",c("d1","d1star")]= c(1.11,7.16) #from Table III for UAH santer["RSS_T2LT",c("d1","d1star")]= c(0.37,2.25) #from Table III for UAH santer["UAH_T2",c("d1","d1star")]= c(1.19,6.78) #from Table III for UAH santer["RSS_T2",c("d1","d1star")]= c(0.44,2.48) #from Table III for UAH #TABLE III CALCULATIONS trend.model=trend.model.T2LT=santer["Multi-model_mean_T2LT","trend"]; # 0.215 # c(trend.model.T2LT , mean(douglass$trend_LT) ) # 0.215 0.208 sd.model=sd.model.T2LT=santer["Inter-model_S.D._T2LT",1]; #c(sd.model.T2LT,sd(douglass$trend_LT) ) #0.0920 0.0913 trend.obs=santer["UAH_T2LT","trend"];trend.obs # 0.06 se.obs=santer["UAH_T2LT","se"];se.obs #0.138 M=19 #from Santer