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

[ilgrad]

Fixes #836.

Problem

With n_bootstrap > 0, CallawaySantAnna clears event_study_vcov, both at fit time and in the post-fit aggregate("event_study") replay. HonestDiD then warns and uses diag(se**2), which drops the covariance between horizons. For the smoothness restriction this is anti-conservative in an ordinary staggered design: in the #836 simulation (parallel trends hold, known effect path), the Δ^SD(M = 0) interval for the average post-period effect covers the truth in about 68% of panels at a nominal 95%. PreTrendsPower and the DiagnosticReport parallel-trends check fall back the same way on bootstrapped CS fits.

Reproducer

main 6874b5d, Python 3.14.7, numpy 2.5.3, pandas 3.0.6, scipy 1.18.1, pure-Python backend (DIFF_DIFF_BACKEND=python). The script is the one in #836.

Script

import warnings

import numpy as np
import pandas as pd
import diff_diff as dd

# Staggered adoption with never-treated controls; parallel trends hold.
T, COHORTS = 10, (4, 6, 8)
rng = np.random.default_rng(3)
n = 300
group = rng.choice([*COHORTS, -1], size=n)
noise = np.empty((n, T))
noise[:, 0] = rng.normal(size=n) / np.sqrt(0.75)
for t in range(1, T):
    noise[:, t] = 0.5 * noise[:, t - 1] + rng.normal(size=n)
y = rng.normal(size=(n, 1)) + 0.3 * np.arange(T) + noise
for i, g in enumerate(group):
    if g >= 0:
        y[i] += [0.0 if t < g else 1.0 + 0.25 * (t - g) for t in range(T)]
df = pd.DataFrame({"unit": np.repeat(np.arange(n), T), "time": np.tile(np.arange(1, T + 1), n),
                   "y": y.ravel(), "first_treat": np.repeat(np.where(group < 0, 0, group + 1), T)})

for n_bootstrap in (0, 999):
    fit = dd.CallawaySantAnna(control_group="never_treated", base_period="universal",
                              n_bootstrap=n_bootstrap, seed=1).fit(
        df, outcome="y", unit="unit", time="time", first_treat="first_treat")
    es = fit.aggregate("event_study")
    with warnings.catch_warnings(record=True) as caught:
        warnings.simplefilter("always")
        h = dd.HonestDiD(method="smoothness", M=0.0).fit(es)
    print(f"n_bootstrap={n_bootstrap}: vcov {None if es.vcov is None else es.vcov.shape}, "
          f"Delta^SD(0) CI ({h.ci_lb:.3f}, {h.ci_ub:.3f}), width {h.ci_ub - h.ci_lb:.3f}")
    for w in caught:
        if "covariance" in str(w.message):
            print("  UserWarning:", w.message)

On main:

n_bootstrap=0: vcov (13, 13), Delta^SD(0) CI (1.520, 2.177), width 0.657
n_bootstrap=999: vcov None, Delta^SD(0) CI (1.573, 1.984), width 0.411
  UserWarning: Event-study container carries no full covariance matrix; using a diagonal approximation from the stored standard errors. Cross-event-time covariance is unavailable on this surface.

On this branch:

n_bootstrap=0: vcov (13, 13), Delta^SD(0) CI (1.520, 2.177), width 0.657
n_bootstrap=999: vcov (13, 13), Delta^SD(0) CI (1.508, 2.178), width 0.670

Cause

diff_diff/staggered.py:2956-2958 (fit time) and diff_diff/staggered_results.py:494-495 (post-fit replay) set event_study_vcov / event_study_vcov_index to None whenever the bootstrap ran. The clearing has a sound reason: the analytical influence-function VCV must not sit beside percentile-bootstrap SEs, because its diagonal would disagree with the reported SEs (_validate_vcov_subblock rejects such a container, results_base.py:842). But the bootstrap already holds a covariance that agrees with its SEs. _run_multiplier_bootstrap builds es_boot_matrix (draws × event times, staggered_bootstrap.py:679), and the event-study SEs are np.std(draws, axis=0, ddof=1) of exactly those columns (bootstrap_utils.py:466; no weighting, no IQR rescaling). So np.cov(es_boot_matrix, rowvar=False, ddof=1) has se**2 on its diagonal by construction. The clearing discarded it.

