Conley standard errors with fastconley

Why It Matters: A Small Example

The quickest way to see the point is to simulate data where both the regressor and the error have spatial structure — the textbook case where conventional standard errors mislead — and compare the inference. We draw 1,500 locations in a box roughly 650 km across and give both x and the error a smooth spatial component with an effective range of about 150 km.

library(fastconley)
set.seed(42)

n <- 1500
d <- data.frame(lat = runif(n, 40, 46), lon = runif(n, 5, 11))

# Smooth spatial fields via an exponential-covariance Gaussian process.
rad <- pi / 180
D <- 6371 * acos(pmin(1,
  sin(d$lat * rad) %o% sin(d$lat * rad) +
  (cos(d$lat * rad) %o% cos(d$lat * rad)) * cos(outer(d$lon, d$lon, "-") * rad)))
L <- chol(exp(-D / 150) + diag(1e-8, n))
d$x <- drop(crossprod(L, rnorm(n))) + 0.5 * rnorm(n)
e   <- drop(crossprod(L, rnorm(n))) + 0.5 * rnorm(n)
d$y <- 1 + 0.5 * d$x + e

library(fixest)
est <- feols(y ~ x, data = d, demeaned = TRUE)

vcov_conley <- function(x) {
  vcovSpHAC(x, lat = "lat", lon = "lon", dist_cutoff = 150,
            kernel = "bartlett", dist_fn = "spherical",
            ncores = 2, data = d)
}

se_tab <- data.frame(
  `standard errors` = c("IID", "Heteroskedasticity-robust", "Conley (150 km)"),
  `se(x)` = sapply(list(vcov(est, vcov = "iid"),
                        vcov(est, vcov = "hetero"),
                        vcov_conley(est)),
                   function(V) sqrt(diag(V))["x"]),
  check.names = FALSE
)
knitr::kable(se_tab, digits = 4, row.names = FALSE)

The coefficient is the same in all three rows; only the uncertainty changes. Ignoring the spatial correlation here understates the standard error by a factor of about 4.2 — enough to turn an insignificant effect “significant”. This is the problem Conley standard errors exist to solve, and everything below is about computing them exactly and quickly.

Quick Start

For fixest::feols(), fit the model with demeaned = TRUE. This keeps the centered design matrix that fastconley needs after estimation. (IV fits work the same way — the stored design is the projected second-stage one, which is exactly what the 2SLS sandwich needs.) GLM fits — fepois() / feglm() — need no flag at all: the stored maximum-likelihood scores and inverse Hessian are used directly, and everything below (including lag_cutoff for panels) applies unchanged.

fixest accepts an external covariance function in its vcov argument. The most convenient pattern is to define a small function that records the coordinate columns, cutoff, kernel, and thread count, then pass that function to summary(), etable(), or even to feols() itself.

library(fixest)
library(fastconley)

est <- feols(
  outcome ~ treatment + x1 + x2 | county + year,
  data = d,
  demeaned = TRUE
)

vcov_fastconley <- function(x) {
  vcovSpHAC(
    x,
    lat = "lat",
    lon = "lon",
    dist_cutoff = 100,     # kilometers
    kernel = "uniform",
    dist_fn = "spherical",
    ncores = 8,
    data = d
  )
}

summary(est, vcov = vcov_fastconley)
etable(est, vcov = vcov_fastconley)

The same wrapper can be used at estimation time. This is useful when you want the fitted object to carry the requested standard errors immediately.

est <- feols(
  outcome ~ treatment + x1 + x2 | county + year,
  data = d,
  demeaned = TRUE,
  vcov = vcov_fastconley
)

You can also compute the covariance matrix once and pass the matrix to fixest. That is preferable when you will print or export the same model many times and do not want the covariance function to be recomputed.

V_conley <- vcov_fastconley(est)
summary(est, vcov = V_conley)
etable(est, vcov = V_conley)

For lfe, use keepCX = TRUE when fitting the model. IV/2SLS fits (the felm multi-part formula with (endog ~ instruments)) work out of the box and agree with the fixest IV path and an independent dense-R sandwich at machine precision. Other k-class estimators, including LIML, are not supported. Explicit unit and time names are matched to absorbed effects or recovered from the model data; when both names are omitted the data are one cross-sectional block (every pair within the cutoff enters), and nothing is inferred from the absorbed effects.

library(lfe)
library(fastconley)

est <- felm(
  outcome ~ treatment + x1 + x2 | county + year,
  data = d,
  keepCX = TRUE
)

V_conley <- vcovSpHAC(
  est,
  lat = "lat",
  lon = "lon",
  dist_cutoff = 100,
  kernel = "uniform",
  dist_fn = "spherical",
  ncores = 8,
  data = d
)

sqrt(diag(V_conley))

For panel data, provide the unit and time variables. Set lag_cutoff > 0 when you also want serial HAC adjustment within units. Supplying time alone still defines the spatial blocks and treats every row as its own unit; serial HAC requires unit.

