[R-sig-Geo] Adjusting for HAC in splm models

[Denys Dukhovnov]

Dear community,

I am trying to interpret and report the results of my spatial panel analysis, but I am running into issues while attempting to accommodate for spatial and serial autocorrelation and clustering in the errors.

What implemented solutions currently exist in R for heteroskedasticity and autocorrelation correction (HAC) and/or standard error clustering for spatial panel models? Beyond the robust that could be run with non-spatial plm, fixest, or lfe packages, I have not been able to find an analogous HAC procedure for spatial panel models. The typical lmtest::coeftest() with sandwich vcov argument do not work directly on splm objects. 

A typical error on application of lmtest:coeftest(splm.obj, .vcov = vcovHAC(splm.obj)) is this:
Error in UseMethod("estfun") : no applicable method for 'estfun' applied to an object of class "c('splm_ML', 'splm', 'splm_ML')"

However, the splm documentation says that splm::vcov.splm() function is able to extract vcov matrix for interoperability with lmtest functions. How would I go about doing it?

I consulted a fairly recent paper on software for spatial panel analysis by Bivand, Millo & Piras (2021), but there is little in the way of discussion on how to adjust the standard errors for spatial and serial autocorrelation and clustering.

I am also aware of Conley spatial errors (Conley (1999)), implemented in fixest package for example, but these tend to require latitude and longitude as inputs (which is fine, if my inputs were points). But my dataset consists of polygons of various sizes where reducing them to a centroid to compute a distance matrix becomes theoretically untenable. I would prefer to stick with adjacency/contiguity weights.

I am running R version 4.2.2 (splm v.1.6.3). Below is a simple reproducible example with similar outcomes:

library(plm)
library(spdep)
library(splm)
library(sandwich)
library(lmtest)

# Read in the data and spatial weights
data("RiceFarms")
data("riceww")

# Convert into panel data frame and convert spatial weights into list form
RiceFarms <- pdata.frame(RiceFarms, index = c("id", "time"))
listw <- mat2listw(riceww)

# Pre-compute lag for variables for use in Durbin specifications
RiceFarms$slag.pesticide <- slag(RiceFarms$pesticide, listw = listw)
RiceFarms$slag.wage <- slag(RiceFarms$wage, listw = listw)

# PLM model
plm.mod <- plm(goutput ~ pesticide + wage,data = RiceFarms, 
               model = "within", 
               effect = "twoways", 
               index = c("id", "time"))

# SDEM model
sdem.mod <- spml(goutput ~ pesticide + wage + slag.pesticide + slag.wage,
                 data = RiceFarms, 
                 listw = listw,
                 spatial.error = "b",
                 lag = FALSE,
                 index = c("id", "time"), 
                 model = "within", 
                 effect = "twoways")

# Apply robust SE corrections
lmtest::coeftest(plm.mod, vcov. = vcovHC(plm.mod, method = "arellano", type = "HC0"))

# Neither works with splm object
lmtest::coeftest(sdem.mod, vcov. = vcovHC(sdem.mod, method = "arellano", type = "HC0"))
lmtest::coeftest(sdem.mod, vcov. = sandwich::vcovHAC(sdem.mod))


A follow-up question is how to deal with HAC standard error corrections for the impacts, given that spatial panel impacts function are not directly implemented in splm as of version 1.6.3, although the coefficients can be computed as per the previous posts on the subject:

https://stat.ethz.ch/pipermail/r-sig-geo/2016-April/024326.html
https://stat.ethz.ch/pipermail/r-sig-geo/2019-July/027513.html

Any help would be much appreciated.


Best,
Denys Dukhovnov