Fix

• CSBootstrapResults gains event_study_vcov and event_study_vcov_index, appended last with None defaults (repr=False on the matrix).
• _run_multiplier_bootstrap stores np.cov(es_boot_matrix[:, kept], rowvar=False, ddof=1), where kept are the horizons with a finite bootstrap SE. A universal base period's reference horizon has all-zero draws and a NaN SE, so it stays out, as it stays out of the analytical index. Nothing is stored, and consumers keep the existing diagonal fallback and its warning, in three cases:
  • a kept horizon has a non-finite draw. Its SE then comes from the finite subset (bootstrap_utils.py:505), so a full-column covariance would not match it;
  • the bootstrap variance is unidentified (fewer than 2 PSUs, the existing NaN branch);
  • n_bootstrap does not exceed the number of kept horizons. The covariance of B draws has rank at most B - 1, so the matrix would be singular by construction. The HonestDiD container route rejects a singular matrix (allow_singular=False), and the native route would use it silently.
• fit() and the post-fit replay hand this matrix on instead of clearing it. Under bootstrap the stored event-study VCV is always the bootstrap one, never the analytical one, on both routes. DMLDiDResults.aggregate("event_study") uses the same replay, so it is fixed too.
• Comments, docstrings and docs that described the clearing: the HonestDiD native-route warning (it still contains "diagonal covariance" and now names the remaining causes), the PreTrendsPower, DiagnosticReport and EventStudyResults.vcov docstrings, the SunAbraham comment that pointed at the CS clearing, REGISTRY (a new CallawaySantAnna **Note:**, the HonestDiD **Note (deviation from R):**, the PreTrendsPower Σ_22 note, the post-fit aggregate() note), REPORTING.md and guides/llms.txt.

Compatibility

• Point estimates, SEs, percentile CIs, p-values and sup-t bands do not change.
• On bootstrapped CallawaySantAnna fits (and DMLDiD replays), event_study_vcov / event_study_vcov_index and the container's vcov / vcov_index are now populated. As a result:
  • HonestDiD bounds change, and its diagonal-fallback warnings no longer fire. In the reproducer, the Δ^SD(0) width goes from 0.411 to 0.670.
  • compute_pretrends_power reports covariance_source="full_pre_period_vcov" instead of "diag_fallback", so power and MDV move.
  • The DiagnosticReport parallel-trends check runs the joint Wald test (joint_wald_event_study, or _survey on survey designs) instead of the Bonferroni fallback. This is a visible change to report output; it is listed under Behavioral Changes in the fragment.
  • Few clusters: both covariances are built from the PSU-level influence-function sums, so with fewer clusters (PSUs) than horizons the bootstrap matrix is singular, like the analytical one. HonestDiD on results.aggregate("event_study") then raises the singular-covariance ValueError it already raises for analytical fits, where main used the diagonal with a warning, and the DiagnosticReport sensitivity check reports that error. HonestDiD on the results object uses the matrix, as for analytical fits. Example: 80 units in 7 PSUs, 9 horizons, base_period="universal", n_bootstrap=199. Both matrices have rank 6. The native relative-magnitude (M = 1) interval is (1.712, 2.443), against (1.706, 2.449) analytically and (1.737, 2.417) with main's diagonal. I kept the analytical path's behavior so that the two paths agree. If you prefer a diagonal fallback for a singular bootstrap matrix, it is a small gate in _run_multiplier_bootstrap. This is also under Behavioral Changes.
