###################################################
##Part III - Relate Seattle level to global level##
###################################################
#note, requires part I to have been run
#now we want to examine the relationshipe between
#global sea level and seattle sea level
#first get the two series to have the same length
#NOTE: following our model, we fit the raw seattle image to the "trended"
#smooth global level.
require(forecast)
yrs=1899:2009
seas.smooth = lowess(church.sealevel,f=1/4)$y
seas.forseattle = seas.smooth[which(is.element(1880:2009,yrs))]
seattle.forseattle = seattle[which(is.element(1899:2012,yrs))]
#now plot the two series against each other
#plot(yrs,seattle.forseattle)
#lines(yrs,seas.forseattle)
#a linear relationship looks good, but we will still
#use a time series ARIMA model
ar.puget.trend = auto.arima(seattle.forseattle,xreg=as.matrix(seas.forseattle))
#now plot these results
#png("seattlevsea.png",height=350,width=700)
#par(mfrow=c(1,2))
#plot(seas.forseattle,seattle.forseattle)
#plot(ar.puget.trend$residuals)
#dev.off()