Example of OLS vs GLS with AR1 residuals for teaching in R

Question

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.

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.)

library(nlme)

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