• CSBootstrapResults is public. The new fields come last with defaults, so keyword and positional construction still work.
• Results pickled before this change keep event_study_vcov=None (the diagonal fallback, as before). Their aggregate("event_study") replays the bootstrap, so the replayed container now carries the matrix.
• StaggeredTripleDifference shares the bootstrap engine. Its CSBootstrapResults now also holds the matrix, but its results expose no event-study covariance, so nothing visible changes.
• Not changed: siblings whose bootstrap also clears an event-study covariance, but whose surface HonestDiD does not read. SunAbraham uses a pairs bootstrap that fills failed draws with the point estimate and computes per-horizon SEs on finite subsets, and HonestDiD does not accept SunAbraham results. TwoStageDiD has a separate multiplier engine, and its HonestDiD container admission is deferred. LWDiD is not admitted either. In all three, only PreTrendsPower (SunAbraham) and the DiagnosticReport parallel-trends check read the matrix; each would need its own derivation and tests. EfficientDiD, ImputationDiD, ContinuousDiD and HAD carry no joint event-study covariance in any mode, and dCDH is diagonal by design. I can follow up on any of these if wanted.

Methodology references (required if estimator / math changes)

• Method name(s): Callaway & Sant'Anna (2021) multiplier bootstrap for the event-study aggregation; Rambachan & Roth (2023) HonestDiD; Roth (2022) pre-trends power.
• Paper / source link(s): https://doi.org/10.1016/j.jeconom.2020.12.001, https://doi.org/10.1093/restud/rdad018, https://doi.org/10.1257/aeri.20210236
• Any intentional deviations from the source (and why): R's honest_did.AGGTEobj builds Σ from the influence-function matrix even when did reports bootstrap SEs. Here, Σ on a bootstrapped fit is the covariance of the bootstrap draws, so its diagonal is the reported SEs squared. It is the multivariate counterpart of those SEs, so it relates to the analytical Σ exactly as the bootstrap SEs relate to the analytical SEs. Documented in the REGISTRY HonestDiD **Note (deviation from R):** and the new CallawaySantAnna **Note:**.

Validation

