Dear R-list,
I would like to show the implications of estimating a linear trend to
time series,
which contain significant serial correlation.
I want to demonstrate this, comparing lm() and an gls() fits, using
the LakeHuron
data set, available in R.
Now in my particular case I would like to draw confidence bands on the plot and
show that there are differences. Unfortunately, I do not know how to
construct those bands
for the gls fit object.
Does anyone know how to do that?
Here is the code:
library(nlme)
library(car)
lm.hu <- lm(LakeHuron~time(LakeHuron))
## use durbin-watson test to estimate the AR(p) model
durbin.watson(lm.hu,max.lag=7)
# I decided to use a AR(2) model
gls.hu <- gls(LakeHuron ~ time(LakeHuron), correlation=corARMA(p=2),
method="ML")
plot(LakeHuron, main="Lake level time series (Lake Huron)")
lines(as.numeric(time(LakeHuron)),predict(lm.hu,interval=c("conf"))[,1],lty=1,col=2)
lines(as.numeric(time(LakeHuron)),predict(lm.hu,interval=c("conf"))[,2],lty=2,col=2)
lines(as.numeric(time(LakeHuron)),predict(lm.hu,interval=c("conf"))[,3],lty=2,col=2)
## predict of the gls class does not contain confidence intervals
lines(as.numeric(time(LakeHuron)),predict(gls.hu,interval=c("conf")),lty=1,col=4)
Thank you,
Maik
--
Maik Renner
Institut f?r Hydrologie und Meteorologie
Technische Universit?t Dresden
