# Spatial Lag Model and Heteroscedasticity

I am using Spatial Lag Models with the form yi = ρWyi + βXi + εi, and am estimating these in R using spdep::lagsarlm. However, Breusch-Pagan-Tests using spdep::bptest.sarlm indicate the presence of heteroscedasticity.

Therefore my question is, whether the asymptotic standard errors of the coefficients are robust against heteroscedasticity? As it is stated in the spdep manual (p. 87f.), spdep::lagsarlm uses Generalized Least Squares, which should be already consistent under heteroscedasticity, right?

But why is it stated in the spdep manual for spdep::bptest.sarlm(p. 17) then:

It is also technically possible to make heteroskedasticity corrections to standard error estimates by using the “lm.target” component of sarlm objects - using functions in the lmtest and sandwich packages.

And if there are indeed corrections for standard errors required, how can I do this in the spdep package?