Skip to content

Commit 5c836c5

Browse files
igerberclaude
andcommitted
Drop overstated Bartlett PSD claim; apply indefiniteness guard to both kernels
Address CI Codex review of PR #411 (P1 + P3): P1 — The implementation called the radial 1-D Bartlett kernel "PSD-guaranteed" by citing Conley 1999 Eq 3.14 + Andrews 1991. Conley's explicit PSD Bartlett formula (Eq 3.14, page 12) is the 2-D separable product window `(1 - |j|/L_M)(1 - |k|/L_N)` indexed on a lattice; the 1-D radial form on pairwise distance that diff-diff and R `conleyreg` implement is a practitioner specialization (Hsiang 2010, Colella et al. 2019) that is not explicitly written in the paper and is therefore not formally PSD-guaranteed. Reframe the kernel docstring around the practitioner-specialization framing, drop the PSD-guaranteed claim, and lift the meat-eigenvalue PSD guard out of the uniform-only branch so it fires for both supported kernels. The warning message now names the active kernel and explicitly states neither radial form is formally PSD. New regression test test_indefinite_meat_warning_fires_for_bartlett locks the lifted guard by patching `_bartlett_kernel` to return an aggressively indefinite matrix and asserting the warning surfaces with the kernel name. P3 — Stale wording cleanup: - conley.py:128-132 missing-coords error message pointed at TwoWayFixedEffects(conley_coords=...) even though TWFE rejects Conley in Phase 1; redirect to LinearRegression / compute_robust_vcov. - TestConleyEstimatorIntegration class docstring claimed panel estimators accept Conley and print a label; rewrote to describe fit-time panel rejection. Doc surfaces (REGISTRY ConleySpatialHAC kernel section, llms-full.txt kernels block, CHANGELOG `conley_kernel` description, conley-1999- review.md PSD-failure note) updated to reflect the both-kernels guard and the radial-specialization framing. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
1 parent 4444375 commit 5c836c5

6 files changed

Lines changed: 130 additions & 31 deletions

File tree

CHANGELOG.md

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
88
## [Unreleased]
99

