Skip to content

Commit a313e05

Browse files
igerberclaude
andcommitted
Address CI AI review round 6: Design 1' support check + empty guard
P1 #1: boundary=0 now enforces a Design 1' support plausibility heuristic: d.min() <= 5% * median(|d|). Samples with d.min() substantially positive (e.g. U(0.5, 1)) are rejected with ValueError directing the caller to boundary=float(d.min()). Threshold chosen at 5% (not REGISTRY's 1%) so the paper's thin-boundary-density DGPs (Beta(2,2), d.min/median ~ 3%) still pass. Reordered so the mass-point check (NotImplementedError, paper Section 3.2.4) fires before the support-check -- mass-point data should be redirected to 2SLS regardless of the boundary the caller picked. P1 #2: Empty-input front-door guard. d.size == 0 raises ValueError with a targeted "must be non-empty" message instead of leaking the NumPy reduction error from d.min(). P3 (docstring sync): _nprobust_port module docstring no longer says weighted data can be handled by the public wrapper -- the wrapper explicitly raises NotImplementedError. Docstring now matches the actual contract. P3 (deferred, same as last round): tri/uni/shifted-boundary golden parity extension. REGISTRY.md Phase 1b note expanded to document the full input contract (nonnegativity, boundary applicability, Design 1' support heuristic, mass-point redirection) so the public API surface is fully specified in the methodology registry. 178 tests pass (up from 177). Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
1 parent 196f910 commit a313e05

4 files changed

Lines changed: 67 additions & 19 deletions

File tree

diff_diff/_nprobust_port.py

Lines changed: 6 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -31,9 +31,12 @@
3131
3232
Deviations from nprobust (documented):
3333
34-
* ``weights=`` is not supported here (nprobust's ``lpbwselect`` has no
35-
weight argument). Weighted data can be handled by the public wrapper
36-
or via the user's own design matrix.
34+
* ``weights=`` is not supported here or in the public wrapper
35+
(nprobust's ``lpbwselect`` has no weight argument, so Phase 1b has
36+
no parity anchor). Weighted-data support is queued for Phase 2+
37+
(survey-design adaptation). The public wrapper
38+
``mse_optimal_bandwidth`` raises ``NotImplementedError`` when a
39+
``weights`` array is passed.
3740
* ``vce="nn"`` is the default and is fully ported. ``vce in
3841
{"hc0", "hc1", "hc2", "hc3"}`` is implemented in ``lprobust_res`` /
3942
``lprobust_vce`` but has not been separately golden-tested; use at

diff_diff/local_linear.py

