##############################################################################################################
##
##	This is R code to compute trends in Temperature data from Hadley & RSS, to determine the pause length.
##	It generates results for McKitrick (2014) "Robust Measurement of the Duration of the Global Warming 
##	Hiatus" Open Journal of Statistics"
##	
## 
##############################################################################################################

rm(list=ls(all=TRUE))		# remove all objects in memory
require(orcutt)			# for AR1 estimation
require(nlme)
require(fUnitRoots)		# for UR testing


##############################################################################################################
# FUNCTION to take 5-col matrix with year, month, global, nh and sh series, and startyear, and returns a table 
# showing, from startyear to T-5, VF CI lower bound, trend, upper bound, p(VFg=0), p(VHn=VFs=0)
#
# The function is called "vfpause"
##############################################################################################################
vfpause = function(x , s) {

	global = x[,3]					# extract columns of input matrix
	hemi   = x[,4:5]
	year   = x[,1]
	month  = x[,2]

	startyear=s					# starting year
	endyear=max(year)				# endyear of sample
	yy=endyear-s+1					# No. years in trend sample
	mm=month[length(month)]				# Final month (number)

#	COMPUTE TRENDS					# starting 5 years back
	table = matrix(0, nrow = (yy-1), ncol = 7)	

	for (y in 5:(yy-1))	{
		start=endyear-y
		Ty=12*y+mm					# length of time trend in subsample
		d = c(1, rep( 0 , times=Ty-1))			# d-vector
		t=c(1:Ty)					# simple trend
		beta1 = lm(t ~ 1)				# partial out the t-terms
		ttilde1 = as.matrix(beta1$residuals)
		eta1 = sum( ttilde1^2 )				# partialling out term (aka inv(X'X) )
#	
		b = lm(global[year>=start]*120 ~ t)		# OLS trend y years back in global series
		hz = matrix(hemi[year>=start]*120, nrow = Ty, ncol = 2)
		bh= lm(hz ~ t)					# OLS trends y years back in NH, SH
		bh_t=bh$coefficients[2,]
		bh_t=as.matrix(bh_t)

#		OMEGA MATRIX & VF SCORE FOR GLOBAL TREND
		u = t(b$residuals)				# transpose resid vector
		S= u[1] 					# S = first col of u 
		OMEGA1 = S*S					# First accum entry for Omega matrix
			for (i in 2:Ty)		{			# Loop to accumulate other Omega=SS' terms
			d[i]=1					# Vector to add 1st i cols of u (sets next element to 1)
			S = u %*% d				# forms vector with sum of 1st i cols of u
			OMEGA1 = OMEGA1 + S*S			# accumulates partial sums
						}
		OMEGA1 = 2*OMEGA1 / (Ty^2)			# Last step
		VF1 = eta1 * ( b$coefficients[2]^2 /OMEGA1 )
		ci1 = 6.482*abs( b$coefficients[2]/sqrt(VF1) )	# 95% Confidence interval width

#		OMEGA MATRIX & VF SCORE FOR HEMISPHERE TRENDS
		u = t(bh$residuals)				# transpose resid vector
		S= u[1:2,1] 					# S = first col of u 
		d = c(1, rep( 0 , times=Ty-1))			# d-vector
		OMEGA2 = S %*% t(S)				# First accum entry for Omega matrix
			for (i in 2:Ty)		{			# Loop to accumulate other Omega=SS' terms
			d[i]=1					# Vector to add 1st i cols of u (sets next element to 1)
			S = u %*% d				# forms Nx1 vector with sum of 1st i cols of u
			OMEGA2 = OMEGA2 + S %*% t(S)		# accumulates partial sums
					}
		OMEGA2 = 2*OMEGA2 / (Ty^2)			# Last step
		R=diag(2)		 			# make R = identity matrix to test joint = 0
		VF2 = eta1 * t(R %*% bh_t) %*% (solve( OMEGA2  ) ) %*% (R %*% bh_t) / 2

#		MONTHLY TREND CI's & VF SCORES
		table[y,1]=y
		table[y,2]=start
		table[y,3]=b$coefficients[2]-ci1
		table[y,4]=b$coefficients[2]
		table[y,5]=b$coefficients[2]+ci1
		table[y,6]=VF1
		table[y,7]=VF2
			}	

	table = round(table,4)
		t_zero = table[,2]
		trend =  table[,4]
		low =    table[,3]
		upp =    table[,5]

#	Display results
	return(table)

#	END OF VF PAUSE FUNCTION
	}



