feat(v4): EventStudyResults full-vcov + per-row df round-trip (M-092 follow-up)

CHANGELOG.md

@@ -20,6 +20,32 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   are unchanged; the renamed section fragments are kept backward-compatible by a
   client-side redirect (`docs/_static/legacy-fragment-redirect.js`) that rewrites old
   `#3.-Fit-Event-Study`-style deep links to their unnumbered targets.
+- **Unified event-study surface: full VCV persistence for StackedDiD/TwoStageDiD +
+  per-row `df` provenance (4.0 program Phase 2, ledger row M-092 follow-up).**
+  (a) `StackedDiDResults` and `TwoStageDiDResults` gain `event_study_vcov` /
+  `event_study_vcov_index`: the full event-study covariance both estimators already
+  computed internally (pooled-regression sub-block / Gardner-GMM matrix) is now
+  persisted and threaded onto `EventStudyResults.vcov`, mirroring the existing
+  CallawaySantAnna/SunAbraham fields. StackedDiD persists in every inference mode (its
+  reported ES SEs are always the matrix diagonal); TwoStageDiD persists on the
+  analytical paths and ships `None` under bootstrap (percentile inference, no
+  covariance) and replicate-weight survey designs (mixed-layout replicate VCV - CS
+  precedent). (b) `EventStudyResults.df` is now PER-ROW (float array aligned to
+  `event_time`, NaN where no df governed the stored p-value/CI; scalar inputs
+  broadcast), and joins the pinned `to_dataframe()` schema between `n` and
+  `is_reference` - an unreleased-schema change, made before `aggregate()` exposes the
+  surface publicly. Three producers genuinely need per-row df: StackedDiD `hc2_bm`
+  (per-event Bell-McCaffrey Satterthwaite df), MultiPeriodDiD `hc2_bm` (per-period BM
+  df - its `inference_df` field stores only the post-average contrast df), and LPDiD
+  (per-horizon realized cluster df). (c) New `event_study_df` provenance channels
+  record the df each event-study row's `safe_inference` actually received:
+  `CallawaySantAnnaResults` (scalar; G-1 on bare-cluster fits, the conservative
+  min-across-horizons effective df on explicit-survey fits - `df_inference` and its
+  narrow HonestDiD contract are untouched), `MultiPeriodDiDResults` (per-period dict),
+  `StackedDiDResults` (per-event dict), `TwoStageDiDResults` (scalar survey df),
+  `LPDiDResults` (per-horizon dict). All channels clear under bootstrap overrides.
+  SunAbraham/dCDH df threading is tracked in TODO.md; their surfaces show df=NaN
+  ("unexposed").
 - **ImputationDiD leave-one-out SE now anchored against Stata `did_imputation, leaveout`
   (no library behavior change).** The Borusyak-Jaravel-Spiess (2024) App. A.9 LOO variance
   (`leave_one_out=True`) has no runnable R reference (R `didimputation` omits LOO), so it
@@ -56,6 +82,22 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   **Estimator equations, weighting, variance, and numerical output are unchanged.**
 
 ### Changed
+- **DiagnosticReport's parallel-trends check upgrades from the Bonferroni fallback to
+  the joint Wald test for StackedDiD and TwoStageDiD event-study fits.** The check's
+  joint-Wald path runs whenever the source fit exposes `event_study_vcov`; now that
+  both estimators persist theirs (see Added), their pre-period tests use the full
+  covariance instead of per-coefficient Bonferroni. All-filtered TwoStageDiD horizons
+  (NaN VCV rows) are excluded from the tested family by the existing undefined-
+  inference collector, so the joint statistic is never computed through NaN entries.
+  Two accompanying guards keep the new path methodology-safe: (1) the joint-Wald
+  path now fail-closes on provenance - if any retained pre-period row carries a
+  non-finite p-value (hc2_bm Bell-McCaffrey DOF failure, collapsed replicate-survey
+  df), the check returns `inconclusive` instead of publishing a finite verdict
+  through inference the estimator refused to produce; (2) `vcov_type="hc2_bm"`
+  source fits skip the generic chi-square joint Wald entirely (it would discard the
+  CR2/BM small-sample correction; the multi-constraint AHT/HTZ CR2 test is tracked
+  in DEFERRED.md) and run Bonferroni over the BM-adjusted per-row p-values. See
+  REPORTING.md "hc2_bm parallel-trends policy".
 - **Diagnostic input validation in the report consumers (4.0 program Phase 2).**
   `BusinessReport`, `practitioner_next_steps`, and `DiagnosticReport` now route
   diagnostics by the new `Diagnostic` marker instead of by result-class name. Passing a

