From 0a40a5ee568e8283ac3bce89d0fc812bb3995e7e Mon Sep 17 00:00:00 2001 From: igerber Date: Mon, 6 Jul 2026 16:51:35 -0400 Subject: [PATCH] test(se-audit): lock estimatr intercept SE, CR2 DOF, LOO df, Yatchew p/sigma2, fixest cluster band Second "fiddly-bits" coverage batch of the SE parity audit. Pins previously- unasserted golden SE/DOF/intercept fields that already matched R but had no CI assertion (the latent-risk pattern that once hid the CallawaySantAnna reg gap): - C4: estimatr HC1/CR1 2SLS intercept SE via an opt-in return_intercept_se hook on the private _fit_mass_point_2sls (surfaces sqrt(V[0,0]); @overload keeps every caller's 3-tuple return byte-unchanged). classical intercept excluded (O(1/n)). - C2: TWFE dof_hc2_bm and MPD per-period dof_per_coef via CI-inversion (no new result fields; reconstruct conf_int from golden DOF, att+se already pinned). - C3: PlaceboTests leave-one-out df via the public leave_one_out_effects count. - C5: Yatchew p_value + sigma2_lin + sigma2_diff (through the N/(N-1) shift). - G2: fixest DiD/TWFE cluster-SE band pin within the documented ~0.25% CR1 band. C1 (StackedDiD intercept SE) dropped after measurement: it diverges ~0.3% from R (nuisance-parameter reference-cell gap, not machine-precision lockable; the event- study interaction SEs already match ~2e-13) - documented in TODO. Remaining deferrals (G2 hetero golden regen, boundary_gap, C6-C8) tracked in TODO. Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_0178r8ZK6VRiXbaBiknrWjQB --- CHANGELOG.md | 14 ++++ TODO.md | 2 +- diff_diff/had.py | 94 +++++++++++++++++++++---- docs/methodology/REGISTRY.md | 1 + tests/test_did_had_parity.py | 40 +++++++++++ tests/test_estimators_vcov_type.py | 32 +++++++-- tests/test_estimatr_iv_robust_parity.py | 66 +++++++++++++++++ tests/test_fixest_did_twfe_parity.py | 8 +++ tests/test_methodology_placebo.py | 6 ++ 9 files changed, 244 insertions(+), 19 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 6681fdb99..5ed25e87c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -22,6 +22,20 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - **StackedDiD** event-study coefficient point estimates; **PlaceboTests** leave-one-out `t_stat`/`p_value`; **DIDHAD** QUG order-statistic test (`t_stat`/`p_value`, previously unasserted against the golden). All at ~1e-10. + - **estimatr 2SLS intercept SE** (HAD mass-point path): the `HC1`/`CR1` intercept SE now + matches `estimatr::iv_robust` `se_intercept` to ~1e-15 (the 2×2 sandwich already computed + `V[0,0]`; an opt-in `return_intercept_se=` hook on the private `_fit_mass_point_2sls` + surfaces it — the default 3-tuple return is byte-unchanged, no public-API change). The + `classical` intercept is excluded (same documented `O(1/n)` projection/DOF deviation as the slope). + - **CR2 Satterthwaite DOF** via CI-inversion (no new result fields): `TwoWayFixedEffects(hc2_bm)` + `dof_hc2_bm` and `MultiPeriodDiD(hc2_bm)` per-period `dof_per_coef` are pinned by reconstructing + `conf_int` from the golden DOF — with ATT+SE already locked, the CI matches iff the DOF matches. + - **PlaceboTests** leave-one-out `df` (recovered from the public `leave_one_out_effects` count); + **DIDHAD** Yatchew linearity test `p_value` and both `sigma2` components (through the documented + `N/(N-1)` sample-vs-population convention shift). All at ~1e-10. + - **`fixest` cluster-robust SE band**: the DiD/TWFE cluster-at-unit SE is pinned within the + documented ~0.25% fixest-CR1 small-sample DOF-convention band (guards an unintended SE-formula + change; the machine-precision hetero/cluster lock is deferred — needs an unbalanced-DGP golden). ### Added - **`ContinuousDiD` lowest-dose-as-control** (`control_group="lowest_dose"`, CGBS 2024 Remark 3.1) for diff --git a/TODO.md b/TODO.md index 52b16feda..4c0b756a2 100644 --- a/TODO.md +++ b/TODO.md @@ -66,7 +66,7 @@ generic sparse-FE, QR+SVD rank-detection redundancy, `check_finite` bypass — m | Issue | Location | Origin | Effort | Priority | |-------|----------|--------|--------|----------| -| SE-audit CI-lock follow-up ("fiddly bits" deferred from the G2+C1-C3/C5 coverage batch). The `fixest_did_twfe_golden.json` stores `cluster_unit` ATT+SE for DiD and TWFE but only the ATT is asserted (the cluster SE carries the documented CR1 DOF-convention difference vs fixest; the hetero SE collapses to iid on these balanced 2-group designs). Also deferred, each needing plumbing or convention work: surface the 2SLS intercept SE for estimatr `se_intercept` parity (C4); surface StackedDiD intercept SEs (`se_cr1/cr2_intercept`, C1); add a `df` field + `boundary_gap` computation for PlaceboTests (C3); assert Yatchew `p`/`sigma2_lin/diff` through the N/(N-1) convention shift + event-map (C5); assert `dof_hc2_bm`/`dof_per_coef` via CI-inversion (C2). Plus the tolerance-tightenings C6-C8. | `tests/test_fixest_did_twfe_parity.py`, `tests/test_estimatr_iv_robust_parity.py`, `tests/test_methodology_stacked_did.py`, `tests/test_methodology_placebo.py`, `tests/test_did_had_parity.py` | SE-audit | Mid | Low | +| SE-audit CI-lock — remaining "fiddly bits" after the second coverage batch landed (that batch pinned C2 `dof_hc2_bm`/`dof_per_coef` via CI-inversion, C3 LOO `df`, C4 estimatr HC1/CR1 intercept SE, C5 Yatchew `p`/`sigma2_lin`/`sigma2_diff`, and the G2 fixest cluster-SE band). Still deferred, each needing a golden regeneration, new computation, or a documented-deviation call: **(a) G2 machine-precision hetero/cluster SE** — needs an unbalanced/heteroskedastic-DGP regen of `fixest_did_twfe_golden.json` (hetero collapses to iid on the current balanced design; the cluster SE is the documented ~0.25% fixest-CR1 DOF-convention deviation, currently band-pinned only); **(b) PlaceboTests `boundary_gap`** — a permutation randomization-inference margin NOT computed anywhere in code (a new feature + result field, not a coverage lock); **(c) StackedDiD intercept SEs** (`se_cr1/cr2_intercept`, C1) — MEASURED to diverge ~0.3% from R: a nuisance-parameter reference-cell/parameterization gap, NOT machine-precision lockable (the event-study interaction SEs already match ~2e-13; surfacing it would add an unasserted, R-divergent public field); **(d) estimatr `classical` intercept SE** — same documented `O(1/n)` projection/DOF deviation as the slope (reference-only, excluded from parity). Plus the tolerance-tightenings C6-C8. | `benchmarks/R/generate_fixest_did_twfe_golden.R`, `tests/test_fixest_did_twfe_parity.py`, `tests/test_methodology_stacked_did.py`, `tests/test_methodology_placebo.py` | SE-audit | Mid | Low | | Render `docs/methodology/REPORTING.md` and `REGISTRY.md` as in-site Sphinx pages so cross-refs can use `:doc:` instead of off-site `blob/main` URLs (stable-docs readers can otherwise land on a different revision than their package version). Two paths: (a) add `myst-parser` to `conf.py` + docs extras and link with `:doc:`, or (b) convert both to `.rst`. **Note:** REGISTRY.md is ~4.5k lines of LaTeX-heavy markdown — high risk under the `-W` (warnings-as-errors) Sphinx build; budget multiple rounds. | `docs/conf.py`, `docs/api/business_report.rst`, `docs/api/diagnostic_report.rst`, tutorials 18 & 19 | follow-up | Mid | Low | | `ImputationDiD` covariate-path variance lacks a dedicated parity anchor — only the no-covariate staggered panel is R-parity'd, though the covariate path shares the same validated projection code. Add a small dense-design **hand-calc** for the covariate projection (no external tooling), or a covariate (time-varying X) R `didimputation` golden asserting overall/ES SE parity (the golden variant needs local R). | `tests/test_methodology_imputation.py`, `benchmarks/R/generate_didimputation_golden.R` | imputation-validation | Mid | Low | | Add true half-sample BRR replicate-weight regressions per estimator family (current tests use Fay-like 0.5/1.5 perturbations; `test_survey_phase6.py` covers true BRR at the helper level). | `tests/test_replicate_weight_expansion.py` | #253 | Mid | Low | diff --git a/diff_diff/had.py b/diff_diff/had.py index e1a525240..ad5e80426 100644 --- a/diff_diff/had.py +++ b/diff_diff/had.py @@ -66,7 +66,7 @@ import warnings from dataclasses import dataclass -from typing import Any, Dict, List, Optional, Tuple, Union +from typing import Any, Dict, List, Literal, Optional, Tuple, Union, overload import numpy as np import pandas as pd @@ -2326,6 +2326,34 @@ def _sup_t_multiplier_bootstrap( # ============================================================================= +@overload +def _fit_mass_point_2sls( + d: np.ndarray, + dy: np.ndarray, + d_lower: float, + cluster: Optional[np.ndarray], + vcov_type: str, + *, + weights: Optional[np.ndarray] = ..., + return_influence: bool = ..., + return_intercept_se: Literal[False] = ..., +) -> Tuple[float, float, Optional[np.ndarray]]: ... + + +@overload +def _fit_mass_point_2sls( + d: np.ndarray, + dy: np.ndarray, + d_lower: float, + cluster: Optional[np.ndarray], + vcov_type: str, + *, + weights: Optional[np.ndarray] = ..., + return_influence: bool = ..., + return_intercept_se: Literal[True], +) -> Tuple[float, float, Optional[np.ndarray], float]: ... + + def _fit_mass_point_2sls( d: np.ndarray, dy: np.ndarray, @@ -2335,7 +2363,11 @@ def _fit_mass_point_2sls( *, weights: Optional[np.ndarray] = None, return_influence: bool = False, -) -> Tuple[float, float, Optional[np.ndarray]]: + return_intercept_se: bool = False, +) -> Union[ + Tuple[float, float, Optional[np.ndarray]], + Tuple[float, float, Optional[np.ndarray], float], +]: """Wald-IV point estimate and structural-residual 2SLS sandwich SE. The just-identified binary instrument ``Z_g = 1{D_{g,2} > d_lower}`` @@ -2390,6 +2422,13 @@ def _fit_mass_point_2sls( ``compute_survey_if_variance(IF, trivial)`` ≈ ``V_HC1[1, 1]`` at ``atol=1e-10`` (PR #359 convention; see the "IF scale convention" section of the Phase 4.5 B plan for derivation). + return_intercept_se : bool + When True, the returned tuple gains a 4th element: the intercept + standard error ``sqrt(V[0, 0])`` (``V[1, 1]`` is the slope + variance that drives ``se_beta``). ``NaN`` whenever ``se_beta`` + is ``NaN`` (early-exit / singular paths). Test-only surfacing hook + for estimatr-parity coverage; production callers leave it False and + the return is the unchanged 3-tuple. Returns ------- @@ -2398,7 +2437,9 @@ def _fit_mass_point_2sls( ``return_influence=True``, else ``None``. NaN for SE when the dose-gap vanishes (``Dbar_{Z=1} == Dbar_{Z=0}``) or the sandwich is singular; in those cases ``psi`` is returned as a - length-``n`` zero array when ``return_influence=True``. + length-``n`` zero array when ``return_influence=True``. When + ``return_intercept_se=True`` a 4th element ``se_intercept`` is + appended. """ d = np.asarray(d, dtype=np.float64) dy = np.asarray(dy, dtype=np.float64) @@ -2443,8 +2484,25 @@ def _fit_mass_point_2sls( _null_psi = np.zeros(n, dtype=np.float64) if return_influence else None + def _pack( + beta: float, + se_b: float, + psi_out: Optional[np.ndarray], + se_int: float = float("nan"), + ) -> Union[ + Tuple[float, float, Optional[np.ndarray]], + Tuple[float, float, Optional[np.ndarray], float], + ]: + # Centralize the optional 4th (intercept-SE) element so every one of + # the function's early-exit and success returns keeps a consistent + # arity. The @overload stubs give callers the precise per-flag type; + # the True branch appends se_int for the test-only surfacing hook. + if return_intercept_se: + return (beta, se_b, psi_out, se_int) + return (beta, se_b, psi_out) + if n_above == 0 or n_at_or_below == 0: - return float("nan"), float("nan"), _null_psi + return _pack(float("nan"), float("nan"), _null_psi) # Point estimate: weighted Wald-IV ratio (reduces to unweighted at w=1). if weighted: @@ -2453,7 +2511,7 @@ def _fit_mass_point_2sls( w_Z1 = float(w_arr[Z1_idx].sum()) w_Z0 = float(w_arr[Z0_idx].sum()) if w_Z1 <= 0.0 or w_Z0 <= 0.0: - return float("nan"), float("nan"), _null_psi + return _pack(float("nan"), float("nan"), _null_psi) dose_gap = float( (w_arr[Z1_idx] * d[Z1_idx]).sum() / w_Z1 - (w_arr[Z0_idx] * d[Z0_idx]).sum() / w_Z0 ) @@ -2470,7 +2528,7 @@ def _fit_mass_point_2sls( if abs(dose_gap) < 1e-12 * max(1.0, abs(d_bar)): # No dose variation around d_lower -> beta undefined. - return float("nan"), float("nan"), _null_psi + return _pack(float("nan"), float("nan"), _null_psi) # dy_gap / dy_bar were computed inside the weighted/unweighted blocks above. beta_hat = float(dy_gap / dose_gap) @@ -2495,7 +2553,7 @@ def _fit_mass_point_2sls( try: ZtWX_inv = np.linalg.inv(ZtWX) except np.linalg.LinAlgError: - return beta_hat, float("nan"), _null_psi + return _pack(beta_hat, float("nan"), _null_psi) vcov_type = vcov_type.lower() if vcov_type in _MASS_POINT_VCOV_UNSUPPORTED: @@ -2525,25 +2583,25 @@ def _fit_mass_point_2sls( n_clusters = len(clusters_unique) if n_clusters < 2: # Cluster-robust SE undefined with a single cluster. - return beta_hat, float("nan"), _null_psi + return _pack(beta_hat, float("nan"), _null_psi) Omega *= (n_clusters / (n_clusters - 1)) * ((n - 1) / (n - k)) elif vcov_type == "classical": if weighted: dof = w_sum - k if dof <= 0: - return beta_hat, float("nan"), _null_psi + return _pack(beta_hat, float("nan"), _null_psi) sigma2 = float((w_arr * w_arr * u * u).sum()) / dof Omega = sigma2 * (Zd.T @ ((w_arr * w_arr)[:, None] * Zd)) else: dof = n - k if dof <= 0: - return beta_hat, float("nan"), _null_psi + return _pack(beta_hat, float("nan"), _null_psi) sigma2 = float((u * u).sum()) / dof Omega = sigma2 * (Zd.T @ Zd) elif vcov_type == "hc1": dof = n - k if dof <= 0: - return beta_hat, float("nan"), _null_psi + return _pack(beta_hat, float("nan"), _null_psi) if weighted: # Pweight HC1: meat = Z' diag(w² u²) Z (Wooldridge 2010 Eq 12.37, # matches linalg.py:1141 convention and estimatr::iv_robust @@ -2561,11 +2619,19 @@ def _fit_mass_point_2sls( V = ZtWX_inv @ Omega @ ZtWX_inv.T var_beta = float(V[1, 1]) if not np.isfinite(var_beta) or var_beta < 0: - return beta_hat, float("nan"), _null_psi + return _pack(beta_hat, float("nan"), _null_psi) se_beta = float(np.sqrt(var_beta)) + # Intercept SE = sqrt(V[0, 0]); surfaced only when return_intercept_se + # (see _pack). NaN if the intercept variance is non-finite/negative. + var_intercept = float(V[0, 0]) + se_intercept = ( + float(np.sqrt(var_intercept)) + if np.isfinite(var_intercept) and var_intercept >= 0 + else float("nan") + ) if not return_influence: - return beta_hat, se_beta, None + return _pack(beta_hat, se_beta, None, se_intercept) # Per-unit influence function on β̂-scale, scaled so that # compute_survey_if_variance(psi, trivial_resolved) ≈ V_HC1[1, 1] @@ -2583,7 +2649,7 @@ def _fit_mass_point_2sls( psi0 = (Zd @ bread_row) * u dof_psi = max(n - k, 1) psi = psi0 * np.sqrt((n - 1) / dof_psi) - return beta_hat, se_beta, psi + return _pack(beta_hat, se_beta, psi, se_intercept) # ============================================================================= diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index 75eb1babc..1918a9003 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -3167,6 +3167,7 @@ Under `survey=SurveyDesign(weights, strata, psu, fpc)`, the variance composes vi - **CR1 pweight-cluster meat**: for each cluster c, `s_c = Z'_c·(w·u)_c`; `Ω_CR1 = (G/(G-1))·((n-1)/(n-k))·Σ_c s_c s_c'` (`w¹` inside cluster score). Bit-exact with `estimatr::iv_robust(..., weights=, clusters=, se_type="stata")` at `atol=1e-10`. - **Classical**: sandwich form `Ω_cl = σ²·Z'·diag(w²)·Z` with `σ² = Σw²u²/(Σw-k)`. Deviates from `estimatr` classical (projection-form + `n-k` DOF) by `O(1/n)` at non-uniform weights; unweighted path is bit-exact by equivalence. Skipped in cross-language parity tests. - **Per-unit IF on β̂-scale** (for Binder-TSL survey composition): `psi_g = [(Z'WX)^{-1} · z_g · w_g · u_g][1] · sqrt((n-1)/(n-k))`. The scaling factor absorbs DOF / small-sample differences so `compute_survey_if_variance(psi, trivial_resolved) ≈ V_HC1[1,1]` at `atol=1e-10` (mirrors PR #359 convention; asserted by `TestIFScaleInvariant` and bit-exact against estimatr HC1 on 4 DGPs). Fixture: `benchmarks/R/generate_estimatr_iv_robust_golden.R` → `benchmarks/data/estimatr_iv_robust_golden.json`. +- **Intercept SE (`return_intercept_se=True`)**: the 2×2 sandwich `V` already carries the intercept variance `V[0,0]`; the opt-in hook surfaces `sqrt(V[0,0])` (default off — the production 3-tuple return is byte-unchanged, so no public-API change). HC1 and CR1 intercept SEs are bit-exact with `estimatr` `se_intercept` at `atol=1e-10` (`tests/test_estimatr_iv_robust_parity.py`); the **classical** intercept carries the same `O(1/n)` projection/DOF deviation as the slope (bullet above) and is likewise excluded from the parity lock. *Event-study survey composition (Phase 4.5 B):* The per-horizon loop in `_fit_event_study` threads `weights_unit_full` + `resolved_survey_unit_full` through to both `_fit_continuous` and `_fit_mass_point_2sls` (the latter with `return_influence=True` under weighted fits). The returned IF matrix `Psi ∈ R^{G × H}` has a shared construction contract across paths — each column on the β̂-scale, such that `compute_survey_if_variance(Psi[:, e], resolved) ≈ V_β[e]`. Per-horizon analytical variance uses Binder-TSL via `compute_survey_if_variance` (survey= path) or the weighted-robust HC1 sandwich (weights= shortcut). `survey_metadata`, `variance_formula` (`"survey_binder_tsl"` / `"survey_binder_tsl_2sls"` / `"pweight"` / `"pweight_2sls"`), and `effective_dose_mean` populate identically to the static path. Pre-PR numerical output is preserved bit-exactly on the unweighted path when `cband=False` (stability invariant; Phase 2b convention unchanged for unweighted fits). diff --git a/tests/test_did_had_parity.py b/tests/test_did_had_parity.py index 6e955e749..466aaef9a 100644 --- a/tests/test_did_had_parity.py +++ b/tests/test_did_had_parity.py @@ -431,6 +431,17 @@ def test_yatchew_t_stat_parity( r_yatchew_t = _as_list(r_result["yatchew_t"]) r_yatchew_n = _as_list(r_result["yatchew_n"]) r_event_ids = _as_list(r_result["event_id"]) + # SE-audit C5: also lock the yatchew p-value and both sigma^2 columns + # (previously only T_hr was asserted). All three carry the SAME + # documented N/(N-1) convention shift as T_hr: R's YatchewTest uses + # base-R `var()` (1/(N-1)) where our estimator uses the population + # (1/G, 1/2G) convention, so R = ours x G/(G-1). Measured to lock at + # ~1e-15 across every row (incl. the literal-2G sigma2_diff divisor). + r_yatchew_p = _as_list(r_result["yatchew_p"]) + r_yatchew_s2_lin = _as_list(r_result["yatchew_sigma2_lin"]) + r_yatchew_s2_diff = _as_list(r_result["yatchew_sigma2_diff"]) + from scipy.stats import norm as _norm + for r_idx, r_id in enumerate(r_event_ids): e = _r_id_to_event_time(int(r_id), trends_lin) if e not in dy_dict: @@ -459,6 +470,35 @@ def test_yatchew_t_stat_parity( f"after N/(N-1) convention shift" ), ) + # p-value: shift T first, THEN evaluate the normal CDF (a p-value + # cannot be rescaled directly). 1 - Phi(T_hr x G/(G-1)) vs golden. + _fac = G_horizon / (G_horizon - 1) + np.testing.assert_allclose( + 1.0 - _norm.cdf(py_t_in_r_convention), + float(r_yatchew_p[r_idx]), + atol=YATCHEW_ATOL, + rtol=0, + err_msg=( + f"{dgp_name}/{combo_name}/Yatchew row {r_idx}: " + f"p-value mismatch after N/(N-1) shift" + ), + ) + # sigma^2 columns carry the same G/(G-1) convention shift as T_hr + # (sigma2_diff holds despite our literal-2G divisor — measured). + np.testing.assert_allclose( + float(r.sigma2_lin) * _fac, + float(r_yatchew_s2_lin[r_idx]), + atol=YATCHEW_ATOL, + rtol=0, + err_msg=(f"{dgp_name}/{combo_name}/Yatchew row {r_idx}: sigma2_lin mismatch"), + ) + np.testing.assert_allclose( + float(r.sigma2_diff) * _fac, + float(r_yatchew_s2_diff[r_idx]), + atol=YATCHEW_ATOL, + rtol=0, + err_msg=(f"{dgp_name}/{combo_name}/Yatchew row {r_idx}: sigma2_diff mismatch"), + ) class TestFixtureMetadata: diff --git a/tests/test_estimators_vcov_type.py b/tests/test_estimators_vcov_type.py index 42ce4dd1a..651e2e257 100644 --- a/tests/test_estimators_vcov_type.py +++ b/tests/test_estimators_vcov_type.py @@ -681,6 +681,24 @@ def _residual(df): "would be ~39, so the contrast-DOF helper may not be wired." ) + # SE-audit C2: lock the per-period BM Satterthwaite DOF (`dof_per_coef`) + # via CI-inversion. Each event-study period effect maps to R's + # `treated_period_p` coefficient; reconstructing that period's CI from + # the golden per-coef DOF and asserting equality with the reported + # `conf_int` pins the DOF (effect+se already fixed → the CI matches iff + # the DOF matches). On this fixture the post-period entries coincide + # with `dof_avg`, so this completes the per-PERIOD path coverage (a + # distinct code path from the avg-contrast DOF pinned above), not a new + # numeric target. + name_to_dof = dict(zip(d["finite_coef_names"], d["dof_per_coef"])) + for p in post_periods: + pe = res.period_effects[p] + gold_dof = name_to_dof.get(f"treated_period_{p}") + assert gold_dof is not None, f"golden dof_per_coef missing treated_period_{p}" + t_crit = float(stats.t.ppf(1.0 - res.alpha / 2.0, gold_dof)) + expected_ci = (pe.effect - t_crit * pe.se, pe.effect + t_crit * pe.se) + np.testing.assert_allclose(pe.conf_int, expected_ci, atol=1e-8, rtol=0) + def test_multi_period_fit_honors_hc2_bm(self): """MultiPeriodDiD.fit with vcov_type='hc2_bm' uses Bell-McCaffrey DOF. @@ -1608,10 +1626,16 @@ def test_absorb_hc2_bm_matches_clubsandwich_singleton_cluster(self): np.testing.assert_allclose(res.se, expected_se_slope, atol=1e-10) # ATT also bit-equal. np.testing.assert_allclose(res.att, float(d["coef"][treat_post_idx]), atol=1e-10) - # Suppress unused-local warning while keeping the constant in scope - # (DOF is exposed indirectly via res.p_value/conf_int but not as a - # standalone field on DiDResults; the SE+ATT parity above suffices). - _ = expected_dof_slope + # SE-audit C2: CI-inversion lock of `dof_hc2_bm`. res.conf_int uses a + # Satterthwaite t-dist with the BM DOF; with att+se already pinned above, + # reconstructing the CI from the golden DOF pins the DOF itself (the CI + # matches iff the estimator's DOF equals the golden DOF). Verified this + # session: matches to ~2e-16, diverging ~0.03 under the wrong n-k DOF. + from scipy import stats + + t_crit = float(stats.t.ppf(1.0 - res.alpha / 2.0, expected_dof_slope)) + expected_ci = (res.att - t_crit * res.se, res.att + t_crit * res.se) + np.testing.assert_allclose(res.conf_int, expected_ci, atol=1e-9, rtol=0) def test_unweighted_cr2_bm_per_coef_dof_no_nonphysical(self): """Unweighted clustered CR2-BM per-coef DOF: physical or NaN, never garbage. diff --git a/tests/test_estimatr_iv_robust_parity.py b/tests/test_estimatr_iv_robust_parity.py index 22ca883e9..973a0d6bf 100644 --- a/tests/test_estimatr_iv_robust_parity.py +++ b/tests/test_estimatr_iv_robust_parity.py @@ -102,6 +102,43 @@ def test_hc1_se_matches_estimatr(self, fixture, dgp_name): err_msg=f"HC1 SE mismatch on DGP {dgp_name}", ) + @pytest.mark.parametrize( + "dgp_name", + ["uniform_n200", "mild_n500", "informative_n500", "heavy_n1000"], + ) + def test_hc1_intercept_se_matches_estimatr(self, fixture, dgp_name): + """Intercept SE = sqrt(V[0, 0]) locks vs estimatr HC1 se_intercept. + + SE-audit C4: the 2SLS sandwich already forms the full 2x2 V but + only the slope diagonal drove the reported SE; surfacing the + intercept diagonal via ``return_intercept_se=True`` pins the + previously-unasserted golden ``se_intercept`` (measured ~1e-15). + Classical is intentionally NOT locked (its DOF / projection + convention diverges O(1/n); see the module docstring / REGISTRY). + """ + dgp = fixture["fixtures"][dgp_name] + d = np.asarray(dgp["d"], dtype=np.float64) + dy = np.asarray(dgp["dy"], dtype=np.float64) + w = np.asarray(dgp["w"], dtype=np.float64) + d_lower = float(dgp["d_lower"]) + + _beta, _se, _psi, se_intercept = _fit_mass_point_2sls( + d, + dy, + d_lower, + None, + "hc1", + weights=w, + return_intercept_se=True, + ) + np.testing.assert_allclose( + se_intercept, + dgp["hc1"]["se_intercept"], + atol=1e-10, + rtol=1e-10, + err_msg=f"HC1 intercept SE mismatch on DGP {dgp_name}", + ) + class TestEstimatrIVRobustCR1Parity: """CR1 pweight-cluster: bit-exact match vs estimatr (se_type='stata').""" @@ -152,6 +189,35 @@ def test_cr1_se_matches_estimatr(self, fixture): rtol=1e-10, ) + def test_cr1_intercept_se_matches_estimatr(self, fixture): + """CR1 intercept SE = sqrt(V[0, 0]) locks vs estimatr (se_type='stata'). + + SE-audit C4 (CR1 arm of the intercept-SE lock; the HC1 arm lives in + TestEstimatrIVRobustHC1Parity). Measured ~1e-15. + """ + dgp = fixture["fixtures"]["informative_cluster_n600"] + d = np.asarray(dgp["d"], dtype=np.float64) + dy = np.asarray(dgp["dy"], dtype=np.float64) + w = np.asarray(dgp["w"], dtype=np.float64) + cluster = np.asarray(dgp["cluster"], dtype=np.int64) + d_lower = float(dgp["d_lower"]) + + _beta, _se, _psi, se_intercept = _fit_mass_point_2sls( + d, + dy, + d_lower, + cluster, + "hc1", + weights=w, + return_intercept_se=True, + ) + np.testing.assert_allclose( + se_intercept, + dgp["cr1"]["se_intercept"], + atol=1e-10, + rtol=1e-10, + ) + class TestIFScaleInvariant: """IF scale invariant: compute_survey_if_variance(psi, trivial) ≈ V_HC1.""" diff --git a/tests/test_fixest_did_twfe_parity.py b/tests/test_fixest_did_twfe_parity.py index 91c53a2cd..a2a074fa5 100644 --- a/tests/test_fixest_did_twfe_parity.py +++ b/tests/test_fixest_did_twfe_parity.py @@ -89,3 +89,11 @@ def test_cluster_att_matches_fixest(self): np.testing.assert_allclose( res.att, golden[key]["cluster_unit"]["att"], atol=1e-10, rtol=0 ) + # SE-audit G2: ratio-band pin on the cluster-robust SE. The exact + # value carries the documented ~0.25% fixest-CR1 small-sample + # DOF-convention deviation (SE_AUDIT.md), so it is not machine- + # precision lockable here; this pins that we never regress BEYOND + # the known band (catches an unintended CR1 SE-formula change). The + # machine-precision hetero/cluster lock is the deferred G2 golden + # regeneration (needs an unbalanced/heteroskedastic DGP). + assert res.se == pytest.approx(golden[key]["cluster_unit"]["se"], rel=0.005) diff --git a/tests/test_methodology_placebo.py b/tests/test_methodology_placebo.py index a3d217465..c914a223b 100644 --- a/tests/test_methodology_placebo.py +++ b/tests/test_methodology_placebo.py @@ -251,6 +251,12 @@ def test_leave_one_out_matches_r(self, golden, panel): assert res.leave_one_out_effects is not None for u, att in res.leave_one_out_effects.items(): assert att == pytest.approx(gl["per_drop_att"][str(int(u))], abs=1e-10) + # SE-audit C3: pin the LOO t-test degrees of freedom (= n_valid - 1). + # The df is baked into t_stat / p_value / conf_int (asserted above) but + # is not surfaced as a result field; recover it from the public + # leave_one_out_effects (dict unit->effect, NaN for a failed drop). + n_valid = sum(1 for a in res.leave_one_out_effects.values() if np.isfinite(a)) + assert n_valid - 1 == gl["df"] def test_fake_group_matches_r(self, golden, panel): fg = golden["fake_group"]