I haven't seen a reply to this, so I will offer a comment:
It looks to me like the correlation you question is the correlation
between the estimated intercept and slope of the gls regression line.
This is different from the correlation between ts.mar and ts.anr. To
understand this, consider the following addition to your script:
plot(ts.mar, ts.anr)
# The intercept and slope are negative
# because the mean of ts.mar is positive.
fit3a <- gls(ts.anr ~ I(ts.mar-1000),correlation =
corARMA(value=c(mod3$coef[1],mod3$coef[2]),p=2))
summary(fit3a)
<snip>
Correlation:
(Intr)
I(ts.mar - 1000) 0.988
# The estimated slope is the same
# in fit3 and fit3a
# but the sign of the correlation between
# intercept and slope has changed.
I hope this response still interests you, even though it comes well
over a month after your post. You are to commended for providing a
good, reasonably complete example. If it had been shorter, it might
have received more comments sooner.
If this does not answer your question or you have another, please let
us know.
spencer graves
Michael Tiemann wrote:
Dear list,
I am trying to re-analyse something. I do have two time series, one
of which (ts.mar) might help explaining the other (ts.anr). In the
original analysis, no-one seems to have cared about the data being
time-series and they just did OLS. This yielded a strong positive
correlation.
I want to know if this correlation is still as strong when the
autocorrelations are taken into account. There are autocorrelations, so
I model the data with arima() to get the parameters and fit it with
gls(). So far, the code seems to work fine, but what puzzles me is that
I get different sings: the gls-fit yields a strong negative correlation.
This shouldn't be so, so I suspect I am doing something wrong.
Here is my code:
# this is my data
ts.mar<-ts(c(431.3,438,389.7,353.3,354.6,371.8,397.7,438.5,467.9,505.7,574.7,644.7,667.8,616.4,509.6,447,413.1,384.1),start=1980,freq=1)
ts.anr<-ts(c(104.1,102.4,97.9,96.2,95.1,95.1,97.9,101.6,105.9,111.1,117.9,121.3,121.8,114.2,107.6,105.1,101.9,98.6),start=1980,freq=1)
# to find autocorrelations via (p)acf's and mle I do:
fun.tsa.mle<-function(x){
par(mfrow=c(3,1))
acf(x)
pacf(x)
# AR model is estimated
m1<- ar.mle(x)
# An estimation of the unexplained portion of variance
m1.1<-m1$var.pred
# plot the function
plot(x)
# Give a printout
print(m1)
print("unexplained portion of variance:")
print(m1.1)
print("Mean:")
print(m1$x.mean)
par(mfrow=c(1,1))
}
#now, the autocorrelations should be consistent with following processes:
fun.tsa.mle(ts.mar) #following DAAG a p=2 AR
fun.tsa.mle(ts.anr) #following DAAG a p=2 AR
#I need to know, wether ts.anr can be explained with ts.mar, so
#according to ar.mle:
mod3<-arima(ts.anr,order=c(2,0,0),xreg=ts.mar,transform.pars=TRUE)
fit3 <- gls(ts.anr ~ ts.mar,correlation =
corARMA(value=c(mod3$coef[1],mod3$coef[2]),p=2))
summary(fit3)
ts.plot(ts.anr,fit3$fitted,col=1:2)
#the puzzling bit is the negative correlation. It ought to be positive,
I think.
#a simple OLS (this is what the people before me have done) yields
test3<-ols(ts.anr~ts.mar)
test3 #with a positive correlation. Why?
Where is the mistake? Up to now, I just thought time-series analyses
would correct parameters and estimations, but simply changing signs?
Appreciating your help and suggestions,
Michael.
------------------------------------------------------------------------
______________________________________________ R-help at stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Spencer Graves, PhD Senior Development Engineer PDF Solutions, Inc. 333 West San Carlos Street Suite 700 San Jose, CA 95110, USA spencer.graves at pdf.com www.pdf.com <http://www.pdf.com> Tel: 408-938-4420 Fax: 408-280-7915