##IPCC TABLE 3.2 TREND CALCULATIONS library(nlme);library(lmtest); ##1. READ CRUTEM, HADSST, HADCRU #function to download annual versions f=function(url){ h<-scan(url) start0<-h[1]; N<-length(h) h<-array(h,dim=c(27,N/27));h<-t(h) f<-ts(h[,14],start=start0,end=2006) #1850 - 2002 f } hadcru=f("http://www.cru.uea.ac.uk/cru/data/temperature/hadcrut3nh.txt") #1850 2007 hadsst=f("http://www.cru.uea.ac.uk/cru/data/temperature/hadsst2nh.txt") crutem=f("http://www.cru.uea.ac.uk/cru/data/temperature/crutem3nh.txt") #reference period 1961-1990 #from 1851:2006 year=1850:2006; N=length(year) Data=data.frame(rep(year,3),c(rep("crutem",N),rep("hadcru",N),rep("hadsst",N) ),c(crutem,hadcru,hadsst) ) names(Data)=c("year","dataset","tt") Data=groupedData(tt~year|dataset,data=Data) ##2. CALCULATE HADSST MODELS #reported 0.038 ± 0.011 deg C per decade #R2=0.54 R1=0.85 DW=2.20 #simple linear model temp=(Data$dataset=="hadsst")&(Data$year<=2005) fm=lm(tt~year,data=Data[temp,]); coef(fm) # (Intercept) year #-6.689188285 0.003400966 summary(fm)$r.squared #[1] 0.4520157 #this does not match reported R2 #best fit with reported slope of 0.038 deg /decade fm2=lm(I(tt-(0.038/10) *year)~1,data=Data[temp,]) #coef divided by 10 to match decade scale coef(fm2) #(Intercept) # -7.458327 test=fm2$residuals+ (0.038/10) *(1850:2005) #lme model with "ML" fm1=lme(tt~year,data=Data[temp,],method="ML"); coef(fm1) # (Intercept) year #hadsst -6.818371 0.003468897 anova(fm1,fm) # Model df AIC BIC logLik Test L.Ratio p-value #fm1 1 6 -99.69398 -81.35651 55.84699 #fm 2 3 -105.69398 -96.52524 55.84699 1 vs 2 2.357855e-08 1 #lme model with REML fm2=update(fm1,method="REML");coef(fm2) # (Intercept) year #hadsst -6.818371 0.003468897 # Model df AIC BIC logLik Test L.Ratio p-value #fm2 1 6 -78.54389 -60.28334 45.27194 #fm 2 3 -84.54389 -75.41361 45.27194 1 vs 2 2.842171e-14 1 ##AR1 correlation structure fm3=update(fm2,corr=corCAR1(form=~year|dataset) );coef(fm3) # (Intercept) year #hadsst -7.152831 0.003646039 ################ summary(fm3) ################## Linear mixed-effects model fit by REML Data: Data[temp, ] AIC BIC logLik -212.4479 -191.1440 113.2240 Random effects: Formula: ~year | dataset Structure: General positive-definite StdDev Corr (Intercept) 3.809902e-02 (Intr) year 1.975549e-05 0 Residual 1.790178e-01 Correlation Structure: Continuous AR(1) Formula: ~year | dataset Parameter estimate(s): Phi 0.7808196 Fixed effects: tt ~ year Value Std.Error DF t-value p-value (Intercept) -7.152831 1.621908 155 -4.410135 0 year 0.003646 0.000841 155 4.335441 0 Correlation: (Intr) year -0.999 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -2.25698742 -0.70827227 -0.05559096 0.52817752 2.53366188 Number of Observations: 157 Number of Groups: 1 #################################### fm3A=update(fm3,correlation = corAR1());coef(fm3A) #(Intercept) year #hadsst -7.152831 0.003646039 fm4=update(fm3,tt~I(year-1928) );coef(fm4) # (Intercept) I(year - 1928) #hadsst -0.1232686 0.003646039 #################### summary(fm4) ###################### Linear mixed-effects model fit by REML Data: Data[temp, ] AIC BIC logLik -212.4479 -191.1440 113.2240 Random effects: Formula: ~I(year - 1928) | dataset Structure: General positive-definite StdDev Corr (Intercept) 0.0380991295 (Intr) I(year - 1928) 0.0008406495 0 Residual 0.1790177714 Correlation Structure: Continuous AR(1) Formula: ~year | dataset Parameter estimate(s): Phi 0.7808196 Fixed effects: tt ~ I(year - 1928) Value Std.Error DF t-value p-value (Intercept) -0.12326856 0.05511832 155 -2.236435 0.0268 I(year - 1928) 0.00364604 0.00118893 155 3.066654 0.0026 Correlation: (Intr) I(year - 1928) 0 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -2.25698742 -0.70827227 -0.05559096 0.52817752 2.53366187 Number of Observations: 157 Number of Groups: 1 ########################################## ################################################ ##DURBIN-WATSON TEST ON TREND RESIDUALS ################################################## dwtest(fm) #with fresh calculation #DW = 0.4852, p-value < 2.2e-16 dwtest(fm2$residuals~1) #DW = 0.4852, p-value < 2.2e-16 dwtest(fm3$residuals[,1]~1) #DW = 0.4841, p-value < 2.2e-16 ########################### ##ARIMA FIT TO RESIDUALS ############################# arima.fit=arima(fm$residuals,order=c(1,0,0));arima.fit #with fresh fitted residuals # ar1 intercept # 0.7571 0.0062 #s.e. 0.0521 0.0358 #sigma^2 estimated as 0.01234: log likelihood = 121.78, aic = -237.56 arima.fit2=arima(fm2$residuals[,1],order=c(1,0,0));arima.fit2 #with fresh fitted residuals #Coefficients: # 0.7571 0.0062 #s.e. 0.0521 0.0358 #sigma^2 estimated as 0.01234: log likelihood = 121.78, aic = -237.56 arima.fit3=arima(fm3$residuals[,1],order=c(1,0,0));arima.fit3 #with fresh fitted residuals # 0.7573 -0.0008 #s.e. 0.0520 0.0359 #sigma^2 estimated as 0.01234: log likelihood = 121.8, aic = -237.6 ###################### ##DURBIN-WATSON TEST ON ARIMA.FIT RESIDUALS ######################### dwtest(arima.fit$residuals~1) # # DW = 1.8888, p-value = 0.2415 dwtest(arima.fit2$residuals~1) #assuming reported trend # DW = 1.8888, p-value = 0.2415 dwtest(arima.fit3$residuals~1) #assuming reported trend # DW = 1.8892, p-value = 0.2423 ########################### ########################## #AN EXPERIMENT ON RANDOM WALK: START FRESH ############################# x=cumsum(rnorm(1055)) y=x[156:1055]-x[1:900] index=(1:900)[y>30] j=sample(index,1) w=x[(j-155):j];plot.ts(w) ##linear model on trend portion of random walk year=1850:2005 Data.rw=data.frame(year,rep("id1",156),w) names(Data.rw)=c("year","dataset","tt") Data.rw=groupedData(tt~year|dataset,data=Data.rw) fm=lm(w~year,data=Data.rw); ############# summary(fm) ############## #Coefficients: # Estimate Std. Error t value Pr(>|t|) #(Intercept) -1.394e+02 9.650e+00 -14.44 <2e-16 *** #year 9.058e-02 5.005e-03 18.10 <2e-16 *** #Residual standard error: 2.815 on 154 degrees of freedom #Multiple R-Squared: 0.6802, Adjusted R-squared: 0.6781 #F-statistic: 327.5 on 1 and 154 DF, p-value: < 2.2e-16 layout(1) plot(year,w,type="l",xlab="",ylab="") lines(year,fm$fitted.values) #lme model with "ML" fm1=lme(tt~year,data=Data.rw,method="REML"); #################### summary(fm1) ############################## Linear mixed-effects model fit by REML Data: Data.rw AIC BIC logLik 785.5467 803.7684 -386.7733 Random effects: Formula: ~year | dataset Structure: General positive-definite StdDev Corr (Intercept) 0.6010743440 (Intr) year 0.0003117564 0 Residual 2.8152810565 Fixed effects: tt ~ year Value Std.Error DF t-value p-value (Intercept) -139.38959 9.669145 154 -14.41592 0 year 0.09058 0.005015 154 18.06206 0 Correlation: (Intr) year -0.996 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -2.3742503 -0.5752575 0.1789003 0.7394398 1.8227061 Number of Observations: 156 Number of Groups: 1 #################################################### ##################### fm3=update(fm1,tt~I(year-1928),correlation = corAR1()); summary(fm3) ##################### Linear mixed-effects model fit by REML Data: Data.rw AIC BIC logLik 426.29 447.5487 -206.145 Random effects: Formula: ~I(year - 1928) | dataset Structure: General positive-definite StdDev Corr (Intercept) 0.88254828 (Intr) I(year - 1928) 0.01959687 0 Residual 4.13363420 Correlation Structure: AR(1) Formula: ~1 | dataset Parameter estimate(s): Phi 0.9759742 Fixed effects: tt ~ I(year - 1928) Value Std.Error DF t-value p-value (Intercept) 36.03912 2.5889472 154 13.920375 0.0000 I(year - 1928) 0.07954 0.0412533 154 1.928117 0.0557 Correlation: (Intr) I(year - 1928) 0.006 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -1.80191863 -0.64351311 -0.04468668 0.35374055 0.92553970 Number of Observations: 156 Number of Groups: 1 ################################### arima.fit=arima(fm3$residuals[,1],order=c(1,0,0)) arima.fit # ar1 intercept # 0.9475 -0.3574 #s.e. 0.0238 1.2385 #sigma^2 estimated as 0.797: log likelihood = -204.79, aic = 415.58 dwtest(arima.fit$residuals~1) #DW = 2.0987, p-value = 0.7325