est <- feols(
  y ~ treatment + x1 + x2 | unit + year,
  data = panel_data,
  demeaned = TRUE
)

vcov_shac <- function(x) {
  vcovSpHAC(
    x,
    unit = "unit",
    time = "year",
    lat = "lat",
    lon = "lon",
    dist_cutoff = 500,
    lag_cutoff = 1,
    kernel = "bartlett",
    dist_fn = "spherical",
    balanced_pnl = TRUE,
    ncores = 8,
    data = panel_data
  )
}

summary(est, vcov = vcov_shac)

Use balanced_pnl = TRUE only when every period contains the same units and each unit has the same coordinates in every period. Otherwise leave it at the default, FALSE.

Choosing Settings

The most important choice is dist_cutoff. It should come from the research design: the distance over which residual correlation is plausible enough to affect inference.

Two conveniences need no argument at all. Weighted regressions (feols(..., weights = ~w) or felm(..., weights = w)) are handled automatically — the scores and the bread pick up the weights from the fit, matching fixest’s own weighted Conley vcov. And if lat / lon are omitted, they are auto-detected from common column names (lat/latitude, lon/long/longitude/lng, case-insensitive), with a message reporting the pick.

Start with pixel = 0. If many observations share the same coordinates, this is exact and often faster. If the data are very large and the application can tolerate approximation, try a positive pixel value that is much smaller than dist_cutoff, then compare the resulting standard errors on a representative subset.

For dist_fn, haversine and spherical are exact arc-distance formulations; chord is the straight-line distance through the sphere. They agree closely at sub-1000 km cutoffs but are not the same geometry, so report which one you used. lag_cutoff is in units of the numeric time variable (periods, not years unless the time variable is years). Numeric-looking character or factor labels are parsed without collapsing gaps ("2000", "2002" stay two units apart); other labels are valid only when lag_cutoff = 0.

cutoffs <- c(50, 100, 150, 300, 600)
se_x <- sapply(cutoffs, function(cut) {
  V <- vcovSpHAC(est, lat = "lat", lon = "lon", dist_cutoff = cut,
                 kernel = "bartlett", dist_fn = "spherical",
                 ncores = 2, data = d)
  sqrt(diag(V))["x"]
})
knitr::kable(data.frame(`cutoff (km)` = cutoffs, `se(x)` = round(se_x, 4),
                        check.names = FALSE), row.names = FALSE)

Matching The Call To The Data

Most users can leave method = "auto". The explicit settings below are useful when you want the call to document the intended data layout, or when you want an early error if the data do not have the structure you expected.

A few concrete templates are often enough.

For ordinary point data:

vcov_points <- function(x) {
  vcovSpHAC(
    x,
    lat = "lat",
    lon = "lon",
    dist_cutoff = 100,
    kernel = "uniform",
    dist_fn = "spherical",
    method = "pairwise",
    ncores = 8,
    data = d
  )
}

For a regular raster:

vcov_raster <- function(x) {
  vcovSpHAC(
    x,
    lat = "cell_lat",
    lon = "cell_lon",
    dist_cutoff = 250,
    kernel = "bartlett",
    dist_fn = "spherical",
    method = "grid",
    ncores = 8,
    data = raster_data
  )
}

For a balanced panel with one serial lag:

vcov_panel <- function(x) {
  vcovSpHAC(
    x,
    unit = "unit",
    time = "year",
    lat = "lat",
    lon = "lon",
    dist_cutoff = 500,
    lag_cutoff = 1,
    kernel = "bartlett",
    dist_fn = "spherical",
    balanced_pnl = TRUE,
    ncores = 8,
    data = panel_data
  )
}

If you are unsure whether a raster is being recognized, run once with verbose = TRUE. If you want to force the regular-raster route, use method = "grid"; a failure then tells you that the coordinates are not on a regular latitude-longitude lattice.

What Is Estimated

Let \(x_i\) be the centered row of regressors for observation \(i\), and let \(\hat u_i\) be the residual. Define the score

\[ s_i = \hat u_i x_i. \]

For a cross-section, the Conley meat is

