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

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.

• What exactly do you believe to be biased when using OLS? – Christoph Hanck Apr 5 '16 at 4:43
• 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. – user111024 Apr 5 '16 at 17:09
• Edited the original post. – user111024 Apr 5 '16 at 17:53
• OK, I see now - see my attempt below. – Christoph Hanck Apr 6 '16 at 12:20

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

• 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. – user111024 Apr 16 '16 at 18:33