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