Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
60,098 characters · 25 sections · 47 citation commands
Bandwidth Selection for Spatial HAC Standard Errors
\setstretch{1.5}
The concept of spatial autocorrelation has served as an organizing principle in spatial statistics for many decades, dating back at least to Student1914's article on “the elimination of spurious correlation due to position in time and space.” Student1914 demonstrated that the application of conventional statistical methods to autocorrelated observations “would lead to altogether misleading values” (p. 179) and presented adjusted inferences as evidence of the importance of such corrections. In the modern econometrics literature, the problems for statistical inference induced by spatial autocorrelation have been recognized since at least AnselinGriffith1988, who showed that spatial effects can substantially distort regression-based inference. Given that virtually all economic data have spatial structure---and that georeferenced data are becoming ubiquitous---this is a growing concern in applied research. The importance of accounting for dependence in inference has been highlighted in many contexts: BertrandDufloMullainathan2004 demonstrated the severity of the problem in panel data with temporal correlation, and casual empiricism suggests that applied researchers are now acutely aware of the potential for distortions. Yet for spatial correlation, no equivalent practical guidance on how to operationalize the correction has emerged.
The most common approach to address spatial autocorrelation in regression analysis in economics and many other social sciences is to rely on spatial heteroskedasticity and autocorrelation consistent (HAC) standard errors, introduced in the seminal contribution of Conley1999. Conley extended the logic of NeweyWest1987 from time-series to the two-dimensional spatial domain, building on the HAC framework of Andrews1991. The core of the technique is a non-parametric estimator for the variance-covariance matrix that uses kernel-weighted averages of cross-products of residuals, where the kernel downweights pairs of observations as a function of the distance between them. A key input for operationalizing this estimator is the choice of a cutoff bandwidth---the distance beyond which residuals are assumed to be uncorrelated. In the time-series analogue, a rich literature on automated and rule-of-thumb bandwidth selection exists NeweyWest1987, Andrews1991, Sun2014. For the spatial case, however, no widely adopted equivalent exists. Instead, researchers typically rely on ad hoc choices guided by intuition or domain knowledge, creating what GelmanLoken2013 term “researcher degrees of freedom.”
This paper makes three contributions. First, I document an empirical regularity with important implications for spatial HAC estimation: the relationship between the bandwidth and the magnitude of spatial HAC standard errors follows an inverse-U shape. This means that both too narrow and too wide bandwidths lead to underestimated standard errors and thus to overrejection of true null hypotheses. This finding contradicts the prevailing conventional wisdom in applied work, which holds that choosing wider bandwidths ensures conservative inference.
Second, I propose a simple, non-parametric, data-driven method to select the appropriate bandwidth. The method is based on the empirical covariogram of regression residuals: the estimated bandwidth corresponds to the distance at which residual covariation first crosses zero. Third, I conduct extensive Monte Carlo simulations to evaluate the method. Using spatial correlation structures calibrated to the geography of the contiguous United States, I show that the proposed bandwidth selector controls the false positive rate at or near the nominal 5% level across a wide range of spatial autocorrelation intensities, sample sizes, and spatial configurations. The simulations also provide a systematic comparison of six commonly used kernel functions---Bartlett, Uniform, Epanechnikov, Gaussian, Parzen, and Quartic Biweight---for spatial HAC estimation, finding that the Bartlett and Epanechnikov kernels deliver the best size control.
Related literature. The most closely related work is KimSun2011, who derive the MSE-optimal bandwidth for spatial HAC estimation and propose a parametric plug-in estimator. Their optimal bandwidth is a function of the sample size (\(n^{\alpha}\) for some \(\alpha > 0\)), which contrasts with the approach proposed here. The plug-in estimator of KimSun2011 relies on parametric pilot estimators of unknown quantities, whereas the covariogram-range method proposed in this paper is fully non-parametric and does not require such assumptions. KelejianPrucha2007 provide a detailed treatment of the spatial HAC estimator's properties and derive consistency results; the present paper builds on their theoretical framework.
Recent contributions have taken a different approach by proposing inference methods that bypass the Conley HAC framework altogether. MullerWatson2022 and MullerWatson2023 develop spatial correlation principal components (SCPC) confidence intervals that are robust to arbitrary spatial correlation without requiring bandwidth selection. MullerWatson2024 show that spatial I(1) processes can induce spuriously significant regression results even with spatial HAC or cluster-robust standard errors, and propose spatial differencing methods. ConleyKelly2025 propose spatial basis regressions combined with placebo tests for persistence studies. DellaVignaImbensKimRitzwoller2025 use multiple outcome variables to estimate cross-sectional dependence, sidestepping the bandwidth problem entirely. While each of these methods has merit, they also require departures from the familiar regression framework that applied researchers are accustomed to. The present paper demonstrates that valid inference can be achieved within the standard Conley HAC framework---provided the bandwidth is chosen appropriately. The method is transparent, intuitive, and computationally efficient.
In the time-series HAC literature, KolokotronesStockWalker2024 show that the Bartlett (Newey--West) kernel delivers the highest power among first-order kernels. The present paper extends the comparison of kernel functions to the spatial domain, where such a systematic evaluation has not been conducted.
The rest of the paper is organized as follows. Section (ref) introduces the econometric framework and the spatial HAC estimator. Section (ref) presents the proposed bandwidth selection method with formal assumptions and a consistency result. Section (ref) describes the simulation design for generating spatially correlated data. Section (ref) documents the inverse-U relationship. Section (ref) presents the Monte Carlo results. Section (ref) provides an empirical application. Section (ref) concludes. Section (ref) contains additional details on spatial HAC estimation, kernel functions, and the geostatistical simulation procedure.
Let \(\mathcal{S} \subset \mathbb{R}^2\) be a subset of two-dimensional Euclidean space and let \(\mathbf{s}_i = (s_{1i}, s_{2i})' \in \mathcal{S}\) denote the spatial location of observation \(i\), with \(s_{1i}\) and \(s_{2i}\) representing the two spatial coordinates.\footnote{For unprojected geographic data, these correspond to longitude and latitude on the sphere \(\mathbb{S}^2\). For projected data, they are coordinates in a planar coordinate system. Without loss of generality, this paper refers to locations in \(\mathbb{R}^2\) and implicitly assumes projected data.} Consider the linear regression model
where \(y_i \equiv y(\mathbf{s}_i)\) is the outcome at location \(\mathbf{s}_i\), \(\mathbf{x}_i\) is a \(p \times 1\) vector of regressors, \(\boldsymbol{\beta}\) is the parameter vector of interest, and \(\varepsilon_i \equiv \varepsilon(\mathbf{s}_i)\) is the error term with \(E[\varepsilon_i] = 0\). The error terms fail to satisfy the independence assumption: \(\text{Cov}(\varepsilon_i, \varepsilon_j) \neq 0\) for nearby locations \(\mathbf{s}_i\) and \(\mathbf{s}_j\). Intuitively, this arises because outcomes are not randomly distributed across space; near things tend to be more related than distant things.\footnote{This paraphrases Tobler's Tobler1970 first law of geography.}
Since the law of large numbers applies to weakly dependent spatial data, the OLS estimator \(\hat{\boldsymbol{\beta}}\) remains consistent. The concern when dealing with spatial autocorrelation is therefore inference: the standard heteroskedasticity-robust variance estimator \[ \hat{V}_{\text{HC}} = \left( \sum_{i=1}^n \mathbf{x}_i \mathbf{x}_i' \right)^{-1} \left( \sum_{i=1}^n \hat{\varepsilon}_i^2 \, \mathbf{x}_i \mathbf{x}_i' \right) \left( \sum_{i=1}^n \mathbf{x}_i \mathbf{x}_i' \right)^{-1} \] ignores cross-products \(\hat{\varepsilon}_i \hat{\varepsilon}_j\) for \(i \neq j\) and thus underestimates the true variance of \(\hat{\boldsymbol{\beta}}\) in the presence of positive spatial autocorrelation. This leads to inflated \(t\)-statistics and rejection rates that exceed the nominal level.
Conley1999 proposes a non-parametric estimator for the variance-covariance matrix that accounts for spatial correlation. The spatial HAC estimator takes the form
where
Here, \(d_{ij} = d(\mathbf{s}_i, \mathbf{s}_j)\) is the distance between observations \(i\) and \(j\), \(\varsigma > 0\) is the bandwidth (cutoff distance), and \(K(\cdot)\) is a kernel function satisfying \(K(0) = 1\), \(K(u) = 0\) for \(|u| > 1\), and \(K(u) = K(-u)\). The kernel downweights the contribution of cross-products as a function of distance, assigning zero weight to pairs separated by more than \(\varsigma\). When \(\varsigma = 0\), the estimator reduces to the standard heteroskedasticity-consistent (HC) estimator.
The choice of bandwidth \(\varsigma\) is the central concern of this paper. It determines the distance within which residuals are allowed to be correlated and thus directly affects the magnitude of the estimated standard errors.
Six kernel functions are considered in this study:
\[
\]
where \(u = d_{ij}/\varsigma\) is the normalized distance. The Uniform kernel assigns equal weight to all pairs within the bandwidth, while the others assign declining weights. The Bartlett kernel corresponds to the spatial analogue of the NeweyWest1987 estimator. The Gaussian kernel used here is truncated at the bandwidth, differing from the standard Gaussian density by omitting the normalization constant; this ensures \(K(0) = 1\) and \(K(u) = 0\) for \(u > 1\), consistent with the spatial HAC framework.
In practice, the spatial relations described above are typically accounted for by a connectivity matrix, often referred to as a weighting matrix---a concept dating back to at least Moran1948. Let \(\mathbf{W}\) be the \(n \times n\) matrix describing a relationship between all locations \(\mathbf{s}_1, \ldots, \mathbf{s}_n\), where the \((i,j)\)-th element, \(w_{ij} \equiv [\mathbf{W}]_{ij}\), refers to the distance between units \(i\) and \(j\). In general, this distance does not have to be spatial but can also refer to, e.g., economic or social distance. If \(w_{ij} \neq 0\), units \(i\) and \(j\) are neighbors. In practice, \(\mathbf{W}\) will be a sparse matrix in most empirical applications, meaning that most pairs \(i\) and \(j\) are not classified as neighbors, i.e., \(w_{ij} = 0\).
To illustrate how a spatially lagged variable is constructed, consider the \(i\)-th row of the connectivity matrix and \(x_i\), the value of the scalar variable \(\mathbf{x}\) for unit \(i\). The spatial lag of \(x_i\) can then be written as \[ [\mathbf{W}\mathbf{x}]_i = w_{i1}x_1 + w_{i2}x_2 + \cdots + w_{in}x_n = \sum_{j=1}^n w_{ij} x_j, \] which is a weighted sum of the values of the same variable across all neighbors. In practice, the entries of \(\mathbf{W}\) are usually row-normalized such that \(\sum_{j=1}^n w_{ij} = 1\) for all rows \(i\). Note also that the main diagonal is set to zero, \(w_{ii} = 0\), meaning that a unit is not its own neighbor. Put differently, \(\mathbf{W}\) can be viewed as the spatial lag operator when operationalized via the matrix product \(\mathbf{W}\mathbf{x}\).
To test each regression specification for residual autocorrelation, I rely on Moran's I Moran1950. In its basic form it can be written as \[ I = \frac{n}{\sum_i \sum_j w_{ij}} \cdot \frac{\sum_{i=1}^n \sum_{j=1}^n w_{ij}(\hat{\varepsilon}_i - \bar{\hat{\varepsilon}})(\hat{\varepsilon}_j - \bar{\hat{\varepsilon}})}{\sum_{i=1}^n (\hat{\varepsilon}_i - \bar{\hat{\varepsilon}})^2}, \] where \(n\) is the number of units and \(\sum_i \sum_j w_{ij}\) is the sum of all spatial weights. The denominator \(\sum_{i=1}^n (\hat{\varepsilon}_i - \bar{\hat{\varepsilon}})^2\) is the variance and the numerator involves the spatial lags discussed above. Note the striking resemblance with Pearson's correlation coefficient: Moran's I quantifies the degree of correlation of \(\hat{\varepsilon}_i\) with its spatially lagged neighbors. The expected value under spatial randomness is \(E[I] = -1/(n-1)\) and inference is based on the standardized statistic \(Z = (I - E[I]) / \sqrt{\text{Var}(I)}\). When applied to regression residuals, the computation accounts for the linear projection following Anselin1988 and CliffOrd1981. All Moran's I tests are conducted using the spdep package BivandAltmanAnselinAssuncaoBerkeEtAl2020.
A drawback when working with connectivity matrices is that they usually have to be defined ex-ante by the researcher. For the present study, this is less of a concern because the matrix is used only for testing, not for estimation. The weights matrix is held constant across all simulations---based on a distance band of 200 km---so that reported Moran's I values are comparable across scenarios. In practice, there are many ways to define a connectivity matrix; this study relies on distance band and \(k\)-nearest neighbor specifications.
The key to the proposed bandwidth selection method is the empirical covariogram of regression residuals. For a set of estimated residuals \(\hat{\varepsilon}_i\), \(i = 1, \ldots, n\), the empirical covariogram at lag distance \(h\) is
where \(N(h) = \{(i,j) : h - \delta \leq d(\mathbf{s}_i, \mathbf{s}_j) < h + \delta\}\) is the set of all location pairs separated by a distance falling in the bin centered at \(h\) with half-width \(\delta\), and \(|N(h)|\) is the cardinality of this set. The bin width \(2\delta\) must be chosen large enough to ensure sufficient observations per bin. In practice, this amounts to specifying a number of equally spaced, non-overlapping distance classes.
The covariogram \(\hat{C}(h)\) estimates the covariance between residuals as a function of the distance separating them. In the presence of positive spatial autocorrelation, \(\hat{C}(h)\) is positive for small \(h\) and declines toward zero as \(h\) increases. Beyond a certain distance, residuals are approximately uncorrelated and \(\hat{C}(h)\) fluctuates around zero.
I propose to use the empirical covariogram to estimate the range of spatial dependence in the regression residuals. This estimated range then serves as the bandwidth for computing spatial HAC standard errors.
In practice, \(\hat{C}(h)\) is evaluated at a fixed set of equally spaced, non-overlapping distance bins. Setting \(\eta = 0\), the estimator selects the first bin center at which the empirical covariogram crosses zero. For covariance functions that decay gradually rather than reaching zero exactly (see Remark 4 below), a small positive \(\eta\) can be used; in the Monte Carlo experiments, the default \(\eta = 0\) performs well because the binned covariogram is sufficiently noisy that it crosses zero near the effective range.
This section provides sufficient conditions under which the covariogram range estimator \(\hat{\varsigma}\) is consistent for the true correlation range and the resulting plug-in spatial HAC estimator is consistent for the variance of \(\hat{\boldsymbol{\beta}}\). The conditions are stated at a level that makes transparent what is required; a fully detailed proof in the increasing-domain asymptotic framework is beyond the scope of this paper. The Monte Carlo evidence in Section (ref) provides the primary validation that the method works in finite-sample settings calibrated to empirically relevant spatial configurations.
The mixing and moment conditions are needed because pairs of residual products \(\hat{\varepsilon}_i \hat{\varepsilon}_j\) sharing common locations are dependent; they ensure a law of large numbers for the binned covariogram averages, which does not follow from the growth in the number of pairs \(|N(h)|\) alone.
The separation condition ensures that the true range \(\varsigma_0\) is well-identified: the covariance is bounded away from zero inside the range and exactly zero outside. This rules out pathological cases where \(C(h)\) is “flat'\,' near \(\varsigma_0\), which would make the threshold crossing unstable under small estimation errors. Note that \(C(h)\) need not be positive for all \(h < \varsigma_0\)---it may be negative at some distances (e.g., due to oscillatory behavior or residualization)---but it must be nonzero.
This is consistent with the asymptotic framework of Conley1999 and KelejianPrucha2007, where spatial HAC consistency is established under increasing-domain asymptotics. The condition that \(|N(h)| \to \infty\) across the relevant range of distances ensures that the binned covariogram averages converge.
The kernel continuity requirement ensures that the mapping from bandwidth to the SHAC variance estimate, \(\varsigma \mapsto \hat{V}_{\text{SHAC}}(\varsigma)\), is continuous: adding or removing pairs near the cutoff boundary produces only a negligible change in the estimate, because their kernel weights are near zero.
Sketch of argument. Let \(\{h_b\}_{b=1}^B\) be the fixed bin centers covering \([0, \varsigma_0 + c]\). Under Assumptions (ref) and (ref), the mixing and moment conditions ensure a law of large numbers for the dependent pairwise averages within each bin, yielding \(\hat{C}(h_b) \to_p C(h_b)\) for each \(b\). Since \(B\) is finite, joint convergence holds, and the maximum over bins converges: \[ \max_{b \leq B} |\hat{C}(h_b) - C(h_b)| \to_p 0. \] Because \(\hat{C}(h)\) is defined as a piecewise-constant function on the bins, this is equivalent to uniform convergence over the interval. Given this and the separation condition in Assumption (ref), the estimator \(\hat{\varsigma}\) converges to \(\varsigma_0\): for any \(\delta > 0\), eventually \(|\hat{C}(h_b)|\) is bounded away from zero (and hence above \(\eta\)) for all bins \(h_b \leq \varsigma_0 - \delta\), and \(|\hat{C}(h_b)|\) is close to zero (and hence below \(\eta\)) for all bins \(h_b \geq \varsigma_0 + \delta\). Thus the first bin at which \(|\hat{C}(h_b)| \leq \eta\) must fall in \((\varsigma_0 - \delta, \varsigma_0 + \delta)\), giving \(\hat{\varsigma} \to_p \varsigma_0\). The second claim follows from the kernel continuity in Assumption (ref)---which ensures \(\hat{V}_{\text{SHAC}}(\hat{\varsigma}) - \hat{V}_{\text{SHAC}}(\varsigma_0) \to_p 0\). Specifically, under increasing-domain sampling, the number of observation pairs in a thin boundary shell \((\varsigma_0 - \epsilon, \varsigma_0 + \epsilon)\) grows at most proportionally to \(\epsilon\), so with a continuous compact-support kernel the marginal contribution of these pairs to \(\hat{\mathbf{\Omega}}\) vanishes as \(\epsilon \downarrow 0\). Combined with the consistency of \(\hat{V}_{\text{SHAC}}(\varsigma_0)\) established by KelejianPrucha2007, this yields the result.
The proposed method is computationally efficient. The empirical covariogram is computed using the variogram() function from gstat Pebesma2004 with the covariogram\ =\ TRUE option. The distance bins are determined by two parameters: width (the bin width \(2\delta\)) and cutoff (the maximum distance). As a default, the bin width is set to ensure approximately 100--200 distance bins up to two-thirds of the maximum inter-point distance. The range is then extracted as the first distance at which the empirical covariogram crosses zero. For a typical dataset with several thousand observations, the entire procedure takes a fraction of a second.
Figure (ref) illustrates the method on four simulated datasets with increasing degrees of spatial autocorrelation.
The covariogram-range bandwidth selector has a useful conceptual connection to the fixed-\(b\) asymptotic framework. In the classical time-series setting, two asymptotic regimes govern the behavior of HAC standard errors. Under the traditional small-\(b\) asymptotics of Andrews1991 and NeweyWest1987, the bandwidth \(M\) grows with the sample size \(T\) but satisfies \(M/T \to 0\), so that the HAC estimator is consistent and the \(t\)-statistic converges to a standard normal. Under the fixed-\(b\) asymptotics of KieferVogelsangBunzel2000 and SunPhillipsJin2008, the bandwidth satisfies \(M/T \to b \in (0, 1]\), so the HAC estimator is inconsistent but the resulting test statistic converges to a non-standard distribution that provides more accurate finite-sample size control. Sun2014 argues that fixed-\(b\) inference is generally preferable, as it accounts for the estimation uncertainty in the HAC estimator itself. BesterConleyHansenVogelsang2016 extend fixed-\(b\) theory to the spatial setting: when the Conley bandwidth grows proportionally with the diameter of the spatial domain---i.e., \(\varsigma / \text{diam}(S_n) \to b > 0\)---the SHAC estimator converges to a nondegenerate random matrix rather than the true variance, and the \(t\)-statistic has a non-standard pivotal limit distribution whose critical values must be obtained by simulation.
While the covariogram-range selector does not implement spatial fixed-\(b\) inference, it shares the key practical intuition that motivates the fixed-\(b\) literature: smoothing should be anchored to the scale of dependence rather than shrinking relative to the sample. The estimator \(\hat{\varsigma}\) targets a fixed physical quantity---the correlation range \(\varsigma_0\)---determined entirely by the spatial dependence structure of the data-generating process. Unlike MSE-optimal rules Andrews1991, where the bandwidth grows with \(T\) but satisfies \(M/T \to 0\) and tends to undersmooth in finite samples, the covariogram-range bandwidth does not depend on \(n\) and is instead pinned to the dependence scale. This fixed-span behavior---anchoring to a physical quantity rather than to sample size---reduces the usual smoothing trade-off that Sun2014 identifies as the primary source of size distortions in HAC inference.
Formally, however, the covariogram-range method operates in a different asymptotic regime. Under the increasing-domain asymptotics used throughout this paper, a fixed physical bandwidth implies \(\varsigma_0 / \text{diam}(S_n) \to 0\), placing the covariogram-range selector in the small-\(b\) regime of BesterConleyHansenVogelsang2016, not the fixed-\(b\) regime. Moreover, under compact support (Assumption (ref)), the covariance is literally zero beyond \(\varsigma_0\), so there is no “missing tail” to truncate---unlike the spatial fixed-\(b\) setting of BesterConleyHansenVogelsang2016, where the bandwidth truncates a non-negligible portion of the dependence structure and the SHAC estimator converges to a nondegenerate random limit. At the true cutoff \(\varsigma_0\), the SHAC estimator captures the entire covariance structure without truncation bias and is therefore consistent. The plug-in bandwidth \(\hat{\varsigma}\) converges to \(\varsigma_0\), standard SHAC consistency obtains, and no non-standard critical values or simulated reference distributions are required. The finite-sample size improvements documented in Section (ref) arise not from the fixed-\(b\) mechanism of embracing estimator randomness, but from targeting the correct dependence scale directly---avoiding the undersmoothing that plagues MSE-optimal bandwidth selectors. Nevertheless, in any finite spatial sample the bandwidth-to-domain ratio is a non-negligible constant, so the resulting smoothing can resemble the larger-smoothing behavior that motivates fixed-\(b\) corrections.
To evaluate the proposed bandwidth selector, I require spatially correlated random variables with precisely controlled correlation ranges. I generate these using unconditional geostatistical simulation (simple kriging) with the gstat package Pebesma2004; see Cressie1993 for a comprehensive treatment of the geostatistical framework. The full details of the simulation procedure---including the kriging system, semivariogram specifications, and implementation---are provided in Section (ref).
Four spatial samples are used throughout the study, all covering the bounding box of the contiguous United States:
The three lattice samples provide controlled variation in sample density, while the county centroid sample tests performance under highly irregular spacing.
In each simulation iteration:
This process is repeated for noise range parameters starting from 0 (no spatial correlation) in 15 sequential increments until extreme degrees of spatial correlation are reached, with 5,000 iterations per configuration. The noise range parameter governs the input to the semivariogram; the realized correlation range varies stochastically across draws. Computations were carried out on the Akropolis cluster at the University of Chicago using Rmpi for parallel processing RCoreTeam2025. To ensure exact replicability, the pseudo-random number seed was set using clusterSetRNGStream() with Pierre L'Ecuyer's RngStreams LEcuyer1999 (seed: 1908).
For each simulation, the following standard errors are computed:
Figure (ref) displays the relationship between the kernel bandwidth and the magnitude of Conley standard errors for a single realization of a spatially correlated random field across the four spatial samples. The pattern is striking and consistent: the standard error first increases as the bandwidth grows from zero (at which point the Conley SE equals the HC1 SE), reaches a maximum near the true correlation range, and then decreases steadily. The method proposed in this paper always picks the peak of this inverse-U curve as the standard error.
This inverse-U shape has immediate practical consequences. Starting from a narrow bandwidth, expanding the range over which residuals are allowed to be correlated initially leads to larger standard errors---as expected, since positive residual correlation adds to the estimated variance. However, beyond the range of actual correlation, additional observations included in the kernel contribute only noise to the variance estimate: their cross-products \(\hat{\varepsilon}_i \hat{\varepsilon}_j\) have expectation zero but dilute the positive contributions from truly correlated pairs. As the bandwidth continues to grow, this dilution effect dominates and the standard error declines.
For very wide bandwidths, the standard error can fall below the HC1 level. This directly contradicts the conventional wisdom that wider bandwidths ensure conservative inference. In fact, choosing a bandwidth much larger than the true correlation range can lead to standard errors that are smaller than even the heteroskedasticity-only estimator, resulting in overrejection that exceeds even the overrejection from ignoring spatial correlation entirely.
The practical implication is clear: the bandwidth should be matched to the actual range of spatial correlation in the residuals. This is precisely what the covariogram range method proposed in Section (ref) achieves.
Figure (ref) visualizes the covariogram range method on three regression residual series with increasing degrees of spatial autocorrelation. In each panel, the empirical covariogram declines from a positive value at short distances and crosses zero at the estimated range \(\hat{\varsigma}\) (marked by the red vertical line). Beyond \(\hat{\varsigma}\), the covariogram fluctuates around zero, indicating the absence of residual correlation. As the degree of spatial autocorrelation increases from panel (a) to (c), the estimated range shifts rightward, correctly tracking the growing extent of spatial dependence.
Table (ref) presents the main results across the three regular lattice samples. The covariogram-range method with the Epanechnikov kernel consistently maintains rejection rates near 5% across all sample sizes and spatial configurations, while HC1 rejection rates rise sharply with spatial autocorrelation. Larger samples (\(n = 10{,}320\)) tend to produce slightly more accurate size control, consistent with the consistency result in Proposition (ref).
Table (ref) presents the detailed results for the county centroid sample (\(n = 3{,}108\)). Several findings stand out.
First, HC1 standard errors lead to substantial overrejection as spatial autocorrelation increases. At the highest noise ranges, the false positive rate exceeds 60%, far above the nominal 5%.
Second, a fixed narrow bandwidth of 25 km offers no improvement over HC1 when the true correlation range exceeds this bandwidth. Conversely, a fixed wide bandwidth of 2,500 km---meant to be “conservative”---also fails to control size, confirming the inverse-U prediction: at very wide bandwidths, the standard error is deflated below the HC1 level.
Third, the proposed covariogram-range method with Bartlett, Uniform, or Epanechnikov kernels controls the rejection rate remarkably well. Across all noise range configurations, the false positive rate remains close to the nominal 5%. At the most extreme degrees of spatial autocorrelation (with noise ranges far above empirically realistic levels), the rejection rate rises to approximately 8--10%, reflecting the inherent difficulty of spatial HAC estimation when the correlation range approaches a substantial fraction of the spatial domain.
Taken together, the results demonstrate that the covariogram-range bandwidth selector provides effective size control across the empirically relevant spectrum of spatial dependence. For the correlation ranges most commonly encountered in applied work---where HC1 rejection rates already exceed 15--30%---the proposed method maintains false positive rates in the range of 5--8%, close to the nominal level. The modest overrejection observed at extreme correlation ranges reflects a regime that pushes far beyond what typical georeferenced datasets exhibit, and even there the covariogram-range method substantially outperforms all fixed-bandwidth alternatives.
Table (ref) reports rejection frequencies for all six kernel functions. The Bartlett and Epanechnikov kernels consistently deliver rejection rates closest to the nominal 5% level. The Uniform kernel performs well at moderate spatial autocorrelation but tends to slightly overreject at high levels. The Gaussian, Parzen, and Quartic Biweight kernels show somewhat more variable performance, with a tendency toward overrejection at higher noise ranges. Based on these results, I recommend the Bartlett or Epanechnikov kernel for applied work, consistent with the time-series finding of KolokotronesStockWalker2024 regarding the optimality of Bartlett among first-order kernels.
Figure (ref) provides a visual summary of the results presented in Table (ref) and Table (ref), plotting every individual simulation run and illustrating how false rejection rates increase with the degree of spatial correlation. Each column corresponds to one of the four spatial samples; each row to a different standard error type. The top row (HC1) shows that \(t\)-statistics grow with the correlation range when both the outcome and regressor are spatially correlated, producing the substantial overrejection documented in the tables. The second row demonstrates the covariogram-range method: the \(t\)-statistics remain centered below the critical value regardless of the correlation range, with approximately 5% of dots exceeding the threshold. The third and fourth rows show the failure modes of fixed bandwidths: a narrow bandwidth (25 km) behaves like HC1 when the true range exceeds 25 km, and a wide bandwidth (2,500 km) deflates standard errors and produces high rejection rates across the board.
To illustrate the practical relevance of the proposed method, I apply it to a regression using U.S. county-level data. The outcome variable is the share of the Black population in each county, and the regressor is the share of the Hispanic population---two variables that exhibit well-known spatial clustering patterns across the United States. The data is from the latest US census, accessed through IPUMS IPUMS2025.
Table (ref) presents the results. The HC1 standard error is the smallest, producing the largest \(t\)-statistic. The Conley standard errors vary substantially with the bandwidth choice. At the estimated covariogram range of approximately 980 km, the standard error is the largest among all choices shown, consistent with the inverse-U pattern. This is the appropriate bandwidth for this regression: the estimated range reflects the spatial extent of residual correlation in these demographic variables. At bandwidths of 500 km and 2,500 km---which might seem “safer” to an applied researcher---the standard errors are actually smaller, illustrating precisely the misleading effect documented in Section (ref).
This example demonstrates a scenario commonly encountered in applied spatial research: widely used variables with strong spatial clustering patterns that extend over large distances. The covariogram range method correctly identifies this extent and produces appropriately sized standard errors.
This paper addresses the unresolved problem of bandwidth selection for spatial HAC standard errors. I make three main contributions.
First, I document that the relationship between the kernel bandwidth and the magnitude of spatial HAC standard errors follows an inverse-U shape. This implies that both too narrow and too wide bandwidths lead to underestimated standard errors. The finding contradicts the common advice that wider bandwidths are more conservative and establishes that the bandwidth choice is not merely a matter of being “generous” with the cutoff distance.
Second, I propose a simple, non-parametric bandwidth selector based on the empirical covariogram of regression residuals. The estimator identifies the distance at which residual covariation first crosses zero---the correlation range---and uses this as the bandwidth for the spatial HAC estimator. Under standard regularity conditions, the estimator is consistent for the true correlation range.
Third, I show through extensive Monte Carlo simulations that the proposed method controls the false positive rate at or near the nominal 5% level across a wide range of spatial autocorrelation intensities and sample configurations. A systematic comparison of six kernel functions finds that the Bartlett and Epanechnikov kernels deliver the best size control, extending the time-series finding of KolokotronesStockWalker2024 to the spatial domain. At extreme levels of spatial correlation, where the correlation range approaches a substantial fraction of the spatial domain, the rejection rate rises to approximately 8--10%. This residual overrejection should be understood as a fundamental limitation of spatial inference under very strong dependence rather than a failure of the proposed selector: when effective degrees of freedom become small, no HAC-type estimator---whether based on fixed bandwidths, cluster-robust methods, or data-driven selection---performs well. The key result is comparative: across all configurations considered, the covariogram-range method substantially improves size control relative to HC1, fixed narrow bandwidths, and fixed wide bandwidths.
The proposed method is implemented in the R package SpatialInference.
During the preparation of this work the author used Claude (Anthropic) in order to assist with manuscript drafting. After using this tool, the author reviewed and edited the content as needed and takes full responsibility for the content of the published article.