• Tests added/updated:

  • New regression tests. All but the single-PSU one fail on main 6874b5d and pass on this branch (18 failed / 1 passed before, 19 passed after):
    • tests/test_staggered.py::TestCallawaySantAnnaBootstrap::test_bootstrap_event_study_vcov_is_draw_covariance (varying base; universal base; universal base with anticipation=1). The VCV is finite, symmetric and positive definite, ordered as its index, with diagonal equal to the bootstrap SEs squared (rtol=1e-12) and non-zero cross-horizon covariance. Reference-only horizons are excluded.
    • tests/test_staggered.py::TestCallawaySantAnnaBootstrap::test_bootstrap_event_study_vcov_needs_more_draws_than_horizons. With n_bootstrap equal to the number of horizons nothing is stored; with one more draw the matrix is stored.
    • tests/test_honest_did.py::TestSurveyVariance::test_bootstrap_vcov_same_on_native_and_container_routes. HonestDiD's Σ equals the stored bootstrap matrix on the native route. The container route gives the same Σ and the same bounds.
    • tests/test_honest_did.py::TestSurveyVariance::test_bootstrap_fit_without_vcov_falls_back_to_diagonal. Below the draw limit both routes keep the diagonal fallback and its warning.
    • tests/test_dml_did.py::TestAggregationBootstrap::test_event_study_replay_carries_bootstrap_covariance. The DMLDiD replay container carries the matrix, with diagonal equal to se**2.
    • tests/test_aggregate_contract.py::TestBootstrapReplayDesigns::test_single_psu_nan_surfaces_and_warning gains assert es.vcov is None. This behavior does not change, so it passes before and after.
  • Updated tests that pinned the clearing (the issue listed these):
    • tests/test_aggregate_contract.py::_assert_es_replay_parity: assert es.vcov is None becomes parity of the replayed vcov / vcov_index with the fit-time event_study_vcov (rtol=atol=1e-13, the helper's existing tolerance). This covers the 7 callers in TestBootstrapReplay / TestBootstrapReplayDesigns (event study, balance_e, bare cluster, stratified survey, FPC, repeated cross-sections, unbalanced panel).
    • tests/test_aggregate_contract.py::TestBootstrapReplayConsumers::test_honest_did_warns_then_uses_diagonal becomes test_honest_did_uses_bootstrap_covariance: no "no full covariance" warning, and a finite original_se.
    • tests/test_aggregate_contract.py::TestBootstrapReplayConsumers::test_pretrends_power_rides_diag_fallback becomes test_pretrends_power_uses_bootstrap_covariance: it now asserts covariance_source == "full_pre_period_vcov" (it previously asserted only power is not None).
    • tests/test_honest_did.py::TestSurveyVariance::test_bootstrap_fit_clears_analytical_vcov becomes test_bootstrap_fit_carries_bootstrap_vcov: the matrix is present, its diagonal equals se**2, and the diagonal warning no longer fires.
    • tests/test_diagnostic_report.py::TestDerivedEventStudySurface::test_bootstrap_plain_cs_derives_and_runs: the method changes from bonferroni to joint_wald_event_study, and covariance_source is full_pre_period_vcov.
  • Commands, run with the pinned tools, the pure-Python backend and DIFF_DIFF_R=skip:
    • ruff check diff_diff tests (0.16.7) and black --check diff_diff tests (26.5.1): clean.
    • mypy diff_diff (2.3.1 with numpy 2.4.5, pandas 3.0.3, scipy 1.17.1, as in lint.yml): no issues in 112 source files.
    • python .claude/scripts/changelog_compile.py check: OK. tests/test_changelog_fragments.py and tests/test_doc_deps_integrity.py: 328 passed.
    • .claude/scripts/premerge_scan.py: 0 pattern findings.
    • Every test file that mentions the CS bootstrap, event_study_vcov, HonestDiD, PreTrendsPower or DiagnosticReport (48 files): 5336 passed, 31 skipped, 0 failed.
    • The remaining files premerge_scan.py resolves, plus the TripleDifference suites (10 files: SunAbraham, StaggeredTripleDifference, TripleDifference, test_staggered_rc.py, test_pretrends_event_study.py and others): 408 passed, 119 skipped (all R-dependent), 0 failed.
    • The 19 regression node ids, rerun on the final code: 19 passed.
    • After the documentation follow-up: the changelog check again, and tests/test_changelog_fragments.py, test_doc_deps_integrity.py, test_docs_ia.py, test_doc_snippets.py, test_guides.py, test_naming_guard.py: 596 passed, 7 skipped.
  • Not run: the full suite (the 48 files above cover every consumer of the changed fields); the Rust backend (not built here); R parity tests (DIFF_DIFF_R=skip); notebooks (none changed).

• Backtest / simulation / notebook evidence (if applicable): the [Bug]: HonestDiD on a bootstrapped CallawaySantAnna event study drops the cross-horizon covariance; Δ^SD(0) covers 68% at 95% #836 coverage script, unchanged (honest_coverage.py <n_units> <n_bootstrap> <panels>). Seeds 0 to N - 1 drive both the panel and the bootstrap, and every arm uses the same seeds. Coverage of the Δ^SD(M = 0) interval for the average post-period effect, nominal 95%:

  n_units	n_bootstrap=0, main = this branch	n_bootstrap=999, main 6874b5d	n_bootstrap=999, this branch
  100	91.8% (median width 1.114)	68.4% (0.675), vcov on 0 of 1000 panels	91.2% (1.112), vcov on 1000 of 1000
  300	93.9% (0.670)	68.2% (0.393), vcov on 0 of 1000 panels	93.8% (0.668), vcov on 1000 of 1000

  1000 panels per cell, seeds 0 to 999. The Monte Carlo SE of a coverage rate is about 1.5 points near 68% and 0.8 to 0.9 points near 92-94%. The analytical arm gives identical numbers on main and on this branch. The bootstrapped arm now matches it within Monte Carlo error on the same panels; the shortfall from 95% at 100 units is in the analytical arm too.

Security / privacy

• Confirm no secrets/PII in this PR: confirmed. Code, tests and docs only; the simulations use synthetic data.

Changelog

• changelog.d/ fragment added (or N/A - no user-visible change): changelog.d/20261005-cs-bootstrap-event-study-vcov.md, with ### Fixed (bootstrapped CS fits carry the bootstrap event-study covariance, so HonestDiD bounds change) and ### Behavioral Changes (PreTrendsPower and DiagnosticReport read it; with few clusters, HonestDiD on the container raises as on analytical fits).