From 6c4ab949cd79f55de7cb8a96813c2fdcdb38a087 Mon Sep 17 00:00:00 2001
From: igerber <isaac.gerber@gmail.com>
Date: Tue, 26 May 2026 13:39:21 -0400
Subject: [PATCH] ConleySpatialHAC methodology-review-tracker promotion +
 Bertanha-Imbens 2014 citation correction
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Three bundled changes:

1. ConleySpatialHAC methodology-review-tracker promotion (In Progress → Complete).
   - METHODOLOGY_REVIEW.md L91 + L1262-1290 flipped Complete with full
     Verified Components / Test Coverage / R Comparison Results inline
     table / Corrections Made / Deviations / Outstanding Concerns
     structure, with Last Review = 2026-05-26.
   - New tests/test_methodology_conley.py (~1600 LoC; 10 classes;
     60 tests, 5 of them @pytest.mark.slow). Paper-equation-numbered
     Verified Components walk-through covering:
       * Eq. 4.2 cross-sectional sandwich (pairwise-distance
         specialization; Eq. 3.13 is the lattice-indexed form reserved
         for grid coordinates)
       * Eq. 4.2 HC0 + rank-1 limits + K(0) = 1 diagonal contribution
       * Andrews (1991) HAC lag truncation matching
         conleyreg::time_dist.cpp
       * Haversine convention with Earth radius 6371.01 km (matches
         conleyreg::haversine_dist)
       * Phase 2 panel block-decomposed sandwich at atol=1e-12 (internal
         cross-check vs hand-coded reference + cluster time-invariance
         contract)
       * Wave A #120 sparse k-d-tree numerical correctness (sparse-vs-
         dense bit-identity at atol=1e-10 on cross-sectional, panel,
         and sparse-with-cluster paths)
       * R conleyreg v0.1.9 parity at atol=1e-6 on 6 fixtures
         (3 cross-sectional + 3 panel) + sparse-forced cross-sectional
         + sparse-forced panel + time-asymmetric kernel literal-matching
     Three dedicated deviations-area classes:
       * TestConleyLibraryExtensions (6 tests, no R correspondence):
         combined spatial+cluster product kernel Wave A #119 (two
         limit-fixture anchors), callable conley_metric validation
         Wave A #123, sparse k-d-tree activation Wave A #120,
         indefiniteness guard on Bartlett + uniform kernels.
       * TestConleyDeviationsFromR (3 tests): 1-D radial Bartlett vs
         paper's 2-D separable Eq. 3.14, time-label normalization,
         independent temporal kernel deferred.
       * TestConleyDeferrals (5 tests): fail-closed NotImplementedError
         / TypeError contracts for LinearRegression + survey_design,
         DiD/MPD/TWFE + survey_design (estimator-level via
         conley.py::_validate_conley_estimator_inputs), Conley +
         weights (rejected for any weight_type — pweight, aweight,
         fweight), SyntheticDiD + Conley (TypeError), wild_bootstrap
         + Conley.
   - tests/test_conley_vcov.py extracted 1135+ LoC out (4248 → 3113);
     defensive surface preserved (input validation, NaN/inf guards,
     dispatch-level validity, estimator-level integration smoke tests,
     set_params atomicity, sparse-path activation thresholds + density-
     gate fallback). Module docstring refreshed to describe its current
     defensive-regression role.

2. Stale priority-queue cleanup at METHODOLOGY_REVIEW.md L1386:
   - PreTrendsPower removed (already Complete since 2026-05-19).
   - ConleySpatialHAC removed (this PR).
   - Substantive-review-blocked renumbered #2-#5 → #1-#4.
   - Consolidation-pass-blocked renumbered #6-#8 → #5-#6.

