diff --git a/.gitignore b/.gitignore index c4dc5889..38c9fefc 100644 --- a/.gitignore +++ b/.gitignore @@ -122,6 +122,6 @@ benchmarks/refresh_2026_07/results/raw/ benchmarks/refresh_2026_07/results/refresh_results_smoke.json # Stata batch mode writes .log to the cwd (repo root) - never committed; -# the golden JSON it produces IS committed (benchmarks/data/lpdid_ra_stata_golden.json). -# Root-anchored to the generator's own log so it doesn't hide logs elsewhere. -/generate_lpdid_ra_golden.log +# the golden JSON each generator produces IS committed (benchmarks/data/*.json). +# Root-anchored to the generators' own logs so it doesn't hide logs elsewhere. +/generate_*_golden.log diff --git a/CHANGELOG.md b/CHANGELOG.md index f13c13c8..2c431a9b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,18 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] ### Added +- **ImputationDiD leave-one-out SE now anchored against Stata `did_imputation, leaveout` + (no library behavior change).** The Borusyak-Jaravel-Spiess (2024) App. A.9 LOO variance + (`leave_one_out=True`) has no runnable R reference (R `didimputation` omits LOO), so it + was validated only by an internal psi-identity + hand-calc + MC coverage. It is now + cross-validated against the authors' own Stata `did_imputation ..., leaveout` via a + committed golden (`benchmarks/stata/generate_imputation_loo_golden.do` → + `didimputation_loo_stata_golden.json`; `tests/test_imputation_loo_stata_parity.py`): + the library LOO SE matches to ~1e-9 (cross-implementation; the point to ~2e-8, gated at + `abs=1e-7`) at the overall ATT and all 6 event-study horizons, with the non-LOO SE + confirmed three-way (library = R = Stata). + Second Stata parity arm and the first SSC-dependent one — `benchmarks/stata/requirements.do` + documents the one-time install; the golden records `ssc_versions` for drift detection. - **LP-DiD regression-adjustment SE now anchored against Stata `teffects` (no library behavior change).** The RA influence-function cluster SE has no runnable R reference (`alexCardazzi/lpdid` does direct covariate inclusion, not RA), so it was previously diff --git a/TODO.md b/TODO.md index e15bd630..a7fbf00c 100644 --- a/TODO.md +++ b/TODO.md @@ -51,3 +51,4 @@ generic sparse-FE, QR+SVD rank-detection redundancy, `check_finite` bypass — m | `worktree-rm` safety via a tested argv helper: the prose rewrite (ask-before-remove confirmation gate, detached-HEAD reachability/rescue, tip-identity force-delete guard) was reverted to the main version because editing the prose repeatedly reintroduced shell-injection (last: sourcing a state file built from a git-derived branch name). Restore those guards in a `worktree_rm.py` that takes the name via file ingress, invokes git through argv arrays, and has a metacharacter-branch/path injection regression test — the pattern that worked for `pr_prepare.py`/`premerge_scan.py`. | `.claude/commands/worktree-rm.md`, `.claude/scripts/` | skill-audit | Heavy | Medium | | `premerge_scan.py` should scan the staged blob (`git show :path`) for staged methodology files, not the working-tree copy — a stage-then-revert-working-copy edit currently reads the safe working version and misses the staged violation. Union staged-index findings with unstaged/untracked filesystem findings. | `.claude/scripts/premerge_scan.py` | skill-audit | Mid | Low | | Re-add committed-range methodology scanning to `/push-pr-update` §3b (clean tree, commits ahead) using `premerge_scan.py --range`, with the comparison ref passed as **data** (resolved into a quoted variable in one Bash call, never a raw ``). It was removed to avoid ref interpolation; the helper already implements and tests `--range`. | `.claude/commands/push-pr-update.md` | skill-audit | Quick | Low | +| Extend the ImputationDiD-LOO Stata anchor to the coarser `aux_partition ∈ {"cohort", "horizon"}` variants (Stata `did_imputation ..., leaveout avgeffectsby(Ei)` / `avgeffectsby(K)`, `K = t - Ei`). These differ from the default only at the *overall* aggregate (a no-op per-horizon) and are the least-validated surface — no R analogue, currently only hand-calc. Extends `benchmarks/stata/generate_imputation_loo_golden.do` + golden. | `benchmarks/stata/`, `tests/test_imputation_loo_stata_parity.py` | stata-arm | Mid | Low | diff --git a/benchmarks/README.md b/benchmarks/README.md index 3a5428be..dc667d55 100644 --- a/benchmarks/README.md +++ b/benchmarks/README.md @@ -51,12 +51,18 @@ and CI never needs Stata. Generators live in `stata/` and are run headless: ```bash # macOS, StataSE 19 (binary is not on PATH by default) -/Applications/Stata/StataSE.app/Contents/MacOS/stata-se -b do \ - benchmarks/stata/generate_lpdid_ra_golden.do +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_imputation_loo_golden.do ``` -Current arm uses only **native** Stata commands (`teffects`) — no SSC packages — -so `version 19` pins behavior. See `stata/README.md`. +The `LPDiD` arm uses only **native** Stata commands (`teffects`), pinned by +`version 19`. The `ImputationDiD` arm depends on SSC packages +(`did_imputation`/`reghdfe`/`ftools`/`require`), which `version 19` does NOT pin +(SSC has no version history) — install them once via `requirements.do` (the +generators do not auto-install) and each golden records `ssc_versions` for drift +detection. See `stata/README.md`. ## Directory Structure @@ -72,8 +78,10 @@ benchmarks/ │ ├── benchmark_honest.R # HonestDiD │ └── benchmark_fixest.R # Basic DiD / TWFE ├── stata/ -│ ├── README.md # Stata arm docs -│ └── generate_lpdid_ra_golden.do # LPDiD RA SE vs teffects ra +│ ├── 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_imputation_loo_golden.do # ImputationDiD LOO SE vs did_imputation leaveout ├── python/ │ ├── utils.py # Common utilities │ ├── benchmark_callaway.py # CallawaySantAnna @@ -96,6 +104,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 | +| `ImputationDiD` (LOO SE) | Stata `did_imputation, leaveout` | Borusyak, Jaravel & Spiess (2024) App. A.9 | ✓ Integrated | | `HonestDiD` | `HonestDiD::createSensitivityResults` | Rambachan & Roth (2023) | Planned | Note: HonestDiD benchmark scripts exist but are not yet integrated into the main runner. diff --git a/benchmarks/data/didimputation_loo_stata_golden.json b/benchmarks/data/didimputation_loo_stata_golden.json new file mode 100644 index 00000000..e3396a14 --- /dev/null +++ b/benchmarks/data/didimputation_loo_stata_golden.json @@ -0,0 +1,28 @@ +{ + "meta": { + "estimator": "ImputationDiD leave-one-out (BJS 2024 App. A.9) SE - Stata did_imputation leaveout", + "generator": "benchmarks/stata/generate_imputation_loo_golden.do", + "source_panel": "benchmarks/data/didimputation_test_panel.csv", + "point_anchor": "benchmarks/data/didimputation_golden.json (overall.att, event_study.att)", + "cmd": "did_imputation y unit time Ei, [horizons(0/5)] leaveout avgeffectsby(Ei t) cluster(unit)", + "avgeffectsby": "Ei t (== library aux_partition=cohort_horizon; pinned explicitly)", + "se_convention": "A.9 leave-one-out finite-sample variance; unit-clustered; se_nonloo is the non-leaveout cluster SE (== R didimputation).", + "ssc_versions": { + "did_imputation": "Version: November 22, 2023", + "reghdfe": "version 6.13.1 10Jan2026", + "ftools": "version 2.50.0 09jan2026", + "require": "version 1.3.1 19sep2023" + }, + "stata_edition": "SE", + "stata_version": 19.0 + }, + "overall": {"att": 2.0456678967434003, "se": 0.021041161682015575, "se_nonloo": 0.020879383566722923, "N": 1440}, + "event_study": { + "0": {"att": 0.99293658500530635, "se": 0.029085929067998187, "se_nonloo": 0.028791608348336495}, + "1": {"att": 1.5129648623314707, "se": 0.02548745171356798, "se_nonloo": 0.025253856015380709}, + "2": {"att": 2.006499293218484, "se": 0.030371068022456742, "se_nonloo": 0.030155388323006894}, + "3": {"att": 2.4998545630646829, "se": 0.029299586939814812, "se_nonloo": 0.029048460851672311}, + "4": {"att": 2.9648169454740363, "se": 0.043460391002896281, "se_nonloo": 0.043064691747540094}, + "5": {"att": 3.4673514576204787, "se": 0.043504637171120686, "se_nonloo": 0.043073537575459654} + } +} diff --git a/benchmarks/stata/README.md b/benchmarks/stata/README.md index b9bc6891..6a5a6828 100644 --- a/benchmarks/stata/README.md +++ b/benchmarks/stata/README.md @@ -115,19 +115,89 @@ sample matches the library's realized size and cluster count — which, together the point (1e-10) and SE (~1e-16) agreement, strongly corroborates the same clean sample. No timestamp — the golden regenerates byte-identically. +--- + +# `did_imputation, leaveout` parity for ImputationDiD LOO SE + +`benchmarks/stata/generate_imputation_loo_golden.do` produces +`benchmarks/data/didimputation_loo_stata_golden.json`, consumed by +`tests/test_imputation_loo_stata_parity.py` to verify that +`ImputationDiD(leave_one_out=True)` — the Borusyak-Jaravel-Spiess (2024) Supplementary +Appendix A.9 finite-sample variance refinement — matches Stata `did_imputation ..., +leaveout` at the overall ATT and all 6 event-study horizons. + +## Why Stata `did_imputation` + +The A.9 leave-one-out (LOO) variance has **no runnable R reference** — R +`didimputation` omits LOO entirely — so the library LOO SE was validated only by an +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. + +## This is the first SSC-dependent arm + +Unlike the native-`teffects` LPDiD arm, `did_imputation` is an SSC package with a +dependency chain `did_imputation → reghdfe → require + ftools`, none pinned by +`version 19`. 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. + +## What the generator does + +No clean-sample reconstruction: `did_imputation` consumes the raw committed R-arm panel +`didimputation_test_panel.csv` (180 units, cohorts {3,5} + never-treated); the only +mapping is `Ei = first_treat` (missing for never-treated). The partition is pinned +`avgeffectsby(Ei t)` **explicitly** (== library `aux_partition="cohort_horizon"`). This +is also `did_imputation`'s current default, but pinning it keeps the validation estimand +self-describing and robust to a future default change. The per-horizon command is: + +```stata +did_imputation y unit time Ei, horizons(0/5) leaveout avgeffectsby(Ei t) cluster(unit) +``` + +Agreement is cross-implementation, not bit-identical (`did_imputation` goes through +`reghdfe`, the library through its own sparse IF solver): the **SE** agrees to ~1e-9 and +the **point** to ~2e-8; the parity test gates at `abs=1e-7`. Both the LOO and the non-LOO +cluster SE are emitted; the non-LOO SE additionally three-way-confirms against the R +golden (`didimputation_golden.json`). + +**Panel-coupling caveat:** `didimputation_test_panel.csv` is owned by the R generator +(`benchmarks/R/generate_didimputation_golden.R`); if it is ever regenerated, the R golden +AND this Stata golden must both be regenerated. + +## Regenerating + +```bash +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_imputation_loo_golden.do +grep -E '^r\([0-9]+\);' generate_imputation_loo_golden.log # must print nothing +``` + +## JSON schema + +```json +{ + "meta": { "...": "...", "avgeffectsby": "Ei t (== cohort_horizon)", + "ssc_versions": {"did_imputation": "...", "reghdfe": "...", "...": "..."} }, + "overall": {"att": , "se": , "se_nonloo": , "N": }, + "event_study": { "0": {"att": , "se": , "se_nonloo": }, "...": "..." } +} +``` + ## 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 generator also runs - an in-`.do` point smoke gate at 1e-8 that surfaces as `r(9);` on a gross port - bug; the Python test's 1e-10 gate is authoritative). + shell exit code — parse the `.log` for `^r\([0-9]+\);`. Each generator also runs an + informational in-`.do` point smoke gate (LPDiD 1e-8, ImputationDiD 1e-6) that + surfaces as `r(9);` on a gross bug; the Python parity test is authoritative (LPDiD + point 1e-10; ImputationDiD `abs=1e-7`). - **`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, with BE by elimination (`c(BE)` is undefined and `cond()` evaluates all branches 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. This arm is exempt because `teffects` is - native (pinned by `version 19`), but any *future* SSC-dependent Stata generator - must record `which ` output verbatim into its golden's `meta` so drift is - at least detectable. + no lockfile / archive to pin against. The LPDiD arm is exempt (`teffects` is native, + pinned by `version 19`); the ImputationDiD arm records its 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_imputation_loo_golden.do b/benchmarks/stata/generate_imputation_loo_golden.do new file mode 100644 index 00000000..bf186a36 --- /dev/null +++ b/benchmarks/stata/generate_imputation_loo_golden.do @@ -0,0 +1,235 @@ +*! Golden generator: ImputationDiD leave-one-out (A.9) SE vs Stata `did_imputation, leaveout`. +*! +*! Purpose +*! Produce an INDEPENDENT external anchor for the ImputationDiD leave-one-out (LOO) +*! standard error - the Borusyak-Jaravel-Spiess (2024) Supplementary Appendix A.9 +*! finite-sample variance refinement (opt-in `leave_one_out=True`). No R package +*! computes it (R `didimputation` omits LOO), so the library LOO SE was validated +*! only by an internal psi-identity + hand-calc + MC coverage. The authors' own +*! Stata `did_imputation` ships the same option (`leaveout`); this generator turns +*! that into a measured cross-implementation anchor. +*! +*! SECOND Stata arm, and the FIRST SSC-dependent one. `did_imputation` is NOT +*! pinned by `version 19` (SSC has no version history), so this generator records +*! the installed package versions in meta.ssc_versions for drift detection. It does +*! NOT install anything - run benchmarks/stata/requirements.do once first. +*! +*! Consuming test +*! tests/test_imputation_loo_stata_parity.py +*! +*! Outputs (checked into the repo) +*! benchmarks/data/didimputation_loo_stata_golden.json +*! +*! Usage (run from the repo root, AFTER benchmarks/stata/requirements.do) +*! /Applications/Stata/StataSE.app/Contents/MacOS/stata-se -b do \ +*! benchmarks/stata/generate_imputation_loo_golden.do +*! Then confirm the log is clean: grep -E '^r\([0-9]+\);' generate_imputation_loo_golden.log +*! (Stata batch mode ALWAYS exits 0, even on error - never trust the exit code.) +*! +*! Notes +*! - Reads (does NOT regenerate) benchmarks/data/didimputation_test_panel.csv, whose +*! sole owner is benchmarks/R/generate_didimputation_golden.R (180 units, cohorts +*! {3,5} + never-treated, t=1..8). If that panel is ever regenerated, the R golden +*! AND this Stata golden must both be regenerated. +*! - No clean-sample reconstruction: did_imputation consumes the raw panel; the only +*! mapping is Ei = first_treat (missing for never-treated). +*! - avgeffectsby(Ei t) is passed EXPLICITLY (== library aux_partition="cohort_horizon"). +*! did_imputation currently defaults to avgeffectsby(Ei t) too; pinning it explicitly +*! makes the validation estimand self-describing and robust to a future default change. +*! - The in-.do point gate (1e-6, informational) checks the Stata coef against the +*! R-anchored points from didimputation_golden.json and aborts early on a gross bug. +*! The AUTHORITATIVE parity gate is tests/test_imputation_loo_stata_parity.py, which +*! compares the recomputed library output against this committed Stata golden. + +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 ".021"/"-.5" (leading dot); JSON requires a leading 0, so we +* patch ".x" -> "0.x" and "-.x" -> "-0.x". Returns r(s). +capture program drop _jnum +program define _jnum, rclass + args x fmt + 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 + +* Extract a JSON-safe version string for an installed ado. Scans the leading `*!` +* header block and PREFERS a line containing "version" (SSC packages put the real +* version/date there - e.g. did_imputation's is on line 2, "Version: ...", after the +* line-1 description); falls back to the first `*!` line. Sanitized (quotes/backslashes +* /tabs stripped, truncated). Returns r(v)="MISSING" if the ado is absent. +capture program drop _adover +program define _adover, rclass + args pkg + capture findfile `pkg'.ado + if _rc { + return local v "MISSING" + exit + } + local fn = r(fn) + tempname vh + local ver "" + 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 + } + } + local ++n + file read `vh' line + } + file close `vh' + if `"`ver'"' == "" local ver `"`first'"' + if `"`ver'"' == "" local ver "unknown" + * JSON-sanitize: drop double-quotes/backslashes/tabs, truncate. + local ver = subinstr(`"`ver'"', `"""', "'", .) + local ver = subinstr(`"`ver'"', "\", "/", .) + local ver = subinstr(`"`ver'"', char(9), " ", .) + if length(`"`ver'"') > 100 local ver = substr(`"`ver'"', 1, 100) + return local v `"`ver'"' +end + +* ------------------------------------------------------------------------------ +* Dependency check + version capture (do NOT install; see requirements.do). +* ------------------------------------------------------------------------------ +foreach p in ftools require reghdfe did_imputation { + capture which `p' + if _rc { + di as error "`p' is not installed. Run: stata -b do benchmarks/stata/requirements.do" + error 111 + } +} +_adover ftools +local v_ftools = r(v) +_adover require +local v_require = r(v) +_adover reghdfe +local v_reghdfe = r(v) +_adover did_imputation +local v_didimp = r(v) +* Each captured identifier must be a real version/date line (not the bare description), +* else drift would go undetected (review: did_imputation's version is on the 2nd *! line). +assert strpos(lower(`"`v_didimp'"'), "version") > 0 +assert strpos(lower(`"`v_reghdfe'"'), "version") > 0 +assert strpos(lower(`"`v_ftools'"'), "version") > 0 +assert strpos(lower(`"`v_require'"'), "version") > 0 + +* ------------------------------------------------------------------------------ +* R-anchored POINT estimates (benchmarks/data/didimputation_golden.json, digits=12). +* Hard-coded for the in-.do smoke gate only; the Python test reads the JSON and is +* authoritative. Points are identical for LOO and non-LOO (LOO only changes the SE). +* ------------------------------------------------------------------------------ +scalar rpt_all = 2.045668026901 +scalar rpt_0 = 0.9929366801146 +scalar rpt_1 = 1.512964957441 +scalar rpt_2 = 2.006499440277 +scalar rpt_3 = 2.499854710124 +scalar rpt_4 = 2.964817090473 +scalar rpt_5 = 3.46735160262 + +* ------------------------------------------------------------------------------ +* Load the committed panel; map first_treat -> Ei (missing for never-treated). +* ------------------------------------------------------------------------------ +import delimited using "benchmarks/data/didimputation_test_panel.csv", clear varnames(1) +confirm numeric variable unit time first_treat y +quietly count +assert r(N) == 1440 // 180 units x 8 periods +gen Ei = first_treat +replace Ei = . if first_treat == 0 + +* ------------------------------------------------------------------------------ +* Overall ATT: LOO + non-LOO. Coefficient is named `tau`. +* ------------------------------------------------------------------------------ +did_imputation y unit time Ei, leaveout avgeffectsby(Ei t) cluster(unit) +scalar att_all = _b[tau] +scalar se_all = _se[tau] +scalar N_all = e(N) +assert reldif(att_all, rpt_all) < 1e-6 + +did_imputation y unit time Ei, avgeffectsby(Ei t) cluster(unit) +scalar se_all_nl = _se[tau] + +* ------------------------------------------------------------------------------ +* Event study horizons 0..5: LOO + non-LOO. Coefficients tau0..tau5. +* ------------------------------------------------------------------------------ +did_imputation y unit time Ei, horizons(0/5) leaveout avgeffectsby(Ei t) cluster(unit) +forvalues h = 0/5 { + scalar att_`h' = _b[tau`h'] + scalar se_`h' = _se[tau`h'] + assert reldif(att_`h', rpt_`h') < 1e-6 +} +did_imputation y unit time Ei, horizons(0/5) avgeffectsby(Ei t) cluster(unit) +forvalues h = 0/5 { + scalar se_nl_`h' = _se[tau`h'] +} + +* ------------------------------------------------------------------------------ +* Emit JSON by hand at %21.17g. Reuses the LPDiD conventions: every STRING field +* carries a trailing comma (the compound-quote `"' delimiter swallows a trailing +* double-quote), so the meta block ENDS with a numeric field. +* ------------------------------------------------------------------------------ +local sver = strtrim(string(c(stata_version), "%4.1f")) +local sedition = cond(c(MP)==1, "MP", cond(c(SE)==1, "SE", "BE")) +tempname fh +file open `fh' using "benchmarks/data/didimputation_loo_stata_golden.json", write replace text +file write `fh' "{" _n +file write `fh' `" "meta": {"' _n +file write `fh' `" "estimator": "ImputationDiD leave-one-out (BJS 2024 App. A.9) SE - Stata did_imputation leaveout","' _n +file write `fh' `" "generator": "benchmarks/stata/generate_imputation_loo_golden.do","' _n +file write `fh' `" "source_panel": "benchmarks/data/didimputation_test_panel.csv","' _n +file write `fh' `" "point_anchor": "benchmarks/data/didimputation_golden.json (overall.att, event_study.att)","' _n +file write `fh' `" "cmd": "did_imputation y unit time Ei, [horizons(0/5)] leaveout avgeffectsby(Ei t) cluster(unit)","' _n +file write `fh' `" "avgeffectsby": "Ei t (== library aux_partition=cohort_horizon; pinned explicitly)","' _n +file write `fh' `" "se_convention": "A.9 leave-one-out finite-sample variance; unit-clustered; se_nonloo is the non-leaveout cluster SE (== R didimputation).","' _n +file write `fh' `" "ssc_versions": {"' _n +file write `fh' `" "did_imputation": "`v_didimp'","' _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' `" },"' _n +file write `fh' `" "stata_edition": "`sedition'","' _n +file write `fh' `" "stata_version": `sver'"' _n +file write `fh' " }," _n + +* overall block +_jnum att_all +local a = r(s) +_jnum se_all +local s = r(s) +_jnum se_all_nl +local snl = r(s) +_jnum N_all "%12.0f" +local nn = r(s) +file write `fh' `" "overall": {"att": `a', "se": `s', "se_nonloo": `snl', "N": `nn'},"' _n + +* event_study block (comma-prefixed; no trailing comma, no blank first line) +file write `fh' `" "event_study": {"' +local sep "" +forvalues h = 0/5 { + _jnum att_`h' + local a = r(s) + _jnum se_`h' + local s = r(s) + _jnum se_nl_`h' + local snl = r(s) + file write `fh' "`sep'" _n `" "`h'": {"att": `a', "se": `s', "se_nonloo": `snl'}"' + local sep "," +} +file write `fh' _n " }" _n +file write `fh' "}" _n +file close `fh' + +display "Wrote benchmarks/data/didimputation_loo_stata_golden.json (overall + 6 horizons)" diff --git a/benchmarks/stata/requirements.do b/benchmarks/stata/requirements.do new file mode 100644 index 00000000..4359c35b --- /dev/null +++ b/benchmarks/stata/requirements.do @@ -0,0 +1,33 @@ +*! Stata SSC requirements for diff-diff parity arms. +*! +*! Run this ONCE before regenerating any SSC-dependent golden. The generators do +*! NOT auto-install: `ssc install` always fetches the LATEST version, which would +*! (a) break byte-identical regeneration across dates, (b) require network at +*! generation time, and (c) defeat the drift-detection that each golden records in +*! its meta.ssc_versions. Installing is therefore a deliberate, separate step. +*! +*! Usage (from the repo root): +*! /Applications/Stata/StataSE.app/Contents/MacOS/stata-se -b do benchmarks/stata/requirements.do +*! +*! SSC has no version history, so there is no way to pin a specific release here; +*! the recorded meta.ssc_versions in each golden is the only drift signal. +*! +*! Packages (dependency order matters - ftools/require before reghdfe): +*! ftools, require - reghdfe infrastructure +*! reghdfe - fast FE regression (did_imputation backend) +*! did_imputation - Borusyak-Jaravel-Spiess (2024) imputation DiD; the only +*! implementation of the App. A.9 leave-one-out variance +*! (consumed by generate_imputation_loo_golden.do) + +version 19 +foreach p in ftools require reghdfe did_imputation { + capture which `p' + if _rc { + di as txt "Installing `p' from SSC ..." + ssc install `p', replace + } + else { + di as txt "`p' already installed." + } +} +di as result "Stata SSC requirements satisfied." diff --git a/docs/benchmarks.rst b/docs/benchmarks.rst index b6b0dcc1..378ab14f 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). Coefficient accuracy, standard error comparison, and performance metrics. + :description: Validation benchmarks comparing diff-diff against R packages (did, synthdid, fixest) and Stata (teffects, did_imputation). 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,9 @@ where no runnable R reference exists: * - ``LPDiD`` (regression-adjustment SE) - Stata ``teffects ra ... atet`` - Dube, Girardi, Jorda & Taylor (2025); no runnable R analogue + * - ``ImputationDiD`` (leave-one-out SE) + - Stata ``did_imputation, leaveout`` + - Borusyak, Jaravel & Spiess (2024) Supp. App. A.9; no runnable R analogue Methodology ----------- @@ -1010,15 +1013,20 @@ Prerequisites pip install -e ".[dev]" -4. (Optional) Stata, only to regenerate the ``LPDiD`` regression-adjustment SE - golden. Uses the **native** ``teffects`` command (no SSC package); the golden - is committed, so this is not needed to run the test suite: +4. (Optional) Stata, only to regenerate the committed Stata goldens (``LPDiD`` + regression-adjustment SE, ``ImputationDiD`` leave-one-out SE). The goldens are + committed, so this is not needed to run the test suite. The ``LPDiD`` arm uses the + **native** ``teffects`` command; the ``ImputationDiD`` arm depends on SSC packages + (``did_imputation``/``reghdfe``/``ftools``/``require``) — install them once via + ``benchmarks/stata/requirements.do`` (the generators do not auto-install): .. code-block:: bash # macOS, StataSE 19 (binary not on PATH by default) - /Applications/Stata/StataSE.app/Contents/MacOS/stata-se -b do \ - benchmarks/stata/generate_lpdid_ra_golden.do + 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_imputation_loo_golden.do Running Benchmarks ~~~~~~~~~~~~~~~~~~ diff --git a/docs/doc-deps.yaml b/docs/doc-deps.yaml index 638e9a76..823fb87e 100644 --- a/docs/doc-deps.yaml +++ b/docs/doc-deps.yaml @@ -232,6 +232,8 @@ sources: type: methodology - path: docs/methodology/papers/borusyak-jaravel-spiess-2024-review.md type: methodology + - path: benchmarks/stata/README.md + type: methodology - path: docs/api/imputation.rst type: api_reference - path: docs/tutorials/11_imputation_did.ipynb diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index 2e27af78..36ff4481 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -1666,7 +1666,7 @@ where `W_it(h) = 1[K_it = h]` are lead indicators, estimated on `Omega_0` only. - **Bootstrap inference:** Uses multiplier bootstrap on the Theorem 3 influence function: `psi_i = sum_t v_it * epsilon_tilde_it`. Cluster-level psi sums are pre-computed for each aggregation target (overall, per-horizon, per-group), then perturbed with multiplier weights (Rademacher by default; configurable via `bootstrap_weights` parameter to use Mammen or Webb weights, matching CallawaySantAnna). This is a library extension (not in the paper) consistent with CallawaySantAnna/SunAbraham bootstrap patterns. - **Auxiliary residuals (Equation 8):** Implements the paper's *unit-clustered* Equation 8 aggregator, `tau_tilde_g = sum_i (sum_{t in G_g,i} v_it)(sum_{t in G_g,i} v_it * tau_hat_it) / sum_i (sum_{t in G_g,i} v_it)^2` (Borusyak-Jaravel-Spiess 2024, eq. 8, p. 3272; minimal-excess-variance derivation in Supplementary Appendix A.8): for each unit form the within-unit weight sum `a_{i,g}` and weighted-effect sum `b_{i,g}` over the unit's observations in group `g`, then combine across units. Groups partition `Omega_1` via `aux_partition` (default `"cohort_horizon"` = cohort × event-time; also `"cohort"` / `"horizon"`). Unimputable (NaN `tau_hat`) and off-target observations carry `v_it = 0` and are excluded from the aggregation — exact for finite `tau_hat` (a zero-weight row adds 0 to both `a` and `b`) and NaN-safe; a group with no contributing observations falls back to the unweighted group mean (a variance no-op, since `psi_g = sum_t v_it * eps_tilde_it = 0` there). - **Note (deviation from R):** R `didimputation::did_imputation` computes the auxiliary aggregator as `sum(v_it^2 * tau_hat_it) / sum(v_it^2)` grouped by cohort × event-time only (no partition control is exposed). At that partition each unit contributes at most one observation per group, so the paper's unit-clustered Equation 8 reduces exactly to `sum(v^2 * tau)/sum(v^2)` — i.e. diff-diff matches R at the default `aux_partition="cohort_horizon"` (pinned in `tests/test_methodology_imputation.py::TestImputationDiDParityR`). diff-diff additionally offers the coarser `aux_partition="cohort"` / `"horizon"` groupings (where a unit may contribute several observations to a group), which have no R analogue and are validated by hand-calculation. Both implement the same paper Equation 8; only the available partition granularity differs. The earlier observation-level mean `sum(v*tau)/sum(v)` (pre-3.5.x) coincided with this only under uniform within-group weights; it was corrected to the exact unit-clustered form during the ImputationDiD methodology validation. -- **Note (leave-one-out variance refinement, Supp. App. A.9):** the opt-in `leave_one_out=True` applies the Borusyak-Jaravel-Spiess (2024) Supplementary Appendix A.9 finite-sample refinement. The non-LOO `tau_tilde_g` (eq. 8) is built from the fitted `tau_hat_it`, which contain the noise `epsilon_it`, so it partially overfits and the auxiliary residuals `epsilon_tilde_it` are too small, biasing the variance **downward**. LOO recomputes each unit's group aggregate excluding that unit, `tau_tilde_it^LO = sum_{j!=i} v_jg^2 T_jg / sum_{j!=i} v_jg^2`, implemented efficiently (A.9) by rescaling each treated residual `epsilon_tilde_it^LO = epsilon_tilde_it / (1 - v_ig^2 / sum_j v_jg^2)` — where `v_ig = sum_{t in G_g,i} v_it` and `sum_j v_jg^2` are already materialized in `_compute_auxiliary_residuals_treated` as `a_{i,g}` and `per_group['den']`. That rescale reproduces the direct-LOO per-unit cluster sum `psi_i = sum_t v_it * epsilon_tilde_it` **exactly** (`psi_i^rescale = a_{i,g}(T_ig - tau_tilde_g) D/(D - a^2) = psi_i^direct-LOO`; a machine-precision-verified identity — **paper-fidelity** to the source `tau_tilde_it^LO`, not merely internal consistency). At the **default unit clustering** LOO gives a larger, less-downward-biased SE (Prop. A8: unbiased for an upper bound on the true variance; an equal-weight K-unit group inflates residuals by exactly `K/(K-1)`). **Default `leave_one_out=False`** preserves R `didimputation` parity (which omits LOO). **Edge (App. A.9 fn. 51):** a group with a single positive-weight unit has an undefined LOO (rescale factor `1/0`); such groups keep the non-LOO residual and the fit emits one consolidated `UserWarning` (a coarser `aux_partition` reduces singletons); a genuinely unit-dominated `>=2`-unit group keeps its large finite factor (intended inflation). **Composition scope:** the rescale operates on `epsilon_tilde`, so it flows through the coarser-`cluster=` CR sum, the analytical survey PSU-TSL variance, and the multiplier bootstrap unchanged — but Prop. A8's upper-bound guarantee and the `LOO >= non-LOO` direction hold at the default unit clustering only (under a coarser `cluster=`, per-unit `psi` inflation can partially cancel), so those compositions are a documented library extension, not paper-derived. **Replicate-weight survey designs** (BRR / Fay / JK1 / JKn / SDR) raise `NotImplementedError` with `leave_one_out=True`: replicate variance is computed by per-replicate point-estimate refits, bypassing the conservative-IF residual path where the LOO rescale lives, so LOO would silently no-op (fail-closed, no-silent-failures). **Effective-singleton guard:** a group's singleton test counts units with *positive* squared weight (not raw rows), and the leave-one-out denominator is the exact sum of the *other* units' squared weights (not `D - v_ig^2`, which can cancel to `<= 0` in float64 for an extremely dominated `>=2`-unit group and would silently revert it to non-LOO). **Reference / validation:** the authors' Stata `did_imputation` ships the same option; there is no CI-runnable anchor (R `didimputation` omits LOO, Stata is not in CI), so validation is the exact psi-identity + hand-calc + MC coverage (`tests/test_methodology_imputation.py::TestB2024AppendixA9LeaveOneOut`). Source: arXiv:2108.12419v5 App. A.9 (the REStud Supplementary Material is canonical); the main-article review's A.9 GAP is now filled (see the `borusyak-jaravel-spiess-2024-review.md` provenance note). The stronger Prop. A8 variant that *also* leaves out for the `delta_hat` covariate estimation (exact-unbiased upper bound) is noted but not implemented — the Stata option is the `tau_tilde` residual rescale only. +- **Note (leave-one-out variance refinement, Supp. App. A.9):** the opt-in `leave_one_out=True` applies the Borusyak-Jaravel-Spiess (2024) Supplementary Appendix A.9 finite-sample refinement. The non-LOO `tau_tilde_g` (eq. 8) is built from the fitted `tau_hat_it`, which contain the noise `epsilon_it`, so it partially overfits and the auxiliary residuals `epsilon_tilde_it` are too small, biasing the variance **downward**. LOO recomputes each unit's group aggregate excluding that unit, `tau_tilde_it^LO = sum_{j!=i} v_jg^2 T_jg / sum_{j!=i} v_jg^2`, implemented efficiently (A.9) by rescaling each treated residual `epsilon_tilde_it^LO = epsilon_tilde_it / (1 - v_ig^2 / sum_j v_jg^2)` — where `v_ig = sum_{t in G_g,i} v_it` and `sum_j v_jg^2` are already materialized in `_compute_auxiliary_residuals_treated` as `a_{i,g}` and `per_group['den']`. That rescale reproduces the direct-LOO per-unit cluster sum `psi_i = sum_t v_it * epsilon_tilde_it` **exactly** (`psi_i^rescale = a_{i,g}(T_ig - tau_tilde_g) D/(D - a^2) = psi_i^direct-LOO`; a machine-precision-verified identity — **paper-fidelity** to the source `tau_tilde_it^LO`, not merely internal consistency). At the **default unit clustering** LOO gives a larger, less-downward-biased SE (Prop. A8: unbiased for an upper bound on the true variance; an equal-weight K-unit group inflates residuals by exactly `K/(K-1)`). **Default `leave_one_out=False`** preserves R `didimputation` parity (which omits LOO). **Edge (App. A.9 fn. 51):** a group with a single positive-weight unit has an undefined LOO (rescale factor `1/0`); such groups keep the non-LOO residual and the fit emits one consolidated `UserWarning` (a coarser `aux_partition` reduces singletons); a genuinely unit-dominated `>=2`-unit group keeps its large finite factor (intended inflation). **Composition scope:** the rescale operates on `epsilon_tilde`, so it flows through the coarser-`cluster=` CR sum, the analytical survey PSU-TSL variance, and the multiplier bootstrap unchanged — but Prop. A8's upper-bound guarantee and the `LOO >= non-LOO` direction hold at the default unit clustering only (under a coarser `cluster=`, per-unit `psi` inflation can partially cancel), so those compositions are a documented library extension, not paper-derived. **Replicate-weight survey designs** (BRR / Fay / JK1 / JKn / SDR) raise `NotImplementedError` with `leave_one_out=True`: replicate variance is computed by per-replicate point-estimate refits, bypassing the conservative-IF residual path where the LOO rescale lives, so LOO would silently no-op (fail-closed, no-silent-failures). **Effective-singleton guard:** a group's singleton test counts units with *positive* squared weight (not raw rows), and the leave-one-out denominator is the exact sum of the *other* units' squared weights (not `D - v_ig^2`, which can cancel to `<= 0` in float64 for an extremely dominated `>=2`-unit group and would silently revert it to non-LOO). **Reference / validation:** the authors' Stata `did_imputation` ships the same option (`leaveout`), and the library LOO SE is now **measured against it** — the LOO SE matches to ~1e-9 (the point to ~2e-8; cross-implementation: `did_imputation` goes through `reghdfe`, the library through its own sparse IF solver; the parity test gates at `abs=1e-7`) across the overall ATT and all event-study horizons on the committed panel (`benchmarks/stata/generate_imputation_loo_golden.do` → `benchmarks/data/didimputation_loo_stata_golden.json`; `tests/test_imputation_loo_stata_parity.py`, which also three-way-confirms the non-LOO SE against R). The partition is pinned explicitly `avgeffectsby(Ei t)` (== `aux_partition="cohort_horizon"`), which is also `did_imputation`'s current default — pinned to keep the validation estimand self-describing and robust to a future default change. Stata is node-locked (not in CI, golden committed, exactly like the R arm), so the CI-runnable guards remain the exact psi-identity + hand-calc + MC coverage (`tests/test_methodology_imputation.py::TestB2024AppendixA9LeaveOneOut`). Source: arXiv:2108.12419v5 App. A.9 (the REStud Supplementary Material is canonical); the main-article review's A.9 GAP is now filled (see the `borusyak-jaravel-spiess-2024-review.md` provenance note). The stronger Prop. A8 variant that *also* leaves out for the `delta_hat` covariate estimation (exact-unbiased upper bound) is noted but not implemented — the Stata option is the `tau_tilde` residual rescale only. - **Note:** The Step-1 iterative FE solver (`_iterative_fe`) routes through the shared bincount Gauss-Seidel helper `diff_diff.utils._iterative_fe_solve`, and the covariate/pre-trend within-transformation routes through the shared MAP engine `diff_diff.utils.demean_by_groups` (factorize-once + `np.bincount`, optional Rust kernel; group order `[time, unit]` preserving the historical time-then-unit sweep) — the same convergence contract, accumulation-order numerics (~1e-10 vs the pre-3.7 pandas loops, not bit-for-bit), and `max_iter=10_000` budget documented under "Absorbed Fixed Effects with Survey Weights". Both surfaces emit `UserWarning` via `diff_diff.utils.warn_if_not_converged` when `max_iter` exhausts without reaching `tol` (the demean warning now carries the shared-engine label naming the affected variables rather than the estimator name). Silent return of the current iterate was classified as a silent failure under the Phase 2 audit and replaced with an explicit signal to match the logistic/Poisson IRLS pattern in `linalg.py`. - **Note:** Zero-total-weight groups (e.g. whole PSUs zeroed by JK1/BRR replicate weights, which reach Step 1 unmasked — `keep_mask` only drops always-treated units): a unit/period whose observations ALL carry zero weight has no identifying contribution and surfaces as `NaN` FE (key retained for the rank-condition membership check; matches the SpilloverDiD `_iterative_fe_subset` REGISTRY contract — never a silent finite `0.0`), and the shared demean engine's inert-row guard leaves those rows un-demeaned instead of NaN-poisoning the column. Before 3.7 the pandas loops divided 0/0 there: the covariate replicate path NaN-poisoned `y_dm`/`X_dm`, failed EVERY replicate refit inside `solve_ols(check_finite=True)`, and returned NaN SEs after a non-convergence warning storm; a main fit with zero-weight rows + covariates raised the same opaque `ValueError`. Both now produce finite results. `TwoStageDiD._mask_nan_ytilde`'s "non-finite imputed outcomes" `UserWarning` is suppressed (via `warn_nan=False`) ONLY inside the replicate-refit closures, where NaN FE for zeroed PSUs is expected mechanics — the main-fit warning is unchanged. - **Note:** `vcov_type` is permanently narrow to `{"hc1"}` per the Theorem 3 IF-based variance decomposition. Analytical-sandwich families `{classical, hc2, hc2_bm}` are rejected at `__init__` — the per-unit influence function aggregation has no equivalent single design matrix on which hat-matrix leverage or Bell-McCaffrey Satterthwaite DOF can be defined. `cluster=` invokes per-cluster IF summation (Theorem 3 equation 7 conservative variance, `sigma_sq = (cluster_psi_sums**2).sum()` — plain CR1 with no Stata-style `(n-1)/(n-p)` finite-sample factor because the IF has no design-matrix `p` in the OLS sense); `cluster=None` (the default) routes the SAME Theorem 3 cluster-summed IF variance with `cluster_var = unit` (the unit column passed to `fit()`), so the summary renders `"CR1 cluster-robust at , G="` rather than the generic `"HC1"` label; `survey_design=` invokes TSL on the combined IF. Under bootstrap (`n_bootstrap > 0`) the analytical variance-family label is suppressed in `summary()` because `fit()` overwrites the reported SE/CI/p-value with bootstrap_results (mirrors the canonical `DiDResults` gate at `results.py:213-226`). `vcov_type='conley'` is deferred to the ImputationDiD Conley follow-up row in DEFERRED.md. diff --git a/docs/methodology/papers/borusyak-jaravel-spiess-2024-review.md b/docs/methodology/papers/borusyak-jaravel-spiess-2024-review.md index 19ad8b7a..4395e65b 100644 --- a/docs/methodology/papers/borusyak-jaravel-spiess-2024-review.md +++ b/docs/methodology/papers/borusyak-jaravel-spiess-2024-review.md @@ -189,4 +189,4 @@ No regularization, bandwidth, factor-count, or cross-validation tuning (the esti - No contradictions were found between the three extraction passes or against the direct page reads; all numbered equations rendered cleanly (no `[UNREADABLE EQUATION]`). - **Resolved in PR-B:** the auxiliary `τ̃_g` aggregator now implements the exact unit-clustered Eq. (8) (`_compute_auxiliary_residuals_treated` + docstring + REGISTRY label), validated by white-box hand-calc + R `didimputation` parity (`tests/test_methodology_imputation.py`). The earlier obs-level mean was a valid conservative simplification but not the variance-minimizing form. - **Done in PR-B:** the **R `didimputation` parity fixture** and **`tests/test_methodology_imputation.py`** are on file. The `v*_it` weights (Supp. Prop. A3, absent from this PDF) were validated *empirically* against the reference: the exact two-way-FE projection `-A₀(A₀'A₀)⁻¹A₁'w` now matches `didimputation` (observed ~1e-10; tests assert SE `abs=1e-7`) — the prior FE-only closed form `-(w_i/n0_i + w_t/n0_t − w/N₀)` was a balanced-panel approximation that biased the SE ~27% on staggered (unbalanced-Ω₀) designs and was corrected in PR-B. The leave-one-out variance refinement (Supp. App. A.9) is now **implemented** as the opt-in `leave_one_out` parameter (default False, preserving R `didimputation` parity) — see the A.9 provenance note below. -- **Supp. App. A.9 provenance (leave-one-out, added post-review):** A.9 was a GAP in this main-article review (the REStud Supplementary Material is not in the reviewed PDF). It was subsequently sourced from the arXiv preprint supplement (**arXiv:2108.12419v5, Appendix A.9 "Leave-Out Conservative Variance Estimation"**) to implement the opt-in `leave_one_out` parameter. The **REStud Supplementary Material is the canonical version**; per this doc's version-pin scope note, the efficient-rescale formula `ε̃^LO = ε̃/(1 − v_ig²/Σ_j v_jg²)`, its exact ψ-level equivalence to the direct leave-one-out `τ̃_it^LO`, Prop. A8 (unbiased for an upper bound), and the footnote-51 single-positive-weight-unit edge are recorded in `docs/methodology/REGISTRY.md` `## ImputationDiD` (NOT transcribed here). Validated by the exact ψ-identity + hand-calc + MC coverage (`tests/test_methodology_imputation.py::TestB2024AppendixA9LeaveOneOut`); the authors' Stata `did_imputation` ships the same option (no CI-runnable anchor). +- **Supp. App. A.9 provenance (leave-one-out, added post-review):** A.9 was a GAP in this main-article review (the REStud Supplementary Material is not in the reviewed PDF). It was subsequently sourced from the arXiv preprint supplement (**arXiv:2108.12419v5, Appendix A.9 "Leave-Out Conservative Variance Estimation"**) to implement the opt-in `leave_one_out` parameter. The **REStud Supplementary Material is the canonical version**; per this doc's version-pin scope note, the efficient-rescale formula `ε̃^LO = ε̃/(1 − v_ig²/Σ_j v_jg²)`, its exact ψ-level equivalence to the direct leave-one-out `τ̃_it^LO`, Prop. A8 (unbiased for an upper bound), and the footnote-51 single-positive-weight-unit edge are recorded in `docs/methodology/REGISTRY.md` `## ImputationDiD` (NOT transcribed here). Validated by the exact ψ-identity + hand-calc + MC coverage (`tests/test_methodology_imputation.py::TestB2024AppendixA9LeaveOneOut`); the authors' Stata `did_imputation` ships the same option (`leaveout`), now a committed Stata parity anchor matching the library LOO SE to ~1e-9 at the overall ATT and all event-study horizons (`benchmarks/stata/generate_imputation_loo_golden.do`; `tests/test_imputation_loo_stata_parity.py`). diff --git a/tests/test_imputation_loo_stata_parity.py b/tests/test_imputation_loo_stata_parity.py new file mode 100644 index 00000000..bf9a1f47 --- /dev/null +++ b/tests/test_imputation_loo_stata_parity.py @@ -0,0 +1,217 @@ +"""External-reference parity: ImputationDiD leave-one-out SE vs Stata `did_imputation`. + +`ImputationDiD(leave_one_out=True)` applies the Borusyak-Jaravel-Spiess (2024) +Supplementary Appendix A.9 finite-sample variance refinement. No R package computes it +(R `didimputation` omits LOO), so the library LOO SE was validated only by an internal +psi-identity + hand-calc + MC coverage. The authors' own Stata `did_imputation` ships +the same option (`leaveout`); this arm turns that into a measured anchor. + +Second Stata parity arm, and the first SSC-dependent one. The golden is produced by +`benchmarks/stata/generate_imputation_loo_golden.do`, which runs `did_imputation ..., +leaveout avgeffectsby(Ei t) cluster(unit)` on the committed R-arm panel +`didimputation_test_panel.csv` (no clean-sample reconstruction - did_imputation takes +the raw panel; the only mapping is Ei = first_treat, missing for never-treated). + +The library uses its own sparse IF solver while did_imputation goes through reghdfe, so +agreement is cross-implementation, not bit-identical: the SE agrees to ~1e-9 and the +point to ~2e-8. Gates use ``abs=1e-7, rel=0`` - the repo's imputation SE-parity tolerance +(`test_methodology_imputation.py::test_overall_se_matches_r` asserts 1e-7 despite +observing ~1e-10) - since the golden is committed and this test recomputes the library +on cross-platform CI. + +Gates (overall + each event-study horizon): +1. **Point** - library ATET vs Stata ATET (same estimand; overall ~6.5e-11, per-horizon + up to ~1.8e-8 locally - both well inside the abs=1e-7 gate). +2. **LOO SE anchor** - library LOO SE vs Stata `leaveout` SE (the point of this arm). +3. **Non-LOO SE** - library non-LOO SE vs Stata non-`leaveout` SE (corroboration). +4. **Non-LOO three-way** - Stata golden `se_nonloo` vs the committed R golden `se` + (library == R == Stata on the same panel; committed-vs-committed). +5. **Warning-cleanliness** - the LOO fit emits no A.9-singleton `UserWarning`, proving + the LOO rescale is genuinely exercised rather than silently falling back to non-LOO. + +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 rather +than fail. Regenerate (after `benchmarks/stata/requirements.do`) with:: + + /Applications/Stata/StataSE.app/Contents/MacOS/stata-se -b do \ + benchmarks/stata/generate_imputation_loo_golden.do +""" + +from __future__ import annotations + +import json +import warnings +from pathlib import Path + +import pandas as pd +import pytest + +from diff_diff import ImputationDiD + +_DATA = Path(__file__).parent.parent / "benchmarks" / "data" +STATA_GOLDEN_PATH = _DATA / "didimputation_loo_stata_golden.json" +R_GOLDEN_PATH = _DATA / "didimputation_golden.json" +PANEL_PATH = _DATA / "didimputation_test_panel.csv" + +_FIXTURE_AVAILABLE = ( + STATA_GOLDEN_PATH.is_file() and R_GOLDEN_PATH.is_file() and PANEL_PATH.is_file() +) + +HORIZONS = [0, 1, 2, 3, 4, 5] + +# Cross-platform tolerance: the library IF variance vs did_imputation (reghdfe) agree +# ~1e-9..1e-10 same-machine; 1e-7 matches the repo's imputation SE-parity convention +# and absorbs cross-platform BLAS variation. rel=0 -> purely absolute bound. +ATOL = 1e-7 + + +def _skip_if_missing() -> None: + if not _FIXTURE_AVAILABLE: + pytest.skip( + "Stata ImputationDiD-LOO parity fixture not present. Regenerate via " + "`stata-se -b do benchmarks/stata/generate_imputation_loo_golden.do` " + "(after benchmarks/stata/requirements.do)." + ) + + +@pytest.fixture(scope="module") +def stata_golden() -> dict: + _skip_if_missing() + with STATA_GOLDEN_PATH.open() as f: + return json.load(f) + + +@pytest.fixture(scope="module") +def r_golden() -> dict: + _skip_if_missing() + with R_GOLDEN_PATH.open() as f: + return json.load(f) + + +@pytest.fixture(scope="module") +def library_fit() -> dict: + """Fit the LOO and non-LOO paths once each; return per-horizon + overall + warnings. + + ImputationDiD exposes ``res.event_study_effects[h]["effect"]`` / ``["se"]`` (a + dict-of-dicts, NOT a ``.event_study`` DataFrame) and scalar ``res.overall_att`` / + ``res.overall_se``. ``aggregate="event_study"`` also computes the overall + unconditionally, so one fit per variance mode suffices. + """ + _skip_if_missing() + panel = pd.read_csv(PANEL_PATH) + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + loo = ImputationDiD(leave_one_out=True).fit( + panel, + outcome="y", + unit="unit", + time="time", + first_treat="first_treat", + aggregate="event_study", + ) + loo_msgs = [str(w.message) for w in caught] + nonloo = ImputationDiD(leave_one_out=False).fit( + panel, + outcome="y", + unit="unit", + time="time", + first_treat="first_treat", + aggregate="event_study", + ) + + def by_h(res): + return {h: res.event_study_effects[h] for h in HORIZONS} + + return { + "loo_overall": (float(loo.overall_att), float(loo.overall_se)), + "loo_es": {h: (float(v["effect"]), float(v["se"])) for h, v in by_h(loo).items()}, + "nonloo_overall_se": float(nonloo.overall_se), + "nonloo_es_se": {h: float(v["se"]) for h, v in by_h(nonloo).items()}, + "loo_warnings": loo_msgs, + } + + +# ----- Gate 1: point (library vs Stata, same estimand) ----- + + +def test_overall_point_matches_stata(stata_golden, library_fit): + lib_att, _ = library_fit["loo_overall"] + st_att = stata_golden["overall"]["att"] + assert lib_att == pytest.approx(st_att, abs=ATOL, rel=0) + + +@pytest.mark.parametrize("h", HORIZONS) +def test_event_study_point_matches_stata(h, stata_golden, library_fit): + lib_att, _ = library_fit["loo_es"][h] + st_att = stata_golden["event_study"][str(h)]["att"] + assert lib_att == pytest.approx(st_att, abs=ATOL, rel=0) + + +# ----- Gate 2: the LOO SE anchor (library vs Stata `leaveout`) ----- + + +def test_overall_loo_se_matches_stata(stata_golden, library_fit): + _, lib_se = library_fit["loo_overall"] + st_se = stata_golden["overall"]["se"] + assert lib_se == pytest.approx( + st_se, abs=ATOL, rel=0 + ), f"overall LOO SE: library {lib_se} != Stata leaveout {st_se}" + + +@pytest.mark.parametrize("h", HORIZONS) +def test_event_study_loo_se_matches_stata(h, stata_golden, library_fit): + _, lib_se = library_fit["loo_es"][h] + st_se = stata_golden["event_study"][str(h)]["se"] + assert lib_se == pytest.approx( + st_se, abs=ATOL, rel=0 + ), f"h={h} LOO SE: library {lib_se} != Stata leaveout {st_se}" + + +# ----- Gate 3: non-LOO SE corroboration (library vs Stata non-`leaveout`) ----- + + +def test_overall_nonloo_se_matches_stata(stata_golden, library_fit): + lib_se = library_fit["nonloo_overall_se"] + st_se = stata_golden["overall"]["se_nonloo"] + assert lib_se == pytest.approx(st_se, abs=ATOL, rel=0) + + +@pytest.mark.parametrize("h", HORIZONS) +def test_event_study_nonloo_se_matches_stata(h, stata_golden, library_fit): + lib_se = library_fit["nonloo_es_se"][h] + st_se = stata_golden["event_study"][str(h)]["se_nonloo"] + assert lib_se == pytest.approx(st_se, abs=ATOL, rel=0) + + +# ----- Gate 4: non-LOO three-way (Stata golden vs committed R golden) ----- + + +def test_overall_nonloo_se_three_way_stata_vs_r(stata_golden, r_golden): + """Committed-vs-committed: Stata non-LOO SE == R didimputation SE (same panel).""" + st_se = stata_golden["overall"]["se_nonloo"] + r_se = r_golden["overall"]["se"] + assert st_se == pytest.approx(r_se, abs=ATOL, rel=0) + + +@pytest.mark.parametrize("h", HORIZONS) +def test_event_study_nonloo_se_three_way_stata_vs_r(h, stata_golden, r_golden): + """Committed-vs-committed, per horizon: Stata non-LOO SE == R didimputation SE.""" + st_se = stata_golden["event_study"][str(h)]["se_nonloo"] + r_by_h = dict(zip(r_golden["event_study"]["horizons"], r_golden["event_study"]["se"])) + assert st_se == pytest.approx(r_by_h[h], abs=ATOL, rel=0) + + +# ----- Gate 5: the LOO rescale is genuinely exercised (no singleton fallback) ----- + + +def test_loo_fit_emits_no_singleton_warning(library_fit): + """App. A.9 fn. 51: a single-positive-weight group has an undefined LOO and falls + back with a UserWarning. On this panel every group has >=2 units, so the LOO fit + must be warning-clean - otherwise the LOO SE would silently equal the non-LOO SE. + """ + singletons = [ + m + for m in library_fit["loo_warnings"] + if "leave-one-out" in m.lower() or "singleton" in m.lower() or "single" in m.lower() + ] + assert not singletons, f"unexpected LOO-singleton warning(s): {singletons}"