Lines changed: 45 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -615,6 +615,11 @@ def mse_optimal_bandwidth(
615615
y = np.asarray(y, dtype=np.float64).ravel()
616616
if d.shape != y.shape:
617617
raise ValueError(f"d and y must have the same shape; got {d.shape} and {y.shape}")
618+
if d.size == 0:
619+
raise ValueError(
620+
"d and y must be non-empty; the selector cannot estimate a "
621+
"bandwidth from zero observations."
622+
)
618623
if not np.all(np.isfinite(d)):
619624
raise ValueError("d contains non-finite values (NaN or Inf)")
620625
if not np.all(np.isfinite(y)):
@@ -668,13 +673,15 @@ def mse_optimal_bandwidth(
668673
)
669674

670675
# Mass-point design check (paper Section 3.2.4, REGISTRY 2% rule).
671-
# When d_min > 0 and there is bunching at d_min, Design 1 requires
672-
# the 2SLS sample-average path (Phase 2), not the CCF nonparametric
673-
# selector. The check applies independently of the boundary the
674-
# user supplied: mass-point data is never appropriate for this
675-
# wrapper. The check explicitly excludes d_min ~ 0, which is the
676-
# Design 1' "untreated units present" subcase that the paper's
677-
# simulations and the Garrett et al. (2020) application accept.
676+
# Must fire BEFORE the Design 1' support check: mass-point data is
677+
# never appropriate for the CCF nonparametric selector regardless
678+
# of the boundary the caller supplied. The correct remediation is
679+
# the 2SLS sample-average path (Phase 2), not a boundary
680+
# reclassification.
681+
#
682+
# The check explicitly excludes d_min ~ 0 (the Design 1'
683+
# "untreated units present" subcase that the paper's simulations
684+
# and the Garrett et al. (2020) application accept).
678685
_MASS_POINT_THRESHOLD = 0.02 # REGISTRY rule: > 2% modal-min
679686
if d_min > _boundary_tol:
680687
eps_eq = 1e-12 * max(1.0, abs(d_min))
@@ -695,6 +702,37 @@ def mse_optimal_bandwidth(
695702
f"applicable."
696703
)
697704

705+
# Design 1' support check: boundary ~ 0 requires the realized
706+
# sample minimum to be compatible with a population support
707+
# infimum at 0. Otherwise the selector calibrates ``h_mse`` at
708+
# an off-support limit.
709+
#
710+
# Rule: when boundary ~ 0 (not also at d.min()), require
711+
# d.min() <= 5% * median(|d|). The 5% threshold is generous
712+
# enough to accept Design 1' samples with vanishing boundary
713+
# density (Beta(2,2): d.min/median ~ 3%) while rejecting samples
714+
# substantially off-support (U(0.5, 1): d.min/median ~ 1.0).
715+
# Samples just between these (e.g. U(0.05, 1), d.min/median ~ 10%)
716+
# are directed to boundary=float(d.min()) for the continuous-
717+
# near-d_lower path.
718+
_DESIGN_1_PRIME_RATIO = 0.05
719+
if _at_zero and not _at_d_min:
720+
d_median_abs = float(np.median(np.abs(d)))
721+
effective_threshold = _DESIGN_1_PRIME_RATIO * max(d_median_abs, 1e-12)
722+
if d_min > effective_threshold:
723+
raise ValueError(
724+
f"boundary ~ 0 selected but d.min()={d_min!r} is not "
725+
f"compatible with a Design 1' support infimum at 0 "
726+
f"(rule: d.min() <= "
727+
f"{_DESIGN_1_PRIME_RATIO} * median(|d|) = "
728+
f"{effective_threshold!r}). This sample is not "
729+
f"Design 1'. Either: (a) pass boundary=float(d.min()) "
730+
f"for the Design 1 continuous-near-d_lower path, or "
731+
f"(b) verify the population support actually has "
732+
f"infimum at 0 (in which case the realized d.min() "
733+
f"would be closer to zero relative to the data scale)."
734+
)
735+
698736
# Defer heavy import to call time to avoid import-cycle risk.
699737
from diff_diff._nprobust_port import lpbwselect_mse_dpi
700738

docs/methodology/REGISTRY.md

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2302,7 +2302,7 @@ Shipped as `did_had_pretest_workflow()` and surfaced via `practitioner_next_step
23022302
- [x] Phase 1a: `vcov_type` enum threaded through `DifferenceInDifferences` (`MultiPeriodDiD`, `TwoWayFixedEffects` inherit); `robust=True` <=> `vcov_type="hc1"`, `robust=False` <=> `vcov_type="classical"`. Conflict detection at `__init__`. Results summary prints the variance-family label.
23032303
- **Note (deviation from the fully-symmetric enum):** `MultiPeriodDiD(cluster=..., vcov_type="hc2_bm")` is intentionally **not supported** and raises `NotImplementedError`. The scalar-coefficient `DifferenceInDifferences` path handles the cluster + CR2 Bell-McCaffrey combination (`_compute_cr2_bm` returns a per-coefficient Satterthwaite DOF that is valid for the single-ATT contrast), but `MultiPeriodDiD` also reports a post-period-average ATT constructed as a *contrast* of the event-study coefficients. The cluster-aware CR2 BM DOF for that contrast (i.e., the Pustejovsky-Tipton 2018 per-cluster adjustment matrices applied to an arbitrary aggregation contrast) is not yet implemented. Pairing CR2 cluster-robust SEs with the one-way Imbens-Kolesar (2016) contrast DOF would be a broken hybrid, so the combination fails fast with a clear workaround message (drop the cluster for one-way HC2+BM, or use `vcov_type="hc1"` with cluster for CR1 Liang-Zeger). Tracked in `TODO.md` under Methodology/Correctness. Applies only to `MultiPeriodDiD`; `DifferenceInDifferences(cluster=..., vcov_type="hc2_bm")` works.
23042304
- [x] Phase 1a: `clubSandwich::vcovCR(..., type="CR2")` parity harness committed: R script at `benchmarks/R/generate_clubsandwich_golden.R` plus a regression-anchor JSON at `benchmarks/data/clubsandwich_cr2_golden.json`. **Note:** the committed JSON currently has `"source": "python_self_reference"` and pins numerical stability only; authoritative R-produced values are generated by running the R script, which the TODO.md row under Methodology/Correctness tracks. The parity test at `tests/test_linalg_hc2_bm.py::TestCR2BMCluster::test_cr2_parity_with_golden` runs at 1e-6 tolerance (Phase 1a plan commits 6-digit parity once R regen completes).
2305-
- [x] Phase 1b: Calonico-Cattaneo-Farrell (2018) MSE-optimal bandwidth selector. In-house port of `nprobust::lpbwselect(bwselect="mse-dpi")` (nprobust 0.5.0, SHA `36e4e53`) as `diff_diff.mse_optimal_bandwidth` and `BandwidthResult`, backed by the private `diff_diff._nprobust_port` module (`kernel_W`, `lprobust_bw`, `lpbwselect_mse_dpi`). Three-stage DPI with four `lprobust.bw` calls at orders `q+1`, `q+2`, `q`, `p`. Parity verified at `0.0000%` on all five stage bandwidths (`c_bw`, `bw_mp2`, `bw_mp3`, `b_mse`, `h_mse`) across three deterministic DGPs (uniform, Beta(2,2), half-normal) via `benchmarks/R/generate_nprobust_golden.R` → `benchmarks/data/nprobust_mse_dpi_golden.json`. **Note:** `weights=` is currently unsupported (raises `NotImplementedError`); nprobust's `lpbwselect` has no weight argument so there is no parity anchor. Weighted-data support deferred to Phase 2 (survey-design adaptation). **Note (public API scope restriction):** the exported wrapper `mse_optimal_bandwidth` hard-codes the HAD Phase 1b configuration (`p=1`, `deriv=0`, `interior=False`, `vce="nn"`, `nnmatch=3`). The underlying port supports a broader surface (`hc0`/`hc1`/`hc2`/`hc3` variance, interior evaluation, higher `p`), but those paths are not parity-tested against `nprobust` and are deferred. Callers needing the broader surface should use `diff_diff._nprobust_port.lpbwselect_mse_dpi` directly and accept that parity has not been verified on non-HAD configurations.
2305+
- [x] Phase 1b: Calonico-Cattaneo-Farrell (2018) MSE-optimal bandwidth selector. In-house port of `nprobust::lpbwselect(bwselect="mse-dpi")` (nprobust 0.5.0, SHA `36e4e53`) as `diff_diff.mse_optimal_bandwidth` and `BandwidthResult`, backed by the private `diff_diff._nprobust_port` module (`kernel_W`, `lprobust_bw`, `lpbwselect_mse_dpi`). Three-stage DPI with four `lprobust.bw` calls at orders `q+1`, `q+2`, `q`, `p`. Parity verified at `0.0000%` on all five stage bandwidths (`c_bw`, `bw_mp2`, `bw_mp3`, `b_mse`, `h_mse`) across three deterministic DGPs (uniform, Beta(2,2), half-normal) via `benchmarks/R/generate_nprobust_golden.R` → `benchmarks/data/nprobust_mse_dpi_golden.json`. **Note:** `weights=` is currently unsupported (raises `NotImplementedError`); nprobust's `lpbwselect` has no weight argument so there is no parity anchor. Weighted-data support deferred to Phase 2 (survey-design adaptation). **Note (public API scope restriction):** the exported wrapper `mse_optimal_bandwidth` hard-codes the HAD Phase 1b configuration (`p=1`, `deriv=0`, `interior=False`, `vce="nn"`, `nnmatch=3`). The underlying port supports a broader surface (`hc0`/`hc1`/`hc2`/`hc3` variance, interior evaluation, higher `p`), but those paths are not parity-tested against `nprobust` and are deferred. Callers needing the broader surface should use `diff_diff._nprobust_port.lpbwselect_mse_dpi` directly and accept that parity has not been verified on non-HAD configurations. **Note (input contract):** the wrapper enforces HAD's support restriction `D_{g,2} >= 0` (front-door `ValueError` on negative doses and empty inputs). `boundary` must equal `0` (Design 1') or `float(d.min())` (Design 1 continuous-near-d_lower) within float tolerance; off-support values raise `ValueError`. When `boundary ~ 0`, the wrapper additionally requires `d.min() <= 0.05 * median(|d|)` as a Design 1' support plausibility heuristic, chosen to pass the paper's thin-boundary-density DGPs (Beta(2,2), d.min/median ~ 3%) while rejecting substantially off-support samples (U(0.5, 1.0), d.min/median ~ 1.0). Detected mass-point designs (`d.min() > 0` with modal fraction at `d.min() > 2%`) raise `NotImplementedError` pointing to the Phase 2 2SLS path per paper Section 3.2.4.
23062306
- [ ] Phase 1c: First-order bias estimator `M̂_{ĥ*_G}` and robust variance `V̂_{ĥ*_G}`.
23072307
- [ ] Phase 1c: Bias-corrected CI (Equation 8) with `nprobust` parity.
23082308
- [ ] Phase 2: `HeterogeneousAdoptionDiD` class with separate code paths for Design 1', Design 1 mass-point, and Design 1 continuous-near-``.

tests/test_bandwidth_selector.py

Lines changed: 15 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -332,18 +332,25 @@ def test_boundary_zero_thin_boundary_density_accepted(self):
332332
assert np.isfinite(h)
333333
assert h > 0.0
334334

335-
def test_boundary_zero_with_data_far_from_zero_fails_gracefully(self):
336-
"""boundary=0 passes the boundary-validation and mass-point
337-
checks but then hits the per-stage count guard deeper in the
338-
selector because the kernel window is empty (d ~ U(0.5, 1.0)
339-
has no data near 0). Must surface a clear ValueError, not an
340-
opaque failure."""
335+
def test_boundary_zero_with_data_far_from_zero_rejected(self):
336+
"""boundary=0 with d.min() substantially positive fails the
337+
Design 1' support check (d.min() > 1% of median(|d|)). The
338+
caller must either pass boundary=float(d.min()) for the
339+
continuous-near-d_lower path or confirm Design 1' applicability."""
341340
rng = np.random.default_rng(2026)
342-
d = rng.uniform(0.5, 1.0, size=1500) # d.min() ~ 0.5, no mass
341+
d = rng.uniform(0.5, 1.0, size=1500) # d.min() ~ 0.5 >> 1% of median
343342
y = d + rng.normal(0, 0.3, size=1500)
344-
with pytest.raises(ValueError, match="lprobust_bw"):
343+
with pytest.raises(ValueError, match="Design 1'"):
345344
mse_optimal_bandwidth(d, y, boundary=0.0)
346345

346+
def test_empty_input_rejected(self):
347+
"""Empty d/y must raise a targeted ValueError up front, not
348+
leak a NumPy reduction error from d.min()."""
349+
d = np.array([], dtype=np.float64)
350+
y = np.array([], dtype=np.float64)
351+
with pytest.raises(ValueError, match="non-empty"):
352+
mse_optimal_bandwidth(d, y)
353+
347354
def test_boundary_zero_with_d_min_mass_point_rejected(self):
348355
"""boundary=0 with d.min() > 0 AND mass at d.min() is a
349356
Design 1 mass-point design and must be redirected to 2SLS."""

0 commit comments

Comments
 (0)