diff --git a/CHANGELOG.md b/CHANGELOG.md index d335b716..3486a2a9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,33 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] ### Added +- **LPDiD non-absorbing SEs now anchored against the authors' Stata `lpdid` package + end-to-end (no library behavior change).** New committed-golden arm + (`benchmarks/stata/generate_lpdid_nonabsorbing_golden.do` → + `benchmarks/data/lpdid_nonabsorbing_stata_golden.json`, consumed by + `tests/test_lpdid_nonabsorbing_stata_parity.py`) running SSC `lpdid` v1.0.3 — + the package builds its own clean samples, unlike the hand-reconstructed + `teffects` RA arm. First external anchors for the non-absorbing REWEIGHTED SE + (previously pinned-only via `RW_SE_PIN` with a documented ~5e-5 feols + weighted-cluster convention gap — the package matches the library convention to + ~1e-9), the non-absorbing POOLED windows (points and SEs), and the Eq. 12 + reweighted point. Scoping: `first_entry` (Eq. 12 == + `nonabsorbing(, firsttreat notyet)`) agrees on every surface of the committed + panel (all event-study horizons incl. placebos + both pooled windows; att + ~1e-14 vw / ~1.3e-8 rw, SE ~6e-17 vw / ~2.3e-9 rw); `effect_stabilization` + (Eq. 13 == `nonabsorbing(L)`) agrees at post horizons + pooled post on a + convention-neutral 47-unit subsample (att ~1.1e-14 / ~1.9e-8, SE ~1.8e-16 / + ~1.6e-9). Three reference-package convention differences were measured, + adjudicated against the paper's Eq. 13 text, and locked by + divergence-documentation gates: the package effectively requires L+1 untreated + periods before re-entry (stricter than Eq. 13's stated levels condition, which + the library implements), admits always-treated units as early-period controls + via missing-lag boundary semantics (paper-silent surface), and builds placebo + samples by recursive lagged intersection rather than the library's backward + window (paper-silent). Eq. 13 full-panel and placebo/pooled-pre inference + remain outside the Stata anchor and are documented as such (REGISTRY `## LPDiD` + Deviation 4; `DEFERRED.md` non-absorbing row narrowed to the RA + + placebo/pooled-pre residuals). - **LWDiD tutorial notebook** (`docs/tutorials/31_lwdid.ipynb`): the replacement for the tutorial withdrawn from PR #588, authored via the numbers-locked workflow (every cited number prototyped in scripts first, notebook assembled and executed once). diff --git a/DEFERRED.md b/DEFERRED.md index 6595a380..5d5047b6 100644 --- a/DEFERRED.md +++ b/DEFERRED.md @@ -65,7 +65,7 @@ exists but parity can't be verified without a local toolchain. | `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 — Stata-`teffects` anchored; optional R follow-up only.** The RA influence-function cluster SE is canonically Stata `teffects ra ... atet vce(cluster)` (no R package computes it — `alexCardazzi/lpdid` does direct covariate inclusion, not RA). **Resolved:** `benchmarks/stata/generate_lpdid_ra_golden.do` reconstructs each horizon's clean sample independently and anchors the SE against `teffects` to ~1e-16 at all 7 horizons (golden `lpdid_ra_stata_golden.json`; `tests/test_lpdid_ra_stata_parity.py`, which also gates the Stata row count vs the library `n_obs`); the point is triple-anchored (Python/R/Stata). Stata is node-locked (not in CI), so the CI-runnable guards remain the regression pin (`RA_SE_PIN`) + MC coverage (`coverage_lpdid_ra.py`). *Optional residual (no longer a blocker):* contribute an RA path to `alexCardazzi/lpdid` so a runnable **R** reference also exists. | `benchmarks/stata/generate_lpdid_ra_golden.do`, `tests/test_lpdid_ra_stata_parity.py` | #B2 · Stata arm | Low | | **`LPDiD` survey scope gaps (PR-D1 deferrals).** Survey support covers the variance-weighted default path only. (a) `survey_design` + `reweight=True` (the equally-weighted / regression-adjustment IF path) is rejected: the weighted RA influence-function variance has **no runnable survey reference** (same class as the RA-SE row above - `survey::svyglm` anchors only the OLS/WLS path). (b) Replicate-weight survey designs (BRR/Fay/JK1/JKn/SDR) and (c) non-pweight (fweight/aweight) types are rejected pending demand. | `lpdid.py`, REGISTRY #8 | PR-D1 | Low | -| **`LPDiD` non-absorbing SE/RA external-reference gap.** The non-absorbing *point* estimands are R-parity-locked (PR-C2: independent `fixest::feols` Eq. 12/13 reconstruction, point+SE ~1e-13/~1e-15 vw; `effect_stabilization` reweighted point + pinned SE), but the authors' canonical non-absorbing SE/RA is Stata `lpdid`/`teffects` only - no faithful R analogue exists (`alexCardazzi/lpdid`'s `nonabsorbing_lag` proved NOT a faithful Eq. 13: off-switch clamp + non-paper boundary/placebo window, diverges ~0.01-0.05 even on a monotone panel; recorded as a divergent reference, not a gate). Revisit if a Stata toolchain or a corrected R package appears - same class as the absorbing RA-SE row above. | `benchmarks/R/generate_lpdid_golden.R`, `tests/test_methodology_lpdid.py` | PR-C2 | Low | +| **`LPDiD` non-absorbing SE/RA external-reference gap.** The authors' canonical non-absorbing SE/RA is Stata `lpdid`/`teffects` only - no faithful R analogue exists (`alexCardazzi/lpdid`'s `nonabsorbing_lag` proved NOT a faithful Eq. 13; recorded as a divergent reference, not a gate). **Resolved (SE half, scoped):** `benchmarks/stata/generate_lpdid_nonabsorbing_golden.do` runs the authors' SSC `lpdid` package end-to-end (golden `lpdid_nonabsorbing_stata_golden.json`; `tests/test_lpdid_nonabsorbing_stata_parity.py`) - Eq. 12 (`nonabsorbing(, firsttreat notyet)` == `first_entry`) anchored on ALL surfaces of the full committed panel (ES incl. placebos + both pooled windows, vw + rw; att ~1e-14/1.3e-8, SE ~6e-17/2.3e-9), Eq. 13 (`nonabsorbing(L)` == `effect_stabilization`) anchored at post horizons + pooled post on a convention-neutral 47-unit subsample (att ~1.1e-14/1.9e-8, SE ~1.8e-16/1.6e-9); first external anchors for the rw SE, the pooled windows, and the Eq. 12 rw point. Three package convention differences measured, paper-adjudicated and divergence-gated (REGISTRY `## LPDiD` Deviation 4). *Residual:* (a) non-absorbing **RA** (`rw` + covariates) - the package's RA engine is `listreg` (not `teffects`) and the committed non-absorbing panel has no covariate; needs a covariate panel + a `listreg` IF-SE convention check. (b) Eq. 13 **placebo/pooled-pre inference** - not anchored by this arm (the package's placebo-window construction differs, convention class 3; a dedicated late-entry/never-treated subsample could in principle align the two placebo constructions - a possible follow-up): the placebo ES SEs remain R-anchored (vw) / `RW_SE_PIN`-pinned (rw, h=-3/-2), and the pooled-pre SE currently has no pin or anchor of any kind (every non-absorbing R-parity fit is event-study-only and the R golden has no pooled block). | `benchmarks/stata/generate_lpdid_nonabsorbing_golden.do`, `tests/test_lpdid_nonabsorbing_stata_parity.py`, `tests/test_methodology_lpdid.py` | PR-C2 · Stata arm | Low | | `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/benchmarks/R/generate_lpdid_golden.R b/benchmarks/R/generate_lpdid_golden.R index b8d71c3a..3ec4906c 100644 --- a/benchmarks/R/generate_lpdid_golden.R +++ b/benchmarks/R/generate_lpdid_golden.R @@ -326,10 +326,12 @@ message(sprintf("Wrote golden: %s", golden_path)) # and ~1e-15 on SE for the variance-weighted variants (reweighted: point ~1e-13, # SE has a small feols-weighted-cluster convention difference, pinned on the # Python side). alexCardazzi::lpdid()'s `nonabsorbing_lag` is NOT a faithful -# Eq.13 (it clamps off-switches via treat_diff[<0]<-0 and uses a non-paper -# boundary/window convention; it diverges ~0.01-0.05 from Eq.13 even on a -# monotone no-off-switch panel) -> recorded in `meta` as a divergent third-party -# reference, NOT a parity gate (the B2 alexCardazzi-pooled precedent). +# Eq.13 (it clamps off-switches via treat_diff[<0]<-0, violating Eq.13's stated +# both-directions control condition; it also holds different boundary/placebo- +# window conventions on surfaces the paper does not specify; it diverges +# ~0.01-0.05 from Eq.13 even on a monotone no-off-switch panel) -> recorded in +# `meta` as a divergent third-party reference, NOT a parity gate (the B2 +# alexCardazzi-pooled precedent). # # APPENDED after the absorbing write_json with its OWN set.seed and distinct # object names; setFixest_ssc is NOT re-called -> the absorbing panel + golden @@ -431,8 +433,9 @@ effect_stab_rw_es <- na_es("eq13", reweight = TRUE) # ---- monotone (no-off-switch) slice: PIN the "alex diverges even without off-switches" # claim with committed evidence. On units whose treatment never decreases, alex's -# off-switch clamp is inert, yet alex still diverges from the paper-faithful Eq.13 (its -# non-paper boundary/window convention), so the recorded max post-horizon |alex - Eq.13| +# off-switch clamp is inert, yet alex still diverges from the library's Eq.13 +# implementation (their boundary/placebo-window conventions differ on surfaces the +# paper does not specify), so the recorded max post-horizon |alex - Eq.13| # is well above 0 -> documents that the divergence is NOT only off-switch handling. na_mono <- na_dt[, if (all(diff(treat) >= 0L)) .SD, by = unit] # drop units with any turn-off na_mode <<- "eq13" diff --git a/benchmarks/README.md b/benchmarks/README.md index 429051e7..e25f0b95 100644 --- a/benchmarks/README.md +++ b/benchmarks/README.md @@ -54,17 +54,19 @@ and CI never needs Stata. Generators live in `stata/` and are run headless: STATA=/Applications/Stata/StataSE.app/Contents/MacOS/stata-se $STATA -b do benchmarks/stata/requirements.do # one-time SSC install $STATA -b do benchmarks/stata/generate_lpdid_ra_golden.do +$STATA -b do benchmarks/stata/generate_lpdid_nonabsorbing_golden.do $STATA -b do benchmarks/stata/generate_imputation_loo_golden.do $STATA -b do benchmarks/stata/generate_etwfe_cs_golden.do $STATA -b do benchmarks/stata/generate_reghdfe_kref_golden.do $STATA -b do benchmarks/stata/generate_lwdid_golden.do # see stata/README.md for its warm-up step ``` -The `LPDiD` arm uses only **native** Stata commands (`teffects`), pinned by +The `LPDiD` RA arm uses only **native** Stata commands (`teffects`), pinned by `version 19`. The `ImputationDiD` arm depends on SSC packages (`did_imputation`/`reghdfe`/`ftools`/`require`), the ETWFE/CS arm on `drdid`/`csdid`/`jwdid`/`hdfe`, the reghdfe K_reference arm on `reghdfe`, -and the LWDiD arm on the authors' `lwdid`; +the LWDiD arm on the authors' `lwdid`, and the LPDiD non-absorbing arm on the +authors' `lpdid` (+ `boottest`/`egenmore`/`listreg`); `version 19` does NOT pin SSC packages (SSC has no version history) — install them once via `requirements.do` (the generators do not auto-install) and each golden records version/checksum @@ -87,6 +89,7 @@ benchmarks/ │ ├── README.md # Stata arm docs │ ├── requirements.do # one-time SSC install (did_imputation etc.) │ ├── generate_lpdid_ra_golden.do # LPDiD RA SE vs teffects ra +│ ├── generate_lpdid_nonabsorbing_golden.do # LPDiD non-absorbing SEs vs the authors' lpdid │ ├── generate_imputation_loo_golden.do # ImputationDiD LOO SE vs did_imputation leaveout │ ├── generate_etwfe_cs_golden.do # ETWFE/CS vs jwdid + csdid (+ subsample ladder) │ ├── generate_reghdfe_kref_golden.do # clustered CR1 K_reference vs reghdfe (disconnected panel) @@ -113,6 +116,7 @@ benchmarks/ | `SyntheticDiD` | `synthdid::synthdid_estimate` | Arkhangelsky et al. (2021) | ✓ Integrated | | `DifferenceInDifferences` | `fixest::feols` | Standard DiD | ✓ Integrated | | `LPDiD` (RA SE) | Stata `teffects ra ... atet` | Dube, Girardi, Jorda & Taylor (2025) | ✓ Integrated | +| `LPDiD` (non-absorbing SEs) | Stata `lpdid, nonabsorbing(...)` (authors' package) | Dube, Girardi, Jorda & Taylor (2025) §4.2 | ✓ Integrated | | `ImputationDiD` (LOO SE) | Stata `did_imputation, leaveout` | Borusyak, Jaravel & Spiess (2024) App. A.9 | ✓ Integrated | | `WooldridgeDiD` / `CallawaySantAnna` | Stata `jwdid` / `csdid` (+ G≈20..500 SE ladder) | Wooldridge (2025) / Callaway & Sant'Anna (2021) | ✓ Integrated | | Clustered CR1 `K_reference` | Stata `reghdfe` + R `fixest` (disconnected-panel arms) | reghdfe/fixest ssc conventions | ✓ Integrated | diff --git a/benchmarks/data/lpdid_nonabsorbing_stata_golden.json b/benchmarks/data/lpdid_nonabsorbing_stata_golden.json new file mode 100644 index 00000000..33e067cc --- /dev/null +++ b/benchmarks/data/lpdid_nonabsorbing_stata_golden.json @@ -0,0 +1,106 @@ +{ + "meta": { + "estimator": "LPDiD non-absorbing SEs - authors' Stata lpdid package, end-to-end", + "generator": "benchmarks/stata/generate_lpdid_nonabsorbing_golden.do", + "source_panel": "benchmarks/data/lpdid_nonabsorbing_panel.csv", + "source_sha256": "71af98c96dec08abc281002b7b88df74033dac7d0ff7707eff9f3ee0c30a3964", + "cmd": "lpdid y, unit(unit) time(time) treat(treat) pre_window(3) post_window(4) nograph + nonabsorbing(, firsttreat notyet) [rw] on the full panel; nonabsorbing(3) [rw] on the subsample and (only_event, vw) on the full panel", + "mapping": "nonabsorbing(, firsttreat notyet) == first_entry; nonabsorbing(L) == effect_stabilization", + "drop_rule": "min(treat)==1 | any(dD==1 & L3.dD==-1)", + "point_anchor": "first_entry vw tau0 gated in-.do vs the committed R golden first_entry[0] att (R-anchored); effect_stab_sub vw tau0 vs the library point (library-anchored, informational); the Python test is authoritative for every surface", + "se_convention": "package reghdfe vce(cluster unit) CR1; the library's cluster SE matches at ~1e-16 (vw) / ~1e-9 (rw) on identical samples", + "convention_notes": "effect_stab full-panel and placebo/pooled-pre rows are DIVERGENCE DOCUMENTATION, not parity: the package admits always-treated units at early t via missing-lag semantics, requires L+1 untreated periods before re-entry (stricter than JAE Eq. 13's levels condition), and builds placebo samples by recursive lagged intersection; see REGISTRY ## LPDiD Deviation 4", + "ssc_versions": { + "lpdid": "lpdid program, version 1.0.3 [checksum:3687664943 len:48365]", + "reghdfe": "version 6.13.1 10Jan2026 [checksum:973401092 len:30978]", + "ftools": "version 2.50.0 09jan2026 [checksum:3739654356 len:1428]", + "require": "version 1.3.1 19sep2023 [checksum:1262842970 len:34415]", + "boottest": "boottest 4.5.3 5 June 2026 [checksum:2334657244 len:63528]", + "listreg": "version 1.0.3 26apr2025 Ben Jann [checksum:1074534106 len:23815]", + "egenmore": "1.1.1 NJC 19 March 2006 [checksum:2753673312 len:1164]", + "egenmore_gclsst": "1.0.0 NJC 11 January 2000 [checksum:885642367 len:916]" + }, + "dropped_units": [24, 25, 27, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40], + "stata_edition": "SE", + "stata_version": 19.0, + "pre_window": 3, + "post_window": 4, + "stabilization_window": 3, + "n_rows_full": 840, + "n_rows_sub": 658 + }, + "first_entry": { + "vw": { + "es": { + "-3": {"att": -0.11058688194906449, "se": 0.083864694854006674, "N": 300}, + "-2": {"att": -0.028508009716748425, "se": 0.059950164863315426, "N": 350}, + "0": {"att": 1.8468153374570833, "se": 0.091536864444450375, "N": 400}, + "1": {"att": 1.9174189184921464, "se": 0.082514757717069953, "N": 350}, + "2": {"att": 1.9541562370690497, "se": 0.095995993613636657, "N": 298}, + "3": {"att": 1.8430002097850104, "se": 0.093215506869608694, "N": 246}, + "4": {"att": 1.8651045890489937, "se": 0.10344600862300475, "N": 211} + }, + "pooled": { + "pre": {"att": -0.083163290352517599, "se": 0.060216134633611737, "N": 300}, + "post": {"att": 1.8463774549223959, "se": 0.082239499731901228, "N": 211} + } + }, + "rw": { + "es": { + "-3": {"att": -0.10982753532939604, "se": 0.085127249305807715, "N": 180}, + "-2": {"att": -0.027473431971800764, "se": 0.059909786016452311, "N": 230}, + "0": {"att": 1.8483385269860748, "se": 0.093169339989145811, "N": 230}, + "1": {"att": 1.9098467525590781, "se": 0.082755294278197722, "N": 203}, + "2": {"att": 1.9542651398293733, "se": 0.095993108882509365, "N": 173}, + "3": {"att": 1.844237043489463, "se": 0.094299729259659446, "N": 123}, + "4": {"att": 1.8646288450757691, "se": 0.10524070978446542, "N": 118} + }, + "pooled": { + "pre": {"att": -0.080380126038013866, "se": 0.061217606891681819, "N": 180}, + "post": {"att": 1.8500925875939953, "se": 0.082980268093888879, "N": 118} + } + } + }, + "effect_stab_sub": { + "n_rows": 658, + "vw": { + "es": { + "-3": {"att": -0.09724719295046115, "N": 342}, + "-2": {"att": -0.007899289624071772, "N": 408}, + "0": {"att": 1.8115854492490506, "se": 0.079386732581726835, "N": 484}, + "1": {"att": 1.9300421146595237, "se": 0.072501471398279477, "N": 414}, + "2": {"att": 1.8965697184052384, "se": 0.084395761100608244, "N": 348}, + "3": {"att": 1.879554866833909, "se": 0.076256728580062952, "N": 279}, + "4": {"att": 1.893416496899663, "se": 0.093768556793994276, "N": 230} + }, + "pooled": { + "pre": {"att": -0.064656623253798173, "N": 342}, + "post": {"att": 1.8538432413410768, "se": 0.072566538322928476, "N": 230} + } + }, + "rw": { + "es": { + "-3": {"att": -0.09552710717227407, "N": 173}, + "-2": {"att": -0.0054487064168838913, "N": 220}, + "0": {"att": 1.813509946055728, "se": 0.081752780459597857, "N": 329}, + "1": {"att": 1.9275161589672816, "se": 0.073720330196303835, "N": 305}, + "2": {"att": 1.9026466911687787, "se": 0.085905795413003486, "N": 280}, + "3": {"att": 1.8767692749265028, "se": 0.078313461669530404, "N": 191}, + "4": {"att": 1.8899064444115865, "se": 0.096189906878782108, "N": 151} + }, + "pooled": { + "pre": {"att": -0.061093488880953395, "N": 173}, + "post": {"att": 1.8560600200442186, "se": 0.074103122657375023, "N": 151} + } + } + }, + "effect_stab_full_vw": { + "es": { + "0": {"att": 1.8496401510611364, "N": 628}, + "1": {"att": 1.9364029248640582, "N": 544}, + "2": {"att": 1.9088763162649005, "N": 463}, + "3": {"att": 1.9438383345066086, "N": 380}, + "4": {"att": 1.9261821895864399, "N": 320} + } + } +} diff --git a/benchmarks/stata/README.md b/benchmarks/stata/README.md index cae1ba14..81d8b04a 100644 --- a/benchmarks/stata/README.md +++ b/benchmarks/stata/README.md @@ -1,7 +1,8 @@ # Stata parity benchmarks -Stata golden generators live here (five arms: LPDiD `teffects ra`, -ImputationDiD leave-one-out, ETWFE/CS, reghdfe K_reference, LWDiD); the +Stata golden generators live here (six arms: LPDiD `teffects ra`, LPDiD +non-absorbing `lpdid`, ImputationDiD leave-one-out, ETWFE/CS, reghdfe +K_reference, LWDiD); the pattern mirrors `benchmarks/R/` (a `generate_*` script writes a committed golden JSON that a skip-guarded test reads, so CI never needs Stata). Stata is node-locked single-user, so — exactly like the R arm — goldens are @@ -118,6 +119,88 @@ sample. No timestamp — the golden regenerates byte-identically. --- +# `lpdid` parity for the LPDiD non-absorbing modes + +`benchmarks/stata/generate_lpdid_nonabsorbing_golden.do` produces +`benchmarks/data/lpdid_nonabsorbing_stata_golden.json`, consumed by +`tests/test_lpdid_nonabsorbing_stata_parity.py`. Unlike the `teffects` arm above, this +arm runs the authors' SSC **`lpdid`** package (Busch & Girardi, in collaboration with +Dube, Jorda & Taylor) **end-to-end** — the package builds its own clean samples — on +the committed `benchmarks/data/lpdid_nonabsorbing_panel.csv` (owner: +`benchmarks/R/generate_lpdid_golden.R`; 60 units x 14 periods, balanced, gap-free). +First external anchors: the non-absorbing **reweighted SE** (previously pinned-only +via `RW_SE_PIN`), the non-absorbing **pooled windows** (points and SEs), and the +Eq. 12 reweighted point. + +## Option-to-mode mapping + +| Stata `lpdid` | diff-diff | +|---|---| +| `nonabsorbing(, firsttreat notyet)` | `LPDiD(non_absorbing="first_entry")` (Eq. 12) | +| `nonabsorbing(L)` | `LPDiD(non_absorbing="effect_stabilization", stabilization_window=L)` (Eq. 13) | + +## Parity scoping and the convention-neutral subsample + +Eq. 12 agrees on EVERY surface of the full panel (all event-study horizons incl. +placebos + both pooled windows, vw and rw; att ~1e-14/1.3e-8, SE ~6e-17/2.3e-9). +Eq. 13 agrees at post horizons + pooled post only, and only on a 47-unit subsample, +because the package's sample construction differs on three measured conventions +(REGISTRY `## LPDiD` Deviation 4): (1) missing-lag boundary semantics admit +always-treated units as early-`t` controls (paper-silent surface); (2) its +switch-free lag window effectively requires `L+1` untreated periods before re-entry +where the paper's Eq. 13 states a levels condition over `[t-L, t-1]` (the package is +stricter than the paper); (3) placebo samples are built by recursive lagged +intersection of its `CCS_0` rather than the library's backward window (paper-silent). +The subsample removes what (1)+(2) admit differently — always-treated units +(`min(treat)==1`: 31–40) plus exact-`L`-respell units (any `t` with +`dD_t==1 & dD_{t-L}==-1`: 24, 25, 27) — making the Eq. 13 post samples row-identical +(658 rows). Convention (3) is not neutralized by THIS subsample (a dedicated +late-entry/never-treated subsample could in principle align the two placebo +constructions — a possible follow-up), so Eq. 13 +placebo/pooled-pre rows and the Eq. 13 full-panel event study are recorded as +measured DIVERGENCE documentation (att + obs only) and gated as such +(divergence-floor + sample-size-mismatch asserts), never as parity. + +## Regenerating + +```bash +# once: install lpdid + deps (boottest, egenmore, listreg) from SSC +/Applications/Stata/StataSE.app/Contents/MacOS/stata-se -b do benchmarks/stata/requirements.do +# from the repo root +/Applications/Stata/StataSE.app/Contents/MacOS/stata-se -b do \ + benchmarks/stata/generate_lpdid_nonabsorbing_golden.do +# batch mode ALWAYS exits 0 — verify the log has no Stata errors: +grep -E '^r\([0-9]+\);' generate_lpdid_nonabsorbing_golden.log +python -m json.tool benchmarks/data/lpdid_nonabsorbing_stata_golden.json > /dev/null +``` + +The generator fails closed (`exit 111`) if any run-time dependency is missing and +records every one in `meta.ssc_versions` (guard set == version set; both guarded +egenmore files carry their own entry — `egenmore` via `_gfilter.ado`, the file the +pooled spec executes, and `egenmore_gclsst` via `_gclsst.ado`, lpdid's startup +which-check). Each version string also embeds the ado file's `checksum:`/`len:` so +a same-version-string upstream replacement still moves the drift metadata. + +## JSON schema + +```json +{"meta": {"...provenance...", "mapping": "...", "drop_rule": "...", + "dropped_units": [...], "source_sha256": "...", "ssc_versions": {...}}, + "first_entry": {"vw": {"es": {"-3".."4": {"att", "se", "N"}}, + "pooled": {"pre": {...}, "post": {...}}}, "rw": {...}}, + "effect_stab_sub": {"n_rows": 658, + "vw": {"es": {"post rows: {att, se, N}; pre rows: {att, N}"}, + "pooled": {"pre": {"att", "N"}, "post": {"att", "se", "N"}}}, + "rw": {...}}, + "effect_stab_full_vw": {"es": {"0".."4": {"att", "N"}}}} +``` + +Divergent rows deliberately omit `se` — every recorded numeric value is covered by a +Python gate, and no gate anchors those SEs. No timestamp — byte-identical +regeneration under fixed Stata + fixed installed SSC versions. + +--- + # `did_imputation, leaveout` parity for ImputationDiD LOO SE `benchmarks/stata/generate_imputation_loo_golden.do` produces @@ -134,11 +217,12 @@ The A.9 leave-one-out (LOO) variance has **no runnable R reference** — R internal psi-identity + hand-calc + MC coverage. The authors' own Stata `did_imputation` (Borusyak) ships the same option (`leaveout`); this arm turns it into a measured anchor. -## SSC dependence (vs the native-`teffects` LPDiD arm) +## SSC dependence (vs the native-`teffects` LPDiD RA arm) -Unlike the native-`teffects` LPDiD arm, `did_imputation` is an SSC package with a +Unlike the native-`teffects` LPDiD RA arm, `did_imputation` is an SSC package with a dependency chain `did_imputation → reghdfe → require + ftools`, none pinned by -`version 19`. (The ETWFE/CS and LWDiD arms are SSC-dependent in the same way.) The generator does **not** install them — run +`version 19`. (The ETWFE/CS, LWDiD and LPDiD non-absorbing arms are SSC-dependent +in the same way.) The generator does **not** install them — run `benchmarks/stata/requirements.do` once first — and it records each package's version in `meta.ssc_versions` (the `*!` ado header line) so drift is detectable. Byte-identical regeneration is therefore scoped to a fixed Stata + fixed installed SSC versions. @@ -392,10 +476,10 @@ explicit missing branch). ## Known constraints - **Batch mode always exits 0**, even on a hard error (`r(NNN);`). Never trust the - shell exit code — parse the `.log` for `^r\([0-9]+\);`. The LPDiD, ImputationDiD - and LWDiD generators also run informational in-`.do` smoke gates (LPDiD 1e-8, - ImputationDiD 1e-6, LWDiD 1e-10/1e-6 on its SMALL-N blocks; the LWDiD Walmart - block gates schema/cardinality/nonmissing only) that surface as `r(9);` on a + shell exit code — parse the `.log` for `^r\([0-9]+\);`. The LPDiD (RA arm 1e-8; + non-absorbing arm 1e-6), ImputationDiD (1e-6) and LWDiD generators (1e-10/1e-6 on + its SMALL-N blocks; the LWDiD Walmart block gates schema/cardinality/nonmissing + only) also run informational in-`.do` smoke gates that surface as `r(9);` on a gross bug; the Python parity tests are authoritative. - **`c(flavor)` misreports the edition** as `IC` on StataSE, and `c(edition)` is unreliable. The generator derives the edition from the `c(MP)` / `c(SE)` 0/1 flags, @@ -403,7 +487,7 @@ explicit missing branch). eagerly, so it must not be referenced). `"SE"` is simply the committed golden's current value. - **SSC has no version history.** `ssc install` always fetches latest and there is - no lockfile / archive to pin against. The LPDiD arm is exempt (`teffects` is native, - pinned by `version 19`); the SSC arms (ImputationDiD, ETWFE/CS, LWDiD) record their - SSC package versions in `meta.ssc_versions` so drift is at least detectable — new - SSC arms should do the same. + no lockfile / archive to pin against. The LPDiD RA arm is exempt (`teffects` is + native, pinned by `version 19`); the SSC arms (ImputationDiD, ETWFE/CS, LWDiD, + LPDiD non-absorbing) record their SSC package versions in `meta.ssc_versions` so + drift is at least detectable — new SSC arms should do the same. diff --git a/benchmarks/stata/generate_lpdid_nonabsorbing_golden.do b/benchmarks/stata/generate_lpdid_nonabsorbing_golden.do new file mode 100644 index 00000000..9852aeeb --- /dev/null +++ b/benchmarks/stata/generate_lpdid_nonabsorbing_golden.do @@ -0,0 +1,465 @@ +*! Golden generator: LPDiD NON-ABSORBING SEs vs the authors' Stata `lpdid` package. +*! +*! Purpose +*! Anchor the library's non-absorbing LP-DiD modes (Dube, Girardi, Jorda & +*! Taylor 2025, JAE Eq. 12 / Eq. 13) against the authors' reference +*! implementation, SSC `lpdid` (Busch & Girardi, in collaboration with Dube, +*! Jorda and Taylor), END-TO-END - the package builds its own clean samples, +*! unlike the sibling `teffects` arm which reconstructs them by hand. This is +*! the first external anchor for the non-absorbing REWEIGHTED SE (previously +*! pinned-only via RW_SE_PIN with a ~5e-5 feols weighted-cluster convention +*! gap), for the non-absorbing POOLED windows (points and SEs), and for the +*! Eq. 12 reweighted point. +*! +*! Mapping (pinned in meta and gated by the consuming test) +*! lpdid, nonabsorbing(, firsttreat notyet) == LPDiD(non_absorbing="first_entry") +*! lpdid, nonabsorbing(L) == LPDiD(non_absorbing="effect_stabilization", +*! stabilization_window=L) +*! +*! Arms +*! A: Eq. 12 (firsttreat notyet), FULL committed panel, vw + rw: +*! all ES horizons (incl. placebos) + pooled Pre/Post. Full parity surface. +*! B: Eq. 13 (nonabsorbing(3)), convention-neutral 47-unit SUBSAMPLE, vw + rw: +*! post horizons + pooled Post are parity surfaces; pre rows + pooled Pre +*! are recorded (att + obs only) as measured DIVERGENCE documentation -- +*! the package builds placebo samples by recursive lagged intersection of +*! CCS_0 while the library uses the backward window [t-max(L,-h), t-1] +*! (a paper-silent surface; see REGISTRY ## LPDiD Deviation 4). +*! C: Eq. 13, FULL panel, vw only, event study only (att + obs): measured +*! divergence documentation for the two sample-admission convention +*! differences (pre-panel boundary handling; exact-L re-entry spells). +*! +*! Subsample drop rule (computable identically here and in pandas; gate 5 of the +*! consuming test asserts both sub-classes): +*! drop units with min(treat)==1 (always-treated: 31..40) +*! drop units with any(dD==1 & L3.dD==-1), dD=D.treat (exact-L respell: 24,25,27) +*! -> 47 units / 658 rows. +*! +*! Consuming test +*! tests/test_lpdid_nonabsorbing_stata_parity.py +*! +*! Outputs (checked into the repo) +*! benchmarks/data/lpdid_nonabsorbing_stata_golden.json +*! +*! Usage (run from the repo root; install SSC deps ONCE via +*! benchmarks/stata/requirements.do first - this generator fails closed) +*! /Applications/Stata/StataSE.app/Contents/MacOS/stata-se -b do \ +*! benchmarks/stata/generate_lpdid_nonabsorbing_golden.do +*! Then confirm the log is clean: +*! grep -E '^r\([0-9]+\);' generate_lpdid_nonabsorbing_golden.log +*! (Stata batch mode ALWAYS exits 0, even on error - never trust the exit code.) +*! +*! Notes +*! - Reads (does NOT regenerate) benchmarks/data/lpdid_nonabsorbing_panel.csv, +*! whose sole owner is benchmarks/R/generate_lpdid_golden.R (Section 5). +*! 60 units x 14 periods, balanced, gap-free. +*! - SSC packages are unpinned; meta.ssc_versions records the installed +*! versions of every fail-closed dependency as the only drift signal. +*! - In-.do point gates are informational; the Python test is authoritative. +*! Arm A vw tau0 is gated against the committed R golden's first_entry["0"] +*! att (R-anchored); Arm B vw tau0 against the library-computed point +*! (library-anchored - no R value exists for the subsample). + +version 19 +clear all +set more off +set type double + +* Format a scalar as a JSON number at round-trip-exact precision. Stata's %21.17g +* renders |x|<1 as ".455"/"-.119" (leading dot); JSON requires a leading 0, so we +* patch ".x" -> "0.x" and "-.x" -> "-0.x". Every value this arm records is +* expected FINITE, so a missing scalar fails the generator loudly (string(.) is +* "." which the leading-dot patch would turn into invalid JSON "0.") instead of +* silently corrupting the golden. Returns r(s). +capture program drop _jnum +program define _jnum, rclass + args x fmt + if missing(`x') { + display as error "_jnum: missing value where a finite scalar was expected" + exit 9 + } + if "`fmt'" == "" local fmt "%21.17g" + local s = strtrim(string(`x', "`fmt'")) + if substr("`s'", 1, 1) == "." local s = "0" + "`s'" + else if substr("`s'", 1, 2) == "-." local s = "-0" + substr("`s'", 2, .) + return local s "`s'" +end + +* Installed-version capture for an unpinned SSC dependency. Widened from the +* ETWFE arm's _adover (generate_etwfe_cs_golden.do): this arm's own packages +* defeat a "*!"-only parse - lpdid.ado's header is a PLAIN-`*` comment +* ("* lpdid program, version 1.0.3"; no "*!" anywhere) and _gfilter.ado's "*!" +* line carries no "version" token ("*! 1.1.1 NJC 19 March 2006"). Preference +* order: a "*!" line containing "version"; else the FIRST "*!" line (boottest's +* "*! boottest 4.5.3 ..." must beat the GPL boilerplate "version 3 of the +* License" plain-comment line below it); else the first PLAIN comment line +* containing "version" (lpdid has no "*!" at all - its header is +* "* lpdid program, version 1.0.3"); else the checksum alone. The file's +* checksum + length are ALWAYS appended to whatever version text was found +* (SSC has no immutable archive, so a same-version-string upstream edit would +* otherwise be invisible in the drift metadata). Takes the FILE basename to +* probe (may differ from the package name, e.g. egenmore -> _gfilter). +capture program drop _adover +program define _adover, rclass + args probefile + capture findfile `probefile'.ado + if _rc { + return local v "MISSING" + exit + } + local fn = r(fn) + tempname vh + local ver "" + local cver "" + local first "" + local n = 0 + file open `vh' using "`fn'", read text + file read `vh' line + while r(eof) == 0 & `n' < 15 { + if substr(`"`macval(line)'"', 1, 2) == "*!" { + local body = strtrim(substr(`"`macval(line)'"', 3, .)) + if `"`first'"' == "" local first `"`body'"' + if strpos(lower(`"`body'"'), "version") > 0 { + local ver `"`body'"' + continue, break + } + } + else if substr(strtrim(`"`macval(line)'"'), 1, 1) == "*" { + local body = strtrim(substr(strtrim(`"`macval(line)'"'), 2, .)) + if `"`cver'"' == "" & strpos(lower(`"`body'"'), "version") > 0 { + local cver `"`body'"' + } + } + local ++n + file read `vh' line + } + file close `vh' + if `"`ver'"' == "" local ver `"`first'"' + if `"`ver'"' == "" local ver `"`cver'"' + * ALWAYS append the file checksum + length: SSC has no immutable archive, so a + * same-version-string replacement of an ado would otherwise be invisible in + * the golden's drift metadata (local review R2). + capture checksum "`fn'" + if _rc { + display as error "_adover: checksum failed for `fn' - the drift signal cannot be recorded" + exit 9 + } + local csum "checksum:`r(checksum)' len:`r(filelen)'" + if `"`ver'"' == "" { + local ver `"`csum'"' + } + else { + local ver `"`ver' [`csum']"' + } + * JSON-sanitize: drop double-quotes/backslashes/tabs, truncate (cap covers the + * version text plus the always-appended [checksum:... len:...] suffix). + local ver = subinstr(`"`ver'"', `"""', "'", .) + local ver = subinstr(`"`ver'"', "\", "/", .) + local ver = subinstr(`"`ver'"', char(9), " ", .) + if length(`"`ver'"') > 160 local ver = substr(`"`ver'"', 1, 160) + return local v `"`ver'"' +end + +* ------------------------------------------------------------------------------ +* Fail closed if ANY run-time dependency is missing: a silent skip would emit a +* truncated golden, and a missing egenmore helper would otherwise surface as an +* opaque mid-run error inside lpdid. (_gfilter is executed by every lpdid pooled +* spec; _gclsst is lpdid's own startup which-check; both ship in SSC egenmore.) +* ------------------------------------------------------------------------------ +foreach p in lpdid reghdfe ftools require boottest listreg _gfilter _gclsst { + capture which `p' + if _rc { + display as error "Missing SSC dependency `p'. Run benchmarks/stata/requirements.do first." + exit 111 + } +} + +* Record installed SSC versions for every fail-closed dependency (guard set == +* version set, matching the ETWFE-arm contract). Both guarded egenmore files are +* recorded: `egenmore` via _gfilter.ado (executed by the pooled spec) and +* `egenmore_gclsst` via _gclsst.ado (lpdid's startup which-check) - a change to +* either file must move the drift metadata. +foreach p in lpdid reghdfe ftools require boottest listreg { + _adover `p' + local v_`p' = r(v) +} +_adover _gfilter +local v_egenmore = r(v) +_adover _gclsst +local v_egenmore_gclsst = r(v) + +* ------------------------------------------------------------------------------ +* Informational in-.do point pins (fatal asserts at reldif < 1e-6, loose enough +* to survive SSC solver drift; the Python test is the authoritative gate). +* rpt_fe : committed R golden first_entry["0"] att (R-anchored; +* benchmarks/data/lpdid_nonabsorbing_golden.json, digits=12) +* rpt_es : library-computed Eq. 13 subsample vw tau0 (library-anchored, +* informational - no R value exists for the subsample) +* ------------------------------------------------------------------------------ +scalar rpt_fe = 1.846815337457 +scalar rpt_es = 1.8115854492 + +* ------------------------------------------------------------------------------ +* Load the committed panel and confirm the expected schema/shape (balanced, +* gap-free; the library's non-absorbing modes require it, so the golden must be +* generated from the same grid). +* ------------------------------------------------------------------------------ +import delimited using "benchmarks/data/lpdid_nonabsorbing_panel.csv", varnames(1) clear +confirm numeric variable unit time treat y +quietly count +assert r(N) == 840 // 60 units x 14 periods +bysort unit: gen _n_obs = _N +assert _n_obs == 14 // balanced: every unit observed at every period +drop _n_obs + +* ------------------------------------------------------------------------------ +* Arm A: Eq. 12 (first_entry), FULL panel, vw + rw. +* e(results) rows: pre3 pre2 pre1 tau0..tau4 (pre1 = zero reference, skipped). +* Columns: coefficient se t p ci_low ci_high obs -> capture 1, 2, 7. +* ------------------------------------------------------------------------------ +foreach mode in vw rw { + local rwopt = cond("`mode'" == "rw", "rw", "") + preserve + lpdid y, unit(unit) time(time) treat(treat) pre_window(3) post_window(4) /// + nonabsorbing(, firsttreat notyet) `rwopt' nograph + matrix J = e(results) + matrix P = e(pooled_results) + foreach h of numlist 3 2 { + scalar fe_`mode'_att_m`h' = J[4 - `h', 1] + scalar fe_`mode'_se_m`h' = J[4 - `h', 2] + scalar fe_`mode'_N_m`h' = J[4 - `h', 7] + } + foreach h of numlist 0 1 2 3 4 { + scalar fe_`mode'_att_`h' = J[4 + `h', 1] + scalar fe_`mode'_se_`h' = J[4 + `h', 2] + scalar fe_`mode'_N_`h' = J[4 + `h', 7] + } + foreach w in Pre Post { + local r = cond("`w'" == "Pre", 1, 2) + local wl = lower("`w'") + scalar fe_`mode'_att_p`wl' = P[`r', 1] + scalar fe_`mode'_se_p`wl' = P[`r', 2] + scalar fe_`mode'_N_p`wl' = P[`r', 7] + } + restore +} +* in-.do smoke gate: Arm A vw tau0 vs the committed R golden value. +assert reldif(fe_vw_att_0, rpt_fe) < 1e-6 + +* ------------------------------------------------------------------------------ +* Arm B: Eq. 13 (effect stabilization, L=3), convention-neutral SUBSAMPLE, +* vw + rw. Drop rule (see header); the variable is named dD, never bare D, so it +* cannot shadow the difference-operator prefix. +* ------------------------------------------------------------------------------ +preserve + tsset unit time + gen dD = D.treat + bysort unit: egen alw = min(treat) + gen respell = (dD == 1 & L3.dD == -1) + bysort unit: egen dropu = max(respell) + drop if alw == 1 | dropu == 1 + quietly count + assert r(N) == 658 // 47 units x 14 periods + drop dD alw respell dropu + foreach mode in vw rw { + local rwopt = cond("`mode'" == "rw", "rw", "") + lpdid y, unit(unit) time(time) treat(treat) pre_window(3) post_window(4) /// + nonabsorbing(3) `rwopt' nograph + matrix J = e(results) + matrix P = e(pooled_results) + * divergent placebo rows: att + obs only (documentation, gate 6b) + foreach h of numlist 3 2 { + scalar es_`mode'_att_m`h' = J[4 - `h', 1] + scalar es_`mode'_N_m`h' = J[4 - `h', 7] + } + * parity-gated post rows: att + se + obs + foreach h of numlist 0 1 2 3 4 { + scalar es_`mode'_att_`h' = J[4 + `h', 1] + scalar es_`mode'_se_`h' = J[4 + `h', 2] + scalar es_`mode'_N_`h' = J[4 + `h', 7] + } + scalar es_`mode'_att_ppre = P[1, 1] + scalar es_`mode'_N_ppre = P[1, 7] + scalar es_`mode'_att_ppost = P[2, 1] + scalar es_`mode'_se_ppost = P[2, 2] + scalar es_`mode'_N_ppost = P[2, 7] + } + * in-.do smoke gate: Arm B vw tau0 vs the library-computed point + * (library-anchored, informational). + assert reldif(es_vw_att_0, rpt_es) < 1e-6 +restore + +* ------------------------------------------------------------------------------ +* Arm C: Eq. 13, FULL panel, vw only, event study only. Divergence documentation +* (att + obs; gate 6a asserts both the att divergence and the obs mismatch that +* locks the boundary/respell sample-admission conventions). +* ------------------------------------------------------------------------------ +quietly count +assert r(N) == 840 // guard: Arm B's restore really returned the full panel +lpdid y, unit(unit) time(time) treat(treat) pre_window(3) post_window(4) /// + nonabsorbing(3) nograph only_event +matrix J = e(results) +foreach h of numlist 0 1 2 3 4 { + scalar ef_att_`h' = J[4 + `h', 1] + scalar ef_N_`h' = J[4 + `h', 7] +} + +* ------------------------------------------------------------------------------ +* Emit JSON by hand (Stata has no jsonlite) at %21.17g (round-trip-exact double). +* No timestamp (byte-identical regeneration). +* ------------------------------------------------------------------------------ +* Stata compound-quote note: the `"' delimiter swallows a trailing double-quote, +* so every STRING field is written with a trailing comma (`...","' -> `...",`); +* the meta block therefore ENDS with numeric fields, which have no trailing quote. +local sver = strtrim(string(c(stata_version), "%4.1f")) +* c(flavor) misreports "IC" on SE and c(edition) is unreliable; c(MP)/c(SE) are +* the authoritative 0/1 flags (c(BE) is UNDEFINED, and cond() evaluates all +* branches eagerly, so BE is resolved by elimination, not referenced). +local sedition = cond(c(MP)==1, "MP", cond(c(SE)==1, "SE", "BE")) + +* Helper: write one ES entry with att/se/N. Assembles trimmed locals first. +tempname fh +file open `fh' using "benchmarks/data/lpdid_nonabsorbing_stata_golden.json", write replace text +file write `fh' "{" _n +file write `fh' `" "meta": {"' _n +file write `fh' `" "estimator": "LPDiD non-absorbing SEs - authors' Stata lpdid package, end-to-end","' _n +file write `fh' `" "generator": "benchmarks/stata/generate_lpdid_nonabsorbing_golden.do","' _n +file write `fh' `" "source_panel": "benchmarks/data/lpdid_nonabsorbing_panel.csv","' _n +file write `fh' `" "source_sha256": "71af98c96dec08abc281002b7b88df74033dac7d0ff7707eff9f3ee0c30a3964","' _n +file write `fh' `" "cmd": "lpdid y, unit(unit) time(time) treat(treat) pre_window(3) post_window(4) nograph + nonabsorbing(, firsttreat notyet) [rw] on the full panel; nonabsorbing(3) [rw] on the subsample and (only_event, vw) on the full panel","' _n +file write `fh' `" "mapping": "nonabsorbing(, firsttreat notyet) == first_entry; nonabsorbing(L) == effect_stabilization","' _n +file write `fh' `" "drop_rule": "min(treat)==1 | any(dD==1 & L3.dD==-1)","' _n +file write `fh' `" "point_anchor": "first_entry vw tau0 gated in-.do vs the committed R golden first_entry[0] att (R-anchored); effect_stab_sub vw tau0 vs the library point (library-anchored, informational); the Python test is authoritative for every surface","' _n +file write `fh' `" "se_convention": "package reghdfe vce(cluster unit) CR1; the library's cluster SE matches at ~1e-16 (vw) / ~1e-9 (rw) on identical samples","' _n +file write `fh' `" "convention_notes": "effect_stab full-panel and placebo/pooled-pre rows are DIVERGENCE DOCUMENTATION, not parity: the package admits always-treated units at early t via missing-lag semantics, requires L+1 untreated periods before re-entry (stricter than JAE Eq. 13's levels condition), and builds placebo samples by recursive lagged intersection; see REGISTRY ## LPDiD Deviation 4","' _n +file write `fh' `" "ssc_versions": {"' _n +file write `fh' `" "lpdid": "`v_lpdid'","' _n +file write `fh' `" "reghdfe": "`v_reghdfe'","' _n +file write `fh' `" "ftools": "`v_ftools'","' _n +file write `fh' `" "require": "`v_require'","' _n +file write `fh' `" "boottest": "`v_boottest'","' _n +file write `fh' `" "listreg": "`v_listreg'","' _n +file write `fh' `" "egenmore": "`v_egenmore'","' _n +file write `fh' `" "egenmore_gclsst": "`v_egenmore_gclsst'""' _n +file write `fh' " }," _n +file write `fh' `" "dropped_units": [24, 25, 27, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40],"' _n +file write `fh' `" "stata_edition": "`sedition'","' _n +file write `fh' `" "stata_version": `sver',"' _n +file write `fh' `" "pre_window": 3,"' _n +file write `fh' `" "post_window": 4,"' _n +file write `fh' `" "stabilization_window": 3,"' _n +file write `fh' `" "n_rows_full": 840,"' _n +file write `fh' `" "n_rows_sub": 658"' _n +file write `fh' " }," _n + +* --- first_entry block: vw + rw, full ES + pooled, all fields --- +file write `fh' `" "first_entry": {"' _n +local firstmode 1 +foreach mode in vw rw { + if !`firstmode' file write `fh' "," _n + local firstmode 0 + file write `fh' `" "`mode'": {"' _n + file write `fh' `" "es": {"' + local sep "" + foreach h of numlist 3 2 { + _jnum fe_`mode'_att_m`h' + local a = r(s) + _jnum fe_`mode'_se_m`h' + local s = r(s) + _jnum fe_`mode'_N_m`h' "%12.0f" + local nn = r(s) + file write `fh' "`sep'" _n `" "-`h'": {"att": `a', "se": `s', "N": `nn'}"' + local sep "," + } + foreach h of numlist 0 1 2 3 4 { + _jnum fe_`mode'_att_`h' + local a = r(s) + _jnum fe_`mode'_se_`h' + local s = r(s) + _jnum fe_`mode'_N_`h' "%12.0f" + local nn = r(s) + file write `fh' "`sep'" _n `" "`h'": {"att": `a', "se": `s', "N": `nn'}"' + local sep "," + } + file write `fh' _n " }," _n + file write `fh' `" "pooled": {"' _n + foreach w in pre post { + _jnum fe_`mode'_att_p`w' + local a = r(s) + _jnum fe_`mode'_se_p`w' + local s = r(s) + _jnum fe_`mode'_N_p`w' "%12.0f" + local nn = r(s) + local wsep = cond("`w'" == "pre", ",", "") + file write `fh' `" "`w'": {"att": `a', "se": `s', "N": `nn'}`wsep'"' _n + } + file write `fh' " }" _n + file write `fh' " }" +} +file write `fh' _n " }," _n + +* --- effect_stab_sub block: vw + rw; post rows full fields, pre rows att+N --- +file write `fh' `" "effect_stab_sub": {"' _n +file write `fh' `" "n_rows": 658,"' _n +local firstmode 1 +foreach mode in vw rw { + if !`firstmode' file write `fh' "," _n + local firstmode 0 + file write `fh' `" "`mode'": {"' _n + file write `fh' `" "es": {"' + local sep "" + foreach h of numlist 3 2 { + _jnum es_`mode'_att_m`h' + local a = r(s) + _jnum es_`mode'_N_m`h' "%12.0f" + local nn = r(s) + file write `fh' "`sep'" _n `" "-`h'": {"att": `a', "N": `nn'}"' + local sep "," + } + foreach h of numlist 0 1 2 3 4 { + _jnum es_`mode'_att_`h' + local a = r(s) + _jnum es_`mode'_se_`h' + local s = r(s) + _jnum es_`mode'_N_`h' "%12.0f" + local nn = r(s) + file write `fh' "`sep'" _n `" "`h'": {"att": `a', "se": `s', "N": `nn'}"' + local sep "," + } + file write `fh' _n " }," _n + file write `fh' `" "pooled": {"' _n + _jnum es_`mode'_att_ppre + local a = r(s) + _jnum es_`mode'_N_ppre "%12.0f" + local nn = r(s) + file write `fh' `" "pre": {"att": `a', "N": `nn'},"' _n + _jnum es_`mode'_att_ppost + local a = r(s) + _jnum es_`mode'_se_ppost + local s = r(s) + _jnum es_`mode'_N_ppost "%12.0f" + local nn = r(s) + file write `fh' `" "post": {"att": `a', "se": `s', "N": `nn'}"' _n + file write `fh' " }" _n + file write `fh' " }" +} +file write `fh' _n " }," _n + +* --- effect_stab_full_vw block: att + N only --- +file write `fh' `" "effect_stab_full_vw": {"' _n +file write `fh' `" "es": {"' +local sep "" +foreach h of numlist 0 1 2 3 4 { + _jnum ef_att_`h' + local a = r(s) + _jnum ef_N_`h' "%12.0f" + local nn = r(s) + file write `fh' "`sep'" _n `" "`h'": {"att": `a', "N": `nn'}"' + local sep "," +} +file write `fh' _n " }" _n +file write `fh' " }" _n +file write `fh' "}" _n +file close `fh' + +display "Wrote benchmarks/data/lpdid_nonabsorbing_stata_golden.json (arms A/B/C)" diff --git a/benchmarks/stata/requirements.do b/benchmarks/stata/requirements.do index 3e436c1e..540e0c97 100644 --- a/benchmarks/stata/requirements.do +++ b/benchmarks/stata/requirements.do @@ -25,11 +25,34 @@ *! (both consumed by generate_etwfe_cs_golden.do) *! lwdid - Lee & Wooldridge rolling DiD, the authors' reference *! implementation (consumed by generate_lwdid_golden.do) +*! lpdid - Dube-Girardi-Jorda-Taylor LP-DiD, the authors' reference +*! implementation (consumed by +*! generate_lpdid_nonabsorbing_golden.do) +*! boottest, egenmore, listreg +*! - lpdid startup dependencies (lpdid `which`-checks all three +*! and exits without them; egenmore's filter() runs in every +*! lpdid pooled spec). egenmore installs no `egenmore.ado` - +*! only `_g*.ado` helpers - so the loop probes BOTH +*! `_gfilter` (executed by the pooled spec) and `_gclsst` +*! (lpdid's own startup which-check) and reinstalls the +*! `egenmore` package if either is missing; requiring both +*! means an incomplete egenmore install is repaired by +*! rerunning this script. version 19 -foreach p in ftools require reghdfe did_imputation drdid csdid hdfe jwdid lwdid { - capture which `p' - if _rc { +foreach p in ftools require reghdfe did_imputation drdid csdid hdfe jwdid lwdid boottest egenmore listreg lpdid { + local missing 0 + if "`p'" == "egenmore" { + foreach probe in _gfilter _gclsst { + capture which `probe' + if _rc local missing 1 + } + } + else { + capture which `p' + if _rc local missing 1 + } + if `missing' { di as txt "Installing `p' from SSC ..." ssc install `p', replace } diff --git a/diff_diff/guides/llms-full.txt b/diff_diff/guides/llms-full.txt index 4b13d380..416a8a20 100644 --- a/diff_diff/guides/llms-full.txt +++ b/diff_diff/guides/llms-full.txt @@ -1196,7 +1196,7 @@ print(grp.to_dataframe()) ### 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 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. Complex-survey designs are supported on the variance-weighted default path via the `survey_design=` argument to `fit()` (probability weights enter the WLS point estimate; the SE is the stratified-PSU Taylor-linearization sandwich with `df = n_PSU - n_strata`, with optional FPC and lonely-PSU handling) — rejected with `reweight=True`, replicate weights, or non-pweight types. +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)` on the Stata-anchored surfaces - non-absorbing parity is scoped because the two implementations' clean-sample conventions differ at panel boundaries and on placebo windows; see the methodology registry, LPDiD Deviation 4); 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. Complex-survey designs are supported on the variance-weighted default path via the `survey_design=` argument to `fit()` (probability weights enter the WLS point estimate; the SE is the stratified-PSU Taylor-linearization sandwich with `df = n_PSU - n_strata`, with optional FPC and lonely-PSU handling) — rejected with `reweight=True`, replicate weights, or non-pweight types. ```python LPDiD( diff --git a/docs/api/lpdid.rst b/docs/api/lpdid.rst index 48fbf725..ec0951c6 100644 --- a/docs/api/lpdid.rst +++ b/docs/api/lpdid.rst @@ -26,9 +26,18 @@ estimand is a strictly non-negatively-weighted average of cohort effects. never-treated units). Non-absorbing modes require a gap-free panel within each unit's observed span and cover the entry-effect estimands. The non-absorbing entry-effect paths are R-parity-validated against an independent - ``fixest::feols`` reconstruction of the paper's Eq. 12/13 (see - ``docs/methodology/REGISTRY.md``); the Appendix-C exit-event dynamics and the - Stata canonical SE remain planned follow-ups. + ``fixest::feols`` reconstruction of the paper's Eq. 12/13 and anchored against + the authors' Stata ``lpdid`` package end-to-end — ``first_entry`` (Eq. 12) on + all surfaces including placebos and both pooled windows; + ``effect_stabilization`` (Eq. 13) at post horizons + pooled post on a + convention-neutral subsample. Outside that anchor sit the Eq. 13 + placebo/pooled-pre inference (the package builds placebo samples with a + different window construction) and full-panel Eq. 13 inference (the two + implementations' boundary/re-entry sample conventions differ there; the + subsample is what makes the samples coincide) - see + ``docs/methodology/REGISTRY.md``, LPDiD Deviation 4. + The Appendix-C exit-event dynamics and the non-absorbing RA path remain + planned follow-ups. Complex-survey designs (probability weights + stratified-PSU Taylor-linearization standard errors with optional finite-population correction and lonely-PSU handling) are supported on the variance-weighted diff --git a/docs/benchmarks.rst b/docs/benchmarks.rst index f1ed32a0..8d10aa7b 100644 --- a/docs/benchmarks.rst +++ b/docs/benchmarks.rst @@ -1,5 +1,5 @@ .. meta:: - :description: Validation benchmarks comparing diff-diff against R packages (did, synthdid, fixest) and Stata (teffects, did_imputation, jwdid/csdid, reghdfe, lwdid). Coefficient accuracy, standard error comparison, and performance metrics. + :description: Validation benchmarks comparing diff-diff against R packages (did, synthdid, fixest) and Stata (teffects, lpdid, did_imputation, jwdid/csdid, reghdfe, lwdid). Coefficient accuracy, standard error comparison, and performance metrics. :keywords: difference-in-differences benchmark, DiD validation R, DiD validation Stata, python econometrics accuracy, did package comparison Benchmarks @@ -46,6 +46,13 @@ where no runnable R reference exists: * - ``LPDiD`` (regression-adjustment SE) - Stata ``teffects ra ... atet`` - Dube, Girardi, Jorda & Taylor (2025); no runnable R analogue + * - ``LPDiD`` (non-absorbing SEs) + - Stata ``lpdid, nonabsorbing(...)`` (the authors' package, end-to-end) + - Dube, Girardi, Jorda & Taylor (2025) §4.2; first external anchor for the + non-absorbing reweighted SE and pooled windows. Eq. 12 on all surfaces; + Eq. 13 at post horizons + pooled post on a convention-neutral subsample + (three package convention differences measured and divergence-gated — see + the methodology registry, LPDiD Deviation 4) * - ``ImputationDiD`` (leave-one-out SE) - Stata ``did_imputation, leaveout`` - Borusyak, Jaravel & Spiess (2024) Supp. App. A.9; no runnable R analogue. @@ -1025,13 +1032,14 @@ Prerequisites pip install -e ".[dev]" 4. (Optional) Stata, only to regenerate the committed Stata goldens (``LPDiD`` - regression-adjustment SE, ``ImputationDiD`` leave-one-out SE, the ETWFE/CS - cross-check, the ``reghdfe`` K_reference convention, and the LWDiD - authors'-package parity). The goldens are committed, so this is not needed to - run the test suite. The ``LPDiD`` arm uses the **native** ``teffects`` command; - the other arms depend on SSC packages + regression-adjustment SE, ``LPDiD`` non-absorbing SEs, ``ImputationDiD`` + leave-one-out SE, the ETWFE/CS cross-check, the ``reghdfe`` K_reference + convention, and the LWDiD authors'-package parity). The goldens are committed, + so this is not needed to run the test suite. The ``LPDiD`` RA arm uses the + **native** ``teffects`` command; the other arms depend on SSC packages (``did_imputation``/``reghdfe``/``ftools``/``require``, - ``drdid``/``csdid``/``jwdid``/``hdfe``, ``lwdid``) — install them once via + ``drdid``/``csdid``/``jwdid``/``hdfe``, ``lwdid``, + ``lpdid``/``boottest``/``egenmore``/``listreg``) — install them once via ``benchmarks/stata/requirements.do`` (the generators do not auto-install): .. code-block:: bash @@ -1040,6 +1048,7 @@ Prerequisites STATA=/Applications/Stata/StataSE.app/Contents/MacOS/stata-se $STATA -b do benchmarks/stata/requirements.do # one-time SSC install $STATA -b do benchmarks/stata/generate_lpdid_ra_golden.do + $STATA -b do benchmarks/stata/generate_lpdid_nonabsorbing_golden.do $STATA -b do benchmarks/stata/generate_imputation_loo_golden.do $STATA -b do benchmarks/stata/generate_etwfe_cs_golden.do $STATA -b do benchmarks/stata/generate_reghdfe_kref_golden.do diff --git a/docs/doc-deps.yaml b/docs/doc-deps.yaml index 3f4f2354..7181b91f 100644 --- a/docs/doc-deps.yaml +++ b/docs/doc-deps.yaml @@ -802,8 +802,12 @@ sources: - path: docs/methodology/REGISTRY.md section: "LPDiD" type: methodology + - path: docs/methodology/papers/dube-2025-review.md + type: methodology - path: benchmarks/stata/README.md type: methodology + - path: docs/benchmarks.rst + type: performance - path: docs/api/lpdid.rst type: api_reference - path: README.md diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index 2d696a9c..b38dfa1f 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -2468,7 +2468,7 @@ Eq. 12 reuses the absorbing clean control and only restricts the treated set (a ### 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). +**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). The paper's materials mention stacking per-horizon regressions (`suest`) for JOINT tests, but the shipped package (v1.0.3, measured by the non-absorbing arm) does NOT obtain pooled inference that way: it constructs a pooled long-difference outcome and fits the pooled windows as separate unit-clustered `reghdfe` regressions - the same construction as Note 6's fixed-composition pooled estimand, which is why the library's pooled SEs match the package to ~1e-16/1e-9 on the anchored surfaces. 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). The absorbing RA SE is measured against Stata `teffects` (Deviation 2); the non-absorbing cluster SEs are measured against the authors' Stata `lpdid` package end-to-end (Eq. 12 full-panel at SE ~6e-17 vw / ~2.3e-9 rw incl. both pooled windows; Eq. 13 post + pooled-post on a convention-neutral subsample at ~1.8e-16 vw / ~1.6e-9 rw - Deviation 4). ### Edge Cases @@ -2476,16 +2476,16 @@ Eq. 12 reuses the absorbing clean control and only restricts the treated set (a - **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):** 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` (PR-C2 documented this as a known divergence from `alexCardazzi/lpdid`, which NA-excludes such first-rows - see Deviation #4). **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). +- **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` (PR-C2 documented this as a known divergence from `alexCardazzi/lpdid`, which NA-excludes such first-rows; the authors' Stata `lpdid` package differs here too, its missing-lag semantics admitting always-treated units as controls at early `t` - the paper does not specify boundary handling, so all three are documented conventions; see Deviation #4). **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 The paper specifies no standard-error formula (Section 1 defers to "standard, well-understood techniques"); the reference Stata `lpdid` uses `vce(cluster unit)`. The entries below document diff-diff's inference and scope choices. -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. **Since 3.9 (M-127) the reference is the `df_convention` knob's `"cluster"` value — LPDiD's DEFAULT, so nothing moved**: `"residual"` opts into t(per-horizon residual df `n_eff − k_kept`; on the RA path `n_total − k0_kept − 1`, the pooled M-estimator's parameter count — nuisance coefficients plus the ATT — a library convention with no external anchor since Stata `teffects ra` reports z), and `"normal"` into z. The G here stays the RAW unique cluster count (deliberately not `effective_cluster_count`: LPDiD's reweights are strictly positive and there is no user weight column, so the two cannot diverge). **Carve-out to the flat `t(G-1)` claim:** the degenerate lanes — the unclustered-refit fallback when the clustered solve fails, the RA `G<=1` lane, and the saturated-design early return (`n_obs <= k`) — keep the literal `df=None` (normal theory, NaN provenance) under ALL knob values, exactly as before 3.9. +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)` - measured on the Stata-anchored surfaces (absorbing via the R/`teffects` arms; non-absorbing per Deviation 4's scoping: Eq. 12 everywhere, Eq. 13 post + pooled-post on the convention-neutral subsample; elsewhere the two implementations' clean-sample conventions differ, see Deviation 4); the paper prescribes no SE. **Since 3.9 (M-127) the reference is the `df_convention` knob's `"cluster"` value — LPDiD's DEFAULT, so nothing moved**: `"residual"` opts into t(per-horizon residual df `n_eff − k_kept`; on the RA path `n_total − k0_kept − 1`, the pooled M-estimator's parameter count — nuisance coefficients plus the ATT — a library convention with no external anchor since Stata `teffects ra` reports z), and `"normal"` into z. The G here stays the RAW unique cluster count (deliberately not `effective_cluster_count`: LPDiD's reweights are strictly positive and there is no user weight column, so the two cannot diverge). **Carve-out to the flat `t(G-1)` claim:** the degenerate lanes — the unclustered-refit fallback when the clustered solve fails, the RA `G<=1` lane, and the saturated-design early return (`n_obs <= k`) — keep the literal `df=None` (normal theory, NaN provenance) under ALL knob values, exactly as before 3.9. 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_reference))` finite-sample factor (see the nested-dummy Note below) - 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)`, while the default path matches `feols`/`reghdfe`. **Originally inferred from the authors' `lpdid_regression_adjustment.do` degrees-of-freedom comments, this is now MEASURED against a runnable `teffects`:** `benchmarks/stata/generate_lpdid_ra_golden.do` runs `teffects ra (Dy x i.time) (tdiff), atet vce(cluster unit)` on an independently reconstructed clean sample (a Stata port of the `generate_lpdid_golden.R` `prep`/`clean_h` recipe), and the library RA IF SE matches it to **~1e-16 at all 7 event-study horizons** - the SE ratio is a flat 1.0 (no `G/(G-1)`, no `(n-1)/(n-k)`), directly confirming the no-finite-sample-factor convention (golden `benchmarks/data/lpdid_ra_stata_golden.json`; test `tests/test_lpdid_ra_stata_parity.py`, which also gates the Stata clean sample's realized size and cluster count `(e(N), e(N_clust))` against the library's `(n_obs, n_clusters)`). The RA *point* estimate is triple-anchored: Python == R full-interaction (`i.dtreat##(i.time c.x)`) == Stata `teffects`, to ~1e-10..1e-13 (`tests/test_methodology_lpdid.py::test_ra_covariate_point`). Because Stata is node-locked and cannot run in CI (goldens are committed, exactly like the R arm), the CI-runnable guards remain the regression pin (`test_ra_covariate_se_regression_pin` / `RA_SE_PIN`) and 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 - non-absorbing treatment (Section 4.2) implements the **entry-effect** estimands (`non_absorbing="first_entry"` / `"effect_stabilization"`, PR-C1). **PR-C2 R-parity-validated both modes against an INDEPENDENT `fixest::feols` reconstruction of the paper's Eq. 12 / Eq. 13 clean-sample restrictions** (point and SE match to ~1e-13/~1e-15 for the variance-weighted variants; the `effect_stabilization` reweighted point matches and its SE is pinned as a regression guard - a small weighted-cluster convention difference vs feols; `tests/test_methodology_lpdid.py::TestLPDiDNonAbsorbingParityR`). The recipe's independence was demonstrated when an earlier draft's Eq. 12 control off-by-one diverged from the already-correct library and was corrected against the paper, plus a hand-computed Python micro-check. **`alexCardazzi/lpdid`'s `nonabsorbing_lag` is NOT a faithful Eq. 13** (it clamps `treat_diff[<0]<-0`, so its clean-control window blocks only treatment turn-*ons*; it reuses a forward placebo window; and it NA-excludes pre-panel-treated rows where the library clamps pre-`min_t` to untreated): it diverges ~0.01-0.05 from Eq. 13 even on a monotone no-off-switch panel, so it is **recorded in the golden `meta` as a divergent third-party reference, not a parity gate** (the alexCardazzi-pooled precedent). The library's "no treatment change" (both directions) and backward placebo window are the more paper-faithful choices. `first_entry` (Eq. 12) has no R-package analogue (anchored on the independent feols recipe only). Appendix-C exit-event dynamics and the Stata canonical SE remain deferred follow-ups. +4. **Deviation from R:** Scope - non-absorbing treatment (Section 4.2) implements the **entry-effect** estimands (`non_absorbing="first_entry"` / `"effect_stabilization"`, PR-C1). **PR-C2 R-parity-validated both modes against an INDEPENDENT `fixest::feols` reconstruction of the paper's Eq. 12 / Eq. 13 clean-sample restrictions** (point and SE match to ~1e-13/~1e-15 for the variance-weighted variants; the `effect_stabilization` reweighted point matches and its SE is pinned as a regression guard - a small weighted-cluster convention difference vs feols; `tests/test_methodology_lpdid.py::TestLPDiDNonAbsorbingParityR`). The recipe's independence was demonstrated when an earlier draft's Eq. 12 control off-by-one diverged from the already-correct library and was corrected against the paper, plus a hand-computed Python micro-check. **`alexCardazzi/lpdid`'s `nonabsorbing_lag` is NOT a faithful Eq. 13** (it clamps `treat_diff[<0]<-0`, so its clean-control window blocks only treatment turn-*ons* - a violation of Eq. 13's stated `Delta-D = 0` both-directions control condition; it also holds two convention differences on surfaces Eq. 13 does not specify: a forward placebo window, and NA-exclusion of pre-panel-treated rows where the library clamps pre-`min_t` to untreated): it diverges ~0.01-0.05 from Eq. 13 even on a monotone no-off-switch panel, so it is **recorded in the golden `meta` as a divergent third-party reference, not a parity gate** (the alexCardazzi-pooled precedent; the golden's `meta.alex_note`/`meta.boundary_note` strings predate this two-tier framing and are generation-time text - this entry is the living statement). `first_entry` (Eq. 12) has no R-package analogue (anchored on the independent feols recipe, and - since the Stata non-absorbing arm - on the authors' Stata `lpdid` package end-to-end via `nonabsorbing(, firsttreat notyet)`). **The Stata canonical SE is now MEASURED against the authors' `lpdid` package (SSC, v1.0.3) end-to-end** - the first external anchor for the non-absorbing reweighted SE, the non-absorbing pooled windows (points and SEs), and the Eq. 12 reweighted point: Eq. 12 agrees on the full committed panel at every surface (all event-study horizons incl. placebos + both pooled windows; att ~1e-14 vw / ~1.3e-8 rw, SE ~6e-17 vw / ~2.3e-9 rw), and Eq. 13 agrees at post horizons + pooled post on a convention-neutral subsample (att ~1.1e-14 vw / ~1.9e-8 rw, SE ~1.8e-16 vw / ~1.6e-9 rw; golden `benchmarks/data/lpdid_nonabsorbing_stata_golden.json`, test `tests/test_lpdid_nonabsorbing_stata_parity.py`). **Three measured convention differences vs the package**, adjudicated against the paper's stated Eq. 13 restriction (verified against the open-access working-paper text, whose Eq. 18 is the published JAE Eq. 13): (a) *exact-L re-entry spells* - Eq. 13's treated condition is a LEVELS condition (`D = 0` on `[t-L, t-1]`), which the library implements; the package's switch-free lag window additionally excludes an exit AT `t-L`, effectively requiring `L+1` untreated periods - **the package is stricter than the paper here**; (b) *pre-panel boundary handling* - the package's missing-lag semantics admit always-treated units as controls at early `t` where the library's clamp excludes them (a paper-silent surface; both are documented conventions); (c) *placebo windows* - the package builds `h<0` clean samples by recursive lagged intersection of its `CCS_0` indicator, the library uses the backward window `[t-max(L,-h), t-1]` (also paper-silent; both are conventions). The full-panel Eq. 13 divergence from (a)+(b) and the placebo/pooled-pre divergence from (c) are recorded in the golden and locked by divergence-documentation gates (att floors + sample-size mismatch asserts), so parity claims stay scoped: Eq. 13 full-panel and placebo/pooled-pre inference remain OUTSIDE the Stata anchor (the full-panel reweighted SEs stay pinned via `RW_SE_PIN`; the pooled-pre SE has no pin or anchor). Appendix-C exit-event dynamics and the non-absorbing RA path (the package's RA engine is `listreg`, and the committed non-absorbing panel has no covariate) remain deferred follow-ups. 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`). @@ -2507,7 +2507,8 @@ The paper specifies no standard-error formula (Section 1 defers to "standard, we - [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) - [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) - [x] Non-absorbing R-parity: both modes vs an independent `fixest::feols` Eq. 12/13 reconstruction (point+SE ~1e-13/~1e-15 vw; reweighted point + pinned SE); `alexCardazzi nonabsorbing_lag` recorded as a divergent reference (not a gate); absorbing B2 goldens byte-identical (PR-C2) -- [ ] Non-absorbing exit-event dynamics (Appendix C `eta_h`) + the Stata canonical RA/SE - deferred +- [x] Non-absorbing Stata canonical SE: authors' `lpdid` package end-to-end - Eq. 12 all surfaces (ES incl. placebos + both pooled windows, vw + rw) and Eq. 13 post + pooled-post on a convention-neutral subsample; three package convention differences measured, adjudicated vs the paper, and divergence-gated (Deviation 4) +- [ ] Non-absorbing exit-event dynamics (Appendix C `eta_h`) + the non-absorbing RA path (package RA engine is `listreg`; needs a covariate panel) + Eq. 13 placebo/pooled-pre inference (not anchored by the current arm - the placebo-window convention differs, and pooled-pre currently has no pin or anchor; a dedicated late-entry/never-treated subsample could in principle align the placebo constructions) - deferred - [x] Survey-design support (PR-D1): pweight + stratified-PSU Taylor-linearization (Binder TSL) variance on the variance-weighted default path; per-sample design re-resolution + unit-as-PSU injection; reweight (incl. RA), replicate-weight, and non-pweight designs rejected; pure-Python invariants (reduction/unit-clustering, FPC-shrinks-SE, stratification, lonely-PSU, NaN-consistency, metadata) (PR-D1) - [x] Survey-design R-parity (PR-D2): all three survey paths (VW strata+PSU+FPC, weights-only inject, direct-covariate) validated end-to-end vs `survey::svyglm` - per-horizon point/SE/df + pooled (point ~1e-6, SE ~1e-5, df exact); clean sample independently `alexCardazzi`-cross-checked (<1e-8); dedicated survey panel keeps the absorbing / non-absorbing goldens byte-identical (PR-D2) diff --git a/docs/methodology/papers/dube-2025-review.md b/docs/methodology/papers/dube-2025-review.md index 1b035c90..4e5a16e4 100644 --- a/docs/methodology/papers/dube-2025-review.md +++ b/docs/methodology/papers/dube-2025-review.md @@ -97,7 +97,7 @@ Single overall ATT over the posttreatment window `h in {0..H}` by using the post *Standard errors (Section 1 - NOT specified by the paper):* - **The paper deliberately gives no SE formula.** Section 1: "all the estimators ... allow for standard statistical inference using well-understood techniques ... For this reason, we do not discuss statistical inference here." - Default in the reference Stata implementation: **cluster-robust at the unit level** (`vce(cluster unit)`, footnote 9). -- Pooled / joint tests across horizons: stack the per-horizon regressions (Stata `suest`). +- Joint tests across horizons: stack the per-horizon regressions (Stata `suest`). (Pooled pre/post inference is a different object - the shipped package v1.0.3 computes it by a separate unit-clustered regression on a pooled long-difference outcome, NOT by stacking; measured by this repo's non-absorbing Stata arm, see REGISTRY `## LPDiD` Standard Errors + Note 6.) - Bootstrap: not discussed in the paper. - => Any analytical SE we ship (and especially the influence-function variance for the RA path) is an **implementation choice to be validated against the reference package**, not against the paper. Document under "Deviations from the paper". @@ -120,6 +120,7 @@ Single overall ATT over the posttreatment window `h in {0..H}` by using the post - R: `alexCardazzi/lpdid` (third-party, A. Cardazzi & Z. Porreca; covers absorbing AND non-absorbing); authors' own R example scripts in `danielegirardi/lpdid`. - Stata RA syntax (footnote 9): `teffects ra (Dhy i.time) (dtreat) if D.treat==1 | Fh.treat==0, atet vce(cluster unit)`. - **RA SE parity (this repo):** the library RA influence-function SE is anchored against this footnote-9 `teffects` recipe (with the covariate added to the outcome model) by `benchmarks/stata/generate_lpdid_ra_golden.do` (golden `benchmarks/data/lpdid_ra_stata_golden.json`, test `tests/test_lpdid_ra_stata_parity.py`) — an independent Stata reconstruction of each horizon's clean sample, matching to ~1e-16 at all 7 event-study horizons and confirming the no-finite-sample-factor convention (REGISTRY LPDiD Deviation 2). +- **Non-absorbing SE parity (this repo):** the Section-4.2 entry-effect modes are anchored against the authors' SSC `lpdid` package end-to-end by `benchmarks/stata/generate_lpdid_nonabsorbing_golden.do` (golden `benchmarks/data/lpdid_nonabsorbing_stata_golden.json`, test `tests/test_lpdid_nonabsorbing_stata_parity.py`) — Eq. 12 on all surfaces of the committed panel, Eq. 13 at post horizons + pooled post on a convention-neutral subsample; three package convention differences measured and documented (REGISTRY LPDiD Deviation 4). **Requirements checklist:** - [ ] Long-difference dependent variable with selectable base period (`t-1` default, PMD `k`) @@ -200,4 +201,4 @@ The paper *proves* numerical equivalences we can exploit for **internal cross-va 3. **Non-absorbing extension (Section 4.2)** is presented as illustrative, not exhaustive ("a comprehensive discussion ... would indeed require a whole article", p.751). The effect-stabilization window `L` (Assumption 9, Equation 13) and the "first-time entry" estimand (Equation 12) are two distinct estimands; a Phase C design will need to choose/expose both deliberately. The contributor's scaffold rejects non-absorbing entirely. 4. **PMD multi-cohort vs BJS** are "very similar (although not identical)" (p.747) - do not assert exact equality in tests except in the single-cohort `k=t-1` case. 5. **`pmd="max"` semantics:** the paper's `k=t-1` "use all pretreatment periods" is per-observation (expanding window). The contributor implemented `pmd="max"` as the expanding mean of all prior periods and integer `pmd=k` as the trailing-`k` mean - verify this matches the package's `pmd` option exactly (the paper fixes the formula but not the option's edge behavior at the panel start). -6. **Pooled SE / joint tests** use stacking (`suest`); the variance of the pooled estimand under clustering should be cross-checked against the package rather than assembled ad hoc. +6. **Joint tests** use stacking (`suest`); pooled inference is separate. ~~The variance of the pooled estimand under clustering should be cross-checked against the package rather than assembled ad hoc.~~ **PARTIALLY RESOLVED** - the shipped package computes pooled pre/post inference as a separate unit-clustered regression on a pooled long-difference outcome, and this repo's non-absorbing Stata arm cross-checked the library's pooled SEs against it end-to-end (Eq. 12 both pooled windows; Eq. 13 pooled post; see REGISTRY `## LPDiD` Standard Errors). Eq. 13 pooled-pre remains outside that anchor (placebo-window convention difference). diff --git a/tests/test_lpdid_nonabsorbing_stata_parity.py b/tests/test_lpdid_nonabsorbing_stata_parity.py new file mode 100644 index 00000000..007cc4b2 --- /dev/null +++ b/tests/test_lpdid_nonabsorbing_stata_parity.py @@ -0,0 +1,502 @@ +"""External-reference parity: LPDiD non-absorbing SEs vs the authors' Stata ``lpdid``. + +The non-absorbing LP-DiD modes (Dube, Girardi, Jorda & Taylor 2025, JAE Eq. 12 / +Eq. 13) were R-parity-locked only on the event-study surfaces the R arm covers +(independent ``fixest::feols`` reconstruction; no Eq. 12 reweighted arm, no pooled +block). The canonical reference is the authors' Stata ``lpdid`` package (SSC), which +this arm runs END-TO-END - the package builds its own clean samples. This anchors, +for the first time externally: the non-absorbing REWEIGHTED SE (previously +pinned-only via ``RW_SE_PIN``), the non-absorbing POOLED windows (points and SEs), +and the Eq. 12 reweighted point. + +Mapping (pinned in the golden ``meta`` and gated below):: + + lpdid, nonabsorbing(, firsttreat notyet) == LPDiD(non_absorbing="first_entry") + lpdid, nonabsorbing(L) == LPDiD(non_absorbing="effect_stabilization", + stabilization_window=L) + +Parity scoping. Eq. 12 (pseudo-absorbing) agrees on EVERY surface of the full +committed panel - all event-study horizons including placebos, plus both pooled +windows, vw and rw. Eq. 13 agrees on POST horizons + pooled post only, and only on a +convention-neutral SUBSAMPLE, because the package's sample construction differs from +the library's on three measured conventions (REGISTRY ``## LPDiD`` Deviation 4): + +1. *Boundary*: the package's missing-lag semantics admit always-treated units as + controls at early t; the library's pre-panel-untreated clamp excludes them until + t >= L+2 (a paper-silent surface - both are conventions). +2. *Exact-L respell*: the package's switch-free lag window effectively requires L+1 + untreated periods before re-entry; the library's levels reading of Eq. 13 + (untreated over [t-L, t-1]) admits re-entry after an exactly-L untreated spell + (the package is stricter than the paper's stated levels condition). +3. *Placebo windows*: the package builds pre-horizon clean samples by recursive + lagged intersection of its CCS_0 indicator; the library uses the backward window + [t-max(L,-h), t-1] (paper-silent; both are conventions). + +The subsample drops the units classes 1+2 admit differently (always-treated +{31..40}; exact-L respell {24, 25, 27}), making the Eq. 13 post-horizon and +pooled-post samples row-identical across implementations. Class 3 is NOT +neutralized by this subsample (a dedicated late-entry/never-treated subsample +could in principle align the two placebo constructions - a possible follow-up, +tracked in DEFERRED.md), so Eq. 13 placebo/pooled-pre rows are recorded as +measured DIVERGENCE documentation (att + obs only) and gated as such (gate 6), never +as parity. Note the per-class attribution is established by unit-set enumeration +(gate 5 asserts each sub-class list separately); gate 6(a) documents the JOINT +classes-1+2 divergence - no separately-measured per-class numeric effect is claimed. + +Gates: + +1. **Stata Eq. 12 vw vs the committed R golden** ``first_entry`` (committed vs + committed, ``abs=1e-9``): the package's end-to-end run reproduces the independent + feols reconstruction. +2. **Library ``first_entry`` vs Stata**, vw + rw: all ES horizons INCLUDING placebos + + pooled pre/post (att ``abs=1e-6``, SE ``abs=1e-7``). +3. **Library ``effect_stabilization`` subsample vs Stata**, vw + rw: post horizons + + pooled post (same tolerances). +4. **Sample shape**: Stata per-row obs == library ``n_obs``, exact, on every parity + surface of gates 2-3. The package's single multi-horizon call returns no + per-horizon cluster count (``e(results)`` columns are + coefficient/se/t/p/ci_low/ci_high/obs), so - unlike the ``teffects`` arm - there + is no per-horizon ``e(N_clust)`` to gate alongside obs; the obs equality plus the + 1e-14-class att/SE agreement carry the same-sample corroboration. +5. **Drop-rule consistency**: golden ``dropped_units``/``n_rows`` vs the pandas + recompute (each convention sub-class asserted separately), and the free-text + ``drop_rule``/``mapping`` meta strings vs hard-coded literals. +6. **Convention-divergence documentation**: every recorded divergence surface stays + measurably divergent - (a) Eq. 13 full-panel att + obs vs the library, (b) + Eq. 13 subsample placebo/pooled-pre att + obs vs the library, both at a 1e-3 att + floor. Threshold rationale: 1e-3 sits 2.35x below the smallest measured + divergence (2.35e-3 at full-panel h=0; every other gated surface clears it by + >3.9x). Goldens regenerate only in deliberate local runs (Stata is not in CI); + the regeneration protocol is to re-measure and refresh these floors if an SSC + release narrows a convention gap. +7. **Warning cleanliness** on every library fit. +8. **Provenance**: ``meta.ssc_versions`` covers every fail-closed generator + dependency with a real identifier that always embeds the ado file's checksum + + length (SSC has no immutable archive, so a same-version-string upstream edit + must still move the metadata); ``source_panel``/``source_sha256`` match the + committed panel; prose meta fields are non-empty. + +Guard per ``feedback_golden_file_pytest_skip``: CI isolated-install jobs copy +``tests/`` only, not ``benchmarks/data/``, so a missing fixture downgrades to +pytest.skip. The guards are SPLIT: the Stata golden + panel carry the module, +while the R golden gates only gate 1 - a missing R fixture cannot silence the +primary Stata-vs-library gates. A PRESENT golden missing the expected top-level +blocks HARD-FAILS +instead (a stale committed golden is a real inconsistency, not an absent fixture). +Regenerate with:: + + /Applications/Stata/StataSE.app/Contents/MacOS/stata-se -b do \ + benchmarks/stata/generate_lpdid_nonabsorbing_golden.do +""" + +from __future__ import annotations + +import hashlib +import json +import re +import warnings +from pathlib import Path + +import pandas as pd +import pytest + +from diff_diff import LPDiD + +_DATA = Path(__file__).parent.parent / "benchmarks" / "data" +STATA_GOLDEN_PATH = _DATA / "lpdid_nonabsorbing_stata_golden.json" +R_GOLDEN_PATH = _DATA / "lpdid_nonabsorbing_golden.json" +PANEL_PATH = _DATA / "lpdid_nonabsorbing_panel.csv" + +# Split availability guards: the Stata golden + panel carry the whole module; the +# R golden is needed only by gate 1 (a missing R fixture must not silence the +# primary Stata-vs-library gates). +_STATA_FIXTURE_AVAILABLE = STATA_GOLDEN_PATH.is_file() and PANEL_PATH.is_file() +_R_FIXTURE_AVAILABLE = R_GOLDEN_PATH.is_file() + +# Event-study horizons pinned by the golden; h = -1 is the omitted reference. +ES_HORIZONS = [-3, -2, 0, 1, 2, 3, 4] +POST_HORIZONS = [0, 1, 2, 3, 4] +PRE_HORIZONS = [-3, -2] +MODES = ["vw", "rw"] +L = 3 + +# committed Stata vs committed R golden (R digits=12); measured ~1e-13 locally. +R_ANCHOR_ATOL = 1e-9 +ATT_ATOL = 1e-6 # library (recomputed) vs committed Stata, cross-platform +SE_ATOL = 1e-7 # the repo's LPDiD cross-platform SE tolerance +DIVERGENCE_FLOOR = 1e-3 # gate 6; see module docstring for the rationale + +# Hard-coded twins of the golden's meta strings (gate 5 / gate 8 desync guards). +DROP_RULE = "min(treat)==1 | any(dD==1 & L3.dD==-1)" +MAPPING = ( + "nonabsorbing(, firsttreat notyet) == first_entry; nonabsorbing(L) == effect_stabilization" +) +ALWAYS_TREATED_UNITS = list(range(31, 41)) # convention class 1 +EXACT_L_RESPELL_UNITS = [24, 25, 27] # convention class 2 +SUBSAMPLE_N_ROWS = 658 +# Every fail-closed generator dependency (guard set == version set; both guarded +# egenmore files carry their own drift signal - _gfilter under "egenmore", +# _gclsst under "egenmore_gclsst"). +SSC_VERSION_KEYS = { + "lpdid", + "reghdfe", + "ftools", + "require", + "boottest", + "listreg", + "egenmore", + "egenmore_gclsst", +} + + +def _skip_if_missing() -> None: + if not _STATA_FIXTURE_AVAILABLE: + pytest.skip( + "Stata LPDiD non-absorbing parity fixture not present. Regenerate via " + "`stata-se -b do benchmarks/stata/generate_lpdid_nonabsorbing_golden.do`." + ) + + +@pytest.fixture(scope="module") +def stata_golden() -> dict: + _skip_if_missing() + with STATA_GOLDEN_PATH.open() as f: + golden = json.load(f) + # Present-but-stale golden = a real inconsistency (forgotten regeneration); + # hard-fail rather than silently skipping the external anchor. + missing = {"meta", "first_entry", "effect_stab_sub", "effect_stab_full_vw"} - set(golden) + assert not missing, ( + f"{STATA_GOLDEN_PATH.name} is present but missing block(s) {sorted(missing)} - " + "stale golden; regenerate via generate_lpdid_nonabsorbing_golden.do." + ) + return golden + + +@pytest.fixture(scope="module") +def r_first_entry() -> dict: + _skip_if_missing() + if not _R_FIXTURE_AVAILABLE: + pytest.skip( + "R non-absorbing golden not present (gate 1 only; the Stata-vs-library " + "gates run without it). Regenerate via benchmarks/R/generate_lpdid_golden.R." + ) + with R_GOLDEN_PATH.open() as f: + return json.load(f)["first_entry"] + + +@pytest.fixture(scope="module") +def panel() -> pd.DataFrame: + _skip_if_missing() + return pd.read_csv(PANEL_PATH) + + +@pytest.fixture(scope="module") +def subsample(panel: pd.DataFrame) -> pd.DataFrame: + """Pandas recompute of the generator's drop rule, sub-class by sub-class.""" + df = panel.sort_values(["unit", "time"]) + grp = df.groupby("unit") + unit_min_treat = grp["treat"].min() + always = sorted(unit_min_treat[unit_min_treat == 1].index) + d_d = grp["treat"].diff() + d_d_lag_l = d_d.groupby(df["unit"]).shift(L) + respell = sorted(df.loc[(d_d == 1) & (d_d_lag_l == -1), "unit"].unique()) + assert ( + always == ALWAYS_TREATED_UNITS + ), f"always-treated recompute {always} != expected {ALWAYS_TREATED_UNITS}" + assert ( + respell == EXACT_L_RESPELL_UNITS + ), f"exact-L respell recompute {respell} != expected {EXACT_L_RESPELL_UNITS}" + sub = panel[~panel["unit"].isin(always + respell)].reset_index(drop=True) + assert len(sub) == SUBSAMPLE_N_ROWS + return sub + + +def _fit(data: pd.DataFrame, mode: str, reweight: bool, only_event: bool = False) -> dict: + """One library fit; returns per-horizon + pooled tables and captured warnings.""" + kwargs = { + "pre_window": 3, + "post_window": 4, + "cluster": "unit", + "non_absorbing": mode, + "reweight": reweight, + } + if mode == "effect_stabilization": + kwargs["stabilization_window"] = L + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + res = LPDiD(**kwargs).fit( + data, + outcome="y", + unit="unit", + time="time", + treatment="treat", + only_event=only_event, + ) + messages = [str(w.message) for w in caught] + es = res.event_study.set_index("horizon") + by_h = { + int(h): { + "att": float(es.loc[h, "coefficient"]), + "se": float(es.loc[h, "se"]), + "n_obs": int(es.loc[h, "n_obs"]), + } + for h in ES_HORIZONS + } + pooled = {} + if not only_event: + for window in ("pre", "post"): + row = res.pooled.loc[res.pooled["window"] == window].iloc[0] + pooled[window] = { + "att": float(row["coefficient"]), + "se": float(row["se"]), + "n_obs": int(row["n_obs"]), + } + return {"by_h": by_h, "pooled": pooled, "warnings": messages} + + +@pytest.fixture(scope="module") +def fits(panel: pd.DataFrame, subsample: pd.DataFrame) -> dict: + """The five library fits (see module docstring): keyed ('fe'|'es'|'es_full', mode).""" + return { + ("fe", "vw"): _fit(panel, "first_entry", reweight=False), + ("fe", "rw"): _fit(panel, "first_entry", reweight=True), + ("es", "vw"): _fit(subsample, "effect_stabilization", reweight=False), + ("es", "rw"): _fit(subsample, "effect_stabilization", reweight=True), + ("es_full", "vw"): _fit(panel, "effect_stabilization", reweight=False, only_event=True), + } + + +def _stata_block(golden: dict, arm: str, mode: str) -> dict: + key = {"fe": "first_entry", "es": "effect_stab_sub"}[arm] + return golden[key][mode] + + +# --------------------------------------------------------------------------- +# Gate 1 +# --------------------------------------------------------------------------- +@pytest.mark.parametrize("h", ES_HORIZONS) +def test_stata_first_entry_vw_matches_r_golden(h, stata_golden, r_first_entry): + """Gate 1: the package's end-to-end Eq. 12 run reproduces the R feols anchor. + + Committed vs committed; ``rel=0`` makes the bound purely absolute. + """ + stata = _stata_block(stata_golden, "fe", "vw")["es"][str(h)] + r_att, r_se = r_first_entry[str(h)] + assert stata["att"] == pytest.approx( + r_att, abs=R_ANCHOR_ATOL, rel=0 + ), f"h={h}: Stata lpdid att {stata['att']} != R anchor {r_att}" + assert stata["se"] == pytest.approx( + r_se, abs=R_ANCHOR_ATOL, rel=0 + ), f"h={h}: Stata lpdid se {stata['se']} != R anchor {r_se}" + + +# --------------------------------------------------------------------------- +# Gates 2-3 (att + SE parity) +# --------------------------------------------------------------------------- +@pytest.mark.parametrize("mode", MODES) +@pytest.mark.parametrize("h", ES_HORIZONS) +def test_library_first_entry_matches_stata(h, mode, stata_golden, fits): + """Gate 2 (ES): library ``first_entry`` == Stata lpdid, every horizon incl. placebos.""" + stata = _stata_block(stata_golden, "fe", mode)["es"][str(h)] + lib = fits[("fe", mode)]["by_h"][h] + assert lib["att"] == pytest.approx( + stata["att"], abs=ATT_ATOL, rel=0 + ), f"first_entry {mode} h={h}: library att {lib['att']} != Stata {stata['att']}" + assert lib["se"] == pytest.approx( + stata["se"], abs=SE_ATOL, rel=0 + ), f"first_entry {mode} h={h}: library se {lib['se']} != Stata {stata['se']}" + + +@pytest.mark.parametrize("mode", MODES) +@pytest.mark.parametrize("window", ["pre", "post"]) +def test_library_first_entry_pooled_matches_stata(window, mode, stata_golden, fits): + """Gate 2 (pooled): both Eq. 12 pooled windows, att + SE - first external anchor.""" + stata = _stata_block(stata_golden, "fe", mode)["pooled"][window] + lib = fits[("fe", mode)]["pooled"][window] + assert lib["att"] == pytest.approx( + stata["att"], abs=ATT_ATOL, rel=0 + ), f"first_entry {mode} pooled {window}: library att {lib['att']} != Stata {stata['att']}" + assert lib["se"] == pytest.approx( + stata["se"], abs=SE_ATOL, rel=0 + ), f"first_entry {mode} pooled {window}: library se {lib['se']} != Stata {stata['se']}" + + +@pytest.mark.parametrize("mode", MODES) +@pytest.mark.parametrize("h", POST_HORIZONS) +def test_library_effect_stab_subsample_matches_stata(h, mode, stata_golden, fits): + """Gate 3 (ES): library Eq. 13 == Stata lpdid on the convention-neutral subsample. + + Post horizons only - the placebo rows diverge by construction (class 3; gate 6b). + The rw SEs here are the previously pinned-only surface, now externally anchored. + """ + stata = _stata_block(stata_golden, "es", mode)["es"][str(h)] + lib = fits[("es", mode)]["by_h"][h] + assert lib["att"] == pytest.approx( + stata["att"], abs=ATT_ATOL, rel=0 + ), f"effect_stab {mode} h={h}: library att {lib['att']} != Stata {stata['att']}" + assert lib["se"] == pytest.approx( + stata["se"], abs=SE_ATOL, rel=0 + ), f"effect_stab {mode} h={h}: library se {lib['se']} != Stata {stata['se']}" + + +@pytest.mark.parametrize("mode", MODES) +def test_library_effect_stab_pooled_post_matches_stata(mode, stata_golden, fits): + """Gate 3 (pooled post): the Eq. 13 pooled-post window, att + SE.""" + stata = _stata_block(stata_golden, "es", mode)["pooled"]["post"] + lib = fits[("es", mode)]["pooled"]["post"] + assert lib["att"] == pytest.approx( + stata["att"], abs=ATT_ATOL, rel=0 + ), f"effect_stab {mode} pooled post: library att {lib['att']} != Stata {stata['att']}" + assert lib["se"] == pytest.approx( + stata["se"], abs=SE_ATOL, rel=0 + ), f"effect_stab {mode} pooled post: library se {lib['se']} != Stata {stata['se']}" + + +# --------------------------------------------------------------------------- +# Gate 4 (sample shape - parity surfaces only; see module docstring on N_clust) +# --------------------------------------------------------------------------- +@pytest.mark.parametrize("mode", MODES) +def test_sample_shape_first_entry(mode, stata_golden, fits): + """Gate 4 (Eq. 12): Stata obs == library n_obs, all horizons + both pooled windows.""" + stata = _stata_block(stata_golden, "fe", mode) + lib = fits[("fe", mode)] + for h in ES_HORIZONS: + assert int(stata["es"][str(h)]["N"]) == lib["by_h"][h]["n_obs"], ( + f"first_entry {mode} h={h}: Stata obs {stata['es'][str(h)]['N']} != " + f"library n_obs {lib['by_h'][h]['n_obs']}" + ) + for window in ("pre", "post"): + assert ( + int(stata["pooled"][window]["N"]) == lib["pooled"][window]["n_obs"] + ), f"first_entry {mode} pooled {window}: sample sizes differ" + + +@pytest.mark.parametrize("mode", MODES) +def test_sample_shape_effect_stab_post(mode, stata_golden, fits): + """Gate 4 (Eq. 13): obs equality on post horizons + pooled post only. + + Pre horizons are excluded for BOTH modes: those samples differ by construction + (class 3), and gate 6(b) asserts the mismatch instead. + """ + stata = _stata_block(stata_golden, "es", mode) + lib = fits[("es", mode)] + for h in POST_HORIZONS: + assert int(stata["es"][str(h)]["N"]) == lib["by_h"][h]["n_obs"], ( + f"effect_stab {mode} h={h}: Stata obs {stata['es'][str(h)]['N']} != " + f"library n_obs {lib['by_h'][h]['n_obs']}" + ) + assert ( + int(stata["pooled"]["post"]["N"]) == lib["pooled"]["post"]["n_obs"] + ), f"effect_stab {mode} pooled post: sample sizes differ" + + +# --------------------------------------------------------------------------- +# Gate 5 (drop-rule consistency) +# --------------------------------------------------------------------------- +def test_drop_rule_consistency(stata_golden, subsample): + """Gate 5: golden provenance == the pandas recompute, sub-class by sub-class. + + The ``subsample`` fixture already asserted the two sub-class unit lists (locking + the per-class attribution) and the 658-row count; this test pins the golden's + record of them plus the free-text rule string against hard-coded twins. + """ + meta = stata_golden["meta"] + assert meta["drop_rule"] == DROP_RULE + assert meta["dropped_units"] == sorted(ALWAYS_TREATED_UNITS + EXACT_L_RESPELL_UNITS) + assert meta["n_rows_sub"] == SUBSAMPLE_N_ROWS == len(subsample) + assert stata_golden["effect_stab_sub"]["n_rows"] == SUBSAMPLE_N_ROWS + + +# --------------------------------------------------------------------------- +# Gate 6 (convention-divergence documentation) +# --------------------------------------------------------------------------- +@pytest.mark.parametrize("h", POST_HORIZONS) +def test_full_panel_divergence_documented(h, stata_golden, fits): + """Gate 6(a): Eq. 13 FULL-panel Stata != library - classes 1+2, jointly. + + Asserts both the att divergence (floor 1e-3; measured 2.35e-3 at h=0 up to + ~4e-2 at h=4) and the obs mismatch (the sample-admission mechanism itself: + measured Stata 628/544/463/380/320 vs library 601/517/436/353/291). + """ + stata = stata_golden["effect_stab_full_vw"]["es"][str(h)] + lib = fits[("es_full", "vw")]["by_h"][h] + att_diff = abs(stata["att"] - lib["att"]) + assert att_diff > DIVERGENCE_FLOOR, ( + f"full-panel h={h}: att divergence {att_diff} fell below the documented floor " + "- the package's boundary/respell conventions may have changed; re-measure " + "and update REGISTRY LPDiD Deviation 4." + ) + assert int(stata["N"]) != lib["n_obs"], ( + f"full-panel h={h}: Stata obs == library n_obs ({lib['n_obs']}) - the " + "classes-1+2 sample-admission difference disappeared; re-measure." + ) + + +@pytest.mark.parametrize("mode", MODES) +def test_subsample_placebo_divergence_documented(mode, stata_golden, fits): + """Gate 6(b): Eq. 13 subsample placebo rows + pooled Pre diverge - class 3. + + Asserts att divergence (floor 1e-3; measured 3.9e-3..7.9e-2) and obs mismatch + (measured ES vw 342/408 vs 390/437, rw 173/220 vs 282/329; pooled Pre vw 342 vs + 390, rw 173 vs 282) for BOTH modes. + """ + stata = _stata_block(stata_golden, "es", mode) + lib = fits[("es", mode)] + for h in PRE_HORIZONS: + s_row = stata["es"][str(h)] + l_row = lib["by_h"][h] + att_diff = abs(s_row["att"] - l_row["att"]) + assert att_diff > DIVERGENCE_FLOOR, ( + f"subsample {mode} h={h}: placebo att divergence {att_diff} fell below " + "the documented floor - the package's placebo-window convention may have " + "changed; re-measure and update REGISTRY LPDiD Deviation 4." + ) + assert ( + int(s_row["N"]) != l_row["n_obs"] + ), f"subsample {mode} h={h}: placebo samples coincide - class 3 disappeared." + s_pre = stata["pooled"]["pre"] + l_pre = lib["pooled"]["pre"] + assert abs(s_pre["att"] - l_pre["att"]) > DIVERGENCE_FLOOR + assert int(s_pre["N"]) != l_pre["n_obs"] + + +# --------------------------------------------------------------------------- +# Gate 7 (warning cleanliness) +# --------------------------------------------------------------------------- +def test_all_fits_warning_clean(fits): + """Gate 7: none of the five library fits emits any warning.""" + noisy = {key: f["warnings"] for key, f in fits.items() if f["warnings"]} + assert not noisy, f"unexpected warning(s): {noisy}" + + +# --------------------------------------------------------------------------- +# Gate 8 (provenance) +# --------------------------------------------------------------------------- +def test_provenance_contract(stata_golden): + """Gate 8: drift signals + load-bearing meta strings. + + ``ssc_versions`` must cover exactly the generator's fail-closed dependency set + with real identifiers - never empty/unknown/MISSING, and each embedding the ado + file's ``checksum:`` marker (the generator always appends checksum + length to + the version text: SSC has no immutable archive, so a same-version-string + upstream edit must still move this metadata). ``mapping`` is the load-bearing + option<->mode claim; the prose fields are asserted non-empty. + """ + meta = stata_golden["meta"] + versions = meta["ssc_versions"] + assert set(versions) == SSC_VERSION_KEYS + for pkg, ver in versions.items(): + assert isinstance(ver, str) and ver.strip(), f"{pkg}: empty version" + assert ver not in {"unknown", "MISSING"}, f"{pkg}: no usable drift signal" + assert re.search( + r"checksum:\d+ len:\d+", ver + ), f"{pkg}: version string lacks the numeric checksum+len drift signal: {ver!r}" + assert meta["mapping"] == MAPPING + assert meta["source_panel"] == "benchmarks/data/lpdid_nonabsorbing_panel.csv" + sha = hashlib.sha256(PANEL_PATH.read_bytes()).hexdigest() + assert meta["source_sha256"] == sha, ( + "golden source_sha256 does not match the committed panel - the panel changed " + "without regenerating the Stata golden." + ) + for field in ("cmd", "point_anchor", "se_convention", "convention_notes"): + assert isinstance(meta[field], str) and meta[field].strip() diff --git a/tests/test_methodology_lpdid.py b/tests/test_methodology_lpdid.py index c5fe91e0..7bf93a0c 100644 --- a/tests/test_methodology_lpdid.py +++ b/tests/test_methodology_lpdid.py @@ -52,7 +52,9 @@ # INDEPENDENT fixest::feols reconstruction of the paper's Eq. 12 (first_entry) and # Eq. 13 (effect_stabilization) clean-sample restrictions; alexCardazzi's # nonabsorbing_lag is recorded in the golden meta as a divergent reference, NOT a -# parity gate (it clamps off-switches + uses a non-paper boundary/window convention). +# parity gate (its off-switch clamp violates Eq. 13's stated both-directions control +# condition, and its boundary/placebo-window conventions differ on surfaces the +# paper does not specify - see REGISTRY ## LPDiD Deviation 4). NA_GOLDEN_PATH = _DATA / "lpdid_nonabsorbing_golden.json" NA_PANEL_PATH = _DATA / "lpdid_nonabsorbing_panel.csv" _NONABSORB_FIXTURE_AVAILABLE = NA_GOLDEN_PATH.is_file() and NA_PANEL_PATH.is_file() @@ -91,6 +93,14 @@ # the reweighted POINT matches the independent feols recipe to ~1e-13, but the # weighted-cluster SE has a small feols convention difference (~5e-5), so the library # value is pinned as a regression guard (the variance-weighted SEs ARE feols-parity). +# The rw SE *convention* is now Stata-anchored (authors' `lpdid` package, ~1e-9) at +# POST horizons + pooled post via the convention-neutral subsample in +# tests/test_lpdid_nonabsorbing_stata_parity.py - but the FULL-panel values pinned +# HERE are not that surface: the full panel diverges from the package by the +# boundary/exact-L-respell sample conventions (its gate 6a), and the h=-3/-2 entries +# sit on the placebo-window convention difference, not anchored by the Stata arm +# (REGISTRY ## LPDiD Deviation 4). This pin therefore stays as the CI-runnable guard +# for the full-panel surface (RA_SE_PIN pattern); Stata is not available in CI. # Refresh from the library if the committed non-absorbing panel changes: fit # LPDiD(pre_window=3, post_window=4, reweight=True, cluster="unit", # non_absorbing="effect_stabilization", stabilization_window=3) on @@ -309,9 +319,14 @@ class TestLPDiDNonAbsorbingParityR: """Pin the C1 non-absorbing modes against the independent feols Eq. 12/13 recipes. The variance-weighted variants match the feols golden on point AND SE; the - reweighted point matches feols while its SE is a pinned regression guard - (documented weighted-cluster convention difference). alexCardazzi's - ``nonabsorbing_lag`` is recorded in the golden meta as a divergent reference only. + reweighted point matches feols while its SE has a documented weighted-cluster + feols convention difference (~5e-5) and is held here by the ``RW_SE_PIN`` + regression guard. The rw SE convention is externally anchored against the + authors' Stata ``lpdid`` package (~1e-9) on its parity surfaces in + ``tests/test_lpdid_nonabsorbing_stata_parity.py`` - see the ``RW_SE_PIN`` + comment for why the full-panel/placebo values pinned here remain pinned-only. + alexCardazzi's ``nonabsorbing_lag`` is recorded in the golden meta as a + divergent reference only. """ def _meta(self, na_golden: dict): @@ -365,7 +380,9 @@ def test_effect_stabilization_reweighted_point( def test_alex_recorded_as_divergent_reference(self, na_golden: dict) -> None: """alexCardazzi nonabsorbing_lag is recorded (meta), NOT a parity gate: it - diverges from the paper-faithful Eq.13 even on this panel.""" + diverges from the library's Eq.13 implementation even on this panel (the + off-switch clamp violates Eq.13's stated control condition; its boundary/ + placebo-window conventions differ on paper-silent surfaces).""" alex = na_golden["alex_nonabsorbing_lag"] assert "error" not in alex, f"alex reference failed to generate: {alex}" es = na_golden["effect_stab"]