Regression - Testing for autocorrelation in the presence of heteroscedasticity
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 &lt;- data.frame(y=rnorm(100), x1=rnorm(100), x2=rnorm(100))
ols &lt;- lm(y~x1+x2, data=DF)
DF$e &lt;- ols$resid
DF$e1 &lt;- c(NA,DF$e[-100])    # are there better ways to lag a variable?
DF$e2 &lt;- c(NA,DF$e1[-100])
aux &lt;- 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 &lt;- 1000
n &lt;- 400

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

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

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

  ols &lt;- lm(y~x1+x2+y1)

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

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


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