diff --git a/CHANGELOG.md b/CHANGELOG.md index c69f0fef..7f7ff8b5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,18 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] ### Added +- **`LPDiD` non-absorbing (reversible) treatment** (Dube, Girardi, Jordà & Taylor 2025, + Section 4.2). New `non_absorbing` parameter: `"first_entry"` (Eq. 12 — the effect of + entering treatment for the first time and staying treated) and `"effect_stabilization"` + (Eq. 13, with `stabilization_window=L` — units whose treatment has been stable for at least + `L` periods serve as clean controls, so estimation is feasible with few or no never-treated + units). The default `non_absorbing=None` is unchanged (absorbing path; still rejects + non-absorbing input, bit-for-bit identical results). Both modes implement the entry-effect + estimands with mode-aware clean-sample masks, a documented "untreated before the first + observed period" boundary convention, and a gap-free-panel requirement; the Appendix-C + exit-event dynamics, R-package parity (PR-C2), and survey-design support remain follow-ups. + Pure-Python validation covers the absorbing reduction, the re-entry mechanism, pre-trend + placebos, non-negative weighting, stabilized-control admission, and DGP recovery. - **`LPDiD` R-parity validation (absorbing).** `tests/test_methodology_lpdid.py` pins the estimator against the method authors' own R recipes (`danielegirardi/lpdid` event-study / reweight / premean / pooled `fixest::feols` specifications) with an `alexCardazzi/lpdid` diff --git a/README.md b/README.md index a33f8be7..d237b050 100644 --- a/README.md +++ b/README.md @@ -117,7 +117,7 @@ Full guide: `diff_diff.get_llm_guide("practitioner")`. - [TROP](https://diff-diff.readthedocs.io/en/stable/api/trop.html) - Triply Robust Panel estimator (Athey et al. 2025) with nuclear norm factor adjustment - [StaggeredTripleDifference](https://diff-diff.readthedocs.io/en/stable/api/staggered.html#staggeredtripledifference) - Ortiz-Villavicencio & Sant'Anna (2025) staggered DDD with group-time ATT - [WooldridgeDiD](https://diff-diff.readthedocs.io/en/stable/api/wooldridge_etwfe.html) - Wooldridge (2023, 2025) ETWFE: saturated OLS, logit/Poisson QMLE (ASF-based ATT). Alias `ETWFE`. -- [LPDiD](https://diff-diff.readthedocs.io/en/stable/api/lpdid.html) - Dube, Girardi, Jorda & Taylor (2025) Local Projections DiD: per-horizon long-difference event study on clean controls (no negative weighting), variance- or equally-weighted ATT, for absorbing treatment +- [LPDiD](https://diff-diff.readthedocs.io/en/stable/api/lpdid.html) - Dube, Girardi, Jorda & Taylor (2025) Local Projections DiD: per-horizon long-difference event study on clean controls (no negative weighting), variance- or equally-weighted ATT, for absorbing or non-absorbing (reversible) treatment - [BaconDecomposition](https://diff-diff.readthedocs.io/en/stable/api/bacon.html) - Goodman-Bacon (2021) decomposition for diagnosing TWFE bias in staggered settings ## Diagnostics & Sensitivity diff --git a/TODO.md b/TODO.md index 938270f9..9be6b18c 100644 --- a/TODO.md +++ b/TODO.md @@ -94,6 +94,8 @@ Not currently actionable. Retained for provenance + AI-review deviation-document | `SpilloverDiD(survey_design=...)` replicate-weight variance (BRR/Fay/JK1/JKn/SDR): Wave E.1 ships Taylor-linearization only. Per Gerber (2026) Appendix A the IF-reweighting shortcut does NOT apply to TwoStageDiD-class estimators (`gamma_hat` is weight-sensitive); correct support needs per-replicate full re-fit of both stages. | `spillover.py`, `survey.py::compute_replicate_refit_variance` | follow-up | Low | | `SpilloverDiD(vcov_type="conley", conley_lag_cutoff>0, survey_design=...)` no-effective-PSU serial Bartlett HAC: weights-only / strata-only designs without a cluster fallback raise `NotImplementedError` (each pseudo-PSU appears in one period, so the serial cross-period loop contributes zero). Needs a unit-level serial fallback derivation or routing through `conley_unit` with documented IF-allocator asymmetry. | `spillover.py`, `two_stage.py::_compute_stratified_serial_bartlett_meat` | Wave E.2 tail | Low | | `SpilloverDiD` data-driven `d_bar` selection (Butts 2021b / 2023 JUE Insight cross-validation). | `spillover.py` | follow-up | Low | +| **`LPDiD` non-absorbing exit-event dynamics** (Dube et al. 2025 online Appendix C `eta_h^{g,n}`): the shipped `non_absorbing` modes estimate the **entry-effect** estimands (Eq. 12/13) only; separate dynamic event-studies for treatment switch-*offs* are not implemented. Needs the exit-event clean-sample derivation + estimand contract. | `lpdid.py`, REGISTRY | PR-C follow-up | Low | +| **`LPDiD` non-absorbing interior-gap support**: non-absorbing modes require a gap-free panel within each unit's observed span and raise on interior time gaps (the `[t-L, t+h]` window conditions can't be verified across a gap). The absorbing path already reindexes interior gaps to the calendar grid; extending that fail-closed handling (per-window gap masking) to non-absorbing is deferred. | `lpdid.py::_prepare_panel` | PR-C follow-up | Low | ### Needs external reference (R / Stata / Julia) @@ -116,6 +118,7 @@ exists but parity can't be verified without a local toolchain. | **`bias_corrected_local_linear` (lprobust) Phase-1c follow-ups:** extend golden parity to `kernel ∈ {triangular, uniform}` (epa-only today); expose `vce ∈ {hc0,hc1,hc2,hc3}` on the public wrapper once R goldens exist (port supports all four; needs a per-mode generator + a hc2/hc3 q-fit-leverage decision); clustered-DGP auto-bandwidth parity is **blocked upstream** on an nprobust singleton-cluster bug in `lpbwselect.mse.dpi` (Phase-1c DGP 4 uses manual `h=b=0.3`). | `_nprobust_port.py`, `local_linear.py`, `generate_nprobust_lprobust_golden.R` | Phase 1c | Low-Med | | `HeterogeneousAdoptionDiD` Stute-family Stata-bridge parity: no public R `Stutetest` package exists; would add `benchmarks/stata/generate_stute_golden.do` + a Stata dependency. | `benchmarks/stata/`, `tests/test_stute_test_parity.py` | follow-up | Low | | **`LPDiD` regression-adjustment SE — no runnable R reference.** The RA influence-function cluster SE is canonically Stata `teffects ra ... atet vce(cluster)` only; no R package computes it (`alexCardazzi/lpdid` does direct covariate inclusion, not RA). Today the RA *point* is R-anchored (~1e-12), the SE is pinned + MC-coverage-validated (`coverage_lpdid_ra.py`). Follow-up: contribute the RA path to `alexCardazzi/lpdid` so a runnable R RA reference exists — only a *trusted* anchor once cross-checked vs Stata `teffects` (else circular). | `tests/test_methodology_lpdid.py`, `benchmarks/python/coverage_lpdid_ra.py` | #B2 follow-up | Low | +| **`LPDiD` non-absorbing R-parity (PR-C2).** The shipped non-absorbing modes (`first_entry` Eq. 12 / `effect_stabilization` Eq. 13) are validated by pure-Python tests (absorbing reduction, re-entry mechanism, placebo, non-negative weighting, DGP recovery) but not yet pinned to `alexCardazzi/lpdid`'s `nonabsorbing` / `nonabsorbing_lag` branches. PR-C2: smoke-gate the option mapping (incl. the equal-weight reweight and the pre-observation boundary convention), then extend `generate_lpdid_golden.R` + `test_methodology_lpdid.py` with non-absorbing variants. | `benchmarks/R/generate_lpdid_golden.R`, `tests/test_methodology_lpdid.py` | PR-C2 | Medium | | `HeterogeneousAdoptionDiD` Phase-3 R-parity: ships coverage-rate validation on synthetic DGPs, not tight point parity vs `chaisemartin::stute_test` / `yatchew_test` (needs bootstrap-seed-semantics + `B` alignment across numpy/R). | `tests/test_had_pretests.py` | Phase 3 | Low | ### Parked — pending user demand / out of scope diff --git a/diff_diff/guides/llms-full.txt b/diff_diff/guides/llms-full.txt index 79349a74..9278bdbd 100644 --- a/diff_diff/guides/llms-full.txt +++ b/diff_diff/guides/llms-full.txt @@ -897,7 +897,7 @@ results.print_summary() ### LPDiD -Local Projections DiD (Dube, Girardi, Jorda & Taylor 2025). Estimates a separate OLS at each event-time horizon of a long difference (`y_{i,t+h} - y_{i,t-1}`) on the treatment-switch indicator plus calendar-time fixed effects (no unit FE), restricted to a flexible "clean control" sample of newly-treated and not-yet-treated units. Excluding already-treated units from the control group removes the negative-weighting bias of naive TWFE, so the default (variance-weighted) estimand has strictly non-negative weights. `reweight=True` yields the equally-weighted ATT (numerically equivalent to Callaway-Sant'Anna); covariates then enter via regression adjustment. Standard errors on the default/weighted path are cluster-robust at the unit level (the paper specifies no SE; matches Stata `lpdid` `vce(cluster unit)`); the regression-adjustment covariate path (`reweight=True`) instead reports an influence-function cluster variance (ImputationDiD/BJS family). Scope: binary, absorbing treatment (rejects panels where treatment turns off). +Local Projections DiD (Dube, Girardi, Jorda & Taylor 2025). Estimates a separate OLS at each event-time horizon of a long difference (`y_{i,t+h} - y_{i,t-1}`) on the treatment-switch indicator plus calendar-time fixed effects (no unit FE), restricted to a flexible "clean control" sample of newly-treated and not-yet-treated units. Excluding already-treated units from the control group removes the negative-weighting bias of naive TWFE, so the default (variance-weighted) estimand has strictly non-negative weights. `reweight=True` yields the equally-weighted ATT (numerically equivalent to Callaway-Sant'Anna); covariates then enter via regression adjustment. Standard errors on the default/weighted path are cluster-robust at the unit level (the paper specifies no SE; matches Stata `lpdid` `vce(cluster unit)`); the regression-adjustment covariate path (`reweight=True`) instead reports an influence-function cluster variance (ImputationDiD/BJS family). Scope: binary treatment; absorbing by default (rejects panels where treatment turns off), with non-absorbing (reversible) treatment available via `non_absorbing` - `"first_entry"` (Dube et al. Eq. 12, the effect of entering for the first time and staying treated) or `"effect_stabilization"` (Eq. 13, requires `stabilization_window=L`; lets units whose treatment has been stable for at least `L` periods act as clean controls, so estimation is feasible with few/no never-treated units). Non-absorbing modes require a gap-free panel within each unit's observed span. ```python LPDiD( @@ -910,6 +910,8 @@ LPDiD( alpha: float = 0.05, cluster: str | None = None, # Cluster column for cluster-robust SEs; defaults to the unit identifier rank_deficient_action: str = "warn", # "warn", "error", or "silent" + non_absorbing: str | None = None, # None=absorbing; "first_entry" (Eq. 12); "effect_stabilization" (Eq. 13) + stabilization_window: int | None = None, # The paper's L; required when non_absorbing="effect_stabilization" ) ``` @@ -921,7 +923,7 @@ lpdid.fit( outcome: str, unit: str, time: str, - treatment: str, # Binary, absorbing treatment indicator (0/1) + treatment: str, # Binary treatment indicator (0/1); absorbing unless non_absorbing is set covariates: list[str] = None, # Direct inclusion (reweight=False) or regression adjustment (reweight=True) ylags: int = 0, # Lagged-outcome controls dylags: int = 0, # Lagged first-difference controls diff --git a/diff_diff/guides/llms.txt b/diff_diff/guides/llms.txt index 92abb8f2..433006bb 100644 --- a/diff_diff/guides/llms.txt +++ b/diff_diff/guides/llms.txt @@ -69,7 +69,7 @@ Full practitioner guide: call `diff_diff.get_llm_guide("practitioner")` - [TROP](https://diff-diff.readthedocs.io/en/stable/api/trop.html): Triply Robust Panel estimator (Athey et al. 2025) with nuclear norm factor adjustment - [StaggeredTripleDifference](https://diff-diff.readthedocs.io/en/stable/api/staggered.html#staggeredtripledifference): Ortiz-Villavicencio & Sant'Anna (2025) staggered DDD with group-time ATT - [WooldridgeDiD](https://diff-diff.readthedocs.io/en/stable/api/wooldridge_etwfe.html): Wooldridge (2023, 2025) ETWFE — saturated OLS, logit/Poisson QMLE (ASF-based ATT). Alias: ETWFE -- [LPDiD](https://diff-diff.readthedocs.io/en/stable/api/lpdid.html): Dube, Girardi, Jorda & Taylor (2025) Local Projections DiD: per-horizon long-difference event study on clean controls (no negative weighting); variance- or equally-weighted ATT, premean differencing, pooled pre/post, fast. Absorbing treatment. +- [LPDiD](https://diff-diff.readthedocs.io/en/stable/api/lpdid.html): Dube, Girardi, Jorda & Taylor (2025) Local Projections DiD: per-horizon long-difference event study on clean controls (no negative weighting); variance- or equally-weighted ATT, premean differencing, pooled pre/post, fast. Absorbing by default; non-absorbing (reversible) treatment via `non_absorbing="first_entry"` (Eq. 12) or `"effect_stabilization"` (Eq. 13, window `L`). - [BaconDecomposition](https://diff-diff.readthedocs.io/en/stable/api/bacon.html): Goodman-Bacon (2021) decomposition for diagnosing TWFE bias in staggered settings ## Diagnostics and Sensitivity Analysis diff --git a/diff_diff/lpdid.py b/diff_diff/lpdid.py index 66db69b6..9878afa0 100644 --- a/diff_diff/lpdid.py +++ b/diff_diff/lpdid.py @@ -23,6 +23,8 @@ def __init__( alpha: float = 0.05, cluster: Optional[str] = None, rank_deficient_action: str = "warn", + non_absorbing: Optional[str] = None, + stabilization_window: Optional[int] = None, ): self.pre_window = pre_window self.post_window = post_window @@ -33,6 +35,8 @@ def __init__( self.alpha = alpha self.cluster = cluster self.rank_deficient_action = rank_deficient_action + self.non_absorbing = non_absorbing + self.stabilization_window = stabilization_window self._validate_params() self.is_fitted_ = False self.results_: Optional[LPDiDResults] = None @@ -56,6 +60,31 @@ def _validate_params(self) -> None: or (isinstance(self.pmd, int) and not isinstance(self.pmd, bool) and self.pmd > 0) ): raise ValueError("pmd must be None, 'max', or a positive integer") + if self.non_absorbing not in (None, "first_entry", "effect_stabilization"): + raise ValueError( + "non_absorbing must be None (absorbing treatment), 'first_entry' " + "(Dube et al. 2025 Eq. 12), or 'effect_stabilization' (Eq. 13)" + ) + if self.non_absorbing == "effect_stabilization": + if ( + isinstance(self.stabilization_window, bool) + or not isinstance(self.stabilization_window, int) + or self.stabilization_window < 1 + ): + raise ValueError( + "stabilization_window (the paper's L) must be a positive integer when " + "non_absorbing='effect_stabilization'" + ) + elif self.stabilization_window is not None: + raise ValueError( + "stabilization_window only applies when non_absorbing='effect_stabilization'; " + "leave it None for absorbing or first_entry modes" + ) + if self.non_absorbing is not None and self.control_group == "never_treated": + raise ValueError( + "control_group='never_treated' is not supported with a non-absorbing mode " + "(the estimand becomes ambiguous); use control_group='clean' (the default)" + ) def _rhs_column_names(self, covariates=None, ylags=0, dylags=0): rhs_columns = list(covariates or []) @@ -105,12 +134,14 @@ def _prepare_panel( # Absorbing-path validation and entry detection run on the OBSERVED rows, # BEFORE any calendar reindex below (the absorbing fill would otherwise make # the monotonicity check trivially pass, and gap rows carry NaN clusters). - treated_cummax = panel.groupby(unit)["_treated"].cummax() - if (treated_cummax > panel["_treated"]).any(): - raise ValueError( - "LPDiD currently requires an absorbing treatment path " - "(once treated, always treated)" - ) + if self.non_absorbing is None: + treated_cummax = panel.groupby(unit)["_treated"].cummax() + if (treated_cummax > panel["_treated"]).any(): + raise ValueError( + "LPDiD requires an absorbing treatment path (once treated, always " + "treated) unless non_absorbing is set to 'first_entry' or " + "'effect_stabilization' (Dube et al. 2025 Section 4.2)" + ) # Entry = first OBSERVED treated period (documented convention; an unobserved # pre-onset gap is unknowable). For an absorbing path this is min(t | D=1). first_treat = panel.loc[panel["_treated"].eq(1)].groupby(unit)[time].min() @@ -125,6 +156,14 @@ def _prepare_panel( # A gap-free panel skips this entirely and is bit-identical to before. span = panel.groupby(unit)[time].agg(["min", "max", "nunique"]) has_gap = bool((span["nunique"] != (span["max"] - span["min"] + 1)).any()) + if self.non_absorbing is not None and has_gap: + raise ValueError( + "LPDiD non-absorbing modes require gap-free panels within each unit's " + "observed span: the [t-L, t+h] window conditions cannot be verified across " + "an interior time gap. Fill or regularize the panel, or use the absorbing " + "path. (Interior-gap support for non-absorbing treatment is a deferred " + "follow-up.)" + ) if has_gap: panel["_observed"] = True grid = pd.concat( @@ -146,6 +185,25 @@ def _prepare_panel( panel["_first_treat"] = panel[unit].map(first_treat).astype(float).fillna(np.inf) panel["_entry"] = (panel[time] == panel["_first_treat"]).astype(float) + if self.non_absorbing is not None: + # Non-absorbing window operators (Dube et al. 2025 Section 4.2 / online + # Appendix C). Non-absorbing requires a gap-free panel (enforced above), so + # each unit's observed rows ARE its complete calendar grid and the groupby + # shift/cumsum below are calendar-correct. `_treated` here is the GENUINE + # treatment (the absorbing fill above is skipped for non-absorbing), so off + # periods are preserved. Boundary convention (extends Deviation 5): periods + # before a unit's first observed period are untreated with no change, so the + # first first-difference uses fill_value=0 and the mask lookups clamp + # pre-`_unit_min_t` offsets to 0. + prev_treated = panel.groupby(unit)["_treated"].shift(1, fill_value=0) + panel["_delta_d"] = panel["_treated"].astype(float) - prev_treated.astype(float) + # `_switch_cum` = cumulative count of treatment CHANGES (for "no change in + # window" conditions); `_d_cum` = cumulative SUM of D (for "D=0 across window" + # level conditions). Both keyed by calendar time for offset lookups. + panel["_switch_cum"] = panel["_delta_d"].abs().groupby(panel[unit]).cumsum() + panel["_d_cum"] = panel["_treated"].astype(float).groupby(panel[unit]).cumsum() + panel["_unit_min_t"] = panel.groupby(unit)[time].transform("min") + # Premean ("max") baseline = mean of all AVAILABLE strictly-prior outcomes. # It must NOT depend on the base row's own outcome y_t: PMD replaces the # t-1 baseline with the premean of prior periods, and the long difference is @@ -206,14 +264,82 @@ def _clean_control_mask(self, panel: pd.DataFrame, *, time: str, horizon: int) - control_mask &= (panel[time] + horizon).lt(panel["_first_treat"]) return control_mask + def _nonabsorbing_masks(self, sample, panel, *, unit, time, horizon): + """Mode-aware (treated_event, clean_control) masks for non-absorbing treatment. + + Dube et al. (2025) Section 4.2 / online Appendix C. Window conditions over + ``[t-L, t+h]`` are evaluated with offset key-lookups against ``panel``'s + cumulative columns (``_switch_cum`` = running count of treatment CHANGES, for + "no change in window"; ``_d_cum`` = running SUM of D, for "D=0 across window"). + Lookups clamp to 0 for periods before a unit's first observed period (boundary + convention; gap-free panels are enforced in ``_prepare_panel`` so every in-range + period is present). Returns boolean Series indexed like ``sample``. + """ + u = sample[unit].to_numpy() + t = sample[time].to_numpy() + min_t = sample["_unit_min_t"].to_numpy() + th = t + horizon + + switch_by_key = panel.set_index([unit, time])["_switch_cum"] + d_by_key = panel.set_index([unit, time])["_d_cum"] + + def _at(series, times): + keys = pd.MultiIndex.from_arrays([u, np.asarray(times)]) + vals = series.reindex(keys).to_numpy(dtype=float) + return np.where(np.asarray(times) < min_t, 0.0, vals) + + if self.non_absorbing == "first_entry": + # Eq. 12: treated = first entry that stays treated through t+h; control = + # untreated from start through t+h (== absorbing clean control, reused). + treated = sample["_entry"].to_numpy(dtype=float) == 1.0 + if horizon >= 0: + # no exit in (t, t+h] + treated = treated & (_at(switch_by_key, th) - _at(switch_by_key, t) == 0.0) + control = self._clean_control_mask(sample, time=time, horizon=horizon).to_numpy() + else: + # Eq. 13 (effect stabilization, window L). + L = self.stabilization_window + delta_d = sample["_delta_d"].to_numpy(dtype=float) + if horizon >= 0: + # treated = fresh entry untreated in [t-L, t-1] (level sum 0) with no other + # change in (t, t+h]; control = no treatment change in [t-L, t+h] (admits + # stabilized already-treated units as controls). + clean_pre = (_at(d_by_key, t - 1) - _at(d_by_key, t - L - 1)) == 0.0 + treated = (delta_d == 1.0) & clean_pre + treated = treated & (_at(switch_by_key, th) - _at(switch_by_key, t) == 0.0) + control = (_at(switch_by_key, th) - _at(switch_by_key, t - L - 1)) == 0.0 + else: + # Placebo horizons: the long-difference y_{t+h} - y_{t-1} reaches back to + # the target t+h, so the clean window must cover the whole pre-span widened + # to the stabilization length: [t-M, t-1] with M = max(L, -h). Treated must + # be untreated across it (a fresh entry contaminated at an earlier period is + # excluded); controls must have no treatment change across it. + m = max(L, -horizon) + clean_pre = (_at(d_by_key, t - 1) - _at(d_by_key, t - m - 1)) == 0.0 + treated = (delta_d == 1.0) & clean_pre + control = (_at(switch_by_key, t - 1) - _at(switch_by_key, t - m - 1)) == 0.0 + + treated_mask = pd.Series(np.asarray(treated, dtype=bool), index=sample.index) + control_mask = pd.Series(np.asarray(control, dtype=bool), index=sample.index) + # Treated events have a change at t, clean controls have none -> disjoint by + # construction; enforce defensively so a row is never classified as both. + control_mask = control_mask & ~treated_mask + return treated_mask, control_mask + def _common_clean_sample_indicator( self, panel: pd.DataFrame, *, unit: str, time: str, outcome: str, max_post_horizon: int ) -> pd.Series: - common_sample = panel["_entry"].eq(1.0) | self._clean_control_mask( - panel, - time=time, - horizon=max_post_horizon, - ) + if self.non_absorbing is None: + common_sample = panel["_entry"].eq(1.0) | self._clean_control_mask( + panel, + time=time, + horizon=max_post_horizon, + ) + else: + _treated_mask, _control_mask = self._nonabsorbing_masks( + panel, panel, unit=unit, time=time, horizon=max_post_horizon + ) + common_sample = _treated_mask | _control_mask # Fixed composition requires the active baseline AND every post-treatment # target outcome (h = 0..max_post_horizon) to be NON-MISSING for each base # observation -- not merely that the target row exists. This keeps the @@ -279,6 +405,7 @@ def _build_horizon_sample( "_first_treat", "_cluster", "_common_event_ok", + *(["_delta_d", "_unit_min_t"] if self.non_absorbing is not None else []), *rhs_columns, *(absorb or []), ] @@ -305,13 +432,26 @@ def _build_horizon_sample( required_columns = [baseline_column, "_target_outcome", *rhs_columns, *(absorb or [])] sample = sample.dropna(subset=required_columns).copy() - treated_mask = sample["_entry"].eq(1.0) - if self.control_group == "never_treated": - control_mask = sample["_entry"].eq(0.0) & np.isinf(sample["_first_treat"]) + if self.non_absorbing is None: + treated_mask = sample["_entry"].eq(1.0) + if self.control_group == "never_treated": + control_mask = sample["_entry"].eq(0.0) & np.isinf(sample["_first_treat"]) + else: + control_mask = self._clean_control_mask(sample, time=time, horizon=horizon) else: - control_mask = self._clean_control_mask(sample, time=time, horizon=horizon) + treated_mask, control_mask = self._nonabsorbing_masks( + sample, panel, unit=unit, time=time, horizon=horizon + ) sample = sample.loc[treated_mask | control_mask].copy() + if self.non_absorbing is not None: + # The per-horizon clean-treated indicator becomes the regression's treatment + # variable and the treated/control key in every downstream path (estimator, + # RA split, reweight denominators, identification check). For absorbing and + # first_entry this equals the original first-entry _entry on the realized + # sample; under effect_stabilization it also marks re-entry events (which have + # _entry==0) as treated. + sample["_entry"] = treated_mask.loc[sample.index].astype(float) # Fixed composition is a POST-treatment contract: apply it only to post # horizons; pre-treatment placebos use whatever pre-period data exists. if self.no_composition and apply_no_composition and horizon >= 0: @@ -784,6 +924,7 @@ def _build_pooled_sample( "_first_treat", "_cluster", "_common_pooled_ok", + *(["_delta_d", "_unit_min_t"] if self.non_absorbing is not None else []), *rhs_columns, *(absorb or []), ] @@ -813,14 +954,28 @@ def _build_pooled_sample( required_columns = [baseline_column, *target_columns, *rhs_columns, *(absorb or [])] sample = sample.dropna(subset=required_columns).copy() - treated_mask = sample["_entry"].eq(1.0) - if self.control_group == "never_treated": - control_mask = sample["_entry"].eq(0.0) & np.isinf(sample["_first_treat"]) + if self.non_absorbing is None: + treated_mask = sample["_entry"].eq(1.0) + if self.control_group == "never_treated": + control_mask = sample["_entry"].eq(0.0) & np.isinf(sample["_first_treat"]) + else: + control_horizon = max(horizons) if kind == "post" else 0 + control_mask = self._clean_control_mask(sample, time=time, horizon=control_horizon) else: - control_horizon = max(horizons) if kind == "post" else 0 - control_mask = self._clean_control_mask(sample, time=time, horizon=control_horizon) + # Pre pooling reaches back to the MOST NEGATIVE horizon, so the + # non-absorbing clean window must cover it (the effect_stabilization + # placebo window widens to [t - max(L, -h), t-1]); using horizon=0 would + # only clean [t-L, t] and leak prior treatment changes inside the pooled + # pre reach-back. (Absorbing pre uses horizon=0 above: not-yet-treated at t + # implies clean across the whole pre-span, so it is unaffected.) + control_horizon = max(horizons) if kind == "post" else min(horizons) + treated_mask, control_mask = self._nonabsorbing_masks( + sample, panel, unit=unit, time=time, horizon=control_horizon + ) sample = sample.loc[treated_mask | control_mask].copy() + if self.non_absorbing is not None: + sample["_entry"] = treated_mask.loc[sample.index].astype(float) if self.no_composition and kind == "post": sample = sample.loc[sample["_common_pooled_ok"]].copy() sample["_event_time"] = sample[time] @@ -1170,6 +1325,8 @@ def fit( absorb=list(absorb) if absorb else None, ylags=ylags, dylags=dylags, + non_absorbing=self.non_absorbing, + stabilization_window=self.stabilization_window, ) self._fit_meta = {"cluster": cluster, "outcome": outcome, "unit": unit, "time": time} self.is_fitted_ = True @@ -1186,6 +1343,8 @@ def get_params(self) -> Dict[str, Any]: "alpha": self.alpha, "cluster": self.cluster, "rank_deficient_action": self.rank_deficient_action, + "non_absorbing": self.non_absorbing, + "stabilization_window": self.stabilization_window, } def set_params(self, **params: Any) -> "LPDiD": diff --git a/diff_diff/lpdid_results.py b/diff_diff/lpdid_results.py index 2ec79733..f0afbd8b 100644 --- a/diff_diff/lpdid_results.py +++ b/diff_diff/lpdid_results.py @@ -41,6 +41,8 @@ class LPDiDResults: absorb: Optional[List[str]] = None ylags: int = 0 dylags: int = 0 + non_absorbing: Optional[str] = None + stabilization_window: Optional[int] = None # ------------------------------------------------------------------ # internal helpers @@ -142,6 +144,9 @@ def to_dict(self) -> Dict[str, Any]: result["cluster_name"] = self.cluster_name if self.n_clusters is not None: result["n_clusters"] = self.n_clusters + if self.non_absorbing is not None: + result["non_absorbing"] = self.non_absorbing + result["stabilization_window"] = self.stabilization_window result["inference_method"] = "cluster_robust" return result @@ -176,6 +181,13 @@ def _fmt(x: Any, nd: int = 4) -> str: f"Estimand: {self.estimand} Control group: {self.control_group}", f"Base period: {self._base_period_label()} No composition: {self.no_composition}", ] + if self.non_absorbing is not None: + mode = ( + f"effect-stabilization (L={self.stabilization_window})" + if self.non_absorbing == "effect_stabilization" + else "first-entry" + ) + lines.append(f"Treatment path: non-absorbing, {mode} (Dube et al. 2025 Sec. 4.2)") if self.covariates or self.absorb or self.ylags or self.dylags: cov_path = "regression-adjustment" if self.reweight else "direct inclusion" lag_bits = [] diff --git a/docs/api/lpdid.rst b/docs/api/lpdid.rst index e6970161..49e99f80 100644 --- a/docs/api/lpdid.rst +++ b/docs/api/lpdid.rst @@ -1,8 +1,9 @@ Local Projections Difference-in-Differences =========================================== -Local Projections DiD (LP-DiD) estimator for staggered, absorbing-treatment -event studies, from Dube, Girardi, Jordà & Taylor (2025). +Local Projections DiD (LP-DiD) estimator for staggered event studies, from +Dube, Girardi, Jordà & Taylor (2025). Absorbing treatment by default, with +optional non-absorbing (reversible) treatment via ``non_absorbing``. LP-DiD estimates a separate regression at each event-time horizon ``h`` of a long difference of the outcome (``y_{i,t+h} - y_{i,t-1}``) on the @@ -14,10 +15,18 @@ estimand is a strictly non-negatively-weighted average of cohort effects. .. note:: - This release implements the **absorbing-treatment main path**: treatment is - binary and, once switched on, stays on. The estimator rejects panels where a - unit's treatment turns off. Non-absorbing (switch on/off) treatment and - survey-design support are planned follow-ups. Covariates and absorbed fixed + Treatment is binary. By default (``non_absorbing=None``) the estimator + follows the **absorbing main path**: once switched on, treatment stays on, + and panels where a unit's treatment turns off are rejected. Non-absorbing + (switch on/off) treatment is supported via ``non_absorbing="first_entry"`` + (Eq. 12, the effect of entering for the first time and staying treated) or + ``non_absorbing="effect_stabilization"`` (Eq. 13, which requires + ``stabilization_window=L`` and lets units whose treatment has been stable for + at least ``L`` periods serve as clean controls — feasible with few or no + never-treated units). Non-absorbing modes require a gap-free panel within + each unit's observed span and cover the entry-effect estimands; the + Appendix-C exit-event dynamics, R-package parity, and survey-design support + are planned follow-ups. Covariates and absorbed fixed effects are supported; under ``reweight=False`` they enter by direct inclusion, which preserves the non-negative weighting result only under homogeneous covariate effects (online Appendix B.2.2) — the @@ -117,6 +126,17 @@ Premean-differenced base period and fixed-composition sample:: results_pmd = lp_pmd.fit(data, outcome="outcome", unit="unit", time="period", treatment="treated") +Non-absorbing (reversible) treatment — units may switch on and off. Use +``effect_stabilization`` (Eq. 13) when there are few or no never-treated units, +so that units whose treatment has been stable for ``L`` periods can serve as +clean controls:: + + # `treated` here is a non-absorbing 0/1 indicator that can turn on and off. + lp_na = LPDiD(pre_window=3, post_window=4, + non_absorbing="effect_stabilization", stabilization_window=5) + results_na = lp_na.fit(panel, outcome="y", unit="unit", + time="t", treatment="treated") + Comparison with Other Staggered Estimators ------------------------------------------ @@ -133,7 +153,7 @@ Comparison with Other Staggered Estimators - Separate 2x2 DiD aggregation - Impute Y(0) via FE model * - Treatment - - Binary, absorbing (this release) + - Binary; absorbing or non-absorbing (``non_absorbing=``) - Binary, absorbing - Binary, absorbing * - Default estimand diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index 40ac7482..97245e5a 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -1828,6 +1828,18 @@ Weights are **always non-negative** (the central result). Via Frisch-Waugh-Lovel **Covariates (paper Section 4.1):** recommended RA path `beta_{h,x}^{RA} = N_TR^{-1} sum_{TR}[(y_{i,t+h}-y_{i,t-1}) - gamma~^h x_i - delta~_t^h]`, with `gamma~^h, delta~_t^h` from a clean-control-only regression. **PMD base period (Section 3.4):** replace `y_{i,t-1}` with the mean of the last `k` pretreatment periods (`k=t-1` = all); single-cohort `k=t-1` == BJS. **Pooled estimand (Section 3.5):** posttreatment-mean long difference `(1/(H+1)) sum_{h=0}^H y_{i,t+h} - y_{i,t-1}` as the dependent variable. +**Non-absorbing treatment (paper Section 4.2; `non_absorbing=`).** Two distinct entry-effect estimands, both fitted by the same per-horizon LP-DiD regression with mode-specific clean-sample restrictions (`Delta_D_it = D_it - D_{i,t-1}`): + + Eq. 12 (first-time entry, "first_entry"): + treated: D_{i,t+j} = 1 for 0 <= j <= h AND untreated before t (stays treated through t+h) + clean control: D_{i,t+j} = 0 for j <= h (== absorbing Eq. 8 control) + + Eq. 13 (effect stabilization, "effect_stabilization", window L; Assumption 9): + treated: fresh entry at t, D = 0 on [t-L, t-1], no other change in (t, t+h] + clean control: no treatment change in [t-L, t+h] (admits stabilized already-treated units) + +Eq. 12 reuses the absorbing clean control and only restricts the treated set (a unit exiting within the horizon drops out); on an absorbing panel it is numerically identical to the absorbing path. Eq. 13 lets units whose treatment has been stable for at least `L` periods serve as clean controls (Assumption 9: dynamic effects stabilize after `L` periods), making estimation feasible when there are few or no never-treated units (e.g. minimum wage). Non-negative weights (online Appendix C, `omega''`) carry over. Placebo horizons (`h < 0`) use the clean window `[t - max(L, -h), t-1]` so the long difference's reach-back to `t+h` cannot be contaminated. + ### Standard Errors **The paper specifies no SE formula** - Section 1 defers to "standard, well-understood techniques." The reference Stata uses **cluster-robust SEs at the unit level** (`vce(cluster unit)`, footnote 9); pooled / joint tests stack the per-horizon regressions (`suest`). No bootstrap is discussed. Any analytical SE the library ships - and in particular an influence-function cluster variance for the RA path - is therefore an **implementation choice validated against the reference package, not against the paper**, and must be documented under Deviations once implemented (PR-B). @@ -1838,7 +1850,7 @@ Weights are **always non-negative** (the central result). Via Frisch-Waugh-Lovel - **Bias-variance (Sections 3.3, 5.3):** variance weighting (default) -> lower variance, some bias; equal weighting (`reweight`) -> unbiased, higher variance. Variance won at short horizons, equal at long horizons in the paper's simulation. - **PMD vs first-lag (Section 3.4):** PMD gains efficiency under low autocorrelation but can amplify bias if PT holds only in some pretreatment periods; first-lag relies on weaker PT (Marcus & Sant'Anna 2021). Choose the base period ex-ante. - **Covariate-weight positivity (online Appendix B.2):** direct covariate inclusion keeps non-negative weights ONLY under linear + homogeneous covariate effects (B.2.1; main-text Assumption 6); in the general case (B.2.2) weights are not guaranteed positive -> prefer the RA covariate path (the direct path should carry a homogeneity-assumption warning). -- **Non-absorbing (Section 4.2, online Appendix C):** few/no never-treated units handled via the effect-stabilization assumption (Assumption 9, window `L`) with a modified clean-control window (Eq. 13). Two distinct estimands - first-time entry (Eq. 12) and effect-stabilization (Eq. 13). Deferred to a later PR; the absorbing main path rejects non-absorbing input. +- **Non-absorbing (Section 4.2, online Appendix C):** implemented via `non_absorbing="first_entry"` (Eq. 12) and `non_absorbing="effect_stabilization"` (Eq. 13, requires `stabilization_window=L`); the default `non_absorbing=None` keeps the absorbing path and still rejects non-absorbing input. Both modes are **entry-effect** estimands; the Appendix-C exit-event dynamics (`eta_h^{g,n}`, separate switch-off event-studies) are a deferred follow-up. **Boundary convention:** periods before a unit's first observed period are treated as untreated with no change (extends Deviation 5), so window conditions clamp pre-`min_t` offsets to 0 - a unit genuinely treated before the panel starts could be misread as a fresh entry under `effect_stabilization` (an assumption to confirm against R in PR-C2). **Interior gaps** make the `[t-L, t+h]` window conditions unverifiable, so non-absorbing modes require gap-free panels within each unit's observed span and raise otherwise (the absorbing path's interior-gap reindex is a deferred follow-up for non-absorbing). ### Deviations from the paper / from R / library extensions @@ -1847,7 +1859,7 @@ The paper specifies no standard-error formula (Section 1 defers to "standard, we 1. **Note:** Standard errors are **cluster-robust at the unit level by default** - `cluster=None` auto-clusters at the unit identifier and the results record `cluster_name`/`n_clusters` - with a `t(G-1)` reference distribution (G = realized clusters in each horizon's clean-control sample). Matches Stata `lpdid` `vce(cluster unit)`; the paper prescribes no SE. 2. **Note:** The regression-adjustment (RA) covariate path (`reweight=True` with covariates/absorb) reports an **influence-function cluster variance** `sum_c (sum_{i in c} psi_i)^2 / n^2`, in the same family as `ImputationDiD`'s Theorem-3 / BJS variance (see "IF-based variance estimators vs analytical-sandwich estimators" above). Its single Gram inversion is routed through `linalg._rank_guarded_inv` (finite SE on the identified subspace under near-collinearity; NaN at rank 0). Unlike the default/weighted `solve_ols` `hc1`-cluster path - which applies the `(G/(G-1))*((n-1)/(n-k))` finite-sample factor - the RA IF variance carries **no finite-sample factor**, while both paths share the `t(G-1)` reference. **PR-B2 validated this asymmetry as faithful to the authors' own tooling**, not a defect: the no-factor RA convention matches the canonical Stata `teffects ra ... atet vce(cluster)` (the authors' `lpdid_regression_adjustment.do` `margins`/`kmatch` degrees-of-freedom comments prove `teffects` applies neither factor), while the default path matches `feols`/`reghdfe`. The RA *point* estimate is R-anchored to ~1e-13 (full-interaction `i.dtreat##(i.time c.x)` == `teffects` point; `tests/test_methodology_lpdid.py::test_ra_covariate_point`). The RA *standard error* itself has **no runnable R reference** (no R package computes the RA IF variance - `alexCardazzi` uses direct covariate inclusion, not RA; the canonical RA SE is Stata `teffects` only), so it is **pinned** as a documented regression value (`test_ra_covariate_se_regression_pin`) and its calibration is validated by the ungated Monte-Carlo coverage study `benchmarks/python/coverage_lpdid_ra.py` (~0.95 empirical coverage of the true effect at cluster counts G in {30, 100, 300}). 3. **Note:** Direct covariate inclusion (`reweight=False` with covariates/absorb) emits a `UserWarning`: per online Appendix B.2.2 it preserves the non-negative LP-DiD weighting result only under linear and homogeneous covariate effects, so the regression-adjustment path (`reweight=True`) is preferred. -4. **Deviation from R:** Scope - this release implements the **absorbing-treatment main path only** (the estimator raises on non-absorbing input). Non-absorbing treatment (Section 4.2) and survey-design support are deferred to later PRs; the reference `lpdid` packages support non-absorbing treatment. +4. **Deviation from R:** Scope - non-absorbing treatment (Section 4.2) is implemented for the **entry-effect** estimands (`non_absorbing="first_entry"` / `"effect_stabilization"`, PR-C1); the Appendix-C exit-event dynamics and survey-design support remain deferred follow-ups. Validation against the reference R packages (`alexCardazzi/lpdid` `nonabsorbing` / `nonabsorbing_lag`) is the planned PR-C2; until then the non-absorbing standard errors are the same cluster-robust family as the absorbing path (an implementation choice, validated against R in PR-C2). 5. **Note:** LP-DiD's per-unit quantities (outcome lags `ylags`, first-difference lags `dylags`, integer-`pmd` premean baselines, treatment-entry detection) are **calendar** quantities (`t-1`, `t-k`), so the estimator requires integer-valued, globally consecutive `time` labels. A unit with an **interior time gap** is handled by reindexing that unit to its complete interior calendar grid `[min_t, max_t]`, computing the features on the grid, then **restricting back to the observed rows** - so a lag/first-difference spanning a gap is NaN and the observation fails closed (never the previous-*observed* row), and no synthetic gap row enters a regression. A gap-free panel skips this entirely and is bit-identical. **Entry = first OBSERVED treated period** (`min(t | D_it=1)`): an unobserved pre-onset gap cannot move a cohort earlier, the only well-defined convention when the true switch falls in an unobserved period. 6. **Note (pooled estimand):** The pooled pre/post ATT (the headline `results.att` is the pooled-post row) is the **unit-equal-weighted average of each unit-event-time's mean long difference** over the window - `mean_h(y_{i,t+h}) - baseline_{i,t}`, one observation per (unit, event-time), regressed on the treatment-switch indicator with event-time fixed effects on the **fixed-composition** sample (only units observing *every* pooled target, with clean controls required through `max(h)`). This equals the mean of the per-horizon event-study coefficients on a balanced panel. **PR-B2 validated it against the authors' runnable R reference**: the pooled estimand matches the authors' own R pooled recipe (`danielegirardi/lpdid`: a `slider` window-mean minus `y_{t-1}` on the clean-through-window-end sample) to ~1e-13 (`tests/test_methodology_lpdid.py::test_pooled`). A prior version of this note speculated the authors used a horizon-**stacked** pooled regression; the authors' R reference in fact uses this same fixed-composition mean-long-difference, so that speculation was incorrect. Unlike the event-study variants (where `alexCardazzi` is a cross-check gate), pooled is anchored to the authors' recipe **only**: `alexCardazzi`'s pooled uses a **laxer** clean-control window, so it differs and is recorded in the golden `meta` for transparency, not as a parity target. 7. **Deviation from R:** `no_composition` is intentionally more faithful to the paper's fixed-composition intent (Section 3.6) than the R packages: it fixes the realized sample across *all* post horizons (every post coefficient shares one sample, even on unbalanced panels) and excludes cohorts with `p_g > T-H`, whereas `alexCardazzi/lpdid` uses a looser per-horizon sample and a stricter `treat_date < T-H` cutoff. It therefore has **no exact R-package anchor** and is validated by the pure-Python tests in `tests/test_lpdid.py` (the R-parity golden omits it; `alexCardazzi`'s looser-semantics value is recorded in the golden `meta`). @@ -1865,7 +1877,8 @@ The paper specifies no standard-error formula (Section 1 defers to "standard, we - [x] doc-deps.yaml mapping for `diff_diff/lpdid.py` + `lpdid_results.py`; llms.txt / llms-full.txt catalog entries (PR-B1, test-enforced) - [x] B1 pure-Python tests: analytical DGPs + cross-estimator equivalence (CS / BJS / DiD; Cengiz-stacked dropped, documented) + unbalanced / interior-gap / RA-overlap / pmd-missing edge cases (PR-B1) - [x] B2: self-generated R-parity (authors' `danielegirardi/lpdid` recipes + `alexCardazzi/lpdid` cross-check; VW / reweight / pmd / direct / pooled / RA-point to ~1e-12; RA SE pinned + MC-coverage-validated; `no_composition` more paper-faithful than R, B1-tested) (PR-B2) -- [ ] Non-absorbing extension (Section 4.2) - deferred to a later PR +- [x] Non-absorbing extension (Section 4.2): entry-effect estimands - first-entry (Eq. 12) + effect-stabilization (Eq. 13, window `L`) via `non_absorbing`; mode-aware clean-sample masks, `C=0`-below-`min_t` boundary convention, gap-free requirement; pure-Python tests (absorbing reduction, re-entry mechanism, placebo, no-negative-weighting, stabilized-control, DGP recovery) (PR-C1) +- [ ] Non-absorbing R-parity vs `alexCardazzi/lpdid` + exit-event dynamics (Appendix C `eta_h`) - deferred (PR-C2) - [ ] Survey-design support - deferred to a later PR --- diff --git a/tests/test_lpdid.py b/tests/test_lpdid.py index c8ba4f62..2b917970 100644 --- a/tests/test_lpdid.py +++ b/tests/test_lpdid.py @@ -32,6 +32,7 @@ LPDiD, LPDiDResults, ) +from tests.conftest import assert_nan_inference # noqa: E402 # --------------------------------------------------------------------------- @@ -762,7 +763,8 @@ def test_ce3_pmd_single_cohort_equals_bjs_imputation(self): # =========================================================================== class TestLPDiDEdgeCases: def test_rejects_non_absorbing_treatment(self): - # Treatment that turns off must raise (absorbing-path scope). + # Treatment that turns off must raise on the DEFAULT (absorbing) path; a + # non_absorbing mode must accept the same panel (Phase C). df = pd.DataFrame( { "unit": [1, 1, 1, 2, 2, 2], @@ -775,6 +777,14 @@ def test_rejects_non_absorbing_treatment(self): LPDiD(pre_window=1, post_window=1).fit( df, outcome="y", unit="unit", time="time", treatment="treat" ) + # first_entry / effect_stabilization accept non-absorbing input (no raise). + for mode_kw in ( + {"non_absorbing": "first_entry"}, + {"non_absorbing": "effect_stabilization", "stabilization_window": 1}, + ): + LPDiD(pre_window=1, post_window=1, **mode_kw).fit( + df, outcome="y", unit="unit", time="time", treatment="treat", only_event=True + ) def test_no_composition_drops_controls(self): df = _make_linear_panel( @@ -1324,3 +1334,378 @@ def test_rejects_fractional_time_labels(self): ) with pytest.raises(ValueError, match="integer-valued"): LPDiD().fit(df, outcome="y", unit="unit", time="time", treatment="treat") + + +# =========================================================================== +# Non-absorbing treatment (Phase C1; Dube, Girardi, Jorda & Taylor 2025 Sec. 4.2) +# =========================================================================== +def _nonabsorbing_homog_panel( + seed=0, n=200, T=12, tau=2.0, p_switch=0.18, error_sd=0.03, trend=0.4 +): + """Non-absorbing panel with a HOMOGENEOUS, instantly-stabilizing effect ``tau``. + + Each unit follows a random on/off path; ``y = unit_fe + trend*t + tau*D + noise``. + The level jump is the same for every (re-)entry and is constant while treated, so + the per-horizon entry effect equals ``tau`` for both estimands and a stabilized + treated unit differences out (``tau`` cancels) -> a valid Eq.13 control. + """ + rng = np.random.default_rng(seed) + rows = [] + for u in range(n): + ufe = rng.normal(0.0, 1.0) + d = 0 + for t in range(T): + if rng.random() < p_switch: + d = 1 - d + y = ufe + trend * t + tau * d + rng.normal(0.0, error_sd) + rows.append((u + 1, t, d, y)) + return pd.DataFrame(rows, columns=["unit", "time", "treat", "y"]) + + +def _nonabsorbing_het_panel( + seed=0, + n_lo=180, + n_hi=60, + T=12, + tau_lo=1.0, + tau_hi=3.0, + p_switch=0.16, + error_sd=0.03, + trend=0.3, +): + """Two effect levels (``tau_lo`` for ``n_lo`` units, ``tau_hi`` for ``n_hi``). + + True per-horizon effects live in ``[tau_lo, tau_hi]`` (no-negative-weighting hull). + The observation-equally-weighted ATT is the unit-count-weighted mean of the two + levels: ``(n_lo*tau_lo + n_hi*tau_hi)/(n_lo+n_hi)``. + """ + rng = np.random.default_rng(seed) + rows = [] + for u in range(n_lo + n_hi): + ufe = rng.normal(0.0, 1.0) + tau = tau_lo if u < n_lo else tau_hi + d = 0 + for t in range(T): + if rng.random() < p_switch: + d = 1 - d + y = ufe + trend * t + tau * d + rng.normal(0.0, error_sd) + rows.append((u + 1, t, d, y)) + return pd.DataFrame(rows, columns=["unit", "time", "treat", "y"]) + + +def _deterministic_panel(specs, tau=2.0, trend=0.5, ufe_step=10.0): + """Noise-free panel from per-unit binary D paths; ``y = unit_fe + trend*t + tau*D``.""" + rows = [] + for i, d_path in enumerate(specs): + ufe = ufe_step * (i + 1) + for t, d in enumerate(d_path): + rows.append((i + 1, t, int(d), ufe + trend * t + tau * int(d))) + return pd.DataFrame(rows, columns=["unit", "time", "treat", "y"]) + + +_FIT_KW = dict(outcome="y", unit="unit", time="time", treatment="treat") + + +class TestLPDiDNonAbsorbingAPI: + def test_get_params_round_trip(self): + est = LPDiD(non_absorbing="effect_stabilization", stabilization_window=4) + p = est.get_params() + assert p["non_absorbing"] == "effect_stabilization" + assert p["stabilization_window"] == 4 + + def test_rejects_bad_mode(self): + with pytest.raises(ValueError, match="non_absorbing"): + LPDiD(non_absorbing="bogus") + + def test_rejects_missing_or_misused_window(self): + with pytest.raises(ValueError, match="stabilization_window"): + LPDiD(non_absorbing="effect_stabilization") # missing L + with pytest.raises(ValueError, match="stabilization_window"): + LPDiD(non_absorbing="effect_stabilization", stabilization_window=0) # not positive + with pytest.raises(ValueError, match="stabilization_window"): + LPDiD(non_absorbing="first_entry", stabilization_window=3) # L without stab mode + + def test_rejects_never_treated_control_group(self): + with pytest.raises(ValueError, match="never_treated"): + LPDiD(non_absorbing="first_entry", control_group="never_treated") + + def test_set_params_atomic_rollback(self): + est = LPDiD(non_absorbing="effect_stabilization", stabilization_window=2) + with pytest.raises(ValueError): + est.set_params(non_absorbing="first_entry") # would leave a dangling window + # rolled back: still effect_stabilization with its window + assert est.non_absorbing == "effect_stabilization" + assert est.stabilization_window == 2 + + +class TestLPDiDNonAbsorbing: + def test_first_entry_reduces_to_absorbing(self): + # Eq.12 clean control == absorbing clean control; on an absorbing panel every + # treated unit trivially "stays treated", so first_entry == absorbing exactly. + df = make_lpdid_panel(cohorts=(4, 8), n_per_cohort=20, n_never=25, n_periods=12, seed=7) + r0 = LPDiD(pre_window=3, post_window=4, cluster="unit").fit(df, **_FIT_KW) + r1 = LPDiD(pre_window=3, post_window=4, cluster="unit", non_absorbing="first_entry").fit( + df, **_FIT_KW + ) + np.testing.assert_allclose( + r1.event_study["coefficient"].to_numpy(dtype=float), + r0.event_study["coefficient"].to_numpy(dtype=float), + atol=1e-10, + equal_nan=True, + ) + np.testing.assert_allclose( + r1.event_study["se"].to_numpy(dtype=float), + r0.event_study["se"].to_numpy(dtype=float), + atol=1e-10, + equal_nan=True, + ) + assert r1.att == pytest.approx(r0.att, abs=1e-10) + assert r1.se == pytest.approx(r0.se, abs=1e-10) + + def test_effect_stabilization_reduces_on_single_cohort(self): + # Single cohort + never-treated controls is the one configuration where Eq.13 + # admits exactly the absorbing clean-control set, so all three modes coincide. + df = make_lpdid_panel(cohorts=(6,), n_per_cohort=30, n_never=40, n_periods=12, seed=11) + r_ab = LPDiD(pre_window=3, post_window=3, cluster="unit").fit(df, **_FIT_KW) + r_fe = LPDiD(pre_window=3, post_window=3, cluster="unit", non_absorbing="first_entry").fit( + df, **_FIT_KW + ) + r_es = LPDiD( + pre_window=3, + post_window=3, + cluster="unit", + non_absorbing="effect_stabilization", + stabilization_window=12, + ).fit(df, **_FIT_KW) + for r in (r_fe, r_es): + np.testing.assert_allclose( + r.event_study["coefficient"].to_numpy(dtype=float), + r_ab.event_study["coefficient"].to_numpy(dtype=float), + atol=1e-10, + equal_nan=True, + ) + assert r.att == pytest.approx(r_ab.att, abs=1e-10) + + def test_effect_stabilization_uses_reentry_events(self): + # Regression guard for the _entry-keyed misclassification (review critical #2): + # first onsets all exit within the horizon (excluded at h=2); only RE-ENTRY events + # (which have _entry==0) stay treated. effect_stabilization must identify h=2 off + # them; first_entry (first onset only) finds no clean treated event -> NaN. + treated_path = [0, 0, 1, 0, 1, 1, 1, 1, 0, 0] # enter 2, exit 3, re-enter 4, stay..7 + specs = [treated_path] * 6 + [[0] * 10] * 6 + [[1] * 10] * 4 # + never + always + df = _deterministic_panel(specs, tau=2.0) + r_es = LPDiD( + pre_window=1, + post_window=2, + cluster="unit", + non_absorbing="effect_stabilization", + stabilization_window=1, + ).fit(df, only_event=True, **_FIT_KW) + r_fe = LPDiD(pre_window=1, post_window=2, cluster="unit", non_absorbing="first_entry").fit( + df, only_event=True, **_FIT_KW + ) + es_h2 = _event_row(r_es, 2) + assert np.isfinite(es_h2["coefficient"]) + assert es_h2["coefficient"] == pytest.approx(2.0, abs=1e-9) + assert not np.isfinite(_event_row(r_fe, 2)["coefficient"]) + + def test_boundary_clamp_retains_early_controls(self): + # An entry at t=2 with L=4 has a window starting at 2-4-1=-3 < min_t=0; the + # C=0-below-min_t convention must keep the never-treated controls (a NaN-from- + # below-min_t bug would drop them and unidentify the horizon). Review critical #3. + specs = [[0, 0, 1, 1, 1, 1, 1, 1]] * 5 + [[0] * 8] * 8 + df = _deterministic_panel(specs, tau=3.0) + r = LPDiD( + pre_window=1, + post_window=2, + cluster="unit", + non_absorbing="effect_stabilization", + stabilization_window=4, + ).fit(df, only_event=True, **_FIT_KW) + h0 = _event_row(r, 0) + assert np.isfinite(h0["coefficient"]) + assert h0["coefficient"] == pytest.approx(3.0, abs=1e-9) + assert h0["n_obs"] > 5 # 5 treated + retained never-treated controls + + def test_negative_horizon_placebos_are_clean(self): + # Parallel trends hold -> pre-trend (h<-1) coefficients ~ 0 for both modes. A + # broken h<0 mask (post-window term spanning the entry) would distort them. + df = _nonabsorbing_homog_panel(seed=3, n=220, T=12, error_sd=0.02) + for mode_kw in ( + {"non_absorbing": "first_entry"}, + {"non_absorbing": "effect_stabilization", "stabilization_window": 2}, + ): + r = LPDiD(pre_window=3, post_window=3, cluster="unit", **mode_kw).fit(df, **_FIT_KW) + pre = r.event_study[r.event_study["horizon"] < -1] + assert (pre["coefficient"].abs() < 0.1).all(), (mode_kw, pre) + + def test_no_negative_weighting(self): + # The central LP-DiD result: non-negative weights -> the estimate stays inside the + # convex hull [tau_lo, tau_hi] of the true effects (TWFE-style negative weighting + # would push it outside). Holds for both modes and the headline ATT. + df = _nonabsorbing_het_panel(seed=5, error_sd=0.03) + for mode_kw in ( + {"non_absorbing": "first_entry"}, + {"non_absorbing": "effect_stabilization", "stabilization_window": 3}, + ): + r = LPDiD(pre_window=2, post_window=3, cluster="unit", **mode_kw).fit(df, **_FIT_KW) + post = r.event_study[r.event_study["horizon"] >= 0]["coefficient"] + assert (post >= 1.0 - 0.25).all() and (post <= 3.0 + 0.25).all(), (mode_kw, post) + assert 1.0 - 0.25 <= r.att <= 3.0 + 0.25 + + def test_effect_stabilization_admits_stabilized_treated_controls(self): + # No never-treated units exist; early units are treated-from-start (stabilized). + # effect_stabilization uses them as clean controls -> identified; first_entry + # (controls must be untreated through t+h) finds none -> NaN. + specs = [[1] * 10] * 8 + [[0, 0, 0, 0, 1, 1, 1, 1, 1, 1]] * 6 # stabilized + entrant@4 + df = _deterministic_panel(specs, tau=2.0) + r_es = LPDiD( + pre_window=1, + post_window=2, + cluster="unit", + non_absorbing="effect_stabilization", + stabilization_window=2, + ).fit(df, only_event=True, **_FIT_KW) + r_fe = LPDiD(pre_window=1, post_window=2, cluster="unit", non_absorbing="first_entry").fit( + df, only_event=True, **_FIT_KW + ) + es0 = _event_row(r_es, 0) + assert np.isfinite(es0["coefficient"]) + assert es0["coefficient"] == pytest.approx(2.0, abs=1e-9) + assert not np.isfinite(_event_row(r_fe, 0)["coefficient"]) + + def test_reweight_first_entry_reduces_to_absorbing(self): + # The reweight (equal-weight) path also reduces: on an absorbing panel, + # reweighted first_entry == reweighted absorbing (== Callaway-Sant'Anna). + df = make_lpdid_panel(cohorts=(4, 8), n_per_cohort=20, n_never=25, n_periods=12, seed=9) + r_ab = LPDiD(pre_window=2, post_window=4, reweight=True, cluster="unit").fit(df, **_FIT_KW) + r_fe = LPDiD( + pre_window=2, post_window=4, reweight=True, cluster="unit", non_absorbing="first_entry" + ).fit(df, **_FIT_KW) + np.testing.assert_allclose( + r_fe.event_study["coefficient"].to_numpy(dtype=float), + r_ab.event_study["coefficient"].to_numpy(dtype=float), + atol=1e-10, + equal_nan=True, + ) + assert r_fe.att == pytest.approx(r_ab.att, abs=1e-10) + + def test_reweight_recovers_equal_weight_nonabsorbing(self): + # Equal weighting on the heterogeneous panel recovers the observation-equally- + # weighted mean effect = (180*1 + 60*3)/240 = 1.5, distinct from the size of the + # groups' effects taken alone. + df = _nonabsorbing_het_panel( + seed=8, n_lo=180, n_hi=60, tau_lo=1.0, tau_hi=3.0, error_sd=0.02 + ) + r = LPDiD( + pre_window=1, + post_window=2, + reweight=True, + cluster="unit", + non_absorbing="effect_stabilization", + stabilization_window=2, + ).fit(df, only_event=True, **_FIT_KW) + assert r.event_study.set_index("horizon").loc[0, "coefficient"] == pytest.approx( + 1.5, abs=0.2 + ) + + def test_dgp_recovery_both_modes(self): + # Homogeneous instant effect tau=2 -> both modes recover ~2 at every post horizon. + df = _nonabsorbing_homog_panel(seed=2, n=250, T=12, tau=2.0, error_sd=0.04) + for mode_kw in ( + {"non_absorbing": "first_entry"}, + {"non_absorbing": "effect_stabilization", "stabilization_window": 2}, + ): + r = LPDiD(pre_window=2, post_window=3, cluster="unit", **mode_kw).fit(df, **_FIT_KW) + post = r.event_study[r.event_study["horizon"] >= 0] + for _, row in post.iterrows(): + assert row["coefficient"] == pytest.approx(2.0, abs=0.2), (mode_kw, row.to_dict()) + + def test_effect_stabilization_pooled_pre_window_clean(self): + # Pooled-pre placebo must be clean over the DEEPEST reach-back, not just + # [t-L, t-1] (codex P1). With L=1 and pre_window=3 the pooled pre window is + # [-3, -2], so a unit with a treated spell at t-3 contaminates the placebo even + # though [t-1] is clean; the pooled mask must use the most-negative horizon. The + # "spell" entrants (treated at t=3, re-enter at t=5) must be excluded from the + # pooled-pre sample, leaving the clean entrants' ~0 placebo. A horizon=0 pooled + # mask would leak the spell rows and bias the estimate to ~tau/2. + clean_entrant = [0, 0, 0, 0, 0, 1, 1, 1] # enter t=5, no prior spell + spell_entrant = [0, 0, 0, 1, 0, 1, 1, 1] # treated spell at t=3, re-enter t=5 + specs = [clean_entrant] * 5 + [spell_entrant] * 5 + [[0] * 8] * 6 + [[1] * 8] * 4 + df = _deterministic_panel(specs, tau=2.0) + r = LPDiD( + pre_window=3, + post_window=1, + cluster="unit", + non_absorbing="effect_stabilization", + stabilization_window=1, + ).fit(df, **_FIT_KW) + pre = r.pooled.loc[r.pooled["window"] == "pre", "coefficient"].iloc[0] + assert abs(pre) < 1e-6, f"pooled pre placebo contaminated by a prior spell: {pre}" + + def test_interior_gap_raises(self): + df = pd.DataFrame( + { + "unit": [1, 1, 1, 2, 2, 2, 2], + "time": [0, 1, 3, 0, 1, 2, 3], # unit 1 missing t=2 (interior gap) + "y": [1.0, 2, 3, 1, 1, 1, 1], + "treat": [0, 1, 0, 0, 0, 0, 0], + } + ) + with pytest.raises(ValueError, match="gap-free"): + LPDiD(pre_window=1, post_window=1, non_absorbing="first_entry").fit(df, **_FIT_KW) + + def test_one_off_treatment_runs(self): + # Natural-disaster style: a single treated period then off forever. + specs = [[0, 0, 1, 0, 0, 0, 0, 0]] * 8 + [[0] * 8] * 8 + df = _deterministic_panel(specs, tau=2.0) + r = LPDiD( + pre_window=1, + post_window=2, + cluster="unit", + non_absorbing="effect_stabilization", + stabilization_window=1, + ).fit(df, only_event=True, **_FIT_KW) + # h=0 (the treated period itself) is identified off the entry jump. + assert np.isfinite(_event_row(r, 0)["coefficient"]) + + def test_no_clean_control_is_nan_consistent(self): + # All units switch on at t=2 -> no clean controls at that event time; the horizon + # is unidentified and ALL inference fields are NaN together. + specs = [[0, 0, 1, 1, 1, 1]] * 6 + df = _deterministic_panel(specs, tau=2.0) + with pytest.warns(UserWarning): + r = LPDiD( + pre_window=1, + post_window=1, + cluster="unit", + non_absorbing="effect_stabilization", + stabilization_window=1, + ).fit(df, only_event=True, **_FIT_KW) + row = _event_row(r, 0) + assert_nan_inference( + { + "se": row["se"], + "t_stat": row["t_stat"], + "p_value": row["p_value"], + "conf_int": (row["conf_low"], row["conf_high"]), + } + ) + + def test_results_metadata_records_mode(self): + df = _nonabsorbing_homog_panel(seed=1, n=60, T=10) + r = LPDiD( + pre_window=1, + post_window=2, + cluster="unit", + non_absorbing="effect_stabilization", + stabilization_window=3, + ).fit(df, only_event=True, **_FIT_KW) + assert r.to_dict()["non_absorbing"] == "effect_stabilization" + assert r.to_dict()["stabilization_window"] == 3 + assert "non-absorbing" in r.summary() + # absorbing results omit the keys + r_ab = LPDiD(pre_window=1, post_window=2, cluster="unit").fit( + make_lpdid_panel(n_periods=8, seed=1), only_event=True, **_FIT_KW + ) + assert "non_absorbing" not in r_ab.to_dict()