From c436dc2326627edd4318b13d78bcd32672e9ed21 Mon Sep 17 00:00:00 2001 From: igerber Date: Sat, 4 Jul 2026 08:32:02 -0400 Subject: [PATCH] feat(had): cluster-robust event-study inference + clustered sup-t band HeterogeneousAdoptionDiD.fit(aggregate="event_study") now honors cluster= end-to-end on both designs: cluster-robust per-horizon pointwise CIs (continuous CCT + mass-point 2SLS) AND a cluster-robust simultaneous sup-t confidence band. Previously cluster= was ignored on the nonparametric event-study path (with a UserWarning) and the weighted mass-point cband=True case raised NotImplementedError. The clustered sup-t band adds a dedicated branch to _sup_t_multiplier_bootstrap: it aggregates the per-unit influence function to cluster level and draws cluster-level Rademacher multipliers, so the perturbation variance is the raw cluster sandwich sum_c s_c^2. This reconciles to each path's analytical cluster-robust SE via a path scalar: 1.0 for continuous (the lprobust cluster meat carries no g/(g-1) correction) and sqrt(G/(G-1)) for mass-point (restoring the CR1 finite-sample factor the returned IF omits). Validated bootstrap-free: sqrt(sum_c (scale*s_c)^2) == se to atol=1e-10 on the real IF for both paths, plus the H=1 -> 1.96 reduction on the clustered branch. cluster= + survey= is rejected (both designs) - route clustering through survey_design=SurveyDesign(psu=). Single cluster -> NaN band + RuntimeWarning. Unclustered fits are byte-unchanged. Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_01LHDijzf8zHXk5T8ahS2mKi --- CHANGELOG.md | 12 ++ TODO.md | 3 +- diff_diff/guides/llms-full.txt | 14 +- diff_diff/had.py | 361 +++++++++++++++++++++----------- docs/api/had.rst | 54 ++--- docs/methodology/REGISTRY.md | 16 +- tests/test_had.py | 363 +++++++++++++++++++++++++++++---- 7 files changed, 621 insertions(+), 202 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index b3bba70b0..9fdf8cb3f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,18 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +### Added +- **`HeterogeneousAdoptionDiD` event-study `cluster=` support** (both designs). On + `aggregate="event_study"`, `cluster=` now produces cluster-robust per-horizon pointwise + confidence intervals AND a cluster-robust simultaneous sup-t confidence band (`cband=True`), + on the continuous (CCT local-linear) and mass-point (2SLS) paths, unweighted or weighted — + previously `cluster=` was ignored on the nonparametric event-study path (with a `UserWarning`) + and the weighted mass-point `cband=True` case raised `NotImplementedError`. The clustered band + draws cluster-level multipliers on the per-unit influence function; the variance family is + reconciled to each path's analytical cluster-robust SE (exact for continuous; `√(G/(G-1))` CR1 + scaling for mass-point). `cluster=` + `survey=` is rejected — route clustering through + `survey_design=SurveyDesign(psu=)`. No behavior change for unclustered fits. + ### Fixed - **Dube, Girardi, Jordà & Taylor (2025) citation corrected to *J. Applied Econometrics* 40(**7**):741-758** (was cited as issue 5) across `docs/references.rst`, diff --git a/TODO.md b/TODO.md index f78938569..0dc210409 100644 --- a/TODO.md +++ b/TODO.md @@ -34,7 +34,6 @@ The `Origin` column (Actionable tables) and the `PR` column (Deferred tables) bo | `TwoWayFixedEffects(vcov_type in {hc2, hc2_bm})` with replicate-weight designs raises `NotImplementedError` (`twfe.py:~233`). The replicate path re-demeans per replicate, which doesn't compose with the full-dummy HC2/HC2-BM build — a correct impl needs per-replicate full-dummy refit. Workaround: `hc1` for replicate-weight CR1. | `twfe.py::fit` | follow-up | Heavy | Low | | TWFE's HC2/HC2-BM inline full-dummy build (`twfe.py:280-315`) duplicates the dummy-construction logic in `DifferenceInDifferences(fixed_effects=...)` (`estimators.py:478-486`). Extract a shared helper, or delegate TWFE's HC2/HC2-BM path to DiD's `fixed_effects=` branch (with TWFE-specific cluster-default threading), to reduce drift risk on FE naming / survey behavior / result-surface conventions. Substantive refactor — touches both estimators. | `twfe.py::fit`, `estimators.py::DifferenceInDifferences.fit` | follow-up | Heavy | Low | | Decide whether to formally deprecate `CallawaySantAnna.cluster=X` in favor of `survey_design=SurveyDesign(psu=X)` (the bare-cluster path already synthesizes a minimal SurveyDesign). Two equivalent paths = redundant surface. Mirrors the question for ImputationDiD / EfficientDiD / TwoStageDiD. | `staggered.py` | follow-up | Mid | Low | -| `HeterogeneousAdoptionDiD` **event-study (Phase 2b)** continuous cluster= threading: Phase 2a static path now threads `cluster=` into `bias_corrected_local_linear` (cluster-robust CCT SE, unweighted + weighted). The per-horizon event-study path still ignores `cluster=` with a `UserWarning` because the `cband` sup-t bootstrap normalizes HC-scale perturbations by the analytical SE and would mix variance families under clustering (mirrors the mass-point `weights= + cluster= + cband=True` `NotImplementedError`). Needs a per-horizon clustered-bootstrap variance-family reconciliation. | `had.py::_fit_event_study` | Phase 2b | Mid | Low | ### Performance @@ -125,7 +124,7 @@ Doable in principle, but no current caller and/or explicitly out of paper scope. | SunAbraham-family fits are SOLVER-bound after the v3.6.x demeaning speedups: on the county-class shape (3.1k units x 60 months, ~100 interaction columns) `solve_ols` = ~60% of fit time (faer SVD 1.24s + pivoted-QR rank detection 0.70s of 3.27s), and the share grows further under the Rust demean kernel; the pivoted QR alone is ~20% of the whole fit (measured attribution in docs/performance-plan.md). Candidate fix is QR-reuse (one factorization for rank detection + solve). **Parked per the 2026-07 correctness-first decision**: the change perturbs every coefficient in the library at the last digits (golden/parity/bit-identity recapture) and removes the SVD's independent singular-value truncation - a second, different near-collinearity detector that is deliberate defense-in-depth for the no-silent-failures contract (the FE-span guard episode showed borderline cases are real). Re-open only on practitioner demand for wide event studies; any implementation gates on fixest/R end-to-end parity, not just internal consistency. | `linalg.py::solve_ols`, `linalg.py::_detect_rank_deficiency` | PR-C attribution | Medium | | dCDH parity-test SE/CI assertions only cover pure-direction scenarios; mixed-direction SE comparison is structurally apples-to-oranges (cell-count vs obs-count weighting). | `test_chaisemartin_dhaultfoeuille_parity.py` | #294 | Low | | `HeterogeneousAdoptionDiD` survey-design API consolidation (**scheduled: next minor bump**): drop the deprecated `survey=` / `weights=` kwargs on all 8 HAD surfaces; only `survey_design=` remains. Also fold the legacy back-end `weights=` routing into the unified `_resolve_survey_for_fit` path. DeprecationWarning has shipped; the removal is ~50 LoC gated on the semver bump. | `had.py`, `had_pretests.py` | next minor | Medium | -| `HeterogeneousAdoptionDiD` joint cross-horizon covariance / sup-t bands: per-horizon SEs use independent sandwiches (paper-faithful pointwise CIs per Pierce-Schott Fig 2). Follow-ups (low demand): IF-based stacking for joint cross-horizon inference; analytical H×H covariance on the weighted ES path; a sup-t band on the unweighted ES path. | `had.py::_fit_event_study` | Phase 2b / 4.5 B | Low | +| `HeterogeneousAdoptionDiD` joint cross-horizon covariance / sup-t bands: per-horizon SEs use independent sandwiches (paper-faithful pointwise CIs per Pierce-Schott Fig 2). Follow-ups (low demand): IF-based stacking for joint cross-horizon inference; analytical H×H covariance on the weighted ES path. (A sup-t band on the unweighted ES path shipped via the clustered band — `cluster=` fires the simultaneous band on the unweighted path.) | `had.py::_fit_event_study` | Phase 2b / 4.5 B | Low | | `HeterogeneousAdoptionDiD` event-study staggered timing beyond the last cohort: Phase 2b auto-filters to the last cohort (paper App B.2); earlier-cohort effects aren't HAD-identified (redirect to dCDH). Full staggered HAD needs a different identification path (out of paper scope). | `had.py::_validate_had_panel_event_study` | Phase 2b | Low | | `HeterogeneousAdoptionDiD` survey-aware support-endpoint test (**research, waits on literature**): needs a calibrated support-infimum test under complex sampling (endpoint EVT × survey-aware functional CLT × tail-empirical-process theory). Permanent `NotImplementedError` on `qug_test(survey=...)`; rationale in REGISTRY § "QUG Null Test" Note (Phase 4.5 C0). | `had_pretests.py::qug_test` | Phase 4.5 C0 | Low | | `HeterogeneousAdoptionDiD` Phase-4.5 weight-aware auto-bandwidth MSE-DPI selector (~300 LoC); users pass `h`/`b` explicitly today. Plus replicate-weight SurveyDesigns on the continuous-dose paths (Rao-Wu-style per-replicate weight-ratio rescaling for the local-linear intercept IF). | `_nprobust_port.py::lpbwselect_mse_dpi`, `had.py::_aggregate_unit_resolved_survey` | Phase 4.5 | Low | diff --git a/diff_diff/guides/llms-full.txt b/diff_diff/guides/llms-full.txt index 8da9890c7..44d2e0b01 100644 --- a/diff_diff/guides/llms-full.txt +++ b/diff_diff/guides/llms-full.txt @@ -738,8 +738,8 @@ HeterogeneousAdoptionDiD( alpha: float = 0.05, vcov_type: str | None = None, # Mass-point only: "classical" (default) or "hc1" robust: bool = False, # Mass-point only: HC1 robust SE shorthand - cluster: str | None = None, # Mass-point only: cluster column for CR1 cluster-robust SE - n_bootstrap: int = 999, # Multiplier-bootstrap iterations for sup-t bands (event-study + weighted) + cluster: str | None = None, # Cluster column for cluster-robust SEs (all designs; event-study adds a cluster-robust sup-t band). cluster= + survey_design= rejected + n_bootstrap: int = 999, # Multiplier-bootstrap iterations for sup-t bands (event-study, when weighted/survey OR clustered) seed: int | None = None, ) ``` @@ -759,7 +759,7 @@ had.fit( aggregate: str = "overall", # "overall" (single scalar WAS) or "event_study" (per-horizon WAS) survey: SurveyDesign | None = None, # DEPRECATED alias of survey_design= weights: np.ndarray | None = None, # DEPRECATED pweight shortcut alias - cband: bool = True, # Simultaneous (sup-t) confidence bands on weighted event-study fits + cband: bool = True, # Simultaneous (sup-t) confidence bands on event-study fits that are weighted/survey OR clustered *, survey_design: SurveyDesign | None = None, # Canonical survey-design kwarg (weights, strata, PSU, FPC) trends_lin: bool = False, # Eq 17 linear-trend detrending. Requires aggregate="event_study"; needs F>=3 (pre-period depth) for the regression; rejects ALL weighting entry paths (survey_design= / survey= / weights= all raise NotImplementedError under trends_lin). @@ -799,7 +799,7 @@ es = est.fit(data_mp, outcome_col='y', unit_col='unit', **Staggered panels.** On multi-cohort panels with `aggregate="event_study"`, `fit()` auto-filters to the last treatment cohort plus never-treated units (paper Appendix B.2) and emits a `UserWarning` naming kept/dropped counts. The estimand is then a **last-cohort-only WAS**, not a multi-cohort average. For full multi-cohort staggered support, see `ChaisemartinDHaultfoeuille`. -**Mass-point + survey constraint.** When fitting `design="mass_point"` with `survey_design=` (or the deprecated `survey=` alias), `vcov_type="hc1"` (or `robust=True`) is required: the survey path composes the standard error via Binder-TSL on the HC1-scale influence function, so the default classical sandwich path raises `NotImplementedError`. The same HC1 requirement also fires on the `weights=` shortcut when `aggregate="event_study"` AND `cband=True`: the per-horizon IF matrix is HC1-scale and the sup-t bootstrap normalizes by it, so mixing in a classical analytical SE would produce inconsistent variance families. Classical vcov is allowed on `weights=` + `aggregate="overall"` and on `weights=` + `aggregate="event_study"` + `cband=False`. Passing `vcov_type="hc1"` is a safe default on weighted survey + sup-t examples since `vcov_type` is unused on the continuous designs (CCT-2014 robust SE is the only formula there). +**Mass-point + survey constraint.** When fitting `design="mass_point"` with `survey_design=` (or the deprecated `survey=` alias), `vcov_type="hc1"` (or `robust=True`) is required: the survey path composes the standard error via Binder-TSL on the HC1-scale influence function, so the default classical sandwich path raises `NotImplementedError`. The same HC1 requirement also fires on the `weights=` shortcut when `aggregate="event_study"` AND `cband=True` — UNLESS `cluster=` is set (a clustered mass-point fit resolves to the CR1 sandwich regardless of `vcov_type` and uses the clustered sup-t band, so there is no classical-vs-HC1 mismatch): the per-horizon IF matrix is HC1-scale and the sup-t bootstrap normalizes by it, so mixing in a classical analytical SE would produce inconsistent variance families. Classical vcov is allowed on `weights=` + `aggregate="overall"` and on `weights=` + `aggregate="event_study"` + `cband=False`. Passing `vcov_type="hc1"` is a safe default on weighted survey + sup-t examples since `vcov_type` is unused on the continuous designs (CCT-2014 robust SE is the only formula there). ### StackedDiD @@ -1545,7 +1545,7 @@ Per-horizon event-study results container for `HeterogeneousAdoptionDiD` with `a | `F` | `object` | First-treatment period label | | `n_units` | `int` | Unique units contributing to the fit (post last-cohort filter) | | `inference_method` | `str` | `"analytical_nonparametric"` or `"analytical_2sls"` | -| `vcov_type` | `str | None` | Mass-point only: `"classical"`, `"hc1"`, or `"cr1"`; `None` on continuous designs | +| `vcov_type` | `str | None` | Mass-point: `"classical"`, `"hc1"`, or `"cr1"` (with `cluster=`). Continuous: `None`, or `"cr1"` with `cluster=` | | `cluster_name` | `str | None` | Cluster column name when CR1 is requested; `None` otherwise | | `survey_metadata` | `SurveyMetadata | None` | Populated on weighted fits | | `bandwidth_diagnostics` | `list[BandwidthResult | None] | None` | Per-horizon MSE-DPI selector output (continuous designs); `None` on `mass_point`; entries can be `None` on degenerate horizons | @@ -1553,10 +1553,10 @@ Per-horizon event-study results container for `HeterogeneousAdoptionDiD` with `a | `filter_info` | `dict | None` | Staggered last-cohort auto-filter metadata (`F_last`, `n_kept`, `n_dropped`, `dropped_cohorts`); `None` when no filter applied | | `variance_formula` | `str | None` | HAD-specific SE label applied UNIFORMLY across all horizons, populated on BOTH continuous and mass-point designs: `"pweight"` (continuous, CCT 2014 weighted-robust on the `weights=` shortcut), `"survey_binder_tsl"` (continuous, Binder 1983 TSL on the `survey_design=` path), `"pweight_2sls"` (mass-point + `weights=`; label applied uniformly across vcov families — classical / HC1 / CR1 — on the weighted 2SLS path, with the actual sandwich resolved via `vcov_type`), or `"survey_binder_tsl_2sls"` (mass-point, Binder 1983 TSL on the `survey_design=` path). `None` on unweighted fits | | `effective_dose_mean` | `float | None` | Weighted denominator used by the β̂-scale rescaling, populated on weighted fits across all designs: weighted `sum(w·d)/sum(w)` (`continuous_at_zero`), weighted `sum(w·(d − d_lower))/sum(w)` (`continuous_near_d_lower`), or weighted Wald-IV dose gap (`mass_point`). Scalar (not per-horizon) because the β̂-scale denominator is computed once on the fit sample. `None` on unweighted fits | -| `cband_low` | `np.ndarray | None` | Simultaneous (sup-t) band lower bounds; `None` on unweighted fits or when `cband=False` | +| `cband_low` | `np.ndarray | None` | Simultaneous (sup-t) band lower bounds; `None` when `cband=False` or on unweighted, unclustered fits (a clustered fit produces the band even when unweighted) | | `cband_high` | `np.ndarray | None` | Simultaneous (sup-t) band upper bounds | | `cband_crit_value` | `float | None` | Sup-t critical value used for the simultaneous band | -| `cband_method` | `str | None` | `"multiplier_bootstrap"` when populated | +| `cband_method` | `str | None` | `"multiplier_bootstrap"` (weighted/survey band) or `"cluster_multiplier_bootstrap"` (clustered band) when populated | | `cband_n_bootstrap` | `int | None` | Bootstrap iterations used for the band | **Methods:** `summary()`, `print_summary()`, `to_dict()`, `to_dataframe()` diff --git a/diff_diff/had.py b/diff_diff/had.py index 8e2b6e08b..e1a525240 100644 --- a/diff_diff/had.py +++ b/diff_diff/had.py @@ -604,9 +604,10 @@ class HeterogeneousAdoptionDiDEventStudyResults: t-inference. Pointwise CIs are always populated; a simultaneous confidence - band is available only on the weighted path via ``cband_*`` - below. Joint cross-horizon analytical covariance is not - computed in this release (tracked in TODO.md). + band (``cband_*`` below) is available on the weighted/survey path + OR whenever ``cluster=`` is set (the clustered band fires even on + an unweighted fit). Joint cross-horizon analytical covariance is + not computed in this release (tracked in TODO.md). t_stat, p_value : np.ndarray, shape (n_horizons,) Per-horizon inference triple element. conf_int_low, conf_int_high : np.ndarray, shape (n_horizons,) @@ -645,9 +646,10 @@ class HeterogeneousAdoptionDiDEventStudyResults: ``"analytical_nonparametric"`` (continuous designs) or ``"analytical_2sls"`` (mass-point). Shared across horizons. vcov_type : str or None - Effective variance-covariance family used on the mass-point path - (``"classical"``, ``"hc1"``, or ``"cr1"`` when cluster supplied). - ``None`` on the continuous paths (they use CCT-2014 robust SE). + Effective variance-covariance family. Mass-point path: + ``"classical"``, ``"hc1"``, or ``"cr1"`` when ``cluster=`` is + supplied. Continuous paths: ``None`` (CCT-2014 robust SE), or + ``"cr1"`` when ``cluster=`` is supplied (cluster-robust CCT SE). cluster_name : str or None Column name of the cluster variable when cluster-robust SE is requested. ``None`` otherwise. @@ -670,10 +672,11 @@ class HeterogeneousAdoptionDiDEventStudyResults: fits. cband_low, cband_high : np.ndarray or None, shape (n_horizons,) Simultaneous confidence-band endpoints constructed by the - multiplier-bootstrap sup-t procedure. ``None`` on unweighted - fits and when ``fit(..., cband=False)`` is passed. Horizons - with ``se <= 0`` or non-finite ``se`` are NaN (matches the - pointwise inference gate from ``safe_inference``). + multiplier-bootstrap sup-t procedure. ``None`` when + ``fit(..., cband=False)`` is passed or on unweighted, unclustered + fits; a clustered fit (``cluster=``) populates the band even when + unweighted. Horizons with ``se <= 0`` or non-finite ``se`` are NaN + (matches the pointwise inference gate from ``safe_inference``). cband_crit_value : float or None Sup-t multiplier-bootstrap critical value at level ``1 - alpha``. Under a trivial resolved design (no strata / @@ -684,8 +687,9 @@ class HeterogeneousAdoptionDiDEventStudyResults: variance matches the analytical Binder-TSL target term-for- term. cband_method : str or None - ``"multiplier_bootstrap"`` on the weighted event-study path - with ``cband=True``, else ``None``. + ``"multiplier_bootstrap"`` on the weighted/survey event-study band, + ``"cluster_multiplier_bootstrap"`` on the clustered band, else + ``None``. cband_n_bootstrap : int or None Number of multiplier-bootstrap replicates used to compute the sup-t critical value. @@ -739,7 +743,9 @@ class HeterogeneousAdoptionDiDEventStudyResults: filter_info: Optional[Dict[str, Any]] # Phase 4.5 B weighted / survey-path extras (optional so unweighted - # fits stay unchanged; all None on unweighted fits). + # fits stay unchanged). ``variance_formula`` / ``effective_dose_mean`` + # are None on unweighted fits; the ``cband_*`` fields are also populated + # on unweighted CLUSTERED fits (Phase 2b clustered band). variance_formula: Optional[str] = None """Per-horizon variance family label (applied uniformly across all horizons in the fit). One of ``"pweight"`` / ``"pweight_2sls"`` (when @@ -756,22 +762,25 @@ class HeterogeneousAdoptionDiDEventStudyResults: mass-point: weighted Wald-IV dose gap. ``None`` on unweighted fits.""" cband_low: Optional[np.ndarray] = None """Simultaneous confidence-band lower endpoints, shape ``(n_horizons,)``. - ``None`` on unweighted fits and when ``cband=False`` on the weighted - event-study path. Derived from multiplier-bootstrap sup-t critical - value: ``cband_low[e] = att[e] − cband_crit_value * se[e]``.""" + ``None`` when ``cband=False``, or on unweighted, unclustered fits (a + clustered fit populates the band even when unweighted). Derived from the + multiplier-bootstrap sup-t critical value: + ``cband_low[e] = att[e] − cband_crit_value * se[e]``.""" cband_high: Optional[np.ndarray] = None """Simultaneous confidence-band upper endpoints, shape ``(n_horizons,)``. See ``cband_low``.""" cband_crit_value: Optional[float] = None """Sup-t multiplier-bootstrap critical value at level ``1 - alpha``. Reduces to ``Φ⁻¹(1 − alpha/2) ≈ 1.96`` at ``H=1`` up to Monte Carlo - error. ``None`` on unweighted fits and when ``cband=False``.""" + error. ``None`` when ``cband=False`` or on unweighted, unclustered fits.""" cband_method: Optional[str] = None - """``"multiplier_bootstrap"`` on the weighted event-study path with - ``cband=True``, else ``None``.""" + """``"multiplier_bootstrap"`` (weighted/survey band) or + ``"cluster_multiplier_bootstrap"`` (clustered band) when populated, else + ``None``.""" cband_n_bootstrap: Optional[int] = None """Number of multiplier-bootstrap replicates used to compute the sup-t - critical value. ``None`` on unweighted fits and when ``cband=False``.""" + critical value. ``None`` when ``cband=False`` or on unweighted, + unclustered fits.""" def __repr__(self) -> str: base = ( @@ -901,10 +910,13 @@ def to_dict(self) -> Dict[str, Any]: "vcov_type": self.vcov_type, "cluster_name": self.cluster_name, "filter_info": _json_safe_filter_info(self.filter_info), - # Phase 4.5 B weighted/survey-path surfaces (None on - # unweighted fits). The full SurveyMetadata dataclass is - # carried as an object, matching the static-path ``to_dict`` - # contract — consumers read attributes uniformly. + # Phase 4.5 B weighted/survey-path surfaces. survey_metadata / + # variance_formula / effective_dose_mean are None on unweighted + # fits; the cband_* fields are also populated on unweighted + # CLUSTERED fits (Phase 2b clustered band). The full + # SurveyMetadata dataclass is carried as an object, matching the + # static-path ``to_dict`` contract — consumers read attributes + # uniformly. "survey_metadata": self.survey_metadata, "variance_formula": self.variance_formula, "effective_dose_mean": self.effective_dose_mean, @@ -919,9 +931,10 @@ def to_dataframe(self) -> pd.DataFrame: """Return a tidy per-horizon DataFrame. Columns: ``event_time, att, se, t_stat, p_value, conf_int_low, - conf_int_high, n_obs``. One row per event-time horizon. On the - weighted event-study path with ``cband=True``, also includes - ``cband_low`` and ``cband_high`` columns. + conf_int_high, n_obs``. One row per event-time horizon. When a + simultaneous band was computed (``cband=True`` with a + weighted/survey OR clustered fit), also includes ``cband_low`` and + ``cband_high`` columns. """ data: Dict[str, Any] = { "event_time": self.event_times, @@ -2043,6 +2056,8 @@ def _sup_t_multiplier_bootstrap( alpha: float, seed: Optional[int], bootstrap_weights: str = "rademacher", + cluster_ids: Optional[np.ndarray] = None, + cluster_if_scale: float = 1.0, ) -> Tuple[float, Optional[np.ndarray], Optional[np.ndarray], int]: """Compute sup-t simultaneous CI via PSU-level multiplier bootstrap. @@ -2094,6 +2109,24 @@ def _sup_t_multiplier_bootstrap( bootstrap_weights : str Passed through to the helper: ``"rademacher"``, ``"mammen"``, or ``"webb"``. Default ``"rademacher"`` (binary ±1 multipliers). + cluster_ids : np.ndarray or None, shape (n_units,) + Per-unit cluster labels aligned to ``influence_matrix`` rows. When + supplied, a dedicated CLUSTER-robust branch fires (bypassing both + the survey and unit-level branches): the per-unit IF is aggregated + to cluster level and cluster-level iid multipliers are drawn, so + ``Var_xi(xi @ Psi_cl) = sum_c (sum_{i in c} psi_i)^2`` — the RAW + cluster sandwich matching the analytical cluster-robust SE (no + stratum-centering / FPC / Bessel; the continuous CCT cluster meat + carries no ``g/(g-1)`` correction). ``resolved_survey`` must be + ``None`` when this is set (``cluster= + survey=`` is rejected + up-front). ``None`` restores the survey / unit-level dispatch. + cluster_if_scale : float + Path-specific finite-sample scalar applied to the cluster-aggregated + IF before the draw: ``1.0`` for the CCT continuous path (exact by + construction), ``sqrt(G/(G-1))`` for the mass-point CR1 path (the + returned mass-point IF carries ``sqrt((n-1)/(n-k))`` but not the CR1 + ``G/(G-1)`` factor, so it is restored here). ``G`` is the full-array + cluster count, identical to the ``n_clusters`` this branch computes. Returns ------- @@ -2124,7 +2157,32 @@ def _sup_t_multiplier_bootstrap( # reserved for the ``weights=`` shortcut (no survey object at all). use_survey_bootstrap = resolved_survey is not None - if use_survey_bootstrap: + if cluster_ids is not None: + # Cluster-robust multiplier bootstrap (had.py clustered event-study + # band). Aggregate the per-unit IF to cluster level and draw + # cluster-level iid multipliers, so + # Var_xi(xi @ Psi_cl) = sum_c (sum_{i in c} psi_i)^2 + # — the RAW cluster sandwich, matching the analytical cluster-robust + # SE. No stratum-centering / FPC / Bessel: the continuous CCT cluster + # meat (lprobust_vce) carries no g/(g-1) correction, so raw Rademacher + # is exact (scale 1.0); the mass-point CR1 path restores its G/(G-1) + # factor via ``cluster_if_scale`` before the draw. ``G`` == this + # branch's ``n_clusters`` == the count the analytical CR uses (full + # array, incl. wholly-zero-weight clusters that contribute 0). + cluster_ids_arr = np.asarray(cluster_ids).ravel() + codes, cluster_labels = pd.factorize(cluster_ids_arr, sort=False) + n_clusters = len(cluster_labels) + if n_clusters < 2: + # Cluster-robust simultaneous band undefined with one cluster. + return float("nan"), None, None, 0 + Psi_cl = np.zeros((n_clusters, n_horizons), dtype=np.float64) + np.add.at(Psi_cl, codes, influence_matrix) + Psi_cl *= float(cluster_if_scale) + perturbations = np.empty((n_bootstrap, n_horizons), dtype=np.float64) + for _cs, _w_block in iter_weight_blocks(n_bootstrap, n_clusters, bootstrap_weights, rng): + with np.errstate(divide="ignore", invalid="ignore", over="ignore"): + perturbations[_cs : _cs + _w_block.shape[0]] = _w_block @ Psi_cl + elif use_survey_bootstrap: # Review R2 P1: lonely_psu="adjust" pools singleton strata into a # pseudo-stratum with NONZERO multipliers in the bootstrap helper, # but the analytical compute_survey_if_variance target for @@ -2622,8 +2680,12 @@ class HeterogeneousAdoptionDiD: composition raises ``NotImplementedError`` (the Binder-TSL survey variance would override the cluster-robust SE — route clustering through ``survey_design=SurveyDesign(psu=)`` instead). - Cluster must be constant within unit. Estimator-level cluster - threading on the event-study path (Phase 2b) remains a follow-up. + Cluster must be constant within unit. On the event-study path + (``aggregate="event_study"``, Phase 2b) ``cluster=`` provides + cluster-robust per-horizon pointwise CIs (both designs) AND a + cluster-robust simultaneous sup-t band (``cband=True``, fires even + on unweighted fits); ``cluster=`` + ``survey_design=`` is rejected + there too. Notes ----- @@ -2729,14 +2791,14 @@ def __init__( self.vcov_type = vcov_type self.robust = robust self.cluster = cluster - # Phase 4.5 B: event-study survey sup-t simultaneous-CI support. - # ``n_bootstrap`` = number of multiplier-bootstrap replicates for - # the sup-t band on the event-study + weighted path. ``seed`` = - # reproducibility seed for the multiplier draws. Both are - # consulted only when ``aggregate="event_study"`` AND a - # ``survey=`` / ``weights=`` is passed to ``fit()`` with - # ``cband=True`` (default). Unweighted event-study skips the - # bootstrap entirely — pre-Phase 4.5 B numerical output preserved. + # Event-study sup-t simultaneous-CI support. ``n_bootstrap`` = + # number of multiplier-bootstrap replicates for the sup-t band; + # ``seed`` = reproducibility seed for the multiplier draws. Both are + # consulted when ``aggregate="event_study"`` + ``cband=True`` + # (default) AND either ``survey=`` / ``weights=`` (weighted/survey + # band, Phase 4.5 B) OR ``cluster=`` (cluster-robust band, Phase 2b + # — fires even unweighted). An unweighted, unclustered event-study + # skips the bootstrap entirely (pre-Phase 4.5 B output preserved). self.n_bootstrap = n_bootstrap self.seed = seed self._validate_constructor_args() @@ -2978,14 +3040,16 @@ def fit( with the per-row → per-unit aggregation invariant intact. Mutually exclusive with ``survey_design=`` and ``survey=``. cband : bool, default True - Phase 4.5 B: controls the multiplier-bootstrap simultaneous - confidence band on the weighted event-study path. When - ``True`` (default) and ``aggregate="event_study"`` AND any of - ``survey_design=`` / ``survey=`` / ``weights=`` is supplied, - the fit populates ``cband_low`` / ``cband_high`` / - ``cband_crit_value`` / ``cband_method`` / ``cband_n_bootstrap`` - on the result. When ``False`` those fields stay ``None``. No - effect on ``aggregate="overall"`` or on unweighted event- + Controls the multiplier-bootstrap simultaneous confidence band + on the event-study path. When ``True`` (default) and + ``aggregate="event_study"`` AND either (a) ``survey_design=`` / + ``survey=`` / ``weights=`` is supplied (weighted/survey band, + Phase 4.5 B) OR (b) ``cluster=`` is supplied (cluster-robust + band — fires even on an unweighted fit, Phase 2b), the fit + populates ``cband_low`` / ``cband_high`` / ``cband_crit_value`` + / ``cband_method`` / ``cband_n_bootstrap`` on the result. When + ``False`` those fields stay ``None``. No effect on + ``aggregate="overall"`` or on unweighted, unclustered event- study. ``n_bootstrap`` and ``seed`` (constructor params) control replicate count and RNG; defaults are 999 / ``None``. trends_lin : bool, default False, keyword-only @@ -4054,17 +4118,21 @@ def _fit_event_study( (``_fit_mass_point_2sls`` HC1 / classical / CR1). Inference is Normal (``df=None``). - The simultaneous confidence band on the weighted path (when - ``cband=True``) is constructed by a shared-PSU multiplier - bootstrap over the stacked per-horizon β̂-scale IF matrix via - :func:`_sup_t_multiplier_bootstrap`. On the ``weights=`` - shortcut, sup-t calibration is routed through a synthetic - trivial ``ResolvedSurveyDesign`` so the centered + + The simultaneous confidence band (when ``cband=True``) is + constructed by a multiplier bootstrap over the stacked per-horizon + β̂-scale IF matrix via :func:`_sup_t_multiplier_bootstrap`. On the + weighted/survey path it draws PSU-level multipliers; when + ``cluster=`` is set it takes the clustered branch (cluster-level + multipliers, raw cluster sandwich; fires even unweighted). On the + ``weights=`` shortcut, sup-t calibration is routed through a + synthetic trivial ``ResolvedSurveyDesign`` so the centered + sqrt(n/(n-1))-corrected survey-aware branch fires uniformly — matches the analytical HC1 variance family at the - compute_survey_if_variance(IF, trivial) ≈ V_HC1 invariant. - Unweighted event-study skips the bootstrap (pre-Phase 4.5 B - numerical output preserved). + compute_survey_if_variance(IF, trivial) ≈ V_HC1 invariant. When + ``cluster=`` is set the clustered branch fires instead (cluster- + level multipliers, raw cluster sandwich), even unweighted. An + unweighted, unclustered event-study skips the bootstrap (pre-Phase + 4.5 B numerical output preserved). """ # ---- Resolve effective fit-time state (local vars only, # feedback_fit_does_not_mutate_config). ---- @@ -4348,46 +4416,33 @@ def _fit_event_study( dose_mean = float(d_arr.mean()) - # ---- Extract cluster IDs on mass-point path only ---- + # ---- Extract cluster IDs (mass-point + continuous paths) ---- + # cluster= threads into the per-horizon CR1 sandwich (mass-point) or + # the cluster-robust CCT SE (continuous, Phase 2a parity), AND into a + # cluster-robust sup-t simultaneous band (both designs) via the + # clustered branch of ``_sup_t_multiplier_bootstrap``. The one + # incompatible composition is cluster= + survey= (either design): the + # survey path composes Binder-TSL variance and would silently override + # the cluster-robust sandwich while metadata still reports the + # cluster-robust vcov. Reject it BEFORE extracting the column (mirrors + # the static-path guard had.py:3361) so the error is predictable even + # for a malformed cluster column. The former weights= + cluster= + + # cband=True mass-point rejection is lifted — the clustered band now + # reconciles the variance family (raw cluster Rademacher; mass-point + # sqrt(G/(G-1)) CR1 scaling). cluster_arr: Optional[np.ndarray] = None - if resolved_design == "mass_point" and cluster_arg is not None: - # Review R4 P1: narrow the cluster+weighted guard (mirrors - # the static-path narrowing). Incompatible cases on the - # event-study path: - # (a) survey= + cluster=: Binder-TSL override would - # silently overwrite CR1. - # (b) weights= shortcut + cluster= + cband=True: the - # sup-t bootstrap normalizes HC1-scale perturbations - # by the CR1 analytical SE, producing an inconsistent - # variance family in the bootstrap t-distribution. - # weights= shortcut + cluster= + cband=False is fine: the - # per-horizon CR1 sandwich is returned as-is and no IF is - # consumed. Unweighted + cluster= also unchanged. - if resolved_survey_unit_full is not None: - raise NotImplementedError( - f"cluster={cluster_arg!r} + survey= on " - f"design='mass_point' event-study is not yet " - f"supported: the survey path composes Binder-TSL " - f"variance per horizon and would silently override " - f"the CR1 cluster-robust sandwich. Pass cluster= " - f"alone (unweighted CR1), or weights= + cluster= " - f"+ cband=False (weighted-CR1 per horizon), or " - f"survey= alone (Binder-TSL). Combined cluster-" - f"robust + survey event-study inference is deferred." - ) - if weights_unit_full is not None and cband: - raise NotImplementedError( - f"cluster={cluster_arg!r} + weights= + cband=True " - f"on design='mass_point' event-study is not yet " - f"supported: the sup-t bootstrap uses an HC1-scale " - f"influence function and normalizes by the CR1 " - f"analytical SE, mixing variance families in the " - f"bootstrap t-distribution. Pass cband=False to " - f"disable the simultaneous band (pointwise CIs " - f"still use the weighted-CR1 sandwich per horizon), " - f"or drop cluster= to use the weighted-HC1 sandwich " - f"with sup-t." - ) + if cluster_arg is not None and resolved_survey_unit_full is not None: + raise NotImplementedError( + f"cluster={cluster_arg!r} + survey= on " + f"design={resolved_design!r} event-study is not supported: " + f"the survey path composes Binder-TSL variance per horizon " + f"and would silently override the cluster-robust sandwich. " + f"Pass cluster= alone (unweighted cluster-robust), weights= " + f"+ cluster= (weighted cluster-robust, pointwise + sup-t " + f"band), or survey_design=SurveyDesign(psu=) to " + f"cluster through the survey (Binder-TSL) path." + ) + if cluster_arg is not None: _, _, cluster_arr, _, _ = _aggregate_multi_period_first_differences( data_filtered, outcome_col, @@ -4427,15 +4482,9 @@ def _fit_event_study( UserWarning, stacklevel=3, ) - if cluster_arg is not None: - warnings.warn( - f"cluster={cluster_arg!r} is ignored on the " - f"'{resolved_design}' path in Phase 2b (estimator-" - f"level cluster threading on the nonparametric path " - f"is queued for a follow-up PR).", - UserWarning, - stacklevel=3, - ) + # cluster= is now threaded into the continuous-path CCT SE and the + # cluster-robust sup-t band (Phase 2b parity with the static path); + # no "cluster ignored" warning. # ---- Resolve vcov label for mass-point ---- if resolved_design == "mass_point": @@ -4455,8 +4504,13 @@ def _fit_event_study( # give wrong t-stats). Matches the static-path rejection — # weighted mass-point paths use the HC1-scale IF # convention uniformly. - _uses_if_matrix = resolved_survey_unit_full is not None or ( - weights_unit_full is not None and cband + # When cluster= is set, the mass-point path computes the CR1 + # cluster-robust sandwich regardless of vcov_type (classical/hc1 + # are ignored), and the clustered sup-t band normalizes by that + # CR1 SE — so there is no classical-vs-HC1 variance-family + # mismatch to reject. Guard on ``cluster_arg is None``. + _uses_if_matrix = cluster_arg is None and ( + resolved_survey_unit_full is not None or (weights_unit_full is not None and cband) ) if vcov_requested == "classical" and _uses_if_matrix: raise NotImplementedError( @@ -4480,8 +4534,11 @@ def _fit_event_study( else: vcov_requested = "" inference_method = "analytical_nonparametric" - vcov_label = None - cluster_label = None + # cluster-robust CCT SE when cluster= is set (Phase 2a static-path + # parity, had.py:3615); labelled "cr1" for surface consistency + # with the mass-point CR1 path. + vcov_label = "cr1" if cluster_arg is not None else None + cluster_label = cluster_arg if cluster_arg is not None else None # ---- Per-horizon loop ---- # On the weighted path, every horizon uses the FULL arrays @@ -4516,12 +4573,16 @@ def _fit_event_study( # compute_survey_if_variance inside _fit_continuous or the # mass-point override below); the STACKED (G, H) IF matrix is # needed only when the sup-t multiplier bootstrap runs - # (``cband=True`` on the weighted path). Splitting them avoids + # (``cband=True`` on a weighted/survey OR clustered fit). Splitting + # them avoids # allocating / filling Psi on the common opt-out path # ``cband=False`` + weights= shortcut, where no IF consumer # exists. + # The clustered sup-t band consumes the stacked IF too, and (unlike + # the survey/weights band) fires even on the UNWEIGHTED path — so the + # stacked-IF flag also switches on when cluster= is set. needs_per_horizon_if = resolved_survey_unit_full is not None or (weighted_es and cband) - needs_stacked_if_matrix = weighted_es and cband + needs_stacked_if_matrix = cband and (weighted_es or cluster_arg is not None) if needs_stacked_if_matrix: Psi = np.full((G_full, n_horizons), np.nan, dtype=np.float64) else: @@ -4542,13 +4603,14 @@ def _fit_event_study( if resolved_survey_unit_full is not None: df_infer = resolved_survey_unit_full.df_survey - # On the weighted event-study path, the sup-t multiplier bootstrap - # operates on the per-horizon IF matrix, so we must force the IF - # computation even on the ``weights=`` shortcut (no survey - # structure → _fit_continuous normally skips IF). Pass through - # the actual ``resolved_survey_unit_full`` (None on shortcut) so + # Whenever the sup-t multiplier bootstrap runs (weighted/survey OR + # clustered — see ``needs_stacked_if_matrix``), it operates on the + # per-horizon IF matrix, so we must force the IF computation even + # when _fit_continuous would normally skip it (``weights=`` shortcut + # or unweighted clustered fit; no survey structure). Pass through + # the actual ``resolved_survey_unit_full`` (None on those paths) so # the per-horizon analytical SE still matches the static-path - # convention (bc_fit.se_robust on shortcut; Binder-TSL on + # convention (bc_fit.se_robust on shortcut/clustered; Binder-TSL on # survey=). IF return is gated on `force_return_influence=True`. # Track the Binder-TSL den for continuous paths so we can @@ -4574,6 +4636,9 @@ def _fit_event_study( force_return_influence=( needs_stacked_if_matrix and resolved_survey_unit_full is None ), + # Phase 2b: cluster-robust CCT SE per horizon (static-path + # parity). None on unclustered fits (byte-unchanged). + cluster_arr=cluster_arr, ) if bc_fits is not None: bc_fits.append(bc_fit_e) @@ -4605,12 +4670,12 @@ def _fit_event_study( cluster_arr, vcov_requested, weights=weights_unit_full, - # Return IF only when a consumer exists: survey= - # path needs it for per-horizon Binder-TSL override; - # weights= shortcut + cband=True needs it for the - # bootstrap. weights= shortcut + cband=False skips - # IF computation entirely (review R6 P2). - return_influence=needs_per_horizon_if, + # Return IF when a consumer exists: survey= path needs it + # for per-horizon Binder-TSL override; the sup-t bootstrap + # (weighted OR clustered) needs it for the stacked matrix. + # weights= shortcut + cband=False + unclustered skips IF + # computation entirely (review R6 P2). + return_influence=(needs_per_horizon_if or needs_stacked_if_matrix), ) # Survey path: override analytical sandwich SE with # Binder-TSL via compute_survey_if_variance (matches @@ -4638,13 +4703,63 @@ def _fit_event_study( ci_lo_arr[i] = float(conf_int_e[0]) ci_hi_arr[i] = float(conf_int_e[1]) - # ---- Sup-t simultaneous confidence band (weighted + cband only) ---- + # ---- Sup-t simultaneous confidence band (weighted/survey OR + # clustered, + cband) ---- cband_low_arr: Optional[np.ndarray] = None cband_high_arr: Optional[np.ndarray] = None cband_crit_value: Optional[float] = None cband_method_label: Optional[str] = None cband_n_bootstrap_eff: Optional[int] = None - if weighted_es and cband and n_horizons >= 1: + if cband and cluster_arg is not None and n_horizons >= 1: + # Cluster-robust simultaneous band (Phase 2b). Route cluster= as + # cluster-level multipliers on the per-unit IF: the raw cluster + # Rademacher variance sum_c (sum_{i in c} psi_i)^2 equals the + # analytical cluster-robust SE² for the continuous CCT path (no + # g/(g-1) correction; scale 1.0) and, for the mass-point CR1 path, + # matches after restoring the CR1 sqrt(G/(G-1)) factor (the + # returned mass-point IF already carries sqrt((n-1)/(n-k)) but not + # G/(G-1)). G is the full-array cluster count — identical to the + # bootstrap branch's own n_clusters, so scale and aggregation + # share one count (a wholly-zero-weight cluster contributes 0 to + # both yet is counted in G by both). Fires on the UNWEIGHTED path + # too (cluster= + survey= is rejected up-front, so no survey + # composition). See REGISTRY "Note (HAD clustered event-study + # sup-t band)". + _cl_G = int(len(pd.unique(np.asarray(cluster_arr)))) + if _cl_G < 2: + warnings.warn( + f"cluster={cluster_arg!r} yields a single cluster; the " + f"cluster-robust sup-t band and pointwise SEs are " + f"undefined (NaN band returned).", + RuntimeWarning, + stacklevel=3, + ) + # Undefined band, NOT skipped: mirror the degenerate + # weighted-band convention (crit=NaN, method + replicate + # count populated, endpoints left None) so callers can tell + # "attempted but undefined" apart from "no band requested". + cband_crit_value = float("nan") + cband_method_label = "cluster_multiplier_bootstrap" + cband_n_bootstrap_eff = n_bootstrap_eff + else: + _cl_scale = ( + float(np.sqrt(_cl_G / (_cl_G - 1))) if resolved_design == "mass_point" else 1.0 + ) + q, cband_low_arr, cband_high_arr, _n_valid = _sup_t_multiplier_bootstrap( + influence_matrix=Psi, + att_per_horizon=att_arr, + se_per_horizon=se_arr, + resolved_survey=None, + n_bootstrap=n_bootstrap_eff, + alpha=float(self.alpha), + seed=seed_eff, + cluster_ids=cluster_arr, + cluster_if_scale=_cl_scale, + ) + cband_crit_value = q + cband_method_label = "cluster_multiplier_bootstrap" + cband_n_bootstrap_eff = n_bootstrap_eff + elif weighted_es and cband and n_horizons >= 1: # Review R7 P0: the per-unit influence function returned by # _fit_continuous / _fit_mass_point_2sls is HC1-scaled per # the PR #359 convention — compute_survey_if_variance(psi, diff --git a/docs/api/had.rst b/docs/api/had.rst index 3e9f08c9c..4070aa0cc 100644 --- a/docs/api/had.rst +++ b/docs/api/had.rst @@ -49,10 +49,9 @@ Unit Remains Untreated" (arXiv:2405.04465v6), which: - **``weights=np.ndarray`` shortcut (deprecated)** - continuous paths reuse the CCT-2014 SE; the mass-point path uses an analytical weighted 2SLS sandwich (``classical`` / ``hc1``; CR1 when - ``cluster=`` is supplied, except ``cluster=`` + - ``aggregate="event_study"`` + ``cband=True`` is rejected outright - regardless of ``vcov_type`` per the cluster-combination deviation - below; ``hc2`` / ``hc2_bm`` raise ``NotImplementedError`` pending a + ``cluster=`` is supplied - including ``aggregate="event_study"`` + + ``cband=True``, which adds a cluster-robust simultaneous band; + ``hc2`` / ``hc2_bm`` raise ``NotImplementedError`` pending a 2SLS-specific leverage derivation). Yields ``variance_formula="pweight"`` / ``"pweight_2sls"``. - **``survey_design=SurveyDesign(weights="col", ...)``** (canonical; @@ -79,16 +78,22 @@ Unit Remains Untreated" (arXiv:2405.04465v6), which: workflow ``did_had_pretest_workflow`` handles this by skipping QUG under survey/weighted dispatch and emitting a ``UserWarning``. - A simultaneous confidence band (sup-t) is available only on the - **weighted event-study path** via ``cband=True``. Joint cross-horizon - analytical covariance is not computed in this release; tracked in - ``TODO.md``. + A simultaneous confidence band (sup-t) is available on the event-study + path via ``cband=True`` whenever ``weights=``/``survey_design=`` **or** + ``cluster=`` is supplied. With ``cluster=`` the band is cluster-robust + (pointwise CIs are cluster-robust too), on both designs and unweighted + or weighted; ``cluster=`` + ``survey_design=`` is rejected (route + clustering through ``survey_design=SurveyDesign(psu=)``). + Joint cross-horizon analytical covariance is not computed in this + release; tracked in ``TODO.md``. **Mass-point ``vcov_type="classical"`` deviation.** The mass-point ``survey_design=SurveyDesign(...)`` paths (static and event-study) and the deprecated ``weights=`` + ``aggregate="event_study"`` + ``cband=True`` path reject ``vcov_type="classical"`` with - ``NotImplementedError``. The per-unit 2SLS influence function returned + ``NotImplementedError`` **when ``cluster=`` is not set** (a ``cluster=`` + fit computes the CR1 sandwich regardless of ``vcov_type``, so no + classical/HC1 mismatch arises). The per-unit 2SLS influence function returned by the mass-point fit is HC1-scaled so that ``compute_survey_if_variance`` and the sup-t bootstrap target ``V_HC1`` consistently; mixing it with a classical analytical SE @@ -98,23 +103,20 @@ Unit Remains Untreated" (arXiv:2405.04465v6), which: ``vcov_type="classical"``, which triggers the guard); a classical-aligned IF derivation is queued for a follow-up PR. - **Mass-point cluster-combination deviation.** On - ``design="mass_point"``, two clustered weighted paths are rejected - outright regardless of ``vcov_type``: - - - ``survey_design=SurveyDesign(...)`` + ``cluster=`` (static and - event-study): the survey path composes Binder-TSL variance, which - would silently override the CR1 cluster-robust sandwich. - Workarounds: ``cluster=`` alone (unweighted CR1), or ``weights=`` - + ``cluster=`` (weighted-CR1 pweight sandwich), or - ``survey_design=`` alone (Binder-TSL). Combined cluster-robust + - survey inference is queued for a follow-up PR. - - Deprecated ``weights=`` shortcut + ``cluster=`` + - ``aggregate="event_study"`` + ``cband=True``: the sup-t bootstrap - normalizes HC1-scale perturbations by the CR1 analytical SE, - mixing variance families. Workarounds: pass ``cband=False`` (keeps - weighted-CR1 per-horizon), or drop ``cluster=`` (keeps - weighted-HC1 sup-t). + **``cluster=`` + ``survey_design=`` deviation.** On both designs, + ``cluster=`` + ``survey_design=SurveyDesign(...)`` (static and + event-study) is rejected outright regardless of ``vcov_type``: the + survey path composes Binder-TSL variance, which would silently + override the cluster-robust sandwich. Workarounds: ``cluster=`` alone + (unweighted cluster-robust), ``weights=`` + ``cluster=`` (weighted + cluster-robust, pointwise + simultaneous band), or route clustering + through ``survey_design=SurveyDesign(psu=)``. All other + ``cluster=`` compositions are supported end-to-end, including the + ``weights=`` + ``cluster=`` + ``aggregate="event_study"`` + + ``cband=True`` mass-point path (formerly rejected): the clustered + sup-t band draws cluster-level multipliers on the per-unit influence + function and normalizes by the CR1 analytical SE (variance families + reconciled via the ``√(G/(G-1))`` CR1 scalar). .. tip:: diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index c7c665a51..cf45ea420 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -3003,7 +3003,7 @@ Under `survey=SurveyDesign(weights, strata, psu, fpc)`, the variance composes vi *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). -*Sup-t multiplier bootstrap (Phase 4.5 B):* Simultaneous confidence band on the weighted event-study path via `_sup_t_multiplier_bootstrap`: +*Sup-t multiplier bootstrap (Phase 4.5 B):* Simultaneous confidence band on the weighted/survey OR clustered event-study path via `_sup_t_multiplier_bootstrap` (the clustered branch is documented in "Note (HAD clustered event-study sup-t band)"): 1. **Multiplier draws**: reuse `diff_diff.bootstrap_utils.generate_survey_multiplier_weights_batch` (survey= path: PSU-level draws with stratum centering, FPC scaling, lonely-PSU handling) or `generate_bootstrap_weights_batch` (weights= shortcut: unit-level Rademacher). Default `n_bootstrap=999` (CS parity); `seed` exposed on `HeterogeneousAdoptionDiD.__init__` for reproducibility. 2. **Perturbations**: `delta = xi @ Psi` — shape `(B, H)` matrix-matrix product, NO `(1/n)` prefactor (matches `staggered_bootstrap.py:373` idiom; `Psi` is already on the β̂-scale). @@ -3011,18 +3011,19 @@ Under `survey=SurveyDesign(weights, strata, psu, fpc)`, the variance composes vi 4. **Sup-t distribution**: `sup_t[b] = max_e |t[b, e]|` with finite-mask filtering of degenerate horizons. 5. **Critical value**: `q = quantile(sup_t[finite], 1 - alpha)`. Simultaneous band: `cband_low[e] = att[e] - q · se[e]`. -**Reduction invariant**: at `H=1`, the sup collapses to the marginal and `q → Φ⁻¹(1 - alpha/2) ≈ 1.96` at `alpha=0.05` up to MC noise. Locked by `TestSupTReducesToNormalAtH1` (G=500, B=5000, seed=42, `atol=0.15` on the quantile) and `TestEventStudySurveyCband::test_trivial_survey_h1_sup_t_matches_analytical` / `test_stratified_h1_sup_t_matches_analytical` for the trivial-survey and stratified cases respectively. +**Reduction invariant**: at `H=1`, the sup collapses to the marginal and `q → Φ⁻¹(1 - alpha/2) ≈ 1.96` at `alpha=0.05` up to MC noise. Locked by `TestSupTReducesToNormalAtH1` (G=500, B=5000, seed=42, `atol=0.15` on the quantile), its clustered variant `test_clustered_sup_t_h1_reduces_to_normal` (both the continuous scale-1.0 and mass-point `√(G/(G-1))` scalars), and `TestEventStudySurveyCband::test_trivial_survey_h1_sup_t_matches_analytical` / `test_stratified_h1_sup_t_matches_analytical` for the trivial-survey and stratified cases respectively. -**Scope**: sup-t bootstrap runs only when `aggregate="event_study"` AND `weights=` or `survey=` is supplied AND `cband=True` (default). Unweighted event-study skips the bootstrap entirely — pre-Phase 4.5 B numerical output bit-exactly preserved. Setting `cband=False` on the weighted path disables the bootstrap (useful for smoke-test bit-parity assertions against the unweighted path at uniform weights). +**Scope**: sup-t bootstrap runs when `aggregate="event_study"` AND `cband=True` (default) AND either (a) `weights=`/`survey=` is supplied (the survey/weights band) OR (b) `cluster=` is supplied (the clustered band — see "Note (HAD clustered event-study sup-t band)"). An unweighted, unclustered event-study skips the bootstrap entirely — pre-Phase 4.5 B numerical output bit-exactly preserved. Setting `cband=False` disables the bootstrap on any path. **`weights=` shortcut ↔ bootstrap routing**: on the `weights=` shortcut (no user-supplied `SurveyDesign`), per-horizon SE stays analytical (CCT-2014 robust for continuous, HC1/classical/CR1 sandwich for mass-point — NOT Binder-TSL), but sup-t calibration is routed through a synthetic trivial `ResolvedSurveyDesign` (pweight, no strata / PSU / FPC, `lonely_psu="remove"`) so the centered + sqrt(n/(n-1))-corrected bootstrap branch fires uniformly. This matches the analytical HC1 variance family — the `compute_survey_if_variance(IF, trivial) ≈ V_HC1` invariant from the IF scale convention — so the bootstrap variance target agrees with the per-horizon SE normalizer. Without this routing, the unit-level bootstrap branch would normalize against raw `sum(ψ²) = ((n-1)/n) · V_HC1` on the HC1-scaled IF and produce silently too-narrow simultaneous bands. Regression-locked by the `weights=`-shortcut / `survey=`-trivial equivalence test in `TestEventStudySurveyCband`. - **Deviation from shared survey-bootstrap contract:** `_sup_t_multiplier_bootstrap` raises `NotImplementedError` on `SurveyDesign(lonely_psu="adjust")` with singleton strata. The shared `generate_survey_multiplier_weights_batch` helper pools singleton PSUs into a pseudo-stratum with NONZERO multipliers, but `compute_survey_if_variance` centers singleton PSU scores at the GLOBAL mean of PSU scores (rather than the pseudo-stratum mean). Matching the two would require a pooled-singleton pseudo-stratum centering transform in the HAD sup-t path that has not been derived. The HAD-specific limitation is scoped to: weighted event-study + `cband=True` + `lonely_psu="adjust"` + at least one singleton stratum. Practitioners can use `lonely_psu="remove"` or `"certainty"` (matches the analytical target bit-exactly on the HAD sup-t path), or pass `cband=False` to skip the simultaneous band. All other survey-bootstrap consumers (CallawaySantAnna, dCDH, SDID) retain full `lonely_psu="adjust"` support through the shared helper. -- **Deviation: weighted mass-point `vcov_type="classical"` on survey/sup-t paths:** `vcov_type="classical"` raises `NotImplementedError` whenever the mass-point IF matrix is consumed downstream — specifically on `design="mass_point"` + `survey=` (static + event-study) and `design="mass_point"` + `weights=` + `aggregate="event_study"` + `cband=True`. The per-unit 2SLS IF returned by `_fit_mass_point_2sls` is scaled (`sqrt((n-1)/(n-k))`) to match V_HC1 via `compute_survey_if_variance`; mixing it with a classical analytical SE would silently return a V_HC1-targeted variance under a classical label. A classical-aligned IF derivation is queued for a follow-up PR. The allowed weighted-mass-point combinations are: `vcov_type="hc1"` on every path; `vcov_type="classical"` on `weights=` + `aggregate="overall"`, and `weights=` + `aggregate="event_study"` + `cband=False` (no IF consumption). +- **Deviation: weighted mass-point `vcov_type="classical"` on survey/sup-t paths:** `vcov_type="classical"` raises `NotImplementedError` whenever the mass-point IF matrix is consumed downstream — specifically on `design="mass_point"` + `survey=` (static + event-study) and `design="mass_point"` + `weights=` + `aggregate="event_study"` + `cband=True` — **only when `cluster=` is NOT set** (with `cluster=` the mass-point path computes the CR1 sandwich regardless of `vcov_type`, so no classical-vs-HC1 mismatch exists and the classical-rejection is guarded by `cluster_arg is None`). The per-unit 2SLS IF returned by `_fit_mass_point_2sls` is scaled (`sqrt((n-1)/(n-k))`) to match V_HC1 via `compute_survey_if_variance`; mixing it with a classical analytical SE would silently return a V_HC1-targeted variance under a classical label. A classical-aligned IF derivation is queued for a follow-up PR. The allowed weighted-mass-point combinations are: `vcov_type="hc1"` on every path; `vcov_type="classical"` on `weights=` + `aggregate="overall"`, and `weights=` + `aggregate="event_study"` + `cband=False` (no IF consumption); and any `cluster=` composition (resolves to CR1). -- **Deviation: mass-point `cluster=` + `survey=` rejected:** the `survey=` path composes Binder-TSL variance via `compute_survey_if_variance`, which would silently overwrite the CR1 cluster-robust sandwich while result metadata still reports `vcov_type="cr1"`. Narrowed (R4) to apply only to: `design="mass_point"` + `survey=` + `cluster=` (static and event-study), and `design="mass_point"` + `weights=` + `cluster=` + `aggregate="event_study"` + `cband=True` (where the sup-t bootstrap normalizes HC1-scale perturbations by the CR1 analytical SE). The `weights=` shortcut + `cluster=` on `aggregate="overall"` or `aggregate="event_study"` with `cband=False` continues to work — returns the weighted-CR1 pweight sandwich bit-exact with `estimatr::iv_robust(..., weights=, clusters=, se_type="stata")`. +- **Deviation: `cluster=` + `survey=` rejected (both designs):** the `survey=` path composes Binder-TSL variance via `compute_survey_if_variance`, which would silently overwrite the cluster-robust sandwich while result metadata still reports `vcov_type="cr1"`. Rejected up-front on `cluster=` + `survey=` for BOTH designs (`continuous_*` and `mass_point`), static and event-study — route clustering through `survey_design=SurveyDesign(psu=)` instead. All other `cluster=` compositions now WORK end-to-end (pointwise CIs AND the simultaneous band): unweighted `cluster=`, and `weights=` + `cluster=` with `cband` either `False` (pointwise cluster-robust only) or `True` (adds the clustered sup-t band — see "Note (HAD clustered event-study sup-t band)" below). The `weights=` shortcut + `cluster=` pweight sandwich remains bit-exact with `estimatr::iv_robust(..., weights=, clusters=, se_type="stata")`. +- **Note (HAD clustered event-study sup-t band):** when `cluster=` is set, `_fit_event_study` produces cluster-robust per-horizon pointwise CIs AND a cluster-robust simultaneous band on BOTH designs (continuous CCT and mass-point 2SLS), unweighted or weighted (`weights=`; `survey=` is rejected — see the deviation above). Pointwise: `cluster_arr` threads into `bias_corrected_local_linear` (continuous, static-path parity) or `_fit_mass_point_2sls` (mass-point CR1). Band: `_sup_t_multiplier_bootstrap` takes a dedicated **clustered branch** — it aggregates the per-unit β̂-scale influence function to cluster level (`s_c = Σ_{i∈c} ψ_i`) and draws cluster-level iid Rademacher multipliers, so the perturbation variance is the RAW cluster sandwich `Σ_c s_c²` (no stratum-centering / FPC / Bessel — distinct from the survey branch). This matches each path's analytical cluster-robust SE via a path scalar: **1.0** for continuous (the `lprobust_vce` cluster meat carries no `g/(g-1)` correction, so `Σ_c s_c² == se_rb²` exactly) and **`√(G/(G-1))`** for mass-point (the returned IF carries `√((n-1)/(n-k))` but not the CR1 `G/(G-1)` factor; `G` is the full-array cluster count, identical to the bootstrap branch's own `n_clusters`, so a wholly-zero-weight cluster contributes `s_c=0` to both the analytical Ω and the bootstrap yet is counted in `G` by both). The variance-family reconciliation is validated bootstrap-free: `sqrt(Σ_c (scale·s_c)²) == se` to `atol=1e-10` on the real IF for both paths (`TestEventStudyClusterBand::test_{continuous,masspoint}_if_reconciliation_deterministic`), plus the `H=1 → 1.96` reduction on the clustered branch (`TestSupTReducesToNormalAtH1::test_clustered_sup_t_h1_reduces_to_normal`). Single cluster (`G<2`) → NaN band + `RuntimeWarning` (CR undefined). `cband_method="cluster_multiplier_bootstrap"` on this path. No R anchor (no reference package computes an HAD clustered sup-t band); the reconciliation identity is the validation. - 2SLS (Design 1 mass-point case): standard 2SLS inference (details not elaborated in the paper). - TWFE with small `G`: HC2 standard errors with Bell-McCaffrey (2002) degrees-of-freedom correction, following Imbens and Kolesar (2016). Used in the Pierce and Schott (2016) application with `G=103`. Added library-wide to `diff_diff/linalg.py` as a new `vcov_type` dispatch (Phase 1a), exposed on `DifferenceInDifferences` and `TwoWayFixedEffects`. - Bootstrap: wild bootstrap with Mammen (1993) two-point weights is used for the Stute test (see Diagnostics below), NOT for the main WAS estimator. Reuses the existing `diff_diff.bootstrap_utils.generate_bootstrap_weights(..., weight_type="mammen")` helper. @@ -3176,7 +3177,7 @@ Shipped in `diff_diff/had_pretests.py` as `stute_joint_pretest()` (residuals-in *Notes #1-#2 lock implementation choices (paper-permitted choices the library codified); Notes #3-#4 document validation-harness work waived in this PR with documented rationale; #5 is a Library extension where the library departs from the paper's prescription toward stricter safety.* - **Note:** Equal-weighting on the continuous path. Paper does not prescribe a unit-weighting scheme on the continuous local-linear paths. Library uses per-unit equal weighting (`w_g = 1` default, matching `diff_diff/_nprobust_port.lprobust`'s default), NOT dose-cell-size weights. Practical consequence: WAS is the population-mean slope from Eq. 3 — `[E(ΔY) − lim_{d↓d̲} E(ΔY | D ≤ d)] / E(D)` (computed as `att = (mean(ΔY) − τ_bc) / mean(D)`), not a cell-size-weighted average; with cell-size weighting, units in less-densely-populated regions of the dose distribution would contribute disproportionately to the boundary slope. User-supplied `weights=` (pweight) overrides the equal-weight default and threads through as `W_combined = k((D − d̲)/h) · w_g`. Lock in `tests/test_methodology_had.py::TestHADDeviations::test_equal_weighting_is_per_row_not_per_dose_cell`. -- **Note:** Sup-t bootstrap gating. Simultaneous-band sup-t multiplier bootstrap runs only when `aggregate="event_study"` AND `(weights= or survey_design= supplied)` AND `cband=True` (default). Unweighted event-study path bit-exactly preserves pre-Phase 4.5 B numerical output (stability invariant). Setting `cband=False` on the weighted event-study path disables the bootstrap (useful for smoke-test bit-parity assertions against the unweighted path at uniform weights). See the algorithmic contract above at `_sup_t_multiplier_bootstrap`. +- **Note:** Sup-t bootstrap gating. Simultaneous-band sup-t multiplier bootstrap runs when `aggregate="event_study"` AND `cband=True` (default) AND either `(weights= or survey_design= supplied)` (weighted/survey band) OR `cluster=` (cluster-robust band — fires even on an unweighted fit, Phase 2b). The unweighted, unclustered event-study path bit-exactly preserves pre-Phase 4.5 B numerical output (stability invariant). Setting `cband=False` disables the bootstrap on any path. See the algorithmic contract above at `_sup_t_multiplier_bootstrap`. - **Note:** Pierce-Schott (2016) Figure 2 replication harness deferred. The paper's empirical application self-acknowledges (Section 5.2; mirrored in `dechaisemartin-2026-review.md:321`) that "NP estimators are too noisy to be informative" on the LBD-restricted PNTR panel. R parity at `atol=1e-8` on 3 DGPs × 5 method combos via `tests/test_did_had_parity.py` (bit-exact, `rtol=0`) is a stronger correctness anchor than reproducing pointwise CIs on LBD-restricted data. **Scope caveat:** R parity locks point estimate, SE, and CI bounds bit-exactly to R's bounds — it does NOT independently verify the asymptotic-coverage properties of the bias-corrected CI in small samples. Paper Table 1 documents under-coverage at small G (89% at G=100 on DGP 1, 93% at G=500, 95% at G=2500); this is inherited from the CCF asymptotic theory itself, and Python is exact-parity with R at the limit-law machinery. - **Note:** Table 1 coverage-rate reproduction deferred. Paper Section 3.1.5 reports 2,000-iter Monte Carlo coverage rates at `G ∈ {100, 500, 2500}` on DGPs 1/2/3. The existing `tests/test_did_had_parity.py` R parity at `atol=1e-8` on the same 3 DGPs reproduces the exact point estimate and SE algorithm to bit-exact tolerance; coverage-rate MC would re-verify the CCF asymptotic coverage already pinned by R parity (Python ≡ R ≡ paper) at the sample-mean level. **Scope caveat (mirrors above):** R parity does NOT re-prove asymptotic-coverage at small G; paper Table 1's 89% / 93% / 95% under-coverage band is valid for both R and Python. - **Library extension:** Staggered-timing fail-closed. Paper Appendix B.2 prescribes "Warn" when staggered treatment timing is detected; library raises `ValueError` at `diff_diff/had.py:1511` when multiple first-treat cohorts are detected without `first_treat_col`. Library extension toward stricter safety: `UserWarning` would let the silent-misuse bug class through (HAD's Appendix B.2 only identifies the LAST cohort under staggered timing); fail-closed forces the user to either supply `first_treat_col` (which activates auto-filter to last-cohort + never-treated per Appendix B.2) or redirect to `ChaisemartinDHaultfoeuille` (`did_multiplegt_dyn`). Lock in `tests/test_methodology_had.py::TestHADDeviations`. @@ -3203,12 +3204,13 @@ Shipped in `diff_diff/had_pretests.py` as `stute_joint_pretest()` (residuals-in - [x] Phase 2a: Panel validator (`diff_diff.had._validate_had_panel`) verifies `D_{g,1} = 0` for all units, rejects negative post-period doses (`D_{g,2} < 0`) front-door on the original (unshifted) scale, rejects `>2` time periods on the `aggregate="overall"` path (multi-period panels must use `aggregate="event_study"`, Phase 2b), and rejects unbalanced panels and NaN in outcome/dose/unit columns. Both Design 1 paths (`continuous_near_d_lower` and `mass_point`) additionally require `d_lower == float(d.min())` within float tolerance; mismatched overrides raise with a pointer to the unsupported (LATE-like / off-support) estimand. - [x] Phase 2a: NaN-propagation tests covering constant-y, degenerate-mass-point, and single-cluster-CR1 inputs. The guaranteed NaN coupling is on the DOWNSTREAM triple (`t_stat`, `p_value`, `conf_int`) via the `safe_inference()` gate, which returns NaN on all three whenever `se` is non-finite, zero, or negative. `att` and `se` themselves are raw estimator outputs: on constant-y / no-dose-variation / divide-by-zero the fit paths return `(att=nan, se=nan)` so all five fields move to NaN together; on the degenerate single-cluster CR1 configuration on the mass-point path, `_fit_mass_point_2sls` returns `(att=beta_hat, se=nan)` - `att` is finite (Wald-IV is well defined) while `se` is NaN, so the downstream triple is NaN while `att` remains the raw 2SLS coefficient. The `assert_nan_inference` fixture in `tests/conftest.py` checks the downstream triple against this contract without requiring `att` to be NaN. - **Note (mass-point SE):** Standard errors on the mass-point path use the structural-residual 2SLS sandwich `[Z'X]^{-1} · Ω · [Z'X]^{-T}` with `Ω` built from the structural residuals `u = ΔY - α̂ - β̂·D` (not the reduced-form residuals from an OLS-on-indicator shortcut). Supported: `classical`, `hc1`, and CR1 (cluster-robust) when `cluster=` is supplied. `hc2` and `hc2_bm` raise `NotImplementedError` pending a 2SLS-specific leverage derivation (the OLS leverage `x_i' (X'X)^{-1} x_i` is wrong for 2SLS; the correct finite-sample correction depends on `(Z'X)^{-1}` rather than `(X'X)^{-1}`) plus a dedicated R parity anchor. Queued for the follow-up PR. - - **Note (continuous cluster-robust SE, Phase 2a):** On the continuous designs (`continuous_at_zero` / `continuous_near_d_lower`), `cluster=` threads the per-unit cluster IDs into `bias_corrected_local_linear`, whose CCT-2014 robust variance then becomes the cluster-robust nonparametric SE; the β̂-scale SE is `se_robust / |den|`. Because the β-scale rescale is a deterministic linear transform, the estimator-level clustered SE equals the direct `bias_corrected_local_linear(cluster=...).se_robust / |den|` to machine precision — this is the in-library validation anchor, and the clustered CCT SE is itself golden-tested against `nprobust` on DGP 4 (see the Phase 1c note above: `atol=1e-14` in first-appearance cluster order). Cluster IDs must be unit-constant (validated up front; a nonexistent column, NaN, or within-unit-varying cluster now raises rather than being silently ignored, mirroring the mass-point path). Cluster-robust inference is **unidentified with fewer than two clusters in the active kernel window** (`eC = cluster[ind]`, the in-bandwidth subset the CCT variance is actually computed on — a stricter condition than the global cluster count, since clusters can be separated from the boundary by the bandwidth): the guard lives in `_nprobust_port.lprobust` and NaNs `se_rb`/`se_cl`, so `bias_corrected_local_linear.se_robust` is NaN and HAD's `safe_inference` NaNs the t-stat / p-value / CI while the point estimate stays finite — the same single-cluster NaN contract as the mass-point CR1 path (`_fit_mass_point_2sls`), applied at the variance-computation site so it also covers the direct `bias_corrected_local_linear` API. A window with ≥2 distinct clusters is bit-identical (DGP-4 golden parity preserved). `cluster=` composes with the `weights=` shortcut (weighted cluster-robust; validated against the direct weighted+clustered local-linear call). The `cluster=` + `survey_design=` composition raises `NotImplementedError`: the Binder (1983) TSL survey variance is composed from the per-unit influence function via `compute_survey_if_variance` and would silently override the cluster-robust SE — route clustering through `survey_design=SurveyDesign(psu=)` instead. Result metadata reports `vcov_type="cr1"` + `cluster_name=` with `inference_method="analytical_nonparametric"` (distinguishing the clustered CCT variance from the mass-point 2SLS CR1 sandwich). Estimator-level cluster threading on the **Phase 2b event-study** path remains a documented follow-up (the per-horizon `cband` sup-t bootstrap normalizes HC-scale perturbations and would mix variance families under clustering); that path still emits the "cluster ignored" `UserWarning`. + - **Note (continuous cluster-robust SE, Phase 2a):** On the continuous designs (`continuous_at_zero` / `continuous_near_d_lower`), `cluster=` threads the per-unit cluster IDs into `bias_corrected_local_linear`, whose CCT-2014 robust variance then becomes the cluster-robust nonparametric SE; the β̂-scale SE is `se_robust / |den|`. Because the β-scale rescale is a deterministic linear transform, the estimator-level clustered SE equals the direct `bias_corrected_local_linear(cluster=...).se_robust / |den|` to machine precision — this is the in-library validation anchor, and the clustered CCT SE is itself golden-tested against `nprobust` on DGP 4 (see the Phase 1c note above: `atol=1e-14` in first-appearance cluster order). Cluster IDs must be unit-constant (validated up front; a nonexistent column, NaN, or within-unit-varying cluster now raises rather than being silently ignored, mirroring the mass-point path). Cluster-robust inference is **unidentified with fewer than two clusters in the active kernel window** (`eC = cluster[ind]`, the in-bandwidth subset the CCT variance is actually computed on — a stricter condition than the global cluster count, since clusters can be separated from the boundary by the bandwidth): the guard lives in `_nprobust_port.lprobust` and NaNs `se_rb`/`se_cl`, so `bias_corrected_local_linear.se_robust` is NaN and HAD's `safe_inference` NaNs the t-stat / p-value / CI while the point estimate stays finite — the same single-cluster NaN contract as the mass-point CR1 path (`_fit_mass_point_2sls`), applied at the variance-computation site so it also covers the direct `bias_corrected_local_linear` API. A window with ≥2 distinct clusters is bit-identical (DGP-4 golden parity preserved). `cluster=` composes with the `weights=` shortcut (weighted cluster-robust; validated against the direct weighted+clustered local-linear call). The `cluster=` + `survey_design=` composition raises `NotImplementedError`: the Binder (1983) TSL survey variance is composed from the per-unit influence function via `compute_survey_if_variance` and would silently override the cluster-robust SE — route clustering through `survey_design=SurveyDesign(psu=)` instead. Result metadata reports `vcov_type="cr1"` + `cluster_name=` with `inference_method="analytical_nonparametric"` (distinguishing the clustered CCT variance from the mass-point 2SLS CR1 sandwich). Estimator-level cluster threading also extends to the **Phase 2b event-study** path (cluster-robust per-horizon pointwise CIs on both designs plus a cluster-robust simultaneous sup-t band; the per-horizon variance-family reconciliation is documented in "Note (HAD clustered event-study sup-t band)" above) — the former "cluster ignored" `UserWarning` is removed. - **Note (Design 1 identification):** `continuous_near_d_lower` and `mass_point` fits emit a `UserWarning` surfacing that `WAS_{d̲}` identification requires Assumption 6 (or Assumption 5 for sign identification only) beyond parallel trends, and that neither is testable via pre-trends. `continuous_at_zero` (Design 1', Assumption 3 only) does not emit this warning. - **Note (CI endpoints):** Because the continuous-path `att` is `(mean(ΔY) - τ̂_bc) / den`, the beta-scale CI endpoints reverse relative to the Phase 1c boundary-limit CI: `CI_lower(β̂) = (mean(ΔY) - CI_upper(τ̂_bc)) / den` and `CI_upper(β̂) = (mean(ΔY) - CI_lower(τ̂_bc)) / den`. The `HeterogeneousAdoptionDiD.fit()` implementation computes `att ± z · se` directly via `safe_inference`, which handles the reversal naturally from the transformed point estimate. - **Note (Phase 2a/2b scope, superseded by Phase 4.5):** Phase 2a ships the single-period `aggregate="overall"` path; Phase 2b lifts `aggregate="event_study"` (Appendix B.2 multi-period extension) which returns a `HeterogeneousAdoptionDiDEventStudyResults` with per-event-time WAS estimates and pointwise CIs. The original Phase 2a/2b release raised `NotImplementedError` on `survey=` and `weights=`, but Phase 4.5 (A/B/C0) lifted both gates with the per-design vcov contract documented above (see L2340-L2379, including the mass-point `vcov_type="classical"` deviation and `cband=True` sup-t restriction); the `survey=` path composes Binder (1983) TSL via `compute_survey_if_variance` on both continuous and mass-point IFs. - **Note (panel-only):** The paper (Section 2) defines HAD on *panel or repeated cross-section* data, but both the overall and event-study paths ship a panel-only implementation: `HeterogeneousAdoptionDiD.fit()` requires a balanced panel with a unit identifier so that unit-level first differences `ΔY_{g,t} = Y_{g,t} - Y_{g,t_anchor}` can be formed. Repeated-cross-section inputs (disjoint unit IDs between periods) are rejected by the balanced-panel validator. RCS support is queued for a follow-up PR (tracked in `TODO.md`); it will need a separate identification path based on pre/post cell means rather than unit-level differences. - [x] Phase 2b: Multi-period event-study extension (Appendix B.2). `aggregate="event_study"` produces per-event-time WAS estimates using a uniform `F-1` baseline (`ΔY_{g,t} = Y_{g,t} - Y_{g,F-1}` for every horizon), reusing the three Phase 2a design paths on per-horizon first differences. Pre-period placebos included for `e <= -2` (the anchor `e = -1` is skipped since `ΔY = 0` trivially). Post-period estimates for `e >= 0`. The joint Stute test (Equation 18) across pre-periods is a SEPARATE diagnostic deferred to a **Phase 3 follow-up patch** (Phase 3 ships the single-horizon Stute test; the joint stacked-residual variant is tracked in `TODO.md`). +- [x] Phase 2b: event-study `cluster=` threading — cluster-robust per-horizon pointwise CIs (both designs) AND a cluster-robust simultaneous band via the clustered branch of `_sup_t_multiplier_bootstrap` (continuous scale 1.0; mass-point `√(G/(G-1))`). Closes the former "cluster ignored on the nonparametric path" deferral (`TODO.md`). See "Note (HAD clustered event-study sup-t band)". `cluster=` + `survey=` remains rejected (route through `survey_design=SurveyDesign(psu=...)`). - **Note (Phase 2b last-cohort filter):** When `first_treat_col` indicates more than one nonzero cohort, the panel is auto-filtered to the last-treatment cohort (`F_last = max(cohorts)`) **plus never-treated units** (`first_treat = 0`), with a `UserWarning` naming kept/dropped unit counts and dropped cohort labels. Paper Appendix B.2 is explicit that HAD "may be used only for the LAST treatment cohort in a staggered design"; the auto-filter implements this prescription, retaining never-treated units per the paper's "there must be an untreated group, at least till the period where the last cohort gets treated" requirement. Only earlier-cohort units (with `first_treat > 0` and `< F_last`) are dropped — never-treated units satisfy the dose invariant at every period (`D = 0` throughout) and preserve Design 1' identifiability (boundary at `0`) when last-cohort doses are uniformly positive. When `first_treat_col` is omitted on a >2-period panel, the validator infers each unit's first-positive-dose period from the dose path; if multiple distinct first-positive-dose cohorts are detected, the estimator raises a front-door `ValueError` directing users to pass `first_treat_col` (which activates the auto-filter) or use `ChaisemartinDHaultfoeuille` for full staggered support — there is no silent acceptance of staggered panels without cohort metadata. Common-adoption panels (single first-positive-dose cohort, or only never-treated + one cohort) pass through unchanged with `F` inferred from the dose invariant, and require dose contiguity (pre-periods < post-periods in natural ordering). Non-contiguous dose sequences (e.g., reverse treatment) raise with a pointer to `ChaisemartinDHaultfoeuille`. - **Note (Phase 2b constant-dose requirement):** The event-study aggregation uses `D_{g, F}` (first-treatment-period dose) as the single regressor for every event-time horizon, per paper Appendix B.2's "once treated, stay treated with the same dose" convention. The validator REJECTS panels where a unit has time-varying dose across post-treatment periods (`D_{g, t} != D_{g, F}` for any `t >= F` within-unit, beyond float tolerance) with a front-door `ValueError`, directing users with genuinely time-varying post-treatment doses to `ChaisemartinDHaultfoeuille` (`did_multiplegt_dyn`). Silent acceptance would misattribute later-horizon treatment-effect heterogeneity to the period-F dose. A follow-up PR could implement a time-varying-dose estimator; tracked in `TODO.md`. - **Note (Phase 2b per-horizon SE):** Each event-time horizon uses an INDEPENDENT sandwich computed on that horizon's first differences: continuous paths use the CCT-2014 robust SE from Phase 1c divided by `|den|`; mass-point path uses the structural-residual 2SLS sandwich from Phase 2a. This produces pointwise CIs per horizon, matching the paper's Pierce-Schott application (Section 5.2, Figure 2: "nonparametric pointwise CIs"). Joint cross-horizon covariance (IF-based stacking or block bootstrap) is NOT computed — the paper does not derive it and all reported CIs are pointwise. Follow-up PRs may add joint covariance for cross-horizon hypothesis tests; current tracking in `TODO.md`. diff --git a/tests/test_had.py b/tests/test_had.py index bfbd39db4..d81ab0c7b 100644 --- a/tests/test_had.py +++ b/tests/test_had.py @@ -3277,7 +3277,7 @@ def test_sklearn_clone_round_trip(self): class TestEventStudyWarnings: - """Continuous-path warnings on event-study mode (vcov/robust/cluster ignored).""" + """Continuous-path warnings on event-study mode (vcov/robust ignored; cluster= is now threaded).""" def _panel(self): rng = np.random.default_rng(0) @@ -3297,17 +3297,30 @@ def test_vcov_type_ignored_on_continuous(self): ] assert len(vcov_warnings) == 1 # ONE per fit, not per horizon - def test_cluster_ignored_on_continuous(self): + def test_cluster_threaded_on_continuous_event_study(self): + # Phase 2b: cluster= is now threaded into the per-horizon CCT SE on + # the continuous event-study path (no longer ignored). No "ignored" + # warning; SE becomes cluster-robust; result surface labels the + # cluster-robust variance. panel = self._panel() panel["state"] = panel["unit"] % 20 - est = HeterogeneousAdoptionDiD(design="auto", cluster="state") with warnings.catch_warnings(record=True) as w: warnings.simplefilter("always") - est.fit(panel, "outcome", "dose", "period", "unit", aggregate="event_study") + r_cl = HeterogeneousAdoptionDiD(design="auto", cluster="state").fit( + panel, "outcome", "dose", "period", "unit", aggregate="event_study" + ) cluster_warnings = [ - msg for msg in w if "cluster=" in str(msg.message) and "ignored" in str(msg.message) + msg + for msg in w + if "cluster=" in str(msg.message) and "ignored" in str(msg.message).lower() ] - assert len(cluster_warnings) == 1 + assert len(cluster_warnings) == 0 + r_un = HeterogeneousAdoptionDiD(design="auto").fit( + panel, "outcome", "dose", "period", "unit", aggregate="event_study" + ) + assert not np.allclose(r_cl.se, r_un.se) + assert r_cl.vcov_type == "cr1" + assert r_cl.cluster_name == "state" class TestEventStudyValidator: @@ -4756,6 +4769,275 @@ def test_sup_t_seed_reproducibility(self): ) assert q1 == q2 + def test_clustered_sup_t_h1_reduces_to_normal(self): + """Clustered branch: at H=1 the sup-t crit reduces to the Normal + quantile for BOTH the continuous scalar (1.0) and the mass-point + CR1 scalar sqrt(G/(G-1)) — the decisive variance-family check.""" + import scipy.stats + + from diff_diff.had import _sup_t_multiplier_bootstrap + + rng = np.random.default_rng(3) + n_units, n_clusters = 300, 30 + cluster_ids = rng.integers(0, n_clusters, size=n_units) + expected = float(scipy.stats.norm.ppf(0.975)) + + def _cluster_se(psi, scale): + labels = pd.unique(cluster_ids) + lab = {v: r for r, v in enumerate(labels)} + pcl = np.zeros((len(labels), psi.shape[1])) + for i in range(n_units): + pcl[lab[cluster_ids[i]]] += psi[i] + return np.sqrt(((scale * pcl) ** 2).sum(axis=0)) + + for scale in (1.0, float(np.sqrt(n_clusters / (n_clusters - 1)))): + psi = rng.standard_normal((n_units, 1)) + se = _cluster_se(psi, scale) + q, _lo, _hi, n_valid = _sup_t_multiplier_bootstrap( + psi, + np.zeros(1), + se, + None, + n_bootstrap=20000, + alpha=0.05, + seed=7, + cluster_ids=cluster_ids, + cluster_if_scale=scale, + ) + assert n_valid >= 18000 + assert abs(q - expected) < 0.10, ( + f"clustered H=1 sup-t (scale={scale:.4f}) should reduce to " + f"Phi^-1(0.975)={expected:.4f}; got {q:.4f} — variance-family " + f"drift in the clustered bootstrap branch." + ) + + def test_clustered_sup_t_single_cluster_nan(self): + """One cluster → NaN crit / None band (CR undefined).""" + from diff_diff.had import _sup_t_multiplier_bootstrap + + rng = np.random.default_rng(0) + psi = rng.standard_normal((100, 1)) + q, lo, hi, n_valid = _sup_t_multiplier_bootstrap( + psi, + np.zeros(1), + np.array([1.0]), + None, + n_bootstrap=500, + alpha=0.05, + seed=1, + cluster_ids=np.zeros(100, dtype=int), + ) + assert np.isnan(q) and lo is None and hi is None and n_valid == 0 + + +class TestEventStudyClusterBand: + """Phase 2b: cluster-robust event-study pointwise CIs + clustered sup-t + simultaneous band (continuous + mass-point). The core reconciliation is + that the cluster-aggregated influence function reproduces the analytical + cluster-robust SE — exactly (continuous, scale 1.0) or after the CR1 + sqrt(G/(G-1)) scalar (mass-point).""" + + @staticmethod + def _clustered_panel(G=240, n_clusters=24, seed=0): + rng = np.random.default_rng(seed) + d = np.where(rng.random(G) < 0.15, 0.0, rng.uniform(0.2, 1.2, size=G)) + d[0] = 0.0 + state = np.repeat(np.arange(n_clusters), G // n_clusters) + panel = _make_multi_period_panel( + d, n_periods=5, F=3, seed=seed, extra_cols={"state": state} + ) + return panel + + def test_continuous_if_reconciliation_deterministic(self): + """sqrt(sum_c (sum_{i in c} IF_i)^2) == se_robust for the cluster- + robust CCT fit (scale 1.0) — bootstrap-free proof of the continuous + reconciliation on the REAL influence function.""" + d, dy, cl = self._make_cluster_dgp(seed=1) + bc = bias_corrected_local_linear(d=d, y=dy, boundary=0.0, cluster=cl, return_influence=True) + recon = self._cluster_agg_norm(bc.influence_function, cl, scale=1.0) + np.testing.assert_allclose(recon, bc.se_robust, rtol=0, atol=1e-10) + + def test_masspoint_if_reconciliation_deterministic(self): + """sqrt(sum_c (sqrt(G/(G-1)) * sum_{i in c} psi_i)^2) == se for the + mass-point CR1 fit — bootstrap-free proof the sqrt(G/(G-1)) scalar + exactly restores the CR1 finite-sample factor on the REAL IF.""" + from diff_diff.had import _fit_mass_point_2sls + + rng = np.random.default_rng(2) + G, d_lower = 240, 0.5 + mass_n = G // 3 + d = np.concatenate([np.full(mass_n, d_lower), rng.uniform(d_lower, 1.0, G - mass_n)]) + rng.shuffle(d) + cl = np.arange(G) % 20 + shock = rng.normal(scale=0.5, size=20)[cl] + dy = 0.3 * d + shock + 0.1 * rng.standard_normal(G) + _beta, se, psi = _fit_mass_point_2sls(d, dy, d_lower, cl, "hc1", return_influence=True) + n_cl = len(pd.unique(cl)) + recon = self._cluster_agg_norm(psi, cl, scale=float(np.sqrt(n_cl / (n_cl - 1)))) + np.testing.assert_allclose(recon, se, rtol=0, atol=1e-10) + + def test_continuous_clustered_band_end_to_end(self): + import scipy.stats + + panel = self._clustered_panel(seed=2) + r = HeterogeneousAdoptionDiD( + design="continuous_at_zero", cluster="state", n_bootstrap=1500, seed=11 + ).fit(panel, "outcome", "dose", "period", "unit", aggregate="event_study", cband=True) + assert r.vcov_type == "cr1" and r.cluster_name == "state" + assert r.cband_low is not None and np.all(np.isfinite(r.cband_low)) + assert r.cband_high is not None and np.all(np.isfinite(r.cband_high)) + assert r.cband_method == "cluster_multiplier_bootstrap" + assert r.cband_crit_value >= float(scipy.stats.norm.ppf(0.975)) - 0.05 + # Simultaneous band is at least as wide as the pointwise CI. + assert np.all(r.cband_low <= r.conf_int_low + 1e-9) + assert np.all(r.cband_high >= r.conf_int_high - 1e-9) + + def test_masspoint_weighted_cluster_cband_no_longer_raises(self): + """Symmetry: weighted mass-point + cluster= + cband=True used to raise + NotImplementedError; the clustered band now reconciles it.""" + rng = np.random.default_rng(4) + G, d_lower = 240, 0.5 + mass_n = G // 3 + d = np.concatenate([np.full(mass_n, d_lower), rng.uniform(d_lower, 1.0, G - mass_n)]) + rng.shuffle(d) + state = np.arange(G) % 24 + panel = _make_multi_period_panel(d, n_periods=5, F=3, seed=5, extra_cols={"state": state}) + w_unit = rng.uniform(0.5, 2.0, G) + w_row = w_unit[panel["unit"].to_numpy()] + r = HeterogeneousAdoptionDiD( + design="mass_point", cluster="state", d_lower=d_lower, n_bootstrap=1500, seed=9 + ).fit( + panel, + "outcome", + "dose", + "period", + "unit", + aggregate="event_study", + weights=w_row, + cband=True, + ) + assert r.vcov_type == "cr1" and r.cluster_name == "state" + assert r.cband_low is not None and np.all(np.isfinite(r.cband_low)) + assert r.cband_method == "cluster_multiplier_bootstrap" + + def test_weighted_continuous_clustered_band_end_to_end(self): + """Weighted continuous event-study + cluster= + cband: finite + cluster-robust band (the weighted arm of the continuous path).""" + rng = np.random.default_rng(6) + G = 240 + d = np.where(rng.random(G) < 0.15, 0.0, rng.uniform(0.2, 1.2, size=G)) + d[0] = 0.0 + state = np.repeat(np.arange(24), G // 24) + panel = _make_multi_period_panel(d, n_periods=5, F=3, seed=6, extra_cols={"state": state}) + w_unit = rng.uniform(0.5, 2.0, G) + w_row = w_unit[panel["unit"].to_numpy()] + r = HeterogeneousAdoptionDiD( + design="continuous_at_zero", cluster="state", n_bootstrap=1500, seed=13 + ).fit( + panel, + "outcome", + "dose", + "period", + "unit", + aggregate="event_study", + weights=w_row, + cband=True, + ) + assert r.vcov_type == "cr1" and r.cluster_name == "state" + assert r.cband_low is not None and np.all(np.isfinite(r.cband_low)) + assert r.cband_method == "cluster_multiplier_bootstrap" + + def test_unweighted_masspoint_clustered_band_end_to_end(self): + """Unweighted mass-point event-study + cluster= + cband: finite + cluster-robust band (the unweighted arm of the mass-point path).""" + rng = np.random.default_rng(7) + G, d_lower = 240, 0.5 + mass_n = G // 3 + d = np.concatenate([np.full(mass_n, d_lower), rng.uniform(d_lower, 1.0, G - mass_n)]) + rng.shuffle(d) + state = np.arange(G) % 24 + panel = _make_multi_period_panel(d, n_periods=5, F=3, seed=7, extra_cols={"state": state}) + r = HeterogeneousAdoptionDiD( + design="mass_point", cluster="state", d_lower=d_lower, n_bootstrap=1500, seed=17 + ).fit(panel, "outcome", "dose", "period", "unit", aggregate="event_study", cband=True) + assert r.vcov_type == "cr1" and r.cluster_name == "state" + assert r.cband_low is not None and np.all(np.isfinite(r.cband_low)) + assert r.cband_method == "cluster_multiplier_bootstrap" + + def test_cluster_survey_event_study_raises(self): + from diff_diff.survey import SurveyDesign + + panel = self._clustered_panel(seed=3) + panel["w"] = 1.0 + with pytest.raises(NotImplementedError, match=r"cluster.*\+ survey="): + HeterogeneousAdoptionDiD(design="continuous_at_zero", cluster="state").fit( + panel, + "outcome", + "dose", + "period", + "unit", + aggregate="event_study", + survey_design=SurveyDesign(weights="w"), + ) + + def test_single_cluster_band_nan_and_warns(self): + # Dense dose (mirrors the static single-cluster test) so the local- + # linear fit resolves gracefully to NaN CR SEs rather than a + # degenerate-bandwidth error; the band guard then fires. + rng = np.random.default_rng(0) + G = 300 + d = rng.uniform(0.0, 1.0, G) + d[0] = 0.0 + panel = _make_multi_period_panel( + d, n_periods=5, F=3, seed=0, extra_cols={"state": np.zeros(G, dtype=int)} + ) + with warnings.catch_warnings(record=True) as w: + warnings.simplefilter("always") + r = HeterogeneousAdoptionDiD( + design="continuous_at_zero", cluster="state", n_bootstrap=500, seed=1 + ).fit(panel, "outcome", "dose", "period", "unit", aggregate="event_study", cband=True) + assert r.cband_low is None and r.cband_high is None + assert any("single cluster" in str(x.message).lower() for x in w) + # "Undefined band" (crit=NaN, method/count populated), NOT "band + # skipped" (all None) — distinguishable by the caller. + assert r.cband_crit_value is not None and np.isnan(r.cband_crit_value) + assert r.cband_method == "cluster_multiplier_bootstrap" + assert r.cband_n_bootstrap == 500 + + def test_clustered_band_determinism(self): + panel = self._clustered_panel(seed=2) + kw = dict(design="continuous_at_zero", cluster="state", n_bootstrap=800, seed=21) + fit_kw = dict(aggregate="event_study", cband=True) + r1 = HeterogeneousAdoptionDiD(**kw).fit( + panel, "outcome", "dose", "period", "unit", **fit_kw + ) + r2 = HeterogeneousAdoptionDiD(**kw).fit( + panel, "outcome", "dose", "period", "unit", **fit_kw + ) + np.testing.assert_array_equal(r1.cband_low, r2.cband_low) + np.testing.assert_array_equal(r1.cband_high, r2.cband_high) + assert r1.cband_crit_value == r2.cband_crit_value + + # ---- helpers ---- + @staticmethod + def _make_cluster_dgp(G=200, n_clusters=20, seed=0): + rng = np.random.default_rng(seed) + d = np.where(rng.random(G) < 0.15, 0.0, rng.uniform(0.2, 1.5, size=G)) + d[0] = 0.0 + cl = np.repeat(np.arange(n_clusters), G // n_clusters) + shock = rng.normal(scale=1.0, size=n_clusters)[cl] + dy = 1.5 * d + shock + rng.normal(scale=0.3, size=G) + return d, dy, cl + + @staticmethod + def _cluster_agg_norm(if_vec, cl, scale): + labels = pd.unique(cl) + lab = {v: r for r, v in enumerate(labels)} + s = np.zeros(len(labels)) + for i in range(len(cl)): + s[lab[cl[i]]] += if_vec[i] + return float(np.sqrt(np.sum((scale * s) ** 2))) + class TestEventStudySurveyCband: """Event-study + weights / survey + sup-t cband scope (Phase 4.5 B).""" @@ -5046,9 +5328,10 @@ def test_mass_point_survey_plus_cluster_rejected_static(self): with pytest.raises(NotImplementedError, match="cluster"): est.fit(panel, "outcome", "dose", "period", "unit", survey=sd) - def test_mass_point_survey_plus_cluster_rejected_event_study(self): - """Review R2 P1 (event-study arm): same rejection must fire on - the multi-period dispatch.""" + def test_mass_point_weights_plus_cluster_event_study_supported(self): + """Phase 2b: mass-point + weights= + cluster= + cband (default True) + is now SUPPORTED (was rejected) — the clustered sup-t band reconciles + the CR1 variance family via the sqrt(G/(G-1)) scalar.""" rng = np.random.default_rng(1) G, T = 150, 4 d_mp = np.concatenate([np.full(30, 0.3), rng.uniform(0.3, 1.0, G - 30)]) @@ -5062,17 +5345,21 @@ def test_mass_point_survey_plus_cluster_rejected_event_study(self): panel = pd.DataFrame(rows, columns=["unit", "period", "dose", "outcome", "state"]) with warnings.catch_warnings(): warnings.simplefilter("ignore", UserWarning) - est = HeterogeneousAdoptionDiD(design="mass_point", vcov_type="hc1", cluster="state") - with pytest.raises(NotImplementedError, match="cluster"): - est.fit( - panel, - "outcome", - "dose", - "period", - "unit", - aggregate="event_study", - weights=np.ones(panel.shape[0]), - ) + est = HeterogeneousAdoptionDiD( + design="mass_point", vcov_type="hc1", cluster="state", seed=0, n_bootstrap=500 + ) + r = est.fit( + panel, + "outcome", + "dose", + "period", + "unit", + aggregate="event_study", + weights=np.ones(panel.shape[0]), + ) + assert r.vcov_type == "cr1" and r.cluster_name == "state" + assert r.cband_low is not None and np.all(np.isfinite(r.cband_low)) + assert r.cband_method == "cluster_multiplier_bootstrap" def test_lonely_psu_adjust_with_singletons_rejected_on_cband(self): """Review R2 P1: sup-t bootstrap rejects lonely_psu='adjust' @@ -5376,11 +5663,11 @@ def test_mass_point_weights_plus_cluster_event_study_cband_false_allowed(self): assert r.variance_formula == "pweight_2sls" assert r.cband_crit_value is None - def test_mass_point_weights_plus_cluster_event_study_cband_true_rejected(self): - """Review R4 P1: event-study + weights= + cluster= + cband=True - IS rejected (HC1-scale bootstrap perturbations normalized by - CR1 analytical SE would mix variance families in the bootstrap - t-distribution).""" + def test_mass_point_weights_plus_cluster_event_study_cband_true_supported(self): + """Phase 2b: event-study + weights= + cluster= + cband=True is now + SUPPORTED (previously rejected). The clustered bootstrap draws + cluster-level multipliers on the per-unit IF and normalizes by the + CR1 analytical SE (variance families reconciled via sqrt(G/(G-1))).""" rng = np.random.default_rng(42) G, T = 180, 4 d_mp = np.concatenate([np.full(36, 0.3), rng.uniform(0.3, 1.0, G - 36)]) @@ -5404,19 +5691,21 @@ def test_mass_point_weights_plus_cluster_event_study_cband_true_rejected(self): vcov_type="hc1", cluster="state", seed=0, - n_bootstrap=100, + n_bootstrap=500, ) - with pytest.raises(NotImplementedError, match="cband=True"): - est.fit( - panel, - "outcome", - "dose", - "period", - "unit", - aggregate="event_study", - weights=w_row, - cband=True, - ) + r = est.fit( + panel, + "outcome", + "dose", + "period", + "unit", + aggregate="event_study", + weights=w_row, + cband=True, + ) + assert r.vcov_type == "cr1" and r.cluster_name == "state" + assert r.cband_low is not None and np.all(np.isfinite(r.cband_low)) + assert r.cband_method == "cluster_multiplier_bootstrap" def test_event_study_zero_weight_units_excluded_from_n_units(self): """Review R4 P2: weighted event-study reports the POSITIVE-WEIGHT