##############################################################################################################
# Define a function that takes 5-col matrix with year, month, global, nh and sh series, and startyear, and returns a table 
# showing, from startyear to T-5, AR1 CI lower bound, trend, upper bound, sh lower, nh lower
#
# The function is called "ar1pause"
##############################################################################################################
ar1pause = function(x , s) {

	global = x[,3]					# extract columns of input matrix
	nh     = x[,4]
	sh     = x[,5]
	year   = x[,1]
	month  = x[,2]

	startyear=s					# starting year
	endyear=max(year)				# endyear of sample
	yy=endyear-s+1					# No. years in trend sample
	mm=month[length(month)]				# Final month (number)

#	COMPUTE TRENDS					# starting 5 years back
	table = matrix(0, nrow = (yy-1), ncol = 7)	

	for (y in 5:(yy-1))	{
		start=endyear-y
		Ty=12*y+mm					# length of time trend in subsample
		t=c(1:Ty)					# simple trend
#	
		b = lm(global[year>=start]*120 ~ t)		# OLS trend y years back in global series
		b1= cochrane.orcutt(b)				# C-O estimate of rho
		r = b1$rho
		se= summary(b)$coef[4]*sqrt((1+r)/(1-r))	# AR1-adjusted SE
		cig = 1.96*se

		bn = lm(nh[year>=start]*120 ~ t)		# OLS trend y years back in global series
		b1= cochrane.orcutt(b)				# C-O estimate of rho
		r = b1$rho
		se= summary(bn)$coef[4]*sqrt((1+r)/(1-r))	# AR1-adjusted SE
		cign = 1.96*se

		bs = lm(global[year>=start]*120 ~ t)		# OLS trend y years back in global series
		b1= cochrane.orcutt(bs)				# C-O estimate of rho
		r = b1$rho
		se= summary(bs)$coef[4]*sqrt((1+r)/(1-r))	# AR1-adjusted SE
		cigs = 1.96*se
		
#		MONTHLY TREND CI's & VF SCORES
		table[y,1]=y
		table[y,2]=start
		table[y,3]=b$coefficients[2]-cig
		table[y,4]=b$coefficients[2]
		table[y,5]=b$coefficients[2]+cig
		table[y,6]=bn$coefficients[2]-cign
		table[y,7]=bn$coefficients[2]-cigs
			}	

	table = round(table,4)

#	Display results
	return(table)

#	END OF FUNCTION
	}









##############################################################################################################
##	UAH GLOBAL LT 
##############################################################################################################


loc="uah.txt"
ht1 = read.table(file = loc, dec=".", skip=2)		# Ignore December 1978 value and start at Jan 1979
	attach(ht1)
	UAH_g = V3
	month = V2
	year = V1
	endyear=max(year)				# endyear of sample
	tt = year+(month-1)/12				# year index for graphing
	yy=endyear-1979+1				# No. years in trend sample
	UAH_nh = V6
	UAH_sh = V9
	yy=endyear-1979+1				# No. years in trend sample
	hemi = as.matrix( data.frame(UAH_nh, UAH_sh) )
	mm=month[length(month)]				# Final month (number)
rm(ht1)	

uah = as.matrix( data.frame(year, month, UAH_g, UAH_nh, UAH_sh) )

	
pause_vfa = vfpause(uah, 1979)				# Compute pause stats using VF method
	t_zeroa = pause_vfa[,2]
	trenda =  pause_vfa[,4]
	lowa =    pause_vfa[,3]
	uppa =    pause_vfa[,5]

pause_ar1a = ar1pause(uah, 1979)			# Compute pause stats using AR1 method


#	SAVE DATA TO FILE
write(t(pause_vfa), file="logfiles/table_vfa.csv", sep=", ", ncolumns=7)

tta = year+(month-1)/12				# post-1979 year index for graphing



##############################################################################################################
##	HADLEY GLOBAL SURFACE 
##############################################################################################################

# startyear=1900		# starting year

loc="http://www.metoffice.gov.uk/hadobs/hadcrut4/data/current/time_series/HadCRUT.4.2.0.0.monthly_ns_avg.txt"
ht1 = read.table(file = loc, dec=".")
	attach(ht1)
	ht=as.matrix(ht1)
	HadCRUT = ht[,2]
	HadCRUT = as.numeric(HadCRUT)

	z = strsplit(ht[,1],"/")			# split date string into year, month
	zm = matrix(unlist(z), ncol=2, byrow=TRUE)	# convert list into 2-col matrix
	year =  as.numeric(zm[,1])			# year as number
	month = as.numeric(zm[,2])			# month as number
	tt = year+(month-1)/12				# year index for graphing
	endyear=max(year)				# endyear of sample
	yy=endyear-1900+1				# No. years in trend sample
	yyh = yy					# Preserve longer value

 