3. Bertanha-Imbens 2014 citation correction across 16 sites:
   - linalg.py × 8, conley.py × 1, llms-full.txt × 2, REGISTRY.md × 4,
     spillover.rst × 1.
   - NBER w20773 is on FRD external validity, NOT weighted spatial-HAC.
   - The boundary is now framed as a tri-part contract:
     * Shipped — SpilloverDiD + Conley + survey via Wave E.1/E.2/E.3
       (PR #468/#474/#482, stratified-Conley sandwich on PSU totals
       with within-PSU serial Bartlett HAC for lag_cutoff > 0);
       TwoStageDiD + Conley + survey via Wave E.3 parity (PR #485).
     * Deferred (generic linalg surface, any weight_type) — DiD/MPD/
       TWFE/LinearRegression generic path + Conley + survey_design;
       LinearRegression / compute_robust_vcov Conley + weights rejected
       for pweight, aweight, AND fweight (weighted Conley is not
       implemented on the generic linalg surface).
     * Open methodological question (subset) — the pweight /
       survey_design portion of the deferral additionally lacks a
       canonical methodological extension of Conley (1999) for weighted
       spatial-HAC under probability sampling.
   - REGISTRY sites use canonical `**Note (open methodological
     question):**` label wrapper per CLAUDE.md "Documenting Deviations".
   - Historical CHANGELOG entries (pre-[Unreleased]) intentionally
     retain Bertanha-Imbens 2014 attribution as accurate records of
     past release claims.

Verification:
- pytest tests/test_methodology_conley.py: 60 tests (55 unit + 5 slow)
  pass.
- pytest tests/test_methodology_conley.py + tests/test_conley_vcov.py:
  176 pass + 5 slow deselected (preserves pre-edit baseline coverage
  176 vs original 175 + 1 new uniform-kernel indefiniteness guard).
- grep -rn "Bertanha" diff_diff/ docs/ benchmarks/: 0 hits.
- black + ruff clean.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 CHANGELOG.md                     |    6 +
 METHODOLOGY_REVIEW.md            |   83 +-
 diff_diff/conley.py              |   12 +-
 diff_diff/guides/llms-full.txt   |    4 +-
 diff_diff/linalg.py              |   87 +-
 docs/api/spillover.rst           |    6 +-
 docs/methodology/REGISTRY.md     |   13 +-
 tests/test_conley_vcov.py        | 1331 ++---------------------
 tests/test_methodology_conley.py | 1709 ++++++++++++++++++++++++++++++
 9 files changed, 1962 insertions(+), 1289 deletions(-)
 create mode 100644 tests/test_methodology_conley.py

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4c34e2d8..a4f5a2b2 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,9 +7,15 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+### Added
+- **ConleySpatialHAC methodology-review-tracker promotion: In Progress → Complete.** Closes the Conley (1999) *Journal of Econometrics* 92(1) primary-source review on the methodology-review tracker. The paper review on file at `docs/methodology/papers/conley-1999-review.md` was previously merged (2026-05-09); this PR is the F.L.I.P. consolidation — new `tests/test_methodology_conley.py` with paper-equation-numbered Verified Components walk-through (~1600 LoC; 10 classes; 60 tests, 5 of them `@pytest.mark.slow`). Coverage: Eq. 4.2 cross-sectional sandwich (pairwise-distance specialization; the project's paper review identifies Eq. 4.2 page 18 as the real-valued/pairwise form, with Eq. 3.13 reserved for the lattice-indexed form), Eq. 4.2 HC0 + rank-1 limits, Andrews (1991) HAC lag truncation matching `conleyreg::time_dist.cpp`, haversine convention with Earth radius 6371.01 km, Phase 2 panel block-decomposed sandwich at `atol=1e-12`, sparse k-d-tree dense-vs-sparse bit-identity (Wave A #120 numerical correctness), and R `conleyreg` v0.1.9 parity at `atol=1e-6` on 6 fixtures (3 cross-sectional + 3 panel) plus the sparse-forced and time-asymmetric kernel parity contracts. Three dedicated deviations-area classes: `TestConleyLibraryExtensions` (Wave A library extensions — combined spatial+cluster product kernel #119, callable conley_metric validation #123, sparse k-d-tree activation #120, indefiniteness guard), `TestConleyDeviationsFromR` (1-D radial Bartlett vs paper's 2-D separable Eq. 3.14, time-label normalization via `np.unique`, independent temporal kernel deferred), and `TestConleyDeferrals` (5 fail-closed `NotImplementedError`/`TypeError` contracts: LinearRegression + survey_design, DiD/MPD/TWFE + survey_design, Conley + weights, SyntheticDiD + Conley, wild_bootstrap + Conley). Methodology-anchored tests extracted from `tests/test_conley_vcov.py`: full classes `TestConleyDirectHelper`, `TestConleyReductions`, `TestConleyReductionsAddendum`, `TestConleyParityR`, `TestConleyParitySpacetime`, `TestConleyPanelHelper`, `TestConleySparseRParityForced`; plus methodology-anchored tests from `TestConleyKernels`, `TestConleyDistanceMetrics`, `TestConleySparse`. File drops 4248 → 3113 lines after extraction. Defensive surface preserved: input validation, NaN/inf guards, dispatch-level validity, estimator-level integration smoke tests, set_params atomicity, sparse-path activation thresholds + density-gate fallback. `METHODOLOGY_REVIEW.md` row L91 promoted to **Complete** with `Last Review = 2026-05-26`; detail block rewritten with Verified Components / Test Coverage / R Comparison Results inline table / Corrections Made / Deviations / Outstanding Concerns. Priority queue at L1386 pruned: PreTrendsPower removed (already Complete since 2026-05-19) and ConleySpatialHAC removed (this PR); substantive-review-blocked renumbered #2-#5 → #1-#4 and consolidation-pass-blocked renumbered #6-#8 → #5-#6.
+
 ### Added / Changed
 - **EfficientDiD `vcov_type` threading + Results metadata harmonization (Phase 1b interstitial #4, permanently narrow).** `EfficientDiD(vcov_type=...)` now accepts `{"hc1"}` only (default). Analytical-sandwich families `{classical, hc2, hc2_bm}` and `conley` are REJECTED at `__init__` / `set_params` with methodology-rooted messages — EfficientDiD uses influence-function-based variance per Chen-Sant'Anna-Xie (2025) achieving the semiparametric efficiency bound; the per-unit EIF aggregation has no single design matrix on which hat-matrix leverage or Bell-McCaffrey Satterthwaite DOF can be defined. `cluster=` (Liang-Zeger CR1 on cluster-aggregated EIF) and `survey_design=` (TSL on combined IF) paths are unchanged. **BC break on `EfficientDiDResults`:** the `cluster` field renamed to `cluster_name`; new `n_clusters` + `vcov_type` fields added; `to_dict()` method added (mirrors TripleDifferenceResults). `DiagnosticReport._pt_hausman` updated to read the renamed `cluster_name` field for the Hausman pretest replay (`diff_diff/diagnostic_report.py:2444`). `EfficientDiD.set_params(vcov_type=bad)` raises immediately rather than deferring to `fit()` — intentional eager-validation pattern matching EfficientDiD's existing handling of `pt_assumption`/`control_group` etc, diverging from `ImputationDiD`/`TripleDifference`/`CallawaySantAnna` (which use sklearn mutate-then-validate-at-use). Survey-PSU bootstrap path returns NaN SE when fewer than 2 independent PSUs are available (was ≈0 SE from BLAS roundoff). New summary block: `Variance estimator: <label>` line rendered after the survey block when not under bootstrap; suppressed under bootstrap (replaced with `Inference method: bootstrap` + `Bootstrap replications: <n>`). Default `cluster=None` (no survey) renders "HC1 heteroskedasticity-robust" — methodologically correct because the per-unit EIF SE `sqrt(mean(EIF²)/n)` is HC1-style (no Liang-Zeger G/(G-1) finite-sample correction); diverges from `ImputationDiD` which auto-clusters at unit per BJS Theorem 3.
 
+### Fixed
+- **Bertanha-Imbens 2014 citation correction (16 sites across 5 files).** A verification spike confirmed the citation across `diff_diff/linalg.py` (×8), `diff_diff/conley.py` (×1), `diff_diff/guides/llms-full.txt` (×2), `docs/methodology/REGISTRY.md` (×4), and `docs/api/spillover.rst` (×1) was incorrect — NBER w20773 *External Validity in Fuzzy Regression Discontinuity Designs* (JBES 2020, 38(3):593-612) by Bertanha & Imbens covers fuzzy RD external validity, NOT weighted spatial-HAC under sampling weights. Replaced across all 16 sites with the open-problem framing: "weighted spatial-HAC under probability sampling is an open methodological question; no canonical extension of Conley (1999) exists for the combination." At the four `REGISTRY.md` sites the replacement is wrapped in the canonical `**Note (open methodological question):**` label per CLAUDE.md "Documenting Deviations (AI Review Compatibility)". REGISTRY ConleySpatialHAC section gains a new `**Note (deferral status, 2026-05-26):**` splitting the boundary into three parts: **Shipped** — SpilloverDiD + Conley + survey via Wave E.1/E.2/E.3 (PR #468/#474/#482), TwoStageDiD + Conley + survey via Wave E.3 parity (PR #485). **Deferred (generic linalg surface, any `weight_type`)** — DiD/MPD/TWFE/LinearRegression generic path + Conley + `survey_design=`; `LinearRegression` / `compute_robust_vcov` Conley + `weights=` rejected for `pweight`, `aweight`, AND `fweight` (weighted Conley is not implemented on the generic linalg surface). **Open methodological question (subset)** — the `pweight` / `survey_design` portion of the deferral additionally lacks a canonical methodological extension of Conley (1999) for weighted spatial-HAC under probability sampling. **No source-code logic changes:** verified via diff-in-diff pytest output before and after the citation strip (175 passed + 14 warnings, bit-identical pass set on `tests/test_conley_vcov.py`). **Historical CHANGELOG entries (pre-[Unreleased]) intentionally retain the Bertanha-Imbens 2014 attribution** as accurate records of what was claimed at the time of each release; the [Unreleased] entry above supersedes those rationales going forward.
+
 ## [3.4.2] - 2026-05-25
 
 ### Fixed
diff --git a/METHODOLOGY_REVIEW.md b/METHODOLOGY_REVIEW.md
index 9f6e99ad..88678d52 100644
--- a/METHODOLOGY_REVIEW.md
+++ b/METHODOLOGY_REVIEW.md
@@ -88,7 +88,7 @@ The catalog grew incrementally over several quarters, so formats vary across the
 
 | Feature | Module | Reference | Status | Last Review |
 |---------|--------|-----------|--------|-------------|
-| ConleySpatialHAC | `conley.py`, `linalg.py` | `conleyreg` (R) / `acreg` (Stata) | **In Progress** | — |
+| ConleySpatialHAC | `conley.py`, `linalg.py` | `conleyreg` (R) / `acreg` (Stata) | **Complete** | 2026-05-26 |
 | Survey Data Support | `survey.py`, `bootstrap_utils.py` | `survey` package (R) | **In Progress** | — |
 
 **Status legend** (matches the contract in [§ What "Complete" means in this tracker](#what-complete-means-in-this-tracker) above):
@@ -1266,26 +1266,57 @@ These are not estimators but variance/inference plumbing used across many estima
 | Module | `conley.py`, `linalg.py` (`_validate_vcov_args`, kernel construction) |
 | Primary Reference | Conley (1999), *GMM Estimation with Cross-Sectional Dependence*, J. Econometrics 92(1), 1-45 |
 | Secondary References | Andrews (1991) HAC theory; Colella, Lalive, Sakalli & Thoenig (2019) for the Stata `acreg` parallel; Düsterhöft (2021) `conleyreg` (CRAN) parity target |
-| Status | **In Progress** |
-| Last Review | — |
+| Status | **Complete** |
+| Last Review | 2026-05-26 |
 
-**Documentation in place:**
-- REGISTRY.md section: `## ConleySpatialHAC` plus three sub-sections (combined spatial + cluster product kernel — Wave A #119; performance/scale — Wave A #120; callable `conley_metric` validation — Wave A #123)
-- **Paper review on file**: `docs/methodology/papers/conley-1999-review.md` (review date 2026-05-09); plus four adjacent paper reviews for the spillover initiative: `butts-2021-review.md`, `butts-2023-review.md` (JUE Insight), `clarke-2017-review.md`, `colella-et-al-2019-review.md`
-- Implementation: 162 tests in `tests/test_conley_vcov.py` (Phase 1 + Phase 2 space-time HAC)
-- Wired through `DifferenceInDifferences`, `MultiPeriodDiD`, `TwoWayFixedEffects` via `vcov_type="conley"` enum
+**Verified Components:**
+- [x] **Eq. 4.2 cross-sectional sandwich (pairwise-distance specialization)** `Var(β) = (X'X)^{-1} (Σ_{i,j} K(d_ij/h) X_i ε_i ε_j X_j') (X'X)^{-1}`. diff-diff implements the real-valued / pairwise-distance form (Conley 1999 Eq. 4.2 plus the "pairwise products at a given distance" remark on page 19); Eq. 3.13 is the lattice-indexed form reserved for grid coordinates — `tests/test_methodology_conley.py::TestConleyEquation42`
+- [x] **Eq. 4.2 limits**: tiny-cutoff → HC0 diagonal; huge-cutoff under uniform → rank-1 correlated limit; K(0) = 1 diagonal contribution always exact — `TestConleyHC0AndRank1Reductions`
+- [x] **Andrews (1991) HAC lag truncation** `(1 - |t-s|/(L+1))` for `0 < |t-s| ≤ L`, matching `conleyreg::time_dist.cpp`; lag 0 excluded to avoid double-counting — `TestConleyAndrewsLagTruncation`
+- [x] **Haversine convention**: Earth radius 6371.01 km matches `conleyreg::haversine_dist`; 1° lat = 111.195 km at equator — `TestConleyHaversineConvention`
+- [x] **Phase 2 panel block-decomposed sandwich** `XeeX = XeeX_spatial + XeeX_serial` matches `conleyreg::time_dist.cpp` at `atol=1e-12`; cluster time-invariance constraint validated on panel path — `TestConleyBlockDecomposition`
+- [x] **Wave A #120 sparse k-d-tree path** produces meat bit-identical to dense at `atol=1e-10` (chord-projection roundoff on haversine absorbed) — `TestConleySparseDenseEquivalence`
 
-**Documentation in place (R parity):**
-- R `conleyreg` goldens committed: `benchmarks/data/r_conleyreg_conley_golden.json`, generator `benchmarks/R/generate_conley_golden.R`
-- Cross-sectional R parity at `atol=1e-6`: `tests/test_conley_vcov.py::TestConleyParityR`
-- Panel (space-time) R parity at `atol=1e-6`: `TestConleyParitySpacetime` (dense path) and `TestConleySparseRParityForced` (sparse path forced)
-- Internal block-decomposition cross-check at machine precision (matches `conleyreg::time_dist.cpp`): `TestConleyParitySpacetime::test_panel_matches_block_decomposed_reference` (inner tolerance `atol=1e-12`)
+**Test Coverage:**
+- `tests/test_methodology_conley.py` (~1600 LoC; 10 paper-anchored / R-parity / deviations classes; 60 tests including 5 `@pytest.mark.slow`)
+- `tests/test_conley_vcov.py` (~3100 LoC remaining after extraction; 11 classes — defensive surface: input validation, NaN/inf guards, dispatch-level validity, estimator-level integration smoke tests, set_params atomicity, sparse-path activation thresholds + density-gate fallback)
 
-**Outstanding for promotion:**
-- Dedicated `tests/test_methodology_conley.py` with paper-equation-numbered Verified Components walk-through (Equation 8 score-covariance, Bartlett kernel, Andrews-style truncation) consolidating the parity tests into a methodology checklist
-- Summary R-parity table in this tracker (currently the parity results are scattered across class-level docstrings in `tests/test_conley_vcov.py`)
-- Document deviation: indefiniteness guard applied to both spatial and cluster kernels (vs. Bartlett's PSD property)
-- Resolution for the Phase 5 spillover-conley dependency on survey-weights interaction (currently raises `NotImplementedError` at the linalg validator)
+**R Comparison Results** (R `conleyreg` v0.1.9, Düsterhöft 2021):
+
+| Fixture | Surface | Tolerance | Test |
+|---------|---------|-----------|------|
+| `small_haversine` | Cross-sectional | `atol=1e-6`, `rtol=1e-6` | `TestConleyParityR::test_parity_small_haversine` |
+| `dense_haversine` | Cross-sectional | `atol=1e-6`, `rtol=1e-6` | `TestConleyParityR::test_parity_dense_haversine` |
+| `lat_lon_realistic` | Cross-sectional | `atol=1e-6`, `rtol=1e-6` | `TestConleyParityR::test_parity_lat_lon_realistic` |
+| `panel_haversine_lag1` | Panel (Phase 2) | `atol=1e-6`, `rtol=1e-6` | `TestConleyParityR::test_parity_panel_haversine_lag1` |
+| `panel_haversine_lag2` | Panel (Phase 2) | `atol=1e-6`, `rtol=1e-6` | `TestConleyParityR::test_parity_panel_haversine_lag2` |
+| `panel_lat_lon_realistic_lag1` | Panel (Phase 2) | `atol=1e-6`, `rtol=1e-6` | `TestConleyParityR::test_parity_panel_lat_lon_realistic_lag1` |
+| Sparse k-d-tree (forced) on cross-sectional fixtures | Cross-sectional | `atol=1e-6`, `rtol=1e-6` | `TestConleyParityR::test_sparse_forced_matches_r_cross_sectional` |
+| Internal block-decomposition cross-check | Phase 2 algebraic identity vs `conleyreg::time_dist.cpp` | `atol=1e-12` | `TestConleyBlockDecomposition::test_panel_matches_block_decomposed_reference` |
+| Time-asymmetric kernel literal-matching (spatial kernel does NOT propagate to temporal component) | Phase 2 contract | `atol=1e-10` | `TestConleyParityR::test_time_asymmetric_kernel_matches_r_literal` |
+
+Goldens at `benchmarks/data/r_conleyreg_conley_golden.json`; generator at `benchmarks/R/generate_conley_golden.R` (Earth radius 6371.01 km matches `conleyreg::haversine_dist`).
+
+**Corrections Made:**
+- **PR (this PR):** Closed Outstanding-for-promotion items via the new methodology test file (`tests/test_methodology_conley.py`), the inline R-parity summary table above, the consolidated deviations-area classes (`TestConleyLibraryExtensions` + `TestConleyDeviationsFromR` + `TestConleyDeferrals`), and the Phase 5 reframing.
+- **Cited-paper correction (this PR):** Stripped Bertanha-Imbens 2014 citation across 16 sites (`linalg.py` × 8, `conley.py` × 1, `llms-full.txt` × 2, `REGISTRY.md` × 4, `spillover.rst` × 1). NBER w20773 is *External Validity in Fuzzy Regression Discontinuity Designs* (JBES 2020), unrelated to weighted spatial-HAC. Replaced with: "weighted spatial-HAC under probability sampling is an open methodological question; no canonical extension of Conley (1999) exists for the combination." At REGISTRY sites the replacement is wrapped in the canonical `**Note (open methodological question):**` label per CLAUDE.md "Documenting Deviations (AI Review Compatibility)".
+
+**Deviations from the paper / from R / library extensions** (consolidated; see REGISTRY `## ConleySpatialHAC` § Note (deviation / source specialization) for full text):
+- **1-D radial Bartlett vs. paper's 2-D separable Eq. 3.14 PSD-guaranteed form** — practitioner specialization matching R `conleyreg`, Stata `acreg`, Hsiang (2010); not formally PSD-guaranteed. Indefiniteness guard at `-1e-12` applies to both spatial and cluster kernels (REGISTRY L3609). `TestConleyDeviationsFromR::test_1d_radial_bartlett_vs_2d_separable_eq314` + `TestConleyLibraryExtensions::test_indefiniteness_guard_fires_on_negative_eigenvalue`.
+- **Combined spatial + cluster product kernel** (Wave A #119, library extension; no R correspondence; two limit-fixture anchors): `K_total[i, j] = K_space(d_ij/h) · 1{c_i = c_j}`. Cluster time-invariance contract enforced on panel path. Anchor 1 (all-unique-clusters → HC0 diagonal): `TestConleyLibraryExtensions::test_combined_spatial_cluster_kernel_wave_a_119_all_unique`. Anchor 2 (huge-cutoff + uniform → CR1 without Liang-Zeger correction): `TestConleyLibraryExtensions::test_combined_spatial_cluster_kernel_wave_a_119_huge_cutoff`.
+- **Callable `conley_metric` validation** (Wave A #123, library extension; no R correspondence): 6 invariants checked (`(n,n)` shape, finite, non-negative, symmetric to `atol=1e-10`, zero diagonal, float64-castable). `TestConleyLibraryExtensions::test_callable_metric_validation_wave_a_123`.
+- **Sparse k-d-tree fast path** (Wave A #120, library extension; auto-activates at `n > 5,000` with Bartlett + haversine/euclidean; density gate at 30% falls back to dense). `TestConleyLibraryExtensions::test_sparse_kd_tree_activation_wave_a_120` (activation contract); `TestConleySparseDenseEquivalence` (numerical equivalence).
+- **Time-label normalization** via `np.unique(return_inverse=True)` (REGISTRY L3592-3603): diff-diff normalizes time labels to dense panel-period codes before lag computation; R `conleyreg` uses raw time values literally. diff-diff is the more robust default on non-dense encodings. `TestConleyDeviationsFromR::test_time_label_normalization_deviation_from_r`.
+- **Time-asymmetric kernel** (matches R `conleyreg` literal; library extension would be follow-up): R's `kernel` argument controls the spatial component only; the temporal kernel is unconditionally Bartlett. diff-diff matches this asymmetry exactly. Parity contract in `TestConleyParityR::test_time_asymmetric_kernel_matches_r_literal`; deferral in `TestConleyDeviationsFromR::test_independent_temporal_kernel_deferred`.
+
+**Outstanding Concerns:**
+- **Spillover-conley dependency: resolved.** SpilloverDiD ships Conley + survey via Wave E.1/E.2/E.3 (PR #468 / #474 / #482, stratified-Conley sandwich on PSU totals with within-PSU serial Bartlett HAC for `lag_cutoff > 0`); TwoStageDiD ships Wave E.3 parity (PR #485, `fdf2cebc`).
+- **Generic LinearRegression / DiD / MPD / TWFE + Conley + survey_design: deferred.** Two fail-closed contracts assert the unsupported combination: the estimator-level gate (DiD / MPD / TWFE) lives in `diff_diff/conley.py::_validate_conley_estimator_inputs` (`TestConleyDeferrals::test_did_mpd_twfe_survey_design_not_implemented`); the generic LinearRegression gate lives in `diff_diff/linalg.py::LinearRegression.fit` (`TestConleyDeferrals::test_linear_regression_survey_design_not_implemented`). The `survey_design=` surface is `LinearRegression` only; `compute_robust_vcov` does not accept `survey_design=`.
+- **Generic LinearRegression / compute_robust_vcov + Conley + `weights=`: deferred for any `weight_type` (`pweight` / `aweight` / `fweight`).** Weighted Conley is not implemented on the generic linalg surface. The `pweight` / `survey_design` subset additionally reflects an **open methodological question** — no canonical extension of Conley (1999) exists for weighted spatial-HAC under probability sampling. `TestConleyDeferrals::test_conley_plus_weights_not_implemented`.
+- **SyntheticDiD + Conley**: raises `TypeError`. SyntheticDiD uses bootstrap / jackknife / placebo variance, not the analytical sandwich. `TestConleyDeferrals::test_synthetic_did_conley_typeerror`.
+- **Wild bootstrap + Conley**: `NotImplementedError`. Wild bootstrap is a separate inference path that does not consume the analytical Conley sandwich (REGISTRY L3547). `TestConleyDeferrals::test_wild_bootstrap_conley_not_implemented`.
+- **No default bandwidth**: Conley 1999 doesn't propose a plug-in selector; REGISTRY recommends `(50, 100, 200, 500)` km sensitivity grid mirroring Conley 1999 § 5. `tests/test_conley_vcov.py::TestConleyValidatorHelpers::test_missing_cutoff_raises` enforces the fail-closed contract.
+- **DiagnosticReport routing** for `SpilloverDiDResults(vcov_type="conley", survey_design=)` is queued for a follow-up Wave F PR — `_APPLICABILITY` / `_PT_METHOD` registration is required before the new combination can be claimed consumable downstream.
 
 ---
 
@@ -1352,23 +1383,21 @@ whereas R's `did::att_gt` would error. This is a defensive enhancement that prov
 more graceful handling of edge cases while still signaling invalid inference to users.
 ```
 
-### Priority Order (2026-05-15)
+### Priority Order (2026-05-26)
 
 Promotion priority for the **In Progress** entries, ordered by what's blocked on substantive review work (top of list = needs review next) vs. consolidation pass (bottom of list = mostly tracker walk-through):
 
 **Substantive-review-blocked (no methodology test file, no paper review, no R parity):**
 
-1. **PreTrendsPower** — small surface, established R package (`pretrends`), Roth (2022) is short.
-2. **PowerAnalysis** — larger surface (MDE / power / sample size / simulation paths); REGISTRY already lists Bloom (1995) and Burlig et al. (2020) as primary sources; least urgent if the library's power-analysis utilities are not heavily used.
-3. **PlaceboTests** — decide first whether to keep standalone or absorb into per-estimator diagnostic sections; methodologically lightweight either way.
-4. **EfficientDiD** — no paper review on file; substantial implementation work (`tests/test_efficient_did.py` + validation tests) needs paper-vs-code audit against Chen, Sant'Anna & Xie (2025).
-5. **ImputationDiD / TwoStageDiD** — natural pair (both single-treatment-effect-imputation methods). Each needs paper review, methodology file, R parity fixture against `didimputation` / `did2s`.
+1. **PowerAnalysis** — larger surface (MDE / power / sample size / simulation paths); REGISTRY already lists Bloom (1995) and Burlig et al. (2020) as primary sources; least urgent if the library's power-analysis utilities are not heavily used.
+2. **PlaceboTests** — decide first whether to keep standalone or absorb into per-estimator diagnostic sections; methodologically lightweight either way.
+3. **EfficientDiD** — no paper review on file; substantial implementation work (`tests/test_efficient_did.py` + validation tests) needs paper-vs-code audit against Chen, Sant'Anna & Xie (2025).
+4. **ImputationDiD / TwoStageDiD** — natural pair (both single-treatment-effect-imputation methods). Each needs paper review, methodology file, R parity fixture against `didimputation` / `did2s`.
 
 **Consolidation-pass-blocked (already has paper review or methodology file or R parity; mostly Verified Components walk-through):**
 
-6. **StaggeredTripleDifference** — shares the primary paper (Ortiz-Villavicencio & Sant'Anna 2025) with TripleDifference, but no dedicated paper review on file yet; needs R parity (R fixtures gitignored — tracked in TODO.md, PR #245).
-7. **ConleySpatialHAC** — paper review + committed R `conleyreg` goldens; needs dedicated methodology test file + summary R-parity table in this tracker.
-8. **Survey Data Support** — cross-cutting feature; promotion requires the per-estimator integration paths to be locked down first.
+5. **StaggeredTripleDifference** — shares the primary paper (Ortiz-Villavicencio & Sant'Anna 2025) with TripleDifference, but no dedicated paper review on file yet; needs R parity (R fixtures gitignored — tracked in TODO.md, PR #245).
+6. **Survey Data Support** — cross-cutting feature; promotion requires the per-estimator integration paths to be locked down first.
 
 ---
 
diff --git a/diff_diff/conley.py b/diff_diff/conley.py
index 4deceff2..855b85a5 100644
--- a/diff_diff/conley.py
+++ b/diff_diff/conley.py
@@ -296,10 +296,14 @@ def _validate_conley_estimator_inputs(
         raise ValueError(f"Cluster column '{cluster}' not found in data")
     if survey_design is not None:
         raise NotImplementedError(
-            f"{estimator_name}(vcov_type='conley') + survey_design is a "
-            "follow-up (Bertanha-Imbens 2014 weighted-Conley). Drop "
-            "survey_design for cross-sectional Conley, or use vcov_type='hc1' "
-            "for survey-aware cluster-robust without spatial HAC."
+            f"{estimator_name}(vcov_type='conley') + survey_design is "
+            "deferred — weighted spatial-HAC under probability sampling "
+            "is an open methodological question; no canonical extension "
+            "of Conley (1999) exists for the combination. Drop "
+            "survey_design for unweighted Conley (cross-sectional or panel "
+            "block-decomposed via conley_lag_cutoff > 0), or use "
+            "vcov_type='hc1' for survey-aware cluster-robust without "
+            "spatial HAC."
         )
     if inference == "wild_bootstrap":
         raise NotImplementedError(
diff --git a/diff_diff/guides/llms-full.txt b/diff_diff/guides/llms-full.txt
index 8307f280..f444a5eb 100644
--- a/diff_diff/guides/llms-full.txt
+++ b/diff_diff/guides/llms-full.txt
@@ -2122,9 +2122,9 @@ SE range.
 - `DifferenceInDifferences` / `MultiPeriodDiD` `+ vcov_type="conley"` without `unit=` at fit-time → `ValueError`.
 - Combining `vcov_type="conley"` with explicit `cluster=<col>` applies the combined spatial + cluster product kernel (Wave A #119). On the panel path the cluster must be constant within each unit across periods (validator raises `ValueError` otherwise). TWFE's default auto-cluster on the Conley path is silently dropped; users opt into the combined kernel explicitly.
 - `DifferenceInDifferences` / `MultiPeriodDiD` / `TwoWayFixedEffects` `(vcov_type="conley", inference="wild_bootstrap")` → `NotImplementedError` (wild bootstrap does not consume the analytical sandwich).
-- `DifferenceInDifferences` / `MultiPeriodDiD` / `TwoWayFixedEffects` `(vcov_type="conley")` + `survey_design=` → `NotImplementedError` (Bertanha-Imbens 2014 follow-up).
+- `DifferenceInDifferences` / `MultiPeriodDiD` / `TwoWayFixedEffects` `(vcov_type="conley")` + `survey_design=` → `NotImplementedError` (estimator-level gate; weighted spatial-HAC under probability sampling is an open methodological question).
 - `SyntheticDiD(vcov_type="conley")` → `TypeError` (uses bootstrap, not analytical sandwich).
-- Any-mode `vcov_type="conley"` `+ weights=` / `survey_design=` on `LinearRegression` / `compute_robust_vcov` → `NotImplementedError` (Bertanha-Imbens 2014 weighted-Conley deferred to follow-up PR).
+- Generic `LinearRegression(vcov_type="conley", survey_design=...)` → `NotImplementedError`. Generic `LinearRegression` / `compute_robust_vcov` Conley rejects `weights=` for any `weight_type` (`pweight` / `aweight` / `fweight`) → `NotImplementedError` (weighted Conley is not implemented on the generic linalg surface; `compute_robust_vcov` does not accept `survey_design=`, the survey-design surface is `LinearRegression` only). The `pweight` / `survey_design` subset additionally reflects an open methodological question — no canonical extension of Conley (1999) exists for weighted spatial-HAC under probability sampling.
 - A sparse k-d-tree fast path auto-activates for `n > 5_000` with `conley_kernel="bartlett"` AND `conley_metric` in `{"haversine", "euclidean"}` (Wave A #120); callable metrics and uniform kernel fall back to the dense path. `n > 20_000` with the dense fallback still emits a memory-OOM `UserWarning`.
 - Callable `conley_metric` is validated at the boundary (shape `(n, n)`, finite, non-negative, symmetric to `atol=1e-10`, AND zero on the diagonal `|d(i, i)| ≤ 1e-10` so the kernel reduces to `K(0) = 1` on the HC0 contribution); each failure raises a targeted `ValueError` (Wave A #123).
 
diff --git a/diff_diff/linalg.py b/diff_diff/linalg.py
index 63e763f2..9098a927 100644
--- a/diff_diff/linalg.py
+++ b/diff_diff/linalg.py
@@ -553,8 +553,12 @@ def solve_ols(
           ``conley_lag_cutoff``. Combining with ``cluster_ids`` applies the
           combined spatial + cluster product kernel (a diff-diff convention;
           see :func:`compute_robust_vcov` for details). Combining with
-          ``weights`` raises ``NotImplementedError`` (Bertanha-Imbens 2014
-          weighted-Conley deferred to a follow-up PR).
+          ``weights`` raises ``NotImplementedError`` regardless of
+          ``weight_type``: weighted Conley is not implemented on the generic
+          linalg surface. For probability-sampling weights (``pweight`` /
+          ``survey_design``) the deferral additionally reflects an open
+          methodological question — no canonical extension of Conley (1999)
+          exists for weighted spatial-HAC under probability sampling.
     conley_coords : ndarray of shape (n, 2), optional
         Required when ``vcov_type="conley"``. Two-column array of
         ``[lat, lon]`` (degrees, for ``conley_metric="haversine"``) or
@@ -1119,13 +1123,16 @@ def _validate_vcov_args(
         combined with a ``vcov_type`` that is one-way only (``classical``,
         ``hc2``).
     NotImplementedError
-        If ``vcov_type == "conley"`` is combined with ``weights`` (the
-        Bertanha-Imbens 2014 weighted-Conley methodology is a separate
-        follow-up). NOT raised here for ``hc2_bm + weights``: that
-        weight-type contract is enforced downstream in
-        ``_compute_robust_vcov_numpy`` (which has access to ``weight_type``
-        and rejects ``aweight`` / ``fweight`` while routing ``pweight``
-        through the clubSandwich WLS-CR2 port).
+        If ``vcov_type == "conley"`` is combined with ``weights`` (regardless
+        of ``weight_type``: weighted Conley is not implemented on the
+        generic linalg surface). For ``pweight`` / probability-sampling
+        designs the deferral additionally reflects an open methodological
+        question — no canonical extension of Conley (1999) exists for
+        weighted spatial-HAC under probability sampling. NOT raised here for
+        ``hc2_bm + weights``: that weight-type contract is enforced
+        downstream in ``_compute_robust_vcov_numpy`` (which has access to
+        ``weight_type`` and rejects ``aweight`` / ``fweight`` while routing
+        ``pweight`` through the clubSandwich WLS-CR2 port).
     """
     if vcov_type not in _VALID_VCOV_TYPES:
         raise ValueError(
@@ -1154,14 +1161,23 @@ def _validate_vcov_args(
         # Conley + cluster_ids is now supported (combined spatial + cluster
         # product kernel; see ``docs/methodology/REGISTRY.md`` § ConleySpatialHAC
         # → "Combined spatial + cluster product kernel"). Conley + weights
-        # (Bertanha-Imbens 2014) is still deferred to a follow-up PR — the
-        # design-based weighted-Conley methodology requires its own
-        # treatment.
+        # remains deferred regardless of weight_type — weighted Conley is
+        # not implemented on the generic linalg surface; for probability-
+        # sampling weights the deferral additionally reflects an open
+        # methodological question with no canonical extension of Conley
+        # (1999) for the combination.
         if weights is not None:
             raise NotImplementedError(
-                "vcov_type='conley' with weights is a Phase 2+ follow-up "
-                "(Bertanha-Imbens 2014). Drop weights for cross-sectional "
-                "Conley, or use vcov_type='hc1' for weighted HC1."
+                "vcov_type='conley' with weights is not implemented on the "
+                "generic linalg surface (any weight_type — pweight, "
+                "aweight, or fweight). For probability-sampling weights "
+                "(pweight / survey_design), the deferral additionally "
+                "reflects an open methodological question: no canonical "
+                "extension of Conley (1999) exists for weighted spatial-"
+                "HAC under probability sampling. Drop weights for "
+                "unweighted Conley (cross-sectional or panel block-"
+                "decomposed via conley_lag_cutoff > 0), or use "
+                "vcov_type='hc1' for weighted HC1."
             )
 
 
@@ -1275,8 +1291,12 @@ def compute_robust_vcov(
       ``K_total[i, j] = K_space(d_ij/h) · 1{c_i = c_j}`` (Wave A #119; on
       the panel path the cluster must be constant within each unit across
       periods). Combining with ``weights`` still raises
-      ``NotImplementedError`` (Bertanha-Imbens 2014 weighted-Conley
-      deferred to a follow-up PR).
+      ``NotImplementedError`` regardless of ``weight_type`` (weighted
+      Conley is not implemented on the generic linalg surface); for
+      probability-sampling weights (``pweight`` / ``survey_design``) the
+      deferral additionally reflects an open methodological question, with
+      no canonical extension of Conley (1999) for weighted spatial-HAC
+      under probability sampling.
 
     Parameters
     ----------
@@ -2992,9 +3012,16 @@ class LinearRegression:
         a positive ``conley_cutoff_km``. Combining ``vcov_type="conley"``
         with ``cluster_ids`` applies the combined spatial + cluster product
         kernel (Wave A #119; cluster must be constant within each unit on
-        the panel path). Combining with ``weights`` / ``survey_design``
-        still raises ``NotImplementedError`` (Bertanha-Imbens 2014
-        weighted-Conley deferred to a follow-up PR). The DiD / MPD / TWFE
+        the panel path). Combining with ``weights`` raises
+        ``NotImplementedError`` regardless of ``weight_type`` (weighted
+        Conley is not implemented on the generic linalg surface);
+        combining with ``survey_design`` (``LinearRegression`` only;
+        ``compute_robust_vcov`` has no survey-design surface) likewise
+        raises ``NotImplementedError``. For probability-sampling weights
+        (``pweight`` / ``survey_design``) the deferral additionally
+        reflects an open methodological question, with no canonical
+        extension of Conley (1999) for weighted spatial-HAC under
+        probability sampling. The DiD / MPD / TWFE
         estimators all support panel Conley by passing ``unit`` at fit-time
         (DiD as a fit-time kwarg; MPD/TWFE via the existing ``unit``
         argument), threading ``conley_time`` / ``conley_unit`` into the
@@ -3232,16 +3259,22 @@ def fit(
         # is set to False on the solve_ols call), and the survey vcov path
         # would silently overwrite the result with a non-Conley variance
         # under a Conley request. Front-door the rejection here so the
-        # contract is enforced uniformly. Phase 5 (Bertanha-Imbens 2014
-        # weighted-Conley) will lift this; Phase 1 supports cross-sectional
-        # unweighted Conley only.
+        # contract is enforced uniformly. Weighted spatial-HAC under
+        # probability sampling is an open methodological question (no
+        # canonical extension of Conley (1999) exists for the combination);
+        # the generic LinearRegression / compute_robust_vcov path supports
+        # unweighted Conley only (cross-sectional or panel block-decomposed
+        # via conley_lag_cutoff > 0).
         if _fit_vcov_type == "conley" and _use_survey_vcov:
             raise NotImplementedError(
                 "LinearRegression(vcov_type='conley', survey_design=...) "
-                "is deferred to Phase 5 (Bertanha-Imbens 2014 weighted-"
-                "Conley). Phase 1 supports cross-sectional unweighted "
-                "Conley only via compute_robust_vcov / LinearRegression "
-                "without a survey design."
+                "is deferred — weighted spatial-HAC under probability "
+                "sampling is an open methodological question; no "
+                "canonical extension of Conley (1999) exists for the "
+                "combination. The generic LinearRegression / "
+                "compute_robust_vcov path supports unweighted Conley only "
+                "(cross-sectional or panel block-decomposed via "
+                "conley_lag_cutoff > 0) without a survey design."
             )
 
         # Resolve effective fit-time weights/weight_type WITHOUT mutating
diff --git a/docs/api/spillover.rst b/docs/api/spillover.rst
index 611f0fb4..9b62c47c 100644
--- a/docs/api/spillover.rst
+++ b/docs/api/spillover.rst
@@ -316,8 +316,10 @@ and planned follow-up enhancements:
     pattern; the Conley cluster product kernel becomes a no-op after
     PSU aggregation.
   - The LinearRegression-side ``vcov_type="conley" + survey_design=``
-    gate at ``diff_diff/linalg.py`` is a separate Bertanha-Imbens 2014
-    weighted-Conley roadmap (not Wave E).
+    gate at ``diff_diff/linalg.py`` remains a separate roadmap (not Wave E)
+    — weighted spatial-HAC under probability sampling is an open
+    methodological question; no canonical extension of Conley (1999) exists
+    for the combination.
   - DiagnosticReport routing for ``SpilloverDiDResults(vcov_type="conley",
     survey_design=)`` is queued for a follow-up (the
     ``_APPLICABILITY`` / ``_PT_METHOD`` wiring must register the new
diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md
index 001d7ffb..101cc3c0 100644
--- a/docs/methodology/REGISTRY.md
+++ b/docs/methodology/REGISTRY.md
@@ -3422,7 +3422,7 @@ Degrees of freedom for the t-distribution lookup use `ResolvedSurveyDesign.df_su
 
 - **Restrictions / out-of-scope (Wave E.2):**
   - Replicate-weight variance (BRR / Fay / JK1 / JKn / SDR) raises `NotImplementedError` (inherits Wave E.1 gate; per-replicate full refit is separate follow-up scope).
-  - LinearRegression-side `vcov_type="conley" + survey_design=` gate at `diff_diff/linalg.py:2853` remains (separate Bertanha-Imbens 2014 weighted-Conley "Phase 5" roadmap, not Wave E).
+  - LinearRegression-side `vcov_type="conley" + survey_design=` gate in `LinearRegression.fit()` (Conley + survey rejection block) remains. **Note (open methodological question):** weighted spatial-HAC under probability sampling is an open methodological question; no canonical extension of Conley (1999) exists for the combination — separate roadmap, not Wave E.
   - DiagnosticReport routing for `SpilloverDiDResults(vcov_type="conley", survey_design=)` is queued for a follow-up Wave F PR — `_APPLICABILITY` / `_PT_METHOD` registration is required before the new combination can be claimed consumable downstream (per `feedback_audit_diagnostic_report_wiring_before_claim`).
 
 **Implementation:** new `_compute_stratified_conley_meat_from_psu_scores` helper in `diff_diff/survey.py` (parallel to existing Binder helper; 3-tuple `(meat, variance_computed, legitimate_zero_count)` return contract; per-stratum loop replaces the inner `centered.T @ centered` with `_compute_conley_meat(centered, coords_h, cutoff, metric, kernel)` cross-sectional mode). New dispatch wrapper `_compute_stratified_conley_meat` in `diff_diff/two_stage.py` (parallel to `_compute_binder_tsl_meat`; per-obs Psi → PSU aggregation via `np.add.at` + PSU centroid derivation via vectorized `np.add.at` sums / `np.bincount` counts + dispatch to survey helper; intentionally no `cluster_ids` parameter). `_compute_gmm_corrected_meat` conley branch extended at `diff_diff/two_stage.py` with `if resolved_survey is not None` routing to the new wrapper; the `resolved_survey is None` branch is bit-identical to Wave D no-survey Conley. Saturation NaN-fail mirrors Wave E.1 (`UserWarning` template "Wave E.2 stratified-Conley sandwich: df_survey = 0..."). Wave E.1 stage-1 weighted FE solver, `finite_mask` survey-array subsetting, `df_survey` threading to `safe_inference` call sites, bread weighting, and `SpilloverDiDResults` survey metadata are all inherited UNCHANGED — Psi construction is bit-identical regardless of `vcov_type`. Tests: `TestSpilloverDiDWaveE2ConleySurveyDesign` + `TestSpilloverDiDWaveE2ConleySurveyDesignEventStudy` at `tests/test_spillover.py`.
@@ -3450,7 +3450,7 @@ Degrees of freedom for the t-distribution lookup use `ResolvedSurveyDesign.df_su
 - **Restrictions (Wave E.2 follow-up):**
   - Requires an **effective PSU**: either an explicit `survey_design.psu` OR a `cluster=<col>` argument that gets injected as the effective PSU per Wave E.1's `_inject_cluster_as_psu` routing. No-effective-PSU survey designs (weights-only / strata-only WITHOUT a cluster fallback) raise `NotImplementedError` at `SpilloverDiD.fit` post-resolution because under the pseudo-PSU = obs-index fallback each pseudo-PSU appears in exactly one period — the per-PSU serial cross-period loop would silently contribute zero. Routing the serial loop to `conley_unit` would mix IF allocators with the spatial term's pseudo-PSU aggregation (the user-confirmed methodology pins a single IF allocator); fail-closed per `feedback_no_silent_failures` until a no-effective-PSU derivation is queued (tracked in TODO.md).
   - Replicate-weight variance (BRR / Fay / JK1 / JKn / SDR) raises `NotImplementedError` (inherits Wave E.1 gate).
-  - LinearRegression-side `vcov_type="conley" + survey_design=` gate at `diff_diff/linalg.py:2853` remains (separate Bertanha-Imbens 2014 weighted-Conley roadmap, not Wave E).
+  - LinearRegression-side `vcov_type="conley" + survey_design=` gate in `LinearRegression.fit()` (Conley + survey rejection block) remains. **Note (open methodological question):** weighted spatial-HAC under probability sampling is an open methodological question; no canonical extension of Conley (1999) exists for the combination — separate roadmap, not Wave E.
   - DiagnosticReport routing for `SpilloverDiDResults(vcov_type="conley", conley_lag_cutoff > 0, survey_design=)` is queued for Wave F follow-up (per `feedback_audit_diagnostic_report_wiring_before_claim` — `_APPLICABILITY` / `_PT_METHOD` registration required before claiming consumability).
 
 **Implementation:** new sibling helper `_compute_stratified_serial_bartlett_meat` in `diff_diff/two_stage.py` (parallel to the Wave E.2 spatial orchestrator; ~200 LoC; three-mode singleton-stratum branching with FPC scaling inside the multi-PSU branch to avoid divide-by-zero; panel-wide dense time codes for the lag math). Orchestrator `_compute_stratified_conley_meat` signature extended with `conley_lag_cutoff: Optional[int] = None`; spatial loop unchanged; serial helper called after spatial loop when `conley_lag_cutoff > 0`; saturation NaN-fail accounting merges both terms' `(variance_computed, legitimate_zero)` flags. Dispatch in `_compute_gmm_corrected_meat` conley branch threads `conley_lag_cutoff` through to the orchestrator. Spillover-side gate at `spillover.py:2210` deleted (Wave E.2 era `NotImplementedError` for lag>0 + survey). Stage-1 weighted FE solver, `finite_mask` survey-array subsetting, `df_survey` threading, bread weighting, and `SpilloverDiDResults` survey metadata are all inherited UNCHANGED — Psi construction is bit-identical regardless of vcov_type or lag. Tests: `TestSpilloverDiDWaveE2FollowupConleySurveyLagCutoff` + `TestSpilloverDiDWaveE2FollowupConleySurveyLagCutoffEventStudy` at `tests/test_spillover.py`.
@@ -3764,9 +3764,9 @@ metrics (`"haversine"`, `"euclidean"`) satisfy this by construction.
 - `TwoWayFixedEffects(vcov_type="conley", cluster=<col>)` is supported (Wave A #119): combined spatial + cluster product kernel applies. The cluster must be time-invariant within each unit on the panel path (validator-enforced). TWFE's default auto-cluster is silently dropped on the Conley path; explicit cluster is required to opt in.
 - `DifferenceInDifferences(vcov_type="conley", cluster=<col>)`: combined kernel applies; same time-invariance contract on the panel path. DiD has no auto-cluster, so the cluster choice is fully explicit.
 - `DifferenceInDifferences` / `MultiPeriodDiD` / `TwoWayFixedEffects` `(vcov_type="conley", inference="wild_bootstrap")` → `NotImplementedError`. (MPD's pre-Conley analytical-fallback `UserWarning` is suppressed when `vcov_type="conley"` so the user gets one consistent error message.)
-- `DifferenceInDifferences` / `MultiPeriodDiD` / `TwoWayFixedEffects` `(vcov_type="conley")` + `survey_design=` → `NotImplementedError` at the estimator level (deferred to Bertanha-Imbens 2014 weighted-Conley follow-up).
+- `DifferenceInDifferences` / `MultiPeriodDiD` / `TwoWayFixedEffects` `(vcov_type="conley")` + `survey_design=` → `NotImplementedError` at the estimator level. **Note (open methodological question):** weighted spatial-HAC under probability sampling is an open methodological question; no canonical extension of Conley (1999) exists for the combination.
 - `SyntheticDiD(vcov_type="conley")` → `TypeError` (SyntheticDiD uses bootstrap/jackknife/placebo variance, not the analytical sandwich; tracked in TODO.md).
-- Any-mode `vcov_type="conley"` `+ weights=` / `survey_design=` → `NotImplementedError` (Bertanha-Imbens 2014 territory; Phase 5 follow-up).
+- Generic `LinearRegression(vcov_type="conley", survey_design=...)` → `NotImplementedError`. Generic `LinearRegression / compute_robust_vcov` Conley rejects `weights=` for any `weight_type` (`pweight` / `aweight` / `fweight`) → `NotImplementedError` (weighted Conley is not implemented on the generic linalg surface). `compute_robust_vcov` does not accept `survey_design=`; the survey-design surface is `LinearRegression` only. **Note (open methodological question):** the `pweight` / `survey_design` subset additionally reflects an open methodological question — no canonical extension of Conley (1999) exists for weighted spatial-HAC under probability sampling. (Estimator-specific shipped surfaces — SpilloverDiD via Wave E.1/E.2/E.3 and TwoStageDiD via Wave E.3 parity — are explicitly excepted; see "Note (deferral status, 2026-05-26)" below.)
 - Panel path: partial `conley_time` / `conley_unit` / `conley_lag_cutoff` (not all three set) → `ValueError` at validator.
 - Panel path with `cluster_ids` that vary across periods within a unit → `ValueError` (time-invariance contract).
 - `n > 20_000` with `conley_kernel='uniform'` or callable metric: emits `UserWarning` about O(n²) distance-matrix memory (sparse fast path doesn't apply; consider switching to bartlett or projecting to euclidean for performance).
@@ -3774,6 +3774,11 @@ metrics (`"haversine"`, `"euclidean"`) satisfy this by construction.
 - Identical coordinates (`d_ij = 0` for `i ≠ j`): `K(0) = 1`, contributing the full HC0 weight per Conley 1999 page 19. Documented behavior; no warning.
 - Callable `conley_metric` returning a non-(n,n)/NaN/inf/negative/asymmetric/non-zero-diagonal matrix raises `ValueError` naming the violated invariant. The zero-diagonal contract (`|d(i, i)| ≤ 1e-10`) is load-bearing for the Conley sandwich's HC0 reduction `K(0) = 1`; see "Callable conley_metric validation" subsection above.
 
+**Note (deferral status, 2026-05-26):** Conley + survey_design / weights surface boundary —
+- **Shipped**: SpilloverDiD + Conley + survey via Wave E.1/E.2/E.3 (PR #468 / #474 / #482, stratified-Conley sandwich on PSU totals with within-PSU serial Bartlett HAC for `lag_cutoff > 0`); TwoStageDiD + Conley + survey via Wave E.3 parity (PR #485 / `fdf2cebc`).
+- **Deferred (generic linalg surface — any weight_type)**: DiD / MPD / TWFE / LinearRegression generic path + Conley + `survey_design=`; `LinearRegression` / `compute_robust_vcov` Conley + `weights=` (rejected for `pweight`, `aweight`, AND `fweight`; weighted Conley is not implemented on the generic linalg surface).
+- **Open methodological question (subset)**: the `pweight` / `survey_design` portion of the deferral additionally has no canonical methodological extension — weighted spatial-HAC under probability sampling lacks a canonical extension of Conley (1999).
+
 **Reference implementations:**
 - R: `conleyreg::conleyreg(...)` (Düsterhöft 2021, CRAN v0.1.9) — **parity benchmark for diff-diff**
 - Stata: `acreg dep indeps, latitude(lat) longitude(lon) spatial dist(km) bartlett` (Colella et al. 2019) — academic spec source for the cluster-flexible variance, not parity-tested here
diff --git a/tests/test_conley_vcov.py b/tests/test_conley_vcov.py
index 2843c064..88602c79 100644
--- a/tests/test_conley_vcov.py
+++ b/tests/test_conley_vcov.py
@@ -1,9 +1,27 @@
-"""Tests for Conley (1999) spatial HAC variance estimator.
-
-Phase 1 scope: pure-numerics helpers (kernels, distance metrics, direct
-sandwich helper) and the dispatch-level validator. Estimator-level
-integration tests, set_params atomicity, and Stata acreg parity land in
-later Phase 1 checkpoints (Steps 4-6 of the plan).
+"""Defensive regression surface for the Conley (1999) spatial HAC
+variance estimator.
+
+Companion to ``tests/test_methodology_conley.py`` (paper-equation-numbered
+methodology walk-through + R `conleyreg` v0.1.9 parity at ``atol=1e-6``
++ library extensions / deviations / deferrals area). This file retains:
+
+- Input validation (callable-metric invariants, lat/lon range bounds,
+  unknown-metric / unknown-kernel raises)
+- Dispatch-level validity (vcov_type membership; cluster + one-way
+  rejections; Conley + weights / survey_design fail-closed gates)
+- Estimator-level integration smoke tests (DiD + Conley, TWFE + Conley,
+  MPD + Conley) and set_params atomicity
+- Sparse-path activation thresholds + density-gate fallback behavior
+  (the methodology file pins numerical equivalence vs dense; activation
+  and fallback are defensive)
+- PSD-guard infrastructure (`_validate_meat_psd` stacklevel attribution,
+  threshold boundary behavior)
+
+Paper-equation-anchored methodology coverage (Eq. 4.2 sandwich, Andrews
+1991 lag truncation, haversine convention, R parity goldens, sparse-vs-
+dense bit-identity, combined-kernel limit anchors, indefiniteness guard
+on both kernels, library-extension / deviation / deferral contracts)
+lives in ``tests/test_methodology_conley.py``.
 """
 
 import warnings
@@ -16,11 +34,7 @@
     _CONLEY_SPARSE_N_THRESHOLD,
     _bartlett_kernel,
     _compute_conley_vcov,
-    _compute_spatial_bartlett_meat_sparse,
-    _haversine_km,
     _pairwise_distance_matrix,
-    _serial_bartlett_kernel_matrix,
-    _uniform_kernel,
     _validate_conley_kwargs,
     _validate_meat_psd,
 )
@@ -56,168 +70,22 @@ def small_ols_with_coords():
 
 
 # ---------------------------------------------------------------------------
-# TestConleyKernels — Bartlett / uniform shape and boundary behavior
+# Note: Methodology-anchored kernel form tests (Bartlett K(u), uniform kernel,
+# serial Bartlett HAC matrix) + paper-anchored haversine/distance-metric
+# numerical contracts have moved to tests/test_methodology_conley.py
+# (TestConleyEquation42, TestConleyAndrewsLagTruncation,
+# TestConleyHaversineConvention). The defensive surface for input validation
+# (callable metric invariant checks, unknown-metric / unknown-kernel raises,
+# lat/lon range validation) remains here.
 # ---------------------------------------------------------------------------
 
 
-class TestConleyKernels:
-    def test_bartlett_at_zero(self):
-        np.testing.assert_allclose(_bartlett_kernel(np.array([0.0])), 1.0)
-
-    def test_bartlett_at_one(self):
-        np.testing.assert_allclose(_bartlett_kernel(np.array([1.0])), 0.0)
-
-    def test_bartlett_above_one_zero(self):
-        u = np.array([1.5, 2.0, 100.0])
-        np.testing.assert_allclose(_bartlett_kernel(u), np.zeros(3))
-
-    def test_bartlett_negative_arg_symmetric(self):
-        """Bartlett uses |u|, so K(-0.3) == K(0.3)."""
-        np.testing.assert_allclose(
-            _bartlett_kernel(np.array([-0.3])), _bartlett_kernel(np.array([0.3]))
-        )
-
-    def test_uniform_kernel_at_boundary(self):
-        """Uniform kernel is closed on the right: K(1) = 1, K(1+eps) = 0."""
-        np.testing.assert_allclose(_uniform_kernel(np.array([1.0])), 1.0)
-
-    def test_uniform_kernel_above_one_zero(self):
-        np.testing.assert_allclose(_uniform_kernel(np.array([1.0001, 2.0, 100.0])), np.zeros(3))
-
-    def test_uniform_kernel_at_zero_one(self):
-        np.testing.assert_allclose(_uniform_kernel(np.array([0.0])), 1.0)
-
-    def test_bartlett_kernel_finite_and_in_unit_interval(self):
-        """Bartlett-weighted kernel matrix on random pairwise distances is
-        finite, symmetric, and bounded in [0, 1]. We do NOT assert PSD here:
-        the radial 1-D Bartlett on pairwise distance is a practitioner
-        specialization of Conley 1999 (matching R conleyreg) and is NOT
-        formally PSD-guaranteed — see REGISTRY ConleySpatialHAC. The
-        runtime path emits a UserWarning if the resulting Conley meat is
-        materially indefinite; that contract is locked separately in
-        ``test_indefinite_meat_warning_fires_for_bartlett``.
-        """
-        rng = np.random.default_rng(seed=11)
-        n = 25
-        coords = rng.uniform(0, 1, size=(n, 2))
-        diff = coords[:, None, :] - coords[None, :, :]
-        D = np.sqrt((diff * diff).sum(axis=-1))
-        K = _bartlett_kernel(D / 0.3)
-        assert K.shape == (n, n)
-        assert np.all(np.isfinite(K))
-        assert np.all(K >= 0.0)
-        assert np.all(K <= 1.0)
-        np.testing.assert_allclose(K, K.T, atol=1e-15)
-
-    def test_serial_bartlett_kernel_matrix_basic(self):
-        """Hand-computed Bartlett HAC kernel for t=[0,1,2,3], L=2.
-
-        K[i,j] = (1 - |i-j|/(L+1)) for 0 < |i-j| <= L, else 0. With L=2
-        and (L+1)=3: lag 1 -> 2/3, lag 2 -> 1/3, lag 3 -> 0 (out of band).
-        """
-        K = _serial_bartlett_kernel_matrix(np.array([0, 1, 2, 3]), L=2)
-        expected = np.array(
-            [
-                [0.0, 2.0 / 3.0, 1.0 / 3.0, 0.0],
-                [2.0 / 3.0, 0.0, 2.0 / 3.0, 1.0 / 3.0],
-                [1.0 / 3.0, 2.0 / 3.0, 0.0, 2.0 / 3.0],
-                [0.0, 1.0 / 3.0, 2.0 / 3.0, 0.0],
-            ]
-        )
-        np.testing.assert_allclose(K, expected, atol=1e-15)
-
-    def test_serial_bartlett_kernel_matrix_l_one(self):
-        """L=1: only adjacent lags survive, with weight (1 - 1/2) = 0.5."""
-        K = _serial_bartlett_kernel_matrix(np.array([0, 1, 2]), L=1)
-        expected = np.array(
-            [
-                [0.0, 0.5, 0.0],
-                [0.5, 0.0, 0.5],
-                [0.0, 0.5, 0.0],
-            ]
-        )
-        np.testing.assert_allclose(K, expected, atol=1e-15)
-
-    def test_serial_bartlett_kernel_matrix_l_zero_returns_zero(self):
-        """L=0 is degenerate-but-callable: every off-diagonal lag fails
-        ``lag <= 0`` (since `lag != 0`), so K is the zero matrix. Callers
-        guard externally (``conley.py`` skips the loop when L == 0)."""
-        K = _serial_bartlett_kernel_matrix(np.array([0, 1, 2]), L=0)
-        np.testing.assert_array_equal(K, np.zeros((3, 3)))
-
-    def test_serial_bartlett_kernel_matrix_single_element(self):
-        """Single-element input yields a 1x1 zero matrix (no off-diagonal
-        lags exist)."""
-        K = _serial_bartlett_kernel_matrix(np.array([7]), L=2)
-        np.testing.assert_array_equal(K, np.zeros((1, 1)))
-
-    def test_serial_bartlett_kernel_matrix_int_input_bit_equal_to_float(self):
-        """Contract test: int64 and float64 inputs must yield bit-equal
-        matrices. The helper does ``astype(np.float64, copy=False)`` and one
-        of the three current call sites (``conley.py`` panel-block branch)
-        passes the result of array slicing on int time codes."""
-        K_int = _serial_bartlett_kernel_matrix(np.array([0, 1, 2, 3], dtype=np.int64), L=2)
-        K_float = _serial_bartlett_kernel_matrix(
-            np.array([0.0, 1.0, 2.0, 3.0], dtype=np.float64), L=2
-        )
-        np.testing.assert_array_equal(K_int, K_float)
-
-
 # ---------------------------------------------------------------------------
-# TestConleyDistanceMetrics — haversine, euclidean, callable
+# TestConleyDistanceMetrics — callable metric validation + unknown-metric raise
 # ---------------------------------------------------------------------------
 
 
 class TestConleyDistanceMetrics:
-    def test_haversine_known_pair_one_degree_equator(self):
-        """1° longitude at the equator = 2π·R/360 ≈ 111.195 km (R=6371)."""
-        d = _haversine_km(np.array(0.0), np.array(0.0), np.array(0.0), np.array(1.0))
-        expected = 2 * np.pi * _CONLEY_EARTH_RADIUS_KM / 360.0
-        np.testing.assert_allclose(d, expected, atol=1e-9)
-
-    def test_haversine_zero_self_distance(self):
-        d = _haversine_km(np.array(45.0), np.array(-122.0), np.array(45.0), np.array(-122.0))
-        np.testing.assert_allclose(d, 0.0, atol=1e-12)
-
-    def test_haversine_symmetric(self):
-        d_ab = _haversine_km(np.array(40.7), np.array(-74.0), np.array(34.0), np.array(-118.2))
-        d_ba = _haversine_km(np.array(34.0), np.array(-118.2), np.array(40.7), np.array(-74.0))
-        np.testing.assert_allclose(d_ab, d_ba, atol=1e-12)
-
-    def test_haversine_pole_to_equator(self):
-        """North pole to equator at any longitude = π/2 · R = ~10007.5 km."""
-        d = _haversine_km(np.array(90.0), np.array(0.0), np.array(0.0), np.array(0.0))
-        expected = np.pi * _CONLEY_EARTH_RADIUS_KM / 2.0
-        np.testing.assert_allclose(d, expected, atol=1e-9)
-
-    def test_haversine_broadcasting_pairwise(self):
-        """Broadcasting (n, 1) vs (1, n) yields the n×n distance matrix."""
-        coords = np.array([[0.0, 0.0], [0.0, 1.0], [0.0, 2.0]])
-        lats = coords[:, 0]
-        lons = coords[:, 1]
-        D = _haversine_km(lats[:, None], lons[:, None], lats[None, :], lons[None, :])
-        assert D.shape == (3, 3)
-        np.testing.assert_allclose(np.diag(D), 0.0, atol=1e-12)
-        # D[0, 2] should be 2 * D[0, 1] for collinear-equator points
-        np.testing.assert_allclose(D[0, 2], 2.0 * D[0, 1], rtol=1e-10)
-
-    def test_pairwise_distance_haversine(self):
-        coords = np.array([[0.0, 0.0], [0.0, 1.0], [10.0, 0.0]])
-        D = _pairwise_distance_matrix(coords, "haversine")
-        assert D.shape == (3, 3)
-        np.testing.assert_allclose(D, D.T, atol=1e-12)
-        np.testing.assert_allclose(np.diag(D), 0.0, atol=1e-12)
-
-    def test_pairwise_distance_euclidean_matches_pdist(self):
-        """Euclidean path matches scipy.spatial.distance squareform exactly."""
-        from scipy.spatial.distance import pdist, squareform
-
-        rng = np.random.default_rng(seed=7)
-        coords = rng.uniform(0, 100, size=(15, 2))
-        D = _pairwise_distance_matrix(coords, "euclidean")
-        D_scipy = squareform(pdist(coords, metric="euclidean"))
-        np.testing.assert_allclose(D, D_scipy, atol=1e-12)
-
     def test_pairwise_distance_callable(self):
         """A user-supplied callable is dispatched and its output preserved.
         Output must satisfy the validator's invariants (zero diagonal, finite,
@@ -963,295 +831,12 @@ def _indefinite(t_codes: np.ndarray, L: int) -> np.ndarray:
 
 
 # ---------------------------------------------------------------------------
-# TestConleyDirectHelper — _compute_conley_vcov correctness
-# ---------------------------------------------------------------------------
-
-
-class TestConleyDirectHelper:
-    def test_returns_psd_with_bartlett(self, small_ols_with_coords):
-        X, residuals, coords, bread = small_ols_with_coords
-        vcov = _compute_conley_vcov(
-            X,
-            residuals,
-            coords,
-            cutoff=2000.0,
-            metric="haversine",
-            kernel="bartlett",
-            bread_matrix=bread,
-        )
-        eigvals = np.linalg.eigvalsh(0.5 * (vcov + vcov.T))
-        assert eigvals.min() > -1e-10
-
-    def test_symmetric_vcov(self, small_ols_with_coords):
-        X, residuals, coords, bread = small_ols_with_coords
-        vcov = _compute_conley_vcov(
-            X,
-            residuals,
-            coords,
-            cutoff=2000.0,
-            metric="haversine",
-            kernel="bartlett",
-            bread_matrix=bread,
-        )
-        np.testing.assert_allclose(vcov, vcov.T, atol=1e-10)
-
-    def test_shape_matches_bread(self, small_ols_with_coords):
-        X, residuals, coords, bread = small_ols_with_coords
-        vcov = _compute_conley_vcov(
-            X,
-            residuals,
-            coords,
-            cutoff=1500.0,
-            metric="haversine",
-            kernel="bartlett",
-            bread_matrix=bread,
-        )
-        k = X.shape[1]
-        assert vcov.shape == (k, k)
-
-    def test_uniform_kernel_negative_eigenvalue_warns(self):
-        """The uniform kernel is documented as not PSD-guaranteed (Conley
-        1999 footnote 11): its spectral window is negative in regions, so
-        the resulting meat can be indefinite. Force the indefinite path
-        deterministically by monkey-patching ``_uniform_kernel`` to return
-        a kernel matrix with aggressive negative off-diagonals (mirroring
-        the bartlett warning test below), and assert the warning surfaces
-        with kernel='uniform' in the message.
-        """
-        from diff_diff import conley as conley_mod
-
-        rng = np.random.default_rng(seed=1)
-        n = 6
-        coords = rng.uniform(0, 1, size=(n, 2))
-        X = np.column_stack([np.ones(n), rng.standard_normal(n)])
-        eps = np.ones(n)
-        bread = X.T @ X
-
-        original = conley_mod._uniform_kernel
-
-        def _indefinite(u: np.ndarray) -> np.ndarray:
-            base = np.eye(u.shape[0])
-            for i in range(u.shape[0]):
-                for j in range(u.shape[0]):
-                    if i != j:
-                        base[i, j] = -10.0
-            return base
-
-        try:
-            conley_mod._uniform_kernel = _indefinite
-            with warnings.catch_warnings(record=True) as w:
-                warnings.simplefilter("always")
-                conley_mod._compute_conley_vcov(
-                    X,
-                    eps,
-                    coords,
-                    cutoff=10.0,
-                    metric="euclidean",
-                    kernel="uniform",
-                    bread_matrix=bread,
-                )
-            psd_warnings = [
-                msg
-                for msg in w
-                if issubclass(msg.category, UserWarning)
-                and "uniform" in str(msg.message)
-                and "negative eigenvalue" in str(msg.message)
-            ]
-            assert len(psd_warnings) >= 1, (
-                f"Expected a UserWarning naming kernel='uniform' and "
-                f"'negative eigenvalue'; got {[str(m.message) for m in w]}"
-            )
-        finally:
-            conley_mod._uniform_kernel = original
-
-    def test_indefinite_meat_warning_fires_for_bartlett(self):
-        """Both kernels (radial 1-D bartlett and uniform) are practitioner
-        specializations of Conley 1999 and are NOT formally PSD-guaranteed
-        (Conley's explicit PSD formula is the 2-D separable product window,
-        Eq 3.14, not the 1-D radial form). The PSD guard must therefore
-        fire for bartlett too, not just uniform.
-
-        Forces the indefinite path by monkey-patching `_bartlett_kernel` to
-        return a kernel matrix with an aggressive negative off-diagonal,
-        making the resulting meat indefinite. Verifies the warning surfaces
-        with the kernel name in the message.
-        """
-        from diff_diff import conley as conley_mod
-
-        rng = np.random.default_rng(seed=11)
-        n = 6
-        coords = rng.uniform(0, 1, size=(n, 2))
-        X = np.column_stack([np.ones(n), rng.standard_normal(n)])
-        eps = np.ones(n)  # non-zero residuals so meat is non-zero
-        bread = X.T @ X
-
-        # Patch the bartlett kernel to return a sign-pattern that DEFINITELY
-        # produces an indefinite meat. The native bartlett is non-negative;
-        # injecting large negative off-diagonals breaks the
-        # PSD-by-non-negative-window heuristic.
-        original = conley_mod._bartlett_kernel
-
-        def _indefinite(u: np.ndarray) -> np.ndarray:
-            base = np.eye(u.shape[0])
-            # Aggressive negative off-diagonals; this kernel is itself
-            # indefinite and so is S.T @ K @ S for generic S.
-            for i in range(u.shape[0]):
-                for j in range(u.shape[0]):
-                    if i != j:
-                        base[i, j] = -10.0
-            return base
-
-        try:
-            conley_mod._bartlett_kernel = _indefinite
-            with warnings.catch_warnings(record=True) as w:
-                warnings.simplefilter("always")
-                conley_mod._compute_conley_vcov(
-                    X,
-                    eps,
-                    coords,
-                    cutoff=10.0,
-                    metric="euclidean",
-                    kernel="bartlett",
-                    bread_matrix=bread,
-                )
-            # Verify a PSD warning fired naming the bartlett kernel
-            psd_warnings = [
-                msg
-                for msg in w
-                if issubclass(msg.category, UserWarning)
-                and "bartlett" in str(msg.message)
-                and "negative eigenvalue" in str(msg.message)
-            ]
-            assert len(psd_warnings) >= 1, (
-                f"Expected a UserWarning naming kernel='bartlett' and "
-                f"'negative eigenvalue'; got {[str(m.message) for m in w]}"
-            )
-        finally:
-            conley_mod._bartlett_kernel = original
-
-
-# ---------------------------------------------------------------------------
-# TestConleyReductions — Bartlett+tiny cutoff → HC0 meat; etc.
+# TestConleyValidationDispatch — dispatch-level validation contracts.
+# (Methodology-anchored Eq. 4.2 PSD / symmetric / shape contracts and
+# HC0+rank-1 reductions moved to tests/test_methodology_conley.py.)
 # ---------------------------------------------------------------------------
 
 
-class TestConleyReductions:
-    def test_tiny_cutoff_distinct_coords_yields_HC0_meat(self):
-        """When the bandwidth is much smaller than the minimum pairwise
-        distance, Conley's kernel is ~0 off-diagonal and 1 on-diagonal, so
-        the meat reduces to Σ x_i² ε_i² x_i x_i' = HC0 meat.
-        """
-        rng = np.random.default_rng(seed=3)
-        n = 15
-        # Distinct coords with min pairwise distance >> 0
-        coords = np.column_stack([np.arange(n) * 100.0, np.arange(n) * 100.0])
-        X = np.column_stack([np.ones(n), rng.standard_normal(n)])
-        eps = rng.standard_normal(n)
-        bread = X.T @ X
-
-        # HC0 meat (Σ x_i x_i' u_i²) — no DOF correction applied
-        meat_hc0 = X.T @ (X * (eps**2)[:, None])
-        bread_inv = np.linalg.solve(bread, np.eye(2))
-        vcov_hc0 = bread_inv @ meat_hc0 @ bread_inv
-
-        vcov_conley = _compute_conley_vcov(
-            X,
-            eps,
-            coords,
-            cutoff=1.0,  # << minimum pairwise distance
-            metric="euclidean",
-            kernel="bartlett",
-            bread_matrix=bread,
-        )
-        np.testing.assert_allclose(vcov_conley, vcov_hc0, atol=1e-12)
-
-    def test_uniform_kernel_tiny_cutoff_yields_HC0_meat(self):
-        """Same reduction with the uniform kernel."""
-        rng = np.random.default_rng(seed=5)
-        n = 12
-        coords = np.column_stack([np.arange(n) * 100.0, np.arange(n) * 100.0])
-        X = np.column_stack([np.ones(n), rng.standard_normal(n)])
-        eps = rng.standard_normal(n)
-        bread = X.T @ X
-
-        meat_hc0 = X.T @ (X * (eps**2)[:, None])
-        bread_inv = np.linalg.solve(bread, np.eye(2))
-        vcov_hc0 = bread_inv @ meat_hc0 @ bread_inv
-
-        vcov_conley = _compute_conley_vcov(
-            X,
-            eps,
-            coords,
-            cutoff=0.5,
-            metric="euclidean",
-            kernel="uniform",
-            bread_matrix=bread,
-        )
-        np.testing.assert_allclose(vcov_conley, vcov_hc0, atol=1e-12)
-
-    def test_huge_cutoff_NOT_HC0(self, small_ols_with_coords):
-        """When cutoff dwarfs all pairwise distances, K -> ones(n, n) and
-        meat = (X·ε)' ones (X·ε) which is the rank-1 outer product of summed
-        scores — NOT HC0. This is the all-correlated limit."""
-        X, residuals, coords, bread = small_ols_with_coords
-        vcov_conley = _compute_conley_vcov(
-            X,
-            residuals,
-            coords,
-            cutoff=1e9,
-            metric="euclidean",
-            kernel="uniform",
-            bread_matrix=bread,
-        )
-        # HC0 for comparison
-        meat_hc0 = X.T @ (X * (residuals**2)[:, None])
-        bread_inv = np.linalg.solve(bread, np.eye(X.shape[1]))
-        vcov_hc0 = bread_inv @ meat_hc0 @ bread_inv
-        # They should differ noticeably
-        assert not np.allclose(vcov_conley, vcov_hc0, atol=1e-6)
-
-    def test_dispatch_matches_direct_helper(self, small_ols_with_coords):
-        """compute_robust_vcov(vcov_type='conley', ...) returns the same
-        vcov as a direct call to _compute_conley_vcov on the same inputs.
-        Atol=1e-14 (bit-equivalence)."""
-        X, residuals, coords, bread = small_ols_with_coords
-        vcov_dispatch = compute_robust_vcov(
-            X,
-            residuals,
-            vcov_type="conley",
-            conley_coords=coords,
-            conley_cutoff_km=2000.0,
-            conley_metric="haversine",
-            conley_kernel="bartlett",
-        )
-        vcov_direct = _compute_conley_vcov(
-            X,
-            residuals,
-            coords,
-            2000.0,
-            "haversine",
-            "bartlett",
-            bread,
-        )
-        np.testing.assert_allclose(vcov_dispatch, vcov_direct, atol=1e-14, rtol=1e-14)
-
-    def test_dispatch_returns_dof_when_requested(self, small_ols_with_coords):
-        """return_dof=True returns (vcov, dof_vec) tuple where dof = n - k."""
-        X, residuals, coords, _ = small_ols_with_coords
-        out = compute_robust_vcov(
-            X,
-            residuals,
-            vcov_type="conley",
-            conley_coords=coords,
-            conley_cutoff_km=2000.0,
-            return_dof=True,
-        )
-        assert isinstance(out, tuple) and len(out) == 2
-        _vcov, dof = out
-        n, k = X.shape
-        np.testing.assert_array_equal(dof, np.full(k, n - k, dtype=np.float64))
-
-
 class TestConleyValidationDispatch:
     """Validation tests at the compute_robust_vcov dispatch level."""
 
@@ -1402,8 +987,10 @@ def test_linear_regression_conley_with_survey_design_raises(self, fit_data):
         rejection: it sets `return_vcov=False` on the solve_ols call when
         survey vcov is needed, skipping the linalg validator, and the survey
         vcov path then overwrites `vcov_` with a non-Conley variance under a
-        Conley request. Phase 5 will lift this rejection (Bertanha-Imbens 2014
-        weighted-Conley); Phase 1 is unweighted only.
+        Conley request. The rejection is permanent on the generic
+        LinearRegression path — weighted spatial-HAC under probability
+        sampling is an open methodological question (no canonical extension
+        of Conley (1999) exists for the combination).
         """
         from diff_diff.survey import make_pweight_design
 
@@ -2606,8 +2193,10 @@ def test_did_conley_combinations(self, df):
         )
 
     def test_did_conley_with_survey_design_raises(self, df):
-        """DiD + Conley + survey_design raises NotImplementedError (deferred
-        to Wave 2 weighted-Conley)."""
+        """DiD + Conley + survey_design raises NotImplementedError — weighted
+        spatial-HAC under probability sampling is an open methodological
+        question (no canonical extension of Conley (1999) exists for the
+        combination)."""
         from diff_diff import DifferenceInDifferences, SurveyDesign
 
         with pytest.raises(NotImplementedError, match="conley.*survey_design"):
@@ -2726,579 +2315,89 @@ def test_valid_kwargs_apply(self):
         assert est.conley_kernel == "uniform"
 
 
-class TestConleyParityR:
-    """R conleyreg parity for the Conley spatial HAC implementation.
+class TestConleySparse:
+    """Sparse k-d-tree fast path for the spatial Bartlett meat.
 
-    Skips when the golden JSON is absent (CI's isolated-install job copies
-    only tests/, not benchmarks/). Local regeneration:
-        cd benchmarks/R && Rscript generate_conley_golden.R
+    The sparse path is gated by three conditions: total n above the
+    threshold, metric in {"haversine", "euclidean"} (no callable), and
+    kernel == "bartlett". Each of these tests exercises one of the
+    gates plus the bit-identity parity claim vs the dense path.
     """
 
-    GOLDEN_PATH = "benchmarks/data/r_conleyreg_conley_golden.json"
-    PARITY_TOL = 1e-6  # Phase 1 success criterion
-
-    @pytest.fixture(scope="class")
-    def golden(self):
-        import json
-        from pathlib import Path
-
-        repo_root = Path(__file__).resolve().parent.parent
-        path = repo_root / self.GOLDEN_PATH
-        if not path.exists():
-            pytest.skip(
-                f"Golden JSON not present at {path}; run "
-                "`cd benchmarks/R && Rscript generate_conley_golden.R` to generate. "
-                "Requires conleyreg R package + sf/lwgeom + system libs gdal/proj/geos/udunits."
-            )
-        return json.loads(path.read_text())
-
-    def _check_fixture(self, golden, name):
-        entry = golden[name]
-        X = np.asarray(entry["x"], dtype=np.float64).reshape(entry["x_shape"])
-        y = np.asarray(entry["y"], dtype=np.float64)
-        coords = np.asarray(entry["coords"], dtype=np.float64).reshape(entry["coords_shape"])
-        vcov_expected = np.asarray(entry["vcov"], dtype=np.float64).reshape(entry["vcov_shape"])
-
+    def _euclidean_fixture(self, n=1000, k=3, cutoff=15.0, seed=11):
+        rng = np.random.default_rng(seed)
+        coords = rng.uniform(0.0, 100.0, size=(n, 2))
+        X = np.column_stack([np.ones(n)] + [rng.standard_normal(n) for _ in range(k - 1)])
+        beta = np.linspace(0.5, 2.0, k)
+        y = X @ beta + rng.standard_normal(n) * 0.5
         coefs, *_ = np.linalg.lstsq(X, y, rcond=None)
         residuals = y - X @ coefs
-        vcov_got = compute_robust_vcov(
-            X,
-            residuals,
-            vcov_type="conley",
-            conley_coords=coords,
-            conley_cutoff_km=entry["cutoff_km"],
-            conley_metric=entry["metric"],
-            conley_kernel=entry["kernel"],
-        )
-        np.testing.assert_allclose(
-            vcov_got, vcov_expected, atol=self.PARITY_TOL, rtol=self.PARITY_TOL
-        )
+        bread = X.T @ X
+        return X, residuals, coords, bread, cutoff
 
-    def test_parity_small_haversine(self, golden):
-        self._check_fixture(golden, "small_haversine")
+    def _haversine_fixture(self, n=1000, k=3, cutoff_km=500.0, seed=13):
+        rng = np.random.default_rng(seed)
+        lats = rng.uniform(-30.0, 30.0, size=n)
+        lons = rng.uniform(-100.0, 100.0, size=n)
+        coords = np.column_stack([lats, lons])
+        X = np.column_stack([np.ones(n)] + [rng.standard_normal(n) for _ in range(k - 1)])
+        beta = np.linspace(0.5, 2.0, k)
+        y = X @ beta + rng.standard_normal(n) * 0.5
+        coefs, *_ = np.linalg.lstsq(X, y, rcond=None)
+        residuals = y - X @ coefs
+        bread = X.T @ X
+        return X, residuals, coords, bread, cutoff_km
 
-    def test_parity_dense_haversine(self, golden):
-        self._check_fixture(golden, "dense_haversine")
+    def test_auto_toggle_above_threshold_uses_sparse(self, monkeypatch):
+        """n > _CONLEY_SPARSE_N_THRESHOLD with bartlett + euclidean must
+        auto-route through the sparse helper. Verified by spying on the
+        sparse helper call count."""
+        import diff_diff.conley as conley_module
 
-    def test_parity_lat_lon_realistic(self, golden):
-        self._check_fixture(golden, "lat_lon_realistic")
+        X, residuals, coords, bread, cutoff = self._euclidean_fixture(
+            n=_CONLEY_SPARSE_N_THRESHOLD + 1
+        )
+        calls = {"n": 0}
+        orig = conley_module._compute_spatial_bartlett_meat_sparse
 
+        def _spy(*args, **kwargs):
+            calls["n"] += 1
+            return orig(*args, **kwargs)
 
-class TestConleyParitySpacetime:
-    """R conleyreg parity on the Phase 2 block-decomposed panel form.
+        monkeypatch.setattr(conley_module, "_compute_spatial_bartlett_meat_sparse", _spy)
+        _compute_conley_vcov(X, residuals, coords, cutoff, "euclidean", "bartlett", bread)
+        assert calls["n"] >= 1, "Sparse helper not called when n > threshold."
 
-    Each fixture has lag_cutoff > 0 and exercises the additive sandwich
-    (within-period spatial + within-unit Bartlett serial). Earth radius
-    6371.01 km. Parity target: atol=1e-6.
-    """
+    def test_auto_toggle_below_threshold_stays_dense(self, monkeypatch):
+        """n <= _CONLEY_SPARSE_N_THRESHOLD must use the dense path even
+        when other sparse conditions (bartlett + euclidean) are met."""
+        import diff_diff.conley as conley_module
 
-    GOLDEN_PATH = "benchmarks/data/r_conleyreg_conley_golden.json"
-    PARITY_TOL = 1e-6
+        X, residuals, coords, bread, cutoff = self._euclidean_fixture(n=_CONLEY_SPARSE_N_THRESHOLD)
+        calls = {"n": 0}
+        orig = conley_module._compute_spatial_bartlett_meat_sparse
 
-    @pytest.fixture(scope="class")
-    def golden(self):
-        import json
-        from pathlib import Path
+        def _spy(*args, **kwargs):
+            calls["n"] += 1
+            return orig(*args, **kwargs)
 
-        repo_root = Path(__file__).resolve().parent.parent
-        path = repo_root / self.GOLDEN_PATH
-        if not path.exists():
-            pytest.skip(
-                f"Golden JSON not present at {path}; run "
-                "`cd benchmarks/R && Rscript generate_conley_golden.R` to generate."
-            )
-        return json.loads(path.read_text())
+        monkeypatch.setattr(conley_module, "_compute_spatial_bartlett_meat_sparse", _spy)
+        _compute_conley_vcov(X, residuals, coords, cutoff, "euclidean", "bartlett", bread)
+        assert calls["n"] == 0, "Sparse helper called below threshold."
 
-    def _check_panel_fixture(self, golden, name):
-        entry = golden[name]
-        X = np.asarray(entry["x"], dtype=np.float64).reshape(entry["x_shape"])
-        y = np.asarray(entry["y"], dtype=np.float64)
-        coords = np.asarray(entry["coords"], dtype=np.float64).reshape(entry["coords_shape"])
-        vcov_expected = np.asarray(entry["vcov"], dtype=np.float64).reshape(entry["vcov_shape"])
-        unit = np.asarray(entry["unit"])
-        time = np.asarray(entry["time"])
+    def test_auto_toggle_callable_metric_stays_dense(self, monkeypatch):
+        """A callable conley_metric forces the dense path even at large n —
+        the kd-tree query needs a vectorizable metric, and callables are
+        not supported via projection."""
+        import diff_diff.conley as conley_module
 
-        coefs, *_ = np.linalg.lstsq(X, y, rcond=None)
-        residuals = y - X @ coefs
-        vcov_got = compute_robust_vcov(
-            X,
-            residuals,
-            vcov_type="conley",
-            conley_coords=coords,
-            conley_cutoff_km=entry["cutoff_km"],
-            conley_metric=entry["metric"],
-            conley_kernel=entry["kernel"],
-            conley_time=time,
-            conley_unit=unit,
-            conley_lag_cutoff=int(entry["lag_cutoff"]),
-        )
-        np.testing.assert_allclose(
-            vcov_got, vcov_expected, atol=self.PARITY_TOL, rtol=self.PARITY_TOL
+        X, residuals, coords, bread, cutoff = self._euclidean_fixture(
+            n=_CONLEY_SPARSE_N_THRESHOLD + 100
         )
 
-    def test_parity_panel_haversine_lag1(self, golden):
-        self._check_panel_fixture(golden, "panel_haversine_lag1")
-
-    def test_parity_panel_haversine_lag2(self, golden):
-        self._check_panel_fixture(golden, "panel_haversine_lag2")
-
-    def test_parity_panel_lat_lon_realistic_lag1(self, golden):
-        self._check_panel_fixture(golden, "panel_lat_lon_realistic_lag1")
-
-
-class TestConleyReductionsAddendum:
-    """Additional reduction tests not covered by the helper-direct class.
-
-    Placeholder: the helper-direct class already covers the essential
-    reductions (HC0 at tiny cutoff, K=ones at huge cutoff, etc.).
-    Kept here so future test expansions have a clear class to attach to.
-    """
-
-    def test_diagonal_of_meat_equals_HC0_contribution(self):
-        """For any kernel, K(0/h) = 1 so the diagonal contribution to the
-        meat is exactly the HC0 term Σ_i X_i ε_i² X_i'."""
-        rng = np.random.default_rng(seed=9)
-        n = 20
-        coords = rng.uniform(0, 1000, size=(n, 2))
-        X = np.column_stack([np.ones(n), rng.standard_normal(n)])
-        eps = rng.standard_normal(n)
-        # Build kernel with cutoff between min and max pairwise distance
-        D = _pairwise_distance_matrix(coords, "euclidean")
-        cutoff = float(D[D > 0].min() * 0.001)  # ensure off-diagonal kernel is 0
-        # With this cutoff, the Bartlett kernel is 1 on the diagonal and 0 off,
-        # so meat == HC0.
-        S = X * eps[:, None]
-        meat_full = S.T @ _bartlett_kernel(D / cutoff) @ S
-        meat_hc0 = X.T @ (X * (eps**2)[:, None])
-        np.testing.assert_allclose(meat_full, meat_hc0, atol=1e-12)
-
-
-# ---------------------------------------------------------------------------
-# TestConleyPanelHelper — _compute_conley_vcov with the block-decomposed
-# panel path (R conleyreg lag_cutoff > 0 form).
-# ---------------------------------------------------------------------------
-
-
-class TestConleyPanelHelper:
-    """The Phase 2 panel block-decomposed form (matches R conleyreg)."""
-
-    def _panel_fixture(self, n_units=5, T=3, k=2, seed=42, cutoff_km=5000):
-        rng = np.random.default_rng(seed)
-        lat_unit = rng.uniform(-30, 30, size=n_units)
-        lon_unit = rng.uniform(-100, 100, size=n_units)
-        unit = np.repeat(np.arange(n_units), T)
-        time = np.tile(np.arange(1, T + 1), n_units)
-        lat = lat_unit[unit]
-        lon = lon_unit[unit]
-        coords = np.column_stack([lat, lon])
-        n = n_units * T
-        X = np.column_stack([np.ones(n)] + [rng.standard_normal(n) for _ in range(k - 1)])
-        beta = np.linspace(0.5, 2.0, k)
-        y = X @ beta + rng.standard_normal(n) * 0.5
-        beta_hat, *_ = np.linalg.lstsq(X, y, rcond=None)
-        residuals = y - X @ beta_hat
-        bread = X.T @ X
-        return X, residuals, coords, time, unit, bread, cutoff_km
-
-    def test_T_eq_1_equals_cross_sectional(self):
-        """Block-decomposed form with T=1 (single time period) should equal
-        the Phase 1 cross-sectional form on the same data."""
-        X, residuals, coords, _, _, bread, cutoff = self._panel_fixture(n_units=8, T=1, k=2)
-        unit_single = np.arange(X.shape[0])
-        time_single = np.ones(X.shape[0], dtype=int)
-        # Phase 1 cross-sectional
-        V_cs = _compute_conley_vcov(X, residuals, coords, cutoff, "haversine", "bartlett", bread)
-        # Phase 2 panel block-decomposed with T=1 and lag_cutoff > 0
-        # (the serial component has nothing to do at T=1 — only one period)
-        V_panel = _compute_conley_vcov(
-            X,
-            residuals,
-            coords,
-            cutoff,
-            "haversine",
-            "bartlett",
-            bread,
-            time=time_single,
-            unit=unit_single,
-            lag_cutoff=2,
-        )
-        np.testing.assert_allclose(V_panel, V_cs, atol=1e-12)
-
-    def test_lag_cutoff_zero_drops_serial(self):
-        """lag_cutoff=0 means the serial component contributes nothing;
-        only the within-period spatial sandwich applies."""
-        X, residuals, coords, time, unit, bread, cutoff = self._panel_fixture()
-        # lag_cutoff=0
-        V0 = _compute_conley_vcov(
-            X,
-            residuals,
-            coords,
-            cutoff,
-            "haversine",
-            "bartlett",
-            bread,
-            time=time,
-            unit=unit,
-            lag_cutoff=0,
-        )
-        # Manually compute the within-period spatial sandwich
-        S = X * residuals[:, None]
-        meat_spatial = np.zeros((X.shape[1], X.shape[1]))
-        for t_val in np.unique(time):
-            mask = time == t_val
-            D_t = _pairwise_distance_matrix(coords[mask], "haversine")
-            K_t = _bartlett_kernel(D_t / cutoff)
-            meat_spatial += S[mask].T @ K_t @ S[mask]
-        V_expected = np.linalg.solve(bread, meat_spatial)
-        V_expected = np.linalg.solve(bread, V_expected.T).T
-        np.testing.assert_allclose(V0, V_expected, atol=1e-12)
-
-    def test_lag_cutoff_positive_adds_serial(self):
-        """lag_cutoff > 0 strictly increases the meat by a positive contribution
-        (off-diagonals from within-unit cross-time pairs)."""
-        X, residuals, coords, time, unit, bread, cutoff = self._panel_fixture()
-        V0 = _compute_conley_vcov(
-            X,
-            residuals,
-            coords,
-            cutoff,
-            "haversine",
-            "bartlett",
-            bread,
-            time=time,
-            unit=unit,
-            lag_cutoff=0,
-        )
-        V1 = _compute_conley_vcov(
-            X,
-            residuals,
-            coords,
-            cutoff,
-            "haversine",
-            "bartlett",
-            bread,
-            time=time,
-            unit=unit,
-            lag_cutoff=1,
-        )
-        # The serial sandwich adds non-trivial off-diagonal contributions
-        # for the panel fixture; V1 should differ from V0
-        assert not np.allclose(
-            V0, V1, atol=1e-8
-        ), "lag_cutoff=1 must differ from lag_cutoff=0 with a serial component"
-
-    def test_panel_matches_block_decomposed_reference(self):
-        """Direct verification that _compute_conley_vcov matches the
-        hand-coded block decomposition from time_dist.cpp at machine precision."""
-        X, residuals, coords, time, unit, bread, cutoff = self._panel_fixture(seed=314)
-        bread_inv = np.linalg.inv(bread)
-        S = X * residuals[:, None]
-        # Hand-coded reference (matches R conleyreg per the spike)
-        for L in (0, 1, 2):
-            meat = np.zeros((X.shape[1], X.shape[1]))
-            for t_val in np.unique(time):
-                mask = time == t_val
-                D_t = _pairwise_distance_matrix(coords[mask], "haversine")
-                K_t = _bartlett_kernel(D_t / cutoff)
-                meat += S[mask].T @ K_t @ S[mask]
-            if L > 0:
-                for u_val in np.unique(unit):
-                    mask = unit == u_val
-                    S_u = S[mask]
-                    t_u = time[mask].astype(np.float64)
-                    lag = np.abs(t_u[:, None] - t_u[None, :])
-                    K_u = ((lag <= L) & (lag != 0)).astype(np.float64) * (1.0 - lag / (L + 1.0))
-                    meat += S_u.T @ K_u @ S_u
-            V_ref = bread_inv @ meat @ bread_inv
-
-            V_helper = _compute_conley_vcov(
-                X,
-                residuals,
-                coords,
-                cutoff,
-                "haversine",
-                "bartlett",
-                bread,
-                time=time,
-                unit=unit,
-                lag_cutoff=L,
-            )
-            np.testing.assert_allclose(V_helper, V_ref, atol=1e-12)
-
-    def test_time_label_normalization_non_unit_spaced_int(self):
-        """Year-like int labels (2020, 2021, 2022) and YYYYMM labels
-        (202011, 202012, 202101) produce the same vcov as the equivalent
-        dense codes (0, 1, 2). Closes Codex P1: `conley_lag_cutoff` is a
-        count of panel periods, not raw label difference."""
-        X, residuals, coords, _, unit, bread, cutoff = self._panel_fixture(n_units=8, T=3, k=2)
-        time_dense = np.tile([1, 2, 3], 8)
-        time_years = np.tile([2020, 2021, 2022], 8)
-        time_yyyymm = np.tile([202011, 202012, 202101], 8)
-        V_dense = _compute_conley_vcov(
-            X,
-            residuals,
-            coords,
-            cutoff,
-            "haversine",
-            "bartlett",
-            bread,
-            time=time_dense,
-            unit=unit,
-            lag_cutoff=1,
-        )
-        for time_alt in (time_years, time_yyyymm):
-            V_alt = _compute_conley_vcov(
-                X,
-                residuals,
-                coords,
-                cutoff,
-                "haversine",
-                "bartlett",
-                bread,
-                time=time_alt,
-                unit=unit,
-                lag_cutoff=1,
-            )
-            np.testing.assert_allclose(V_alt, V_dense, atol=1e-12)
-
-    def test_time_label_normalization_datetime64(self):
-        """datetime64 time labels normalize to dense codes via np.unique."""
-        X, residuals, coords, _, unit, bread, cutoff = self._panel_fixture(n_units=6, T=3, k=2)
-        time_dense = np.tile([0, 1, 2], 6)
-        time_dt = np.tile(
-            np.array(["2024-01-01", "2024-04-01", "2024-08-01"], dtype="datetime64[D]"),
-            6,
-        )
-        V_dense = _compute_conley_vcov(
-            X,
-            residuals,
-            coords,
-            cutoff,
-            "haversine",
-            "bartlett",
-            bread,
-            time=time_dense,
-            unit=unit,
-            lag_cutoff=1,
-        )
-        V_dt = _compute_conley_vcov(
-            X,
-            residuals,
-            coords,
-            cutoff,
-            "haversine",
-            "bartlett",
-            bread,
-            time=time_dt,
-            unit=unit,
-            lag_cutoff=1,
-        )
-        np.testing.assert_allclose(V_dt, V_dense, atol=1e-12)
-
-    def test_serial_kernel_bartlett_hardcoded_even_when_kernel_uniform(self):
-        """conleyreg::time_dist hardcodes Bartlett-style temporal kernel
-        regardless of the user's `kernel` choice. We mirror that asymmetry."""
-        # Two panels: same data, one bartlett spatial, one uniform spatial.
-        # The serial contribution should be IDENTICAL because the temporal
-        # kernel is Bartlett-hardcoded.
-        X, residuals, coords, time, unit, bread, cutoff = self._panel_fixture()
-        V_bartlett_L0 = _compute_conley_vcov(
-            X,
-            residuals,
-            coords,
-            cutoff,
-            "haversine",
-            "bartlett",
-            bread,
-            time=time,
-            unit=unit,
-            lag_cutoff=0,
-        )
-        V_bartlett_L2 = _compute_conley_vcov(
-            X,
-            residuals,
-            coords,
-            cutoff,
-            "haversine",
-            "bartlett",
-            bread,
-            time=time,
-            unit=unit,
-            lag_cutoff=2,
-        )
-        V_uniform_L0 = _compute_conley_vcov(
-            X,
-            residuals,
-            coords,
-            cutoff,
-            "haversine",
-            "uniform",
-            bread,
-            time=time,
-            unit=unit,
-            lag_cutoff=0,
-        )
-        V_uniform_L2 = _compute_conley_vcov(
-            X,
-            residuals,
-            coords,
-            cutoff,
-            "haversine",
-            "uniform",
-            bread,
-            time=time,
-            unit=unit,
-            lag_cutoff=2,
-        )
-        # The serial delta should be the same regardless of spatial kernel.
-        # Convert vcov back to meat: meat = bread @ V @ bread
-        delta_bartlett = bread @ (V_bartlett_L2 - V_bartlett_L0) @ bread
-        delta_uniform = bread @ (V_uniform_L2 - V_uniform_L0) @ bread
-        np.testing.assert_allclose(delta_bartlett, delta_uniform, atol=1e-10)
-
-
-# ---------------------------------------------------------------------------
-# TestConleySparse — sparse k-d-tree fast path (Wave A item #120).
-# ---------------------------------------------------------------------------
-
-
-class TestConleySparse:
-    """Sparse k-d-tree fast path for the spatial Bartlett meat.
-
-    The sparse path is gated by three conditions: total n above the
-    threshold, metric in {"haversine", "euclidean"} (no callable), and
-    kernel == "bartlett". Each of these tests exercises one of the
-    gates plus the bit-identity parity claim vs the dense path.
-    """
-
-    def _euclidean_fixture(self, n=1000, k=3, cutoff=15.0, seed=11):
-        rng = np.random.default_rng(seed)
-        coords = rng.uniform(0.0, 100.0, size=(n, 2))
-        X = np.column_stack([np.ones(n)] + [rng.standard_normal(n) for _ in range(k - 1)])
-        beta = np.linspace(0.5, 2.0, k)
-        y = X @ beta + rng.standard_normal(n) * 0.5
-        coefs, *_ = np.linalg.lstsq(X, y, rcond=None)
-        residuals = y - X @ coefs
-        bread = X.T @ X
-        return X, residuals, coords, bread, cutoff
-
-    def _haversine_fixture(self, n=1000, k=3, cutoff_km=500.0, seed=13):
-        rng = np.random.default_rng(seed)
-        lats = rng.uniform(-30.0, 30.0, size=n)
-        lons = rng.uniform(-100.0, 100.0, size=n)
-        coords = np.column_stack([lats, lons])
-        X = np.column_stack([np.ones(n)] + [rng.standard_normal(n) for _ in range(k - 1)])
-        beta = np.linspace(0.5, 2.0, k)
-        y = X @ beta + rng.standard_normal(n) * 0.5
-        coefs, *_ = np.linalg.lstsq(X, y, rcond=None)
-        residuals = y - X @ coefs
-        bread = X.T @ X
-        return X, residuals, coords, bread, cutoff_km
-
-    def test_sparse_vs_dense_bit_identity_euclidean_cross_sectional(self):
-        """Sparse and dense paths produce the same meat on a 1000-row
-        euclidean+bartlett fixture (atol=1e-10). Headroom over the ~1e-14
-        roundoff in chord-projection (haversine) and matmul ordering."""
-        X, residuals, coords, bread, cutoff = self._euclidean_fixture(n=1000)
-        V_dense = _compute_conley_vcov(
-            X,
-            residuals,
-            coords,
-            cutoff,
-            "euclidean",
-            "bartlett",
-            bread,
-            _conley_sparse=False,
-        )
-        V_sparse = _compute_conley_vcov(
-            X,
-            residuals,
-            coords,
-            cutoff,
-            "euclidean",
-            "bartlett",
-            bread,
-            _conley_sparse=True,
-        )
-        np.testing.assert_allclose(V_sparse, V_dense, atol=1e-10, rtol=1e-10)
-
-    def test_sparse_vs_dense_bit_identity_haversine_cross_sectional(self):
-        """Sparse and dense paths produce the same meat on a 1000-row
-        haversine+bartlett fixture (atol=1e-10). Haversine adds the
-        chord-projection roundoff that the tolerance must absorb."""
-        X, residuals, coords, bread, cutoff = self._haversine_fixture(n=1000)
-        V_dense = _compute_conley_vcov(
-            X,
-            residuals,
-            coords,
-            cutoff,
-            "haversine",
-            "bartlett",
-            bread,
-            _conley_sparse=False,
-        )
-        V_sparse = _compute_conley_vcov(
-            X,
-            residuals,
-            coords,
-            cutoff,
-            "haversine",
-            "bartlett",
-            bread,
-            _conley_sparse=True,
-        )
-        np.testing.assert_allclose(V_sparse, V_dense, atol=1e-10, rtol=1e-10)
-
-    def test_auto_toggle_above_threshold_uses_sparse(self, monkeypatch):
-        """n > _CONLEY_SPARSE_N_THRESHOLD with bartlett + euclidean must
-        auto-route through the sparse helper. Verified by spying on the
-        sparse helper call count."""
-        import diff_diff.conley as conley_module
-
-        X, residuals, coords, bread, cutoff = self._euclidean_fixture(
-            n=_CONLEY_SPARSE_N_THRESHOLD + 1
-        )
-        calls = {"n": 0}
-        orig = conley_module._compute_spatial_bartlett_meat_sparse
-
-        def _spy(*args, **kwargs):
-            calls["n"] += 1
-            return orig(*args, **kwargs)
-
-        monkeypatch.setattr(conley_module, "_compute_spatial_bartlett_meat_sparse", _spy)
-        _compute_conley_vcov(X, residuals, coords, cutoff, "euclidean", "bartlett", bread)
-        assert calls["n"] >= 1, "Sparse helper not called when n > threshold."
-
-    def test_auto_toggle_below_threshold_stays_dense(self, monkeypatch):
-        """n <= _CONLEY_SPARSE_N_THRESHOLD must use the dense path even
-        when other sparse conditions (bartlett + euclidean) are met."""
-        import diff_diff.conley as conley_module
-
-        X, residuals, coords, bread, cutoff = self._euclidean_fixture(n=_CONLEY_SPARSE_N_THRESHOLD)
-        calls = {"n": 0}
-        orig = conley_module._compute_spatial_bartlett_meat_sparse
-
-        def _spy(*args, **kwargs):
-            calls["n"] += 1
-            return orig(*args, **kwargs)
-
-        monkeypatch.setattr(conley_module, "_compute_spatial_bartlett_meat_sparse", _spy)
-        _compute_conley_vcov(X, residuals, coords, cutoff, "euclidean", "bartlett", bread)
-        assert calls["n"] == 0, "Sparse helper called below threshold."
-
-    def test_auto_toggle_callable_metric_stays_dense(self, monkeypatch):
-        """A callable conley_metric forces the dense path even at large n —
-        the kd-tree query needs a vectorizable metric, and callables are
-        not supported via projection."""
-        import diff_diff.conley as conley_module
-
-        X, residuals, coords, bread, cutoff = self._euclidean_fixture(
-            n=_CONLEY_SPARSE_N_THRESHOLD + 100
-        )
-
-        def callable_metric(c1, c2):
-            diff = c1[:, None, :] - c2[None, :, :]
-            return np.sqrt(np.sum(diff * diff, axis=-1))
+        def callable_metric(c1, c2):
+            diff = c1[:, None, :] - c2[None, :, :]
+            return np.sqrt(np.sum(diff * diff, axis=-1))
 
         calls = {"n": 0}
         orig = conley_module._compute_spatial_bartlett_meat_sparse
@@ -3400,28 +2499,6 @@ def _spy(*args, **kwargs):
             conley_module._compute_spatial_bartlett_meat_sparse = orig
         assert calls["n"] == 0, "Sparse helper called when _conley_sparse=False."
 
-    def test_helper_direct_matches_dense_meat_euclidean(self):
-        """Call _compute_spatial_bartlett_meat_sparse directly and compare
-        to the dense matmul S' K S on the same data."""
-        X, residuals, coords, _, cutoff = self._euclidean_fixture(n=500)
-        S = X * residuals[:, None]
-        D = _pairwise_distance_matrix(coords, "euclidean")
-        K = _bartlett_kernel(D / cutoff)
-        meat_dense = S.T @ K @ S
-        meat_sparse = _compute_spatial_bartlett_meat_sparse(S, coords, cutoff, "euclidean")
-        np.testing.assert_allclose(meat_sparse, meat_dense, atol=1e-10, rtol=1e-10)
-
-    def test_helper_direct_matches_dense_meat_haversine(self):
-        """Helper direct match for haversine — exercises the chord-projection
-        + exact-distance refinement path."""
-        X, residuals, coords, _, cutoff = self._haversine_fixture(n=500)
-        S = X * residuals[:, None]
-        D = _pairwise_distance_matrix(coords, "haversine")
-        K = _bartlett_kernel(D / cutoff)
-        meat_dense = S.T @ K @ S
-        meat_sparse = _compute_spatial_bartlett_meat_sparse(S, coords, cutoff, "haversine")
-        np.testing.assert_allclose(meat_sparse, meat_dense, atol=1e-10, rtol=1e-10)
-
     def test_sparse_haversine_cutoff_above_half_earth_circumference(self):
         """Sparse haversine path with conley_cutoff_km > π·R_earth (~20,015 km)
         must include all geometrically eligible pairs. Without the arc-radians
@@ -3539,93 +2616,6 @@ def test_sparse_density_gate_does_not_trigger_below_threshold(self):
                 len(density_warnings) == 0
             ), "Density gate should not have triggered at low density"
 
-    def test_sparse_with_cluster_matches_dense(self):
-        """Sparse + combined cluster kernel matches the dense path bit-for-bit
-        at atol=1e-10 (cross-sectional). Codex CI R8 P1: load-bearing
-        sparse+cluster regression that was missing prior."""
-        rng = np.random.default_rng(seed=181)
-        n = 600
-        coords = rng.uniform(0.0, 100.0, size=(n, 2))
-        cluster_ids = rng.integers(0, 8, size=n)  # 8 clusters
-        X = np.column_stack([np.ones(n), rng.standard_normal(n), rng.standard_normal(n)])
-        y = X @ np.array([1.0, 1.5, -0.5]) + rng.standard_normal(n) * 0.4
-        coefs, *_ = np.linalg.lstsq(X, y, rcond=None)
-        residuals = y - X @ coefs
-        bread = X.T @ X
-        V_dense = _compute_conley_vcov(
-            X,
-            residuals,
-            coords,
-            15.0,
-            "euclidean",
-            "bartlett",
-            bread,
-            cluster_ids=cluster_ids,
-            _conley_sparse=False,
-        )
-        V_sparse = _compute_conley_vcov(
-            X,
-            residuals,
-            coords,
-            15.0,
-            "euclidean",
-            "bartlett",
-            bread,
-            cluster_ids=cluster_ids,
-            _conley_sparse=True,
-        )
-        np.testing.assert_allclose(V_sparse, V_dense, atol=1e-10, rtol=1e-10)
-
-    def test_sparse_with_cluster_panel_matches_dense(self):
-        """Sparse + cluster + panel block-decomposed sandwich matches the
-        dense path on a 3-period panel with time-invariant cluster.
-        Codex CI R8 P1."""
-        rng = np.random.default_rng(seed=191)
-        n_units = 200
-        T = 3
-        unit = np.repeat(np.arange(n_units), T)
-        time = np.tile(np.arange(T), n_units)
-        n = n_units * T
-        # Time-invariant coords + cluster (per-unit)
-        unit_coords = rng.uniform(0.0, 100.0, size=(n_units, 2))
-        coords = unit_coords[unit]
-        cluster_per_unit = rng.integers(0, 5, size=n_units)
-        cluster_ids = cluster_per_unit[unit]
-        X = np.column_stack([np.ones(n), rng.standard_normal(n), rng.standard_normal(n)])
-        y = X @ np.array([1.0, 1.5, -0.3]) + rng.standard_normal(n) * 0.4
-        coefs, *_ = np.linalg.lstsq(X, y, rcond=None)
-        residuals = y - X @ coefs
-        bread = X.T @ X
-        V_dense = _compute_conley_vcov(
-            X,
-            residuals,
-            coords,
-            15.0,
-            "euclidean",
-            "bartlett",
-            bread,
-            time=time,
-            unit=unit,
-            lag_cutoff=1,
-            cluster_ids=cluster_ids,
-            _conley_sparse=False,
-        )
-        V_sparse = _compute_conley_vcov(
-            X,
-            residuals,
-            coords,
-            15.0,
-            "euclidean",
-            "bartlett",
-            bread,
-            time=time,
-            unit=unit,
-            lag_cutoff=1,
-            cluster_ids=cluster_ids,
-            _conley_sparse=True,
-        )
-        np.testing.assert_allclose(V_sparse, V_dense, atol=1e-10, rtol=1e-10)
-
     def test_sparse_density_gate_cluster_aware(self):
         """High global spatial density + many small clusters: within-cluster
         density is low, so the sparse path should NOT spuriously fall back
@@ -3716,111 +2706,6 @@ def test_sparse_haversine_cutoff_at_exactly_half_earth_circumference(self):
         )
         np.testing.assert_allclose(V_sparse, V_dense, atol=1e-10, rtol=1e-10)
 
-    def test_panel_block_decomposed_sparse_matches_dense(self):
-        """Panel block-decomposed sandwich produces the same vcov whether
-        the spatial component is computed dense or sparse. The serial
-        component is always dense regardless of the flag."""
-        X, residuals, coords, _, cutoff = self._euclidean_fixture(n=900, seed=21)
-        # Synthetic 3-period panel with 300 units per period
-        time = np.repeat(np.arange(3), 300)
-        unit = np.tile(np.arange(300), 3)
-        bread = X.T @ X
-        V_dense = _compute_conley_vcov(
-            X,
-            residuals,
-            coords,
-            cutoff,
-            "euclidean",
-            "bartlett",
-            bread,
-            time=time,
-            unit=unit,
-            lag_cutoff=1,
-            _conley_sparse=False,
-        )
-        V_sparse = _compute_conley_vcov(
-            X,
-            residuals,
-            coords,
-            cutoff,
-            "euclidean",
-            "bartlett",
-            bread,
-            time=time,
-            unit=unit,
-            lag_cutoff=1,
-            _conley_sparse=True,
-        )
-        np.testing.assert_allclose(V_sparse, V_dense, atol=1e-10, rtol=1e-10)
-
-
-class TestConleySparseRParityForced:
-    """R conleyreg parity at atol=1e-6 with the sparse path FORCED on the
-    three panel R fixtures (bartlett kernel, haversine metric)."""
-
-    GOLDEN_PATH = "benchmarks/data/r_conleyreg_conley_golden.json"
-    PARITY_TOL = 1e-6
-
-    @pytest.fixture(scope="class")
-    def golden(self):
-        import json
-        from pathlib import Path
-
-        repo_root = Path(__file__).resolve().parent.parent
-        path = repo_root / self.GOLDEN_PATH
-        if not path.exists():
-            pytest.skip(
-                f"Golden JSON not present at {path}; run "
-                "`cd benchmarks/R && Rscript generate_conley_golden.R` to generate."
-            )
-        return json.loads(path.read_text())
-
-    def _check_panel_forced_sparse(self, golden, name):
-        entry = golden[name]
-        # Sparse path requires bartlett kernel; skip if fixture is uniform.
-        if entry["kernel"] != "bartlett":
-            pytest.skip(f"Fixture {name!r} is not bartlett; skipped for sparse parity.")
-        X = np.asarray(entry["x"], dtype=np.float64).reshape(entry["x_shape"])
-        y = np.asarray(entry["y"], dtype=np.float64)
-        coords = np.asarray(entry["coords"], dtype=np.float64).reshape(entry["coords_shape"])
-        vcov_expected = np.asarray(entry["vcov"], dtype=np.float64).reshape(entry["vcov_shape"])
-        unit = np.asarray(entry["unit"])
-        time = np.asarray(entry["time"])
-
-        coefs, *_ = np.linalg.lstsq(X, y, rcond=None)
-        residuals = y - X @ coefs
-        bread = X.T @ X
-        vcov_got = _compute_conley_vcov(
-            X,
-            residuals,
-            coords,
-            entry["cutoff_km"],
-            entry["metric"],
-            entry["kernel"],
-            bread,
-            time=time,
-            unit=unit,
-            lag_cutoff=int(entry["lag_cutoff"]),
-            _conley_sparse=True,
-        )
-        np.testing.assert_allclose(
-            vcov_got, vcov_expected, atol=self.PARITY_TOL, rtol=self.PARITY_TOL
-        )
-
-    def test_sparse_parity_panel_haversine_lag1(self, golden):
-        self._check_panel_forced_sparse(golden, "panel_haversine_lag1")
-
-    def test_sparse_parity_panel_haversine_lag2(self, golden):
-        self._check_panel_forced_sparse(golden, "panel_haversine_lag2")
-
-    def test_sparse_parity_panel_lat_lon_realistic_lag1(self, golden):
-        self._check_panel_forced_sparse(golden, "panel_lat_lon_realistic_lag1")
-
-
-# ---------------------------------------------------------------------------
-# TestConleyCluster — combined spatial + cluster product kernel (Wave A #119).
-# ---------------------------------------------------------------------------
-
 
 class TestConleyCluster:
     """Combined spatial + cluster product kernel: K(d_ij/h) * 1{c_i = c_j}.
diff --git a/tests/test_methodology_conley.py b/tests/test_methodology_conley.py
new file mode 100644
index 00000000..06d57b6e
--- /dev/null
+++ b/tests/test_methodology_conley.py
@@ -0,0 +1,1709 @@
+"""ConleySpatialHAC methodology verification tests.
+
+Targets Conley, T. G. (1999), *GMM Estimation with Cross-Sectional Dependence*,
+Journal of Econometrics 92(1), 1-45. DOI: 10.1016/S0304-4076(98)00084-0.
+
+Secondary sources:
+- Andrews, D. W. K. (1991). Heteroskedasticity and autocorrelation consistent
+  covariance matrix estimation. *Econometrica* 59(3), 817-858.
+- Duesterhoeft, C. (2021). conleyreg: Estimations using Conley Standard Errors.
+  CRAN R package. Parity target for the panel block-decomposed sandwich.
+
+Paper review on file: ``docs/methodology/papers/conley-1999-review.md``.
+
+Equation / section walk-through:
+
+- **Eq. 4.2 cross-sectional sandwich (pairwise-distance specialization)**:
+  `Var(beta) = (X'X)^{-1} (Sum_{i,j} K(d_ij/h) X_i eps_i eps_j X_j') (X'X)^{-1}`.
+  diff-diff implements the real-valued / pairwise-distance form (Conley 1999
+  Equation 4.2 plus the "pairwise products at a given distance" remark on
+  page 19); Equation 3.13 is the lattice-indexed form reserved for grid
+  coordinates. Verified by ``TestConleyEquation42`` (Bartlett kernel form,
+  PSD/symmetric/shape contracts, uniform kernel boundary). The 1-D radial
+  Bartlett implementation matches R ``conleyreg``, Stata ``acreg`` (Colella
+  et al. 2019), and Hsiang (2010) as a practitioner specialization of Eq. 4.2.
+- **Eq. 4.2 limits**: tiny-cutoff reduces to HC0 diagonal `Sum_i X_i eps_i^2
+  X_i'`; huge-cutoff under uniform reduces to rank-1 correlated limit. Verified
+  by ``TestConleyHC0AndRank1Reductions``.
+- **Eq. 3.14**: paper's 2-D separable product window
+  `K(j, k) = (1 - |j|/L_M)(1 - |k|/L_N)` on a lattice is the explicitly
+  PSD-guaranteed form. The 1-D radial Bartlett implementation is a practitioner
+  specialization that is **not formally PSD-guaranteed**. The indefiniteness
+  guard at `< -1e-12` applies to both kernels (REGISTRY L3609); locked by
+  ``TestConleyLibraryExtensions::test_indefiniteness_guard_fires_on_negative_eigenvalue``.
+- **Andrews (1991) HAC truncation**: panel space-time Bartlett kernel
+  `(1 - |t-s|/(L+1))` for `0 < |t-s| <= L`, matching
+  ``conleyreg::time_dist.cpp``. Lag = 0 is excluded to avoid double-counting the
+  diagonal already in the spatial component. Verified by
+  ``TestConleyAndrewsLagTruncation``.
+- **Haversine convention**: Earth radius 6371.01 km (matches
+  ``conleyreg::haversine_dist``); 1 degree latitude = 111.195 km at the
+  equator. Verified by ``TestConleyHaversineConvention``.
+- **Phase 2 panel block-decomposed sandwich**: additive decomposition
+  `XeeX = XeeX_spatial + XeeX_serial` (matches ``conleyreg::time_dist.cpp``
+  at `atol=1e-12`). Verified by ``TestConleyBlockDecomposition``.
+- **Wave A #120 sparse k-d-tree path**: numerical correctness of sparse vs
+  dense (bit-identical to `atol=1e-10`). Verified by
+  ``TestConleySparseDenseEquivalence``. Activation thresholds + density-gate
+  fallback are defensive surface (kept in ``tests/test_conley_vcov.py``).
+- **R ``conleyreg`` v0.1.9 parity at `atol=1e-6`**: 6 fixtures (3 cross-sectional
+  + 3 panel) plus the time-asymmetric kernel literal-matching contract. Verified
+  by ``TestConleyParityR``. Skips when the golden JSON is absent (isolated CI).
+
+Library extensions, deviations from R, and deferrals (3 dedicated classes):
+
+- ``TestConleyLibraryExtensions`` (6 tests): combined spatial + cluster product
+  kernel (Wave A #119) — two limit-fixture anchors (all-unique-clusters and
+  huge-cutoff + uniform), callable conley_metric validation (Wave A #123),
+  sparse k-d-tree activation (Wave A #120), and indefiniteness guard
+  separately on the Bartlett and uniform kernel paths (REGISTRY L3609; both
+  kernels are practitioner specializations not formally PSD-guaranteed).
+  All have no R correspondence.
+- ``TestConleyDeviationsFromR`` (3 tests): 1-D radial Bartlett vs paper's 2-D
+  separable Eq. 3.14, time-label normalization via `np.unique`, independent
+  temporal kernel deferred (R's `kernel` arg controls spatial only).
+- ``TestConleyDeferrals`` (5 tests): fail-closed `NotImplementedError` /
+  `TypeError` contracts for LinearRegression + survey_design, DiD/MPD/TWFE +
+  survey_design, Conley + weights, SyntheticDiD + Conley, wild_bootstrap +
+  Conley.
+
+**Generic-path deferral (open methodological question)**: weighted
+spatial-HAC under probability sampling (LinearRegression / DiD / MPD /
+TWFE generic path + Conley + survey_design) is an open methodological
+question; no canonical extension of Conley (1999) exists for the
+combination. The estimator-specific shipped surfaces are: SpilloverDiD
++ Conley + survey via Wave E.1/E.2/E.3 (PR #468/#474/#482, stratified-
+Conley sandwich on PSU totals); TwoStageDiD + Conley + survey via Wave
+E.3 parity (PR #485).
+
+Companion to ``tests/test_conley_vcov.py`` (defensive surface: input
+validation, NaN/inf guards, dispatch-level validity, estimator-level
+integration smoke tests, set_params atomicity). The methodology file
+consolidates paper-anchored R-parity tests into a checklist; the defensive
+file retains the high-volume input-validation + regression surface.
+"""
+
+import json
+import warnings
+from pathlib import Path
+
+import numpy as np
+import pytest
+
+from diff_diff.conley import (
+    _CONLEY_EARTH_RADIUS_KM,
+    _CONLEY_SPARSE_N_THRESHOLD,
+    _bartlett_kernel,
+    _compute_conley_vcov,
+    _haversine_km,
+    _pairwise_distance_matrix,
+    _serial_bartlett_kernel_matrix,
+    _uniform_kernel,
+)
+from diff_diff.linalg import (
+    LinearRegression,
+    compute_robust_vcov,
+)
+
+# ---------------------------------------------------------------------------
+# Shared fixtures and helpers
+# ---------------------------------------------------------------------------
+
+GOLDEN_PATH = "benchmarks/data/r_conleyreg_conley_golden.json"
+PARITY_TOL = 1e-6  # R conleyreg v0.1.9 parity success criterion (Phase 1+2)
+
+
+def _load_r_conley_goldens():
+    """Load the R conleyreg golden fixtures.
+
+    Returns a dict keyed by fixture name; each value is the per-fixture entry
+    (``x``, ``x_shape``, ``y``, ``coords``, ``coords_shape``, ``cutoff_km``,
+    ``metric``, ``kernel``, ``vcov``, ``vcov_shape``; panel fixtures add
+    ``unit``, ``time``, ``lag_cutoff``).
+    """
+    repo_root = Path(__file__).resolve().parent.parent
+    path = repo_root / GOLDEN_PATH
+    if not path.exists():
+        pytest.skip(
+            f"Golden JSON not present at {path}; run "
+            "`cd benchmarks/R && Rscript generate_conley_golden.R` to "
+            "generate. Requires conleyreg R package + sf/lwgeom + system "
+            "libs gdal/proj/geos/udunits."
+        )
+    return json.loads(path.read_text())
+
+
+def _make_conley_cross_section(*, n=20, seed=42):
+    """20-row synthetic OLS dataset with lat/lon coords.
+
+    Used across `TestConleyEquation42`, `TestConleyHC0AndRank1Reductions`,
+    and `TestConleyHaversineConvention` (Eq. 4.2 cross-sectional surface).
+    """
+    rng = np.random.default_rng(seed=seed)
+    X = np.column_stack([np.ones(n), rng.standard_normal(n)])
+    eps = rng.standard_normal(n)
+    y = X @ np.array([1.0, 2.0]) + eps
+    coefs, *_ = np.linalg.lstsq(X, y, rcond=None)
+    residuals = y - X @ coefs
+    bread = X.T @ X
+    coords = np.column_stack(
+        [
+            rng.uniform(-30, 30, n),  # lat
+            rng.uniform(-100, 100, n),  # lon
+        ]
+    )
+    return X, residuals, coords, bread
+
+
+def _make_conley_panel(*, n_units=5, T=3, k=2, seed=42, cutoff_km=5000):
+    """Balanced panel fixture for the Phase 2 block-decomposed surface.
+
+    Used across `TestConleyAndrewsLagTruncation` and
+    `TestConleyBlockDecomposition`. Mirrors the helper in
+    `tests/test_conley_vcov.py::TestConleyPanelHelper`.
+    """
+    rng = np.random.default_rng(seed)
+    lat_unit = rng.uniform(-30, 30, size=n_units)
+    lon_unit = rng.uniform(-100, 100, size=n_units)
+    unit = np.repeat(np.arange(n_units), T)
+    time = np.tile(np.arange(1, T + 1), n_units)
+    lat = lat_unit[unit]
+    lon = lon_unit[unit]
+    coords = np.column_stack([lat, lon])
+    n = n_units * T
+    X = np.column_stack([np.ones(n)] + [rng.standard_normal(n) for _ in range(k - 1)])
+    beta = np.linspace(0.5, 2.0, k)
+    y = X @ beta + rng.standard_normal(n) * 0.5
+    beta_hat, *_ = np.linalg.lstsq(X, y, rcond=None)
+    residuals = y - X @ beta_hat
+    bread = X.T @ X
+    return X, residuals, coords, time, unit, bread, cutoff_km
+
+
+# ---------------------------------------------------------------------------
+# TestConleyEquation42 — Conley 1999 Eq. 4.2 cross-sectional sandwich
+# ---------------------------------------------------------------------------
+
+
+class TestConleyEquation42:
+    """Conley 1999 Eq. 4.2 cross-sectional Bartlett-kernel sandwich.
+
+    `Var(beta) = (X'X)^{-1} (Sum_{i,j} K(d_ij/h) X_i eps_i eps_j X_j') (X'X)^{-1}`
+
+    Verifies the Bartlett kernel form K(u) = max(0, 1 - |u|), the uniform
+    kernel K(u) = 1{|u| <= 1}, and the sandwich PSD/symmetric/shape contracts
+    on a small synthetic OLS dataset.
+    """
+
+    def test_bartlett_at_zero(self):
+        """K(0) = 1: same-location pair contributes the full HC0 term."""
+        np.testing.assert_allclose(_bartlett_kernel(np.array([0.0])), 1.0)
+
+    def test_bartlett_at_one(self):
+        """K(1) = 0: pair at exactly the bandwidth contributes zero."""
+        np.testing.assert_allclose(_bartlett_kernel(np.array([1.0])), 0.0)
+
+    def test_bartlett_above_one_zero(self):
+        """K(u) = 0 for u > 1 (compact support of the Bartlett kernel)."""
+        u = np.array([1.5, 2.0, 100.0])
+        np.testing.assert_allclose(_bartlett_kernel(u), np.zeros(3))
+
+    def test_bartlett_negative_arg_symmetric(self):
+        """K(-u) = K(u): Bartlett kernel uses |u| (Conley 1999 p. 11)."""
+        np.testing.assert_allclose(
+            _bartlett_kernel(np.array([-0.3])),
+            _bartlett_kernel(np.array([0.3])),
+        )
+
+    def test_bartlett_kernel_finite_and_in_unit_interval(self):
+        """Bartlett-weighted kernel matrix on random pairwise distances is
+        finite, symmetric, and bounded in [0, 1]. We do NOT assert PSD here:
+        the radial 1-D Bartlett on pairwise distance is a practitioner
+        specialization of Conley 1999 and is NOT formally PSD-guaranteed —
+        see REGISTRY ConleySpatialHAC § Note (deviation / source
+        specialization). The runtime path emits a UserWarning if the
+        resulting Conley meat is materially indefinite; that contract is
+        locked in
+        ``TestConleyLibraryExtensions::test_indefiniteness_guard_fires_on_negative_eigenvalue``.
+        """
+        rng = np.random.default_rng(seed=11)
+        n = 25
+        coords = rng.uniform(0, 1, size=(n, 2))
+        diff = coords[:, None, :] - coords[None, :, :]
+        D = np.sqrt((diff * diff).sum(axis=-1))
+        K = _bartlett_kernel(D / 0.3)
+        assert K.shape == (n, n)
+        assert np.all(np.isfinite(K))
+        assert np.all(K >= 0.0)
+        assert np.all(K <= 1.0)
+        np.testing.assert_allclose(K, K.T, atol=1e-15)
+
+    def test_uniform_kernel_at_boundary(self):
+        """Uniform kernel is closed on the right: K(1) = 1, K(1+eps) = 0
+        (Conley 1999 p. 11; spectral window is negative in regions per
+        footnote 11)."""
+        np.testing.assert_allclose(_uniform_kernel(np.array([1.0])), 1.0)
+
+    def test_uniform_kernel_above_one_zero(self):
+        np.testing.assert_allclose(_uniform_kernel(np.array([1.0001, 2.0, 100.0])), np.zeros(3))
+
+    def test_uniform_kernel_at_zero(self):
+        np.testing.assert_allclose(_uniform_kernel(np.array([0.0])), 1.0)
+
+    def test_sandwich_returns_psd_with_bartlett(self):
+        """Eq. 4.2 sandwich is PSD on well-behaved residuals (no
+        indefiniteness warning fires)."""
+        X, residuals, coords, bread = _make_conley_cross_section()
+        vcov = _compute_conley_vcov(
+            X,
+            residuals,
+            coords,
+            cutoff=2000.0,
+            metric="haversine",
+            kernel="bartlett",
+            bread_matrix=bread,
+        )
+        eigvals = np.linalg.eigvalsh(0.5 * (vcov + vcov.T))
+        assert eigvals.min() > -1e-10
+
+    def test_sandwich_symmetric_vcov(self):
+        """Eq. 4.2 sandwich is symmetric to machine precision."""
+        X, residuals, coords, bread = _make_conley_cross_section()
+        vcov = _compute_conley_vcov(
+            X,
+            residuals,
+            coords,
+            cutoff=2000.0,
+            metric="haversine",
+            kernel="bartlett",
+            bread_matrix=bread,
+        )
+        np.testing.assert_allclose(vcov, vcov.T, atol=1e-10)
+
+    def test_sandwich_shape_matches_bread(self):
+        """Eq. 4.2 sandwich has shape (k, k) matching the bread matrix."""
+        X, residuals, coords, bread = _make_conley_cross_section()
+        vcov = _compute_conley_vcov(
+            X,
+            residuals,
+            coords,
+            cutoff=1500.0,
+            metric="haversine",
+            kernel="bartlett",
+            bread_matrix=bread,
+        )
+        k = X.shape[1]
+        assert vcov.shape == (k, k)
+
+
+# ---------------------------------------------------------------------------
+# TestConleyHC0AndRank1Reductions — Eq. 4.2 limits
+# ---------------------------------------------------------------------------
+
+
+class TestConleyHC0AndRank1Reductions:
+    """Conley 1999 Eq. 4.2 limits.
+
+    - Tiny-cutoff: K is ~0 off-diagonal and 1 on-diagonal, so the meat
+      reduces to HC0 `Sum_i X_i eps_i^2 X_i'`.
+    - Huge-cutoff under uniform: K -> ones(n, n) and meat = (X*eps)' ones
+      (X*eps) which is rank-1 (NOT HC0). This is the all-correlated limit.
+    - Diagonal contribution: K(0) = 1 so the diagonal always contributes the
+      full HC0 term exactly.
+    """
+
+    def test_tiny_cutoff_bartlett_yields_HC0_meat(self):
+        """Bartlett + bandwidth << min pairwise distance reduces to HC0."""
+        rng = np.random.default_rng(seed=3)
+        n = 15
+        coords = np.column_stack([np.arange(n) * 100.0, np.arange(n) * 100.0])
+        X = np.column_stack([np.ones(n), rng.standard_normal(n)])
+        eps = rng.standard_normal(n)
+        bread = X.T @ X
+
+        meat_hc0 = X.T @ (X * (eps**2)[:, None])
+        bread_inv = np.linalg.solve(bread, np.eye(2))
+        vcov_hc0 = bread_inv @ meat_hc0 @ bread_inv
+
+        vcov_conley = _compute_conley_vcov(
+            X,
+            eps,
+            coords,
+            cutoff=1.0,  # << minimum pairwise distance
+            metric="euclidean",
+            kernel="bartlett",
+            bread_matrix=bread,
+        )
+        np.testing.assert_allclose(vcov_conley, vcov_hc0, atol=1e-12)
+
+    def test_tiny_cutoff_uniform_yields_HC0_meat(self):
+        """Uniform kernel + tiny cutoff also reduces to HC0."""
+        rng = np.random.default_rng(seed=5)
+        n = 12
+        coords = np.column_stack([np.arange(n) * 100.0, np.arange(n) * 100.0])
+        X = np.column_stack([np.ones(n), rng.standard_normal(n)])
+        eps = rng.standard_normal(n)
+        bread = X.T @ X
+
+        meat_hc0 = X.T @ (X * (eps**2)[:, None])
+        bread_inv = np.linalg.solve(bread, np.eye(2))
+        vcov_hc0 = bread_inv @ meat_hc0 @ bread_inv
+
+        vcov_conley = _compute_conley_vcov(
+            X,
+            eps,
+            coords,
+            cutoff=0.5,
+            metric="euclidean",
+            kernel="uniform",
+            bread_matrix=bread,
+        )
+        np.testing.assert_allclose(vcov_conley, vcov_hc0, atol=1e-12)
+
+    def test_huge_cutoff_uniform_NOT_HC0(self):
+        """Huge cutoff under uniform: K -> ones(n, n), meat is rank-1
+        all-correlated limit — NOT HC0."""
+        X, residuals, coords, bread = _make_conley_cross_section()
+        vcov_conley = _compute_conley_vcov(
+            X,
+            residuals,
+            coords,
+            cutoff=1e9,
+            metric="euclidean",
+            kernel="uniform",
+            bread_matrix=bread,
+        )
+        meat_hc0 = X.T @ (X * (residuals**2)[:, None])
+        bread_inv = np.linalg.solve(bread, np.eye(X.shape[1]))
+        vcov_hc0 = bread_inv @ meat_hc0 @ bread_inv
+        assert not np.allclose(vcov_conley, vcov_hc0, atol=1e-6)
+
+    def test_diagonal_of_meat_equals_HC0_contribution(self):
+        """K(0/h) = 1 always, so the diagonal contribution to the meat is
+        exactly the HC0 term `Sum_i X_i eps_i^2 X_i'` regardless of cutoff
+        choice. Extracted from `test_conley_vcov.py::TestConleyReductionsAddendum`."""
+        rng = np.random.default_rng(seed=9)
+        n = 20
+        coords = rng.uniform(0, 1000, size=(n, 2))
+        X = np.column_stack([np.ones(n), rng.standard_normal(n)])
+        eps = rng.standard_normal(n)
+        D = _pairwise_distance_matrix(coords, "euclidean")
+        cutoff = float(D[D > 0].min() * 0.001)
+        S = X * eps[:, None]
+        meat_full = S.T @ _bartlett_kernel(D / cutoff) @ S
+        meat_hc0 = X.T @ (X * (eps**2)[:, None])
+        np.testing.assert_allclose(meat_full, meat_hc0, atol=1e-12)
+
+
+# ---------------------------------------------------------------------------
+# TestConleyAndrewsLagTruncation — Andrews (1991) HAC truncation
+# ---------------------------------------------------------------------------
+
+
+class TestConleyAndrewsLagTruncation:
+    """Andrews (1991) HAC lag truncation for the panel space-time path.
+
+    The within-unit serial Bartlett kernel is `(1 - |t-s|/(L+1))` for
+    `0 < |t-s| <= L`, zero otherwise, matching ``conleyreg::time_dist.cpp``.
+    `lag = 0` is excluded to avoid double-counting the diagonal already in
+    the within-period spatial component.
+    """
+
+    def test_serial_bartlett_kernel_matrix_basic(self):
+        """Hand-computed Bartlett HAC kernel for t=[0,1,2,3], L=2.
+
+        K[i,j] = (1 - |i-j|/(L+1)) for 0 < |i-j| <= L, else 0. With L=2
+        and (L+1)=3: lag 1 -> 2/3, lag 2 -> 1/3, lag 3 -> 0 (out of band).
+        """
+        K = _serial_bartlett_kernel_matrix(np.array([0, 1, 2, 3]), L=2)
+        expected = np.array(
+            [
+                [0.0, 2.0 / 3.0, 1.0 / 3.0, 0.0],
+                [2.0 / 3.0, 0.0, 2.0 / 3.0, 1.0 / 3.0],
+                [1.0 / 3.0, 2.0 / 3.0, 0.0, 2.0 / 3.0],
+                [0.0, 1.0 / 3.0, 2.0 / 3.0, 0.0],
+            ]
+        )
+        np.testing.assert_allclose(K, expected, atol=1e-15)
+
+    def test_serial_bartlett_kernel_matrix_l_one(self):
+        """L=1: only adjacent lags survive, with weight (1 - 1/2) = 0.5."""
+        K = _serial_bartlett_kernel_matrix(np.array([0, 1, 2]), L=1)
+        expected = np.array(
+            [
+                [0.0, 0.5, 0.0],
+                [0.5, 0.0, 0.5],
+                [0.0, 0.5, 0.0],
+            ]
+        )
+        np.testing.assert_allclose(K, expected, atol=1e-15)
+
+    def test_serial_bartlett_kernel_matrix_l_zero_returns_zero(self):
+        """L=0 is degenerate-but-callable: every off-diagonal lag fails the
+        `lag <= 0` band test (since `lag != 0`), so K is the zero matrix.
+        Callers guard externally — ``conley.py`` skips the serial loop when
+        L == 0."""
+        K = _serial_bartlett_kernel_matrix(np.array([0, 1, 2]), L=0)
+        np.testing.assert_array_equal(K, np.zeros((3, 3)))
+
+    def test_serial_bartlett_kernel_matrix_single_element(self):
+        """Single-element input yields a 1x1 zero matrix (no off-diagonal
+        lags exist)."""
+        K = _serial_bartlett_kernel_matrix(np.array([7]), L=2)
+        np.testing.assert_array_equal(K, np.zeros((1, 1)))
+
+    def test_serial_bartlett_kernel_matrix_int_input_bit_equal_to_float(self):
+        """Contract test: int64 and float64 inputs must yield bit-equal
+        matrices. The helper does ``astype(np.float64, copy=False)`` and one
+        of the call sites (`conley.py` panel-block branch) passes the result
+        of array slicing on int time codes."""
+        K_int = _serial_bartlett_kernel_matrix(np.array([0, 1, 2, 3], dtype=np.int64), L=2)
+        K_float = _serial_bartlett_kernel_matrix(
+            np.array([0.0, 1.0, 2.0, 3.0], dtype=np.float64), L=2
+        )
+        np.testing.assert_array_equal(K_int, K_float)
+
+    def test_T_eq_1_equals_cross_sectional(self):
+        """Block-decomposed form with T=1 (single time period) equals the
+        Phase 1 cross-sectional form on the same data (serial component
+        contributes nothing at T=1)."""
+        X, residuals, coords, _, _, bread, cutoff = _make_conley_panel(n_units=8, T=1, k=2)
+        unit_single = np.arange(X.shape[0])
+        time_single = np.ones(X.shape[0], dtype=int)
+        V_cs = _compute_conley_vcov(X, residuals, coords, cutoff, "haversine", "bartlett", bread)
+        V_panel = _compute_conley_vcov(
+            X,
+            residuals,
+            coords,
+            cutoff,
+            "haversine",
+            "bartlett",
+            bread,
+            time=time_single,
+            unit=unit_single,
+            lag_cutoff=2,
+        )
+        np.testing.assert_allclose(V_panel, V_cs, atol=1e-12)
+
+    def test_lag_cutoff_zero_drops_serial(self):
+        """L=0 means the serial component contributes nothing; only the
+        within-period spatial sandwich applies."""
+        X, residuals, coords, time, unit, bread, cutoff = _make_conley_panel()
+        V0 = _compute_conley_vcov(
+            X,
+            residuals,
+            coords,
+            cutoff,
+            "haversine",
+            "bartlett",
+            bread,
+            time=time,
+            unit=unit,
+            lag_cutoff=0,
+        )
+        S = X * residuals[:, None]
+        meat_spatial = np.zeros((X.shape[1], X.shape[1]))
+        for t_val in np.unique(time):
+            mask = time == t_val
+            D_t = _pairwise_distance_matrix(coords[mask], "haversine")
+            K_t = _bartlett_kernel(D_t / cutoff)
+            meat_spatial += S[mask].T @ K_t @ S[mask]
+        V_expected = np.linalg.solve(bread, meat_spatial)
+        V_expected = np.linalg.solve(bread, V_expected.T).T
+        np.testing.assert_allclose(V0, V_expected, atol=1e-12)
+
+    def test_lag_cutoff_positive_adds_serial(self):
+        """L > 0 strictly increases the meat by a positive serial
+        contribution (within-unit cross-time pairs)."""
+        X, residuals, coords, time, unit, bread, cutoff = _make_conley_panel()
+        V0 = _compute_conley_vcov(
+            X,
+            residuals,
+            coords,
+            cutoff,
+            "haversine",
+            "bartlett",
+            bread,
+            time=time,
+            unit=unit,
+            lag_cutoff=0,
+        )
+        V1 = _compute_conley_vcov(
+            X,
+            residuals,
+            coords,
+            cutoff,
+            "haversine",
+            "bartlett",
+            bread,
+            time=time,
+            unit=unit,
+            lag_cutoff=1,
+        )
+        assert not np.allclose(V0, V1, atol=1e-8), (
+            "lag_cutoff=1 must differ from lag_cutoff=0 with a serial " "component"
+        )
+
+
+# ---------------------------------------------------------------------------
+# TestConleyHaversineConvention — Earth radius + degree conversion
+# ---------------------------------------------------------------------------
+
+
+class TestConleyHaversineConvention:
+    """Great-circle distance convention.
+
+    diff-diff uses Earth's mean radius 6371.01 km, matching
+    ``conleyreg::haversine_dist`` (per Duesterhoeft 2021 CRAN v0.1.9 src).
+    At the equator, 1 degree longitude = 2*pi*R/360 = 111.195 km.
+    """
+
+    def test_haversine_known_pair_one_degree_equator(self):
+        """1 degree longitude at the equator = 2*pi*R/360 ~ 111.195 km."""
+        d = _haversine_km(np.array(0.0), np.array(0.0), np.array(0.0), np.array(1.0))
+        expected = 2 * np.pi * _CONLEY_EARTH_RADIUS_KM / 360.0
+        np.testing.assert_allclose(d, expected, atol=1e-9)
+
+    def test_haversine_zero_self_distance(self):
+        """d(p, p) = 0 (load-bearing for the K(0) = 1 identity in
+        Eq. 4.2's diagonal contribution)."""
+        d = _haversine_km(
+            np.array(45.0),
+            np.array(-122.0),
+            np.array(45.0),
+            np.array(-122.0),
+        )
+        np.testing.assert_allclose(d, 0.0, atol=1e-12)
+
+    def test_haversine_symmetric(self):
+        """d(a, b) = d(b, a) (load-bearing for symmetric kernel matrix)."""
+        d_ab = _haversine_km(
+            np.array(40.7),
+            np.array(-74.0),
+            np.array(34.0),
+            np.array(-118.2),
+        )
+        d_ba = _haversine_km(
+            np.array(34.0),
+            np.array(-118.2),
+            np.array(40.7),
+            np.array(-74.0),
+        )
+        np.testing.assert_allclose(d_ab, d_ba, atol=1e-12)
+
+    def test_haversine_pole_to_equator(self):
+        """North pole to equator at any longitude = pi/2 * R ~ 10007.5 km."""
+        d = _haversine_km(
+            np.array(90.0),
+            np.array(0.0),
+            np.array(0.0),
+            np.array(0.0),
+        )
+        expected = np.pi * _CONLEY_EARTH_RADIUS_KM / 2.0
+        np.testing.assert_allclose(d, expected, atol=1e-9)
+
+    def test_haversine_broadcasting_pairwise(self):
+        """Broadcasting (n, 1) vs (1, n) yields the n*n distance matrix."""
+        coords = np.array([[0.0, 0.0], [0.0, 1.0], [0.0, 2.0]])
+        lats = coords[:, 0]
+        lons = coords[:, 1]
+        D = _haversine_km(lats[:, None], lons[:, None], lats[None, :], lons[None, :])
+        assert D.shape == (3, 3)
+        np.testing.assert_allclose(np.diag(D), 0.0, atol=1e-12)
+        np.testing.assert_allclose(D[0, 2], 2.0 * D[0, 1], rtol=1e-10)
+
+    def test_pairwise_distance_haversine_dispatch(self):
+        """Pairwise distance helper dispatches haversine correctly."""
+        coords = np.array([[0.0, 0.0], [0.0, 1.0], [10.0, 0.0]])
+        D = _pairwise_distance_matrix(coords, "haversine")
+        assert D.shape == (3, 3)
+        np.testing.assert_allclose(D, D.T, atol=1e-12)
+        np.testing.assert_allclose(np.diag(D), 0.0, atol=1e-12)
+
+    def test_pairwise_distance_euclidean_matches_pdist(self):
+        """Euclidean path matches scipy.spatial.distance squareform exactly
+        (sanity check for the alternative metric path)."""
+        from scipy.spatial.distance import pdist, squareform
+
+        rng = np.random.default_rng(seed=7)
+        coords = rng.uniform(0, 100, size=(15, 2))
+        D = _pairwise_distance_matrix(coords, "euclidean")
+        D_scipy = squareform(pdist(coords, metric="euclidean"))
+        np.testing.assert_allclose(D, D_scipy, atol=1e-12)
+
+
+# ---------------------------------------------------------------------------
+# TestConleyBlockDecomposition — Phase 2 panel block-decomposed sandwich
+# ---------------------------------------------------------------------------
+
+
+class TestConleyBlockDecomposition:
+    """Phase 2 panel block-decomposed sandwich `XeeX = XeeX_spatial +
+    XeeX_serial` matching ``conleyreg::time_dist.cpp``.
+
+    Internal cross-check at `atol=1e-12` against a hand-coded reference
+    that re-implements the block decomposition directly from the C++ form.
+    """
+
+    def test_panel_matches_block_decomposed_reference(self):
+        """Direct verification that ``_compute_conley_vcov`` matches a
+        hand-coded block decomposition (matches R ``conleyreg::time_dist.cpp``)
+        at machine precision across L in (0, 1, 2)."""
+        X, residuals, coords, time, unit, bread, cutoff = _make_conley_panel(seed=314)
+        bread_inv = np.linalg.inv(bread)
+        S = X * residuals[:, None]
+        for L in (0, 1, 2):
+            meat = np.zeros((X.shape[1], X.shape[1]))
+            for t_val in np.unique(time):
+                mask = time == t_val
+                D_t = _pairwise_distance_matrix(coords[mask], "haversine")
+                K_t = _bartlett_kernel(D_t / cutoff)
+                meat += S[mask].T @ K_t @ S[mask]
+            if L > 0:
+                for u_val in np.unique(unit):
+                    mask = unit == u_val
+                    S_u = S[mask]
+                    t_u = time[mask].astype(np.float64)
+                    lag = np.abs(t_u[:, None] - t_u[None, :])
+                    K_u = ((lag <= L) & (lag != 0)).astype(np.float64) * (1.0 - lag / (L + 1.0))
+                    meat += S_u.T @ K_u @ S_u
+            V_ref = bread_inv @ meat @ bread_inv
+
+            V_helper = _compute_conley_vcov(
+                X,
+                residuals,
+                coords,
+                cutoff,
+                "haversine",
+                "bartlett",
+                bread,
+                time=time,
+                unit=unit,
+                lag_cutoff=L,
+            )
+            np.testing.assert_allclose(V_helper, V_ref, atol=1e-12)
+
+    def test_cluster_time_invariance_validated_on_panel_path(self):
+        """REGISTRY § ConleySpatialHAC § Combined spatial + cluster product
+        kernel: on the panel path the validator REQUIRES cluster_ids to be
+        constant within each unit across periods, raising ValueError with
+        the violating unit named otherwise. Anchors the contract that
+        guarantees the within-unit serial sandwich's cluster mask is
+        trivially all-ones, simplifying the math to bare serial Bartlett HAC
+        weighted by the spatial mask only."""
+        rng = np.random.default_rng(seed=17)
+        n_units = 4
+        T = 3
+        n = n_units * T
+        unit = np.repeat(np.arange(n_units), T)
+        time = np.tile(np.arange(T), n_units)
+        X = np.column_stack([np.ones(n), rng.standard_normal(n)])
+        eps = rng.standard_normal(n)
+        coords = rng.uniform(0, 1, size=(n, 2))
+        # Time-varying cluster for unit 0
+        cluster = unit.copy()
+        cluster[0] = 99  # unit 0 at time 0 has cluster 99 instead of 0
+        with pytest.raises(ValueError, match=r"(?i)cluster|time-invariant|unit"):
+            compute_robust_vcov(
+                X,
+                eps,
+                vcov_type="conley",
+                conley_coords=coords,
+                conley_cutoff_km=10.0,
+                conley_metric="euclidean",
+                conley_time=time,
+                conley_unit=unit,
+                conley_lag_cutoff=1,
+                cluster_ids=cluster,
+            )
+
+
+# ---------------------------------------------------------------------------
+# TestConleySparseDenseEquivalence — Wave A #120 numerical correctness
+# ---------------------------------------------------------------------------
+
+
+class TestConleySparseDenseEquivalence:
+    """Sparse k-d-tree path produces meat bit-identical (to atol=1e-10) to
+    the dense path.
+
+    The methodology anchor for Wave A #120 is *numerical correctness*: the
+    sparse path is an optimization, not a methodology change, so its output
+    must agree with the dense reference at full precision (modulo
+    chord-projection roundoff on haversine, which is absorbed by atol=1e-10).
+    Sparse-path activation thresholds, density-gate fallback, callable-
+    metric / uniform-kernel fallbacks are defensive surface and remain in
+    ``tests/test_conley_vcov.py``.
+    """
+
+    def _euclidean_fixture(self, *, n=1000, k=3, cutoff=15.0, seed=11):
+        rng = np.random.default_rng(seed)
+        coords = rng.uniform(0.0, 100.0, size=(n, 2))
+        X = np.column_stack([np.ones(n)] + [rng.standard_normal(n) for _ in range(k - 1)])
+        beta = np.linspace(0.5, 2.0, k)
+        y = X @ beta + rng.standard_normal(n) * 0.5
+        coefs, *_ = np.linalg.lstsq(X, y, rcond=None)
+        residuals = y - X @ coefs
+        bread = X.T @ X
+        return X, residuals, coords, bread, cutoff
+
+    def _haversine_fixture(self, *, n=1000, k=3, cutoff_km=500.0, seed=13):
+        rng = np.random.default_rng(seed)
+        lats = rng.uniform(-30.0, 30.0, size=n)
+        lons = rng.uniform(-100.0, 100.0, size=n)
+        coords = np.column_stack([lats, lons])
+        X = np.column_stack([np.ones(n)] + [rng.standard_normal(n) for _ in range(k - 1)])
+        beta = np.linspace(0.5, 2.0, k)
+        y = X @ beta + rng.standard_normal(n) * 0.5
+        coefs, *_ = np.linalg.lstsq(X, y, rcond=None)
+        residuals = y - X @ coefs
+        bread = X.T @ X
+        return X, residuals, coords, bread, cutoff_km
+
+    @pytest.mark.slow
+    def test_sparse_vs_dense_euclidean_cross_sectional(self):
+        """Sparse and dense paths produce the same meat on a 1000-row
+        euclidean+bartlett fixture at `atol=1e-10`."""
+        X, residuals, coords, bread, cutoff = self._euclidean_fixture(n=1000)
+        V_dense = _compute_conley_vcov(
+            X,
+            residuals,
+            coords,
+            cutoff,
+            "euclidean",
+            "bartlett",
+            bread,
+            _conley_sparse=False,
+        )
+        V_sparse = _compute_conley_vcov(
+            X,
+            residuals,
+            coords,
+            cutoff,
+            "euclidean",
+            "bartlett",
+            bread,
+            _conley_sparse=True,
+        )
+        np.testing.assert_allclose(V_sparse, V_dense, atol=1e-10, rtol=1e-10)
+
+    @pytest.mark.slow
+    def test_sparse_vs_dense_haversine_cross_sectional(self):
+        """Sparse and dense paths produce the same meat on a 1000-row
+        haversine+bartlett fixture at `atol=1e-10`. Haversine adds the
+        chord-projection roundoff that the tolerance must absorb."""
+        X, residuals, coords, bread, cutoff = self._haversine_fixture(n=1000)
+        V_dense = _compute_conley_vcov(
+            X,
+            residuals,
+            coords,
+            cutoff,
+            "haversine",
+            "bartlett",
+            bread,
+            _conley_sparse=False,
+        )
+        V_sparse = _compute_conley_vcov(
+            X,
+            residuals,
+            coords,
+            cutoff,
+            "haversine",
+            "bartlett",
+            bread,
+            _conley_sparse=True,
+        )
+        np.testing.assert_allclose(V_sparse, V_dense, atol=1e-10, rtol=1e-10)
+
+    @pytest.mark.slow
+    def test_sparse_vs_dense_panel_block_decomposed(self):
+        """Panel block-decomposed sandwich produces the same vcov whether the
+        spatial component is computed dense or sparse. The serial component
+        is always dense regardless of the flag. Locks Wave A #120 numerical
+        correctness on the Phase 2 path."""
+        X, residuals, coords, bread, cutoff = self._euclidean_fixture(n=900, seed=21)
+        time = np.repeat(np.arange(3), 300)
+        unit = np.tile(np.arange(300), 3)
+        V_dense = _compute_conley_vcov(
+            X,
+            residuals,
+            coords,
+            cutoff,
+            "euclidean",
+            "bartlett",
+            bread,
+            time=time,
+            unit=unit,
+            lag_cutoff=1,
+            _conley_sparse=False,
+        )
+        V_sparse = _compute_conley_vcov(
+            X,
+            residuals,
+            coords,
+            cutoff,
+            "euclidean",
+            "bartlett",
+            bread,
+            time=time,
+            unit=unit,
+            lag_cutoff=1,
+            _conley_sparse=True,
+        )
+        np.testing.assert_allclose(V_sparse, V_dense, atol=1e-10, rtol=1e-10)
+
+    @pytest.mark.slow
+    def test_sparse_with_cluster_cross_sectional(self):
+        """Sparse + combined cluster kernel (Wave A #119) matches the dense
+        path bit-for-bit at `atol=1e-10` on a cross-sectional fixture.
+        Load-bearing for the sparse + Wave A #119 product-kernel composition."""
+        rng = np.random.default_rng(seed=181)
+        n = 600
+        coords = rng.uniform(0.0, 100.0, size=(n, 2))
+        cluster_ids = rng.integers(0, 8, size=n)
+        X = np.column_stack([np.ones(n), rng.standard_normal(n), rng.standard_normal(n)])
+        y = X @ np.array([1.0, 1.5, -0.5]) + rng.standard_normal(n) * 0.4
+        coefs, *_ = np.linalg.lstsq(X, y, rcond=None)
+        residuals = y - X @ coefs
+        bread = X.T @ X
+        V_dense = _compute_conley_vcov(
+            X,
+            residuals,
+            coords,
+            15.0,
+            "euclidean",
+            "bartlett",
+            bread,
+            cluster_ids=cluster_ids,
+            _conley_sparse=False,
+        )
+        V_sparse = _compute_conley_vcov(
+            X,
+            residuals,
+            coords,
+            15.0,
+            "euclidean",
+            "bartlett",
+            bread,
+            cluster_ids=cluster_ids,
+            _conley_sparse=True,
+        )
+        np.testing.assert_allclose(V_sparse, V_dense, atol=1e-10, rtol=1e-10)
+
+    @pytest.mark.slow
+    def test_sparse_with_cluster_panel(self):
+        """Sparse + Wave A #119 cluster product kernel + Phase 2 panel
+        block-decomposed sandwich matches the dense path on a 3-period
+        panel with time-invariant cluster. Cluster-time-invariance contract
+        guarantees the within-unit serial mask is trivially all-ones."""
+        rng = np.random.default_rng(seed=191)
+        n_units = 200
+        T = 3
+        unit = np.repeat(np.arange(n_units), T)
+        time = np.tile(np.arange(T), n_units)
+        n = n_units * T
+        unit_coords = rng.uniform(0.0, 100.0, size=(n_units, 2))
+        coords = unit_coords[unit]
+        cluster_per_unit = rng.integers(0, 5, size=n_units)
+        cluster_ids = cluster_per_unit[unit]
+        X = np.column_stack([np.ones(n), rng.standard_normal(n), rng.standard_normal(n)])
+        y = X @ np.array([1.0, 1.5, -0.3]) + rng.standard_normal(n) * 0.4
+        coefs, *_ = np.linalg.lstsq(X, y, rcond=None)
+        residuals = y - X @ coefs
+        bread = X.T @ X
+        V_dense = _compute_conley_vcov(
+            X,
+            residuals,
+            coords,
+            15.0,
+            "euclidean",
+            "bartlett",
+            bread,
+            time=time,
+            unit=unit,
+            lag_cutoff=1,
+            cluster_ids=cluster_ids,
+            _conley_sparse=False,
+        )
+        V_sparse = _compute_conley_vcov(
+            X,
+            residuals,
+            coords,
+            15.0,
+            "euclidean",
+            "bartlett",
+            bread,
+            time=time,
+            unit=unit,
+            lag_cutoff=1,
+            cluster_ids=cluster_ids,
+            _conley_sparse=True,
+        )
+        np.testing.assert_allclose(V_sparse, V_dense, atol=1e-10, rtol=1e-10)
+
+
+# ---------------------------------------------------------------------------
+# TestConleyParityR — R conleyreg v0.1.9 parity at atol=1e-6
+# ---------------------------------------------------------------------------
+
+
+class TestConleyParityR:
+    """R ``conleyreg`` v0.1.9 parity (Duesterhoeft 2021 CRAN).
+
+    Goldens at ``benchmarks/data/r_conleyreg_conley_golden.json``; generator
+    at ``benchmarks/R/generate_conley_golden.R``. 6 fixtures:
+
+    - 3 cross-sectional (Phase 1): ``small_haversine``, ``dense_haversine``,
+      ``lat_lon_realistic``.
+    - 3 panel (Phase 2 block-decomposed with `lag_cutoff > 0`):
+      ``panel_haversine_lag1``, ``panel_haversine_lag2``,
+      ``panel_lat_lon_realistic_lag1``.
+
+    Parity tolerance: `atol=1e-6` and `rtol=1e-6`. Earth radius 6371.01 km
+    matches ``conleyreg::haversine_dist``. The sparse k-d-tree path forced
+    on the same goldens is exercised at ``test_sparse_forced_matches_r_*``
+    to verify the optimization path also matches R parity.
+
+    Also verifies the time-asymmetric kernel parity contract: R
+    ``conleyreg``'s `kernel` argument controls the spatial component only;
+    the temporal kernel is unconditionally Bartlett (`(1 - |t-s|/(L+1))`)
+    regardless of `kernel` (see ``conleyreg::time_dist.cpp`` for the
+    hardcoding). diff-diff matches this asymmetry exactly. This is a
+    parity contract, not a deviation — the independent-temporal-kernel
+    follow-up is documented as a deferral in
+    ``TestConleyDeviationsFromR::test_independent_temporal_kernel_deferred``.
+    """
+
+    def _check_cross_section_fixture(self, golden, name):
+        entry = golden[name]
+        X = np.asarray(entry["x"], dtype=np.float64).reshape(entry["x_shape"])
+        y = np.asarray(entry["y"], dtype=np.float64)
+        coords = np.asarray(entry["coords"], dtype=np.float64).reshape(entry["coords_shape"])
+        vcov_expected = np.asarray(entry["vcov"], dtype=np.float64).reshape(entry["vcov_shape"])
+
+        coefs, *_ = np.linalg.lstsq(X, y, rcond=None)
+        residuals = y - X @ coefs
+        vcov_got = compute_robust_vcov(
+            X,
+            residuals,
+            vcov_type="conley",
+            conley_coords=coords,
+            conley_cutoff_km=entry["cutoff_km"],
+            conley_metric=entry["metric"],
+            conley_kernel=entry["kernel"],
+        )
+        np.testing.assert_allclose(vcov_got, vcov_expected, atol=PARITY_TOL, rtol=PARITY_TOL)
+
+    def _check_panel_fixture(self, golden, name):
+        entry = golden[name]
+        X = np.asarray(entry["x"], dtype=np.float64).reshape(entry["x_shape"])
+        y = np.asarray(entry["y"], dtype=np.float64)
+        coords = np.asarray(entry["coords"], dtype=np.float64).reshape(entry["coords_shape"])
+        vcov_expected = np.asarray(entry["vcov"], dtype=np.float64).reshape(entry["vcov_shape"])
+        unit = np.asarray(entry["unit"])
+        time = np.asarray(entry["time"])
+
+        coefs, *_ = np.linalg.lstsq(X, y, rcond=None)
+        residuals = y - X @ coefs
+        vcov_got = compute_robust_vcov(
+            X,
+            residuals,
+            vcov_type="conley",
+            conley_coords=coords,
+            conley_cutoff_km=entry["cutoff_km"],
+            conley_metric=entry["metric"],
+            conley_kernel=entry["kernel"],
+            conley_time=time,
+            conley_unit=unit,
+            conley_lag_cutoff=int(entry["lag_cutoff"]),
+        )
+        np.testing.assert_allclose(vcov_got, vcov_expected, atol=PARITY_TOL, rtol=PARITY_TOL)
+
+    def test_parity_small_haversine(self):
+        self._check_cross_section_fixture(_load_r_conley_goldens(), "small_haversine")
+
+    def test_parity_dense_haversine(self):
+        self._check_cross_section_fixture(_load_r_conley_goldens(), "dense_haversine")
+
+    def test_parity_lat_lon_realistic(self):
+        self._check_cross_section_fixture(_load_r_conley_goldens(), "lat_lon_realistic")
+
+    def test_parity_panel_haversine_lag1(self):
+        self._check_panel_fixture(_load_r_conley_goldens(), "panel_haversine_lag1")
+
+    def test_parity_panel_haversine_lag2(self):
+        self._check_panel_fixture(_load_r_conley_goldens(), "panel_haversine_lag2")
+
+    def test_parity_panel_lat_lon_realistic_lag1(self):
+        self._check_panel_fixture(_load_r_conley_goldens(), "panel_lat_lon_realistic_lag1")
+
+    def test_sparse_forced_matches_r_cross_sectional(self):
+        """Sparse path forced on the R conleyreg cross-sectional fixtures
+        also passes parity at `atol=1e-6`. Verifies Wave A #120 sparse
+        optimization preserves R parity, not just dense-vs-sparse internal
+        equivalence."""
+        golden = _load_r_conley_goldens()
+        for name in ("small_haversine", "dense_haversine", "lat_lon_realistic"):
+            entry = golden[name]
+            if entry["kernel"] != "bartlett":
+                continue  # sparse only supports bartlett
+            X = np.asarray(entry["x"], dtype=np.float64).reshape(entry["x_shape"])
+            y = np.asarray(entry["y"], dtype=np.float64)
+            coords = np.asarray(entry["coords"], dtype=np.float64).reshape(entry["coords_shape"])
+            vcov_expected = np.asarray(entry["vcov"], dtype=np.float64).reshape(entry["vcov_shape"])
+            coefs, *_ = np.linalg.lstsq(X, y, rcond=None)
+            residuals = y - X @ coefs
+            bread = X.T @ X
+            vcov_got = _compute_conley_vcov(
+                X,
+                residuals,
+                coords,
+                entry["cutoff_km"],
+                entry["metric"],
+                entry["kernel"],
+                bread,
+                _conley_sparse=True,
+            )
+            np.testing.assert_allclose(vcov_got, vcov_expected, atol=PARITY_TOL, rtol=PARITY_TOL)
+
+    def test_sparse_forced_matches_r_panel(self):
+        """Sparse path forced on the R conleyreg panel fixtures (Phase 2
+        block-decomposed sandwich) also passes parity at `atol=1e-6`.
+        Verifies Wave A #120 sparse optimization preserves R parity on the
+        spatial component of the panel path; the serial component is always
+        dense regardless of the flag."""
+        golden = _load_r_conley_goldens()
+        for name in (
+            "panel_haversine_lag1",
+            "panel_haversine_lag2",
+            "panel_lat_lon_realistic_lag1",
+        ):
+            entry = golden[name]
+            if entry["kernel"] != "bartlett":
+                continue  # sparse only supports bartlett
+            X = np.asarray(entry["x"], dtype=np.float64).reshape(entry["x_shape"])
+            y = np.asarray(entry["y"], dtype=np.float64)
+            coords = np.asarray(entry["coords"], dtype=np.float64).reshape(entry["coords_shape"])
+            vcov_expected = np.asarray(entry["vcov"], dtype=np.float64).reshape(entry["vcov_shape"])
+            unit = np.asarray(entry["unit"])
+            time = np.asarray(entry["time"])
+            coefs, *_ = np.linalg.lstsq(X, y, rcond=None)
+            residuals = y - X @ coefs
+            bread = X.T @ X
+            vcov_got = _compute_conley_vcov(
+                X,
+                residuals,
+                coords,
+                entry["cutoff_km"],
+                entry["metric"],
+                entry["kernel"],
+                bread,
+                time=time,
+                unit=unit,
+                lag_cutoff=int(entry["lag_cutoff"]),
+                _conley_sparse=True,
+            )
+            np.testing.assert_allclose(vcov_got, vcov_expected, atol=PARITY_TOL, rtol=PARITY_TOL)
+
+    def test_time_asymmetric_kernel_matches_r_literal(self):
+        """R ``conleyreg::time_dist.cpp`` hardcodes the Bartlett temporal
+        kernel regardless of the user's `kernel` argument — the spatial
+        kernel choice does NOT propagate to the serial component. diff-diff
+        matches this asymmetry exactly: the serial-contribution delta
+        between two `lag_cutoff` values is IDENTICAL whether the spatial
+        kernel is bartlett or uniform.
+
+        Moved from `TestConleyDeviationsFromR` to `TestConleyParityR` per
+        plan review MEDIUM #1 — this is a parity contract with R
+        ``conleyreg`` literal, not a deviation. The independent-temporal-
+        kernel API extension would be a deviation and is documented as a
+        deferral in
+        ``TestConleyDeviationsFromR::test_independent_temporal_kernel_deferred``.
+        """
+        X, residuals, coords, time, unit, bread, cutoff = _make_conley_panel()
+        V_bartlett_L0 = _compute_conley_vcov(
+            X,
+            residuals,
+            coords,
+            cutoff,
+            "haversine",
+            "bartlett",
+            bread,
+            time=time,
+            unit=unit,
+            lag_cutoff=0,
+        )
+        V_bartlett_L2 = _compute_conley_vcov(
+            X,
+            residuals,
+            coords,
+            cutoff,
+            "haversine",
+            "bartlett",
+            bread,
+            time=time,
+            unit=unit,
+            lag_cutoff=2,
+        )
+        V_uniform_L0 = _compute_conley_vcov(
+            X,
+            residuals,
+            coords,
+            cutoff,
+            "haversine",
+            "uniform",
+            bread,
+            time=time,
+            unit=unit,
+            lag_cutoff=0,
+        )
+        V_uniform_L2 = _compute_conley_vcov(
+            X,
+            residuals,
+            coords,
+            cutoff,
+            "haversine",
+            "uniform",
+            bread,
+            time=time,
+            unit=unit,
+            lag_cutoff=2,
+        )
+        # The serial-contribution delta must match across spatial kernels
+        # because the temporal kernel is hardcoded Bartlett in both
+        # diff-diff and R conleyreg.
+        delta_bartlett = bread @ (V_bartlett_L2 - V_bartlett_L0) @ bread
+        delta_uniform = bread @ (V_uniform_L2 - V_uniform_L0) @ bread
+        np.testing.assert_allclose(delta_bartlett, delta_uniform, atol=1e-10)
+
+
+# ---------------------------------------------------------------------------
+# TestConleyLibraryExtensions — Wave A library extensions (no R correspondence)
+# ---------------------------------------------------------------------------
+
+
+class TestConleyLibraryExtensions:
+    """Wave A library extensions documented as library-side additions with
+    no paper or R correspondence (6 tests).
+
+    [x] Combined spatial + cluster product kernel (Wave A #119) — Anchor 1:
+        all-unique-clusters → HC0 diagonal
+    [x] Combined spatial + cluster product kernel (Wave A #119) — Anchor 2:
+        huge-cutoff + uniform → within-cluster CR1 plain
+    [x] Callable conley_metric validation (Wave A #123)
+    [x] Sparse k-d-tree fast path (Wave A #120) — activation thresholds and
+        boundary behavior
+    [x] Indefiniteness guard on Bartlett kernel
+    [x] Indefiniteness guard on uniform kernel
+    """
+
+    def test_combined_spatial_cluster_kernel_wave_a_119_all_unique(self):
+        """Wave A #119 library extension: combined kernel
+        `K_total[i, j] = K_space(d_ij/h) * 1{cluster_i = cluster_j}`.
+
+        Anchor 1 (REGISTRY § Combined spatial + cluster product kernel):
+        all-unique-clusters reduction — cluster mask is identity, meat
+        reduces to HC0 diagonal `Sum_i X_i eps_i^2 X_i'`.
+        """
+        rng = np.random.default_rng(seed=21)
+        n = 30
+        X = np.column_stack([np.ones(n), rng.standard_normal(n)])
+        eps = rng.standard_normal(n)
+        coords = rng.uniform(0, 1, size=(n, 2))
+        cluster_unique = np.arange(n)
+        bread = X.T @ X
+        bread_inv = np.linalg.solve(bread, np.eye(2))
+        meat_hc0 = X.T @ (X * (eps**2)[:, None])
+        vcov_hc0 = bread_inv @ meat_hc0 @ bread_inv
+        vcov_combined = compute_robust_vcov(
+            X,
+            eps,
+            vcov_type="conley",
+            conley_coords=coords,
+            conley_cutoff_km=10.0,  # large; full kernel reduces to ones
+            conley_metric="euclidean",
+            conley_kernel="uniform",
+            cluster_ids=cluster_unique,
+        )
+        np.testing.assert_allclose(vcov_combined, vcov_hc0, atol=1e-12)
+
+    def test_combined_spatial_cluster_kernel_wave_a_119_huge_cutoff(self):
+        """Wave A #119 library extension: combined kernel
+        `K_total[i, j] = K_space(d_ij/h) * 1{cluster_i = cluster_j}`.
+
+        Anchor 2 (REGISTRY § Combined spatial + cluster product kernel):
+        huge-cutoff under uniform — `K_space = 1` on every pair, so the
+        meat reduces to the pure within-cluster sum
+        `Sum_g X_g' eps_g eps_g' X_g` (CR1 without the Liang-Zeger
+        small-sample correction). This exact identity holds only for
+        `conley_kernel="uniform"` (K_uniform(u) = 1 for |u| <= 1); the
+        Bartlett kernel gives `K(u) = 1 - |u| < 1` for u > 0, so the
+        Bartlett huge-cutoff limit is asymptotic, not exact at any finite
+        cutoff. The anchor uses uniform for an exact identity check.
+        """
+        rng = np.random.default_rng(seed=22)
+        n = 24
+        n_clusters = 6
+        X = np.column_stack([np.ones(n), rng.standard_normal(n), rng.standard_normal(n)])
+        eps = rng.standard_normal(n)
+        coords = rng.uniform(0, 1, size=(n, 2))
+        cluster_ids = np.repeat(np.arange(n_clusters), n // n_clusters)
+        bread = X.T @ X
+        bread_inv = np.linalg.solve(bread, np.eye(X.shape[1]))
+        # Pure within-cluster CR1 sandwich (no Liang-Zeger correction)
+        S = X * eps[:, None]
+        meat_cr1_plain = np.zeros((X.shape[1], X.shape[1]))
+        for c in np.unique(cluster_ids):
+            S_c = S[cluster_ids == c]
+            S_c_sum = S_c.sum(axis=0)
+            meat_cr1_plain += np.outer(S_c_sum, S_c_sum)
+        vcov_cr1_plain = bread_inv @ meat_cr1_plain @ bread_inv
+        # Huge-cutoff under uniform reduces combined kernel to within-cluster
+        vcov_combined = compute_robust_vcov(
+            X,
+            eps,
+            vcov_type="conley",
+            conley_coords=coords,
+            conley_cutoff_km=1e9,  # >> all pairwise distances
+            conley_metric="euclidean",
+            conley_kernel="uniform",
+            cluster_ids=cluster_ids,
+        )
+        np.testing.assert_allclose(vcov_combined, vcov_cr1_plain, atol=1e-12)
+
+    def test_callable_metric_validation_wave_a_123(self):
+        """Wave A #123 library extension: callable conley_metric is
+        validated at the boundary with 6 invariants — float64-castable,
+        `(n, n)` shape, finite, non-negative, symmetric to `atol=1e-10`,
+        and zero diagonal `|d(i, i)| <= 1e-10` (so K(0) = 1 reduces to
+        the HC0 diagonal contribution). Asymmetric callable raises
+        `ValueError`."""
+        rng = np.random.default_rng(seed=31)
+        n = 10
+        X = np.column_stack([np.ones(n), rng.standard_normal(n)])
+        eps = rng.standard_normal(n)
+        coords = rng.uniform(0, 1, size=(n, 2))
+
+        def asymmetric_metric(c1, c2):
+            d = np.zeros((len(c1), len(c2)))
+            for i in range(len(c1)):
+                for j in range(len(c2)):
+                    d[i, j] = float(i + 2 * j)  # asymmetric: d[i, j] != d[j, i]
+            return d
+
+        with pytest.raises(ValueError, match=r"(?i)symmetric|asymmetric"):
+            compute_robust_vcov(
+                X,
+                eps,
+                vcov_type="conley",
+                conley_coords=coords,
+                conley_cutoff_km=10.0,
+                conley_metric=asymmetric_metric,
+            )
+
+    def test_sparse_kd_tree_activation_wave_a_120(self):
+        """Wave A #120 library extension: sparse k-d-tree path auto-activates
+        when `n > _CONLEY_SPARSE_N_THRESHOLD` AND `conley_kernel == "bartlett"`
+        AND `conley_metric` is haversine or euclidean (no callable). At
+        `n + 1` above threshold the sparse helper is invoked; at threshold
+        the dense path stays."""
+        import diff_diff.conley as conley_module
+
+        rng = np.random.default_rng(seed=41)
+        n = _CONLEY_SPARSE_N_THRESHOLD + 1
+        coords = rng.uniform(0.0, 100.0, size=(n, 2))
+        X = np.column_stack([np.ones(n), rng.standard_normal(n)])
+        eps = rng.standard_normal(n)
+        bread = X.T @ X
+
+        calls = {"n": 0}
+        orig = conley_module._compute_spatial_bartlett_meat_sparse
+
+        def _spy(*args, **kwargs):
+            calls["n"] += 1
+            return orig(*args, **kwargs)
+
+        conley_module._compute_spatial_bartlett_meat_sparse = _spy
+        try:
+            _compute_conley_vcov(X, eps, coords, 15.0, "euclidean", "bartlett", bread)
+        finally:
+            conley_module._compute_spatial_bartlett_meat_sparse = orig
+        assert calls["n"] >= 1, "Sparse helper not called when n > threshold."
+
+    def test_indefiniteness_guard_fires_on_negative_eigenvalue(self):
+        """REGISTRY L3609: the 1-D radial Bartlett implementation and the
+        uniform kernel are practitioner specializations of Conley 1999
+        (not the paper's explicitly PSD-guaranteed 2-D separable Eq. 3.14
+        form). Library hardening: the indefiniteness guard emits a
+        `UserWarning` when any meat eigenvalue is below `-1e-12`. Forces
+        the indefinite path by monkey-patching the kernel to return a
+        sign-pattern that produces an indefinite meat, then asserts the
+        warning surfaces with the kernel name."""
+        from diff_diff import conley as conley_mod
+
+        rng = np.random.default_rng(seed=11)
+        n = 6
+        coords = rng.uniform(0, 1, size=(n, 2))
+        X = np.column_stack([np.ones(n), rng.standard_normal(n)])
+        eps = np.ones(n)
+        bread = X.T @ X
+
+        original = conley_mod._bartlett_kernel
+
+        def _indefinite(u: np.ndarray) -> np.ndarray:
+            base = np.eye(u.shape[0])
+            for i in range(u.shape[0]):
+                for j in range(u.shape[0]):
+                    if i != j:
+                        base[i, j] = -10.0
+            return base
+
+        try:
+            conley_mod._bartlett_kernel = _indefinite
+            with warnings.catch_warnings(record=True) as w:
+                warnings.simplefilter("always")
+                conley_mod._compute_conley_vcov(
+                    X,
+                    eps,
+                    coords,
+                    cutoff=10.0,
+                    metric="euclidean",
+                    kernel="bartlett",
+                    bread_matrix=bread,
+                )
+            psd_warnings = [
+                msg
+                for msg in w
+                if issubclass(msg.category, UserWarning)
+                and "bartlett" in str(msg.message)
+                and "negative eigenvalue" in str(msg.message)
+            ]
+            assert len(psd_warnings) >= 1, (
+                f"Expected UserWarning naming kernel='bartlett' and "
+                f"'negative eigenvalue'; got {[str(m.message) for m in w]}"
+            )
+        finally:
+            conley_mod._bartlett_kernel = original
+
+    def test_indefiniteness_guard_fires_on_negative_eigenvalue_uniform(self):
+        """REGISTRY L3609 + the production indefiniteness guard at
+        ``diff_diff/linalg.py:L3024-L3030`` covers BOTH kernels (uniform is
+        documented as not PSD-guaranteed per Conley 1999 footnote 11 —
+        spectral window is negative in regions). Library hardening: the
+        guard emits a ``UserWarning`` naming the kernel when any meat
+        eigenvalue is below ``-1e-12``. Forces the indefinite path by
+        monkey-patching ``_uniform_kernel`` to a sign-pattern that produces
+        an indefinite meat, then asserts the warning surfaces with the
+        uniform-kernel name (regression for the uniform-kernel branch of
+        the same contract)."""
+        from diff_diff import conley as conley_mod
+
+        rng = np.random.default_rng(seed=12)
+        n = 6
+        coords = rng.uniform(0, 1, size=(n, 2))
+        X = np.column_stack([np.ones(n), rng.standard_normal(n)])
+        eps = np.ones(n)
+        bread = X.T @ X
+
+        original = conley_mod._uniform_kernel
+
+        def _indefinite(u: np.ndarray) -> np.ndarray:
+            base = np.eye(u.shape[0])
+            for i in range(u.shape[0]):
+                for j in range(u.shape[0]):
+                    if i != j:
+                        base[i, j] = -10.0
+            return base
+
+        try:
+            conley_mod._uniform_kernel = _indefinite
+            with warnings.catch_warnings(record=True) as w:
+                warnings.simplefilter("always")
+                conley_mod._compute_conley_vcov(
+                    X,
+                    eps,
+                    coords,
+                    cutoff=10.0,
+                    metric="euclidean",
+                    kernel="uniform",
+                    bread_matrix=bread,
+                )
+            psd_warnings = [
+                msg
+                for msg in w
+                if issubclass(msg.category, UserWarning)
+                and "uniform" in str(msg.message)
+                and "negative eigenvalue" in str(msg.message)
+            ]
+            assert len(psd_warnings) >= 1, (
+                f"Expected UserWarning naming kernel='uniform' and "
+                f"'negative eigenvalue'; got {[str(m.message) for m in w]}"
+            )
+        finally:
+            conley_mod._uniform_kernel = original
+
+
+# ---------------------------------------------------------------------------
+# TestConleyDeviationsFromR — Deviations from R conleyreg literal
+# ---------------------------------------------------------------------------
+
+
+class TestConleyDeviationsFromR:
+    """Deviations from the R ``conleyreg`` literal that are intentional,
+    documented in REGISTRY, and exercised here.
+
+    [x] 1-D radial Bartlett vs paper's 2-D separable Eq. 3.14 (matches
+        ``conleyreg``, Stata ``acreg``, Hsiang 2010 — not formally
+        PSD-guaranteed)
+    [x] Time-label normalization via `np.unique(return_inverse=True)`
+        (diff-diff normalizes to dense panel-period codes; R uses raw
+        time values literally)
+    [x] Independent temporal kernel deferred (R `kernel` arg controls
+        spatial only; the literal-matching parity assertion is in
+        ``TestConleyParityR::test_time_asymmetric_kernel_matches_r_literal``)
+    """
+
+    def test_1d_radial_bartlett_vs_2d_separable_eq314(self):
+        """The 1-D radial Bartlett K(u) = max(0, 1 - |u|) on pairwise
+        distance is a practitioner specialization that diff-diff,
+        ``conleyreg``, Stata ``acreg``, and Hsiang (2010) all share. It
+        differs from Conley 1999's explicitly PSD-guaranteed 2-D separable
+        Eq. 3.14 form `K(j, k) = (1 - |j|/L_M)(1 - |k|/L_N)` indexed on a
+        lattice. The 1-D form is NOT formally PSD-guaranteed, hence the
+        indefiniteness guard at `< -1e-12` documented in REGISTRY L3609.
+
+        This test pins the 1-D form numerically (NOT the 2-D form). The
+        2-D form would require a lattice-indexed implementation which
+        diff-diff does not provide."""
+        u = np.array([0.0, 0.25, 0.5, 0.75, 1.0, 1.25])
+        expected = np.array([1.0, 0.75, 0.5, 0.25, 0.0, 0.0])
+        np.testing.assert_allclose(_bartlett_kernel(u), expected, atol=1e-15)
+
+    def test_time_label_normalization_deviation_from_r(self):
+        """REGISTRY L3592-3603: diff-diff normalizes time labels via
+        `np.unique(return_inverse=True)` to dense panel-period codes
+        before lag computation; R ``conleyreg`` uses raw time values
+        literally (`time_dist.cpp`'s `t_diff = abs(times - times[i])`).
+
+        On dense integer labels (the parity-test convention), the two
+        paths produce bit-identical results. For non-dense encodings —
+        e.g., `time = 202012, 202101` for monthly panels — the raw R
+        difference is 89, so a `lag_cutoff=1` request silently drops
+        valid lag-1 serial pairs in R. diff-diff is the more robust
+        default; this test verifies bit-equality between dense codes
+        and YYYYMM labels for the same underlying panel."""
+        X, residuals, coords, _, unit, bread, cutoff = _make_conley_panel(n_units=8, T=3, k=2)
+        time_dense = np.tile([1, 2, 3], 8)
+        time_yyyymm = np.tile([202011, 202012, 202101], 8)
+        V_dense = _compute_conley_vcov(
+            X,
+            residuals,
+            coords,
+            cutoff,
+            "haversine",
+            "bartlett",
+            bread,
+            time=time_dense,
+            unit=unit,
+            lag_cutoff=1,
+        )
+        V_yyyymm = _compute_conley_vcov(
+            X,
+            residuals,
+            coords,
+            cutoff,
+            "haversine",
+            "bartlett",
+            bread,
+            time=time_yyyymm,
+            unit=unit,
+            lag_cutoff=1,
+        )
+        np.testing.assert_allclose(V_yyyymm, V_dense, atol=1e-12)
+
+    def test_independent_temporal_kernel_deferred(self):
+        """REGISTRY L3584-3591: R ``conleyreg``'s `kernel` argument controls
+        the spatial component only; the temporal kernel is unconditionally
+        Bartlett. diff-diff matches this asymmetry exactly for parity (see
+        ``TestConleyParityR::test_time_asymmetric_kernel_matches_r_literal``
+        for the literal-matching assertion).
+
+        Independent temporal-kernel choice would be a follow-up API
+        extension if user demand emerges. This test pins the current
+        contract by asserting `conley_kernel` is the only kernel knob in
+        the supported API surface — there is no `conley_temporal_kernel`
+        kwarg yet."""
+        import inspect
+
+        sig = inspect.signature(compute_robust_vcov)
+        param_names = set(sig.parameters.keys())
+        assert "conley_kernel" in param_names
+        # No independent temporal kernel parameter is supported.
+        assert "conley_temporal_kernel" not in param_names
+        assert "conley_time_kernel" not in param_names
+
+
+# ---------------------------------------------------------------------------
+# TestConleyDeferrals — Fail-closed NotImplementedError / TypeError contracts
+# ---------------------------------------------------------------------------
+
+
+class TestConleyDeferrals:
+    """Fail-closed contracts asserting unsupported combinations raise
+    explicit `NotImplementedError` / `TypeError`.
+
+    [x] `LinearRegression(vcov_type="conley", survey_design=...)` (generic-
+        path deferral; open methodological question)
+    [x] DiD / MPD / TWFE + Conley + survey_design (estimator-level)
+    [x] Conley + weights (cross-sectional sampling weights without survey
+        design)
+    [x] SyntheticDiD + Conley (uses bootstrap/jackknife/placebo variance)
+    [x] wild_bootstrap + Conley (separate inference path)
+
+    Generic-path deferral rationale: weighted spatial-HAC under probability
+    sampling is an open methodological question; no canonical extension of
+    Conley (1999) exists for the combination. Shipped estimator-specific
+    surfaces: SpilloverDiD + Conley + survey via Wave E.1/E.2/E.3 (PR
+    #468/#474/#482); TwoStageDiD + Conley + survey via Wave E.3 parity
+    (PR #485).
+    """
+
+    def test_linear_regression_survey_design_not_implemented(self):
+        """`LinearRegression(vcov_type='conley', survey_design=...)` raises
+        `NotImplementedError` at the LinearRegression entry (linalg.py
+        front-door rejection). Uses `inspect.signature` to confirm the
+        `survey_design` kwarg exists (so the rejection path is reachable)
+        before asserting the error."""
+        import inspect
+
+        import pandas as pd
+
+        from diff_diff.survey import SurveyDesign
+
+        sig = inspect.signature(LinearRegression.__init__)
+        assert "survey_design" in sig.parameters, (
+            "LinearRegression must accept survey_design kwarg for the "
+            "generic-path Conley+survey rejection contract to be reachable"
+        )
+
+        rng = np.random.default_rng(seed=51)
+        n = 30
+        y = rng.standard_normal(n)
+        X = np.column_stack([np.ones(n), rng.standard_normal(n)])
+        coords = np.column_stack([rng.uniform(-30, 30, n), rng.uniform(-100, 100, n)])
+        df = pd.DataFrame(
+            {
+                "psu": np.repeat(np.arange(5), n // 5),
+                "weights": np.ones(n),
+            }
+        )
+        sd = SurveyDesign(psu="psu", weights="weights").resolve(df)
+        with pytest.raises(NotImplementedError, match=r"(?i)conley.*survey|open methodological"):
+            LinearRegression(
+                vcov_type="conley",
+                conley_coords=coords,
+                conley_cutoff_km=1000.0,
+                survey_design=sd,
+            ).fit(X, y)
+
+    def test_did_mpd_twfe_survey_design_not_implemented(self):
+        """REGISTRY L3760: DiD / MPD / TWFE + Conley + survey_design raises
+        `NotImplementedError` at the estimator level. Tests on DiD as the
+        primary surface; MPD/TWFE share the same fail-closed contract."""
+        import pandas as pd
+
+        from diff_diff.estimators import DifferenceInDifferences
+        from diff_diff.survey import SurveyDesign
+
+        rng = np.random.default_rng(seed=61)
+        n_units = 6
+        T = 2
+        n = n_units * T
+        df = pd.DataFrame(
+            {
+                "y": rng.standard_normal(n),
+                "treated": np.tile([0, 0, 1, 1, 0, 1], 2),
+                "post": np.repeat([0, 1], n_units),
+                "unit": np.tile(np.arange(n_units), T),
+                "psu": np.tile(np.repeat([0, 1, 2], 2), T),
+                "lat": np.tile(rng.uniform(-30, 30, n_units), T),
+                "lon": np.tile(rng.uniform(-100, 100, n_units), T),
+                "weights": np.ones(n),
+            }
+        )
+        sd = SurveyDesign(psu="psu", weights="weights")
+        with pytest.raises(NotImplementedError, match=r"(?i)conley.*survey|open methodological"):
+            DifferenceInDifferences(
+                vcov_type="conley",
+                conley_coords=("lat", "lon"),
+                conley_cutoff_km=1000.0,
+                conley_lag_cutoff=0,
+            ).fit(
+                df,
+                formula="y ~ treated * post",
+                unit="unit",
+                survey_design=sd,
+            )
+
+    def test_conley_plus_weights_not_implemented(self):
+        """REGISTRY L3762: cross-sectional sampling weights without survey
+        design also raise `NotImplementedError`. Same open methodological
+        question — weighted spatial-HAC under probability sampling has no
+        canonical extension of Conley (1999)."""
+        rng = np.random.default_rng(seed=71)
+        n = 20
+        X = np.column_stack([np.ones(n), rng.standard_normal(n)])
+        eps = rng.standard_normal(n)
+        coords = rng.uniform(0, 1, size=(n, 2))
+        weights = rng.uniform(0.5, 1.5, size=n)
+        with pytest.raises(NotImplementedError, match=r"(?i)conley.*weight|open methodological"):
+            compute_robust_vcov(
+                X,
+                eps,
+                vcov_type="conley",
+                conley_coords=coords,
+                conley_cutoff_km=10.0,
+                conley_metric="euclidean",
+                weights=weights,
+            )
+
+    def test_synthetic_did_conley_typeerror(self):
+        """`SyntheticDiD(vcov_type='conley')` raises `TypeError` —
+        SyntheticDiD uses bootstrap / jackknife / placebo variance, not
+        the analytical sandwich. Tracked in TODO.md."""
+        from diff_diff.synthetic_did import SyntheticDiD
+
+        with pytest.raises(TypeError, match=r"(?i)conley|vcov_type"):
+            SyntheticDiD(vcov_type="conley")
+
+    def test_wild_bootstrap_conley_not_implemented(self):
+        """REGISTRY L3547: wild bootstrap is a separate inference path that
+        does not consume the analytical Conley sandwich. Combining
+        `vcov_type='conley'` with `inference='wild_bootstrap'` raises
+        `NotImplementedError`."""
+        import pandas as pd
+
+        from diff_diff.estimators import DifferenceInDifferences
+
+        rng = np.random.default_rng(seed=81)
+        n_units = 6
+        T = 2
+        n = n_units * T
+        df = pd.DataFrame(
+            {
+                "y": rng.standard_normal(n),
+                "treated": np.tile([0, 0, 1, 1, 0, 1], 2),
+                "post": np.repeat([0, 1], n_units),
+                "unit": np.tile(np.arange(n_units), T),
+                "lat": np.tile(rng.uniform(-30, 30, n_units), T),
+                "lon": np.tile(rng.uniform(-100, 100, n_units), T),
+            }
+        )
+        with pytest.raises(NotImplementedError, match=r"(?i)wild.bootstrap|conley.*inference"):
+            DifferenceInDifferences(
+                vcov_type="conley",
+                conley_coords=("lat", "lon"),
+                conley_cutoff_km=1000.0,
+                conley_lag_cutoff=0,
+                inference="wild_bootstrap",
+            ).fit(df, formula="y ~ treated * post", unit="unit")