--- title: "Conley standard errors with fastconley" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Conley standard errors with fastconley} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>") has_fixest <- requireNamespace("fixest", quietly = TRUE) read_csv <- function(file) { paths <- c( system.file("benchmarks", file, package = "fastconley"), file.path("..", "inst", "benchmarks", file), file.path("inst", "benchmarks", file) ) paths <- paths[nzchar(paths) & file.exists(paths)] if (!length(paths)) stop("Cannot find benchmark file: ", file) read.csv(paths[1], stringsAsFactors = FALSE) } bench <- read_csv("fastconley-benchmark-results.csv") large <- read_csv("fastconley-large-data-results.csv") fmt_sec <- function(x) sprintf("%.3f", x) fmt_x <- function(x) sprintf("%.1fx", x) fmt_e <- function(x) formatC(x, format = "e", digits = 2) fmt_int <- function(x) format(x, big.mark = ",", scientific = FALSE, trim = TRUE) fmt_pairs <- function(x) formatC(x, format = "e", digits = 2) session_path <- system.file("benchmarks", "fastconley-large-data-session.txt", package = "fastconley") if (session_path == "") { session_path <- file.path("..", "inst", "benchmarks", "fastconley-large-data-session.txt") } session_lines <- if (file.exists(session_path)) readLines(session_path, warn = FALSE) else "" line_one <- function(pattern) { out <- grep(pattern, session_lines, value = TRUE) if (length(out)) out[1] else NA_character_ } strip_value <- function(pattern, prefix) sub(prefix, "", line_one(pattern)) cpu_model <- strip_value("Model name:", "Model name:\\s+") cpu_count <- strip_value("^CPU\\(s\\):", "CPU\\(s\\):\\s+") core_count <- strip_value("Core\\(s\\) per socket:", "Core\\(s\\) per socket:\\s+") thread_count <- strip_value("Thread\\(s\\) per core:", "Thread\\(s\\) per core:\\s+") mem_line <- line_one("^Mem:") r_line <- line_one("^R version") os_line <- strip_value("^Running under:", "Running under:\\s+") run_date <- sub("\\..*$", "", unique(large$run_date)[1]) ``` When regression errors are spatially correlated — nearby observations hit by the same local shocks — conventional standard errors understate the uncertainty of the estimates, often badly. Conley (1999) standard errors fix this the same way Newey-West fixes serial correlation: covariance between observations is admitted up to a chosen distance cutoff and zeroed beyond it. If your regressors have spatial structure (almost everything geographic does: climate, infrastructure, institutions, treatments rolled out by region), this is the difference between a significant and an insignificant result far more often than is comfortable. `fastconley` computes Conley spatial HAC standard errors for models estimated with `fixest` (`feols()` including IV, `feglm()`, `fepois()`) or `lfe::felm()` (OLS and IV/2SLS) — for cross-sections and for panels, where it implements the full Conley (1999) spatial-*temporal* HAC. For GLM fits the variance is the M-estimation sandwich built from the scores and inverse Hessian `fixest` stores on the fit, so Poisson panels get the spatial + serial correction too. Four properties distinguish it: - **Exact.** Distances use the selected exact geometry (haversine or spherical arc distance, or straight-line chord distance), and every within-cutoff pair enters with its exact kernel weight — no dense-matrix shortcuts, no approximate distance formulas. Against a brute-force reference, results agree to near machine precision (the dense-baseline benchmark below shows this directly). - **Complete for panels.** With `lag_cutoff > 0` the spatial adjustment within each period is combined with a Newey-West serial HAC term within each unit — the spatial-temporal correction of Conley (1999) as popularized by Hsiang (2010). `fixest`'s `vcov_conley()` is spatial-only; it makes you choose between the Conley and the Driscoll-Kraay/Newey-West corrections, while here they enter one covariance matrix together. - **Reproducible.** Results are bit-identical across thread counts: the same data and settings give the same covariance matrix on a laptop with `ncores = 2` and a server with `ncores = 64`. - **Fast at any size.** Small problems are instant; a 1,000,000-observation cross-section takes about a second, and regular rasters with hundreds of billions of within-cutoff pairs take seconds rather than hours (the benchmark sections quantify this). The usual workflow is: 1. Estimate the model. 2. Pick a spatial cutoff that is meaningful for the application. 3. Pass the fitted model and coordinate columns to `vcovSpHAC()` (when the columns have standard names like `lat` / `lon`, they are auto-detected and can be omitted). 4. Use the returned covariance matrix in `summary()`, `etable()`, or your preferred table-making code. The examples below use `fixest`, because many users already rely on its formula syntax and reporting tools. The same `vcovSpHAC()` call also works with `felm` models when the fit was created with `keepCX = TRUE`. ## 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. ```{r motivation-sim, eval=has_fixest} 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 ``` ```{r motivation-fit, eval=has_fixest} 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 `r if (has_fixest) sprintf("about %.1f", se_tab[3, 2] / se_tab[1, 2]) else "several"` — 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. ```{r quick-start-fixest, eval=FALSE} 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. ```{r quick-start-fixest-estimation-vcov, eval=FALSE} 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. ```{r quick-start-fixest-matrix, eval=FALSE} 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. ```{r quick-start-lfe, eval=FALSE} 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`. ```{r quick-start-panel, eval=FALSE} 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. | Argument | Practical meaning | |---|---| | `dist_cutoff` | Spatial range, in kilometers. Observations farther apart do not contribute to each other's Conley adjustment. | | `kernel = "uniform"` | All observations inside the cutoff receive equal weight. | | `kernel = "bartlett"` | Observations inside the cutoff are downweighted with distance. | | `lag_cutoff` | Number of time lags for serial HAC in panel data. Use `0` for spatial-only inference. | | `ncores` | Number of threads used for the covariance calculation. The default uses `getOption("fastconley.ncores")` when set, otherwise all detected logical cores; `_R_CHECK_LIMIT_CORES_=TRUE` caps it at two. | | `pixel = 0` | Exact deduplication of identical coordinates. Positive values trade a little spatial precision for more aggregation. | | `ssc = TRUE` | Small-sample correction. The default scales the matrix by $n/(n-K)$, counting absorbed fixed-effect levels in $K$ — exactly `fixest`'s default Conley correction. $K$ is the fitting package's own parameter count; with three or more absorbed fixed effects `lfe` and `fixest` can count the estimable levels differently, so their default-`ssc` results differ by that factor while `ssc = FALSE` results are identical. `FALSE` applies none, reproducing `rbluhm/conley` and fastconley before 0.9.0. | | `psd_fix = TRUE` | Conley kernels do not formally guarantee a positive semi-definite matrix. The default clamps negative eigenvalues (same semantics as `fixest`'s `vcov_fix`); either way a warning reports a noticeably non-PSD result. | 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`. ### Cutoff Sensitivity Is Cheap — Report It The honest answer to "what cutoff should I use?" is usually a range, and referees increasingly ask for it. Because the covariance call is fast, looping over cutoffs costs almost nothing — there is no reason to report a single number when you can report the curve: ```{r cutoff-sensitivity, eval=has_fixest} 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) ``` In this simulation the errors were built with a ~150 km correlation range, and the standard error climbs steeply until the cutoff covers that range, then levels off — exactly the pattern to look for in real data. If the standard errors keep growing with the cutoff, the spatial correlation is longer-ranged than the design assumed (or there is a spatial trend the model should absorb). ## 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. | Data layout | Recommended `vcovSpHAC()` settings | What it means | |---|---|---| | Scattered cross-section | Omit `unit` and `time`; use `lag_cutoff = 0`; leave `method = "auto"` or set `method = "pairwise"`. | Exact spatial Conley for point locations. Under `auto`, ordinary scattered data use the pairwise engine. | | Repeated cross-sections | Supply `time`, omit `unit`, and keep `lag_cutoff = 0`. | Spatial pairs are limited to each period; every row is its own unit. | | Repeated locations | Same as scattered data; start with `pixel = 0`. | Rows with identical coordinates are combined exactly before the covariance is computed. | | Approximate location snapping | Set a positive `pixel`, such as `pixel = 5`, only after checking sensitivity. | Nearby locations are snapped to a kilometer grid. This can be faster, but it is an approximation. | | Regular raster or grid | Use `method = "auto"`, or `method = "grid"` if you want to require a regular latitude-longitude lattice. | Exact Conley for complete or partly occupied regular rasters. If `method = "grid"` and the data are not a lattice, the call errors instead of falling back. | | Balanced panel with fixed locations | Provide `unit` and `time`; set `lag_cutoff` as needed; set `balanced_pnl = TRUE`. | Spatial correlation is computed within each period; serial HAC is computed within unit. The balanced setting is faster when every period has the same units and locations do not move. | | Unbalanced panel or moving locations | Provide `unit` and `time`; set `lag_cutoff` as needed; leave `balanced_pnl = FALSE`. | Uses the general panel route. This is safer when units enter, exit, duplicate within a period, or change coordinates. | A few concrete templates are often enough. For ordinary point data: ```{r template-point-data, eval=FALSE} 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: ```{r template-raster-data, eval=FALSE} 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: ```{r template-panel-data, eval=FALSE} 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. ### Machine ```{r machine-table, echo=FALSE} machine_tab <- data.frame( item = c("run date", "CPU", "logical CPUs", "cores per socket", "threads per core", "memory", "operating system", "R", "fastconley", "fixest", "lfe"), value = c(run_date, cpu_model, cpu_count, core_count, thread_count, mem_line, os_line, r_line, unique(large$fastconley_version)[1], unique(large$fixest_version)[1], unique(large$lfe_version)[1]), check.names = FALSE ) knitr::kable(machine_tab) ``` ### One Million Observations This cross-section has 1,000,000 observations, 10 regressors, a 100 km cutoff, the uniform kernel, and spherical distance. A dense spatial weight matrix would be about 7.3 TiB. ```{r large-xsection-table, echo=FALSE} lx <- large[large$section == "large_xsection", ] lx <- lx[order(lx$threads), ] base <- lx$seconds[lx$threads == min(lx$threads)][1] lx_tab <- data.frame( observations = fmt_int(lx$n_obs), regressors = lx$k, cutoff_km = lx$cutoff_km, threads = lx$threads, seconds = fmt_sec(lx$seconds), speedup_vs_2_threads = fmt_x(base / lx$seconds), dense_matrix_size = lx$implied_dense_weight, check.names = FALSE ) knitr::kable(lx_tab) ``` On this machine, the 16-thread run computes the covariance matrix in about `r fmt_sec(lx$seconds[lx$threads == 16][1])` seconds. ### Comparison With fixest `fixest` has its own Conley covariance implementation through `vcov_conley()`. The comparison below uses smaller cross-sections where both methods are comfortable to run locally. Both models are fitted once outside the timer; only the post-estimation covariance call is timed. ```{r fixest-comparison-table, echo=FALSE} xs <- bench[bench$section == "cross_section", ] xs_fc <- xs[xs$method == "fastconley", ] xs_fx <- xs[xs$method == "fixest", ] xs_w <- merge(xs_fc, xs_fx, by = c("benchmark", "n_obs", "k", "cutoff_km", "threads"), suffixes = c("_fastconley", "_fixest")) xs_w <- xs_w[order(xs_w$n_obs, xs_w$cutoff_km, xs_w$threads), ] xs_tab <- data.frame( observations = fmt_int(xs_w$n_obs), cutoff_km = xs_w$cutoff_km, threads = xs_w$threads, fastconley_seconds = fmt_sec(xs_w$seconds_fastconley), fixest_seconds = fmt_sec(xs_w$seconds_fixest), speedup = fmt_x(xs_w$seconds_fixest / xs_w$seconds_fastconley), max_abs_difference = fmt_e(xs_w$max_abs_diff_fastconley), check.names = FALSE ) knitr::kable(xs_tab) ``` The last column needs an explanation, because the two packages do not return the same matrix: on these configurations they disagree by roughly one percent in relative terms. The disagreement is **not** a fastconley approximation. A brute-force reference — forming the full great-circle distance matrix and the exact double sum, feasible at a few thousand observations — matches `fastconley` to near machine precision (about 1e-12; see the dense-baseline table below, which performs exactly this check). The gap comes from `fixest` 0.14.x's internal `distance = "spherical"` formula, which flips a small share of pairs near the cutoff boundary relative to the exact great-circle distance. The workload being timed is the same; `fastconley` is the numerically exact one. #### Matching fixest's numbers If you are migrating from `vcov_conley()` and want to reconcile results, two things matter: 1. **Small-sample corrections.** Since v0.9.0 the defaults already agree: both packages scale by $n/(n-K)$, where $K$ counts all estimated parameters including absorbed fixed-effect levels (`fixest`'s cluster adjustment is a no-op for Conley vcovs). To compare *uncorrected* matrices instead, pass `ssc = FALSE` to `vcovSpHAC()` and `ssc(adj = FALSE, cluster.adj = FALSE)` to `fixest`. 2. **Distances.** After aligning the corrections, the remaining ~1e-2 relative difference is fixest's approximate spherical distance, as described above. It shrinks with the cutoff's distance from the mass of pairs and is not something a setting can remove. On panels the two packages compute structurally different objects (see the Balanced Panel section below), so coefficients' standard errors should not be expected to match even after these alignments. ### Balanced Panel The panel scenario has 50,000 units observed for 10 periods, for 500,000 rows and 10 regressors. The covariance includes within-period spatial Conley adjustment and one serial HAC lag within unit. ```{r panel-table, echo=FALSE} pn <- large[large$section == "large_panel", ] pn_tab <- data.frame( units = fmt_int(pn$n_unit), periods = pn$n_time, observations = fmt_int(pn$n_obs), regressors = pn$k, cutoff_km = pn$cutoff_km, lag_cutoff = 1, threads = pn$threads, seconds = fmt_sec(pn$seconds), dense_one_period_size = pn$implied_dense_weight, check.names = FALSE ) knitr::kable(pn_tab) ``` This is not the same calculation as applying a cross-sectional Conley covariance to all stacked panel rows. Spatial correlation is restricted to observations in the same period; serial correlation is handled separately within unit. ### Regular Raster Regular rasters are common in climate, remote-sensing, and gridded economic data. The benchmark below uses a 1500 by 1500 latitude-longitude lattice, or 2.25 million cells, with a 250 km cutoff. ```{r raster-table, echo=FALSE} rg <- large[large$section == "large_grid", ] rg <- rg[order(rg$kernel), ] rg_tab <- data.frame( kernel = rg$kernel, cells = fmt_int(rg$n_obs), cutoff_km = rg$cutoff_km, threads = rg$threads, seconds = fmt_sec(rg$seconds), dense_matrix_size = rg$implied_dense_weight, implied_within_cutoff_pairs = fmt_pairs(rg$implied_pairs), check.names = FALSE ) knitr::kable(rg_tab) ``` The important point for users is that this remains an exact Conley calculation for the raster. It is just using the regular grid structure instead of treating the data as millions of unrelated points. For scale: treating the same lattice as scattered points and running the pairwise engine on all ~1.8e11 within-cutoff pairs takes minutes, not seconds — measured during v0.8.0 development on this machine at 16 threads, about 244 s for the uniform kernel and 710 s for bartlett, against the sub-five-second grid-engine times in the table. That two-orders-of-magnitude gap is why `method = "auto"` exists: gridded data should never pay the point-pair price. ### Repeated Locations Repeated coordinates appear in panels, administrative data, firm or household microdata linked to towns, and gridded data with multiple observations per cell. At `pixel = 0`, `fastconley` combines only exactly identical coordinates, so the result is unchanged apart from floating-point summation order. ```{r pixel-table, echo=FALSE} pix <- bench[bench$section == "pixel", ] pix <- pix[order(pix$cases, decreasing = TRUE), ] pixel0 <- pix$seconds[pix$method == "pixel=0"][1] pix_tab <- data.frame( setting = pix$method, raw_rows = fmt_int(pix$raw_rows), coordinate_cases = fmt_int(pix$cases), seconds = fmt_sec(pix$seconds), speed_vs_pixel0 = fmt_x(pixel0 / pix$seconds), max_abs_difference_vs_pixel0 = fmt_e(pix$max_abs_diff), check.names = FALSE ) knitr::kable(pix_tab) ``` Positive `pixel` values are optional. They can be useful for very large data sets with dense nearby coordinates, but they should be treated as an approximation and checked against `pixel = 0` on a smaller sample. Rounding the two coordinate axes can move one point by roughly `pixel / sqrt(2)` and can change a two-point distance by as much as roughly `sqrt(2) * pixel`. ### Dense Baseline For intuition, the small benchmark below actually forms the dense comparison object. This is feasible at a few thousand rows and already shows why the dense approach does not scale. ```{r dense-table, echo=FALSE} dense <- bench[bench$section == "dense", ] dense_fast <- dense[dense$method == "fastconley", ] dense_base <- dense[dense$method == "dense", ] dense_fast <- dense_fast[order(dense_fast$n_obs), ] dense_base <- dense_base[order(dense_base$n_obs), ] dense_tab <- data.frame( observations = fmt_int(dense_fast$n_obs), regressors = dense_fast$k, cutoff_km = dense_fast$cutoff_km, dense_seconds = fmt_sec(dense_base$seconds), fastconley_seconds = fmt_sec(dense_fast$seconds), speedup = fmt_x(dense_base$seconds / dense_fast$seconds), max_abs_difference = fmt_e(dense_fast$max_abs_diff), check.names = FALSE ) knitr::kable(dense_tab) ``` ## Reproducing The Benchmarks The quick local benchmark suite is: ```bash FASTCONLEY_BENCH_THREADS=8 FASTCONLEY_BENCH_REPS=2 \ Rscript inst/benchmarks/run-fastconley-benchmarks.R ``` The large-data benchmark suite is: ```bash 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.)