#	GET NH and SH 
loc="http://www.metoffice.gov.uk/hadobs/hadcrut4/data/current/time_series/HadCRUT.4.2.0.0.monthly_nh.txt"
ht1 = read.table(file = loc, dec=".")
	attach(ht1)
	ht=as.matrix(ht1)
	HadNH = ht[,2]
	HadNH = as.numeric(HadNH)
loc="http://www.metoffice.gov.uk/hadobs/hadcrut4/data/current/time_series/HadCRUT.4.2.0.0.monthly_sh.txt"
ht1 = read.table(file = loc, dec=".")
	attach(ht1)
	ht=as.matrix(ht1)
	HadSH = ht[,2]
	HadSH = as.numeric(HadSH)
rm(z,zm,ht1,ht)	

hadley = as.matrix( data.frame(year, month, HadCRUT, HadNH, HadSH) )

pause_vfh = vfpause(hadley, 1900)			# Compute Pause stats using VF method

t_zeroh = pause_vfh[,2]
trendh =  pause_vfh[,4]
lowh =    pause_vfh[,3]
upph =    pause_vfh[,5]

#	SAVE DATA TO FILE
write(t(pause_vfh), file="logfiles/table_vfh.csv", sep=", ", ncolumns=7)


pause_ar1h = ar1pause(hadley, 1900)			# Compute Pause stats using AR1


tth = year+(month-1)/12				# post-1900 year index for graphing






##############################################################################################################
##	RSS GLOBAL LT 
##############################################################################################################

loc="http://data.remss.com/msu/monthly_time_series/RSS_Monthly_MSU_AMSU_Channel_TLT_Anomalies_Land_and_Ocean_v03_3.txt"
ht1 = read.table(file = loc, dec=".", skip=3)
	attach(ht1)
	RSS_g = V3
	month = V2
	year = V1
	endyear=max(year)				# endyear of sample
	tt = year+(month-1)/12				# year index for graphing
	yy=endyear-1979+1				# No. years in trend sample
	RSS_nh = V10
	RSS_sh = V11
	yy=endyear-1979+1				# No. years in trend sample
	hemi = as.matrix( data.frame(RSS_nh, RSS_sh) )
	mm=month[length(month)]				# Final month (number)
rm(ht1)	

rss = as.matrix( data.frame(year, month, RSS_g, RSS_nh, RSS_sh) )
	
pause_vfr = vfpause(rss, 1979)				# Compute pause stats using VF method
pause_ar1r = ar1pause(rss, 1979)			# Compute pause stats using AR1 method

ttr = year+(month-1)/12				# post-1900 year index for graphing

t_zeror = pause_vfr[,2]
trendr =  pause_vfr[,4]
lowr =    pause_vfr[,3]
uppr =    pause_vfr[,5]

#	SAVE DATA TO FILE
write(t(pause_vfr), file="logfiles/table_vfr.csv", sep=", ", ncolumns=7)




##############################################################################################################
##	DRAW GRAPHS 
##############################################################################################################

#	TREND GRAPHS

#	HADLEY
#postscript('graphics/fig1.eps', horizontal=F)			# write figure to EPS file
	par(mar=c(10,7,4,2) )
	plot(HadCRUT~tth, xlab="Year", cex.axis=2, cex.lab=2, ylab="HadCRUT deg C anomaly", type="p", cex=.5, col=gray(.5), ylim=c(-1, 1))
	par(new="T")
		sh=lowess(y=HadCRUT, x=tth, f=0.09)
		plot(sh$y~tth, xlab="", ylab="", type="p", cex=.5, ylim=c(-1, 1), xaxt="n", yaxt="n")
#dev.off()

#	UAH
#postscript('graphics/fig2.eps', horizontal=F)			# write figure to EPS file
	par(mar=c(10,7,4,2) )
#	par(mar=c(13,6,20,1) )
	plot(UAH_g ~ tta, xlab="Year", cex.axis=2, cex.lab=2, ylab="UAH LT deg C anomaly", type="p", cex=.5, col=gray(.5),  xaxs="r", ylim=c(-1, 1))
	par(new="T")
		su=lowess(y=UAH_g, x=tta, f=0.11)
		plot(su$y ~ tta, xlab="", ylab="", type="p", cex=.5,  ylim=c(-1, 1), xaxt="n", yaxt="n")
