############################################# ##Part II - Relate Sea level to temperature## ############################################# require(forecast) #our model is rate of sea level rise ~ beta_0 + beta_1*temperature #from Rahmstorf 2007 #first we smooth our temperature data and our sea level data temps.smooth = lowess(giss.temp,f=1/4)$y seas.smooth = lowess(church.sealevel,f=1/4)$y s1999 = seas.smooth[120] #this is used when we integrate #now get change in sea level rate.seas = diff(seas.smooth) #the differencing procedure leaves us with one fewer observation so: temps.smooth2 = temps.smooth[2:130] #first create a linear model lm1 = lm(rate.seas~temps.smooth2) #now check the residual structure #auto.arima(lm1$residuals) #should be order (2,0,1) #now fit a time series model using auto.arima armodel = auto.arima(rate.seas,xreg=as.matrix(temps.smooth2)) #armodel gives us the model relating temperature and change in sea level #we integrate change in sea level to get sea level #how good is our model? There is some strangeness in the #residual analysis, as seen here: #png("DiagSea.png") #par(mfrow=c(2,2)) #plot(armodel$residuals) #plot(temps.smooth2,rate.seas) #plot(rate.seas) #plot(temps.smooth2) #dev.off()