##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")


  ##2. CALCULATE HADSST TREND

	#reported 0.038  ± 0.011 deg C per decade
	#R2=0.54 R1=0.85 DW=2.20

	library(nlme)
	Data=groupedData(tt~year|dataset,data=Data)
	temp=(Data$dataset=="hadsst")&(Data$year<=2005)

     #simple linear model
	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

     #lme model with "ML"	
	fm1=lme(tt~year,data=Data[temp,],method="REML"); 
	coef(fm1)
    	#      (Intercept)        year
	#hadsst   -6.689188 0.003400966


     #assume reported fit
	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)


  ##DURBIN-WATSON TEST ON TREND RESIDUALS
	dwtest(fm)  #with fresh calculation
		#DW = 0.4926, p-value < 2.2e-16

	dwtest(test~1)  #assuming reported trend
		#DW = 0.2704, p-value < 2.2e-16


  ##ARIMA FIT TO RESIDUALS
	arima.fit=arima(fm$residuals,order=c(1,0,1));arima.fit  #with fresh fitted residuals

	#Coefficients:
        #	 ar1     ma1  intercept
	#      0.6571  0.2165     0.0049
	#s.e.  0.0964  0.1475     0.0311
	#sigma^2 estimated as 0.01228:  log likelihood = 121.4,  aic = -234.81


	arima.fit2=arima(fm2$residuals,order=c(1,0,1));arima.fit2 #with fresh fitted residuals
	#Coefficients:
 	#       	 ar1     ma1  intercept
	#      0.6624  0.2113     0.0050
	#s.e.  0.0954  0.1481     0.0315
	#sigma^2 estimated as 0.01228:  log likelihood = 121.37,  aic = -234.73

   ##DURBIN-WATSON TEST ON ARIMA.FIT RESIDUALS
	dwtest(arima.fit$residuals~1)  #with fresh calculation
	#	DW = 2.0433, p-value = 0.6072
	#alternative hypothesis: true autocorrelation is greater than 0 

	dwtest(arima.fit2$residuals~1)  #assuming reported trend
	#DW = 2.0435, p-value = 0.6078
	#alternative hypothesis: true autocorrelation is greater than 0 




##geoR IMPLementation

	library(geoR)

	fm.ml <- likfit(Data[temp,], trend="1st", ini=c(1.5,1.2), fix.nug=TRUE)   # default is ML estimation