\[ \widehat\Omega = \sum_{i=1}^{n} \sum_{j=1}^{n} K\left(\frac{d_{ij}}{c}\right) 1\{d_{ij} \le c\} s_i' s_j, \]

where \(d_{ij}\) is the selected spatial distance between observations \(i\) and \(j\), \(c\) is dist_cutoff, and \(K(\cdot)\) is the chosen kernel:

\[ K_{\text{uniform}}(u) = 1, \qquad K_{\text{bartlett}}(u) = 1 - u, \qquad u = d_{ij}/c \le 1 . \]

The uniform kernel weights every within-cutoff pair equally; the bartlett kernel downweights linearly with distance, which makes a non-positive semi-definite estimate much less likely (neither spatial kernel formally guarantees one — by default negative eigenvalues are clamped, with the same semantics as fixest’s vcov_fix; pass psd_fix = FALSE to return the matrix as computed). The covariance matrix is the usual sandwich:

\[ \widehat V = (X'X)^{-1}\widehat\Omega(X'X)^{-1}. \]

For weighted (WLS) fits the scores become \(s_i = w_i \hat u_i x_i\) and the bread uses \(X'WX\), matching fixest’s weighted Conley vcov.

For panel data, the spatial part is applied within time period. If lag_cutoff is positive, fastconley adds a serial HAC term within unit:

\[ \widehat\Omega = \sum_t \sum_{i:t_i=t} \sum_{j:t_j=t} K\left(\frac{d_{ij}}{c}\right)1\{d_{ij}\le c\}s_i's_j + \sum_g \sum_t \sum_{\ell=1}^{L} \left(1-\frac{\ell}{L+1}\right) \left(s_{g,t}'s_{g,t-\ell} + s_{g,t-\ell}'s_{g,t}\right), \]

where \(L\) is lag_cutoff. In words: spatial correlation is handled among units observed in the same period, and serial correlation is handled along the same unit over time. Note the serial sum starts at \(\ell = 1\): the contemporaneous (\(\ell = 0\)) covariance already lives in the spatial term, which includes each observation’s own score product. This is why the panel object is not the same as applying a cross-sectional Conley formula to the stacked rows.

Two computational guarantees are worth stating precisely. First, the sums above are evaluated exactly: every pair with \(d_{ij} \le c\) enters with its exact kernel weight under the selected distance geometry — there is no sampling, no tapering beyond the kernel itself, and no approximate distance formula. Second, the result is deterministic: summation is organized in fixed-size chunks, so the same data and settings produce a bit-identical matrix regardless of ncores. Both properties hold for every engine (pairwise and grid), which is why method is purely a speed choice.

By default \(\widehat V\) is scaled by the \(n/(n-K)\) small-sample correction — the same one fixest applies to Conley vcovs — so out of the box the two packages’ defaults line up. Pass ssc = FALSE to report \(\widehat V\) exactly as defined above (the convention of rbluhm/conley and fastconley versions before 0.9.0).

Why It Scales

A direct Conley calculation is often described as if it required an \(n \times n\) spatial weight matrix. That becomes impossible quickly. With 1,000,000 observations, a dense double-precision matrix would require about 7.3 TiB of memory before any covariance algebra.

fastconley avoids that object. It works with the score rows that enter the sandwich estimator and only spends time on observations that can affect the answer. It also uses three common data features directly:

• If rows share coordinates, their scores can be combined exactly.
• If panel locations are fixed over time, the same spatial relationships can be reused across periods.
• If the data are a regular latitude-longitude raster, the calculation is done as a grid problem instead of a point-pair problem.

You normally do not have to choose among these paths. The default method = "auto" chooses between the pairwise and grid engines; ordinary point data use pairwise while qualifying regular rasters can use grid.

Benchmarks

The tables below report post-estimation covariance time. The model is fitted once outside the timer, so the numbers answer a specific question: after the regression is estimated, how long does it take to compute the Conley or spatial HAC covariance matrix?

The benchmark files are shipped with the package. They are not rerun when this vignette is built.

Reproducing The Benchmarks

The quick local benchmark suite is:

FASTCONLEY_BENCH_THREADS=8 FASTCONLEY_BENCH_REPS=2 \
  Rscript inst/benchmarks/run-fastconley-benchmarks.R

The large-data benchmark suite is:

FASTCONLEY_LARGE_THREADS=2,4,8,16 FASTCONLEY_LARGE_REPS=1 \
  Rscript inst/benchmarks/run-large-data-benchmarks.R

The large optional fixest comparison is disabled by default because it can take several minutes. To run it, add FASTCONLEY_LARGE_FIXEST=1.

References

Conley, T. G. (1999). GMM estimation with cross sectional dependence. Journal of Econometrics, 92(1), 1–45.

Conley, T. G. (2010). Spatial econometrics. In S. N. Durlauf & L. E. Blume (Eds.), Microeconometrics (pp. 303–313). Palgrave Macmillan.

Hsiang, S. M. (2010). Temperatures and cyclones strongly associated with economic production in the Caribbean and Central America. Proceedings of the National Academy of Sciences, 107(35), 15367–15372. (Source of the within-period spatial + within-unit serial HAC panel estimator this package implements.)

Newey, W. K., & West, K. D. (1987). A simple, positive semi-definite, heteroskedasticity and autocorrelation consistent covariance matrix. Econometrica, 55(3), 703–708.

Bergé, L. (2018). Efficient estimation of maximum likelihood models with multiple fixed-effects: the R package FENmlm. CREA Discussion Papers, 13. (The fixest package, whose vcov_conley() is benchmarked above and whose “On standard-errors” vignette is the model for the small-sample-correction discussion.)