#dev.off()

#	RSS
#postscript('graphics/fig3.eps', horizontal=F)			# write figure to EPS file
	par(mar=c(10,7,4,2) )
#	par(mar=c(13,6,20,1) )
	plot(RSS_g ~ ttr, xlab="Year", cex.axis=2, cex.lab=2, ylab="RSS LT deg C anomaly", type="p", cex=.5, col=gray(.5),  xaxs="r", ylim=c(-1, 1))
	par(new="T")
		sr=lowess(y=RSS_g, x=ttr, f=0.11)
		plot(sr$y ~ ttr, xlab="", ylab="", type="p", cex=.5,  ylim=c(-1, 1), xaxt="n", yaxt="n")
#dev.off()



#	TRUMPET GRAPHS

#	HADLEY
postscript('graphics/fig4.eps', horizontal=F)			# write figure to EPS file
#par(mar=c(5,7,2,2) )
	par(mar=c(13,6,20,1) )
	plot(trendh[5:yyh]~t_zeroh[5:yyh], cex=1, cex.axis=2, cex.lab=2, col="black", pch=19, type="p", xlab="Start Year of Trend", 
		ylab="HadCRUT Trend (C/decade)", ylim=c(-.4,.4))
		abline(h=0, lty=2)
	par(new="T")
		plot(lowh[5:yyh]~t_zeroh[5:yyh], type="l", xlab="", ylab="", ylim=c(-.4,.4), lty=1, xaxt="n", yaxt="n")
	par(new="T")
		plot(upph[5:yyh]~t_zeroh[5:yyh], type="l", xlab="", ylab="", ylim=c(-.4,.4), lty=1, xaxt="n", yaxt="n")
dev.off()

#	UAH
postscript('graphics/fig5.eps', horizontal=F)			# write figure to EPS file
#par(mar=c(5,7,2,2) )
	par(mar=c(13,6,20,2) )
	plot(trenda[5:yy]~t_zeroa[5:yy], cex=1.5, cex.axis=2, cex.lab=2, col="black", pch=19, type="p", xlab="Start Year of Trend", 
		ylab="UAH LT Trend (C/decade)", ylim=c(-.4,.4))
		abline(h=0, lty=2)
	par(new="T")
		plot(lowa[5:yy]~t_zeroa[5:yy], type="l", xlab="", ylab="", ylim=c(-.4,.4), lty=1, lwd=2, xaxt="n", yaxt="n")
	par(new="T")
		plot(uppa[5:yy]~t_zeroa[5:yy], type="l", xlab="", ylab="", ylim=c(-.4,.4), lty=1, lwd=2, xaxt="n", yaxt="n")
dev.off()
#	RSS
postscript('graphics/fig6.eps', horizontal=F)			# write figure to EPS file
#par(mar=c(5,7,2,2) )
	par(mar=c(13,6,20,2) )
	plot(trendr[5:yy]~t_zeror[5:yy], cex=1.5, cex.axis=2, cex.lab=2, col="black", pch=19, type="p", xlab="Start Year of Trend", 
		ylab="RSS LT Trend (C/decade)", ylim=c(-.4,.4))
		abline(h=0, lty=2)
	par(new="T")
		plot(lowr[5:yy]~t_zeror[5:yy], type="l", xlab="", ylab="", ylim=c(-.4,.4), lty=1, lwd=2, xaxt="n", yaxt="n")
	par(new="T")
		plot(uppr[5:yy]~t_zeror[5:yy], type="l", xlab="", ylab="", ylim=c(-.4,.4), lty=1, lwd=2, xaxt="n", yaxt="n")

dev.off()


Hadley	= data.frame(tth, HadCRUT, sh$y)
  write.table(Hadley, file = "fi/Figure1.csv", sep=",", row.names=FALSE)
RSS 	= data.frame(ttr, RSS_g, sr$y)
  write.table(RSS, file = "fi/Figure2.csv", sep=",", row.names=FALSE)

Hadley	= data.frame(lowh[5:yyh], trendh[5:yyh], upph[5:yyh], t_zeroh[5:yyh])
  write.table(Hadley, file = "fi/Figure3.csv", sep=",", row.names=FALSE)
RSS 	= data.frame(lowr[5:yy], trendr[5:yy], uppr[5:yy], t_zeror[5:yy])
  write.table(RSS, file = "fi/Figure4.csv", sep=",", row.names=FALSE)