# Peter Guttorp 5/14/2013, using and modifying to analyze VR09
# All code by Nicolas Nierenberg 4/6/2009 except where remarked by PG
# Updated 4/13/2009 comments to reflect changes to use new SSA
# Also added new analysis to the end.
# Analyzing Rahmstorf 2007 (R07), and taking a look at an alternative
# model of sea level rise.
# In his response to comments on his original paper Dr. Rahmstorf
# included the data file that he used as well as the matlab script.
# Unfortunately the matlab script uses a function called ssatrend which
# I can't find anywhere on the internet, nor does there appear to be an
# equivalent in R.
# Based on a comment from Dr. Rahmstorf on the internet I wrote Aslak Grinsted
# and he sent me an unpublished version of ssatrend.m which I have now translated
# into R
# In his response to comments Dr. Rahmstorf also supplied the Church and White
# data set he used in his paper. I didn't find this exact data set at http://www.cmar.csiro.au/sealevel/sl_data_cmar.html
# which appears to be their primary data source. It is an annual data set, but
# it appears to be very similar to the monthly data set which appears on that site.
# In the following code I try to stick as closely as possible to the Rahmstorf
# matlab script "sealevel.m"
# Install ssa package
install.packages("simsalabim",repos="http://R-Forge.R-project.org")
library("simsalabim")
install.packages("matlab")
library("matlab")
library("nnet")
# This function is an R rewrite of ssatrend by Aslak Grinsted, but not as complete
ssatrend = function(x,window) {
x.decomp = decompSSA(x,window)
# Now we search the eigenvectors for the one representing the trend.
e = x.decomp$U
indx = which.is.max(abs(sum(e))/sum(abs(e)))
# Just need the one vector representing the trend
e = e[,indx]
# Now create a filter using this vector
efilt = convolve(e,flipud(e),type="o")
# Normalize
efilt = efilt/sum(efilt)
# Now we have to pad the original series so that we can get a full trend.
# We do this by adding values to the start and end based on the linear trend
# of the first window of values.
n = length(x)
df = data.frame(x=x[1:window],time=seq(1,window))
res = lm(x~time,data=df)
left.df = data.frame(time=seq(1-window,0))
x.pad = c(predict.lm(res,left.df),x)
df = data.frame(x=x[(n-window+1):n],time=seq(1,window))
res = lm(x~time,data=df)
right.df = data.frame(time=seq(window+1,2*window))
x.pad = c(x.pad,predict.lm(res,right.df))
# Now just run it through the filter
result = filter(x.pad,efilt)
result = result[(window+1):(length(result)-window)]
result
}
# Let's begin by loading the Rahmstorf supplied data
church_13221 = read.table("church_13221.txt")
seayear = church_13221$V1
sealevel = church_13221$V2/10
giss_landocean = read.table("giss_landocean.txt")
gissyear=giss_landocean$V1
gisstemp=giss_landocean$V2 * .01
# Use ssatrend to do the smoothing to as in the original paper
window = 15
smoothsealevel = ssatrend(sealevel,window)
smoothgisstemp = ssatrend(gisstemp,window)
# For some reason he rescales to 1990 sea level.
g1990=(smoothsealevel[120]+smoothsealevel[121])/120
sealevel = sealevel-g1990
smoothsealevel = smoothsealevel-g1990
rateofrise = diff(smoothsealevel)
#PG: Taking out the binning
# # Now he bins the data rather than just running a regression against the smoothed data
# # He uses the indexes to sort of align the years between Sea level which begins in 1870
# # and temperature which begins in 1880
# # mtemp = numeric()
# mrate = numeric()
# myear = numeric()
# for(i in 1:24) {
# mtemp[i] = mean(smoothgisstemp[(5*i-2):(5*i+2)])
# mrate[i] = mean(rateofrise[(5*i+7):(5*i+11)])
# myear[i] = mean(gissyear[(5*i-2):(5*i+2)]) # I guess this is for plotting
# }
# # Rahmstorf reports the correlation, which is .88, note that with a seven year embedding this
# # drops to .68, at 17 years correlation is at .9. Above 17 years the correlation stays about the same
# # but the rate starts to drop.
#
# # Now we are ready to compute the regression, note that the p-value is very similar to R07
# # at 1.3 x 10^-8
#PG:Doing a regular OLS regression to start with
#PG aligning series
mrate=rateofrise[10:132]
mtemp=smoothgisstemp[1:122]
temprate=diff(smoothgisstemp[1:122])
temprate=c(-0.036,temprate) #extrapolated value
r07.ols = lm(mrate~mtemp)
summary(r07.ols)
plot(r07.ols$residuals,temprate)#PG: Look at relation to dT//dt
# Plot equivalent to Fig. 2 in Rahmstorf 2007
plot(NULL,xlim=c(-.5,.5),ylim=c(0,3.5),type="n",col="black",xlab="Warming above 1951-1980 mean",ylab="Rate of Sea Level Change (mm/year)")
points(mtemp,mrate*10,col="red")
# I cheated and just put in the values from the regression rescaled for the graph
# to draw the line.
lines(c(-.5,.5),c(-.001721,.5*3.39748+1.68153),col="black")
# Now just regress the full delta sea level with full smooth temp, the results are basically the
# same, but it sets up the plot.
res = lm(rateofrise[12:131]~smoothgisstemp[3:122])#PG:Why start at year 3?
summary(res)
plot(NULL,xlim=c(1870,2001),ylim=c(0,4),type="n",col="black",xlab="Year",ylab="Rate of Change (mm/yr)")
# PG: I am now going on to implement VR09 and skipping the rest of Nierenberg's analysis
# lines(gissyear[3:122],10*fitted(res))
# lines(gissyear[3:122],10*ssatrend(rateofrise[12:131],window),col="red")
# # Naively compare a model based on the first 12 bins with a model based on the second 12 bins.
# res = lm(mrate[1:12]~mtemp[1:12]) # coefficient is .44 versus .34 in the full model
# summary(res)
# res = lm(mrate[13:24]~mtemp[13:24])
# summary(res)
# # The issue with the simple test above is that due to smoothing the first 12 bins are based
# # partially on data in the second period. (Likewise for the second 12 bins). This was
# # noted in a technical correction posted by Dr. Rahmstorf in Science in October 2008, more than a year
# # after his original response.
# # So to do this correctly we need to break the data sets up in to two pieces, and then process each
# # piece separately. This is where the date ranges get interesting. It turns out it matters
# # that you start with 1882 instead of 1880. This range is different was noted in the technical correction
# # but it is different than the 1880 to 1940 stated in the original response. I note however that the
# # way that the binning was done, I believe that both the original paper and the response actually start
# # in 1882. That made no difference to the full period results.
# # Another important point is that the reported result in the technical correction can only be achieved
# # by reducing the window size to 10 years. It makes sense to reduce the window with the smaller data
# # set, but this reduction was not mentioned in the tehcnical correction.
# # The technical correction also doesn't report on the result for the second half of the data, using the
# # split data sets, the result is very different.
# # These choices yield the same result as the techincal correction
# window=10
# gissrange.per1=3:62
# gissrange.per2=63:122
# sealevelrange.per1=12:72
# sealevelrange.per2=72:132
# # Moving the first period back one year changes the coefficient to .26
# #gissrange.per1=2:61
# #sealevelrange.per1=11:71
# # Moving the first period forward one year changes the coefficient to .41
# #gissrange.per1=4:62
# #sealevelrange.per1=13:73
# smoothgisstemp.per1=ssatrend(gisstemp[3:62],10)
# smoothgisstemp.per2=ssatrend(gisstemp[63:122],10)
# smoothsealevel.per1=ssatrend(sealevel[12:72],10)
# smoothsealevel.per2=ssatrend(sealevel[72:132],10)
# rateofrise.per1=diff(smoothsealevel.per1)
# rateofrise.per2=diff(smoothsealevel.per2)
# res = lm(rateofrise.per1~smoothgisstemp.per1)
# summary(res) # The first period coefficient is .35 As reported
# res = lm(rateofrise.per2~smoothgisstemp.per2[2:61])
# summary(res) # The second period coefficient is .25 As not reported
# # Just to show that binning makes no difference once again.
# for(i in 1:12) {
# mtemp[i] = mean(smoothgisstemp.per1[(5*i-4):(5*i)])
# mrate[i] = mean(rateofrise.per1[(5*i-4):(5*i)])
# }
# res = lm(mrate[1:12]~mtemp[1:12])
# summary(res) # coefficient is .35, as reported.
# # Let's try that for the second period
# for(i in 1:12) {
# mtemp[i] = mean(smoothgisstemp.per2[(5*i-4):(5*i)])
# mrate[i] = mean(rateofrise.per2[(5*i-4):(5*i)])
# }
# res = lm(mrate[1:12]~mtemp[1:12])
# summary(res) # coefficient is .25, not reported.
# # Make a quick estimate of sea level rise given a 3C temperature increase
# # define f as the linear equation from the model (times ten to make it millimeters/year)
# f = function(x)(x*.339748+.168153)*10
# # The beginning rate is defined by the temperature in the last five year period, which is close enough
# # to the satellite observed rate at 3.1. The ending temperature is defined as the starting temperature
# # plus 3. I average the rate and multiply by 100 for the number of years (just multiply by 50)
# (f(3.31)+f(mrate[24]))*50 # This comes out to 78 cm
# # If you want to use 4C then
# (f(4.31)+f(mrate[24]))*50 # It is 95cm
#PG: Back again. Here we go. Everything from now on is PG.
vr.ols=lm(mrate~mtemp+temprate)
summary(vr.ols)
#First plot the smooths, the data, and the VR-fit
plot(gissyear[1:122],smoothsealevel[11:132],type="l")
lines(gissyear[1:122],7.477612+church.smooth,col="red")#line up with smoothsealevel[1]
lines(gissyear[1:122],smoothsealevel[11]+cumsum(0.177011+0.316575*mtemp-1.464321*temprate),col="darkgreen")
lines(gissyear[1:122],7.477612+church[-(123:128),2])
text(1890,10,"SSA smooth")
text(1890,9,"Lowess",col="red")
text(1890,8,"VR reg",col="darkgreen")
library(forecast)
auto.arima(vr.ols$residuals)
#We get an ARMA(3,2)-model
library(nlme)
vr.gls=gls(mrate~mtemp+temprate,correlation=corARMA(p=3,q=2),method="ML")
summary(vr.gls)
#temprate has changed sign relative to OLS
#and is no longer significant
lines(gissyear[1:122],smoothsealevel[11]+cumsum(0.1630910+0.2724587*mtemp+0.3688486*temprate),col="darkblue")
text(1890,7,"VR GLS",col="darkblue")
plot(gissyear[1:122],smoothsealevel[11:132]-smoothsealevel[11]-cumsum(0.177011+0.316575*mtemp-1.464321*temprate),type="l")
lines(gissyear[1:122],smoothsealevel[11:132]-7.477612-church.smooth,col="red")
abline(h=0)