DEFERRED.md

@@ -48,6 +48,7 @@ exists but parity can't be verified without a local toolchain.
 | `StaggeredTripleDifference` per-cohort group-effect SEs include WIF (conservative vs R's `wif=NULL`); documented in REGISTRY. Could override the mixin for an exact R match (verification needs R `triplediff`). | `staggered_triple_diff.py` | #245 | Low |
 | **WooldridgeDiD follow-up cluster** (PR-B Stage D/E fail-closed surfaces; re-enable after R/Stata validation): <br>• QMLE sandwich uses `aweight` cluster adjustment `(G/(G-1))·(n-1)/(n-k)` vs Stata's `G/(G-1)` (conservative); add a `qmle` weight type if Stata goldens confirm a material difference (`wooldridge.py`, `linalg.py`). <br>• response-scale APE / log-link coefficient bridge for R `etwfe(family=poisson|logit)` cell-level parity — needs `emfx()` APE extraction or link-inversion with baseline-mean adjustment (`generate_wooldridge_golden.R`, `test_methodology_wooldridge.py`). <br>• `aggregate(weights="cohort_share")` on survey-weighted fits: `_n_g_per_cohort` uses raw `unit.nunique()`; implement design-weighted unit totals per cohort (paper W2025 §7) and lift the `ValueError` gate (`wooldridge.py`, `wooldridge_results.py`). <br>• unconditional inference for `cohort_share` accounting for ω̂_g sampling uncertainty (W2025 §7.5); currently NaN-closed (`wooldridge_results.py`). <br>• `cohort_trends=True × survey_design` and `× control_group="never_treated"` raise `NotImplementedError` (unvalidated TSL variance / trend columns spanned by placebo cell-dummies) (`wooldridge.py`). <br>• Stata `jwdid` golden-value `TestReferenceValues` (`tests/test_wooldridge.py`). | `wooldridge.py`, `wooldridge_results.py`, `linalg.py`, benchmarks | #216 · PR-B | Med-Low |
 | Extend `WooldridgeDiD` `method ∈ {logit, poisson}` with `vcov_type ∈ {classical, hc2, hc2_bm}`: composing HC2 leverage + Bell-McCaffrey DOF with the QMLE pseudo-residual sandwich needs derivation + R parity vs `clubSandwich::vcovCR(glm, type="CR2")`. Rejected at `__init__`. | `wooldridge.py` | follow-up | Medium |
+| Multi-constraint CR2 parallel-trends test (AHT/HTZ) for `hc2_bm` fits: DiagnosticReport's PT check routes `vcov_type="hc2_bm"` sources to Bonferroni over the BM-adjusted per-row p-values because the generic chi-square joint Wald would discard the CR2 small-sample correction (see REPORTING.md "hc2_bm parallel-trends policy"). The proper joint test is the AHT/HTZ Wald with a Satterthwaite-style denominator df over the pre-period contrast block; needs derivation for the stacked/pooled WLS-CR2 layout + parity vs `clubSandwich::Wald_test(..., test="HTZ")`. | `diagnostic_report.py`, `linalg.py` | vcov/df round-trip PR | Low |
 | `PreTrendsPower` CS/SA `anticipation=1` R-parity fixture: R `pretrends` has no anticipation parameter, so the Python `_extract_pre_period_params` anticipation filter isn't R-parity-locked. Build a synthetic CS/SA result with `anticipation=1` and assert γ_p matches R's `slope_for_power()`. (Mechanism already covered by MC + full-VCV tests.) | `tests/test_methodology_pretrends.py`, `generate_pretrends_golden.R` | PR-C | Low |
 | Harmonize SunAbraham's HC1 within-transform finite-sample correction with `fixest::sunab()` — SA applies `n/(n-k_dm)`, fixest applies `n/(n-k_total)` (counts absorbed FE); ~1-2% SE difference, documented as a "Deviation from R" and pinned at `atol=5e-3`. Either thread `df_adjustment` or keep as an intentional, R-verified difference. | `sun_abraham.py`, `linalg.py` | follow-up | Low |
 | Absorbed-FE **clustered** CR1 with *non-nested* FE: for `absorb=[FE1,FE2], cluster=FE1` (e.g. `absorb=["unit","time"], cluster="unit"`), `fixest` counts the non-nested FE (time) in the CR1 `(n-1)/(n-k)` finite-sample denominator, but the clustered path uses only `k_visible`. D4 harmonized the *non-clustered* classical/hc1 full-K scale (`_absorbed_fe_vcov_scale`) and left the clustered path unchanged — correct for FE nested in the cluster, a small deviation for non-nested FE (documented in REGISTRY within-transform note). Thread a non-nested `df_adjustment` into the clustered CR1 factor; verify vs `fixest::feols(..., cluster=)`. | `linalg.py`, `estimators.py` | SE-audit D4 | Low |

TODO.md

@@ -23,7 +23,7 @@ Related tracking surfaces:
 |-------|----------|--------|--------|----------|
 | `SyntheticControl` conformal (CWZ 2021) AR / innovation-permutation path (Lemmas 5-7) for time-series proxies — the residual-permutation shortcut is only valid for time-permutation-invariant proxies (SC/Lasso/DiD); an AR proxy needs innovation permutation. | `diff_diff/conformal.py`, `diff_diff/synthetic_control_results.py` | CWZ-2021 | Heavy | Low |
 | `ContinuousDiD` CGBS-2024 remaining extensions (earlier phases — `covariates=` reg/dr, `treatment_type="discrete"`, single-cohort `control_group="lowest_dose"` with estimand `ATT(d)−ATT(d_L)` — are already supported; see REGISTRY Note #7). Remaining (all deferred `NotImplementedError`, documented): `estimation_method="ipw"` on the dose curve (scalar-adjustment / degenerate); `covariates=` × `survey_design=` (weighted OR + weighted nuisance IF); multi-cohort **heterogeneous-support** discrete aggregation (support-aware: average each dose only over the cohorts that observe it); **multi-cohort `lowest_dose`** (within-cohort `d_L` reference + support-aware cross-cohort aggregation); and **`covariates=` × `lowest_dose`** (conditional-PT-relative-to-`d_L` estimand). Single-cohort / 2-period / shared-support multi-cohort are supported. | `continuous_did.py` | CGBS-2024 | Heavy | Low |
-| `EventStudyResults` full-vcov + df round-trip for more producers (4.0 program Phase 2 follow-up). StackedDiD and TwoStageDiD compute a full event-study VCV (regression / Gardner GMM) internally but retain only marginal SEs on their result containers, so their unified surfaces carry `vcov=None`; persist the event-time submatrix + ordered horizon index on those containers and thread it in. Likewise thread producer-specific inference df into `EventStudyResults.df` (MPD `inference_df`, explicit-survey CS `survey_metadata.df_survey`, HC2-BM horizon-specific contrast df). The adapter is faithful to what containers expose today; this is additive container work, not a defect. | `diff_diff/results_base.py`, `diff_diff/stacked_did_results.py`, `diff_diff/two_stage_results.py` | #M-092 | Mid | Low |
+| `event_study_df` provenance for the remaining producers (4.0 program Phase 2 follow-up; the Stacked/TwoStage vcov + CS/MPD/Stacked/TwoStage/LPDiD df round-trip landed with the per-row `EventStudyResults.df` schema). Remaining: SunAbraham (per-event BM Satterthwaite df under `hc2_bm` at `sun_abraham.py:1206`, scalar survey df at `:1256` - mechanism identical to StackedDiD's recorded-at-use dict); dCDH (its event-study inference receives a `df_inference` argument, `chaisemartin_dhaultfoeuille.py:6311`, never stored - needs tracing to the stored rows); LPDiD pooled-window df (popped and discarded at fit assembly). All purely additive: their surfaces show df=NaN ("unexposed") until threaded. | `diff_diff/sun_abraham.py`, `diff_diff/chaisemartin_dhaultfoeuille.py`, `diff_diff/lpdid.py` | #M-092 | Quick | Low |
 
 ### Performance
 

diff_diff/diagnostic_report.py

@@ -1398,24 +1398,72 @@ def _pt_event_study(self) -> Dict[str, Any]:
         #   2. ``event_study_vcov_index`` + ``event_study_vcov`` — the
         #      CallawaySantAnnaResults convention, where the event-study
         #      covariance is stored separately from the full regression vcov.
+        # Fail-closed inference guard for the JOINT-WALD path. The Wald
+        # statistic consumes only effects + covariance, so with a persisted
+        # vcov it would silently bypass the estimator's own fail-closed
+        # inference: a StackedDiD hc2_bm row whose Bell-McCaffrey DOF is
+        # unavailable (or a replicate-survey row whose effective df
+        # collapsed) carries a FINITE effect/SE and finite covariance
+        # entries but deliberately-NaN t/p/CI. The retained-row collector
+        # only drops non-finite effect/SE rows, and the nan-p guard below
+        # runs only on the Bonferroni fallback - so without this check a
+        # covariance-backed fit publishes a finite joint p-value for
+        # inference the estimator explicitly refused to compute. Same
+        # contract as the round-33/42 Bonferroni guards, applied BEFORE any
+        # covariance path is selected.
+        _n_nonfinite_p = sum(
+            1 for (_, _, _, p) in pre_coefs if not (isinstance(p, (int, float)) and np.isfinite(p))
+        )
+        if _n_nonfinite_p > 0:
+            return {
+                "status": "ran",
+                "method": "inconclusive",
+                "joint_p_value": None,
+                "test_statistic": None,
+                "df": len(pre_coefs),
+                "n_pre_periods": len(pre_coefs),
+                "n_dropped_undefined": _n_nonfinite_p,
+                "per_period": per_period,
+                "verdict": "inconclusive",
+                "reason": (
+                    f"{_n_nonfinite_p} retained pre-period coefficient(s) "
+                    "have non-finite per-period p-value: the source "
+                    "estimator's inference failed closed (e.g. hc2_bm "
+                    "Bell-McCaffrey DOF unavailable, or collapsed "
+                    "replicate-survey df). A joint Wald over the persisted "
+                    "covariance would silently convert that undefined "
+                    "inference into a finite verdict; the PT test is "
+                    "inconclusive on this fit."
+                ),
+            }
         vcov_for_wald: Optional[Any] = None
         idx_map_for_wald: Optional[Any] = None
         vcov_method_tag = "joint_wald"
-        if vcov is not None and interaction_indices is not None:
-            vcov_for_wald = vcov
-            idx_map_for_wald = interaction_indices
-        else:
-            es_vcov = getattr(r, "event_study_vcov", None)
-            es_vcov_index = getattr(r, "event_study_vcov_index", None)
-            if es_vcov is not None and es_vcov_index is not None:
-                vcov_for_wald = es_vcov
-                # ``event_study_vcov_index`` is an ordered list of relative-time
-                # keys; convert it into a dict mapping key -> position.
-                try:
-                    idx_map_for_wald = {k: i for i, k in enumerate(es_vcov_index)}
-                    vcov_method_tag = "joint_wald_event_study"
-                except TypeError:
-                    idx_map_for_wald = None
+        # hc2_bm-sourced fits skip the covariance-based joint Wald entirely:
+        # the generic chi-square (or F(k, df_survey)) reference ignores the
+        # CR2 / Bell-McCaffrey small-sample correction the estimator's own
+        # per-row inference applies, and the proper multi-constraint CR2
+        # test (AHT/HTZ, clubSandwich::Wald_test(test="HTZ")) is not
+        # implemented (tracked in DEFERRED.md). Bonferroni over the
+        # BM-adjusted per-row p-values keeps the small-sample correction.
+        # See REGISTRY.md StackedDiD / REPORTING.md PT notes.
+        if getattr(r, "vcov_type", None) != "hc2_bm":
+            if vcov is not None and interaction_indices is not None:
+                vcov_for_wald = vcov
+                idx_map_for_wald = interaction_indices
+            else:
+                es_vcov = getattr(r, "event_study_vcov", None)
+                es_vcov_index = getattr(r, "event_study_vcov_index", None)
+                if es_vcov is not None and es_vcov_index is not None:
+                    vcov_for_wald = es_vcov
+                    # ``event_study_vcov_index`` is an ordered list of
+                    # relative-time keys; convert it into a dict mapping
+                    # key -> position.
+                    try:
+                        idx_map_for_wald = {k: i for i, k in enumerate(es_vcov_index)}
+                        vcov_method_tag = "joint_wald_event_study"
+                    except TypeError:
+                        idx_map_for_wald = None
         df_denom: Optional[float] = None
         if vcov_for_wald is not None and idx_map_for_wald is not None and df > 0:
             try:
@@ -1425,7 +1473,33 @@ def _pt_event_study(self) -> Dict[str, Any]:
                     beta_map = {k: eff for (k, eff, _, _) in pre_coefs}
                     beta = np.array([beta_map[k] for k in keys_in_vcov], dtype=float)
                     v_sub = np.asarray(vcov_for_wald)[np.ix_(idx, idx)]
-                    stat = float(beta @ np.linalg.solve(v_sub, beta))
+                    # Rank/PSD guard. A cluster-sandwich
+                    # covariance has rank bounded by the cluster-score
+                    # dimension, so a singular pre-period block is a NORMAL
+                    # edge case whenever the number of tested leads
+                    # approaches the number of clusters. Solving through a
+                    # singular / non-PSD block produces a garbage quadratic
+                    # form (observed: stat ~ -1e15 -> chi2 p=1.0 -> a false
+                    # stakeholder-facing "no violation" verdict). Require a
+                    # finite block that is numerically full-rank and PSD
+                    # (matrix_rank's standard tolerance: k * eps * largest
+                    # singular value); otherwise raise into the Bonferroni
+                    # fallback over the (finite, per Guard above) per-row
+                    # p-values. Belt-and-braces: a Wald statistic is a PSD
+                    # quadratic form, so any non-finite or negative result
+                    # is invalid regardless of how it arose.
+                    if not np.isfinite(v_sub).all():
+                        raise ValueError("non-finite covariance block")
+                    v_sym = 0.5 * (v_sub + v_sub.T)
+                    if np.linalg.matrix_rank(v_sym) < df:
+                        raise ValueError("rank-deficient covariance block")
+                    _eigs = np.linalg.eigvalsh(v_sym)
+                    _scale = float(np.max(np.abs(_eigs))) if _eigs.size else 0.0
+                    if _scale <= 0.0 or float(_eigs.min()) < -df * np.finfo(float).eps * _scale:
+                        raise ValueError("covariance block is not positive semi-definite")
+                    stat = float(beta @ np.linalg.solve(v_sym, beta))
+                    if not np.isfinite(stat) or stat < 0.0:
+                        raise ValueError("invalid Wald statistic")
 
                     # Round-27 P1 CI review on PR #318: survey-backed
                     # fits carry a finite ``df_survey`` on
@@ -1698,8 +1772,10 @@ def _infer_cov_source(source_fit: Any) -> str:
           this fallback entirely.
         - ``"diag_fallback"`` — event-study result types with
           ``event_study_vcov is None`` (bootstrap or replicate-weight
-          CS / SA fits, plus ImputationDiD / Stacked / EfficientDiD /
-          TwoStageDiD / etc. which don't yet expose ``event_study_vcov``);
+          CS / SA / TwoStageDiD fits, plus ImputationDiD / EfficientDiD /
+          etc. which don't yet expose ``event_study_vcov``; StackedDiD
+          persists its VCV in every inference mode, so it reaches this
+          fallback only when no event study was requested);
           OR ``MultiPeriodDiDResults`` without ``interaction_indices``
           (genuine diag-only path inside ``pretrends.py:_extract_pre_period_params``,
           no "available but unused" concern, so no downgrade applies).