I'm looking for an example to show my class. We are covering OLS vs GLS with autocorrelated errors -- I've got the class to the point where they understand (some of them) why the the standard errors on the coefficient estimates tend to be underestimated (and the t-scores overestimated) when the autocorrelations of the errors at low lags are positive. But it would be nice to have an example where they can see how inference might change based on fitting a ``better'' model with gls.

Thus, I'd like to show an example where the residuals from an OLS model y~x have an AR1 structure that would give spurious results (std errors on the estimates, p-values, etc.). Then I would show them a gls implementation that accounts for correlated residuals. E.g.,

# example of what'd I'd like given the right x and y 
ols1 <- lm(y~x) #ols model
summary(ols1) # would show a biased model
acf(residuals(ols1)) # would show AR1 structure
# vs
gls1 <- gls(y~x,correlation = corAR1())

Any ideas out there?

Edited for clarity.

  • 1
    $\begingroup$ What exactly do you believe to be biased when using OLS? $\endgroup$ Commented Apr 5, 2016 at 4:43
  • $\begingroup$ Ooops. Sorry. Not biased -- wrong wording. But rather with high AR wikipedia tells me that ``the standard errors tend to be underestimated (and the t-scores overestimated) when the autocorrelations of the errors at low lags are positive.'' I guess I'm looking for an example where that would be the case that might impact how the regression would be interpreted. $\endgroup$
    – user111024
    Commented Apr 5, 2016 at 17:09
  • $\begingroup$ Edited the original post. $\endgroup$
    – user111024
    Commented Apr 5, 2016 at 17:53
  • $\begingroup$ OK, I see now - see my attempt below. $\endgroup$ Commented Apr 6, 2016 at 12:20

1 Answer 1


Here is a little Monte Carlo simulation that illustrates that you will reject the null of a mean of zero way too often if you neglect the AR(1) autocorrelation in the time series generated , i.e., if you estimate the variance of the mean as $s^2/n$. GLS produces rejection rates closer to the nominal level. (I am not so familiar with the nlme package, so I hope I have used it wisely.)


reps <- 1000
n <- 300

ols.rejections <- gls.rejections <- rep(NA,reps)

for (i in 1:reps){
  y <- arima.sim(n=n,list(ar=0.5))
  ols.reg <- lm(y~1)
  gls.reg <- gls(y~1,correlation = corAR1(form = ~ 1))
  ols.rejections[i] <- summary(ols.reg)$coefficients[1,4] < 0.05
      gls.rejections[i] <- summary(gls.reg)$tTable[,"p-value"] < 0.05

The result:

> mean(ols.rejections)
[1] 0.264

> mean(gls.rejections)
[1] 0.046
  • $\begingroup$ This is helpful. Thanks much. From the simulations I've done the SE on the model coef lm(y~x+w) where w is AR(1) are pretty stable when OLS and GLS are compared. It's interesting. I have to think some more about it all. $\endgroup$
    – user111024
    Commented Apr 16, 2016 at 18:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.