fix(staggered): carry the bootstrap event-study covariance instead of clearing it (#839)

* fix(staggered): carry the bootstrap event-study covariance instead of clearing it

With n_bootstrap > 0, CallawaySantAnna cleared event_study_vcov at fit time
and in the post-fit aggregate("event_study") replay, so HonestDiD fell back
to diag(se**2) and its Delta^SD intervals under-covered (#836). The clearing
kept the analytical VCV away from percentile-bootstrap SEs, but the
multiplier bootstrap already holds a consistent covariance: the event-study
SEs are np.std(draws, ddof=1) of es_boot_matrix, so the covariance of the
same draws has se**2 on its diagonal by construction.

CSBootstrapResults now carries that covariance over the horizons with a
finite bootstrap SE, and fit() and the replay (CallawaySantAnna and DMLDiD)
hand it on. It is not stored when n_bootstrap does not exceed the number of
horizons (B draws give rank <= B - 1), when a draw is non-finite, or with
fewer than two PSUs; consumers then keep the diagonal fallback. PreTrendsPower
and the DiagnosticReport parallel-trends check now read the same matrix.

* docs: state the few-cluster HonestDiD error on bootstrapped CS fits

With fewer PSUs than event-study horizons the bootstrap covariance is
singular, as the analytical one is, so HonestDiD's container route now raises
its singular-covariance ValueError where it used to fall back to the diagonal.
State it under Behavioral Changes and in the CallawaySantAnna REGISTRY note.

* fix(honest_did): validate the native-route event-study covariance like the container route

HonestDiD.fit(results) on a CallawaySantAnnaResults object assigned the
stored event_study_vcov unvalidated, while the aggregate('event_study')
container route ran _validate_vcov_subblock(..., allow_singular=False).
With fewer PSUs than horizons the stored matrix (analytical or bootstrap)
is singular, so the native route returned finite bounds where the
container route raised. Run the same check on the native route.

- tests: singular few-PSU regression on both routes (analytical and
  bootstrap); the survey-df test moves from 4 PSUs (df=2, rank-2 VCV
  over 6 horizons - singular) to 10 nested PSUs (df=8).
- docs: REGISTRY CS/HonestDiD notes, changelog fragment, and the
  practitioner guide line that still described the diagonal fallback on
  bootstrapped CS containers.

* docs(ledger): M-020 note - bootstrap replay carries the draw covariance, not vcov=None

---------

Co-authored-by: igerber <isaac.gerber@gmail.com>

Author: ilgrad

changelog.d/20261005-cs-bootstrap-event-study-vcov.md

@@ -0,0 +1,35 @@
+### Fixed
+- **Bootstrapped CallawaySantAnna fits now carry an event-study covariance**
+  (#836): with `n_bootstrap > 0`, `event_study_vcov` (fit time) and the `vcov`
+  of `aggregate("event_study")` (post-fit replay, CallawaySantAnna and DMLDiD)
+  are the covariance (`ddof=1`) of the event-study bootstrap draws, indexed by
+  the horizons with a finite bootstrap SE. Its diagonal is the bootstrap SEs
+  squared, because both come from the same draws. Previously the bootstrap
+  path cleared the covariance and HonestDiD fell back to `diag(se**2)` with a
+  warning, which made its Delta^SD intervals under-cover in staggered designs
+  with serially correlated errors. HonestDiD bounds on bootstrapped CS fits
+  therefore change, and the diagonal-fallback warnings no longer fire. The
+  matrix is still not stored, and the diagonal fallback still applies, when
+  `n_bootstrap` does not exceed the number of event-study horizons (the
+  covariance would be singular), when a draw is non-finite, or with fewer
+  than two PSUs.
+
+### Behavioral Changes
+- **Pre-trend diagnostics on bootstrapped CallawaySantAnna fits read the
+  bootstrap covariance** (#836): `compute_pretrends_power` now reports
+  `covariance_source="full_pre_period_vcov"` instead of `"diag_fallback"`,
+  and the `DiagnosticReport` parallel-trends check runs the joint Wald test
+  on the event-study covariance (`method="joint_wald_event_study"`, or its
+  `_survey` variant) instead of the Bonferroni fallback.
+- **HonestDiD on a CallawaySantAnna fit with few clusters now raises on
+  both input routes** (#836): when the clusters (PSUs) are too few for the
+  event-study horizons, the stored event-study covariance is singular
+  (analytical and bootstrap alike). HonestDiD on
+  `results.aggregate("event_study")` raised the singular-covariance
+  `ValueError` on analytical fits already and now does so on bootstrapped
+  fits too, instead of falling back to `diag(se**2)` with a warning, and the
+  `DiagnosticReport` sensitivity check reports that error. HonestDiD on the
+  native `CallawaySantAnnaResults` object now runs the same
+  `_validate_vcov_subblock(..., allow_singular=False)` check before using
+  `event_study_vcov` (previously it used the stored matrix unvalidated and
+  returned finite bounds on a singular covariance), so the two routes agree.

diff_diff/diagnostic_report.py

@@ -2206,8 +2206,10 @@ def _infer_cov_source(source_fit: Any) -> str:
           New PR-B fits set ``covariance_source`` directly and bypass
           this fallback entirely.
         - ``"diag_fallback"`` — event-study result types with
-          ``event_study_vcov is None`` (bootstrap or replicate-weight
-          CS / SA / TwoStageDiD fits, plus ImputationDiD / EfficientDiD /
+          ``event_study_vcov is None`` (replicate-weight CS fits,
+          bootstrap CS fits with ``n_bootstrap`` at or below the number
+          of event times, bootstrap or replicate-weight SA / TwoStageDiD
+          fits, plus ImputationDiD / EfficientDiD /
           etc. which don't yet expose ``event_study_vcov``; StackedDiD
           persists its VCV in every inference mode and, since 3.9
           (row M-024), always materializes the event-study surface, so

diff_diff/guides/llms-practitioner.txt

@@ -409,8 +409,9 @@ estimated pre-periods exist).
   `aggregate=` route is deprecated since 3.9). BOOTSTRAPPED fits
   (`n_bootstrap > 0`) included: the recompute levels REPLAY the
   fit-time multiplier bootstrap from the kit-retained RNG state
-  (percentile inference; the container carries no joint covariance, so
-  HonestDiD uses its diagonal approximation), while
+  (percentile inference; the container's vcov is the covariance of the
+  bootstrap draws, so HonestDiD reads a full covariance whose diagonal is
+  the bootstrap se**2), while
   `aggregate('simple')`/`aggregate('total')` relay the stored
   inference; only pre-replay legacy pickles and cross-weight-backend
   artifacts fail closed with a refit message.

diff_diff/guides/llms.txt

@@ -21,7 +21,7 @@ diagnostic steps produces unreliable results.
 4. **Choose estimator** — staggered adoption → CS/SA/BJS (NOT plain TWFE); few treated units → SDiD; factor confounding → TROP; simple 2x2 → DiD. Run `BaconDecomposition` to diagnose TWFE bias.
 5. **Estimate** — `estimator.fit(data, ...)`. Always print the cluster count first and choose inference method based on the result (cluster-robust if >= 50 clusters, wild bootstrap if fewer — for DifferenceInDifferences pass `cluster=`; TwoWayFixedEffects auto-clusters at unit level).
 6. **Sensitivity analysis** — `compute_honest_did(results)` for bounds under PT violations (MultiPeriodDiD, CS, or dCDH natively; the TwoWayFixedEffects `event_study=True` surface and a StackedDiD `results.aggregate('event_study')` container also admit - Stacked needs `kappa_pre >= 2` so estimated pre-periods exist), `run_all_placebo_tests()` for 2x2 falsification, specification comparisons for staggered designs.
-7. **Heterogeneity** — CS: `results.aggregate('group')`/`.aggregate('event_study')` post-fit, no refit (fit-time `aggregate=`/`balance_e=` are deprecated since 3.9, removed in 4.0; `compute_honest_did` / `compute_pretrends_power` / `plot_event_study` all accept the post-fit `results.aggregate('event_study')` container directly; BOOTSTRAPPED CS fits included: the recompute levels `'event_study'`/`'group'` REPLAY the fit-time multiplier bootstrap from the kit's retained RNG state (percentile inference matching a fit-time aggregation; container carries vcov=None), while `.aggregate('simple')` and, where supported, `.aggregate('total')` relay the stored bootstrap inference (NaN df column); only pre-replay legacy pickles and cross-weight-backend artifacts fail closed with a refit message). NEW in 3.10: `.aggregate('total')` on CS/EfficientDiD/ImputationDiD/TwoStageDiD - the estimator-owned total incremental outcome (exact relay C x overall over the finite-masked complete-case support; single target='total' row; bootstrap-safe RELAY; panel non-survey fits only - repeated-cross-section and declared survey_design fits raise NotImplementedError with the reason); dCDH: `results.aggregate('event_study')`/`.aggregate('simple')` post-fit views (bootstrap fits included — pure views); SA: `results.event_study_effects`/`to_dataframe(level='cohort')`; StackedDiD: `results.aggregate('event_study')`/`.aggregate('simple')` post-fit views (the surface is ALWAYS computed at fit since 3.9 - row M-024 - and the container admits into `compute_honest_did`/`compute_pretrends_power` with `kappa_pre >= 2`); EDiD: `results.aggregate('event_study')`/`.aggregate('group')`/`.aggregate('simple')` post-fit, RECOMPUTED from retained EIFs (3.9, row M-023; fit-time `aggregate=`/`balance_e=` deprecated; on bootstrapped EDiD fits the recompute levels REPLAY the fit-time multiplier bootstrap from the kit's retained RNG state (percentile inference matching a fit-time aggregation; only pre-replay legacy pickles and cross-weight-backend artifacts fail closed with a refit message) while `.aggregate('simple')` relays the stored bootstrap inference; EDiD containers are NOT admitted into honest/pretrends - no joint ES covariance); BJS/TwoStageDiD: `results.aggregate('event_study')`/`.aggregate('group')`/`.aggregate('simple')` post-fit on ImputationDiD and TwoStageDiD too (3.9, rows M-021/M-022; recomputed from panel-backed kits, `balance_e=` on `aggregate('event_study')`; on bootstrapped fits the recompute levels raise while `.aggregate('simple')` relays the stored bootstrap inference - use the deprecated fit-time aggregation for a bootstrapped ES/group surface; their containers are not admitted into honest/pretrends - Imputation by design, TwoStage deferred pending a normalization derivation); CGBS continuous: ContinuousDiD is a MIXED adopter (3.9, row M-025) - `results.aggregate('dose')` (ATT(d)+ACRT(d) rows) and `.aggregate('simple')` (att+acrt rows) are views over the always-computed curves and work on ANY fit incl. bootstrapped, while `.aggregate('event_study')` recomputes the binarized event study from a pruned per-cell IF kit and raises on bootstrapped fits (use the deprecated fit-time `aggregate='eventstudy'` there until 4.0; its container is not admitted into honest/pretrends - no joint ES covariance and no reference normalization); HAD: `results.aggregate('simple')` (overall two-period fits; the target column carries the WAS estimand label) / `.aggregate('event_study')` (multi-period fits) - pure views, work on any fit (3.9, rows M-027/M-139; fit() selects the mode from the panel shape; HAD containers are not admitted into honest/pretrends - no joint cross-horizon covariance, deferred); subgroup re-estimation.
+7. **Heterogeneity** — CS: `results.aggregate('group')`/`.aggregate('event_study')` post-fit, no refit (fit-time `aggregate=`/`balance_e=` are deprecated since 3.9, removed in 4.0; `compute_honest_did` / `compute_pretrends_power` / `plot_event_study` all accept the post-fit `results.aggregate('event_study')` container directly; BOOTSTRAPPED CS fits included: the recompute levels `'event_study'`/`'group'` REPLAY the fit-time multiplier bootstrap from the kit's retained RNG state (percentile inference matching a fit-time aggregation; the container's vcov is the covariance of the bootstrap draws, so its diagonal is the bootstrap se**2), while `.aggregate('simple')` and, where supported, `.aggregate('total')` relay the stored bootstrap inference (NaN df column); only pre-replay legacy pickles and cross-weight-backend artifacts fail closed with a refit message). NEW in 3.10: `.aggregate('total')` on CS/EfficientDiD/ImputationDiD/TwoStageDiD - the estimator-owned total incremental outcome (exact relay C x overall over the finite-masked complete-case support; single target='total' row; bootstrap-safe RELAY; panel non-survey fits only - repeated-cross-section and declared survey_design fits raise NotImplementedError with the reason); dCDH: `results.aggregate('event_study')`/`.aggregate('simple')` post-fit views (bootstrap fits included — pure views); SA: `results.event_study_effects`/`to_dataframe(level='cohort')`; StackedDiD: `results.aggregate('event_study')`/`.aggregate('simple')` post-fit views (the surface is ALWAYS computed at fit since 3.9 - row M-024 - and the container admits into `compute_honest_did`/`compute_pretrends_power` with `kappa_pre >= 2`); EDiD: `results.aggregate('event_study')`/`.aggregate('group')`/`.aggregate('simple')` post-fit, RECOMPUTED from retained EIFs (3.9, row M-023; fit-time `aggregate=`/`balance_e=` deprecated; on bootstrapped EDiD fits the recompute levels REPLAY the fit-time multiplier bootstrap from the kit's retained RNG state (percentile inference matching a fit-time aggregation; only pre-replay legacy pickles and cross-weight-backend artifacts fail closed with a refit message) while `.aggregate('simple')` relays the stored bootstrap inference; EDiD containers are NOT admitted into honest/pretrends - no joint ES covariance); BJS/TwoStageDiD: `results.aggregate('event_study')`/`.aggregate('group')`/`.aggregate('simple')` post-fit on ImputationDiD and TwoStageDiD too (3.9, rows M-021/M-022; recomputed from panel-backed kits, `balance_e=` on `aggregate('event_study')`; on bootstrapped fits the recompute levels raise while `.aggregate('simple')` relays the stored bootstrap inference - use the deprecated fit-time aggregation for a bootstrapped ES/group surface; their containers are not admitted into honest/pretrends - Imputation by design, TwoStage deferred pending a normalization derivation); CGBS continuous: ContinuousDiD is a MIXED adopter (3.9, row M-025) - `results.aggregate('dose')` (ATT(d)+ACRT(d) rows) and `.aggregate('simple')` (att+acrt rows) are views over the always-computed curves and work on ANY fit incl. bootstrapped, while `.aggregate('event_study')` recomputes the binarized event study from a pruned per-cell IF kit and raises on bootstrapped fits (use the deprecated fit-time `aggregate='eventstudy'` there until 4.0; its container is not admitted into honest/pretrends - no joint ES covariance and no reference normalization); HAD: `results.aggregate('simple')` (overall two-period fits; the target column carries the WAS estimand label) / `.aggregate('event_study')` (multi-period fits) - pure views, work on any fit (3.9, rows M-027/M-139; fit() selects the mode from the panel shape; HAD containers are not admitted into honest/pretrends - no joint cross-horizon covariance, deferred); subgroup re-estimation.
 8. **Robustness** — compare 2-3 estimators (CS vs SA vs BJS), MUST report with and without covariates (shows whether conditioning drives identification), present pre-trends and sensitivity bounds.
 
 After estimation, call `practitioner_next_steps(results)` for context-aware

diff_diff/honest_did.py

@@ -1457,9 +1457,23 @@ def _extract_event_study_params(
                         sigma = vcov
                     else:
                         sigma = np.diag(np.array(ses) ** 2)
+                    # Same integrity contract as the container route
+                    # (_extract_container_params): finite, symmetric,
+                    # diagonal == se**2, and allow_singular=False because
+                    # Rambachan-Roth inference assumes covariance eigenvalues
+                    # bounded away from zero. A few-PSU fit (fewer clusters
+                    # than horizons) stores a singular matrix on both the
+                    # analytical and the bootstrap path, so it fails loud
+                    # here exactly as results.aggregate('event_study') does.
+                    sigma = _validate_vcov_subblock(
+                        sigma, np.asarray(ses, dtype=float), "HonestDiD", allow_singular=False
+                    )
                 else:
-                    # No full VCV available. Check if this is a bootstrap fit
-                    # (VCV was cleared to prevent mixing analytical/bootstrap).
+                    # No full VCV available. A bootstrap fit carries the
+                    # covariance of its event-study draws instead of the
+                    # analytical VCV, except when n_bootstrap does not exceed
+                    # the number of event times (singular by construction) or
+                    # a draw is non-finite.
                     if (
                         hasattr(results, "bootstrap_results")
                         and results.bootstrap_results is not None
@@ -1468,9 +1482,11 @@ def _extract_event_study_params(
 
                         warnings.warn(
                             "HonestDiD on bootstrap-fitted CallawaySantAnna results "
-                            "uses a diagonal covariance matrix (cross-event-time "
-                            "covariance is not available from bootstrap). For full "
-                            "covariance structure, use analytical SEs (n_bootstrap=0).",
+                            "uses a diagonal covariance matrix (no covariance of the "
+                            "event-study bootstrap draws is available: n_bootstrap "
+                            "does not exceed the number of event times, or a draw is "
+                            "non-finite). For full covariance structure, increase "
+                            "n_bootstrap or use analytical SEs (n_bootstrap=0).",
                             UserWarning,
                             stacklevel=4,
                         )

diff_diff/pretrends.py

@@ -185,12 +185,14 @@ def _extract_event_study_vcov_subblock(
 
     This is the canonical Σ_22 routing path for ``compute_pretrends_power``
     when the event-study result type exposes a full event-study covariance
-    matrix (CallawaySantAnnaResults non-bootstrap fits at
-    ``staggered_results.py:126-128`` and SunAbrahamResults non-bootstrap
-    fits via the W-matrix construction added in PR-B Step 3). Bootstrap
-    fits and replicate-weight survey fits clear ``event_study_vcov`` so
-    the analytical VCV is not mixed with bootstrap / replicate SE
-    overrides — those cases naturally fall through to the diag fallback.
+    matrix (CallawaySantAnnaResults fits at ``staggered_results.py:126-128``
+    and SunAbrahamResults non-bootstrap fits via the W-matrix construction
+    added in PR-B Step 3). A bootstrapped CallawaySantAnna fit carries the
+    covariance of its event-study bootstrap draws, whose diagonal is the
+    bootstrap SEs squared. Bootstrap SunAbraham fits and replicate-weight
+    survey fits clear ``event_study_vcov`` so the analytical VCV is not
+    mixed with bootstrap / replicate SE overrides — those cases naturally
+    fall through to the diag fallback.
 
     Parameters
     ----------
@@ -216,7 +218,9 @@ def _extract_event_study_vcov_subblock(
         ``"full_pre_period_vcov"`` when the full event-study sub-block
         was used (no off-diagonal information was discarded), or
         ``"diag_fallback"`` when ``event_study_vcov`` was missing /
-        cleared (bootstrap / replicate-weight CS or SA paths).
+        cleared (replicate-weight CS or SA paths, bootstrap SA paths, and
+        bootstrap CS fits with ``n_bootstrap`` at or below the number of
+        event times).
     """
     es_vcov = getattr(results, "event_study_vcov", None)
     es_vcov_index = getattr(results, "event_study_vcov_index", None)
@@ -1310,9 +1314,9 @@ def _extract_pre_period_params(
                 ses = np.array([pre_effects[t]["se"] for t in pre_periods])
 
                 # Route through full event_study_vcov when available
-                # (non-bootstrap CS fits at staggered_results.py:126-128).
-                # Bootstrap CS fits clear event_study_vcov at
-                # staggered.py:2032-2036, falling through to diag.
+                # (CS fits at staggered_results.py:126-128; bootstrap CS
+                # fits carry the covariance of their bootstrap draws).
+                # Replicate-weight CS fits clear it, falling through to diag.
                 vcov, covariance_source = _extract_event_study_vcov_subblock(
                     results, pre_periods, ses
                 )

diff_diff/results_base.py

@@ -281,7 +281,9 @@ class EventStudyResults(BaseResults):
         ordered by ``vcov_index``. None when the producer records no matrix
         or when the stored SEs are no longer its diagonal (bootstrap and
         replicate-weight overrides clear it rather than ship an
-        inconsistent matrix).
+        inconsistent matrix). A bootstrapped CallawaySantAnna (or DMLDiD)
+        surface instead carries the covariance of its event-study bootstrap
+        draws, whose diagonal is the stored bootstrap SEs squared.
     vcov_index : np.ndarray or None
         ``event_time`` labels labelling ``vcov``'s rows/columns (explicit
         ordering for HonestDiD / PreTrendsPower consumption).

diff_diff/staggered.py

@@ -2949,13 +2949,14 @@ def fit(
 
         # Store results
         # Retrieve event-study VCV from aggregation mixin (Phase 7d).
-        # Clear it when bootstrap overwrites event-study SEs to prevent
-        # HonestDiD from mixing analytical VCV with bootstrap SEs.
+        # When bootstrap overwrites the event-study SEs, the VCV is the
+        # covariance of the same bootstrap draws (diagonal = bootstrap SE**2),
+        # so HonestDiD never mixes the analytical VCV with bootstrap SEs.
         event_study_vcov = es_aggregation.vcov if es_aggregation is not None else None
         event_study_vcov_index = es_aggregation.vcov_index if es_aggregation is not None else None
-        if bootstrap_results is not None and event_study_vcov is not None:
-            event_study_vcov = None
-            event_study_vcov_index = None
+        if bootstrap_results is not None:
+            event_study_vcov = bootstrap_results.event_study_vcov
+            event_study_vcov_index = bootstrap_results.event_study_vcov_index
         # Same clearing rule for the ES df provenance: bootstrap replaces the
         # stored ES se/p/CI with percentile values that never used the
         # analytical df, so surfacing it would be false provenance.