diff --git a/CHANGELOG.md b/CHANGELOG.md index 121fafb0..c7ba25a5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,18 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +### Added +- **Chang (2020) §4.2.2 RCS simulation-DGP replication** (`tests/test_methodology_dml_did.py`, + DML PR-B2): the paper's own kernel-design repeated-cross-section DGP as + maintainer validation fixtures for `DMLDiD(panel=False)` — a DGP-shape pin + (distributions, all three innovation scales, the design's built-in confounded + contrast → θ₀+1, and both correct-specification facts), seed-pinned recovery at + both paper sample sizes with a discriminating comparison against the unadjusted + contrast, and a slow Monte Carlo coverage lane. The §4.2 parameterizations are + extracted into the paper review; the §4.2.1 ML design is documented as not + replicable with the bundled unpenalized learners (narrowed TODO row) and the + REGISTRY carries the replication-scope Note. + ### Changed - **`n_bootstrap` type guards aligned onto `utils.validate_n_bootstrap`** for the estimators the M-081 sweep deliberately left out — diff --git a/TODO.md b/TODO.md index 00502909..b8152d2d 100644 --- a/TODO.md +++ b/TODO.md @@ -78,7 +78,7 @@ generic sparse-FE, QR+SVD rank-detection redundancy, `check_finite` bypass — m | Issue | Location | Origin | Effort | Priority | |-------|----------|--------|--------|----------| -| Replicate Chang (2020) §4's own RCS simulation DGPs (pp. 17-21, "fully specified" per the paper review) as recovery/coverage fixtures for the `DMLDiD(panel=False)` lane — the shipped tests use a library-authored RCS design (documented in the REGISTRY checklist caveat); needs the paper PDF to extract the parameterization | `tests/test_methodology_dml_did.py`, `docs/methodology/papers/chang-2020-review.md` | DML PR-B2 | Mid | Low | +| Replicate Chang (2020) §4.2.1 ML-design RCS DGP — requires a penalized (Lasso-class) propensity learner or a maintainer fixture via the duck-typed `propensity_learner=` object route; native logit+linear verified pure noise at (N=500, p=100) across ~20 review seeds (att ~ −23..+16, SE ~4-10, EPV ~1.4-1.6 vs threshold 10, fitted out-of-fold clipping ~5-29%, mean ~14%); at (N=200, p=100) `outcome_learner="linear"` fails closed on control-fold rank deficiency even under an oracle propensity (use ridge/sieve there); §4.2.2 is replicated | `tests/test_methodology_dml_did.py`, `diff_diff/_learners.py` | DML PR-B2 | Mid | Low | | Optional scheduled end-to-end execution gate for the MMM tutorials (29/30): a cron-only workflow (or extension of `mmm-interop.yml`) that executes both notebooks in isolated exact-pin environments, so a stale/invalid committed posterior cannot stay green indefinitely - today the hybrid posture (deliberate: notebooks execute locally with committed outputs; CI smoke-tests the exporters without sampling; drift tests pin source + committed-output needles) leaves the MCMC claims un-re-executed in CI | `.github/workflows/mmm-interop.yml`, `docs/tutorials/29_mmm_calibration_pymc.ipynb`, `docs/tutorials/30_mmm_calibration_meridian.ipynb` | mmm-interop | Mid | Low | | Committed `fixest::feols` event-study golden for TWFE `event_study=True` (within + pooled specs, unbalanced + covariate panels, matched CR1 cluster convention, per-period effects + vcov block) - the in-suite gates are shared-core cross-checks (TWFE-within == MPD-absorb, pooled == MPD bit-exact), so a defect common to the shared core would pass; the live-R harness (`benchmarks/R/benchmark_multiperiod.R`, `feols(y ~ treated * time_f \| unit)`) validated the within design in `docs/benchmarks.rst` but is not a committed regression test - follow the `fixest_did_twfe_golden.json` committed-golden pattern (pytest.skip when absent) | `tests/test_fixest_did_twfe_parity.py`, `benchmarks/R/` | 3(a) R2 | Mid | Medium | | Evaluate adding the `BaseEstimator` param surface (get_params/set_params) to the exported classes that never had it - `PowerAnalysis`, `LinearRegression`, `BusinessReport`, `DiagnosticReport`, `TWFEWeightsResult` (a NEW public surface, deliberately out of the 2(c)-i pure-refactor scope; `LinearRegression` is the one `fit`-bearing class excluded from the contract suite's roster-completeness test). | `diff_diff/linalg.py`, `diff_diff/power.py` | mixin PR | Mid | Low | diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index 629d6b99..654bd878 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -3129,6 +3129,8 @@ the finite-dimensional `p_0` is handled by the variance correction below. at 0.0 observed diff under identical folds and sklearn learners). - R / Stata: none — the paper ships no companion package. +- **Note:** Chang §4.2 simulation-DGP replication scope (PR-B2). The §4.2.2 kernel-design RCS DGP (arXiv v3 p. 19) is replicated verbatim in `tests/test_methodology_dml_did.py`, estimated with the library's NATIVE learners rather than the paper's Gaussian-kernel first stages (no kernel first stage is bundled); the fixture is correctly specified on both nuisances (true propensity `sigmoid(X − 1/2)` by Bayes' rule; true `ℓ₂₀ = 0.25·X`), so it validates paper-DGP faithfulness — the discriminating assertion compares against the design's own confounded unadjusted contrast (→ θ₀ + 1), not against a misspecified learner. The design's unbounded X means the a.s. strict-overlap Assumption 3.2(a) holds for no fixed κ > 0 (mass outside [0.01, 0.99] ≈ 2.1e-5 — immaterial at the fixture Ns; the same quirk applies more strongly to §4.2.1, tail ≈ 0.06%). The §4.2.1 ML design is NOT replicable with the bundled unpenalized learners: at (N=500, p=100) the native logit+linear stack is order-of-magnitude noise across seeds (SE ~4-10; `solve_logit` EPV ~1.4-1.6 vs threshold 10; fitted out-of-fold clipping ~5-29% of rows — estimated propensities, not true ones, whose out-of-[0.01, 0.99] mass is ≈0.06%), and at (N=200, p=100) the linear outcome learner fails closed on control-fold rank deficiency even under an oracle propensity (ridge/sieve fit there) — tracked in TODO.md (penalized propensity learner + p-robust outcome learner). + **Requirements checklist (shipped Case 1 panel + Case 2 RCS staggered lanes):** - [x] Neyman-orthogonal Case 1 score (3.1) implemented exactly (Abadie score + mean-zero adjustment; `chang_panel_score`) - [x] DML2 cross-fitting: per-cell K-fold partition, nuisances fit on fold complements, never on the evaluation fold @@ -3139,7 +3141,7 @@ the finite-dimensional `p_0` is handled by the variance correction below. - [x] Per-cell degenerate guards (zero treated/control, cell < K, singleton stratum, empty untreated complement) — closed skip vocabulary, consolidated warning - [x] Normal-approximation inference via `safe_inference()` - [x] Validation: 2-period DoubleMLDID + staggered per-cell DoubleMLDIDBinary parity spikes (version-pinned, committed, golden literals consumed in-tests) + oracle-nuisance closed-form equivalence + degenerate-cell hand-pipeline equivalence (rtol 1e-14 — BLAS reduces differently-laid-out inputs in platform-dependent order, so bit identity does not hold cross-platform) + Monte Carlo coverage sanity -- [x] Case 2 (repeated cross sections): Equation 3.2 score + λ-corrected Theorem 2 variance — SHIPPED as `DMLDiD(panel=False)` (`chang_rcs_score` / `chang_rcs_lambda_slope` / `chang_rcs_score_augmented`; equation-level fixtures, oracle closed forms, derivative-identity checks, DR both directions, characterization spike, MC coverage). The paper's own §4 RCS simulation DGPs are NOT replicated (tracked TODO row — needs the paper PDF pp. 17-21); the shipped recovery tests use a library-authored RCS design. +- [x] Case 2 (repeated cross sections): Equation 3.2 score + λ-corrected Theorem 2 variance — SHIPPED as `DMLDiD(panel=False)` (`chang_rcs_score` / `chang_rcs_lambda_slope` / `chang_rcs_score_augmented`; equation-level fixtures, oracle closed forms, derivative-identity checks, DR both directions, characterization spike, MC coverage). The paper's own §4.2.2 kernel-design RCS DGP is replicated (`tests/test_methodology_dml_did.py` "Chang Sec. 4.2.2" section: DGP-shape pin, seed-pinned recovery at both paper Ns with a discriminating comparison against the design's confounded unadjusted contrast, slow MC coverage — see the §4.2 replication Note below); the §4.2.1 ML design is not replicable with the bundled learners (narrowed TODO.md row). - [ ] Case 3 (multilevel treatment): deferred (`DEFERRED.md`; implementation-required overlap conditions per the paper review's Case 3 caution) --- diff --git a/docs/methodology/papers/chang-2020-review.md b/docs/methodology/papers/chang-2020-review.md index a554dbc4..49d9c06c 100644 --- a/docs/methodology/papers/chang-2020-review.md +++ b/docs/methodology/papers/chang-2020-review.md @@ -192,7 +192,7 @@ items remain open, tracked in DEFERRED.md): - [x] Normal-approximation inference via `safe_inference()` (no-design fits; survey/bare-cluster fits use finite-df t inference via `df=df_survey` — library extension, not from the paper) - [ ] Multilevel treatment (Case 3): open — DEFERRED - [ ] Case 3 guards: open — DEFERRED (see the Case 3 caution above) -- [x] Validation (Case 1): `doubleml.DoubleMLDID` (2-period) + `DoubleMLDIDBinary` (staggered per-cell, end-to-end public fit) parity spikes, doubleml==0.11.4 pinned, golden literals in-repo; [x] Case 2 equation-level fixtures — SHIPPED (closed-form/oracle fixtures, derivative-identity checks, DR both directions, `DoubleMLDIDCSBinary` characterization spike — no parity oracle exists). CAVEAT: the paper's own §4 RCS simulation DGPs (pp. 17-21) are NOT replicated — the shipped recovery/coverage tests use a library-authored RCS design; replication is a tracked TODO.md row (needs the paper PDF); [ ] Case 3 fixtures — open +- [x] Validation (Case 1): `doubleml.DoubleMLDID` (2-period) + `DoubleMLDIDBinary` (staggered per-cell, end-to-end public fit) parity spikes, doubleml==0.11.4 pinned, golden literals in-repo; [x] Case 2 equation-level fixtures — SHIPPED (closed-form/oracle fixtures, derivative-identity checks, DR both directions, `DoubleMLDIDCSBinary` characterization spike — no parity oracle exists). [x] §4.2.2 kernel-design RCS DGP replicated (`tests/test_methodology_dml_did.py`, "Chang Sec. 4.2.2" section: DGP-shape pin, seed-pinned recovery at both paper Ns with a discriminating comparison against the design's own confounded unadjusted contrast, slow MC coverage; the library's native learners stand in for the paper's Gaussian-kernel first stages, and the fixture is correctly specified on BOTH nuisances, so it adds paper-DGP faithfulness, not plug-in-bias reproduction). The §4.2.1 ML-design DGP is extracted below but NOT replicable with the bundled unpenalized learners — narrowed TODO.md row (penalized propensity learner + ridge/sieve outcome learner required); [ ] Case 3 fixtures — open --- @@ -244,6 +244,8 @@ items remain open, tracked in DEFERRED.md): 7. **Case 2 response in the Lasso penalty-loading recipe (appendix "The Lasso Penalty", p. 26):** the recipe literally prints "Let `y_i` denote `Y_i(1) - Y_i(0)` or `(T_i - λ̂_k)`" — the Case 2 response omits the `Y_i` multiplier (verified against the PDF page 2026-08-22), while the modified-Lasso objective it feeds (p. 12) minimizes over `(T_i - λ̂_k)·Y_i - q_i'β`. Loadings computed on a bare `(T_i - λ̂_k)` would produce a different first-stage estimator; implement loadings on the objective's response, `(T_i - λ̂_k)·Y_i`, pending comparison with the published article. 8. **`ȳ_k` population/denominator inconsistency (p. 26):** `ȳ_k = M^{-1} Σ_{i∈I_k^c} y_i` is defined as the mean over the FULL auxiliary sample (size `M`), while both loading sums run over the untreated subset `I_kz^c` with normalizer `1/M_k` (verified against the PDF page 2026-08-22). Whether the centering mean was intended over `I_k^c` or `I_kz^c` is not resolvable from the arXiv text; record the choice made at implementation time and compare with the published article. +9. **§4.2.1's printed-but-unused `β_0` (p. 19):** the RCS ML simulation section prints "where `β_0 = γ_0 + 0.5`" although `β_0` appears in none of its equations (`Y⁰(0) = 1 + ε_1` as printed), while §4.1.1 (p. 18) and §4.3.1 (p. 20) print `Y⁰(0) = X'β_0 + ε_1` in the parallel position — a suspected dropped `X'β_0` term; not resolvable from the arXiv text. Either reading leaves §4.2.1's parallel trends intact (the level term cancels in the recursive trend `Y⁰(1) = Y⁰(0) + 1 + ε_2`) and `ℓ_20` constant in X, so the library's §4.2.1 non-replicability conclusion (above) is robust to the ambiguity; check the published version. + **Contradiction — where `p̂_k` (and `λ̂_k`) is computed (algorithms vs proofs):** - Both printed algorithm statements agree with each other and were verified against the PDF pages directly (2026-08-22): main-text Algorithm 1 (p. 11) prints `p̂_k = (1/n) Σ_{i∈I_k^c} D_i` and the appendix multilevel algorithm (p. 25, step 2) prints `p̂_w = (1/n) Σ_{i∈I_k^c} D_i` — in both cases a `1/n` normalizer (n = |I_k|) over the **auxiliary sample `I_k^c`** (size `M = N - n`), which is not a valid mean of anything as printed (typo 1 above). - The proofs of Theorems 1-2, however, use `p̂_k = E_{n,k}[D]` (fold mean over the **main fold `I_k`**, pp. 31, 33, 37), i.e. `p̂_k - p_0 = E_{n,k}[D - p_0]`, which is what generates the `G_1p0 (D - p_0)` influence-correction term in the variance algebra. @@ -258,6 +260,8 @@ items remain open, tracked in DEFERRED.md): - Multilevel treatment: variance estimator and full theorem statements are not printed in the reviewed text (Theorems 2/4 cover Cases 1-2 only; multilevel results "can be proven using the same arguments", p. 13). The analogy `Ĝ_wp = -θ̃_w/p̂_w` is plausible but NOT stated in the paper. - Theorem 2 proof for the repeated-cross-sections case concludes within the reviewed range (pp. 51-56); no gaps remain in the proofs of Theorems 1-4 as reviewed, but the appendix's stated Lasso-penalty constants `c`, `γ`, `B` have no defaults (p. 26). - No clustering, no bootstrap, no missing-data handling, no aggregation across multilevel treatment arms. +- The §4 simulation designs violate the paper's own regularity conditions in two symmetric ways: the Gaussian kernel used for the kernel-estimation variants violates Assumption 3.3's compact-support requirement (see the Tuning Parameters kernel row), and BOTH §4.2 propensities are unbounded-Gaussian-index logistics, so the a.s. strict-overlap Assumption 3.2(a) (`Pr(κ ≤ g_0(X) ≤ 1−κ) = 1` for fixed `κ > 0`) holds for NO fixed κ — §4.2.1's tail mass outside [0.01, 0.99] is ≈0.06% (`X'γ_0 ~ N(0.685, 1.21²)`), §4.2.2's ≈2.1e-5 (immaterial at the simulated Ns, but a theory/design inconsistency worth recording). - Empirical application (Section 5, pp. 21-22, repeated cross sections via Eqs. 2.2/3.2): Sequeira (2016) tariff/bribery data, N = 1084. Table 1 exact estimate (SE) pairs (verified against the PDF 2026-08-22): Sequeira (2016) TWFE -2.928 (0.944); Abadie kernel -7.986 (3.028); orthogonal θ̃ kernel -8.670 (3.643); Abadie Lasso -7.499 (2.746); orthogonal θ̃ Lasso -9.191 (4.854) — usable as a rough replication target only if the Sequeira data is obtainable. -- Simulation DGPs (Section 4, pp. 17-21) are fully specified and are the recommended validation fixtures: e.g. 4.1.1 repeated outcomes ML: `N ∈ {200, 500}`, `p ∈ {100, 300}`, `X ~ N(0, I_p)`, `γ_0 = (1, 1/2, 1/3, 1/4, 1/5, 0, ...)`, logistic PS, `β_0 = γ_0 + 0.5`, `θ_0 = 3`, errors N(0, 0.1) (variance/SD as printed "N(0,0.1)" — ambiguous); 4.3.1 multilevel: `W ∈ {0,1,2}`, shares (0.3, 0.3, 0.4), `θ_10 = 3`, `θ_20 = 6`. +- Simulation DGPs (Section 4, pp. 17-21) are fully specified and are the recommended validation fixtures: e.g. 4.1.1 repeated outcomes ML: `N ∈ {200, 500}`, `p ∈ {100, 300}`, `X ~ N(0, I_p)`, `γ_0 = (1, 1/2, 1/3, 1/4, 1/5, 0, ...)`, logistic PS, `β_0 = γ_0 + 0.5`, `θ_0 = 3`, errors N(0, 0.1) (variance/SD as printed "N(0,0.1)" — ambiguous as printed; the shipped §4.2.2 fixtures ADOPT the variance-0.1 reading, `σ = √0.1`, stated in the fixture comment); 4.3.1 multilevel: `W ∈ {0,1,2}`, shares (0.3, 0.3, 0.4), `θ_10 = 3`, `θ_20 = 6`. +- §4.2 repeated-cross-section DGPs (pp. 18-20, extracted 2026-08-29 for the PR-B2 fixtures). **§4.2.1 (ML estimation):** `N ∈ {200, 500}`, `p ∈ {100, 300}`, `X_i ~ N(0.3·1, I_p)` (note the 0.3 mean, unlike §4.1.1), same `γ_0`, logistic PS `P(D=1|X) = sigmoid(X'γ_0)`; as printed `Y⁰(0) = 1 + ε_1`, `Y⁰(1) = Y⁰(0) + 1 + ε_2`, `Y¹(1) = θ_0 + Y⁰(1) + ε_3`, `θ_0 = 3`, `T ~ Bernoulli(0.5)`, observed `Y = Y(0) + T(Y(1) − Y(0))` (see suspected-typo 9 on the printed-but-unused `β_0`). Under our native unpenalized learners this design is NOT estimable: at (N=500, p=100) the logit+linear stack is pure noise across ~20 review seeds (att ~ −23..+16, SE ~4-10; `solve_logit`'s EPV guard fires every fold at EPV ~1.4-1.6 vs threshold 10; FITTED out-of-fold clipping ~5-29% of rows across seeds, mean ~14% — the estimated-propensity leg of strict overlap; TRUE-propensity mass outside [0.01, 0.99] is only ≈0.06%, `X'γ_0 ~ N(0.685, 1.21²)`), and at (N=200, p=100) `outcome_learner="linear"` fails closed on control-fold rank deficiency even under an oracle propensity (ridge/sieve fit there). Reproduce by generating the equations above at (500, 100) with `DMLDiD(panel=False)` and native learners over any batch of ~20 seeds — the qualitative behavior is seed-robust; the ranges are deliberately qualitative, not seed-pinned. **§4.2.2 (kernel estimation):** `N ∈ {200, 500}`, `D ~ Bernoulli(0.5)`, scalar `X | D ~ N(D, 1)`, `Y⁰(0) = ε_1`, `Y⁰(1) = Y⁰(0) + X + ε_2`, `Y¹(1) = θ_0 + Y⁰(1) + ε_3`, `θ_0 = 3`, `T ~ Bernoulli(0.5)`, same observation rule. True propensity `= sigmoid(X − 1/2)` by Bayes' rule; true `ℓ_20(X) = λ_0(1−λ_0)·X = 0.25·X`. The trend `X + ε_2` with group-imbalanced X violates UNCONDITIONAL parallel trends by construction (the unadjusted 2×2 contrast converges to `θ_0 + 1 = 4`) — the design's covariate-adjustment point. Replicated verbatim in `tests/test_methodology_dml_did.py`. - Finding to preserve for docs: in the repeated-cross-section simulations the orthogonal estimator appears well centered on the truth but visibly NOISIER than Abadie's plug-in at small N (Figures 9-14) — expect larger SEs. Orthogonality removes first-order nuisance-estimation bias (an asymptotic property under the Theorem 1 rates); it does NOT guarantee finite-sample unbiasedness, and histogram centering in simulations cannot establish it. diff --git a/tests/test_methodology_dml_did.py b/tests/test_methodology_dml_did.py index fe966b9f..848cbe11 100644 --- a/tests/test_methodology_dml_did.py +++ b/tests/test_methodology_dml_did.py @@ -343,9 +343,11 @@ def test_misspecified_outcome_correct_propensity(self): # =========================================================================== # # Shared fixed RCS DGP (library-authored, in the spirit of Chang Sec. 4 — -# NOT the paper's own high-dimensional Sec. 4 RCS parameterization, whose -# replication is a tracked TODO row): X ~ N(0, I_2), D ~ -# Bernoulli(sigmoid(0.5 X1 - 0.5 X2)), T ~ Bernoulli(0.5), levels +# NOT the paper's own Sec. 4 RCS parameterization; the paper's Sec. 4.2.2 +# kernel-design DGP is replicated in the "Chang Sec. 4.2.2" section below, +# and the Sec. 4.2.1 ML design is tracked in the narrowed TODO.md row): +# X ~ N(0, I_2), D ~ Bernoulli(sigmoid(0.5 X1 - 0.5 X2)), T ~ +# Bernoulli(0.5), levels # Y = 1 + X1 + 0.5 X2 + T*(0.5 + 0.4 X1) + D + T*D*theta0 + eps. RCS_THETA0 = 3.0 @@ -691,3 +693,197 @@ def test_coverage_sanity(self, ci_params): coverage = hits / n_reps lo_band, hi_band = (0.90, 0.99) if n_reps >= 100 else (0.80, 1.00) assert lo_band <= coverage <= hi_band, coverage + + +# =========================================================================== +# Chang Sec. 4.2.2 kernel-design RCS DGP (paper replication, DML PR-B2) +# =========================================================================== +# +# The paper's OWN repeated-cross-section simulation design (arXiv:1812.10846v3 +# p. 19, Sec. 4.2.2 "Kernel Estimation"), replicated verbatim: D ~ +# Bernoulli(0.5), scalar X | D ~ N(D, 1), Y0(0) = e1, Y0(1) = Y0(0) + X + e2, +# Y1(1) = theta0 + Y0(1) + e3, theta0 = 3, T ~ Bernoulli(0.5), observed +# Y = Y(0) + T*(Y(1) - Y(0)). "N(0, 0.1)" is read as VARIANCE 0.1 +# (sigma = sqrt(0.1)); the printed notation is ambiguous and the adopted +# reading is recorded in docs/methodology/papers/chang-2020-review.md. +# +# The design's trend (X + e2) with group-imbalanced X | D ~ N(D, 1) violates +# UNCONDITIONAL parallel trends by construction: the unadjusted 2x2 group-mean +# DiD converges to theta0 + (E[X|D=1] - E[X|D=0]) = 4, which is what makes it +# a genuine covariate-adjustment fixture. Under our native learners both +# nuisances are correctly specified (true propensity = sigmoid(X - 1/2) by +# Bayes' rule; true outcome nuisance l20(X) = E[(T - lam0) Y | X, D=0] = +# lam0(1 - lam0) X = 0.25 X, from the identity E[(T - lam0) T] = +# lam0(1 - lam0) — the l20 definition is in diff_diff/_dr_scores.py), so this +# section adds paper-DGP faithfulness, not new discriminating power vs the +# incumbent _rcs_frame; the discriminating assertion below compares against +# the design's own confounded contrast instead. The paper estimates this +# design with Gaussian-KERNEL first stages; the fixtures use the library's +# native learners (documented REGISTRY Note). The Sec. 4.2.1 ML design is +# NOT replicable with the bundled unpenalized learners (narrowed TODO.md row). +# The paper publishes histograms only, so every assertion is recovery / +# centering / self-coverage — never a published-number pin. + +CHANG_S422_THETA0 = 3.0 + + +def chang_s422_kernel_frame(n, seed): + """Sec. 4.2.2 verbatim; the EXACT draw order below is load-bearing (the + seed-pinned tolerances in TestChangS422Recovery were measured under it).""" + rng = np.random.default_rng(seed) + D = (rng.uniform(size=n) < 0.5).astype(int) + X = D + rng.standard_normal(n) # X | D ~ N(D, 1) + e1, e2, e3 = rng.normal(0, np.sqrt(0.1), (3, n)) + y00 = e1 + y01 = y00 + X + e2 + y11 = CHANG_S422_THETA0 + y01 + e3 + T = (rng.uniform(size=n) < 0.5).astype(int) + y = np.where(T == 1, np.where(D == 1, y11, y01), y00) + return pd.DataFrame( + { + "unit": np.arange(n), + "time": T + 1, # waves {1, 2}; cohort 2 = treated at period 2 + "first_treat": D * 2, + "y": y, + "x1": X, + } + ) + + +def _fit_s422_expecting_only_a23(est, df): + """Fit recording ALL warnings; tolerate only the deliberate panel=False + Assumption 2.3 lane warning and the occasional propensity-trimming + notice (the design's unbounded-X tails put rare fitted propensities + outside [trim, 1-trim] — the documented overlap quirk; fires on some MC + replicates, not at the pinned recovery seeds). Any OTHER warning + (learner, degenerate cell, inference) fails the fixture instead of + being blanket-suppressed (review round: simplefilter("ignore") weakened + the regression guard).""" + tolerated = ("Assumption 2.3", "will be trimmed") + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + res = est.fit(df, **FIT_KW, covariates=["x1"]) + unexpected = [ + f"{w.category.__name__}: {w.message}" + for w in caught + # both tolerated messages are documented UserWarnings - a matching + # text under a different category is still unexpected + if not (w.category is UserWarning and any(pat in str(w.message) for pat in tolerated)) + ] + assert not unexpected, f"unexpected warnings: {unexpected}" + return res + + +def _unadjusted_did(df): + """The covariate-blind 2x2 group-mean DiD contrast on a Sec.-4.2.2 frame + (shared by the shape pin and every recovery test's discriminating + comparison). Converges to theta0 + 1 on this design.""" + T = (df["time"] == 2).to_numpy() + D = (df["first_treat"] > 0).to_numpy() + m = lambda mask: float(df.loc[mask, "y"].mean()) # noqa: E731 + return (m(T & D) - m(~T & D)) - (m(T & ~D) - m(~T & ~D)) + + +class TestChangS422FrameShape: + """Pin the generator to the paper's distributions: recovery alone passes + under plausible mis-codings (dropped mean shift, dropped trend, SD-vs- + variance misread), so the defining moments are asserted directly.""" + + def test_frame_matches_paper_dgp(self): + df = chang_s422_kernel_frame(200_000, seed=40_000) + x = df["x1"].to_numpy() + T = (df["time"] == 2).to_numpy() + D = (df["first_treat"] > 0).to_numpy() + y = df["y"].to_numpy() + + # X | D ~ N(D, 1); shares 0.5 + assert abs(x[D].mean() - 1.0) < 0.02 and abs(x[~D].mean() - 0.0) < 0.02 + assert abs(x[D].var() - 1.0) < 0.03 and abs(x[~D].var() - 1.0) < 0.03 + assert abs(D.mean() - 0.5) < 0.01 and abs(T.mean() - 0.5) < 0.01 + + # Stationary RCS sampling (Assumption 2.3): the wave draw is + # independent of (D, X) — composition stable across waves. + assert abs(T[D].mean() - T[~D].mean()) < 0.01 + assert abs(x[T].mean() - x[~T].mean()) < 0.02 + assert abs(np.corrcoef(T.astype(float), x)[0, 1]) < 0.01 + + # Control-post regression on X: slope 1, residual var = Var(e1+e2) = 0.2 + mask = T & ~D + slope, intercept = np.polyfit(x[mask], y[mask], 1) + resid = y[mask] - (slope * x[mask] + intercept) + assert abs(slope - 1.0) < 0.03 + assert abs(resid.var() - 0.2) < 0.01 + + # Pre-wave outcome is pure e1: Var = 0.1 + assert abs(y[~T].var() - 0.1) < 0.005 + + # Treated-post residual (OLS of y[T=1,D=1] on x1) = Var(e1+e2+e3) = 0.3 + mask = T & D + slope, intercept = np.polyfit(x[mask], y[mask], 1) + resid = y[mask] - (slope * x[mask] + intercept) + assert abs(resid.var() - 0.3) < 0.015 + + # The design's built-in confounding: the UNADJUSTED contrast converges + # to theta0 + (E[X|D=1] - E[X|D=0]) = 4, NOT theta0. + assert abs(_unadjusted_did(df) - 4.0) < 0.05 + + # Correct-specification facts the recovery narrative rests on: + # true propensity = sigmoid(X - 1/2) (logistic in X)... + logit = LogitLearner().fit(x.reshape(-1, 1), D.astype(float)) + assert abs(logit.intercept_ - (-0.5)) < 0.05 + assert abs(float(np.asarray(logit.coef_).ravel()[0]) - 1.0) < 0.05 + # ...and l20(X) = 0.25 X (OLS of (T - 0.5) y on x1 over controls). + mask = ~D + slope, _ = np.polyfit(x[mask], (T[mask] - 0.5) * y[mask], 1) + assert abs(slope - 0.25) < 0.02 + + +class TestChangS422Recovery: + """Seed-pinned recovery on the paper's design. The absolute bounds are + SEED-PINS at the named data seeds (measured margins 0.242 / 0.095 / + 0.442), NOT distributional bounds (~0.5% of N=500 replicates exceed + 1.25); a failure under the exact chang_s422_kernel_frame draw order is a + frame-implementation discrepancy, not a methodology signal — do not + retune. Each test also asserts the DISCRIMINATING comparison: the + adjusted estimate must land strictly closer to theta0 than the design's + own confounded unadjusted contrast (4.008 / 3.978 / 3.617 at these + seeds), which a covariate-blind estimator cannot do.""" + + @pytest.mark.parametrize( + "n,learner,seed,abs_bound", + [ + (500, "sieve", 40_001, 1.25), + (500, "linear", 40_002, 1.25), + # the paper's small cell; centered but noisier at small N (the + # review's Figures 9-14 finding), hence the wider pin + (200, "linear", 40_003, 2.0), + ], + ) + def test_theta_recovery(self, n, learner, seed, abs_bound): + df = chang_s422_kernel_frame(n, seed=seed) + res = _fit_s422_expecting_only_a23(DMLDiD(outcome_learner=learner, seed=0, panel=False), df) + err = abs(res.overall_att - CHANG_S422_THETA0) + assert err < 4 * res.overall_se, (res.overall_att, res.overall_se) + assert err < abs_bound, (res.overall_att, res.overall_se) + confounded_err = abs(_unadjusted_did(df) - CHANG_S422_THETA0) + assert err < confounded_err, (err, confounded_err) + + +@pytest.mark.slow +class TestChangS422MonteCarloCoverage: + def test_coverage_sanity(self, ci_params): + # Acceptance band CONDITIONAL on the scaled rep count (at ~22 reps the + # tight band fails ~42% of the time at nominal coverage). + n_reps = ci_params.bootstrap(200) + n = 500 + hits = 0 + for rep in range(n_reps): + df = chang_s422_kernel_frame(n, seed=41_000 + rep) + res = _fit_s422_expecting_only_a23( + DMLDiD(outcome_learner="linear", seed=rep, panel=False), df + ) + lo, hi = res.conf_int + hits += int(lo <= CHANG_S422_THETA0 <= hi) + coverage = hits / n_reps + lo_band, hi_band = (0.90, 0.99) if n_reps >= 100 else (0.80, 1.00) + assert lo_band <= coverage <= hi_band, coverage diff --git a/tests/test_methodology_dr_scores.py b/tests/test_methodology_dr_scores.py index 3c4507a9..23cc5b04 100644 --- a/tests/test_methodology_dr_scores.py +++ b/tests/test_methodology_dr_scores.py @@ -198,8 +198,11 @@ def test_native_chang_estimator_matches_doubleml_goldens(self): # # LIBRARY-AUTHORED low-dimensional RCS design in the spirit of Chang Sec. 4 # (Gaussian X, logistic PS, theta0 = 3) — NOT the paper's own Sec. 4 RCS -# parameterization (a high-dimensional p in {100, 300} ML design, not -# extracted into the paper review; replication is a tracked TODO row). +# parameterization. The paper's Sec. 4.2 RCS DGPs are now extracted into +# docs/methodology/papers/chang-2020-review.md; the Sec. 4.2.2 kernel design +# is replicated in tests/test_methodology_dml_did.py ("Chang Sec. 4.2.2" +# section), and the Sec. 4.2.1 ML design (p in {100, 300}) is tracked in the +# narrowed TODO.md row (needs a penalized propensity learner). # Rows i.i.d.: X ~ N(0, I_2); D ~ Bernoulli(sigmoid(0.5 X1 - 0.5 X2)); # T ~ Bernoulli(lam0) independent; levels # Y = 1 + X1 + 0.5 X2 + T*(0.5 + 0.4 X1) + D*1.0 + T*D*theta0 + eps.