diff --git a/CHANGELOG.md b/CHANGELOG.md index af7f45c0..91bf6f6a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,30 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] ### Added +- **LWDiD Stata parity arm + acceptance-suite re-anchor (ahead of PR #588's merge).** + New golden generator `benchmarks/stata/generate_lwdid_golden.do` runs the + authors' own SSC `lwdid` package (v2.4.2, version line pinned in the + golden's meta) over the loader-cached Prop 99 / Walmart panels and the + committed castle subset, emitting `benchmarks/data/lwdid_stata_golden.json`: + full-precision small-N ATT/SE (Prop 99 both rollings + castle `tau_omega`), + Prop 99 randomization-inference p at 100k reps, and the six Walmart + event-study configs' per-r WATT points + multiplier-bootstrap SEs at + B=9,999. The maintainer validation suite `tests/test_methodology_lwdid.py` + is rewritten against the agreed post-rename LWDiD API as the acceptance + bar for the contribution's final round: the RI golden re-anchors from the + paper's printed 0.020 to the authors'-package inclusive-ties convention + (~0.051 — the printed value is not reproducible with the package; REGISTRY + LWDiD RI note), and the Walmart SE xfail scaffolding is retired: the RA + config's bootstrap SEs become a passing Monte-Carlo-tolerance comparison + against the Stata golden, while the four IPWRA SE columns are committed as + provenance only — they diverge systematically (~15%) pending the E.3 + influence-function adjudication in the contribution's final round (REGISTRY + LWDiD IPWRA-variance note). A new ungated + `tests/test_lwdid_stata_golden_schema.py` validates the golden's schema on + main (the methodology suite import-skips until `diff_diff.lwdid` lands). + REGISTRY's LWDiD section gains the maintainer implementation notes (RI + convention, dual overall-ATT conventions, Sec 4.3 rejection, inference + default, API conformance). - **`results.aggregate('total')` - the estimator-owned total incremental outcome** on CallawaySantAnna, EfficientDiD, ImputationDiD, and TwoStageDiD, promoted into the library-wide aggregation vocabulary. The single diff --git a/DEFERRED.md b/DEFERRED.md index 5048e077..32ec54b2 100644 --- a/DEFERRED.md +++ b/DEFERRED.md @@ -75,6 +75,7 @@ For survey-specific limitations (`NotImplementedError` paths), see the | Issue | Location | PR | Priority | |-------|----------|----|----------| | RDDensityTest public bandwidth-selector helper (an `rdbwdensity`-equivalent exposed function) and a Tutorial-28 executed density-test demo cell (Act 4d; needs notebook re-execution + drift-suite updates): the selector ships internal-only and the tutorial carries a prose pointer for now. | `diff_diff/rddensity.py`, `docs/tutorials/28_rdd_scholarship_illusion.ipynb` | rddensity PR-B | Low | +| LW 2025 Sec 4.3 all-eventually-treated mode for LWDiD (drop `D_infinity`, effects relative to the last cohort as reference, last cohort's own effect not estimable). Current behavior (by decision, see the REGISTRY LWDiD Sec 4.3 note): such designs raise ValueError under both control strategies rather than silently truncating | `diff_diff/lwdid_staggered.py` | #588 | Low | | EfficientDiD/ImputationDiD/TwoStageDiD kit retention has no memory opt-out: EfficientDiD's kit always retains the per-(g,t) EIF dict (M-023), ImputationDiD's references the `_fit_data` panel objects (M-021, pre-existing retention), and TwoStageDiD's retains a column-subset working-frame copy (M-022, the first new panel retention); a `store_kit=False` constructor knob would restore an opt-out at the cost of post-fit aggregation on that result | `diff_diff/efficient_did.py`, `diff_diff/imputation.py`, `diff_diff/two_stage.py` | 2(b) PR-3a | Low | | MMM interop, StackedDiD totals (remainder of the estimator-owned-totals delivery - `aggregate('total')` shipped for CS/EfficientDiD/ImputationDiD/TwoStageDiD on panel non-survey fits): StackedDiD is STAGED OUT because the total estimand is ill-defined under `weighting="population"`/`"sample_share"` (a weighted estimand for which raw treated exposure is not the right mass), and its results class is view-only (a total would be its first computed level). Needs an estimand decision per weighting variant before any implementation. | `diff_diff/mmm.py`, `diff_diff/stacked.py` | mmm-interop | Low | | `SyntheticControl` fit-snapshot residency (`_SyntheticControlFitSnapshot`) — **investigated 2026-07-07, parked**: the snapshot ALIASES the fit's own working pivots (zero extra construction cost); the retained residency implements the documented freeze contract (post-fit mutation of estimator inputs must not change `in_space_placebo()` / `leave_one_out()` / conformal output on an already-returned results object, and `__getstate__` already excludes it from pickles). A compact array representation saves only pandas overhead (the float panel dominates); releasing residency needs new API surface (`release`/opt-out flag) or a freeze-contract change. Revisit on user demand for very large donor panels. | `synthetic_control.py`, `synthetic_control_results.py` | follow-up | Low | @@ -131,6 +132,7 @@ decisions (refactor waivers, perf trade-offs, test-infrastructure calls) are rec | Decision | Location | Verified | |----------|----------|----------| +| **Fixed-B Stata SE-golden comparisons in the LWDiD suite run unmarked in default CI** (no `ci_params.bootstrap()` scaling, no `slow` marker): the SE tolerance derives from BOTH fixed rep counts (ours B=999, Stata's R from the golden meta) so scaling either side would invalidate the committed-golden comparison — the ci_params convention governs convergence-style tests, not fixed-B golden comparisons — and these tests are the PR #588 acceptance bar, which must run in default CI; module-scoped fit memoization bounds the cost (~1-2 s per fit locally) | `tests/test_methodology_lwdid.py` | 2026-08-16 | | **The 4.0 migration guide's code blocks are not snippet-executed.** `tests/test_doc_snippets.py` discovers a hardcoded list of `.rst` files and only `.. code-block:: python` / RST `::` bodies, so `docs/migration-4.0.md` gets no coverage. Deliberate: the guide is a MIXED document - most "after" examples (the renames, `results.att`) run on the current release, but the `field-flip` and `df-convention-flip` examples describe 4.0 behaviour that cannot run until 4.0, so a blanket execution lane would fail by construction. Closing the gap fully means a markdown-fence extractor plus a skip-marker convention for the future-API blocks - a harness change, out of scope for a docs PR. What IS gated: the appendix's ledger parity (`test_migration_guide_*`), which pins the row set and every mechanically checkable cell; and, since the first local review found all three merge examples carrying invalid keywords, `test_migration_guide_examples_bind_to_real_signatures`, which ast-parses the guide's python blocks and asserts every constructor/`fit()` keyword exists on the target signature. That is signature binding, NOT execution - it deliberately skips calls whose owner it cannot resolve (e.g. `results.aggregate(...)`), and it cannot catch a wrong VALUE or a wrong sequence of calls. The hand-written `Fix` prose remains unverifiable by any available means. | `docs/migration-4.0.md`, `tests/test_doc_snippets.py` | Phase 4 / 2026-08-09 | | **MultiPeriodDiD deprecation shim loses static constructor-arg checking (3.9 window).** The M-010 shim is `__init__(*args, **kwargs)` + an import-time `__signature__` mirror of DiD's constructor: runtime introspection (get_params/set_params, `inspect.signature`) and eager validation are fully preserved, but static type checkers / IDEs cannot check constructor arguments for the deprecated class until its 4.0 removal. Accepted: the alternative (hand-mirroring ~20 parameters) is a drift magnet on a class with one minor version of remaining life. | `diff_diff/estimators.py` | 3(a) / 2026-08-07 | | **DCDH `sklearn.base.clone` param-identity failure won't-fix.** `ChaisemartinDHaultfoeuille._validate_paths_of_interest` unconditionally canonicalizes `paths_of_interest` into a fresh `List[Tuple[int, ...]]`, so sklearn `clone()`'s post-construction `param1 is param2` identity check fails for configured instances - a pre-existing normalization the BaseEstimator mixin PR documented rather than changed (get_params/set_params signatures are clone-compatible; the dependency-free `cls(**est.get_params())` config-equality contract is the enforced one, `tests/test_base_estimator.py`). Fixing would mean returning the caller's raw object from a validator whose job is canonicalization. | `chaisemartin_dhaultfoeuille.py` | mixin PR / 2026-08-01 | diff --git a/TODO.md b/TODO.md index 81cbad4b..2ef1d085 100644 --- a/TODO.md +++ b/TODO.md @@ -21,6 +21,7 @@ Related tracking surfaces: | Issue | Location | Origin | Effort | Priority | |-------|----------|--------|--------|----------| +| Expose cell-mass overall ATT (Stata `Post_avg` convention; = CS-simple on balanced panels) as an aggregate extra on LWDiD results — the fit's `.att` is the paper's `tau_omega` (cohort-mean-then-treated-weight, eq. 7.18); the authors' large-N display uses cell-mass weighting instead, and both are legitimate estimands (see the REGISTRY LWDiD Aggregation note). Lands only after PR #588 merges | `diff_diff/lwdid_results.py` | #588 | Quick | Low | | Post-fit `aggregate()` for the staggered DDD container: `StaggeredTripleDiffResults` carries no `AggregationMixin`, which is why the phase-3(b) merge had to carry fit-time `aggregate=`/`balance_e=` onto the surviving `TripleDifference` (rows M-140/M-141) as the ONE documented exception to the section-6 aggregate-postfit program. Porting the container onto the M-122 aggregation contract retires both rows; note the bootstrapped-fit recompute levels will need draw retention or a fail-closed relay, the same problem tracked for CS/EfficientDiD/ImputationDiD. Until it lands, the DDD docs deliberately keep teaching the fit-time kwarg (the canonical route there) | `diff_diff/staggered_triple_diff_results.py`, `diff_diff/aggregation.py`, `docs/api/triple_diff.rst`, `docs/tutorials/08_triple_diff.ipynb` | 3(b) | Heavy | Medium | | Staggered-DDD power support: `simulate_power`/`simulate_mde`/`simulate_sample_size` now REJECT a staggered-configured `TripleDifference` (both registered DDD generators emit 2x2x2 data and fit with `(group, partition, post)`, so a staggered config would be simulated under the wrong design). Support needs a staggered DDD DGP profile plus fit-kwargs builder, and a decision on whether the mode is selected by profile or by the estimator's own config | `diff_diff/power.py` | 3(b) | Mid | Low | | Bootstrap-`seed` provenance on multiplier-bootstrap results containers: neither `StaggeredTripleDiffResults` nor `CallawaySantAnnaResults` carries the `seed` that generated its bootstrap SEs / p-values / sup-t bands, so a serialized result cannot report the random configuration behind its inference. NOT a 3(b) regression - `seed` reaches the engine and `get_params()` correctly (same seed reproduces the SE bit-exactly, a different seed moves it), the gap is results-object observability only, it predates the merge, and both containers inherit it from the shared `CallawaySantAnnaBootstrapMixin`. Add `seed` (and consider `n_bootstrap`/`bootstrap_weights`/`cband`) to BOTH containers plus `to_dict()`, with seeded and unseeded pins; sequence it with the M-014 container unification rather than schema-changing one container mid-merge. Precedent for exposing it: `ContinuousDiDResults`, `EfficientDiDResults`, `SyntheticDiDResults` already do | `diff_diff/staggered_triple_diff_results.py`, `diff_diff/staggered_results.py` | 3(b) | Quick | Low | diff --git a/benchmarks/README.md b/benchmarks/README.md index f7604893..429051e7 100644 --- a/benchmarks/README.md +++ b/benchmarks/README.md @@ -57,12 +57,14 @@ $STATA -b do benchmarks/stata/generate_lpdid_ra_golden.do $STATA -b do benchmarks/stata/generate_imputation_loo_golden.do $STATA -b do benchmarks/stata/generate_etwfe_cs_golden.do $STATA -b do benchmarks/stata/generate_reghdfe_kref_golden.do +$STATA -b do benchmarks/stata/generate_lwdid_golden.do # see stata/README.md for its warm-up step ``` The `LPDiD` arm uses only **native** Stata commands (`teffects`), pinned by `version 19`. The `ImputationDiD` arm depends on SSC packages (`did_imputation`/`reghdfe`/`ftools`/`require`), the ETWFE/CS arm on -`drdid`/`csdid`/`jwdid`/`hdfe`, and the reghdfe K_reference arm on `reghdfe`; +`drdid`/`csdid`/`jwdid`/`hdfe`, the reghdfe K_reference arm on `reghdfe`, +and the LWDiD arm on the authors' `lwdid`; `version 19` does NOT pin SSC packages (SSC has no version history) — install them once via `requirements.do` (the generators do not auto-install) and each golden records version/checksum @@ -87,7 +89,8 @@ benchmarks/ │ ├── generate_lpdid_ra_golden.do # LPDiD RA SE vs teffects ra │ ├── generate_imputation_loo_golden.do # ImputationDiD LOO SE vs did_imputation leaveout │ ├── generate_etwfe_cs_golden.do # ETWFE/CS vs jwdid + csdid (+ subsample ladder) -│ └── generate_reghdfe_kref_golden.do # clustered CR1 K_reference vs reghdfe (disconnected panel) +│ ├── generate_reghdfe_kref_golden.do # clustered CR1 K_reference vs reghdfe (disconnected panel) +│ └── generate_lwdid_golden.do # LWDiD vs the authors' lwdid (small-N, RI, event-study bootstrap) ├── python/ │ ├── utils.py # Common utilities │ ├── benchmark_callaway.py # CallawaySantAnna diff --git a/benchmarks/data/lwdid_stata_golden.json b/benchmarks/data/lwdid_stata_golden.json new file mode 100644 index 00000000..cb91de60 --- /dev/null +++ b/benchmarks/data/lwdid_stata_golden.json @@ -0,0 +1,290 @@ +{ + "meta": { + "estimator": "LWDiD (Lee & Wooldridge 2025, 2026) - authors' Stata lwdid package parity", + "generator": "benchmarks/stata/generate_lwdid_golden.do", + "ssc_versions": {"lwdid": "version 2.4.2 15 June 2026"}, + "stata_edition": "SE", + "bootstrap_scheme": "lwdid large-N multiplier bootstrap (package default; reps() draws with set seed); compared against diff-diff's unit-level Rademacher multiplier bootstrap at B=999", + "control_pool": {"small_N": "composite regression, never-treated-based by construction (no control option exists); equivalent to Python control_group='never_treated'", "large_N": "never-treated + not-yet-treated default (the 'never' option is NOT passed); matches Python default control_group='not_yet_treated'"}, + "datasets": {"prop99": {"url": "http://fmwww.bc.edu/repec/bocode/l/lw_smoking.dta", "sha256": "16c3ac1da351788817433fc890ec2f502a8bdfcb46cbc8d693653330e71d5a65"}, "walmart": {"url": "http://fmwww.bc.edu/repec/bocode/l/lw_walmart.dta", "sha256": "410885572143dceb9daa643a8097768f1bc3493f9437451a9e4d1d5dc1e18d14"}, "castle": {"path": "benchmarks/data/real/castle_lw_subset.csv"}}, + "stata_version": 19.0, + "rireps": 100000, + "riseed": 20260815, + "bootstrap_reps": 9999, + "bootstrap_seed": 20260815 + }, + "prop99": { + "demean": {"att": -0.42217461502012654, "se": 0.12079952386677335, "p_ri": 0.026749732502674973}, + "detrend": {"att": -0.22698869955616763, "se": 0.094068893843888307, "p_ri": 0.052639473605263949} + }, + "castle": { + "demean": {"att": 0.091745387139613596, "se": 0.057102695603325655}, + "detrend": {"att": 0.066550335128826035, "se": 0.056012386378555774} + }, + "walmart": { + "detrend_ra__log_retail_emp": { + "watt": { + "-22": [-0.032869808983039699, 0.017624710036477202], + "-21": [-0.021320017034415301, 0.010151158896656401], + "-20": [-0.0116210174606169, 0.0083639469796090004], + "-19": [-0.00061104820844, 0.0073885953707725997], + "-18": [-0.0023026261440382998, 0.0059356869986891998], + "-17": [-0.0028910304069036, 0.0042525688546915996], + "-16": [0.0052021302742854998, 0.0033687438708670998], + "-15": [0.0048990622495043001, 0.0030726745105532], + "-14": [0.0026436326348103002, 0.0026767184052707001], + "-13": [-0.0060195910334129004, 0.0027390010384942998], + "-12": [-0.0061429921219214, 0.0021660988380691002], + "-11": [-0.0036769404905967998, 0.0020471717626491001], + "-10": [0.0012255892680701, 0.0022008323968882998], + "-9": [0.0011690810389403999, 0.0022047461508291001], + "-8": [0.0038259121091685001, 0.0021612168484445998], + "-7": [0.0059497199978954, 0.0022031017003438002], + "-6": [0.0038542023917781999, 0.0020452847439453], + "-5": [0.0047282745924016002, 0.0020240572667549], + "-4": [0.0053559210065129998, 0.0019114316013613001], + "-3": [-0.0000178099765863, 0.0016779961906666001], + "-2": [0, 0], + "-1": [0, 0], + "0": [0.0072331135173352001, 0.0036635874455351], + "1": [0.032212750139948398, 0.0051302591263951], + "2": [0.0240114130515417, 0.0063392958928369996], + "3": [0.018115740657803001, 0.0078007498337440999], + "4": [0.0150672686437532, 0.0091322647377050006], + "5": [0.0164467023898466, 0.010918123423304799], + "6": [0.0197611344107241, 0.012594966119080101], + "7": [0.035660153543117398, 0.014526941082212899], + "8": [0.040341176071916798, 0.016391654044761499], + "9": [0.037549751765288797, 0.019155799882757101], + "10": [0.023958717335560601, 0.023376790404146801], + "11": [0.0000960360982631, 0.028761490367063901], + "12": [-0.0023635346661423002, 0.037012184957851198], + "13": [0.027538329065794301, 0.054609122586710197] + }, + "overall": { + "Pre_avg": [0.0012741669959059, 0.0003094184778805], + "Post_avg": [0.0219775794324495, 0.0101621100306762] + } + }, + "detrend_ipwra__log_retail_emp": { + "watt": { + "-22": [-0.036379516584030401, 0.016859957226860001], + "-21": [-0.0233149790840733, 0.0094037484122564997], + "-20": [-0.011027271225191201, 0.0079153126461820999], + "-19": [-0.0007168158562135, 0.0069641339415532997], + "-18": [-0.0046787480270942998, 0.0056844647941504001], + "-17": [-0.0041303179839107997, 0.0039897592895552], + "-16": [0.0048505478178386, 0.0031644442670286], + "-15": [0.0044040650363834999, 0.0028938431529255001], + "-14": [0.0028112497717054001, 0.0026002741184215999], + "-13": [-0.0054165655796852997, 0.0025964617540279], + "-12": [-0.0060993819029491001, 0.0020947192022076002], + "-11": [-0.0033211951703885998, 0.0020048639833157001], + "-10": [0.0013240077646047999, 0.0020874992799277], + "-9": [0.0016277734465345, 0.0020768297213930001], + "-8": [0.0044241819270507002, 0.0020625660031432002], + "-7": [0.0064614313755062, 0.0021503070285278001], + "-6": [0.0033130263195181998, 0.0019875453399713001], + "-5": [0.0048859348344042, 0.0019502418228817], + "-4": [0.0050924505736281997, 0.0018357464510820999], + "-3": [0.0002513706389189, 0.0015873934773629999], + "-2": [0, 0], + "-1": [0, 0], + "0": [0.0066771286223030997, 0.0035706718994881], + "1": [0.031907034937737, 0.0050501392557210003], + "2": [0.0240026554302009, 0.0062605125627464], + "3": [0.019460903515924902, 0.0077190353580456001], + "4": [0.0159784513647818, 0.0091123640934158995], + "5": [0.015692360835684399, 0.010845762946476101], + "6": [0.016981627212872299, 0.012628491673815599], + "7": [0.033304463585542197, 0.0143921022930162], + "8": [0.038500545923522998, 0.016065688849282701], + "9": [0.039180179428989999, 0.018379749829394901], + "10": [0.0350779602807312, 0.021983038422147098], + "11": [0.016344964680412199, 0.026545678089310899], + "12": [0.0149337945856171, 0.033085893907452597], + "13": [0.045984500886364098, 0.043122591591231003] + }, + "overall": { + "Pre_avg": [0.0013561892258513999, 0.0003054220826937], + "Post_avg": [0.0228723314177301, 0.0097836143482841007] + } + }, + "demean_ipwra__log_retail_emp": { + "watt": { + "-22": [-0.1265667598306377, 0.032931268156581303], + "-21": [-0.1026851197654966, 0.0235327295408422], + "-20": [-0.10130695840767839, 0.017533610163888699], + "-19": [-0.087464578552356995, 0.016178725705569399], + "-18": [-0.063348840778495105, 0.011697670645966799], + "-17": [-0.057257902370084802, 0.0086675828188626999], + "-16": [-0.041830671140123898, 0.0071151289773139], + "-15": [-0.029885368181652501, 0.0062773935142942999], + "-14": [-0.022554327102101399, 0.0054487247995867001], + "-13": [-0.019603637734096301, 0.0048876585476156], + "-12": [-0.0158875541467235, 0.0039518581921481003], + "-11": [-0.0117357190832646, 0.0033140256228865001], + "-10": [-0.0051802742616998999, 0.0030223805764341], + "-9": [-0.0019979266912907002, 0.0026944764689437999], + "-8": [0.0049509384066898999, 0.0024339708510338], + "-7": [0.0111406444726097, 0.0025189166384971999], + "-6": [0.012144696034085901, 0.0026004190187776002], + "-5": [0.017870061166436101, 0.0030428568217146998], + "-4": [0.0222290335231243, 0.0035116674900589999], + "-3": [0.0215404102058794, 0.0038897869498251002], + "-2": [0.020046028089758701, 0.0043963628970100001], + "-1": [0, 0], + "0": [0.0404235380416565, 0.0054313798840825997], + "1": [0.072211736332771406, 0.0064235331362959003], + "2": [0.071909394356022102, 0.0074688921648073999], + "3": [0.074169561112794394, 0.0086548448387444004], + "4": [0.079111442339598395, 0.0097528159992984996], + "5": [0.090527687992410094, 0.0108525912787115], + "6": [0.1011768365401205, 0.012077602732584], + "7": [0.1216548754381906, 0.0131949734659286], + "8": [0.1366012794805587, 0.0143858897904555], + "9": [0.14437865372909131, 0.016216201104351301], + "10": [0.1681346396438935, 0.019142795455582701], + "11": [0.1783024590143957, 0.023221750323480499], + "12": [0.181936603107721, 0.029684161737235201], + "13": [0.23051619104812621, 0.042992346638219901] + }, + "overall": { + "Pre_avg": [-0.0014088984066793001, 0.0003735591081783], + "Post_avg": [0.094983398621604501, 0.0099759478833273008] + } + }, + "detrend_ra__log_wholesale_emp": { + "watt": { + "-22": [-0.0166781630345583, 0.0249492193583177], + "-21": [-0.0134691682024187, 0.0265168102485593], + "-20": [0.026035636445363899, 0.019027857265189501], + "-19": [0.0049602149724297002, 0.0178840558847141], + "-18": [0.017197316166776399, 0.0133094001814444], + "-17": [0.0074859656420785004, 0.011136686024101301], + "-16": [0.0078825140645386992, 0.0073619761132143998], + "-15": [-0.0093029431286250006, 0.0062763506963301003], + "-14": [-0.0093539419100747007, 0.0055937395036290997], + "-13": [-0.0092153264224708997, 0.0047957942232543004], + "-12": [0.0011165858906435999, 0.0047207752574429003], + "-11": [0.0050323342119105996, 0.0044338702678121998], + "-10": [0.0017875518835511999, 0.0050205976612657001], + "-9": [-0.0013154102176849999, 0.0047132142516927], + "-8": [-0.0043346277557406998, 0.0048052649147116997], + "-7": [0.0010069632358766, 0.0046482366121568001], + "-6": [0.0049587667225499, 0.0044095144398822997], + "-5": [-0.0014905350715235001, 0.0041819164090552002], + "-4": [0.0029289676004634999, 0.0039579143747065001], + "-3": [0.0017760208016357, 0.0039843286246427998], + "-2": [0, 0], + "-1": [0, 0], + "0": [-0.0005069292346066, 0.0085100172113640998], + "1": [-0.016846810748210699, 0.0106210615627813], + "2": [-0.0284552384106747, 0.0131526737742133], + "3": [-0.038859293296893901, 0.016022745703066301], + "4": [-0.045758518493014502, 0.019108628831135299], + "5": [-0.044836384811121503, 0.022751982472146201], + "6": [-0.042960975866873599, 0.026176340371541498], + "7": [-0.039660094925869298, 0.0305793973475528], + "8": [-0.055117198740565101, 0.034935840023438601], + "9": [-0.058070644974873097, 0.040108491818666901], + "10": [-0.054003601550404401, 0.0467318924104691], + "11": [-0.1007609253610636, 0.057564991726528299], + "12": [-0.13516963107884869, 0.067835632181930103], + "13": [-0.087657900312403295, 0.085535360168670596] + }, + "overall": { + "Pre_avg": [0.0002711857038081, 0.0007320396867235], + "Post_avg": [-0.039034281629469897, 0.020773624819511401] + } + }, + "detrend_ipwra__log_wholesale_emp": { + "watt": { + "-22": [-0.019764386787799799, 0.024676054681120399], + "-21": [-0.0178217797296569, 0.023173086008546701], + "-20": [0.031714915614891398, 0.016677016561166901], + "-19": [0.0069114735873840004, 0.015295246143585999], + "-18": [0.017449180256758599, 0.012236442798624301], + "-17": [0.0070166242146020002, 0.010557203322134501], + "-16": [0.0099031820172495004, 0.0067424007414697996], + "-15": [-0.0079105817650971001, 0.0058643022278788001], + "-14": [-0.0088156964435695002, 0.0052675055015838001], + "-13": [-0.0093060294268772006, 0.0045874291334756998], + "-12": [-0.0002291244656504, 0.0045129153339050999], + "-11": [0.0034128484150541, 0.0043174467222244], + "-10": [-0.0003849064017842, 0.0049519772121045004], + "-9": [-0.0016674045939451001, 0.0045915474697652003], + "-8": [-0.0034144455258635002, 0.004687198428712], + "-7": [0.0015998753558063, 0.0045234990360592999], + "-6": [0.0062285881148214997, 0.0043255624349487999], + "-5": [-0.0005767423630826, 0.0041321879149434004], + "-4": [0.0032423471381093998, 0.0039079361235291], + "-3": [0.001073622195384, 0.003886495169969], + "-2": [0, 0], + "-1": [0, 0], + "0": [-0.0023122056658507, 0.0083989966830585999], + "1": [-0.019535535524106602, 0.010631007939586001], + "2": [-0.031630167724447399, 0.013346508707818399], + "3": [-0.043944027661363801, 0.0163909226165036], + "4": [-0.052026092956788103, 0.019778157754672999], + "5": [-0.052411492632748601, 0.0236629535155604], + "6": [-0.056548554036921897, 0.027385820791360298], + "7": [-0.051844321886251502, 0.031936723097364401], + "8": [-0.066694517861813005, 0.0362014503817518], + "9": [-0.062490021296021901, 0.041323806556528397], + "10": [-0.052771518159207699, 0.0492058332120279], + "11": [-0.089105382421253695, 0.060360236212299898], + "12": [-0.1135089103074004, 0.068194934789793202], + "13": [-0.061573519276753301, 0.082660130560822498] + }, + "overall": { + "Pre_avg": [0.0002982690648296, 0.0007306411115489], + "Post_avg": [-0.043900207742629802, 0.021319152953250101] + } + }, + "demean_ipwra__log_wholesale_emp": { + "watt": { + "-22": [-0.1224387461546648, 0.051695933924493398], + "-21": [-0.1068198066158309, 0.036705055844681603], + "-20": [-0.060619633753150402, 0.029255774675624398], + "-19": [-0.0691612907620026, 0.0240491673625358], + "-18": [-0.0547882983631367, 0.021421805650571901], + "-17": [-0.056256032382420898, 0.018214931997429299], + "-16": [-0.053089300567689901, 0.0140020410815971], + "-15": [-0.046324385105713801, 0.0119291628424739], + "-14": [-0.040605676168759, 0.0096631056067872992], + "-13": [-0.026961042565563299, 0.0088116890007322], + "-12": [-0.011828204095582599, 0.0075501428084263], + "-11": [-0.0071423471463131997, 0.0061998992881783996], + "-10": [-0.0106360540950203, 0.0059664044858119998], + "-9": [-0.0083068767967495998, 0.0052774799011478003], + "-8": [-0.0044867848784985002, 0.0050076079287419001], + "-7": [0.0060946688533405, 0.0049877637594234999], + "-6": [0.016290514462525101, 0.0053627900334861003], + "-5": [0.0150523168347907, 0.0058464119529952003], + "-4": [0.024438539186151999, 0.0064296929075360002], + "-3": [0.027836947093596399, 0.0071566348765411004], + "-2": [0.029583682017902499, 0.0082127438446350995], + "-1": [0, 0], + "0": [0.0411525177828699, 0.0101974189955248], + "1": [0.032204764567774199, 0.011771862647732299], + "2": [0.030447600531221902, 0.0135185772846647], + "3": [0.027855923205975502, 0.015639708941682599], + "4": [0.033149171196456502, 0.017594275172949699], + "5": [0.042306878508044897, 0.019873895569967898], + "6": [0.051454325677815702, 0.021936091218538602], + "7": [0.055119227393426, 0.024157470780753799], + "8": [0.064081530924753605, 0.026358482538956], + "9": [0.074537256919657902, 0.028794405758284399], + "10": [0.12898874442364389, 0.033292969172122003], + "11": [0.1382786833158752, 0.039470959438951098], + "12": [0.1266235757436317, 0.047822421642295797], + "13": [0.16462473590902241, 0.060072341831702] + }, + "overall": { + "Pre_avg": [-0.0028165493242277, 0.0006962843765859], + "Post_avg": [0.050389056036868503, 0.017231313376858199] + } + } + } +} diff --git a/benchmarks/stata/README.md b/benchmarks/stata/README.md index 50c58be1..80f093dc 100644 --- a/benchmarks/stata/README.md +++ b/benchmarks/stata/README.md @@ -1,10 +1,11 @@ # 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. +Stata golden generators live here (five arms: LPDiD `teffects ra`, +ImputationDiD leave-one-out, ETWFE/CS, reghdfe K_reference, LWDiD); the +pattern mirrors `benchmarks/R/` (a `generate_*` script writes a committed +golden JSON that a skip-guarded test reads, so CI never needs Stata). Stata +is node-locked single-user, so — exactly like the R arm — goldens are +committed and only regenerated locally. Locating the binary (macOS, StataSE 19; **not** on `PATH` by default): @@ -133,11 +134,11 @@ The A.9 leave-one-out (LOO) variance has **no runnable R reference** — R internal psi-identity + hand-calc + MC coverage. The authors' own Stata `did_imputation` (Borusyak) ships the same option (`leaveout`); this arm turns it into a measured anchor. -## This is the first SSC-dependent arm +## SSC dependence (vs the native-`teffects` LPDiD 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 +`version 19`. (The ETWFE/CS and LWDiD arms are SSC-dependent in the same way.) The generator does **not** install them — run `benchmarks/stata/requirements.do` once first — and it records each package's version in `meta.ssc_versions` (the `*!` ado header line) so drift is detectable. Byte-identical regeneration is therefore scoped to a fixed Stata + fixed installed SSC versions. @@ -281,13 +282,89 @@ Two arms: sqrt((N−10)/(N−11))`. On CONNECTED designs the two coincide — the jwdid subsample ladder above pins machine-precision agreement at every G. +# LWDiD parity vs the authors' `lwdid` package + +`benchmarks/stata/generate_lwdid_golden.do` produces +`benchmarks/data/lwdid_stata_golden.json`, consumed by +`tests/test_methodology_lwdid.py` (import-skip-gated until `diff_diff.lwdid` +lands via PR #588) and schema-validated on main by the ungated +`tests/test_lwdid_stata_golden_schema.py`. + +## Why the authors' `lwdid` + +Printed-table goldens (LW 2026 Table 3, LW 2025 Tables A4/A5) only reach +three decimals and don't cover the randomization-inference convention or the +multiplier-bootstrap SEs at all. The authors' SSC `lwdid` (Hur, Lee & +Wooldridge) is the reference implementation; this arm pins it at full +precision: Prop 99 small-N ATT/SE + RI p at 100k reps (both rollings), +castle-doctrine `tau_omega` ATT/SE (both rollings), and the six Walmart +large-N event-study configs (per-r WATT points + B=9,999 multiplier-bootstrap +SEs, `set seed` pinned). It is an SSC arm: `meta.ssc_versions` records the +verbatim `lwdid` version line (fail-closed capture — the generator aborts if +no version line is found), and `capture which lwdid` exits 111 if the package +is missing (`requirements.do` installs it). + +## Inputs and their integrity + +Prop 99 / Walmart come from the loader cache +(`~/.cache/diff_diff/datasets/{prop99,walmart}.dta` — the authors' SSC +ancillary files). Loader success does NOT guarantee the on-disk bytes (cache +writes are best-effort), so the regeneration recipe hashes the cache files +directly. Castle uses the committed +`benchmarks/data/real/castle_lw_subset.csv`; note `set type double` makes +`import delimited` read it at double precision (a float import shifts the +8th decimal — the committed golden is the double-precision read, matching +Python). In-`.do` gates: `confirm numeric variable` + row-count asserts per +dataset. Measured-value smoke gates cover the SMALL-N blocks only (prop99 +att + detrend RI p; castle att/se); the Walmart block's gates are +schema/cardinality/nonmissing (one row per event time, no missing cells) - +its value anchoring lives in the Python parity tests. + +## Regenerating + +```bash +# 1. fail-closed warm-up: loaders genuine AND on-disk cache bytes SHA-verified +python -c "from diff_diff import load_prop99, load_walmart; \ + import hashlib, pathlib; \ + dfs = {'prop99': load_prop99(), 'walmart': load_walmart()}; \ + assert all(df.attrs.get('source') != 'synthetic_fallback' for df in dfs.values()); \ + pins = {'prop99': '16c3ac1da351788817433fc890ec2f502a8bdfcb46cbc8d693653330e71d5a65', \ + 'walmart': '410885572143dceb9daa643a8097768f1bc3493f9437451a9e4d1d5dc1e18d14'}; \ + cache = pathlib.Path.home() / '.cache' / 'diff_diff' / 'datasets'; \ + [1/0 for n, w in pins.items() if hashlib.sha256((cache / (n + '.dta')).read_bytes()).hexdigest() != w]" +# 2. generate (about 15 minutes; six B=9,999 bootstrap runs dominate) +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_lwdid_golden.do +grep -E '^r\([0-9]+\);' generate_lwdid_golden.log # must print nothing +``` + +## JSON schema + +```json +{ + "meta": {"ssc_versions": {"lwdid": "version 2.4.2 ..."}, "bootstrap_scheme": "...", + "control_pool": {...}, "datasets": {...}, "stata_version": 19, + "rireps": 100000, "riseed": ..., "bootstrap_reps": 9999, "bootstrap_seed": ...}, + "prop99": {"demean": {"att": ..., "se": ..., "p_ri": ...}, "detrend": {...}}, + "castle": {"demean": {"att": ..., "se": ...}, "detrend": {...}}, + "walmart": {"___": { + "watt": {"-22": [point, se], ..., "13": [point, se]}, + "overall": {"Pre_avg": [point, se], "Post_avg": [point, se]}}, ...} +} +``` + +Stata missing values are emitted as JSON `null` (the `_jnum` writer has an +explicit missing branch). + ## Known constraints - **Batch mode always exits 0**, even on a hard error (`r(NNN);`). Never trust the - shell exit code — parse the `.log` for `^r\([0-9]+\);`. 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`). + shell exit code — parse the `.log` for `^r\([0-9]+\);`. The LPDiD, ImputationDiD + and LWDiD generators also run informational in-`.do` smoke gates (LPDiD 1e-8, + ImputationDiD 1e-6, LWDiD 1e-10/1e-6 on its SMALL-N blocks; the LWDiD Walmart + block gates schema/cardinality/nonmissing only) that surface as `r(9);` on a + gross bug; the Python parity tests are authoritative. - **`c(flavor)` misreports the edition** as `IC` on StataSE, and `c(edition)` is unreliable. The generator derives the edition from the `c(MP)` / `c(SE)` 0/1 flags, with BE by elimination (`c(BE)` is undefined and `cond()` evaluates all branches @@ -295,5 +372,6 @@ Two arms: current value. - **SSC has no version history.** `ssc install` always fetches latest and there is no lockfile / archive to pin against. The LPDiD arm is exempt (`teffects` is native, - pinned by `version 19`); the ImputationDiD arm records its SSC package versions in - `meta.ssc_versions` so drift is at least detectable — new SSC arms should do the same. + pinned by `version 19`); the SSC arms (ImputationDiD, ETWFE/CS, LWDiD) record their + SSC package versions in `meta.ssc_versions` so drift is at least detectable — new + SSC arms should do the same. diff --git a/benchmarks/stata/generate_lwdid_golden.do b/benchmarks/stata/generate_lwdid_golden.do new file mode 100644 index 00000000..9d8ade00 --- /dev/null +++ b/benchmarks/stata/generate_lwdid_golden.do @@ -0,0 +1,334 @@ +*! Golden generator: LWDiD parity vs the authors' Stata `lwdid` package (SSC). +*! +*! Purpose +*! Produce the INDEPENDENT external anchor for the LWDiD estimator +*! (Lee & Wooldridge 2025, 2026) ahead of PR #588's merge: full-precision +*! small-N results on Prop 99 (incl. randomization inference at high reps) +*! and Castle Doctrine (tau_omega), plus the large-N Walmart event-study +*! WATT points and multiplier-bootstrap SEs for the six suite configs. +*! Printed-table goldens only reach three decimals; this arm anchors the +*! unprinted paths (RI convention, bootstrap SEs) at full precision. +*! +*! Consuming tests +*! tests/test_methodology_lwdid.py (gated on diff_diff.lwdid; activates +*! when PR #588 lands) and tests/test_lwdid_stata_golden_schema.py +*! (ungated schema/cardinality check that runs on main). +*! +*! Outputs (checked into the repo) +*! benchmarks/data/lwdid_stata_golden.json +*! +*! Usage (run from the repo root) +*! 1. FAIL-CLOSED warm-up - verifies the loader cache on disk (loader +*! success alone does not guarantee the cached bytes; cache writes are +*! best-effort): +*! python -c "from diff_diff import load_prop99, load_walmart; \ +*! import hashlib, pathlib; \ +*! dfs = {'prop99': load_prop99(), 'walmart': load_walmart()}; \ +*! assert all(df.attrs.get('source') != 'synthetic_fallback' for df in dfs.values()); \ +*! pins = {'prop99': '16c3ac1da351788817433fc890ec2f502a8bdfcb46cbc8d693653330e71d5a65', \ +*! 'walmart': '410885572143dceb9daa643a8097768f1bc3493f9437451a9e4d1d5dc1e18d14'}; \ +*! cache = pathlib.Path.home() / '.cache' / 'diff_diff' / 'datasets'; \ +*! [1/0 for n, w in pins.items() if hashlib.sha256((cache / (n + '.dta')).read_bytes()).hexdigest() != w]" +*! 2. /Applications/Stata/StataSE.app/Contents/MacOS/stata-se -b do \ +*! benchmarks/stata/generate_lwdid_golden.do +*! 3. Confirm the log is clean: grep -E '^r\([0-9]+\);' generate_lwdid_golden.log +*! (Stata batch mode ALWAYS exits 0, even on error - never trust the +*! exit code.) +*! +*! Notes +*! - `lwdid` is an SSC package (no version pinning possible); meta records +*! the verbatim `which lwdid` version line under `ssc_versions` so drift +*! is detectable, and the capture FAILS CLOSED if no version line is found. +*! - Small-N mode has NO control-group option: its composite regression is +*! never-treated-based by construction - empirically equivalent to the +*! Python `control_group='never_treated'` fits the parity tests use. +*! - Large-N default pool = never-treated + not-yet-treated (no `never` +*! option passed), matching the Python default +*! `control_group='not_yet_treated'`. +*! - In-.do smoke gates are informational early-abort checks against the +*! 2026-08-15/16 measured values; the Python tests are authoritative. + +version 19 +clear all +set more off +set type double + +* ------------------------------------------------------------------------------ +* Fail-closed dependency guard (SSC package; cannot be pinned by `version`). +* ------------------------------------------------------------------------------ +capture which lwdid +if _rc { + display as error "lwdid not installed - run benchmarks/stata/requirements.do" + exit 111 +} + +* ------------------------------------------------------------------------------ +* _jnum: format a scalar as a JSON number at round-trip-exact precision. +* Stata's %21.17g renders |x|<1 as ".455"/"-.455" (leading dot); JSON requires +* a leading 0, so patch ".x" -> "0.x" and "-.x" -> "-0.x". Stata MISSING is +* emitted as JSON null (the parity tables can carry missing cells). Returns r(s). +* ------------------------------------------------------------------------------ +capture program drop _jnum +program define _jnum, rclass + args x fmt + if "`fmt'" == "" local fmt "%21.17g" + if missing(`x') { + return local s "null" + exit + } + 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 + +* ------------------------------------------------------------------------------ +* _adover: capture an ado-file's version header for meta.ssc_versions. +* FAILS CLOSED: aborts if the captured line does not look like a version line +* (an "unknown" capture would silently disable drift detection). +* ------------------------------------------------------------------------------ +capture program drop _adover +program define _adover, rclass + args pkg + quietly findfile `pkg'.ado + local path "`r(fn)'" + tempname fh + file open `fh' using "`path'", read text + file read `fh' line + local ver "" + local tries 0 + while r(eof) == 0 & `tries' < 10 { + if strpos("`line'", "*!") == 1 & strpos(lower("`line'"), "version") > 0 { + local ver = strtrim(substr("`line'", 3, .)) + continue, break + } + file read `fh' line + local tries = `tries' + 1 + } + file close `fh' + if strpos(lower("`ver'"), "version") == 0 { + display as error "no version line found in `pkg'.ado - refusing to emit undetectable provenance" + exit 498 + } + return local s "`ver'" +end + +_adover lwdid +local lwdid_ver "`r(s)'" +display "lwdid version line: `lwdid_ver'" + +local RIREPS 100000 +local RISEED 20260815 +local BREPS 9999 +local BSEED 20260815 +local home : env HOME +local cache "`home'/.cache/diff_diff/datasets" + +* ============================================================================== +* Block 1: Prop 99 (small-N), demean + detrend, RI at high reps. +* ============================================================================== +use "`cache'/prop99.dta", clear +confirm numeric variable state year first_year lcigsale +quietly count +assert r(N) == 1209 + +foreach r in demean detrend { + lwdid lcigsale, small ivar(state) tvar(year) gvar(first_year) /// + rolling(`r') ri rireps(`RIREPS') riseed(`RISEED') + scalar p99_att_`r' = e(att) + scalar p99_se_`r' = e(se_att) + scalar p99_pri_`r' = e(p_ri) +} +* informational smoke gates (measured 2026-08-15) +assert reldif(p99_att_demean, -0.4221746150201265) < 1e-10 +assert reldif(p99_att_detrend, -0.2269886995561676) < 1e-10 +assert abs(p99_pri_detrend - 0.0508) < 0.005 + +* ============================================================================== +* Block 2: Castle Doctrine (small-N staggered), demean + detrend. +* ============================================================================== +import delimited using "benchmarks/data/real/castle_lw_subset.csv", clear case(preserve) +replace effyear = 0 if missing(effyear) +egen sid = group(state) +confirm numeric variable sid year effyear lhomicide +quietly count +assert r(N) == 550 + +foreach r in demean detrend { + lwdid lhomicide, small ivar(sid) tvar(year) gvar(effyear) rolling(`r') + scalar cas_att_`r' = e(att) + scalar cas_se_`r' = e(se_att) +} +* informational smoke gates (measured 2026-08-16 under `set type double`: +* `import delimited` reads the CSV as doubles here, unlike the exploratory +* smoke whose float import gave 0.09174538052... - an 8th-decimal float +* artifact; the double values below match the Python double-precision read) +assert reldif(cas_att_demean, 0.091745387139613596) < 1e-10 +assert reldif(cas_se_demean, 0.0571027) < 1e-6 +assert reldif(cas_att_detrend, 0.066550335128826035) < 1e-10 +assert reldif(cas_se_detrend, 0.0560124) < 1e-6 + +* ============================================================================== +* Block 3: Walmart (large-N), 3 configs x 2 outcomes, multiplier bootstrap. +* save() results go to tempfiles (auto-deleted; nothing written under the repo). +* Emitter runs AFTER all six fits so the panel-replacement hazard of reading a +* result .dta never touches a pending fit. +* ============================================================================== +use "`cache'/walmart.dta", clear +confirm numeric variable cid year first_year log_retail_emp log_wholesale_emp x1 x2 x3 +quietly count +assert r(N) == 29371 + +* config list: key|rolling|method|covars|outcome +local c1 "detrend_ra__log_retail_emp|detrend|ra||log_retail_emp" +local c2 "detrend_ipwra__log_retail_emp|detrend|ipwra|x1 x2 x3|log_retail_emp" +local c3 "demean_ipwra__log_retail_emp|demean|ipwra|x1 x2 x3|log_retail_emp" +local c4 "detrend_ra__log_wholesale_emp|detrend|ra||log_wholesale_emp" +local c5 "detrend_ipwra__log_wholesale_emp|detrend|ipwra|x1 x2 x3|log_wholesale_emp" +local c6 "demean_ipwra__log_wholesale_emp|demean|ipwra|x1 x2 x3|log_wholesale_emp" + +forvalues i = 1/6 { + * gettoken with parse("|") returns delimiter tokens; an EMPTY covars + * field surfaces as a second "|" token, handled below. + gettoken key rest : c`i', parse("|") + gettoken bar rest : rest, parse("|") + gettoken roll rest : rest, parse("|") + gettoken bar rest : rest, parse("|") + gettoken method rest : rest, parse("|") + gettoken bar rest : rest, parse("|") + gettoken covs rest : rest, parse("|") + if "`covs'" == "|" { + local covs "" + gettoken outc rest : rest, parse("|") + } + else { + gettoken bar rest : rest, parse("|") + gettoken outc rest : rest, parse("|") + } + display as txt "config `i': key=`key' rolling=`roll' method=`method' covs=[`covs'] outcome=`outc'" + + tempfile res`i' + set seed `BSEED' + preserve + lwdid `outc' `covs', ivar(cid) tvar(year) gvar(first_year) /// + rolling(`roll') method(`method') reps(`BREPS') save(`res`i'') + restore + local key`i' "`key'" +} + +* ------------------------------------------------------------------------------ +* Emit JSON. +* Compound-quote note: the `"' delimiter swallows a trailing double-quote, so +* every STRING field is written with a trailing comma; numeric meta fields have +* no trailing quote hazard. No timestamps (byte-identical regeneration). +* ------------------------------------------------------------------------------ +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/lwdid_stata_golden.json", write replace text +file write `fh' "{" _n +file write `fh' `" "meta": {"' _n +file write `fh' `" "estimator": "LWDiD (Lee & Wooldridge 2025, 2026) - authors' Stata lwdid package parity","' _n +file write `fh' `" "generator": "benchmarks/stata/generate_lwdid_golden.do","' _n +file write `fh' `" "ssc_versions": {"lwdid": "`lwdid_ver'"},"' _n +file write `fh' `" "stata_edition": "`sedition'","' _n +file write `fh' `" "bootstrap_scheme": "lwdid large-N multiplier bootstrap (package default; reps() draws with set seed); compared against diff-diff's unit-level Rademacher multiplier bootstrap at B=999","' _n +file write `fh' `" "control_pool": {"small_N": "composite regression, never-treated-based by construction (no control option exists); equivalent to Python control_group='never_treated'", "large_N": "never-treated + not-yet-treated default (the 'never' option is NOT passed); matches Python default control_group='not_yet_treated'"},"' _n +file write `fh' `" "datasets": {"prop99": {"url": "http://fmwww.bc.edu/repec/bocode/l/lw_smoking.dta", "sha256": "16c3ac1da351788817433fc890ec2f502a8bdfcb46cbc8d693653330e71d5a65"}, "walmart": {"url": "http://fmwww.bc.edu/repec/bocode/l/lw_walmart.dta", "sha256": "410885572143dceb9daa643a8097768f1bc3493f9437451a9e4d1d5dc1e18d14"}, "castle": {"path": "benchmarks/data/real/castle_lw_subset.csv"}},"' _n +file write `fh' `" "stata_version": `sver',"' _n +file write `fh' `" "rireps": `RIREPS',"' _n +file write `fh' `" "riseed": `RISEED',"' _n +file write `fh' `" "bootstrap_reps": `BREPS',"' _n +file write `fh' `" "bootstrap_seed": `BSEED'"' _n +file write `fh' " }," _n + +* --- prop99 block: {"demean": {"att":..,"se":..,"p_ri":..}, "detrend": {...}} +file write `fh' `" "prop99": {"' +local sep "" +foreach r in demean detrend { + _jnum p99_att_`r' + local a = r(s) + _jnum p99_se_`r' + local s = r(s) + _jnum p99_pri_`r' + local p = r(s) + file write `fh' "`sep'" _n `" "`r'": {"att": `a', "se": `s', "p_ri": `p'}"' + local sep "," +} +file write `fh' _n " }," _n + +* --- castle block: {"demean": {"att":..,"se":..}, "detrend": {...}} +file write `fh' `" "castle": {"' +local sep "" +foreach r in demean detrend { + _jnum cas_att_`r' + local a = r(s) + _jnum cas_se_`r' + local s = r(s) + file write `fh' "`sep'" _n `" "`r'": {"att": `a', "se": `s'}"' + local sep "," +} +file write `fh' _n " }," _n + +* --- walmart block: one key per (config, outcome); wrapper {"watt": {...}, "overall": {...}} +file write `fh' `" "walmart": {"' +local csep "" +forvalues i = 1/6 { + preserve + use "`res`i''", clear + * schema gate (from the save() spike): effect ryear watt se + aggregates + confirm string variable effect + confirm numeric variable ryear watt se + quietly count + local nrows = r(N) + quietly count if !missing(ryear) + local nwatt = r(N) + assert `nrows' == `nwatt' + 2 // exactly Pre_avg + Post_avg extra + * fail closed on incomplete regeneration: every emitted cell nonmissing, + * exactly one row per event time and per aggregate label + assert !missing(watt) & !missing(se) + quietly duplicates report ryear if !missing(ryear) + assert r(unique_value) == r(N) + file write `fh' "`csep'" _n `" "`key`i''": {"' _n + file write `fh' `" "watt": {"' + local rsep "" + quietly levelsof ryear if !missing(ryear), local(rs) + foreach rv of local rs { + quietly summarize watt if ryear == `rv', meanonly + scalar w_pt = r(mean) + quietly summarize se if ryear == `rv', meanonly + scalar w_se = r(mean) + _jnum w_pt + local a = r(s) + _jnum w_se + local s = r(s) + local rint = int(`rv') + file write `fh' "`rsep'" _n `" "`rint'": [`a', `s']"' + local rsep "," + } + file write `fh' _n " }," _n + file write `fh' `" "overall": {"' + local osep "" + foreach agg in Pre_avg Post_avg { + quietly summarize watt if effect == "`agg'", meanonly + scalar o_pt = r(mean) + quietly summarize se if effect == "`agg'", meanonly + scalar o_se = r(mean) + _jnum o_pt + local a = r(s) + _jnum o_se + local s = r(s) + file write `fh' "`osep'" _n `" "`agg'": [`a', `s']"' + local osep "," + } + file write `fh' _n " }" _n + file write `fh' " }" + local csep "," + restore +} +file write `fh' _n " }" _n +file write `fh' "}" _n +file close `fh' + +display "Wrote benchmarks/data/lwdid_stata_golden.json" diff --git a/benchmarks/stata/requirements.do b/benchmarks/stata/requirements.do index 8fd64d7a..3e436c1e 100644 --- a/benchmarks/stata/requirements.do +++ b/benchmarks/stata/requirements.do @@ -23,9 +23,11 @@ *! jwdid dependency; jwdid errors with "You need to install *! hdfe from SSC" without it) *! (both consumed by generate_etwfe_cs_golden.do) +*! lwdid - Lee & Wooldridge rolling DiD, the authors' reference +*! implementation (consumed by generate_lwdid_golden.do) version 19 -foreach p in ftools require reghdfe did_imputation drdid csdid hdfe jwdid { +foreach p in ftools require reghdfe did_imputation drdid csdid hdfe jwdid lwdid { capture which `p' if _rc { di as txt "Installing `p' from SSC ..." diff --git a/docs/benchmarks.rst b/docs/benchmarks.rst index 378ab14f..dd08bfc2 100644 --- a/docs/benchmarks.rst +++ b/docs/benchmarks.rst @@ -1,5 +1,5 @@ .. meta:: - :description: Validation benchmarks comparing diff-diff against R packages (did, synthdid, fixest) and Stata (teffects, did_imputation). 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, jwdid/csdid, reghdfe, lwdid). Coefficient accuracy, standard error comparison, and performance metrics. :keywords: difference-in-differences benchmark, DiD validation R, DiD validation Stata, python econometrics accuracy, did package comparison Benchmarks @@ -49,6 +49,15 @@ where no runnable R reference exists: * - ``ImputationDiD`` (leave-one-out SE) - Stata ``did_imputation, leaveout`` - Borusyak, Jaravel & Spiess (2024) Supp. App. A.9; no runnable R analogue + * - ``ETWFE`` / ``CallawaySantAnna`` (cross-check) + - Stata ``jwdid`` / ``csdid`` + - Wooldridge (2021) ETWFE and CS via the authors' Stata implementations + * - Clustered CR1 ``K_reference`` (shared linalg) + - Stata ``reghdfe`` + - Disconnected-panel absorbed-FE rank convention + * - ``LWDiD`` (validation arm precommitted; estimator in review, PR #588) + - Stata ``lwdid`` (authors' package) + - Lee & Wooldridge (2025, 2026); small-N exact + RI + event-study bootstrap Methodology ----------- @@ -1014,10 +1023,13 @@ Prerequisites pip install -e ".[dev]" 4. (Optional) Stata, only to regenerate the committed Stata goldens (``LPDiD`` - regression-adjustment SE, ``ImputationDiD`` leave-one-out SE). The 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 + regression-adjustment SE, ``ImputationDiD`` leave-one-out SE, the ETWFE/CS + cross-check, the ``reghdfe`` K_reference convention, and the LWDiD + authors'-package parity). The goldens are committed, so this is not needed to + run the test suite. The ``LPDiD`` arm uses the **native** ``teffects`` command; + the other arms depend on SSC packages + (``did_imputation``/``reghdfe``/``ftools``/``require``, + ``drdid``/``csdid``/``jwdid``/``hdfe``, ``lwdid``) — install them once via ``benchmarks/stata/requirements.do`` (the generators do not auto-install): .. code-block:: bash @@ -1027,6 +1039,9 @@ Prerequisites $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 + $STATA -b do benchmarks/stata/generate_etwfe_cs_golden.do + $STATA -b do benchmarks/stata/generate_reghdfe_kref_golden.do + $STATA -b do benchmarks/stata/generate_lwdid_golden.do # warm-up step: see benchmarks/stata/README.md Running Benchmarks ~~~~~~~~~~~~~~~~~~ diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index c4c9154c..917dbf84 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -2522,8 +2522,9 @@ Exact reviewed artifacts (PDF SHA-256) and live-verified SSRN metadata are pinne Full maintainer paper reviews (equation-level detail, replication targets): `docs/methodology/papers/lee-wooldridge-2025-review.md`, `docs/methodology/papers/lee-wooldridge-2026-review.md`. -- **Note:** Registry entry authored with the paper reviews ahead of the implementation (PR #588, third-party contribution under maintainer revision). Checklist boxes are unchecked until the implementation lands; implementation-specific edge-case notes will be finalized in that PR. -- **Note:** Maintainer validation suite: `tests/test_methodology_lwdid.py` (import-skip-gated until `diff_diff.lwdid` exists; xfail markers encode PR #588's outstanding acceptance criteria; replication goldens in `benchmarks/data/lwdid_walmart_eventstudy_golden.json` and `benchmarks/data/real/castle_lw_subset.csv`). +- **Note:** Registry entry authored with the paper reviews ahead of the implementation (PR #588, third-party contribution under maintainer revision). Checklist boxes are unchecked until the implementation lands; the implementation-specific Notes in this section (control pools, RI convention, overall conventions, Sec 4.3 rejection, API conformance) are the maintainer pre-pass of that finalization — remaining edge-case notes close with the merge. +- **Note:** Maintainer validation suite: `tests/test_methodology_lwdid.py` (import-skip-gated until `diff_diff.lwdid` exists) is the ACCEPTANCE BAR for PR #588's final round: it is written against the agreed post-rename API and intentionally fails on the pre-rename head; the round is complete when it passes unmodified. The earlier xfail scaffolding is retired. Goldens: `benchmarks/data/lwdid_walmart_eventstudy_golden.json` (printed Tables A4/A5; point halves consumed, printed SE halves superseded), `benchmarks/data/real/castle_lw_subset.csv`, and `benchmarks/data/lwdid_stata_golden.json` (full-precision authors'-package parity, `benchmarks/stata/generate_lwdid_golden.do`; schema checked ungated by `tests/test_lwdid_stata_golden_schema.py`). +- **Note:** Pre-release API-conformance decisions (locked 2026-08-15/16; no deprecation-ledger rows are owed because LWDiD never shipped on main): `estimation_method=` with house-CS values `reg`/`ipw`/`dr` plus the LWDiD-only `psm`; `vcov_type=` (name locked; value set finalized in the contributor's round within the house `linalg.py` vocabulary — `hc3` retained as the paper-mandated extra via a shared linalg extension, the `"cluster"` mode value retired in favor of clustering-by-`cluster=`, `hc0`/`hc4`/`conley` not offered); `seed=` with default `None`; `pscore_trim=`; `cluster=` on the constructor; the `LW` alias, module-level `lwdid()` wrapper, `lwdid_trend_diagnostics` module, and `overall_att`/`period_effects`/`period_specific` surfaces retired (per-period effects live on the unified post-fit event-study surface). **Key implementation requirements:** @@ -2546,20 +2547,28 @@ Event-study/placebo transformations over ALL periods (Appendix D): demeaning (D. *Estimation (per (g, t) cell; LW 2025 eq. (E.1)):* - Cross-sectional RA regression on the cell sample `A_{g,t}`: `Y_dot on 1, D_g, X, D_g (X - Xbar_g)`; ATT(g,t) = coefficient on `D_g`. No-covariate case reduces to plain DiD (eq. (3.4)); Theorem 3.1: common-timing per-period regressions are numerically equivalent to pooled OLS (3.6) (r = g reproduces the ETWFE estimand; `r > g` does not). - IPWRA (workhorse): logit propensity score per cell + WLS with weights `w = D + (1-D) p/(1-p)`; IPW = special case without the outcome-regression component. +- **Note:** API vocabulary (as specified for the PR #588 implementation, in review): `estimation_method="reg"` = the paper's RA, `"ipw"` = IPW, `"dr"` = IPWRA (house CallawaySantAnna values; the doubly-robust option), `"psm"` = propensity-score matching, an LWDiD-only extra with no CS analog. - Control pool at (g, r): `A_{r+1} = 1` (never-treated + not-yet-treated) by default; NT-only optional. Pre-treatment placebo cells use the Appendix D.3 rule `A_{g,t} = {G = g} ∪ {G = 0} ∪ {G > max(g,t)}`. *Aggregation:* - Event-study: `WATT(r) = sum_{g in G_r} omega_{g,r} ATT(g, g+r)` with `omega_{g,r}` = (treated units of cohort g contributing at event time r) / (total treated units contributing at event time r) - the operative definition per LW 2025 Appendix E.1, required under unbalanced panels where a cohort's contributing count at r can differ from `N_g`. In balanced panels this simplifies to `N_g / N_{G_r}` (Sec. 6.2/D.3). Aggregated influence function `IF_{i,r} = sum_g omega_{g,r} IF_{i,g,g+r}`. - Overall: composite-outcome single regression (LW 2026 eqs. (7.18)-(7.19)) — `tau_omega` with cohort-share weights `omega_g = N_g / N_treat`; automatically accounts for correlation among per-cohort effects and supports exact small-N inference. +- **Note:** The authors use TWO overall conventions across their own modes (measured 2026-08-15 against Stata `lwdid` v2.4.2): the small-N composite regression targets `tau_omega` (cohort-mean-then-treated-weight, eq. (7.18)) while the large-N display's `Post_avg` is the cell-mass (treated-count-per-cell) weighted average of the post ATT(g,t) — which equals the house CallawaySantAnna `"simple"` convention ON BALANCED PANELS (verified to the last digit on Walmart; CS-simple prefers fixed unit-cohort mass under unbalanced panels, so the equivalence carries that qualification). The implementation's `.att` is `tau_omega` — what the papers print, and the estimand the composite regression's inference is built for; a `vcov_type` selection must never move the point estimate (cross-path pins in the validation suite, including an unbalanced-panel pin). `aggregate("simple")` view-relays the fit per the house aggregation contract; exposing the cell-mass overall as an aggregate extra is a TODO row. *Standard errors:* -- Large-N default: influence-function **multiplier bootstrap** (LW 2025 Algorithm 1): IFs per (g,t) from E.2 (RA, finite-sample exact), E.3 (IPWRA, stacked M-estimator with logit-score correction), E.4 (IPW, `psi - Gamma' IF_gamma` correction); centered IFs; **unit-level Rademacher multipliers** (one draw per unit across all cells — unit clustering by construction); sup-t simultaneous bands over the event-study path; B = 999 in the paper's application; anchor periods excluded. +- Large-N default: influence-function **multiplier bootstrap** (LW 2025 Algorithm 1): IFs per (g,t) from E.2 (RA, finite-sample exact), E.3 (IPWRA, stacked M-estimator with logit-score correction), E.4 (IPW, `psi - Gamma' IF_gamma` correction); centered IFs; **unit-level Rademacher multipliers** (one draw per unit across all cells — unit clustering by construction); sup-t simultaneous bands over the event-study path; B = 999 in the paper's application; anchor periods excluded from the raw effect set (the public `EventStudyResults` surface still emits OBSERVED anchors as zero-valued `is_reference` rows — that is a display convention, not a contradiction of the exclusion). +- **Note:** Inference DEFAULT (as specified for the PR #588 implementation, in review): analytical vcov with the bootstrap as opt-in (`n_bootstrap=999` for the paper's setting) — the house convention (CallawaySantAnna ships the same analytical default), while this bullet describes the paper's recommended large-N inference. - Small-N exact (LW 2026): usual OLS SE on the collapsed cross-sectional regression with exact `T_{N-2}` / `T_{N-K-2}` reference distribution; valid down to `N = 3` and a single treated unit (`N1 = 1` — the t statistic is the studentized residual; same for `N_g = 1` per cohort in (7.8)/(7.10)). +- **Note (controlled exact inference is design-coherent):** LW 2026 Section 2 sanctions two controlled collapsed regressions — plain controls `(1, D, X)` with exact reference `T_{N-K-2}`, and interacted controls `(1, D, X, D(X - Xbar_1))` whose design rank implies `T_{N-2K-2}` (with the per-group guards `N0 > K+1` and `N1 > K+1`). Either is admissible for `vcov_type="classical"`; mixing them (fitting one design while reporting the other design's df) yields anti-conservative exact p-values and is a defect. The validation suite pins coherence: the reported p-value must use the residual df of whichever design reproduces the fitted point estimate. - Alternatives: HC3 when there are "at least a handful" of treated units; randomization inference for the sharp null (two-sided p = c / #permutations; Stata `lwdid` `ri` option); higher-level clustering and Conley SHAC SEs for larger cross sections (LW 2026 Sec. 8.2, citing Abadie-Athey-Imbens-Wooldridge 2023). +- **Note (IPWRA variance forms, measured divergence):** the PR #588 implementation's IPW/IPWRA influence functions are AIPW/Lunceford-Davidian-style, NOT the papers' E.2-E.4 stacked forms that the authors' Stata package implements. Measured 2026-08-16 on the Walmart application: the IPWRA multiplier-bootstrap SEs diverge SYSTEMATICALLY from the package's (~15% at event-time level, far beyond Monte-Carlo bounds), while the RA config's SEs agree within the MC bound and IPWRA POINT estimates agree to ~1e-3 (logit-optimizer paths). The validation suite therefore gates bootstrap-SE parity on the RA config only; the Stata IPWRA SE columns are committed as provenance in `lwdid_stata_golden.json`, and the E.3-form adjudication (implement the stacked IF, or document-and-anchor the AIPW alternative) is a required item of the contribution's final round — the checklist's E.2/E.3/E.4 box stays unchecked until it resolves. +- **Note (RI convention):** the paper states `p = c / #permutations`, but the authors' own package (Stata `lwdid` v2.4.2, measured 2026-08-15) implements the INCLUSIVE Phipson-Smyth rule — Monte-Carlo shuffles of the treatment vector with `p = (#{|coef| >= |b0|} + 1) / (reps + 1)`, ties counted as extreme — converging to ~0.0508 on the Prop 99 detrend application at 100k reps. The paper's printed RI p = 0.020 is NOT reproducible with the package (~4.5 binomial SD away) and is recorded as an as-printed discrepancy; the implementation and the validation goldens follow the package convention. The paper reviews (`docs/methodology/papers/lee-wooldridge-2026-review.md`) remain paper-faithful and state the paper's c/N convention as printed. +- **Note:** Conley SHAC SEs (listed above as a paper alternative) are NOT offered by the implementation — LWDiD exposes no spatial-coordinate inputs; the vcov design keeps to the house `linalg.py` vocabulary. *Edge cases:* -- Anchor periods: event-study omits `r = -1` (demeaning) / `r = -2, -1` (detrending); bootstrap excludes them. +- Anchor periods: event-study omits `r = -1` (demeaning) / `r = -2, -1` (detrending); bootstrap excludes them. (Raw-effect exclusion; the public results surface emits observed anchors as `is_reference` rows — see the Standard errors note.) - All units eventually treated (LW 2025 Sec. 4.3): drop `D_infinity`; effects defined relative to the last cohort; no effect estimable for the last cohort. +- **Note:** The PR #588 implementation (in review) REJECTS all-eventually-treated designs with a ValueError under BOTH control strategies, via mutually exclusive branches on `control_group` (the default `not_yet_treated` strategy has an explicit guard; the optional `never_treated` strategy rejects via its fewer-than-two-never-treated-units guard). Sec 4.3's last-cohort-reference mode is deferred — see the DEFERRED.md row. - Unbalanced panels (Sec. 4.4): transform observed data; per-cell observability requirements (1 pre-period demeaning / 2 detrending + outcome observed at r); selection may correlate with unit heterogeneity (and trends, under detrending) but not with shocks. - Anticipation: drop periods just prior to the intervention from the pre-average/trend (LW 2025 Sec. 4.4; LW 2026 eq. (2.22) anchor `Ybar_{i,S0}`, `S0 < S-1`). - Periods with no newly-treated units: no effects estimated there. @@ -2575,11 +2584,11 @@ Event-study/placebo transformations over ALL periods (Appendix D): demeaning (D. - LW 2025 Algorithm 1 (multiplier bootstrap, 7 steps): compute per-cell IFs -> aggregate WATT(r) -> center IFs -> B unit-level Rademacher draws -> bootstrap SEs -> sup statistics -> sup-t bands. **Reference implementation(s):** -- Stata: user-written `lwdid` (Hur, Lee and Wooldridge 2026; SSC) — implements the full procedure, multiplier-bootstrap inference, and randomization inference (`ri` option); ancillary datasets `lw_smoking.dta` / `lw_walmart.dta` (MIT) are the sources for `load_prop99()` / `load_walmart()`. +- Stata: user-written `lwdid` (Hur, Lee and Wooldridge 2026; SSC) — implements the full procedure, multiplier-bootstrap inference, and randomization inference (`ri` option); ancillary datasets `lw_smoking.dta` / `lw_walmart.dta` (MIT) are the sources for `load_prop99()` / `load_walmart()`. Full-precision parity goldens against this package (v2.4.2 recorded in the golden meta): `benchmarks/stata/generate_lwdid_golden.do` -> `benchmarks/data/lwdid_stata_golden.json`. - R: none. **Replication targets (from the papers; datasets available via `diff_diff.datasets`):** -- Prop 99 (LW 2026 Table 3, 38-state donor pool): demeaning ATT = -0.422 (SE 0.121); detrending ATT = -0.227 (SE 0.094), exact p = 0.021 vs randomization-inference p = 0.020. +- Prop 99 (LW 2026 Table 3, 38-state donor pool): demeaning ATT = -0.422 (SE 0.121); detrending ATT = -0.227 (SE 0.094), exact p = 0.021 vs randomization-inference p = 0.020 as printed. **Note:** the printed RI p is anchored INSTEAD to the authors'-package value (~0.051) — see the RI-convention note under Standard errors. - Castle laws (LW 2026 Sec. 7): `tau_omega` = 0.092 (demeaning; OLS SE 0.057), 0.067 (detrending). - Walmart entry (LW 2025 Tables A4/A5, 1,277 counties): per-relative-period WATT(r) with SEs for r = 0..13. @@ -2593,9 +2602,10 @@ Event-study/placebo transformations over ALL periods (Appendix D): demeaning (D. - [ ] Algorithm 1 multiplier bootstrap: unit-level Rademacher, sup-t simultaneous bands - [ ] Composite-outcome overall aggregation (7.18)/(7.19) with cohort-share weights - [ ] Exact-t inference: `T_{N-2}` / `T_{N-K-2}`, valid to `N = 3`, `N1 = 1`, `N_g = 1`; sample-size guards enforced -- [ ] HC3 alternative; randomization inference (two-sided p = c / #permutations); higher-level clustering / SHAC for larger N +- [ ] HC3 alternative; randomization inference (paper convention p = c / #permutations — implemented per the authors'-package inclusive convention, see the RI note); higher-level clustering per the vcov design (SHAC/Conley not offered, see note) - [ ] Anticipation-robustness period dropping; seasonal dummies in the transformation step -- [ ] All-eventually-treated (Sec. 4.3) and unbalanced-panel (Sec. 4.4) support +- [ ] All-eventually-treated (Sec. 4.3) support — deferred by decision (implementation rejects; see the Edge cases note and the DEFERRED.md row) +- [ ] Unbalanced-panel (Sec. 4.4) support - [ ] Common-timing no-covariate case reproduces plain DiD (3.4); Theorem 3.1 pooled-OLS equivalence (cross-estimator test vs `DifferenceInDifferences` / ETWFE at r = g) - [ ] Prop 99 / castle-laws / Walmart replication targets pinned as tests diff --git a/tests/test_lwdid_stata_golden_schema.py b/tests/test_lwdid_stata_golden_schema.py new file mode 100644 index 00000000..49158978 --- /dev/null +++ b/tests/test_lwdid_stata_golden_schema.py @@ -0,0 +1,118 @@ +"""Schema/cardinality validation for the LWDiD Stata parity golden. + +This module is deliberately UNGATED (no ``diff_diff.lwdid`` import): the +methodology suite (``tests/test_methodology_lwdid.py``) skips entirely until +the LWDiD estimator lands, so without this file main's CI would never open +``benchmarks/data/lwdid_stata_golden.json`` and a malformed or regenerated +golden could merge unvalidated. + +It also HOSTS the expected key-set constants shared with the methodology +suite - the import must run in this direction (the gated module imports from +here), because importing the gated module raises ``Skipped`` at collection +on main. + +The golden is produced by ``benchmarks/stata/generate_lwdid_golden.do`` +against the authors' SSC ``lwdid`` package; see that file and +``benchmarks/stata/README.md`` for provenance and regeneration. +""" + +import json +from pathlib import Path + +import pytest + +STATA_GOLDEN_PATH = ( + Path(__file__).resolve().parent.parent / "benchmarks" / "data" / "lwdid_stata_golden.json" +) + +#: The six Walmart (config, outcome) table keys: "___" +#: (method in the authors' Stata vocabulary: ra / ipwra). +WALMART_CONFIG_KEYS = ( + "detrend_ra__log_retail_emp", + "detrend_ipwra__log_retail_emp", + "demean_ipwra__log_retail_emp", + "detrend_ra__log_wholesale_emp", + "detrend_ipwra__log_wholesale_emp", + "demean_ipwra__log_wholesale_emp", +) + +#: Measured per-config event-time label sets of the Stata WATT tables +#: (surface (b) of the cardinality pins): r = -22..13 INCLUDING the anchor +#: rows, identical across all six configs (measured from the generator run; +#: the walmart panel is 1977-1999 with cohorts 1986-1999, so the reachable +#: event-time range is data-determined and config-independent). +WALMART_WATT_LABELS = {key: tuple(range(-22, 14)) for key in WALMART_CONFIG_KEYS} + +#: Fixed meta keys the consuming tests read. +REQUIRED_META_KEYS = ( + "ssc_versions", + "bootstrap_scheme", + "control_pool", + "datasets", + "rireps", + "riseed", + "bootstrap_reps", + "bootstrap_seed", +) + + +def _load(): + if not STATA_GOLDEN_PATH.exists(): + # The Rust CI matrix copies only tests/ to a temp dir; skip (never + # error) when benchmarks/data is absent, like every Stata arm test. + pytest.skip(f"{STATA_GOLDEN_PATH.name} not committed (partial checkout)") + return json.loads(STATA_GOLDEN_PATH.read_text()) + + +def _is_finite_number(x): + import math + + return isinstance(x, (int, float)) and not isinstance(x, bool) and math.isfinite(x) + + +def test_top_level_blocks(): + golden = _load() + assert sorted(golden.keys()) == ["castle", "meta", "prop99", "walmart"] + + +def test_meta_keys_and_version_line(): + meta = _load()["meta"] + for key in REQUIRED_META_KEYS: + assert key in meta, key + # ssc drift detection depends on a real version line, not a fallback + assert "version" in meta["ssc_versions"]["lwdid"].lower() + assert isinstance(meta["bootstrap_reps"], int) and meta["bootstrap_reps"] >= 999 + assert isinstance(meta["rireps"], int) and meta["rireps"] >= 10000 + + +def test_prop99_block(): + block = _load()["prop99"] + assert sorted(block.keys()) == ["demean", "detrend"] + for rolling, cell in block.items(): + assert sorted(cell.keys()) == ["att", "p_ri", "se"], rolling + assert all(_is_finite_number(v) for v in cell.values()), rolling + + +def test_castle_block(): + block = _load()["castle"] + assert sorted(block.keys()) == ["demean", "detrend"] + for rolling, cell in block.items(): + assert sorted(cell.keys()) == ["att", "se"], rolling + assert all(_is_finite_number(v) for v in cell.values()), rolling + + +def test_walmart_block_cardinality(): + block = _load()["walmart"] + assert sorted(block.keys()) == sorted(WALMART_CONFIG_KEYS) + for key, table in block.items(): + assert sorted(table.keys()) == ["overall", "watt"], key + labels = sorted(int(r) for r in table["watt"]) + assert tuple(labels) == tuple(sorted(WALMART_WATT_LABELS[key])), key + for r, cell in table["watt"].items(): + assert isinstance(cell, list) and len(cell) == 2, f"{key} r={r}" + # fail closed: a regeneration with missing cells must not pass + assert all(_is_finite_number(v) for v in cell), f"{key} r={r}" + assert sorted(table["overall"].keys()) == ["Post_avg", "Pre_avg"], key + for agg, cell in table["overall"].items(): + assert isinstance(cell, list) and len(cell) == 2, f"{key} {agg}" + assert all(_is_finite_number(v) for v in cell), f"{key} {agg}" diff --git a/tests/test_methodology_lwdid.py b/tests/test_methodology_lwdid.py index d5ed9d07..5116738f 100644 --- a/tests/test_methodology_lwdid.py +++ b/tests/test_methodology_lwdid.py @@ -1,37 +1,43 @@ """Independent methodology-validation suite for the LWDiD estimator (PR #588). -This module is the maintainer-side acceptance suite for the third-party LWDiD +This module is the maintainer-side ACCEPTANCE BAR for the third-party LWDiD (Lee & Wooldridge rolling-transformation DiD) contribution. It is merged to main BEFORE the estimator exists: the module-level ``pytest.importorskip`` -makes it skip cleanly until ``diff_diff.lwdid`` lands, at which point every -test activates automatically on the estimator branch. +makes it skip cleanly until ``diff_diff.lwdid`` lands. + +The suite is written against the POST-RENAME canonical API agreed for the +final PR #588 round (``estimation_method=`` with values ``reg/ipw/dr/psm``, +``vcov_type=``, ``seed=``, ``pscore_trim=``, constructor ``cluster=``, +fit-time ``first_treat=``/``covariates=``, post-fit +``results.aggregate("event_study")``). It therefore INTENTIONALLY FAILS on +the pre-rename contribution head; the contributor's round is complete when +this suite passes unmodified. This module verifies that the LWDiD implementation matches: 1. The published replication targets of Lee & Wooldridge (2026), Tables 3, 4 and A1 (California Prop 99, three donor pools) including exact-inference - and randomization-inference p-values. + p-values, with randomization inference anchored to the authors' Stata + package convention (see the REGISTRY RI note: the paper's printed 0.020 + is not reproducible with the authors' package; the package's inclusive + Phipson-Smyth convention yields ~0.051 and is the anchored behavior). 2. The castle-doctrine staggered application of LW (2026) Section 7.2 - (tau_omega via the composite-outcome regression (7.18)/(7.19)). + (tau_omega via the composite-outcome regression (7.18)/(7.19)), including + cross-path estimand consistency on balanced AND unbalanced panels. 3. The event-study replication targets of Lee & Wooldridge (2025), Appendix F - Tables A4/A5 (Walmart entry), via the normative event-study API specified - in ``TestEventStudySpec`` (xfail until PR #588's Appendix D work lands). + Tables A4/A5 (Walmart entry) for point estimates, with bootstrap SEs + anchored to the full-precision Stata golden + (``benchmarks/data/lwdid_stata_golden.json``) at a Monte-Carlo-informed + tolerance. 4. Estimator-independent properties: translation invariance of SEs, - cross-estimator equivalences (plain DiD; per-period panel-DiD identity of - LW 2026 eq. (2.20); the T=3 detrending closed form of LW 2025 eq. (5.7)), - from-scratch reference implementations of Procedures 2.1/3.1 and the - staggered per-cohort demeaning, exact small-sample t inference - (T_{N-2} / T_{N-K-2}), and the Monte Carlo bias ordering of LW 2026 - Section 5 under heterogeneous trends. -5. REGISTRY.md edge cases: minimum pre-treatment periods (>= 1 demeaning, - >= 2 detrending). - -xfail semantics (first use of xfail in this codebase): an ``xfail`` marker -encodes an agreed, outstanding work item of PR #588 - the marker's ``reason`` -names the item. ``strict=True`` markers MUST be removed by the commit that -fixes the item (the test then passing turns XPASS into a hard failure, -forcing explicit acceptance). ``strict=False`` is used only where numerical -fragility across environments is plausible. + cross-estimator equivalences, from-scratch reference implementations, + exact small-sample t inference (T_{N-2} / T_{N-K-2}), and the Monte Carlo + bias ordering of LW 2026 Section 5 under heterogeneous trends. +5. REGISTRY.md edge cases and input contracts: minimum pre-treatment + periods, all-eventually-treated rejection under both control strategies, + unit-constant covariate/cluster validation on every path, mixed + time-family rejection, and the retirement of the pre-release API surfaces + (``TestRetiredSurfaces``). Data sources: @@ -45,11 +51,19 @@ year, effyear, lhomicide, homicide, population). Originally committed because the ``load_castle_doctrine`` upstream source was dead; that source is now pinned and verified, but this subset stays as the pinned artifact - these goldens were captured against. Consolidating onto the loader is a - separate decision, not a consequence of the source being reachable again. -- Walmart event-study goldens: + these goldens were captured against. +- Walmart event-study PRINT goldens: ``benchmarks/data/lwdid_walmart_eventstudy_golden.json`` (Tables A4/A5, - provenance embedded in the file). + provenance embedded in the file). NOTE: since the Stata arm landed, only + the point halves of its ``[att, se]`` cells are consumed - the + three-decimal printed SEs are superseded by the full-precision seeded + Stata golden below and stay in the file as provenance. +- Stata parity golden: ``benchmarks/data/lwdid_stata_golden.json`` + (generated by ``benchmarks/stata/generate_lwdid_golden.do`` against the + authors' SSC ``lwdid`` package; version recorded in its meta). Expected + key-set constants shared with ``tests/test_lwdid_stata_golden_schema.py`` + (the ungated schema test), imported FROM there because this module's + importorskip would break the reverse direction. References: @@ -64,7 +78,9 @@ - docs/methodology/papers/lee-wooldridge-{2025,2026}-review.md. """ +import importlib import json +import warnings from pathlib import Path import numpy as np @@ -72,6 +88,12 @@ import pytest from scipy import stats +from tests.test_lwdid_stata_golden_schema import ( + STATA_GOLDEN_PATH, + WALMART_CONFIG_KEYS, + WALMART_WATT_LABELS, +) + pytest.importorskip( "diff_diff.lwdid", reason="LWDiD estimator not yet on main (arrives via PR #588)", @@ -96,7 +118,7 @@ "detrend": {1989: (-0.043, 0.059), 1995: (-0.282, 0.112), 2000: (-0.403, 0.152)}, } TABLE3_DETREND_EXACT_P = 0.021 -TABLE3_DETREND_RI_P = 0.020 +PROP99_TREAT_START = 1989 # g: first treated year (event time r = year - g) # Table 4 (Southern pool) and Table A1 (Midwestern pool): Average rows SOUTHERN_POOL = ["Alabama", "Arkansas", "Louisiana", "Mississippi"] @@ -112,6 +134,34 @@ # Printed-precision tolerance for 3-decimal published tables PRINTED_ATOL = 1e-3 +# Full-precision Stata parity tolerances: observed cross-BLAS agreement is +# ~1e-9..3e-9 ABSOLUTE while some parity cells are near zero, so an atol +# floor is mandatory (rtol alone would demand ~1e-11 on ~1e-5 cells). +PARITY_RTOL = 1e-6 +PARITY_ATOL = 1e-8 +# IPWRA configs go through a logit propensity fit whose optimizer path +# differs across implementations: measured cross-implementation agreement vs +# Stata is ~9e-4 relative on well-sized points (Post_avg ~2e-4) and ~2e-5 +# ABSOLUTE on near-zero placebo cells - consistent with the contribution's +# own documented "within 1%" vs lwdid-py - while the RA config agrees to +# ~1e-9. The IPWRA gates therefore use documented looser tolerances (~5x +# headroom over the measured divergence, still far below any structural +# break). +IPWRA_PARITY_RTOL = 5e-3 +IPWRA_PARITY_ATOL = 1e-4 + +#: Stata golden keys use the authors' method vocabulary (ra/ipwra); the +#: canonical API uses reg/dr. ONE decoder so no test drifts. +STATA_TO_CANONICAL_METHOD = {"ra": "reg", "ipwra": "dr"} + + +def _decode_walmart_key(key): + """ "___" -> (rolling, canonical_method, outcome).""" + rolling, method_outcome = key.split("_", 1) + stata_method, outcome = method_outcome.split("__") + return rolling, STATA_TO_CANONICAL_METHOD[stata_method], outcome + + _CASTLE_CSV = ( Path(__file__).resolve().parent.parent / "benchmarks" / "data" / "real" / "castle_lw_subset.csv" ) @@ -122,33 +172,18 @@ / "lwdid_walmart_eventstudy_golden.json" ) -XFAIL_IPW_CENTERING = pytest.mark.xfail( - strict=True, - reason="PR #588 step-2 item 1: IPW influence function is un-centered, " - "making the IPW SE translation-variant. Remove this marker in the " - "commit that centers the IPW IF.", -) -XFAIL_EVENT_STUDY = pytest.mark.xfail( - strict=True, - reason="PR #588 Option A: Appendix D event study + Algorithm 1 " - "multiplier bootstrap not yet implemented. Remove this marker in the " - "commit that implements the event study (deterministic spec tests).", -) -XFAIL_EVENT_STUDY_GOLDENS = pytest.mark.xfail( - strict=False, - reason="PR #588 Option A: Appendix D event study + Algorithm 1 not yet " - "implemented. Non-strict (numerical fragility): the golden SEs are the " - "paper's printed B=999 multiplier-bootstrap draws; a re-seeded bootstrap " - "can sit near the printed-precision tolerance boundary across " - "platforms. Re-calibrate the SE tolerance when the event study lands, " - "then remove the marker.", -) - # --------------------------------------------------------------------------- # Helpers # --------------------------------------------------------------------------- +def _load_stata_golden(): + """Load the Stata parity golden or skip (partial checkout).""" + if not STATA_GOLDEN_PATH.exists(): + pytest.skip(f"{STATA_GOLDEN_PATH.name} not committed (partial checkout)") + return json.loads(STATA_GOLDEN_PATH.read_text()) + + def _real_prop99(): """Load the real Prop 99 panel or skip (offline / synthetic fallback). @@ -172,7 +207,7 @@ def _real_walmart(): def _fit_prop99(df, rolling, **kwargs): - est = LWDiD(rolling=rolling, estimator="ra", vce="classical", **kwargs) + est = LWDiD(rolling=rolling, estimation_method="reg", vcov_type="classical", **kwargs) return est.fit(df, outcome="lcigsale", unit="state", time="year", treatment="treated") @@ -273,19 +308,100 @@ def _synthetic_staggered(n_units=120, t_max=10, cohorts=(5, 7), nt_share=0.4, se return pd.DataFrame(rows) +def _synthetic_unbalanced_staggered(): + """Deterministic UNBALANCED never-treated staggered panel. + + Exact recipe of the 2026-08-16 verification run that demonstrated the + composite-vs-joint-IF point divergence (2.1318 vs 2.1136) on the + contribution head: draws happen inside the unit/time loops in this exact + order (alpha per unit first, then one noise draw per RETAINED row - the + unbalancing `continue` fires BEFORE the row's noise draw, so the stream + depends on the drop pattern exactly as below). + """ + rng = np.random.default_rng(7) + rows = [] + uid = 0 + for g, n in [(0, 8), (4, 5), (6, 5)]: + for _ in range(n): + alpha = rng.normal() + for t in range(1, 9): + if uid % 3 == 1 and t in (2, 8): + continue + d = int(g > 0 and t >= g) + y = alpha + 0.1 * t + rng.normal(scale=0.3) + (2.0 + 0.5 * (g == 6)) * d + rows.append(dict(unit=uid, time=t, first=g, treat=d, y=y)) + uid += 1 + return pd.DataFrame(rows) + + # --------------------------------------------------------------------------- # 1. Prop 99, Table 3 (38-state donor pool) # --------------------------------------------------------------------------- +@pytest.fixture(scope="module") +def prop99(): + return _real_prop99() + + +@pytest.fixture(scope="module") +def walmart(): + return _real_walmart() + + +@pytest.fixture(scope="module") +def print_golden(): + """The printed Tables A4/A5 golden (point halves consumed only).""" + if not _WALMART_ES_GOLDEN.exists(): + pytest.skip(f"{_WALMART_ES_GOLDEN.name} not committed (partial checkout)") + return json.loads(_WALMART_ES_GOLDEN.read_text()) + + +@pytest.fixture(scope="module") +def stata_golden(): + return _load_stata_golden() + + +_ES_FIT_CACHE = {} + + +@pytest.fixture(scope="module") +def es_fit(walmart): + """Memoizing factory: each Walmart event-study config fits exactly ONCE + across all consuming classes (B=999 bootstrap fits are the expensive + unit of this module).""" + + def _fit(rolling, estimation_method, outcome="log_retail_emp"): + key = (rolling, estimation_method, outcome) + if key not in _ES_FIT_CACHE: + # The dr/ipw configs REQUIRE covariates: without them the paths + # degrade to reg (with a warning), which would invalidate the + # Stata parity comparison (the Stata runs pass x1-x3). + covariates = ["x1", "x2", "x3"] if estimation_method in ("dr", "ipw") else None + est = LWDiD( + rolling=rolling, + estimation_method=estimation_method, + n_bootstrap=999, + seed=42, + ) + _ES_FIT_CACHE[key] = est.fit( + walmart, + outcome=outcome, + unit="cid", + time="year", + treatment="treated", + first_treat="first_year", + covariates=covariates, + ) + return _ES_FIT_CACHE[key] + + return _fit + + @pytest.mark.realdata class TestProp99Table3Goldens: """LW (2026) Table 3: the authors' Stata `lwdid` output, frozen in print.""" - @pytest.fixture(scope="class") - def prop99(self): - return _real_prop99() - @pytest.mark.parametrize("rolling", ["demean", "detrend"]) def test_average_att_and_se(self, prop99, rolling): res = _fit_prop99(prop99, rolling) @@ -295,38 +411,66 @@ def test_average_att_and_se(self, prop99, rolling): @pytest.mark.parametrize("rolling", ["demean", "detrend"]) def test_per_period_atts(self, prop99, rolling): - res = _fit_prop99(prop99, rolling, period_specific=True) - assert res.period_effects, "period_specific=True should populate period_effects" + """Common-timing per-period effects via the unified post-fit + event-study surface (the `period_specific=`/`period_effects` pair is + retired). Event time r = year - g with g = 1989; Table 3 prints both + the point and the usual OLS SE per year.""" + res = _fit_prop99(prop99, rolling) + es = res.aggregate("event_study") + labels = {int(e): i for i, e in enumerate(np.asarray(es.event_time))} for year, (att, se) in TABLE3_PER_PERIOD[rolling].items(): - eff = res.period_effects[year] - np.testing.assert_allclose(eff["att"], att, atol=PRINTED_ATOL) - np.testing.assert_allclose(eff["se"], se, atol=PRINTED_ATOL) + r = year - PROP99_TREAT_START + assert r in labels, f"event time {r} (year {year}) missing from surface" + i = labels[r] + assert not bool(np.asarray(es.is_reference)[i]) + np.testing.assert_allclose(np.asarray(es.att)[i], att, atol=PRINTED_ATOL) + np.testing.assert_allclose(np.asarray(es.se)[i], se, atol=PRINTED_ATOL) def test_detrend_exact_inference_p_value(self, prop99): res = _fit_prop99(prop99, "detrend") np.testing.assert_allclose(res.p_value, TABLE3_DETREND_EXACT_P, atol=PRINTED_ATOL) - @pytest.mark.xfail( - strict=False, - reason="PR #588 step-2 discussion: RI p-value convention diverges " - "from LW 2026 Table 3 Note 2 (implementation gives the seed-stable " - "~2/39 two-sided exact permutation atom for N1=1 among 39 states - " - "arguably the standard exact answer - vs the paper's 0.020, whose " - "permutation scheme is under-documented; see the maintainer review " - "doc Gaps section). Reconcile against the authors' Stata " - "`lwdid, ri` behavior.", - ) - def test_detrend_randomization_inference_p_value(self, prop99): + def test_detrend_randomization_inference_p_value(self, prop99, stata_golden): + """RI anchored to the AUTHORS' PACKAGE convention, not the paper's + printed 0.020. + + Adjudicated 2026-08-15 against Stata `lwdid` v2.4.2: its RI counts + permutations at least as extreme (inclusive ties) under the + Phipson-Smyth plus-one rule, converging to ~0.0508 at 100k reps - + the paper's printed 0.020 is NOT reproducible (~4.5 binomial SD + away) and is recorded as an as-printed discrepancy in the REGISTRY + RI note. Tolerance covers ~3 binomial SEs at p ~= 0.05 with 1,000 + replications. + """ from diff_diff.lwdid_randomization import randomization_inference + anchor = stata_golden["prop99"]["detrend"]["p_ri"] ybar = _detrend_reference(prop99, "state", "year", "lcigsale", pre_end=1988) treated = prop99.groupby("state")["first_year"].first().loc[ybar.index] > 0 ri = randomization_inference( ybar.to_numpy(), treated.to_numpy(dtype=float), n_reps=1000, seed=2026 ) - # Paper: RI p = 0.020 from 1,000 permutations; tolerance covers ~3 - # binomial standard errors at p ~= 0.02 with 1,000 replications. - np.testing.assert_allclose(ri.pvalue, TABLE3_DETREND_RI_P, atol=0.015) + n_reps = 1000 + atol = 3 * np.sqrt(anchor * (1 - anchor) / n_reps) + np.testing.assert_allclose(ri.pvalue, anchor, atol=atol) + + def test_randomization_inference_plus_one_formula(self, prop99): + """Deterministic pin of the inclusive Phipson-Smyth formula: the + reported p equals (c + 1) / (B + 1) EXACTLY, with c the count of + permuted statistics at least as extreme as observed, recomputed here + from the returned draw distribution (the MC-tolerance anchor above + cannot distinguish (c+1)/(B+1) from c/B).""" + from diff_diff.lwdid_randomization import randomization_inference + + ybar = _detrend_reference(prop99, "state", "year", "lcigsale", pre_end=1988) + treated = prop99.groupby("state")["first_year"].first().loc[ybar.index] > 0 + ri = randomization_inference( + ybar.to_numpy(), treated.to_numpy(dtype=float), n_reps=499, seed=7 + ) + dist = np.asarray(ri.att_distribution, dtype=float) + assert dist.shape == (499,) + c = int(np.sum(np.abs(dist) >= abs(ri.att_observed))) + np.testing.assert_allclose(ri.pvalue, (c + 1) / (499 + 1), rtol=0, atol=0) # --------------------------------------------------------------------------- @@ -338,10 +482,6 @@ def test_detrend_randomization_inference_p_value(self, prop99): class TestDonorPoolVariants: """LW (2026) Tables 4 / A1: Southern and Midwestern donor pools.""" - @pytest.fixture(scope="class") - def prop99(self): - return _real_prop99() - @pytest.mark.parametrize( "pool,targets", [(SOUTHERN_POOL, TABLE4_AVERAGE), (MIDWEST_POOL, TABLEA1_AVERAGE)], @@ -362,6 +502,43 @@ def test_average_att_and_se(self, prop99, pool, targets, rolling): # --------------------------------------------------------------------------- +@pytest.fixture(scope="module") +def castle(): + if not _CASTLE_CSV.exists(): + pytest.skip(f"{_CASTLE_CSV.name} not committed (partial checkout)") + df = pd.read_csv(_CASTLE_CSV) + # Data-integrity assertion: this must be the paper's exact sample. + fy = df.groupby("state")["effyear"].first() + counts = fy.dropna().astype(int).value_counts().sort_index().to_dict() + assert counts == CASTLE_COHORTS, f"castle cohorts {counts} != paper {CASTLE_COHORTS}" + assert df["state"].nunique() == 50 + assert int(fy.isna().sum()) == 29 + df = df.copy() + df["first_year"] = df["effyear"].fillna(0).astype(int) + df["treat"] = ((df["first_year"] > 0) & (df["year"] >= df["first_year"])).astype(int) + return df + + +def _fit_castle(castle, rolling, vcov_type="classical"): + # The paper's per-cohort effects (7.10) use never-treated controls; + # PR #588's default control pool is not-yet-treated, which moves the + # castle point estimate from 0.092 to 0.074 (calibrated 2026-07-13). + est = LWDiD( + rolling=rolling, + estimation_method="reg", + vcov_type=vcov_type, + control_group="never_treated", + ) + return est.fit( + castle, + outcome="lhomicide", + unit="state", + time="year", + treatment="treat", + first_treat="first_year", + ) + + @pytest.mark.realdata class TestCastleTauOmegaAdjudicator: """LW (2026) Section 7.2: staggered tau_omega via composite regression. @@ -372,37 +549,6 @@ class TestCastleTauOmegaAdjudicator: the paper's point estimate and OLS SE simultaneously. """ - @pytest.fixture(scope="class") - def castle(self): - if not _CASTLE_CSV.exists(): - pytest.skip(f"{_CASTLE_CSV.name} not committed (partial checkout)") - df = pd.read_csv(_CASTLE_CSV) - # Data-integrity assertion: this must be the paper's exact sample. - fy = df.groupby("state")["effyear"].first() - counts = fy.dropna().astype(int).value_counts().sort_index().to_dict() - assert counts == CASTLE_COHORTS, f"castle cohorts {counts} != paper {CASTLE_COHORTS}" - assert df["state"].nunique() == 50 - assert int(fy.isna().sum()) == 29 - df = df.copy() - df["first_year"] = df["effyear"].fillna(0).astype(int) - df["treat"] = ((df["first_year"] > 0) & (df["year"] >= df["first_year"])).astype(int) - return df - - @staticmethod - def _fit(castle, rolling): - # The paper's per-cohort effects (7.10) use never-treated controls; - # PR #588's default control pool is not-yet-treated, which moves the - # castle point estimate from 0.092 to 0.074 (calibrated 2026-07-13). - est = LWDiD(rolling=rolling, estimator="ra", vce="classical", control_group="never_treated") - return est.fit( - castle, - outcome="lhomicide", - unit="state", - time="year", - treatment="treat", - cohort="first_year", - ) - @staticmethod def _composite_reference(castle, rolling): """From-scratch LW 2026 (7.18)/(7.19): composite outcome + single @@ -436,25 +582,58 @@ def test_composite_regression_reference_reproduces_paper(self, castle): np.testing.assert_allclose(tau_dt, CASTLE_TAU_DETREND, atol=PRINTED_ATOL) def test_demean_tau_omega_point(self, castle): - res = self._fit(castle, "demean") + res = _fit_castle(castle, "demean") np.testing.assert_allclose(res.att, CASTLE_TAU_DEMEAN[0], atol=PRINTED_ATOL) - @pytest.mark.xfail( - strict=True, - reason="PR #588 step-2 aggregation: the overall SE must come from the " - "composite-outcome regression (7.18)/(7.19) (paper OLS SE 0.057; " - "implementation's independence-across-cohorts SE gives 0.051). " - "Remove this marker in the commit that adopts the composite " - "regression.", - ) def test_demean_tau_omega_ols_se(self, castle): - res = self._fit(castle, "demean") + """The overall SE comes from the composite-outcome regression + (7.18)/(7.19) - the paper's OLS SE (0.057).""" + res = _fit_castle(castle, "demean") np.testing.assert_allclose(res.se, CASTLE_TAU_DEMEAN[1], atol=PRINTED_ATOL) def test_detrend_tau_omega_point(self, castle): - res = self._fit(castle, "detrend") + res = _fit_castle(castle, "detrend") np.testing.assert_allclose(res.att, CASTLE_TAU_DETREND, atol=PRINTED_ATOL) + def test_balanced_cross_path_estimand_consistency(self, castle): + """The composite path (classical gate) and the joint-IF path (hc1) + must estimate the SAME tau_omega: a variance option must never move + the point estimate.""" + rc = _fit_castle(castle, "demean", vcov_type="classical") + rh = _fit_castle(castle, "demean", vcov_type="hc1") + np.testing.assert_allclose(rc.att, rh.att, rtol=0, atol=1e-10) + + +class TestUnbalancedEstimandConsistency: + """Cross-path estimand consistency on an UNBALANCED panel. + + Verified by execution 2026-08-16 on the contribution head: the composite + path and the joint-IF path returned DIFFERENT points (2.1318 vs 2.1136, + diff 1.8e-2) on this exact fixture, switching purely on ``vcov_type`` - + the two weightings coincide under balance, so the balanced castle check + cannot see it. The joint-IF aggregation weights must implement the SAME + (7.18) tau_omega estimand as the composite regression. + """ + + def test_unbalanced_cross_path_estimand_consistency(self): + df = _synthetic_unbalanced_staggered() + kw = dict(outcome="y", unit="unit", time="time", treatment="treat", first_treat="first") + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + rc = LWDiD( + rolling="demean", + estimation_method="reg", + vcov_type="classical", + control_group="never_treated", + ).fit(df, **kw) + rh = LWDiD( + rolling="demean", + estimation_method="reg", + vcov_type="hc1", + control_group="never_treated", + ).fit(df, **kw) + np.testing.assert_allclose(rc.att, rh.att, rtol=0, atol=1e-10) + # --------------------------------------------------------------------------- # 4. Translation invariance of SEs (property test) @@ -467,28 +646,27 @@ class TestTranslationInvariance: SHIFT = 100.0 - def _fit_pair(self, estimator): + def _fit_pair(self, estimation_method): df = _synthetic_common_timing() df2 = df.copy() df2["y"] = df2["y"] + self.SHIFT * df2["post"] - kw = dict(outcome="y", unit="unit", time="time", treatment="treat", controls=["x"]) - r1 = LWDiD(rolling="demean", estimator=estimator, vce="hc1").fit(df, **kw) - r2 = LWDiD(rolling="demean", estimator=estimator, vce="hc1").fit(df2, **kw) + kw = dict(outcome="y", unit="unit", time="time", treatment="treat", covariates=["x"]) + r1 = LWDiD(rolling="demean", estimation_method=estimation_method, vcov_type="hc1").fit( + df, **kw + ) + r2 = LWDiD(rolling="demean", estimation_method=estimation_method, vcov_type="hc1").fit( + df2, **kw + ) return r1, r2 - @pytest.mark.parametrize("estimator", ["ra", "ipwra", "ipw"]) - def test_att_translation_invariant(self, estimator): - r1, r2 = self._fit_pair(estimator) + @pytest.mark.parametrize("estimation_method", ["reg", "dr", "ipw"]) + def test_att_translation_invariant(self, estimation_method): + r1, r2 = self._fit_pair(estimation_method) np.testing.assert_allclose(r1.att, r2.att, rtol=0, atol=1e-10) - @pytest.mark.parametrize("estimator", ["ra", "ipwra"]) - def test_se_translation_invariant(self, estimator): - r1, r2 = self._fit_pair(estimator) - np.testing.assert_allclose(r1.se, r2.se, rtol=0, atol=1e-10) - - @XFAIL_IPW_CENTERING - def test_ipw_se_translation_invariant(self): - r1, r2 = self._fit_pair("ipw") + @pytest.mark.parametrize("estimation_method", ["reg", "dr", "ipw"]) + def test_se_translation_invariant(self, estimation_method): + r1, r2 = self._fit_pair(estimation_method) np.testing.assert_allclose(r1.se, r2.se, rtol=0, atol=1e-10) @@ -500,28 +678,34 @@ def test_ipw_se_translation_invariant(self): class TestCrossEstimatorEquivalence: """Theory-mandated numerical identities against already-validated code.""" - def test_demean_ra_equals_plain_did(self): + def test_demean_reg_equals_plain_did(self): """Common timing, no covariates: rolling demeaning + RA reproduces the standard DiD estimator exactly (LW 2026, eq. (2.5) / Section 9).""" df = _synthetic_common_timing() - lw = LWDiD(rolling="demean", estimator="ra", vce="classical").fit( + lw = LWDiD(rolling="demean", estimation_method="reg", vcov_type="classical").fit( df, outcome="y", unit="unit", time="time", treatment="treat" ) dd = DifferenceInDifferences().fit(df, outcome="y", treatment="treated_group", post="post") np.testing.assert_allclose(lw.att, dd.att, rtol=0, atol=1e-10) - @pytest.mark.parametrize("r", [5, 7]) - def test_per_period_equals_subset_panel_did(self, r): + @pytest.mark.parametrize("r_calendar", [5, 7]) + def test_per_period_equals_subset_panel_did(self, r_calendar): """LW 2026 eq. (2.20): the per-period effect tau_hat_{t,DD} is numerically identical to a standard panel DiD run on periods - {1, ..., S-1, t}.""" + {1, ..., S-1, t}. Read through the unified event-study surface + (calendar period r_calendar with treatment start s=5 is event time + r_calendar - 5).""" df = _synthetic_common_timing(t_max=8, s=5) - res = LWDiD(rolling="demean", estimator="ra", vce="classical", period_specific=True).fit( + res = LWDiD(rolling="demean", estimation_method="reg", vcov_type="classical").fit( df, outcome="y", unit="unit", time="time", treatment="treat" ) - sub = df[(df["time"] < 5) | (df["time"] == r)] + es = res.aggregate("event_study") + labels = {int(e): i for i, e in enumerate(np.asarray(es.event_time))} + r_event = r_calendar - 5 + assert r_event in labels + sub = df[(df["time"] < 5) | (df["time"] == r_calendar)] dd = DifferenceInDifferences().fit(sub, outcome="y", treatment="treated_group", post="post") - np.testing.assert_allclose(res.period_effects[r]["att"], dd.att, rtol=0, atol=1e-10) + np.testing.assert_allclose(np.asarray(es.att)[labels[r_event]], dd.att, rtol=0, atol=1e-10) def test_detrend_t3_closed_form(self): """LW 2025 eq. (5.7): with T=3, S=3, no covariates, the detrending @@ -536,7 +720,7 @@ def test_detrend_t3_closed_form(self): y = a + b * t + rng.normal(0, 0.3) + (0.8 if treated and t == 3 else 0.0) rows.append({"unit": i, "time": t, "y": y, "treat": int(treated and t == 3)}) df = pd.DataFrame(rows) - res = LWDiD(rolling="detrend", estimator="ra", vce="classical").fit( + res = LWDiD(rolling="detrend", estimation_method="reg", vcov_type="classical").fit( df, outcome="y", unit="unit", time="time", treatment="treat" ) m = df.assign(g=(df["unit"] < 30).astype(int)).groupby(["g", "time"])["y"].mean() @@ -559,7 +743,7 @@ def test_procedure_2_1_demeaning(self): ybar = _demean_reference(df, "unit", "time", "y", pre_end=4) treated = df.groupby("unit")["treated_group"].first().loc[ybar.index] tau, se, dof = _cross_section_did(ybar, treated) - res = LWDiD(rolling="demean", estimator="ra", vce="classical").fit( + res = LWDiD(rolling="demean", estimation_method="reg", vcov_type="classical").fit( df, outcome="y", unit="unit", time="time", treatment="treat" ) np.testing.assert_allclose(res.att, tau, rtol=0, atol=1e-10) @@ -570,7 +754,7 @@ def test_procedure_3_1_detrending(self): ybar = _detrend_reference(df, "unit", "time", "y", pre_end=4) treated = df.groupby("unit")["treated_group"].first().loc[ybar.index] tau, se, dof = _cross_section_did(ybar, treated) - res = LWDiD(rolling="detrend", estimator="ra", vce="classical").fit( + res = LWDiD(rolling="detrend", estimation_method="reg", vcov_type="classical").fit( df, outcome="y", unit="unit", time="time", treatment="treat" ) np.testing.assert_allclose(res.att, tau, rtol=0, atol=1e-10) @@ -582,8 +766,13 @@ def test_staggered_per_cohort_demeaning(self): (control_group='never_treated' matches the (7.10) sample).""" df = _synthetic_staggered() res = LWDiD( - rolling="demean", estimator="ra", vce="classical", control_group="never_treated" - ).fit(df, outcome="y", unit="unit", time="time", treatment="treat", cohort="first_year") + rolling="demean", + estimation_method="reg", + vcov_type="classical", + control_group="never_treated", + ).fit( + df, outcome="y", unit="unit", time="time", treatment="treat", first_treat="first_year" + ) assert res.cohort_effects, "staggered fit should populate cohort_effects" fy = df.groupby("unit")["first_year"].first() for g in sorted(set(fy[fy > 0])): @@ -606,17 +795,55 @@ class TestExactSmallSampleInference: def test_classical_p_value_uses_t_n_minus_2(self): df = _synthetic_common_timing(n_treat=6, n_control=10, t_max=6, s=4) - res = LWDiD(rolling="demean", estimator="ra", vce="classical").fit( + res = LWDiD(rolling="demean", estimation_method="reg", vcov_type="classical").fit( df, outcome="y", unit="unit", time="time", treatment="treat" ) n = df["unit"].nunique() p_expected = 2 * stats.t.sf(abs(res.t_stat), n - 2) np.testing.assert_allclose(res.p_value, p_expected, rtol=1e-10) + def test_covariate_exact_inference_is_design_coherent(self): + """Covariate-adjusted exact inference must be DESIGN-COHERENT. + + LW 2026 Section 2 sanctions two controlled collapsed regressions: + the plain-controls design ``(1, D, X)`` with exact reference + T_{N-K-2}, and the interacted design ``(1, D, X, D(X - Xbar_1))`` + whose design rank implies T_{N-2K-2} (with the per-group sample + guards). Either is admissible; what is NOT is mixing them - fitting + one design while reporting the other design's df gives + anti-conservative exact p-values. MEASURED on the contribution head + 2026-08-16: the fitted att reproduces the INTERACTED design exactly + while the p-value uses the plain design's T_{N-K-2} (df 19 vs the + coherent 18 on a N=22, K=1 panel) - a required fix. This test + identifies which design produced the point estimate and requires the + p-value to use THAT design's residual df. + """ + df = _synthetic_common_timing(n_treat=8, n_control=14, t_max=6, s=4) + res = LWDiD(rolling="demean", estimation_method="reg", vcov_type="classical").fit( + df, outcome="y", unit="unit", time="time", treatment="treat", covariates=["x"] + ) + ybar = _demean_reference(df, "unit", "time", "y", pre_end=3) + d = df.groupby("unit")["treated_group"].first().loc[ybar.index].to_numpy(dtype=float) + xv = df.groupby("unit")["x"].first().loc[ybar.index].to_numpy(dtype=float) + y = ybar.to_numpy(dtype=float) + plain = np.column_stack([np.ones_like(d), d, xv]) + xbar_treated = xv[d == 1].mean() + interacted = np.column_stack([np.ones_like(d), d, xv, d * (xv - xbar_treated)]) + n = len(y) + matches = {} + for name, X in (("plain", plain), ("interacted", interacted)): + beta, *_ = np.linalg.lstsq(X, y, rcond=None) + matches[name] = (abs(float(beta[1]) - res.att) < 1e-10, X.shape[1]) + matched = [name for name, (hit, _) in matches.items() if hit] + assert len(matched) == 1, f"att matches designs: {matched or 'none'}" + n_params = matches[matched[0]][1] + p_expected = 2 * stats.t.sf(abs(res.t_stat), n - n_params) + np.testing.assert_allclose(res.p_value, p_expected, rtol=1e-10) + def test_single_treated_unit_inference_is_finite(self): """N1 = 1: exact inference remains valid (studentized-residual t).""" df = _synthetic_common_timing(n_treat=1, n_control=12, t_max=6, s=4) - res = LWDiD(rolling="demean", estimator="ra", vce="classical").fit( + res = LWDiD(rolling="demean", estimation_method="reg", vcov_type="classical").fit( df, outcome="y", unit="unit", time="time", treatment="treat" ) assert np.isfinite(res.att) @@ -626,73 +853,90 @@ def test_single_treated_unit_inference_is_finite(self): p_expected = 2 * stats.t.sf(abs(res.t_stat), n - 2) np.testing.assert_allclose(res.p_value, p_expected, rtol=1e-10) - @pytest.mark.xfail( - strict=True, - reason="PR #588 step-2: no N_infinity >= 2 guard exists for the " - "never-treated-only staggered control strategy (LW 2026, p26).", - ) def test_never_treated_pool_of_one_is_rejected(self): + """The N_infinity >= 2 requirement applies to the NT-only control + strategy (LW 2026, p26); NYT controls are exempt.""" df = _synthetic_staggered(n_units=30, nt_share=0.0, seed=5) # Force exactly one never-treated unit first_unit = df["unit"] == 0 df.loc[first_unit, "first_year"] = 0 df.loc[first_unit, "treat"] = 0 with pytest.raises(ValueError, match="[Nn]ever[- ]treated"): - # The N_infinity >= 2 requirement applies to the NT-only control - # strategy (LW 2026, p26); NYT controls are exempt. LWDiD( rolling="demean", - estimator="ra", - vce="classical", + estimation_method="reg", + vcov_type="classical", control_group="never_treated", - ).fit(df, outcome="y", unit="unit", time="time", treatment="treat", cohort="first_year") + ).fit( + df, + outcome="y", + unit="unit", + time="time", + treatment="treat", + first_treat="first_year", + ) + + @pytest.mark.parametrize("control_group", ["not_yet_treated", "never_treated"]) + def test_all_eventually_treated_is_rejected(self, control_group): + """All-eventually-treated designs are rejected under BOTH control + strategies (REGISTRY Sec 4.3 note; the last-cohort-reference mode is + deferred, see DEFERRED.md).""" + rows = [] + rng = np.random.default_rng(12) + for unit in range(12): + g = 3 if unit < 6 else 5 + for t in range(1, 7): + rows.append( + { + "unit": unit, + "time": t, + "first_year": g, + "treat": int(t >= g), + "y": float(rng.normal()), + } + ) + df = pd.DataFrame(rows) + with pytest.raises(ValueError): + LWDiD(rolling="demean", estimation_method="reg", control_group=control_group).fit( + df, + outcome="y", + unit="unit", + time="time", + treatment="treat", + first_treat="first_year", + ) # --------------------------------------------------------------------------- -# 8. Event-study specification (normative API; Option A work) +# 8. Event-study surface (post-fit aggregate; Walmart replication) # --------------------------------------------------------------------------- @pytest.mark.realdata class TestEventStudySpec: - """Normative event-study API for LWDiD (maintainer-specified). - - Invocation: ``fit(..., aggregate="event_study")`` (CallawaySantAnna - precedent). Results must expose ``event_study_effects: Dict[int, Dict]`` - keyed by relative period r, each with keys ``effect, se, t_stat, - p_value, conf_int`` and (when simultaneous bands are computed via - Algorithm 1) ``cband_conf_int``; result-level metadata ``cband_method, - cband_crit_value, cband_n_bootstrap``. Anchor periods are excluded: - r = -1 (demeaning); r = -2, -1 (detrending). - - All tests xfail until PR #588's Appendix D + Algorithm 1 work lands. + """Normative event-study surface for LWDiD (maintainer-specified). + + Per-relative-period effects are computed jointly at fit and read back + via post-fit ``results.aggregate("event_study")`` (fit-time + ``aggregate=`` is forbidden for new estimators, docs/v4-design.md). The + raw ``event_study_effects: Dict[int, Dict]`` dict excludes the anchor + periods (r = -1 demeaning; r = -2, -1 detrending); the public + ``EventStudyResults`` surface emits observed anchors as zero-valued + ``is_reference`` rows. Bootstrap SEs come from the Algorithm 1 + unit-level Rademacher multiplier bootstrap (B = 999 here). """ - @pytest.fixture(scope="class") - def walmart(self): - return _real_walmart() - - @pytest.fixture(scope="class") - def golden(self): - if not _WALMART_ES_GOLDEN.exists(): - pytest.skip(f"{_WALMART_ES_GOLDEN.name} not committed (partial checkout)") - return json.loads(_WALMART_ES_GOLDEN.read_text()) - - def _fit_es(self, walmart, rolling, estimator, outcome="log_retail_emp"): - # The golden SEs are Algorithm 1 multiplier-bootstrap SEs (B = 999): - # the spec requires the bootstrap path, not analytical vce. - est = LWDiD(rolling=rolling, estimator=estimator, n_bootstrap=999, bootstrap_seed=42) - return est.fit( - walmart, - outcome=outcome, - unit="cid", - time="year", - treatment="treated", - cohort="first_year", - aggregate="event_study", - ) + def test_walmart_raw_label_sets_match_measured_constants(self, es_fit): + """Fail-closed coverage: the raw event_study_effects key set equals + the measured per-config constant (shared with the schema test) so a + silently shrunken surface cannot pass by iterating fewer rows.""" + for key, expected in WALMART_WATT_LABELS.items(): + rolling, method, outcome = _decode_walmart_key(key) + res = es_fit(rolling, method, outcome=outcome) + anchors = {-1} if rolling == "demean" else {-2, -1} + expected_raw = set(expected) - anchors + assert set(res.event_study_effects) == expected_raw, key - @XFAIL_EVENT_STUDY @pytest.mark.parametrize( "outcome,table_key", [ @@ -702,63 +946,83 @@ def _fit_es(self, walmart, rolling, estimator, outcome="log_retail_emp"): ids=["a4_retail", "a5_wholesale"], ) @pytest.mark.parametrize( - "rolling,estimator,column", + "rolling,method,column", [ - ("detrend", "ra", "rolling_ra_detrend"), - ("detrend", "ipwra", "rolling_ipwra_detrend"), - ("demean", "ipwra", "rolling_ipwra_demean"), + ("detrend", "reg", "rolling_ra_detrend"), + ("detrend", "dr", "rolling_ipwra_detrend"), + ("demean", "dr", "rolling_ipwra_demean"), ], ) def test_walmart_eventstudy_point_goldens( - self, walmart, golden, rolling, estimator, column, outcome, table_key + self, print_golden, es_fit, rolling, method, column, outcome, table_key ): - """Deterministic WATT(r) point estimates vs Tables A4/A5 (strict).""" - res = self._fit_es(walmart, rolling, estimator, outcome=outcome) - table = golden[table_key] + """Deterministic WATT(r) point estimates vs printed Tables A4/A5.""" + res = es_fit(rolling, method, outcome=outcome) + table = print_golden[table_key] + # printed tables cover r = 0..13, a subset of the computed surface + assert {int(r) for r in table} <= set(res.event_study_effects) for r_str, cols in table.items(): r = int(r_str) - att, _se = cols[column] + att, _se = cols[column] # printed SE halves retired (see docstring) eff = res.event_study_effects[r] np.testing.assert_allclose(eff["effect"], att, atol=PRINTED_ATOL) - @XFAIL_EVENT_STUDY_GOLDENS @pytest.mark.parametrize( - "outcome,table_key", - [ - ("log_retail_emp", "table_a4_log_retail"), - ("log_wholesale_emp", "table_a5_log_wholesale"), - ], - ids=["a4_retail", "a5_wholesale"], + "outcome", + ["log_retail_emp", "log_wholesale_emp"], + ids=["retail", "wholesale"], ) - @pytest.mark.parametrize( - "rolling,estimator,column", - [ - ("detrend", "ra", "rolling_ra_detrend"), - ("detrend", "ipwra", "rolling_ipwra_detrend"), - ("demean", "ipwra", "rolling_ipwra_demean"), - ], - ) - def test_walmart_eventstudy_se_goldens( - self, walmart, golden, rolling, estimator, column, outcome, table_key - ): - """Bootstrap SEs vs the paper's printed B=999 draws (non-strict: - re-seeded bootstrap noise can sit near printed precision).""" - res = self._fit_es(walmart, rolling, estimator, outcome=outcome) - table = golden[table_key] - for r_str, cols in table.items(): - r = int(r_str) - _att, se = cols[column] - eff = res.event_study_effects[r] - np.testing.assert_allclose(eff["se"], se, atol=PRINTED_ATOL) + @pytest.mark.parametrize("rolling,method", [("detrend", "reg")]) + def test_walmart_eventstudy_se_vs_stata(self, stata_golden, es_fit, rolling, method, outcome): + """Bootstrap SEs vs the authors' Stata high-B multiplier bootstrap + (RA config ONLY - the verified-equivalent scheme). + + Both sides are FIXED-rep Monte Carlo estimates of the same bootstrap + SD, so the tolerance is the 4-sigma SE-of-SE bound + ``4*sqrt(1/(2*B_py) + 1/(2*B_stata))`` (4-sigma because the bound is + applied elementwise; at 3-sigma the family-wise chance of one MC + excursion across all rows is ~40%). Any failure = STOP and + investigate a scheme divergence; never widen this tolerance. + + The IPWRA configs are deliberately NOT gated here: measured + 2026-08-16, their bootstrap SEs diverge SYSTEMATICALLY from Stata's + (~15% at row level, far beyond the MC bound) - the known SE-form + difference (the contribution's IPWRA influence functions are + AIPW-style, not the papers' E.3 stacked form that Stata implements). + The Stata IPWRA SE columns stay in the golden as provenance; gating + them awaits the E.3 IF-form adjudication on the PR #588 work list. + """ + stata_method = {"reg": "ra", "dr": "ipwra"}[method] + key = f"{rolling}_{stata_method}__{outcome}" + table = stata_golden["walmart"][key]["watt"] + r_stata = stata_golden["meta"]["bootstrap_reps"] + rtol = 4 * np.sqrt(1 / (2 * 999) + 1 / (2 * r_stata)) + res = es_fit(rolling, method, outcome=outcome) + for r in sorted(res.event_study_effects): + cell = table.get(str(r)) + assert cell is not None and cell[1] is not None, f"{key} r={r} missing" + np.testing.assert_allclose( + res.event_study_effects[r]["se"], cell[1], rtol=rtol, err_msg=f"{key} r={r}" + ) - @XFAIL_EVENT_STUDY - def test_anchor_periods_excluded(self, walmart): - res_dm = self._fit_es(walmart, "demean", "ra") + def test_anchor_periods_excluded_from_raw_dict(self, es_fit): + res_dm = es_fit("demean", "reg") assert -1 not in res_dm.event_study_effects - res_dt = self._fit_es(walmart, "detrend", "ra") + res_dt = es_fit("detrend", "reg") assert -1 not in res_dt.event_study_effects assert -2 not in res_dt.event_study_effects + def test_public_surface_emits_observed_anchors_as_reference(self, es_fit): + res = es_fit("detrend", "reg") + es = res.aggregate("event_study") + ev = np.asarray(es.event_time) + ref = np.asarray(es.is_reference) + for anchor in (-2, -1): + idx = np.nonzero(ev == anchor)[0] + assert idx.size == 1, f"anchor {anchor} missing from public surface" + assert bool(ref[idx[0]]) + np.testing.assert_allclose(np.asarray(es.att)[idx[0]], 0.0) + def test_detrend_insample_residuals_sum_to_zero(self): """LW 2026 (pp20-21): per-unit detrended residuals sum to zero over the fitted pre-window - a property of OLS-with-intercept residuals @@ -776,9 +1040,8 @@ def test_detrend_insample_residuals_sum_to_zero(self): resid = pre["y"].to_numpy(dtype=float) - X @ beta np.testing.assert_allclose(resid.sum(), 0.0, atol=1e-9) - @XFAIL_EVENT_STUDY - def test_simultaneous_band_metadata(self, walmart): - res = self._fit_es(walmart, "detrend", "ra") + def test_simultaneous_band_metadata(self, es_fit): + res = es_fit("detrend", "reg") assert res.cband_method is not None assert res.cband_n_bootstrap >= 999 any_r = next(iter(res.event_study_effects)) @@ -786,7 +1049,545 @@ def test_simultaneous_band_metadata(self, walmart): # --------------------------------------------------------------------------- -# 9. Monte Carlo bias ordering (LW 2026, Section 5) +# 9. Full-precision Stata parity (authors' package, benchmarks/stata arm) +# --------------------------------------------------------------------------- + + +@pytest.mark.realdata +class TestStataParity: + """Full-precision anchors vs benchmarks/data/lwdid_stata_golden.json. + + The golden is generated by benchmarks/stata/generate_lwdid_golden.do + against the authors' SSC ``lwdid`` package (version pinned in its meta). + All comparisons use PARITY_RTOL/PARITY_ATOL (see module constants). + """ + + def test_prop99_att_se(self, prop99, stata_golden): + for rolling in ("demean", "detrend"): + res = _fit_prop99(prop99, rolling) + g = stata_golden["prop99"][rolling] + np.testing.assert_allclose(res.att, g["att"], rtol=PARITY_RTOL, atol=PARITY_ATOL) + np.testing.assert_allclose(res.se, g["se"], rtol=PARITY_RTOL, atol=PARITY_ATOL) + + def test_castle_att_se(self, castle, stata_golden): + """Stata small-N staggered mode is never-treated-based by + construction (no control option exists); the Python fits pass + control_group='never_treated' accordingly - the golden meta records + this equivalence.""" + for rolling in ("demean", "detrend"): + res = _fit_castle(castle, rolling) + g = stata_golden["castle"][rolling] + np.testing.assert_allclose(res.att, g["att"], rtol=PARITY_RTOL, atol=PARITY_ATOL) + np.testing.assert_allclose(res.se, g["se"], rtol=PARITY_RTOL, atol=PARITY_ATOL) + + def test_walmart_watt_points(self, stata_golden, es_fit): + """Per-r WATT point parity at full precision, over the measured raw + label set (any golden row absent from the raw dict must be + anchor-or-nonfinite; a finite Stata value missing from our surface is + a coverage hole, asserted in the label-set test).""" + for key in WALMART_CONFIG_KEYS: + rolling, method, outcome = _decode_walmart_key(key) + res = es_fit(rolling, method, outcome=outcome) + table = stata_golden["walmart"][key]["watt"] + rtol = PARITY_RTOL if method == "reg" else IPWRA_PARITY_RTOL + atol = PARITY_ATOL if method == "reg" else IPWRA_PARITY_ATOL + for r in sorted(res.event_study_effects): + cell = table.get(str(r)) + assert cell is not None and cell[0] is not None, f"{key} r={r} missing" + np.testing.assert_allclose( + res.event_study_effects[r]["effect"], + cell[0], + rtol=rtol, + atol=atol, + err_msg=f"{key} r={r}", + ) + + def test_walmart_post_avg_convention(self, stata_golden, es_fit): + """Stata's Post_avg equals the cell-mass (n_treated-weighted) average + of the finite post cells - the convention pin behind the REGISTRY + Aggregation note (verified to the last digit in the 2026-08-15 + smoke). Computed over finite-att cells with relative_time >= 0 (the + container includes non-estimable NaN cells).""" + for key in WALMART_CONFIG_KEYS: + rolling, method, outcome = _decode_walmart_key(key) + res = es_fit(rolling, method, outcome=outcome) + cells = [ + v + for v in res.cohort_time_effects.values() + if v["relative_time"] >= 0 and np.isfinite(v["att"]) + ] + att = np.array([v["att"] for v in cells]) + wt = np.array([v["n_treated"] for v in cells], dtype=float) + post_avg = float(np.average(att, weights=wt)) + golden = stata_golden["walmart"][key]["overall"]["Post_avg"][0] + rtol = PARITY_RTOL if method == "reg" else IPWRA_PARITY_RTOL + np.testing.assert_allclose(post_avg, golden, rtol=rtol, atol=PARITY_ATOL, err_msg=key) + + +# --------------------------------------------------------------------------- +# 10. Common-timing aggregation surface + degenerate/bootstrap contracts +# --------------------------------------------------------------------------- + + +class TestCommonTimingAggregationSurface: + """The guard relaxation covers event_study AND simple; group still raises.""" + + @pytest.fixture(scope="class") + def fitted(self): + df = _synthetic_common_timing(t_max=8, s=5) + return LWDiD(rolling="demean", estimation_method="reg", vcov_type="classical").fit( + df, outcome="y", unit="unit", time="time", treatment="treat" + ) + + def test_simple_relays_the_fit(self, fitted): + agg = fitted.aggregate("simple") + df_out = agg.to_dataframe() + np.testing.assert_allclose(df_out["att"].iloc[0], fitted.att, rtol=0, atol=0) + np.testing.assert_allclose(df_out["se"].iloc[0], fitted.se, rtol=0, atol=0) + + def test_group_raises_on_common_timing(self, fitted): + with pytest.raises(ValueError): + fitted.aggregate("group") + + def test_event_study_bootstrap_regime(self): + """Under n_bootstrap>0 the common-timing ES rows carry the bootstrap + regime exactly like the staggered surface: the full quintet is + finite and regime-consistent, and the band metadata is present.""" + df = _synthetic_common_timing(t_max=8, s=5) + res = LWDiD(rolling="demean", estimation_method="reg", n_bootstrap=199, seed=42).fit( + df, outcome="y", unit="unit", time="time", treatment="treat" + ) + es = res.aggregate("event_study") + ref = np.asarray(es.is_reference) + for arr in (es.att, es.se, es.t_stat, es.p_value, es.conf_int_lower, es.conf_int_upper): + vals = np.asarray(arr, dtype=float)[~ref] + assert np.isfinite(vals).all() + assert res.cband_method is not None + + def test_degenerate_fit_follows_house_empty_contract(self): + """A fit with no finite effects (single-pre-period detrend) yields + the house requested-but-empty event-study surface: zero + non-reference rows, any anchor row only if genuinely observed - and + never NaN non-reference rows or a ValueError.""" + df = TestMinimumPrePeriods._panel(first_period_treated=2) + with pytest.warns(UserWarning): + res = LWDiD(rolling="detrend", estimation_method="reg", vcov_type="classical").fit( + df, outcome="y", unit="unit", time="time", treatment="treat" + ) + es = res.aggregate("event_study") + ref = np.asarray(es.is_reference, dtype=bool) + assert (~ref).sum() == 0 + # single pre-period: the detrend anchor r=-2 is NOT observed and + # must not be synthesized + ev = set(int(e) for e in np.asarray(es.event_time)) + assert -2 not in ev + + +# --------------------------------------------------------------------------- +# 11. Shared input-validation contracts (every path) +# --------------------------------------------------------------------------- + + +class TestSharedValidationContracts: + """Unit-constancy and time-family contracts enforced in shared fit() + validation for every LWDiD path (common-timing AND staggered).""" + + @staticmethod + def _panel_and_kwargs(mode): + """(df, fit_kwargs, extra_init) for each design path.""" + if mode == "common": + df = _synthetic_common_timing(t_max=6, s=4) + kw = dict(outcome="y", unit="unit", time="time", treatment="treat") + extra = {} + else: + df = _synthetic_staggered() + kw = dict( + outcome="y", unit="unit", time="time", treatment="treat", first_treat="first_year" + ) + extra = dict(control_group="never_treated") + return df, kw, extra + + @pytest.mark.parametrize("mode", ["common", "staggered"]) + def test_time_varying_covariate_rejected(self, mode): + df, kw, extra = self._panel_and_kwargs(mode) + rng = np.random.default_rng(5) + df["xt"] = rng.normal(size=len(df)) # time-varying within unit + with pytest.raises(ValueError): + LWDiD(rolling="demean", estimation_method="reg", **extra).fit( + df, covariates=["xt"], **kw + ) + + @pytest.mark.parametrize("mode", ["common", "staggered"]) + def test_time_varying_cluster_rejected(self, mode): + df, kw, extra = self._panel_and_kwargs(mode) + rng = np.random.default_rng(6) + df["cl"] = rng.integers(0, 5, size=len(df)) # time-varying within unit + with pytest.raises(ValueError): + LWDiD(rolling="demean", estimation_method="reg", cluster="cl", **extra).fit(df, **kw) + + @pytest.mark.parametrize("mode", ["common", "staggered"]) + def test_row_order_invariance(self, mode): + """Shuffling row order within units leaves att/se identical (with + unit-CONSTANT covariates; the time-varying case is covered by the + rejection tests above, since post-fix it raises).""" + if mode == "common": + df = _synthetic_common_timing(t_max=6, s=4) + kw = dict(outcome="y", unit="unit", time="time", treatment="treat", covariates=["x"]) + est = dict(rolling="demean", estimation_method="reg", vcov_type="hc1") + else: + df = _synthetic_staggered() + # unit-constant covariate and cluster so both selection paths are + # exercised under reordering (time-varying cases raise, above) + units = df["unit"].unique() + rng = np.random.default_rng(41) + xmap = dict(zip(units, rng.normal(size=len(units)))) + df["xu"] = df["unit"].map(xmap) + kw = dict( + outcome="y", + unit="unit", + time="time", + treatment="treat", + first_treat="first_year", + covariates=["xu"], + ) + # cluster selection under reordering is covered by the + # clustered-inference oracle + the time-varying rejection above; + # keeping cluster= out here keeps this invariance lane executable + # on the pre-rename head too. + est = dict( + rolling="demean", + estimation_method="reg", + vcov_type="hc1", + control_group="never_treated", + ) + shuffled = df.sample(frac=1.0, random_state=99).reset_index(drop=True) + r1 = LWDiD(**est).fit(df, **kw) + r2 = LWDiD(**est).fit(shuffled, **kw) + np.testing.assert_allclose(r1.att, r2.att, rtol=0, atol=1e-12) + np.testing.assert_allclose(r1.se, r2.se, rtol=0, atol=1e-12) + + @pytest.mark.parametrize( + "direction", ["datetime_time_period_cohort", "period_time_datetime_cohort"] + ) + def test_mixed_time_families_rejected(self, direction): + """Mixed datetime/Period time families raise the estimator's + documented ValueError, not an internal pandas TypeError (both + directions pass the boolean date-like check today).""" + periods = pd.date_range("2020-01-01", periods=6, freq="YS") + rows = [] + for unit, cohort_i in [("t1", 2), ("t2", 2), ("c1", None), ("c2", None)]: + for i, p in enumerate(periods): + treat = int(cohort_i is not None and i >= cohort_i) + rows.append({"unit": unit, "i": i, "time": p, "treat": treat, "y": float(i)}) + df = pd.DataFrame(rows) + if direction == "datetime_time_period_cohort": + df["cohort"] = df["unit"].map( + {"t1": pd.Period("2022", freq="Y"), "t2": pd.Period("2022", freq="Y")} + ) + else: + df["time"] = df["time"].dt.to_period("Y") + df["cohort"] = df["unit"].map( + {"t1": pd.Timestamp("2022-01-01"), "t2": pd.Timestamp("2022-01-01")} + ) + with pytest.raises(ValueError): + LWDiD(rolling="demean", estimation_method="reg", control_group="never_treated").fit( + df, outcome="y", unit="unit", time="time", treatment="treat", first_treat="cohort" + ) + + def test_dr_without_covariates_warns_and_equals_reg(self): + """no-silent-failures: a dr fit with no covariates degrades to reg + and must SAY so (the head's IPWRA branch is silent against its own + routing docstring; ipw and psm already warn). Warning contract per + CONTRIBUTING: message assertion + behavioral assertion.""" + df = _synthetic_common_timing(t_max=6, s=4) + kw = dict(outcome="y", unit="unit", time="time", treatment="treat") + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + r_dr = LWDiD(rolling="demean", estimation_method="dr").fit(df, **kw) + assert any( + any(s in str(w.message).lower() for s in ("covariate", "control", "fall")) + for w in caught + ), "dr-without-covariates fallback must warn" + r_reg = LWDiD(rolling="demean", estimation_method="reg").fit(df, **kw) + np.testing.assert_allclose(r_dr.att, r_reg.att, rtol=0, atol=1e-12) + np.testing.assert_allclose(r_dr.se, r_reg.se, rtol=0, atol=1e-12) + + def test_ipw_without_covariates_warns(self): + """The ipw branch already warns at the head (expected-pass lane).""" + df = _synthetic_common_timing(t_max=6, s=4) + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + LWDiD(rolling="demean", estimation_method="ipw").fit( + df, outcome="y", unit="unit", time="time", treatment="treat" + ) + assert any( + any(s in str(w.message).lower() for s in ("covariate", "control", "fall")) + for w in caught + ) + + +# --------------------------------------------------------------------------- +# 12. Sensitivity module contracts +# --------------------------------------------------------------------------- + + +class TestSensitivityContracts: + """lwdid_sensitivity failure contracts + rename sweep.""" + + @staticmethod + def _df(): + return _synthetic_common_timing(n_treat=20, n_control=30, t_max=8, s=5, effect=0.0) + + def test_zero_baseline_is_not_estimable(self): + """A near-zero baseline ATT must classify as not_estimable with NaN + ratio, never highly_robust (the relative ratio is undefined).""" + from diff_diff.lwdid_sensitivity import robustness_pre_periods + + df = self._df() + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + r = robustness_pre_periods(df, outcome="y", unit="unit", time="time", treatment="treat") + if abs(getattr(r, "baseline_att", 1.0)) < 1e-10: + assert r.robustness_level == "not_estimable" + assert np.isnan(r.sensitivity_ratio) + else: + # zero-effect DGP but baseline not exactly ~0: exercise the + # classifier directly on a forced zero baseline + from diff_diff.lwdid_sensitivity import _compute_sensitivity_ratio + + ratio = _compute_sensitivity_ratio(0.0, [0.0, 1.0, -1.0]) + assert np.isnan(ratio) + + def test_missing_column_propagates(self): + """A missing-column programming error must PROPAGATE as ValueError, + not be converted into a silently failed specification (narrowed + except contract).""" + from diff_diff.lwdid_sensitivity import robustness_pre_periods + + with pytest.raises(ValueError): + robustness_pre_periods( + self._df(), outcome="missing_outcome", unit="unit", time="time", treatment="treat" + ) + + def test_public_entry_points_use_canonical_kwargs(self): + """The sensitivity API carries the post-rename vocabulary; the old + kwargs are rejected.""" + import inspect + + from diff_diff.lwdid_sensitivity import robustness_pre_periods, sensitivity_no_anticipation + + for fn in (robustness_pre_periods, sensitivity_no_anticipation): + params = inspect.signature(fn).parameters + assert "estimation_method" in params, fn.__name__ + assert "estimator" not in params, fn.__name__ + assert "vce" not in params, fn.__name__ + + def test_specification_result_serializes_canonical_key(self): + from diff_diff.lwdid_sensitivity import robustness_pre_periods + + df = _synthetic_common_timing(n_treat=20, n_control=30, t_max=8, s=5) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + r = robustness_pre_periods(df, outcome="y", unit="unit", time="time", treatment="treat") + payload = json.dumps(r.to_dict()) if hasattr(r, "to_dict") else "{}" + assert '"estimator"' not in payload + assert '"vce"' not in payload + + +# --------------------------------------------------------------------------- +# 13. Retired pre-release surfaces (negative assertions) +# --------------------------------------------------------------------------- + + +class TestRetiredSurfaces: + """The locked pre-release renames/removals, enforced. + + Versioning rationale: NO deprecation shims or ledger rows - LWDiD has + never shipped on main, so these are pre-release renames, not breaking + changes (see the REGISTRY validation-suite note). + """ + + @staticmethod + def _small_fit(**kwargs): + df = _synthetic_common_timing(n_treat=10, n_control=15, t_max=6, s=4) + est = LWDiD(rolling="demean", estimation_method="reg", **kwargs) + return est.fit(df, outcome="y", unit="unit", time="time", treatment="treat") + + def test_lw_alias_gone(self): + import diff_diff + + assert not hasattr(diff_diff, "LW") + assert "LW" not in diff_diff.__all__ + + def test_canonical_export(self): + import diff_diff + + assert "LWDiD" in diff_diff.__all__ + + def test_module_level_wrapper_gone(self): + mod = importlib.import_module("diff_diff.lwdid") + assert not hasattr(mod, "lwdid") + + def test_trend_diagnostics_module_deleted(self): + with pytest.raises(ModuleNotFoundError): + importlib.import_module("diff_diff.lwdid_trend_diagnostics") + + def test_retired_result_fields_gone(self): + res = self._small_fit() + assert not hasattr(res, "overall_att") + assert not hasattr(res, "period_effects") + payload = res.to_dict() + assert "overall_att" not in payload + assert "period_effects" not in payload + + def test_result_field_mirrors_renamed(self): + res = self._small_fit() + assert not hasattr(res, "estimator") or res.estimator == "LWDiD" # house NAME field only + assert not hasattr(res, "vce_type") + assert res.estimation_method == "reg" + assert res.vcov_type == "hc1" or res.vcov_type == "classical" + d = res.to_dict() + assert "estimation_method" in d and "vcov_type" in d + assert "vce_type" not in d + frame_cols_or_payload = json.dumps(d) + assert '"vce_type"' not in frame_cols_or_payload + text = res.summary() or "" + rep = repr(res) + for banned in ("Estimator:", "estimator=", "vce="): + assert banned not in text + assert banned not in rep + assert "estimation_method" in rep or "reg" in text + + def test_psm_routes(self): + """psm actually routes (with covariates - the covariate-free PSM + branch falls back to reg) and differs from reg on the same data.""" + df = _synthetic_common_timing(t_max=6, s=4) + kw = dict(outcome="y", unit="unit", time="time", treatment="treat", covariates=["x"]) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + r_psm = LWDiD(rolling="demean", estimation_method="psm").fit(df, **kw) + r_reg = LWDiD(rolling="demean", estimation_method="reg").fit(df, **kw) + assert np.isfinite(r_psm.att) + assert r_psm.estimation_method == "psm" + assert abs(r_psm.att - r_reg.att) > 1e-12 + + def test_legacy_values_rejected(self): + for legacy in ("ra", "ipwra"): + with pytest.raises(ValueError): + LWDiD(estimation_method=legacy) + + def test_legacy_kwargs_rejected(self): + for kwargs in ( + {"estimator": "ra"}, + {"vce": "hc1"}, + {"bootstrap_seed": 42}, + {"period_specific": True}, + {"trim_threshold": 0.05}, + ): + with pytest.raises(TypeError): + LWDiD(**kwargs) + + def test_cluster_placement(self): + """cluster= lives on the constructor (CS/LPDiD placement); fit() + rejects it.""" + LWDiD(cluster="cl") # accepted + df = _synthetic_common_timing(t_max=6, s=4) + with pytest.raises(TypeError): + LWDiD(rolling="demean", estimation_method="reg").fit( + df, outcome="y", unit="unit", time="time", treatment="treat", cluster="unit" + ) + + def test_post_rename_defaults(self): + params = LWDiD().get_params() + assert params["seed"] is None + assert params["estimation_method"] == "reg" + assert params["vcov_type"] == "hc1" + assert params["pscore_trim"] == 0.01 + + def test_vcov_cluster_value_retired(self): + """The 'cluster' MODE value is retired AT THE CONSTRUCTOR (clustering + activates via the cluster= param); a fit-time probe would pass at the + pre-rename head for the wrong reason.""" + with pytest.raises(ValueError, match="cluster"): + LWDiD(vcov_type="cluster") + + def test_cluster_param_produces_clustered_inference(self): + """Behavioral oracle: with induced within-cluster correlation the + clustered SE differs from the unclustered hc1 SE (metadata alone + would pass a silent-HC1 implementation), and the house + cluster-metadata contract is honored.""" + rng = np.random.default_rng(17) + # cluster-level POST-period shocks: a unit-level or time-invariant + # shock is removed by the pre-period demeaning, leaving nothing for + # clustering to pick up - the collapsed outcomes must be genuinely + # correlated WITHIN clusters for this oracle to bite. + cluster_post_shock = {c: rng.normal(0, 1.5) for c in range(10)} + rows = [] + for i in range(80): + cl = i // 8 # 10 clusters + is_treated = i < 30 + for t in range(1, 7): + post = t >= 4 + y = ( + (cluster_post_shock[cl] if post else 0.0) + + rng.normal(0, 0.4) + + (1.0 if is_treated and post else 0.0) + ) + rows.append( + {"unit": i, "cl": cl, "time": t, "y": y, "treat": int(is_treated and post)} + ) + df = pd.DataFrame(rows) + r_cl = LWDiD(rolling="demean", estimation_method="reg", cluster="cl").fit( + df, outcome="y", unit="unit", time="time", treatment="treat" + ) + r_un = LWDiD(rolling="demean", estimation_method="reg", vcov_type="hc1").fit( + df, outcome="y", unit="unit", time="time", treatment="treat" + ) + assert r_cl.cluster_name == "cl" + assert r_cl.n_clusters == 10 + # strong within-cluster correlation must MATERIALLY widen the SE, + # not merely change it (a silent-HC1 implementation with metadata + # would fail here; measured ratio ~2.7 on this fixture) + assert r_cl.se > 1.15 * r_un.se, (r_cl.se, r_un.se) + # ... and the value itself must equal the CR1 sandwich on the + # collapsed cross-section, small-sample factor G/(G-1)*(N-1)/(N-K) + # (verified identical to the contribution head's clustered path) + ybar = _demean_reference(df, "unit", "time", "y", pre_end=3) + d = df.groupby("unit")["treat"].max().loc[ybar.index].to_numpy(dtype=float) + cl_ids = df.groupby("unit")["cl"].first().loc[ybar.index].to_numpy() + y = ybar.to_numpy(dtype=float) + X = np.column_stack([np.ones_like(d), d]) + beta, *_ = np.linalg.lstsq(X, y, rcond=None) + e = y - X @ beta + n_units, k = X.shape + n_cl = len(np.unique(cl_ids)) + bread = np.linalg.inv(X.T @ X) + meat = np.zeros((k, k)) + for g_id in np.unique(cl_ids): + m = cl_ids == g_id + s = X[m].T @ e[m] + meat += np.outer(s, s) + factor = (n_cl / (n_cl - 1)) * ((n_units - 1) / (n_units - k)) + se_cr1 = float(np.sqrt((factor * bread @ meat @ bread)[1, 1])) + np.testing.assert_allclose(r_cl.se, se_cr1, rtol=1e-10) + + def test_pscore_trim_affects_behavior(self): + """set_params(pscore_trim=...) actually changes an ipw fit (the + generic roster's mutation lane picks alpha, so nothing else + exercises this parameter).""" + df = _synthetic_common_timing(t_max=6, s=4) + kw = dict(outcome="y", unit="unit", time="time", treatment="treat", covariates=["x"]) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + r1 = LWDiD(rolling="demean", estimation_method="ipw").fit(df, **kw) + est = LWDiD(rolling="demean", estimation_method="ipw") + est.set_params(pscore_trim=0.4) + r2 = est.fit(df, **kw) + assert abs(r1.att - r2.att) > 0 or abs(r1.se - r2.se) > 0 + + +# --------------------------------------------------------------------------- +# 14. Monte Carlo bias ordering (LW 2026, Section 5) # --------------------------------------------------------------------------- @@ -832,7 +1633,7 @@ def _one_rep(self, rng): ) out = {} for rolling in ("demean", "detrend"): - res = LWDiD(rolling=rolling, estimator="ra", vce="classical").fit( + res = LWDiD(rolling=rolling, estimation_method="reg", vcov_type="classical").fit( df, outcome="y", unit="unit", time="time", treatment="treat" ) out[rolling] = res.att - sample_att @@ -858,7 +1659,7 @@ def test_bias_ordering_under_heterogeneous_trends(self): # --------------------------------------------------------------------------- -# 10. Minimum pre-treatment periods (REGISTRY edge case) +# 15. Minimum pre-treatment periods (REGISTRY edge case) # --------------------------------------------------------------------------- @@ -886,7 +1687,7 @@ def _panel(first_period_treated): def test_demeaning_with_single_pre_period_works(self): df = self._panel(first_period_treated=2) # exactly one pre-period - res = LWDiD(rolling="demean", estimator="ra", vce="classical").fit( + res = LWDiD(rolling="demean", estimation_method="reg", vcov_type="classical").fit( df, outcome="y", unit="unit", time="time", treatment="treat" ) assert np.isfinite(res.att) @@ -897,7 +1698,7 @@ def test_detrending_with_single_pre_period_warns_and_nans(self): NaN inference - loud, house-compatible (warn + NaN, never silent).""" df = self._panel(first_period_treated=2) # one pre-period: rank-deficient with pytest.warns(UserWarning, match="at least 2 pre-treatment periods"): - res = LWDiD(rolling="detrend", estimator="ra", vce="classical").fit( + res = LWDiD(rolling="detrend", estimation_method="reg", vcov_type="classical").fit( df, outcome="y", unit="unit", time="time", treatment="treat" ) assert np.isnan(res.att) @@ -909,7 +1710,7 @@ def test_detrending_with_single_pre_period_warns_and_nans(self): def test_detrending_with_two_pre_periods_works(self): df = self._panel(first_period_treated=3) # two pre-periods: minimum - res = LWDiD(rolling="detrend", estimator="ra", vce="classical").fit( + res = LWDiD(rolling="detrend", estimation_method="reg", vcov_type="classical").fit( df, outcome="y", unit="unit", time="time", treatment="treat" ) assert np.isfinite(res.att) diff --git a/tests/test_naming_guard.py b/tests/test_naming_guard.py index 2b50716c..bc41306b 100644 --- a/tests/test_naming_guard.py +++ b/tests/test_naming_guard.py @@ -1020,6 +1020,13 @@ def _token_family_code_refs(tok): "canonical successor method; HAD-section fit-time mode mentions " "were migrated with M-027" ), + ("period_effects", "docs/methodology/REGISTRY.md"): ( + "the LWDiD API-conformance/aggregation Notes name PR #588's " + "pre-release `period_effects` surface being RETIRED before it ever " + "ships - not a read of MultiPeriodDiDResults.period_effects (M-016), " + "whose reader inventory this registry entry documents rather than " + "consumes" + ), ("aggregate", "docs/methodology/REPORTING.md"): ( "canonical post-fit aggregate() mention only (zero fit-time " "kwarg sites; the fit-time-population clause was reworded with "