1010
### Added
11-
- **Conley (1999) spatial-HAC standard errors via `vcov_type="conley"`** on cross-sectional `LinearRegression` / `compute_robust_vcov` (Phase 1 of the spillover-conley initiative). Keyword arguments: `conley_coords` (n × 2 array of lat/lon or projected coords), `conley_cutoff_km=<float>` (positive finite bandwidth in km for haversine, or coord units for euclidean — REQUIRED, no default per the no-silent-failures contract), `conley_metric="haversine"|"euclidean"|callable` (default `"haversine"`; great-circle uses Earth's mean radius 6371.01 km matching R `conleyreg`), `conley_kernel="bartlett"|"uniform"` (default `"bartlett"` is PSD-guaranteed; `"uniform"` emits `UserWarning` if the meat has a materially negative eigenvalue per Conley 1999 footnote 11). Variance estimator `Var̂(β) = (X'X)^{-1} · ( Σ_{i,j} K(d_ij/h) · X_i ε_i ε_j X_j' ) · (X'X)^{-1}` (Conley 1999 Eq 4.2). **Panel estimators (`DifferenceInDifferences`, `TwoWayFixedEffects`, `MultiPeriodDiD`) reject `vcov_type="conley"` at fit-time with `NotImplementedError`** — Phase 1's cross-sectional Conley does not handle the time dimension. Applying it over (unit, time) rows would treat same-unit cross-time pairs as `d_ij = 0 → K = 1`, mishandling the space-time HAC. Practitioners needing Conley with a panel design should pre-collapse to per-unit first-differences and call `compute_robust_vcov` directly on a single-period regression. Phase 2 will add the space-time product kernel (Driscoll-Kraay) for full panel support. `SyntheticDiD(vcov_type="conley")` raises `TypeError` (uses bootstrap variance, not analytical sandwich); `set_params` mirrors the constructor rejection. `vcov_type="conley"` + `cluster_ids=` / `weights=` / `survey_design=` raises `NotImplementedError` (combined product kernel + Bertanha-Imbens 2014 weighted-Conley deferred to follow-up phases). `n > 20_000` emits a `UserWarning` about the dense O(n²) distance-matrix memory; sparse k-d-tree fast path is queued for Phase 2. Helpers live in new module `diff_diff/conley.py` (`_haversine_km`, `_pairwise_distance_matrix`, `_bartlett_kernel`, `_uniform_kernel`, `_validate_conley_kwargs`, `_compute_conley_vcov`); `compute_robust_vcov` in `diff_diff/linalg.py` imports the dispatch helpers. R `conleyreg` parity (Düsterhöft 2021, CRAN v0.1.9) on three benchmark fixtures (`benchmarks/data/r_conleyreg_conley_golden.json`, regenerable via `benchmarks/R/generate_conley_golden.R`); observed max abs diff 5.7e-16. Earth radius 6371.01 km matches `conleyreg::haversine_dist`. Test file `tests/test_conley_vcov.py` skips parity cleanly when the JSON is absent. New REGISTRY section `## ConleySpatialHAC`. Tracked on `BRIEFING.md` as Phase 1 of the 6-phase initiative (Phase 2: space-time product kernel + sparse fast path + panel-estimator support; Phase 3: ring-indicator spillover-aware DiD per Butts 2021; Phase 4a/4b: mechanical extension to IF-aggregation and sandwich-derived estimators; Phase 5: survey design support).
11+
- **Conley (1999) spatial-HAC standard errors via `vcov_type="conley"`** on cross-sectional `LinearRegression` / `compute_robust_vcov` (Phase 1 of the spillover-conley initiative). Keyword arguments: `conley_coords` (n × 2 array of lat/lon or projected coords), `conley_cutoff_km=<float>` (positive finite bandwidth in km for haversine, or coord units for euclidean — REQUIRED, no default per the no-silent-failures contract), `conley_metric="haversine"|"euclidean"|callable` (default `"haversine"`; great-circle uses Earth's mean radius 6371.01 km matching R `conleyreg`), `conley_kernel="bartlett"|"uniform"` (default `"bartlett"` evaluated on pairwise distance `d_ij/h`, matching R `conleyreg`; both kernels emit a `UserWarning` if the resulting meat has a materially negative eigenvalue. Conley 1999's explicit PSD Bartlett formula is the 2-D separable product window on a lattice (Eq 3.14); the 1-D radial pairwise specialization that diff-diff and R `conleyreg` implement is a practitioner convention that is not formally PSD-guaranteed). Variance estimator `Var̂(β) = (X'X)^{-1} · ( Σ_{i,j} K(d_ij/h) · X_i ε_i ε_j X_j' ) · (X'X)^{-1}` (Conley 1999 Eq 4.2). **Panel estimators (`DifferenceInDifferences`, `TwoWayFixedEffects`, `MultiPeriodDiD`) reject `vcov_type="conley"` at fit-time with `NotImplementedError`** — Phase 1's cross-sectional Conley does not handle the time dimension. Applying it over (unit, time) rows would treat same-unit cross-time pairs as `d_ij = 0 → K = 1`, mishandling the space-time HAC. Practitioners needing Conley with a panel design should pre-collapse to per-unit first-differences and call `compute_robust_vcov` directly on a single-period regression. Phase 2 will add the space-time product kernel (Driscoll-Kraay) for full panel support. `SyntheticDiD(vcov_type="conley")` raises `TypeError` (uses bootstrap variance, not analytical sandwich); `set_params` mirrors the constructor rejection. `vcov_type="conley"` + `cluster_ids=` / `weights=` / `survey_design=` raises `NotImplementedError` (combined product kernel + Bertanha-Imbens 2014 weighted-Conley deferred to follow-up phases). `n > 20_000` emits a `UserWarning` about the dense O(n²) distance-matrix memory; sparse k-d-tree fast path is queued for Phase 2. Helpers live in new module `diff_diff/conley.py` (`_haversine_km`, `_pairwise_distance_matrix`, `_bartlett_kernel`, `_uniform_kernel`, `_validate_conley_kwargs`, `_compute_conley_vcov`); `compute_robust_vcov` in `diff_diff/linalg.py` imports the dispatch helpers. R `conleyreg` parity (Düsterhöft 2021, CRAN v0.1.9) on three benchmark fixtures (`benchmarks/data/r_conleyreg_conley_golden.json`, regenerable via `benchmarks/R/generate_conley_golden.R`); observed max abs diff 5.7e-16. Earth radius 6371.01 km matches `conleyreg::haversine_dist`. Test file `tests/test_conley_vcov.py` skips parity cleanly when the JSON is absent. New REGISTRY section `## ConleySpatialHAC`. Tracked on `BRIEFING.md` as Phase 1 of the 6-phase initiative (Phase 2: space-time product kernel + sparse fast path + panel-estimator support; Phase 3: ring-indicator spillover-aware DiD per Butts 2021; Phase 4a/4b: mechanical extension to IF-aggregation and sandwich-derived estimators; Phase 5: survey design support).
1212
- **Tutorial 21: HAD Pre-test Workflow** (`docs/tutorials/21_had_pretest_workflow.ipynb`) — composite pre-test walkthrough for `HeterogeneousAdoptionDiD` building on Tutorial 20's brand-campaign framing. Uses a 60-DMA × 8-week panel close in shape to T20's but with the dose distribution drawn from `Uniform[$0.01K, $50K]` (vs T20's `[$5K, $50K]`); the true support is strictly positive but very near zero, chosen so the QUG step in `did_had_pretest_workflow` fails-to-reject `H0: d_lower = 0` in this finite sample and the verdict text fires the load-bearing "Assumption 7 deferred" pivot for the upgrade-arc narrative. (HAD's `design="auto"` selector — a separate min/median heuristic at `had.py::_detect_design`, NOT the QUG p-value — independently lands on the `continuous_at_zero` identification path with target `WAS` on this panel because `d.min() < 0.01 * median(|d|)`. The QUG test and the design selector are independent rules that point to the same identification path here.) Walks through three surfaces: (a) `did_had_pretest_workflow(aggregate="overall")` on a two-period collapse, where the verdict explicitly flags Step 2 (Assumption 7 pre-trends) as not run because a single pre-period structurally cannot support a pre-trends test, and the structural fields `pretrends_joint` / `homogeneity_joint` are both `None`; (b) `did_had_pretest_workflow(aggregate="event_study")` on the full multi-period panel, where the verdict reads "TWFE admissible under Section 4 assumptions" because all three testable diagnostics (QUG + joint pre-trends Stute over 3 horizons + joint homogeneity Stute over 4 horizons) fail-to-reject — non-rejection evidence under finite-sample power and test specification, not proof that the identifying assumptions hold; and (c) a side panel exercising both `yatchew_hr_test` null modes — `null="linearity"` (default, paper Theorem 7) vs `null="mean_independence"` (Phase 4 R-parity with R `YatchewTest::yatchew_test(order=0)`) — on the within-pre-period first-difference paired with post-period dose, illustrating the stricter null's larger residual variance (`sigma2_lin` 7.01 vs 6.53) and smaller p-value (0.29 vs 0.49). Companion drift-test file `tests/test_t21_had_pretest_workflow_drift.py` (16 tests pinning panel composition, both verdict pivots, structural anchors on both paths, deterministic QUG / Yatchew statistics, bootstrap p-value tolerance bands per `feedback_bootstrap_drift_tests_need_backend_tolerance`, and `HAD(design="auto")` resolution to `continuous_at_zero` on this panel). T20's "Composite pretest workflow" Extensions bullet updated with a forward-pointer to T21. T22 weighted/survey HAD tutorial remains queued as a separate notebook PR.
1313
- **`ChaisemartinDHaultfoeuille.by_path` and `paths_of_interest` now compose with `survey_design`** for analytical Binder TSL SE and replicate-weight bootstrap variance. The `NotImplementedError` gate at `chaisemartin_dhaultfoeuille.py:1233-1239` is replaced by a per-path multiplier-bootstrap-only gate (`survey_design + n_bootstrap > 0` under by_path / paths_of_interest still raises, since the survey-aware perturbation pivot for path-restricted IFs is methodologically underived). Per-path SE routes through the existing `_survey_se_from_group_if` cell-period allocator: the per-period IF (`U_pp_l_path`) is built with non-path switcher-side contributions skipped (control contributions are unchanged, matching the joiners/leavers IF convention; preserves the row-sum identity `U_pp.sum(axis=1) == U`), cohort-recentered via `_cohort_recenter_per_period`, then expanded to observations as `psi_i = U_pp[g_i, t_i] · (w_i / W_{g_i, t_i})`. Replicate-weight designs unconditionally use the cell allocator (Class A contract from PR #323). New `_refresh_path_inference` helper post-call refreshes `safe_inference` on every populated entry across `multi_horizon_inference`, `placebo_horizon_inference`, `path_effects`, and `path_placebos` so all four surfaces use the same final `df_survey` after per-path replicate fits append `n_valid` to the shared accumulator. Path-enumeration ranking under `survey_design` remains unweighted (group-cardinality, not population-weight mass). Lonely-PSU policy stays sample-wide, not per-path. Telescope invariant: on a single-path panel, per-path SE matches the global non-by_path survey SE bit-exactly. **No R parity** — R `did_multiplegt_dyn` does not support survey weighting; this is a Python-only methodology extension. The global non-by_path TSL multiplier-bootstrap path is unaffected (anti-regression test `tests/test_chaisemartin_dhaultfoeuille.py::TestByPathSurveyDesignAnalytical::test_global_survey_plus_n_bootstrap_still_works` locks the per-path-only scope of the new gate). Cross-surface invariants regression-tested at `TestByPathSurveyDesignAnalytical` (~17 tests across gate / dispatch / analytical SE / replicate-weight SE / per-path placebos / `trends_linear` composition / unobserved-path warnings / final-df refresh regressions) and `TestByPathSurveyDesignTelescope`. See `docs/methodology/REGISTRY.md` §`ChaisemartinDHaultfoeuille` `Note (Phase 3 by_path ...)` → "Per-path survey-design SE" for the full contract.
1414
- **Inference-field aliases on staggered result classes** for adapter / external-consumer compatibility. Read-only `@property` aliases expose the flat `att` / `se` / `conf_int` / `p_value` / `t_stat` names (matching `DiDResults` / `TROPResults` / `SyntheticDiDResults` / `HeterogeneousAdoptionDiDResults`) on every result class that previously only carried prefixed canonical fields: `CallawaySantAnnaResults`, `StackedDiDResults`, `EfficientDiDResults`, `ChaisemartinDHaultfoeuilleResults`, `StaggeredTripleDiffResults`, `WooldridgeDiDResults`, `SunAbrahamResults`, `ImputationDiDResults`, `TwoStageDiDResults` (mapping to `overall_*`); `ContinuousDiDResults` (mapping to `overall_att_*`, ATT-side as the headline, ACRT-side accessible unchanged via `overall_acrt_*`); `MultiPeriodDiDResults` (mapping to `avg_*`). `ContinuousDiDResults` additionally exposes `overall_se` / `overall_conf_int` / `overall_p_value` / `overall_t_stat` aliases for naming consistency with the rest of the staggered family. Aliases are pure read-throughs over the canonical fields — no recomputation, no behavior change — so the `safe_inference()` joint-NaN contract (per CLAUDE.md "Inference computation") is inherited automatically (NaN canonical → NaN alias, locked at `tests/test_result_aliases.py::test_pattern_b_aliases_propagate_nan`). The native `overall_*` / `overall_att_*` / `avg_*` fields remain canonical for documentation and computation. Motivated by the `balance.interop.diff_diff.as_balance_diagnostic()` adapter (`facebookresearch/balance` PR #465) which calls `getattr(res, "se", None)` / `getattr(res, "conf_int", None)` without a fallback chain — pre-alias, every staggered result class returned `None` on those keys, silently dropping `se` and `conf_int` from the adapter's diagnostic dict. 23 alias-mechanic + balance-adapter regression tests at `tests/test_result_aliases.py`. Patch-level (additive on stable surfaces).

diff_diff/conley.py

Lines changed: 43 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -86,10 +86,18 @@ def _pairwise_distance_matrix(coords: np.ndarray, metric) -> np.ndarray:
8686

8787

8888
def _bartlett_kernel(u: np.ndarray) -> np.ndarray:
89-
"""Bartlett (linear taper) kernel: K(u) = max(0, 1 - |u|).
90-
91-
Conley (1999) Eq 3.14 + Andrews (1991). PSD-guaranteed (non-negative
92-
spectral window), so the resulting Conley meat is PSD.
89+
"""Bartlett (linear taper) kernel on pairwise distance: K(u) = max(0, 1 - |u|).
90+
91+
This is the radial 1-D specialization of Conley (1999)'s Bartlett window
92+
that R ``conleyreg`` (Düsterhöft 2021), Stata ``acreg`` (Colella et al.
93+
2019), and Hsiang (2010) all use as their Bartlett path. Conley's
94+
explicit PSD formula (Eq 3.14, page 12) is the **2-D separable product
95+
window** ``K(j, k) = (1 - |j|/L_M)(1 - |k|/L_N)`` indexed on a lattice;
96+
the 1-D radial form on pairwise distance is a practitioner specialization
97+
that is not explicitly written in the paper and is therefore **not
98+
PSD-guaranteed**. The caller checks the resulting meat for indefiniteness
99+
and emits a ``UserWarning`` if the smallest eigenvalue is materially
100+
negative (regardless of kernel).
93101
"""
94102
return np.maximum(0.0, 1.0 - np.abs(u))
95103

@@ -99,7 +107,7 @@ def _uniform_kernel(u: np.ndarray) -> np.ndarray:
99107
100108
Cited as White (1980) truncated estimator; Conley (1999) page 11. Easier
101109
to interpret than Bartlett but the spectral window is negative in regions
102-
(Conley 1999 footnote 11), so the resulting meat is NOT guaranteed PSD.
110+
(Conley 1999 footnote 11), so the resulting meat is not guaranteed PSD.
103111
Caller emits ``UserWarning`` if any meat eigenvalue is materially negative.
104112
"""
105113
return (np.abs(u) <= 1.0).astype(np.float64)
@@ -128,8 +136,10 @@ def _validate_conley_kwargs(
128136
if coords is None:
129137
raise ValueError(
130138
"vcov_type='conley' requires conley_coords (n×2 array of [lat, lon] "
131-
"or projected coords). Pass via TwoWayFixedEffects(conley_coords=...) "
132-
"or compute_robust_vcov(conley_coords=...)."
139+
"or projected coords). Pass via LinearRegression(conley_coords=...) "
140+
"or compute_robust_vcov(conley_coords=...) on a cross-sectional "
141+
"design (Phase 1 supports cross-sectional Conley only; panel "
142+
"estimators are deferred to Phase 2)."
133143
)
134144
coords_arr = np.asarray(coords, dtype=np.float64)
135145
if coords_arr.ndim != 2 or coords_arr.shape[1] != 2:
@@ -209,9 +219,13 @@ def _compute_conley_vcov(
209219
210220
Notes
211221
-----
212-
For ``kernel == "uniform"`` the meat is not guaranteed PSD (Conley 1999
213-
footnote 11); a ``UserWarning`` is emitted if the smallest meat eigenvalue
214-
is materially negative (< -1e-12).
222+
Neither the uniform kernel (negative spectral regions, Conley 1999
223+
footnote 11) nor the **radial 1-D Bartlett** specialization implemented
224+
here is PSD-guaranteed. Conley's explicit PSD formula (Eq 3.14) is the
225+
2-D separable product window on a lattice; the radial pairwise form is
226+
a practitioner specialization (R ``conleyreg``, Stata ``acreg``, Hsiang
227+
2010) that is not formally PSD. We emit a ``UserWarning`` if the smallest
228+
meat eigenvalue is materially negative (< -1e-12) regardless of kernel.
215229
"""
216230
coords_arr = np.asarray(coords, dtype=np.float64)
217231
D = _pairwise_distance_matrix(coords_arr, metric)
@@ -243,18 +257,25 @@ def _compute_conley_vcov(
243257
"score matrix for NaN/Inf."
244258
)
245259

246-
# PSD guard for uniform kernel (Conley 1999 fn 11)
247-
if kernel == "uniform":
248-
eigvals = np.linalg.eigvalsh(meat)
249-
if eigvals.size and eigvals.min() < -1e-12:
250-
warnings.warn(
251-
f"Conley meat with uniform kernel has a negative eigenvalue "
252-
f"({eigvals.min():.2e}); the variance estimator is not "
253-
"guaranteed PSD. Consider conley_kernel='bartlett' (PSD by "
254-
"construction).",
255-
UserWarning,
256-
stacklevel=3,
257-
)
260+
# PSD guard. Neither the uniform kernel (Conley 1999 fn 11) nor the
261+
# radial 1-D Bartlett specialization is formally PSD-guaranteed —
262+
# Conley's explicit PSD Bartlett formula (Eq 3.14) is the 2-D separable
263+
# product window, not the 1-D radial pairwise form that R `conleyreg`,
264+
# Stata `acreg`, and this implementation use. Check both kernels.
265+
eigvals = np.linalg.eigvalsh(meat)
266+
if eigvals.size and eigvals.min() < -1e-12:
267+
warnings.warn(
268+
f"Conley meat with conley_kernel={kernel!r} has a materially "
269+
f"negative eigenvalue ({eigvals.min():.2e}); the variance "
270+
"estimator is not guaranteed PSD on this design. Both "
271+
"supported kernels (radial bartlett and uniform) are "
272+
"practitioner specializations of Conley 1999 and are not "
273+
"formally PSD-guaranteed; consider varying conley_cutoff_km "
274+
"or reviewing the design for collinearity / degenerate "
275+
"residual structure.",
276+
UserWarning,
277+
stacklevel=3,
278+
)
258279

259280
# Sandwich via two solves (mirrors _compute_cr2_bm pattern in linalg.py)
260281
try:

diff_diff/guides/llms-full.txt

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1918,8 +1918,10 @@ se = np.sqrt(np.diag(reg.vcov_))
19181918
Var̂(β) = (X'X)^{-1} · ( Σ_{i,j} K(d_ij / h) · X_i ε_i ε_j X_j' ) · (X'X)^{-1}
19191919

19201920
**Kernels:**
1921-
- `"bartlett"` (default): `K(u) = max(0, 1 - |u|)`. PSD-guaranteed.
1922-
- `"uniform"`: `K(u) = 1{|u| ≤ 1}`. Easier to interpret; emits `UserWarning` if the resulting meat has a materially negative eigenvalue.
1921+
- `"bartlett"` (default): `K(u) = max(0, 1 - |u|)` on pairwise distance `d_ij/h`. The radial 1-D form, matching R `conleyreg` / Stata `acreg`. Conley 1999's explicit PSD-guaranteed Bartlett formula (Eq 3.14) is the 2-D **separable product** window on a lattice; the 1-D radial specialization that diff-diff implements is a practitioner convention and is not formally PSD-guaranteed.
1922+
- `"uniform"`: `K(u) = 1{|u| ≤ 1}`. Easier to interpret.
1923+
1924+
Both kernels: `UserWarning` is emitted if the resulting meat has a materially negative eigenvalue (< -1e-12) — neither kernel is formally PSD-guaranteed in the radial 1-D pairwise-distance form.
19231925

19241926
**Distance metrics:**
19251927
- `"haversine"` (default): great-circle in km, Earth's mean radius 6371.01 km (matching R `conleyreg`). Validates `lat ∈ [-90, 90]`, `lon ∈ [-180, 180]`.

0 commit comments

Comments
 (0)