diff --git a/.gitignore b/.gitignore index 6a62748f2..c4dc58890 100644 --- a/.gitignore +++ b/.gitignore @@ -120,3 +120,8 @@ tools/reviewer-eval/runs/ benchmarks/refresh_2026_07/venvs/ 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 diff --git a/CHANGELOG.md b/CHANGELOG.md index 1a13b6d16..573289d0f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,18 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +### Added +- **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 + pinned against itself + MC-coverage-validated. It is now cross-validated against Stata + `teffects ra ... atet vce(cluster)` via a committed golden + (`benchmarks/stata/generate_lpdid_ra_golden.do` → `lpdid_ra_stata_golden.json`; + `tests/test_lpdid_ra_stata_parity.py`): the library matches `teffects` to ~1e-16 at all + 7 event-study horizons, confirming the no-finite-sample-factor convention (REGISTRY + LPDiD Deviation 2). This is the repo's first Stata parity arm — `teffects` is native + (no SSC dependency) and goldens are committed, so CI needs no Stata. + ### Changed - **Internal tracking docs reorganized (no library behavior change).** `TODO.md` is now the actionable backlog only; blocked / parked work and won't-fix decisions moved to the diff --git a/DEFERRED.md b/DEFERRED.md index 9146a0331..d0cbf57ac 100644 --- a/DEFERRED.md +++ b/DEFERRED.md @@ -57,7 +57,7 @@ exists but parity can't be verified without a local toolchain. | `CallawaySantAnna` bootstrap: align p-value computation with R `did`'s symmetric-percentile method (former "CallawaySantAnna Bootstrap Improvements" section). | `staggered.py` | — | Low | | **`bias_corrected_local_linear` (lprobust) Phase-1c follow-ups:** extend golden parity to `kernel ∈ {triangular, uniform}` (epa-only today); expose `vce ∈ {hc0,hc1,hc2,hc3}` on the public wrapper once R goldens exist (port supports all four; needs a per-mode generator + a hc2/hc3 q-fit-leverage decision); clustered-DGP auto-bandwidth parity is **blocked upstream** on an nprobust singleton-cluster bug in `lpbwselect.mse.dpi` (Phase-1c DGP 4 uses manual `h=b=0.3`). | `_nprobust_port.py`, `local_linear.py`, `generate_nprobust_lprobust_golden.R` | Phase 1c | Low-Med | | `HeterogeneousAdoptionDiD` Stute-family Stata-bridge parity: no public R `Stutetest` package exists; would add `benchmarks/stata/generate_stute_golden.do` + a Stata dependency. | `benchmarks/stata/`, `tests/test_stute_test_parity.py` | follow-up | Low | -| **`LPDiD` regression-adjustment SE — no runnable R reference.** The RA influence-function cluster SE is canonically Stata `teffects ra ... atet vce(cluster)` only; no R package computes it (`alexCardazzi/lpdid` does direct covariate inclusion, not RA). Today the RA *point* is R-anchored (~1e-12), the SE is pinned + MC-coverage-validated (`coverage_lpdid_ra.py`). Follow-up: contribute the RA path to `alexCardazzi/lpdid` so a runnable R RA reference exists — only a *trusted* anchor once cross-checked vs Stata `teffects` (else circular). | `tests/test_methodology_lpdid.py`, `benchmarks/python/coverage_lpdid_ra.py` | #B2 follow-up | Low | +| **`LPDiD` 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 | | `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 | diff --git a/benchmarks/README.md b/benchmarks/README.md index b03d270d2..3a5428bee 100644 --- a/benchmarks/README.md +++ b/benchmarks/README.md @@ -43,6 +43,21 @@ Required R packages: - `jsonlite` - JSON interchange - `data.table` - Fast data manipulation +### Stata + +A few goldens are anchored against Stata where no runnable R reference exists. +Stata is node-locked single-user, so — like the R arm — the goldens are committed +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 +``` + +Current arm uses only **native** Stata commands (`teffects`) — no SSC packages — +so `version 19` pins behavior. See `stata/README.md`. + ## Directory Structure ``` @@ -56,6 +71,9 @@ benchmarks/ │ ├── benchmark_synthdid.R # Synthetic DiD │ ├── 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 ├── python/ │ ├── utils.py # Common utilities │ ├── benchmark_callaway.py # CallawaySantAnna @@ -72,11 +90,12 @@ benchmarks/ ## Estimator Comparisons -| diff-diff | R Package | Reference | Status | +| diff-diff | Reference Package | Reference | Status | |-----------|-----------|-----------|--------| | `CallawaySantAnna` | `did::att_gt` | Callaway & Sant'Anna (2021) | ✓ Integrated | | `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 | | `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/lpdid_ra_stata_golden.json b/benchmarks/data/lpdid_ra_stata_golden.json new file mode 100644 index 000000000..991e74096 --- /dev/null +++ b/benchmarks/data/lpdid_ra_stata_golden.json @@ -0,0 +1,24 @@ +{ + "meta": { + "estimator": "LPDiD regression-adjustment (RA) SE - Stata teffects ra atet", + "generator": "benchmarks/stata/generate_lpdid_ra_golden.do", + "source_panel": "benchmarks/data/lpdid_test_panel.csv", + "point_anchor": "benchmarks/data/lpdid_golden.json ra_cov[h][0]", + "stata_edition": "SE", + "cmd": "teffects ra (Dy x i.time) (tdiff), atet vce(cluster unit)", + "se_convention": "cluster-robust IF at unit, NO finite-sample factor (teffects convention); t(G-1) reference on the library side", + "note": "Independent Stata reconstruction of the per-horizon clean sample (ports generate_lpdid_golden.R prep/clean_h). N and G are emitted to power the Python row-count sample-integrity gate.", + "stata_version": 19.0, + "pre_window": 3, + "post_window": 4 + }, + "ra_se": { + "0": {"att": 2.888871838355342, "se": 0.45516589747218866, "N": 522, "G": 60}, + "1": {"att": 3.225268910679945, "se": 0.44522904812910252, "N": 461, "G": 60}, + "2": {"att": 3.5779502216280248, "se": 0.46179396866595729, "N": 401, "G": 60}, + "3": {"att": 3.973278128236954, "se": 0.39295354484873846, "N": 353, "G": 60}, + "4": {"att": 4.2651037814454096, "se": 0.39398602601018534, "N": 305, "G": 60}, + "-2": {"att": 0.1204444745695894, "se": 0.43592937799131937, "N": 461, "G": 60}, + "-3": {"att": -0.11970460945194859, "se": 0.43296260860283592, "N": 401, "G": 60} + } +} diff --git a/benchmarks/stata/README.md b/benchmarks/stata/README.md new file mode 100644 index 000000000..b9bc6891e --- /dev/null +++ b/benchmarks/stata/README.md @@ -0,0 +1,133 @@ +# Stata parity benchmarks + +Stata golden generators live here. They are the **first Stata arm** in the repo; +the pattern mirrors `benchmarks/R/` (a `generate_*` script writes a committed +golden JSON that a skip-guarded `tests/test_*_parity.py` reads, so CI never needs +Stata). Stata is node-locked single-user, so — exactly like the R arm — goldens +are committed and only regenerated locally. + +Locating the binary (macOS, StataSE 19; **not** on `PATH` by default): + +``` +/Applications/Stata/StataSE.app/Contents/MacOS/stata-se +``` + +--- + +# `teffects ra` parity for the LP-DiD regression-adjustment SE + +`benchmarks/stata/generate_lpdid_ra_golden.do` produces +`benchmarks/data/lpdid_ra_stata_golden.json`, consumed by +`tests/test_lpdid_ra_stata_parity.py` to verify that diff-diff's LP-DiD +regression-adjustment (RA) standard error (Dube, Girardi, Jorda & Taylor 2025) +matches Stata `teffects ra ... atet vce(cluster)` across the 7 event-study +horizons. + +## Why Stata `teffects` + +The RA covariate path reports an influence-function cluster variance with **no +finite-sample factor**. **No R package computes it** — `alexCardazzi/lpdid` uses +direct covariate inclusion, not RA — so there is no runnable R analogue. The +canonical reference is Stata `teffects ra ... atet vce(cluster)`, which the LP-DiD +authors themselves invoke (paper footnote 9). `teffects` is a **native** Stata +command (no SSC package), so `version 19` in the generator fully pins its +numerical behavior; there is no third-party package version to track. + +This arm converts REGISTRY `## LPDiD` Deviation 2 from an *inference* (the +no-finite-sample-factor convention was argued from degrees-of-freedom comments in +the authors' `.do` files) into a *measurement*: the library RA IF SE matches +`teffects` to ~1e-16 at every horizon (same machine / BLAS). + +## What the generator does + +It reads (does **not** regenerate) `benchmarks/data/lpdid_test_panel.csv` — whose +sole owner is `benchmarks/R/generate_lpdid_golden.R` — and **independently +reconstructs** each horizon's clean sample by porting the R `prep` + `clean_h` +recipe (`generate_lpdid_golden.R:97-123`) to `tsset` + `L.`/`F.` time-series +operators. This gives three independent sample constructions (Python / R / Stata); +a port bug is very likely to be caught before it could masquerade as an SE finding, +because the Python test's point gate (Stata ATET vs the R-anchored `ra_cov[h][0]`) and +sample-shape gate (`(e(N), e(N_clust))` vs the library's `(n_obs, n_clusters)`) would +fail first. + +`tsfill` is load-bearing: it mirrors R `fill_gaps()` so the absorbing `treat` +recompute and the long differences see the completed calendar grid across the +deliberate interior gap at `(unit==60, time==7)`. + +The per-horizon command is: + +```stata +teffects ra (Dy x i.time) (tdiff), atet vce(cluster unit) +``` + +(The authors' footnote-9 syntax omits the covariate; the RA *covariate* path adds +`x` to the outcome model, as here.) + +## Regenerating + +```bash +# from the repo root +/Applications/Stata/StataSE.app/Contents/MacOS/stata-se -b do \ + benchmarks/stata/generate_lpdid_ra_golden.do + +# batch mode ALWAYS exits 0 — verify the log has no Stata errors: +grep -E '^r\([0-9]+\);' generate_lpdid_ra_golden.log # must print nothing +``` + +The output JSON is **committed** so CI doesn't need Stata. Only re-run when the +committed panel changes (then also refresh `RA_SE_PIN` in +`tests/test_methodology_lpdid.py`) or the horizon set changes. + +## Skip behavior + +`tests/test_lpdid_ra_stata_parity.py` calls `pytest.skip(...)` when the golden (or +the shared panel / R golden it cross-checks against) is absent, so CI passes +without Stata. The estimator itself is exercised independently by +`tests/test_methodology_lpdid.py` and `tests/test_lpdid.py`. + +## JSON schema + +```json +{ + "meta": { + "estimator": "LPDiD regression-adjustment (RA) SE - Stata teffects ra atet", + "generator": "benchmarks/stata/generate_lpdid_ra_golden.do", + "source_panel": "benchmarks/data/lpdid_test_panel.csv", + "point_anchor": "benchmarks/data/lpdid_golden.json ra_cov[h][0]", + "stata_edition": "SE", + "cmd": "teffects ra (Dy x i.time) (tdiff), atet vce(cluster unit)", + "se_convention": "...", + "note": "...", + "stata_version": 19.0, + "pre_window": 3, + "post_window": 4 + }, + "ra_se": { + "0": {"att": , "se": , "N": , "G": }, + "...": "post {0..4} then pre {-2,-3}; h=-1 is the omitted reference" + } +} +``` + +Values are written at `%21.17g` (round-trip-exact IEEE-754 double). `N` and `G` +(`e(N)` / `e(N_clust)`) are emitted so the Python gate can assert the Stata clean +sample matches the library's realized size and cluster count — which, together with +the point (1e-10) and SE (~1e-16) agreement, strongly corroborates the same clean +sample. No timestamp — the golden regenerates byte-identically. + +## 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). +- **`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. diff --git a/benchmarks/stata/generate_lpdid_ra_golden.do b/benchmarks/stata/generate_lpdid_ra_golden.do new file mode 100644 index 000000000..4aa5e2e07 --- /dev/null +++ b/benchmarks/stata/generate_lpdid_ra_golden.do @@ -0,0 +1,195 @@ +*! Golden generator: LPDiD regression-adjustment (RA) SE vs Stata `teffects ra`. +*! +*! Purpose +*! Produce an INDEPENDENT external anchor for the LP-DiD regression-adjustment +*! standard error (Dube, Girardi, Jorda & Taylor 2025). The canonical RA SE is +*! Stata `teffects ra ... atet vce(cluster)` only - no R package computes it +*! (`alexCardazzi/lpdid` does direct covariate inclusion, not RA). Until now the +*! library RA IF-cluster SE was pinned against itself (`RA_SE_PIN` in +*! tests/test_methodology_lpdid.py) and calibration-validated only by a +*! Monte-Carlo coverage study. This generator converts that self-pin into a +*! measured cross-implementation anchor. +*! +*! This is the FIRST Stata arm in the repo. `teffects` is NATIVE to Stata (no +*! SSC dependency), so `version 19` fully pins the numerical behavior. +*! +*! Consuming test +*! tests/test_lpdid_ra_stata_parity.py +*! +*! Outputs (checked into the repo) +*! benchmarks/data/lpdid_ra_stata_golden.json +*! +*! Usage (run from the repo root) +*! /Applications/Stata/StataSE.app/Contents/MacOS/stata-se -b do \ +*! benchmarks/stata/generate_lpdid_ra_golden.do +*! Then confirm the log is clean: grep -E '^r\([0-9]+\);' generate_lpdid_ra_golden.log +*! (Stata batch mode ALWAYS exits 0, even on error - never trust the exit code.) +*! +*! Notes +*! - Reads (does NOT regenerate) benchmarks/data/lpdid_test_panel.csv, whose sole +*! owner is benchmarks/R/generate_lpdid_golden.R. 60 units x 12 periods, one +*! deliberate interior gap at (unit==60, time==7) on a never-treated unit. +*! - Independently reconstructs the per-horizon clean sample by porting the R +*! `prep` + `clean_h` recipe (generate_lpdid_golden.R:97-123). `tsfill` is +*! load-bearing: it mirrors R `fill_gaps()` so the absorbing `treat` recompute +*! and the F./L. long differences see the completed calendar grid (a lead +*! across the unit-60 gap must yield a real value / a missing outcome, not a +*! spurious neighbor). +*! - The in-.do point gate (1e-8, informational) aborts early on a gross port +*! bug; the AUTHORITATIVE point gate is the Python test asserting Stata ATET +*! vs benchmarks/data/lpdid_golden.json `ra_cov[h][0]` at 1e-10. + +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". 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 + +* ------------------------------------------------------------------------------ +* R-anchored RA POINT estimates (benchmarks/data/lpdid_golden.json `ra_cov[h][0]`, +* R generator digits=12). Hard-coded here only for the in-.do smoke gate; the +* Python test reads the JSON directly and is the authoritative check. +* ------------------------------------------------------------------------------ +scalar rpt_0 = 2.888871838355 +scalar rpt_1 = 3.22526891068 +scalar rpt_2 = 3.577950221628 +scalar rpt_3 = 3.973278128237 +scalar rpt_4 = 4.265103781445 +scalar rpt_m2 = 0.1204444745696 +scalar rpt_m3 = -0.1197046094519 + +* ------------------------------------------------------------------------------ +* Load the committed panel and confirm the expected schema/types. +* ------------------------------------------------------------------------------ +import delimited using "benchmarks/data/lpdid_test_panel.csv", clear varnames(1) +confirm numeric variable unit time treat y x +quietly count +assert r(N) == 719 // 60*12 - 1 interior gap + +* ------------------------------------------------------------------------------ +* prep: port of generate_lpdid_golden.R:97-110. +* - tsfill materializes the interior-gap row (mirrors R fill_gaps()). +* - treat is recomputed as an ABSORBING fill on the completed grid. +* ------------------------------------------------------------------------------ +tsset unit time +tsfill // <-- mirrors fill_gaps() + +* treat_date = first period the imported treat==1, per unit (never-treated -> .). +* tsfill rows carry treat==. so cond(treat==1,.)-> . and egen min ignores them. +bysort unit (time): egen treat_date = min(cond(treat==1, time, .)) + +replace treat = (!missing(treat_date) & time >= treat_date) // absorbing on grid +gen double Ly = L.y // base for long diff +gen tdiff = treat - L.treat +replace tdiff = 0 if missing(tdiff) | tdiff < 0 +gen byte obs = !missing(y) + +* ------------------------------------------------------------------------------ +* Per-horizon: rebuild the clean sample (clean_h) inside preserve/restore, run +* teffects ra ... atet vce(cluster unit), capture ATET/SE/N/G into scalars. +* Post h in {0,1,2,3,4}; pre h in {2,3} (key = -h; h=-1 is the omitted reference). +* ------------------------------------------------------------------------------ +tempname fh +foreach h of numlist 0 1 2 3 4 { + preserve + * clean_h(post): Dy = lead(y,h) - Ly ; Fh = lead(treat,h) + gen double Dy = F`h'.y - Ly + gen Fh = F`h'.treat + keep if obs & !missing(Dy) & !missing(tdiff) & !missing(Fh) & (tdiff==1 | Fh==0) + teffects ra (Dy x i.time) (tdiff), atet vce(cluster unit) + scalar att_`h' = _b[r1vs0.tdiff] + scalar se_`h' = _se[r1vs0.tdiff] + scalar N_`h' = e(N) + scalar G_`h' = e(N_clust) + * in-.do smoke gate (informational; Python test is authoritative) + assert reldif(att_`h', rpt_`h') < 1e-8 + restore +} +foreach h of numlist 2 3 { + preserve + * clean_h(pre): Dy = lag(y,h) - Ly ; filter on treat, not Fh + gen double Dy = L`h'.y - Ly + keep if obs & !missing(Dy) & !missing(tdiff) & !missing(treat) & (tdiff==1 | treat==0) + teffects ra (Dy x i.time) (tdiff), atet vce(cluster unit) + scalar att_m`h' = _b[r1vs0.tdiff] + scalar se_m`h' = _se[r1vs0.tdiff] + scalar N_m`h' = e(N) + scalar G_m`h' = e(N_clust) + assert reldif(att_m`h', rpt_m`h') < 1e-8 + restore +} + +* ------------------------------------------------------------------------------ +* Emit JSON by hand (Stata has no jsonlite) at %21.17g (round-trip-exact double). +* Per horizon: {att, se, N, G}. 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")) +* Derive the running edition so provenance is correct if regenerated under a +* different edition. 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")) +file open `fh' using "benchmarks/data/lpdid_ra_stata_golden.json", write replace text +file write `fh' "{" _n +file write `fh' `" "meta": {"' _n +file write `fh' `" "estimator": "LPDiD regression-adjustment (RA) SE - Stata teffects ra atet","' _n +file write `fh' `" "generator": "benchmarks/stata/generate_lpdid_ra_golden.do","' _n +file write `fh' `" "source_panel": "benchmarks/data/lpdid_test_panel.csv","' _n +file write `fh' `" "point_anchor": "benchmarks/data/lpdid_golden.json ra_cov[h][0]","' _n +file write `fh' `" "stata_edition": "`sedition'","' _n +file write `fh' `" "cmd": "teffects ra (Dy x i.time) (tdiff), atet vce(cluster unit)","' _n +file write `fh' `" "se_convention": "cluster-robust IF at unit, NO finite-sample factor (teffects convention); t(G-1) reference on the library side","' _n +file write `fh' `" "note": "Independent Stata reconstruction of the per-horizon clean sample (ports generate_lpdid_golden.R prep/clean_h). N and G are emitted to power the Python row-count sample-integrity gate.","' _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' " }," _n +file write `fh' `" "ra_se": {"' + +* post entries (0..4) then pre (-2,-3); comma-PREFIXED so there is no trailing +* comma and no blank first line. Each entry assembled into trimmed locals first. +local sep "" +foreach h of numlist 0 1 2 3 4 { + _jnum att_`h' + local a = r(s) + _jnum se_`h' + local s = r(s) + _jnum N_`h' "%12.0f" + local nn = r(s) + _jnum G_`h' "%12.0f" + local gg = r(s) + file write `fh' "`sep'" _n `" "`h'": {"att": `a', "se": `s', "N": `nn', "G": `gg'}"' + local sep "," +} +foreach h of numlist 2 3 { + _jnum att_m`h' + local a = r(s) + _jnum se_m`h' + local s = r(s) + _jnum N_m`h' "%12.0f" + local nn = r(s) + _jnum G_m`h' "%12.0f" + local gg = r(s) + file write `fh' "`sep'" _n `" "-`h'": {"att": `a', "se": `s', "N": `nn', "G": `gg'}"' + local sep "," +} +file write `fh' _n " }" _n +file write `fh' "}" _n +file close `fh' + +display "Wrote benchmarks/data/lpdid_ra_stata_golden.json (7 horizons)" diff --git a/docs/benchmarks.rst b/docs/benchmarks.rst index 6215b4c1e..b6b0dcc16 100644 --- a/docs/benchmarks.rst +++ b/docs/benchmarks.rst @@ -1,6 +1,6 @@ .. meta:: - :description: Validation benchmarks comparing diff-diff against R packages (did, synthdid, fixest). Coefficient accuracy, standard error comparison, and performance metrics. - :keywords: difference-in-differences benchmark, DiD validation R, python econometrics accuracy, did package comparison + :description: Validation benchmarks comparing diff-diff against R packages (did, synthdid, fixest) and Stata (teffects). 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 ========== @@ -18,14 +18,15 @@ the tables below runs the same wheel with the backend disabled Overview -------- -diff-diff is validated against the following R packages: +diff-diff is validated primarily against established R packages, plus Stata +where no runnable R reference exists: .. list-table:: :header-rows: 1 :widths: 30 30 40 * - diff-diff Estimator - - R Package + - Reference Package - Reference * - ``DifferenceInDifferences`` - ``fixest::feols`` @@ -42,6 +43,9 @@ diff-diff is validated against the following R packages: * - ``SyntheticDiD`` - ``synthdid::synthdid_estimate`` - Arkhangelsky et al. (2021) + * - ``LPDiD`` (regression-adjustment SE) + - Stata ``teffects ra ... atet`` + - Dube, Girardi, Jorda & Taylor (2025); no runnable R analogue Methodology ----------- @@ -51,8 +55,9 @@ Validation Approach 1. **Synthetic Data**: Generate data with known true effects using ``generate_did_data()`` from diff_diff.prep -2. **Identical Inputs**: Both Python and R estimators receive the same CSV data -3. **JSON Interchange**: R scripts output JSON for comparison +2. **Identical Inputs**: Both Python and the reference implementation (R, or Stata + where no runnable R analogue exists) receive the same data +3. **JSON Interchange**: reference generators (R / Stata) output JSON for comparison 4. **Automated Comparison**: Python script validates numerical equivalence 5. **Multiple Scales**: Test at small (200-400 obs), 1K, 5K, 10K, and 20K unit scales 6. **Replicated Timing**: multiple fresh-subprocess replications per benchmark @@ -1005,6 +1010,16 @@ 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: + + .. 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 + Running Benchmarks ~~~~~~~~~~~~~~~~~~ diff --git a/docs/doc-deps.yaml b/docs/doc-deps.yaml index 36554e505..e034e059c 100644 --- a/docs/doc-deps.yaml +++ b/docs/doc-deps.yaml @@ -657,6 +657,8 @@ sources: - path: docs/methodology/REGISTRY.md section: "LPDiD" type: methodology + - path: benchmarks/stata/README.md + type: methodology - path: docs/api/lpdid.rst type: api_reference - path: README.md diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index 36ce13205..d1feb82f3 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -2127,7 +2127,7 @@ Eq. 12 reuses the absorbing clean control and only restricts the treated set (a 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. -2. **Note:** The regression-adjustment (RA) covariate path (`reweight=True` with covariates/absorb) reports an **influence-function cluster variance** `sum_c (sum_{i in c} psi_i)^2 / n^2`, in the same family as `ImputationDiD`'s Theorem-3 / BJS variance (see "IF-based variance estimators vs analytical-sandwich estimators" above). Its single Gram inversion is routed through `linalg._rank_guarded_inv` (finite SE on the identified subspace under near-collinearity; NaN at rank 0). Unlike the default/weighted `solve_ols` `hc1`-cluster path - which applies the `(G/(G-1))*((n-1)/(n-k))` finite-sample factor - the RA IF variance carries **no finite-sample factor**, while both paths share the `t(G-1)` reference. **PR-B2 validated this asymmetry as faithful to the authors' own tooling**, not a defect: the no-factor RA convention matches the canonical Stata `teffects ra ... atet vce(cluster)` (the authors' `lpdid_regression_adjustment.do` `margins`/`kmatch` degrees-of-freedom comments prove `teffects` applies neither factor), while the default path matches `feols`/`reghdfe`. The RA *point* estimate is R-anchored to ~1e-13 (full-interaction `i.dtreat##(i.time c.x)` == `teffects` point; `tests/test_methodology_lpdid.py::test_ra_covariate_point`). The RA *standard error* itself has **no runnable R reference** (no R package computes the RA IF variance - `alexCardazzi` uses direct covariate inclusion, not RA; the canonical RA SE is Stata `teffects` only), so it is **pinned** as a documented regression value (`test_ra_covariate_se_regression_pin`) and its calibration is validated by the ungated Monte-Carlo coverage study `benchmarks/python/coverage_lpdid_ra.py` (~0.95 empirical coverage of the true effect at cluster counts G in {30, 100, 300}). +2. **Note:** The regression-adjustment (RA) covariate path (`reweight=True` with covariates/absorb) reports an **influence-function cluster variance** `sum_c (sum_{i in c} psi_i)^2 / n^2`, in the same family as `ImputationDiD`'s Theorem-3 / BJS variance (see "IF-based variance estimators vs analytical-sandwich estimators" above). Its single Gram inversion is routed through `linalg._rank_guarded_inv` (finite SE on the identified subspace under near-collinearity; NaN at rank 0). Unlike the default/weighted `solve_ols` `hc1`-cluster path - which applies the `(G/(G-1))*((n-1)/(n-k))` finite-sample factor - the RA IF variance carries **no finite-sample factor**, while both paths share the `t(G-1)` reference. **PR-B2 validated this asymmetry as faithful to the authors' own tooling**, not a defect: the no-factor RA convention matches the canonical Stata `teffects ra ... atet vce(cluster)`, 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. 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. diff --git a/docs/methodology/papers/dube-2025-review.md b/docs/methodology/papers/dube-2025-review.md index 0c02f0002..1b035c904 100644 --- a/docs/methodology/papers/dube-2025-review.md +++ b/docs/methodology/papers/dube-2025-review.md @@ -119,6 +119,7 @@ Single overall ATT over the posttreatment window `h in {0..H}` by using the post - Stata: `lpdid` (SSC `s459273`) - the authors' reference implementation. - 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). **Requirements checklist:** - [ ] Long-difference dependent variable with selectable base period (`t-1` default, PMD `k`) diff --git a/tests/test_lpdid_ra_stata_parity.py b/tests/test_lpdid_ra_stata_parity.py new file mode 100644 index 000000000..0419f5a64 --- /dev/null +++ b/tests/test_lpdid_ra_stata_parity.py @@ -0,0 +1,198 @@ +"""External-reference parity: LPDiD regression-adjustment SE vs Stata ``teffects ra``. + +The LP-DiD (Dube, Girardi, Jorda & Taylor 2025) regression-adjustment (RA) covariate +path reports an influence-function cluster variance with **no finite-sample factor**. +No R package computes it (``alexCardazzi/lpdid`` uses direct covariate inclusion, not +RA), so the canonical reference is Stata ``teffects ra ... atet vce(cluster)``. This is +the repo's FIRST Stata parity arm; ``teffects`` is native to Stata (no SSC dependency), +so the generator's ``version 19`` fully pins the numerical behavior. + +The golden is produced by ``benchmarks/stata/generate_lpdid_ra_golden.do``, which +INDEPENDENTLY reconstructs each horizon's clean sample (porting the R +``generate_lpdid_golden.R`` prep/clean_h recipe) and runs ``teffects``. Three gates: + +1. **Point** - Stata ATET vs the R-anchored ``ra_cov[h][0]`` (both committed) at + ``atol=1e-10``, plus a *direct* library-ATET vs Stata-ATET check at the repo's + cross-platform RA-point tolerance: strongly corroborates that the Stata + reconstruction used the same clean sample as R and the library (three-way point + agreement, no longer transitive through R). +2. **Sample shape** - Stata ``(e(N), e(N_clust))`` vs the library's realized + ``(n_obs, n_clusters)`` (exact): narrows the "point matches on a *different* sample" + hole - equal size + cluster count alongside the 1e-10 point and ~1e-16 SE agreement + strongly corroborate the same realized clean sample (strong corroboration, not a + formal membership proof). +3. **SE anchor** - the library RA IF SE vs the Stata ``teffects`` SE. Same-machine + agreement is ~1e-16 (the two compute the identical estimator - IF cluster variance, + no finite-sample factor), but the library side is recomputed on CI and its RA SE is + held only to ``abs=1e-6`` cross-platform (see ``test_ra_covariate_se_regression_pin`` + / ``RA_SE_PIN``), so the assertion uses ``atol=1e-7`` - the repo's LPDiD R-parity SE + tolerance - for BLAS/OS robustness. + +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 with:: + + /Applications/Stata/StataSE.app/Contents/MacOS/stata-se -b do \ + benchmarks/stata/generate_lpdid_ra_golden.do +""" + +from __future__ import annotations + +import json +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_ra_stata_golden.json" +R_GOLDEN_PATH = _DATA / "lpdid_golden.json" +PANEL_PATH = _DATA / "lpdid_test_panel.csv" + +# Both goldens + the shared panel must be present (the point gate reads the R golden). +_FIXTURE_AVAILABLE = ( + STATA_GOLDEN_PATH.is_file() and R_GOLDEN_PATH.is_file() and PANEL_PATH.is_file() +) + +# Horizons pinned by the Stata golden: post {0..4} plus pre placebos {-2, -3}; +# h = -1 is the omitted (zero) reference. +HORIZONS = [0, 1, 2, 3, 4, -2, -3] + +POINT_ATOL = 1e-10 # committed Stata ATET vs committed R ra_cov point (R digits=12) +SE_ATOL = 1e-7 # library RA SE (recomputed) vs committed Stata SE; cross-platform +# library point (recomputed) vs committed Stata ATET; matches the repo's LPDiD +# RA-point cross-platform tolerance (test_methodology_lpdid.py::test_ra_covariate_point). +POINT_XPLATFORM_ATOL = 1e-6 + + +def _skip_if_missing() -> None: + if not _FIXTURE_AVAILABLE: + pytest.skip( + "Stata LPDiD-RA parity fixture not present. Regenerate via " + "`stata-se -b do benchmarks/stata/generate_lpdid_ra_golden.do`." + ) + + +@pytest.fixture(scope="module") +def stata_golden() -> dict: + _skip_if_missing() + with STATA_GOLDEN_PATH.open() as f: + return json.load(f)["ra_se"] + + +@pytest.fixture(scope="module") +def r_ra_points() -> dict: + _skip_if_missing() + with R_GOLDEN_PATH.open() as f: + return json.load(f)["ra_cov"] + + +@pytest.fixture(scope="module") +def library_fit() -> dict: + """Fit the RA path once; return {horizon: {coef, se, n_obs, n_clusters}} + warnings. + + Mirrors the generator's configuration exactly: + ``LPDiD(pre_window=3, post_window=4, reweight=True, cluster="unit")`` with + ``covariates=["x"]``. + """ + _skip_if_missing() + panel = pd.read_csv(PANEL_PATH) + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + res = LPDiD(pre_window=3, post_window=4, reweight=True, cluster="unit").fit( + panel, + covariates=["x"], + outcome="y", + unit="unit", + time="time", + treatment="treat", + ) + messages = [str(w.message) for w in caught] + es = res.event_study.set_index("horizon") + by_h = { + int(h): { + "coef": float(es.loc[h, "coefficient"]), + "se": float(es.loc[h, "se"]), + "n_obs": int(es.loc[h, "n_obs"]), + "n_clusters": int(es.loc[h, "n_clusters"]), + } + for h in HORIZONS + } + return {"by_h": by_h, "warnings": messages} + + +@pytest.mark.parametrize("h", HORIZONS) +def test_stata_point_matches_r_anchor(h, stata_golden, r_ra_points): + """Gate 1: the independent Stata reconstruction reproduces the R-anchored ATET. + + ``rel=0`` makes this a purely absolute ``abs=POINT_ATOL`` bound, independent of + pytest's default relative tolerance. + """ + stata_att = stata_golden[str(h)]["att"] + r_att = r_ra_points[str(h)][0] # ra_cov[h] = [att, conditional_CR0_se_ref_only, null] + assert stata_att == pytest.approx(r_att, abs=POINT_ATOL, rel=0), ( + f"h={h}: Stata ATET {stata_att} != R anchor {r_att} - the Stata clean-sample " + "reconstruction diverged from R; the SE anchor below would be meaningless." + ) + + +@pytest.mark.parametrize("h", HORIZONS) +def test_stata_sample_shape_matches_library(h, stata_golden, library_fit): + """Gate 2: Stata (e(N), e(N_clust)) == library (n_obs, n_clusters), exact. + + Equal size AND cluster count do not by themselves *prove* identical (unit, time) + membership, but combined with the 1e-10 point gate and the ~1e-16 SE match they + provide strong corroboration of the same realized clean sample. + """ + stata_n = int(stata_golden[str(h)]["N"]) + stata_g = int(stata_golden[str(h)]["G"]) + lib = library_fit["by_h"][h] + assert (stata_n, stata_g) == (lib["n_obs"], lib["n_clusters"]), ( + f"h={h}: Stata (N={stata_n}, G={stata_g}) != library " + f"(n_obs={lib['n_obs']}, n_clusters={lib['n_clusters']}) - the clean samples differ." + ) + + +@pytest.mark.parametrize("h", HORIZONS) +def test_library_point_matches_stata_teffects(h, stata_golden, library_fit): + """Direct Python<->Stata ATET agreement (not transitive through the R golden). + + Gate 1 pins Stata vs R (committed values); this closes the loop by comparing the + recomputed library ATET directly to the Stata ATET, at the repo's cross-platform + RA-point tolerance (``rel=0`` -> purely absolute). + """ + stata_att = stata_golden[str(h)]["att"] + lib_coef = library_fit["by_h"][h]["coef"] + assert lib_coef == pytest.approx( + stata_att, abs=POINT_XPLATFORM_ATOL, rel=0 + ), f"h={h}: library ATET {lib_coef} != Stata teffects ATET {stata_att}" + + +def test_library_fit_emits_no_ra_drop_warning(library_fit): + """Gate 2 (companion): no treated observation was dropped for non-identification. + + A drop would mean the library and Stata clean samples could diverge silently; on + the committed panel every horizon is identified (G=60, all treated cells have a + clean control), so the fit must be warning-clean. + """ + drops = [m for m in library_fit["warnings"] if "regression adjustment: dropped" in m] + assert not drops, f"unexpected RA drop warning(s): {drops}" + + +@pytest.mark.parametrize("h", HORIZONS) +def test_library_ra_se_matches_stata_teffects(h, stata_golden, library_fit): + """Gate 3 (the anchor): library RA IF SE == Stata teffects atet cluster SE. + + Confirms the no-finite-sample-factor RA convention (REGISTRY LPDiD Deviation 2) + against an independent implementation. Local agreement is ~1e-16; SE_ATOL is the + repo's cross-platform LPDiD SE tolerance (see module docstring). ``rel=0`` makes + this a purely absolute ``abs=SE_ATOL`` bound, independent of pytest's default rel. + """ + stata_se = stata_golden[str(h)]["se"] + lib_se = library_fit["by_h"][h]["se"] + assert lib_se == pytest.approx( + stata_se, abs=SE_ATOL, rel=0 + ), f"h={h}: library RA SE {lib_se} != Stata teffects SE {stata_se}"