Regression - Testing for autocorrelation in the presence of heteroscedasticity

Question

I have constructed a linear time series regression model and estimated the parameters by applying OLS. I now want to test wether the assumptions for proper large sample inference (asymptotic Gauß Markov assumptions) are fulfilled.

Now, I am not sure how to test wether the residuals are autocorrelated or not. Since my model contains lagged dependent variables I can not use the Durbin-Watson test (since my independent variables are not strictly exogenous). Following Wooldridge I decided to apply the Breusch-Godfrey test. But the residuals are heteroscedastic, which I tested for via applying the Bresuch-Pagan test.

Wooldridge says that in the case of heteroscedasticity, one can not apply the usual Breusch-Godfrey test. How can I test for autocorrelation in the presence of heteroscedasticity? Is there any robust method? If that is of any interest - I am using R, so it would be helpful if there would be an implementation of the method (if there is one) in R.

EDIT: I have found a quite interesting paper that proposes a method of dealing with the topic: The modified Breusch-Godfrey test. Link: http://www.naun.org/main/NAUN/mcs/17-542.pdf.

Yet, I did not find any practical implementations of this test. As I am (just) an undergrad student, my possibilities regarding implementing such methods on my own are rather limited. So I am still looking for a general approach/test or method. (And I assume there has to to be a method, because the problem I am having strikes me as a rather common one.) Thank you!

Answer 1

General remarks: The BG test under homoskedasticity can be done using the bgtest command in the lmtest package of R. The $(n-p)R_{aux}^2$ version mentioned in link works only under homoskedasticity. In the presence of heteroskedasticity, Wooldridge (1991, JoE) gives a discussion (as noted in the Wooldridge textbook you mentioned).

What I think: I guess that what Wooldridge does is to use a heteroskedasticity-robust variance estimator. For this, (i) get the OLS residuals, (ii) regress e(t) on e(t-1), ..., e(t-p) and X, and test the joint significance of e(t-1), ..., e(t-p) using a heteroskedasticity-robust covariance estimate. If you want to use R, do the following for AR(2):

DF <- data.frame(y=rnorm(100), x1=rnorm(100), x2=rnorm(100))
ols <- lm(y~x1+x2, data=DF)
DF$e <- ols$resid
DF$e1 <- c(NA,DF$e[-100])    # are there better ways to lag a variable?
DF$e2 <- c(NA,DF$e1[-100])
aux <- lm(e~e1+e2+x1+x2, data=DF)
library(car)
lht(aux, c('e1','e2'), white.adjust='hc3')

Discussions: That said, there is the generated regressor problem in the aux regression above, that is, the some of RHS variables (e1 and e2) are generated using the OLS regression results. This usually causes trouble. However some tests can be done even if the regressors are generated. I guess this is one, but I haven't checked it formally.

Further discussions: The BG test is an LM test, while lht does a Wald test. The difference should be minor.

Simulation results: I did simulations. The ordinary BG test seems to fail. The robustified version seems working.

library(car)
iterate <- 1000
n <- 400

ans <- data.frame(ord=rep(NA,iterate), rob=rep(NA,iterate))
set.seed(1)
for (iter in seq_len(iterate)) {
  x1 <- rnorm(n+1)
  x2 <- rnorm(n+1)
  u0 <- rnorm(n+1)
  u0[1:floor(n/2)] <- 2*u0[1:floor(n/2)]
  u <- sqrt(1+abs(x1+x2))*u0

  ## y(t) = 1+x1(t)+x2(t)+0.5*y(t-1)+u(t)
  y <- filter(1+x1+x2+u, 0.5, method='recursive')

  y1 <- y[-(n+1)]      # y(1), ..., y(n)
  y <- y[-1]           # y(2), ..., y(n+1)
  x1 <- x1[-1]         # x1(2), ..., x1(n+1)
  x2 <- x2[-1]         # x2(2), ..., x2(n+1)

  ols <- lm(y~x1+x2+y1)

  e <- ols$resid
  e1 <- c(NA,e[-n])    # NA, e(1), ..., e(n-1)
  e2 <- c(NA,e1[-n])   # NA, NA, e(1), ..., e(n-2)
  aux <- lm(e~e1+e2+x1+x2+y1)
  tst0 <- lht(aux, c('e1','e2'), white.adjust=FALSE)
  ans$ord[iter] <- as.numeric(tst0$`Pr(>F)`[2] < 0.05)
  tst1 <- lht(aux, c('e1','e2'), white.adjust=TRUE)
  ans$rob[iter] <- as.numeric(tst1$`Pr(>F)`[2] < 0.05)
}

print(colMeans(ans))
##   ord   rob
## 0.090 0.046

Note: I edited the simulation a lot. I hope this one is OK.