diff --git a/CHANGELOG.md b/CHANGELOG.md index 1c0422ca..9922cb00 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,19 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] ### Added +- **`SyntheticControl` ADH-2015 §4 tail diagnostics** (two opt-in `SyntheticControlResults` + methods, closing the last two ADH-2015 §4 checklist items). `regression_weights()` reports the + implied donor weights `W^reg = X0a'(X0a X0a')^{-1} X1a` of the regression counterfactual + (intercept-augmented so they sum to 1 at full row rank) and flags donors outside `[0, 1]` — the + extrapolation an OLS counterfactual incurs but the simplex-constrained synthetic control cannot; + pure linear algebra, min-norm least-squares with a rank-deficient warning. `sparse_synthetic_control()` + exhaustively searches `C(J, l)` size-`l` donor subsets holding `V` fixed at the baseline fit, + reporting how fit / ATT degrade as the synthetic is forced sparse (per-size winner table + + `get_sparse_synthetic_control_gaps()` overlay); the default `sizes=[1,2,3]` sweep skips over-cap + sizes with a warning while an explicitly requested oversized `l` raises (`max_subsets` guard). Both + are purely additive, surface under `estimator_native_diagnostics`, and leave the analytical + inference contract unchanged (`se`/`t_stat`/`p_value`/`conf_int` stay NaN). No behavior change to + any existing fit. - **Rust `demean_map` kernel for FE absorption** (optional backend acceleration). When the Rust extension is available, the method-of-alternating-projections sweeps in `demean_by_groups`/`within_transform` run in a compiled kernel, rayon-parallel across the diff --git a/TODO.md b/TODO.md index ca0e0fcf..fa310747 100644 --- a/TODO.md +++ b/TODO.md @@ -28,7 +28,6 @@ The `Origin` column (Actionable tables) and the `PR` column (Deferred tables) bo | Issue | Location | Origin | Effort | Priority | |-------|----------|--------|--------|----------| -| `SyntheticControl` remaining ADH-2015 §4 items: the regression-weight `W^reg = X_0'(X_0 X_0')^{-1} X_1` extrapolation diagnostic (flag implied OLS weights outside `[0,1]`) and sparse-SC subset search (`l < J`, holding `V` fixed). LOO, in-time placebo, CV `V`-selection, and inverse-variance `V` have landed; these two are the deferred tail. | `synthetic_control.py`, `synthetic_control_results.py` | ADH-2015 | Mid | Low | | `SyntheticControl` conformal (CWZ 2021) extensions: (a) one-sided / signed-`t` variants (§7); (b) covariates in the conformal proxy (`X_jt`, eqs 4/6 — current proxy is outcomes-only); (c) AR / innovation-permutation path (Lemmas 5-7) for time-series proxies. The joint test, pointwise CIs, and average-effect CI have landed. | `conformal.py`, `synthetic_control_results.py` | CWZ-2021 | Heavy | Low | | `ContinuousDiD` CGBS-2024 extensions (matches R `contdid` v0.1.0 deferral set): (a) `covariates=` kwarg; (b) discrete-treatment saturated regression (integer dose currently warned, not routed to per-level coefficients); (c) lowest-dose-as-control per Remark 3.1 when `P(D=0)=0`. REGISTRY `## ContinuousDiD` → Implementation Checklist marks these `[ ]`. | `continuous_did.py` | CGBS-2024 | Heavy | Low | | `ImputationDiD` LOO conservative-variance refinement (BJS 2024 Supp. Appendix A.9) — a finite-sample improvement to the auxiliary-model residuals reducing overfit of `tau_tilde_g` to `epsilon`. Asymptotic Theorem-3 variance is implemented and matches R `didimputation` (which also omits LOO by default). | `imputation.py` | imputation-validation | Mid | Low | diff --git a/diff_diff/diagnostic_report.py b/diff_diff/diagnostic_report.py index a3674f53..b91e7936 100644 --- a/diff_diff/diagnostic_report.py +++ b/diff_diff/diagnostic_report.py @@ -2532,6 +2532,107 @@ def _scm_native(self, r: Any) -> Dict[str, Any]: ), } + # Regression-weight extrapolation diagnostic (ADH 2015 §4): opt-in, surfaced once + # the user has run results.regression_weights() (pure linear algebra — no refits). + if getattr(r, "_regw_df", None) is not None: + regw_status = getattr(r, "_regw_status", None) + if regw_status == "ran": + out["regression_weights"] = { + "status": "ran", + # Headline: how many donors an OLS/regression counterfactual would push + # outside [0,1] (extrapolate). The SC simplex never does. + "n_extrapolating": _to_python_scalar(getattr(r, "_regw_n_extrapolating", None)), + # True if W^reg is a non-unique min-norm solution (not full row rank); + # weight_sum then need not equal 1. + "rank_deficient": bool(getattr(r, "_regw_rank_deficient", False)), + "weight_sum": _to_python_float(getattr(r, "_regw_weight_sum", None)), + } + else: + _regw_reasons = { + "treated_fit_nonconverged": ( + "regression_weights() was run but the treated unit's own SCM fit " + "did not converge at fit time, so the synthetic control it is " + "compared against is not a valid optimum." + ), + "too_few_donors": ( + "regression_weights() was run but fewer than 2 donors are available " + "(W^reg is trivially [1] with a single donor)." + ), + } + out["regression_weights"] = { + "status": "infeasible", + "reason_code": regw_status, + "reason": _regw_reasons.get( + regw_status, "regression_weights() produced no table." + ), + } + else: + out["regression_weights"] = { + "status": "not_run", + "reason": ( + "Call results.regression_weights() to flag donors an OLS/regression " + "counterfactual would weight outside [0,1] (ADH 2015 §4 extrapolation " + "diagnostic; opt-in, no refit)." + ), + } + + # Sparse-SC subset search (ADH 2015 §4): opt-in, surfaced once the user has run + # results.sparse_synthetic_control() (exhaustive subset search, V held fixed). + if getattr(r, "_sparse_df", None) is not None: + sparse_status = getattr(r, "_sparse_status", None) + if sparse_status == "ran": + # Compact per-size summary: how far the ATT / fit move as the synthetic is + # forced sparse (the baseline row is dropped — it is the full-fit reference). + records = r._sparse_df.to_dict("records") + per_size = [ + { + "size": _to_python_scalar(rec["size"]), + "att": _to_python_float(rec["att"]), + "delta_att": _to_python_float(rec["delta_att"]), + "pre_rmspe": _to_python_float(rec["pre_rmspe"]), + "n_failed": _to_python_scalar(rec["n_failed"]), + "status": rec["status"], + } + for rec in records + if rec["status"] != "baseline" + ] + out["sparse_synthetic_control"] = { + "status": "ran", + # Headline: the largest baseline-relative ATT swing across searched sizes. + "max_abs_delta_att": _to_python_float( + getattr(r, "_sparse_max_abs_delta_att", None) + ), + "sizes": per_size, + } + else: + _sparse_reasons = { + "treated_fit_nonconverged": ( + "sparse_synthetic_control() was run but the treated unit's own SCM " + "fit did not converge at fit time, so the baseline ATT is not a " + "valid reference for the sparse deltas." + ), + "too_few_donors": ( + "sparse_synthetic_control() was run but fewer than 2 donors are " + "available (a sparse subset must be smaller than the pool)." + ), + } + out["sparse_synthetic_control"] = { + "status": "infeasible", + "reason_code": sparse_status, + "reason": _sparse_reasons.get( + sparse_status, "sparse_synthetic_control() produced no valid subsets." + ), + } + else: + out["sparse_synthetic_control"] = { + "status": "not_run", + "reason": ( + "Call results.sparse_synthetic_control() to test how the fit / ATT " + "degrade when the synthetic is forced to a few donors (ADH 2015 §4 " + "sparse subset search; opt-in, V held fixed)." + ), + } + # Test-inversion confidence set (Firpo & Possebom 2018 §4): opt-in, surfaced once # the user has run results.confidence_set() (it reuses the in-space placebo # reference set — no refits). The analytical conf_int stays NaN; this is a SEPARATE diff --git a/diff_diff/guides/llms-full.txt b/diff_diff/guides/llms-full.txt index 281c2ef3..5ae5511b 100644 --- a/diff_diff/guides/llms-full.txt +++ b/diff_diff/guides/llms-full.txt @@ -617,7 +617,7 @@ scm.fit( ) -> SyntheticControlResults ``` -**Inference:** NONE analytical — `se`/`t_stat`/`p_value`/`conf_int` are always NaN. `att` is the mean post-period gap. Significance via in-space placebo permutation inference: `results.in_space_placebo()` reassigns treatment to each donor, refits against the other J-1 donors (the real treated unit is excluded from every placebo pool), and sets `placebo_p_value = rank/(n_placebos+1)` from the post/pre RMSPE-ratio. The permutation `placebo_p_value` is a SEPARATE field from the (NaN) `p_value`; `is_significant` stays bound to `p_value`. **ADH-2015 §4 robustness (opt-in, analytical inference unchanged):** `results.leave_one_out()` drops each reportably-weighted donor (weight > 1e-6) and re-fits (per-drop ATT/`delta_att` — large `delta_att` ⇒ single-donor dependence); `results.in_time_placebo()` backdates the intervention and checks for a spurious pre-period gap (TRUNCATE windowing — predictor windows in the held-out region are dropped). **Confidence sets by test inversion (Firpo-Possebom 2018 §4, opt-in, also non-analytical):** `results.test_sharp_null(effect, gamma=0.1)` tests `H_0: α_1t = f(t)` by re-ranking the in-space placebo gaps (no refits; `test_sharp_null(0)` == `placebo_p_value`), and `results.confidence_set(family="constant"|"linear", gamma=0.1)` inverts it into a confidence set for the effect path (constant-effect interval / linear-slope set, strict `p>gamma`), surfaced on `effect_confidence_set` / `get_confidence_set_df()` — the analytical `conf_int` still stays NaN. **Conformal inference (Chernozhukov-Wüthrich-Zhu 2021, opt-in, also non-analytical):** unlike the Firpo path (which re-ranks the cross-unit placebo gaps), the conformal layer fits its OWN time-permutation-invariant constrained-LS proxy (eqs 3–4, no V-matrix) under the null on ALL periods and permutes residuals over time for the single treated unit. `results.conformal_test(effect, q=1, scheme="moving_block")` gives a joint sharp-null p-value (eqs 1–2; `q∈{1,2,∞}`); `results.conformal_confidence_intervals(alpha=0.1)` gives pointwise per-period CIs (Algorithm 1 — each period `t` uses `Z=(pre, t)`, the other post-periods dropped); `results.conformal_average_effect(alpha=0.1)` gives an average-effect CI by collapsing into `T*`-blocks (Appendix A.1). `scheme="moving_block"` (default, serial-dependence-robust) or `"iid"` (finer). Surfaced on `conformal_inference` / `get_conformal_grid_df()`; `conf_int` stays NaN. Predictor periods must lie within the pre window; `post_periods` must be a contiguous suffix cross-checked against `D` (no anticipation). +**Inference:** NONE analytical — `se`/`t_stat`/`p_value`/`conf_int` are always NaN. `att` is the mean post-period gap. Significance via in-space placebo permutation inference: `results.in_space_placebo()` reassigns treatment to each donor, refits against the other J-1 donors (the real treated unit is excluded from every placebo pool), and sets `placebo_p_value = rank/(n_placebos+1)` from the post/pre RMSPE-ratio. The permutation `placebo_p_value` is a SEPARATE field from the (NaN) `p_value`; `is_significant` stays bound to `p_value`. **ADH-2015 §4 robustness (opt-in, analytical inference unchanged):** `results.leave_one_out()` drops each reportably-weighted donor (weight > 1e-6) and re-fits (per-drop ATT/`delta_att` — large `delta_att` ⇒ single-donor dependence); `results.in_time_placebo()` backdates the intervention and checks for a spurious pre-period gap (TRUNCATE windowing — predictor windows in the held-out region are dropped); `results.regression_weights()` computes the implied regression-counterfactual weights `W^reg = X0a'(X0a X0a')^{-1}X1a` (intercept-augmented, sums to 1 at full row rank) and flags donors outside `[0,1]` — the extrapolation an OLS counterfactual incurs but the simplex-constrained SC cannot (pure linear algebra, no refit; min-norm + `_regw_rank_deficient` when `k+1>J`); `results.sparse_synthetic_control(sizes=None, max_subsets=50000)` exhaustively searches `C(J,l)` size-`l` donor subsets holding `V` FIXED at the baseline (default sweep `[1,2,3]` skips over-cap sizes with a warning; an explicitly-requested oversize `l` raises), reporting how fit/ATT degrade as the synthetic is forced sparse (+ `get_sparse_synthetic_control_gaps()`). **Confidence sets by test inversion (Firpo-Possebom 2018 §4, opt-in, also non-analytical):** `results.test_sharp_null(effect, gamma=0.1)` tests `H_0: α_1t = f(t)` by re-ranking the in-space placebo gaps (no refits; `test_sharp_null(0)` == `placebo_p_value`), and `results.confidence_set(family="constant"|"linear", gamma=0.1)` inverts it into a confidence set for the effect path (constant-effect interval / linear-slope set, strict `p>gamma`), surfaced on `effect_confidence_set` / `get_confidence_set_df()` — the analytical `conf_int` still stays NaN. **Conformal inference (Chernozhukov-Wüthrich-Zhu 2021, opt-in, also non-analytical):** unlike the Firpo path (which re-ranks the cross-unit placebo gaps), the conformal layer fits its OWN time-permutation-invariant constrained-LS proxy (eqs 3–4, no V-matrix) under the null on ALL periods and permutes residuals over time for the single treated unit. `results.conformal_test(effect, q=1, scheme="moving_block")` gives a joint sharp-null p-value (eqs 1–2; `q∈{1,2,∞}`); `results.conformal_confidence_intervals(alpha=0.1)` gives pointwise per-period CIs (Algorithm 1 — each period `t` uses `Z=(pre, t)`, the other post-periods dropped); `results.conformal_average_effect(alpha=0.1)` gives an average-effect CI by collapsing into `T*`-blocks (Appendix A.1). `scheme="moving_block"` (default, serial-dependence-robust) or `"iid"` (finer). Surfaced on `conformal_inference` / `get_conformal_grid_df()`; `conf_int` stays NaN. Predictor periods must lie within the pre window; `post_periods` must be a contiguous suffix cross-checked against `D` (no anticipation). **Usage:** @@ -1360,7 +1360,7 @@ Returned by `SyntheticControl.fit()`. | `effect_confidence_set` | `dict \| None` | Test-inversion confidence-set summary (Firpo-Possebom 2018 §4): `{family, parameter, gamma, lower, upper, contiguous, status ("ran"/"empty"/"unbounded"), ...}`; None until `confidence_set()` runs. SEPARATE from the always-NaN analytical `conf_int` (a permutation set at level 1−gamma, possibly a set/unbounded). | | `conformal_inference` | `dict \| None` | CWZ (2021) conformal-inference summary of the most recent run: `{kind ("joint"/"pointwise"/"average"), scheme, status, n_perms, ...}` (joint adds `joint_p_value`/`proxy_converged`; average/pointwise add `alpha` + CI fields); None until a `conformal_*()` method runs. SEPARATE from the always-NaN analytical `conf_int`. | -**Methods:** `in_space_placebo()` (opt-in permutation inference; refits one synthetic control per donor), `get_placebo_df()` (per-unit RMSPE-ratio table incl. the treated row), `leave_one_out()` (ADH-2015 §4 donor robustness; drops each reportably-weighted donor (weight > 1e-6) → per-drop ATT/`delta_att` table) + `get_leave_one_out_df()`/`get_leave_one_out_gaps()`, `in_time_placebo()` (ADH-2015 §4 backdating placebo; reassigns the intervention earlier, TRUNCATE windowing, placebo ATT ~0 if no real pre-effect) + `get_in_time_placebo_df()`/`get_in_time_placebo_gaps()`, `test_sharp_null(effect, gamma=0.1)` (Firpo-Possebom 2018 §4: test a sharp null α_1t=f(t) by re-ranking the in-space placebo gaps — `effect` is a scalar or a post-period array; `test_sharp_null(0)` is identically `placebo_p_value`), `confidence_set(family="constant"|"linear", gamma=0.1, bounds=None, n_grid=200)` (invert that test for a confidence set of the effect path — a constant-effect interval / linear-slope set; strict p>gamma membership; exact piecewise-constant breakpoint inversion when bounds=None, else a fixed grid; `conf_int` stays NaN) + `get_confidence_set_df()`, `conformal_test(effect, q=1, scheme="moving_block", n_iid=10000, seed=None)` (CWZ 2021 joint sharp-null conformal p-value — fits its own constrained-LS proxy under the null on all periods, permutes residuals over time; `q∈{1,2,∞}`), `conformal_confidence_intervals(alpha=0.1, scheme="moving_block", bounds=None, n_grid=100, seed=None)` (pointwise per-period CIs, Algorithm 1 — each period uses Z=(pre, t)), `conformal_average_effect(alpha=0.1, scheme="moving_block", bounds=None, n_grid=200, seed=None)` (average-effect CI by T*-block collapse, Appendix A.1) + `get_conformal_grid_df()`, `summary()`, `print_summary()`, `to_dict()`, `to_dataframe()`, `get_gap_df()`, `get_weights_df()` +**Methods:** `in_space_placebo()` (opt-in permutation inference; refits one synthetic control per donor), `get_placebo_df()` (per-unit RMSPE-ratio table incl. the treated row), `leave_one_out()` (ADH-2015 §4 donor robustness; drops each reportably-weighted donor (weight > 1e-6) → per-drop ATT/`delta_att` table) + `get_leave_one_out_df()`/`get_leave_one_out_gaps()`, `in_time_placebo()` (ADH-2015 §4 backdating placebo; reassigns the intervention earlier, TRUNCATE windowing, placebo ATT ~0 if no real pre-effect) + `get_in_time_placebo_df()`/`get_in_time_placebo_gaps()`, `regression_weights()` (ADH-2015 §4 regression-weight extrapolation diagnostic; intercept-augmented `W^reg`, flags donors outside [0,1]; pure linear algebra) + `get_regression_weights_df()`, `sparse_synthetic_control(sizes=None, max_subsets=50000)` (ADH-2015 §4 sparse subset search; exhaustive `C(J,l)` subsets holding V fixed; default-skip vs explicit-raise cap) + `get_sparse_synthetic_control_df()`/`get_sparse_synthetic_control_gaps()`, `test_sharp_null(effect, gamma=0.1)` (Firpo-Possebom 2018 §4: test a sharp null α_1t=f(t) by re-ranking the in-space placebo gaps — `effect` is a scalar or a post-period array; `test_sharp_null(0)` is identically `placebo_p_value`), `confidence_set(family="constant"|"linear", gamma=0.1, bounds=None, n_grid=200)` (invert that test for a confidence set of the effect path — a constant-effect interval / linear-slope set; strict p>gamma membership; exact piecewise-constant breakpoint inversion when bounds=None, else a fixed grid; `conf_int` stays NaN) + `get_confidence_set_df()`, `conformal_test(effect, q=1, scheme="moving_block", n_iid=10000, seed=None)` (CWZ 2021 joint sharp-null conformal p-value — fits its own constrained-LS proxy under the null on all periods, permutes residuals over time; `q∈{1,2,∞}`), `conformal_confidence_intervals(alpha=0.1, scheme="moving_block", bounds=None, n_grid=100, seed=None)` (pointwise per-period CIs, Algorithm 1 — each period uses Z=(pre, t)), `conformal_average_effect(alpha=0.1, scheme="moving_block", bounds=None, n_grid=200, seed=None)` (average-effect CI by T*-block collapse, Appendix A.1) + `get_conformal_grid_df()`, `summary()`, `print_summary()`, `to_dict()`, `to_dataframe()`, `get_gap_df()`, `get_weights_df()` ### TripleDifferenceResults diff --git a/diff_diff/synthetic_control.py b/diff_diff/synthetic_control.py index 1cfc13f2..731c3a75 100644 --- a/diff_diff/synthetic_control.py +++ b/diff_diff/synthetic_control.py @@ -31,6 +31,7 @@ §SyntheticControl for the deviation/Note labels. """ +import itertools import warnings from dataclasses import dataclass, replace from typing import Any, Dict, List, Optional, Tuple, cast @@ -676,6 +677,21 @@ def fit( # estimator's inputs cannot silently change in_space_placebo() output on an # already-returned results object. (pivots are freshly pivoted here, and the # period/id lists are re-listed below, so those are not caller-aliased.) + # + # Capture the EXACT (X1s, X0s, V) triple that produced the treated unit's final + # donor weights, so the ADH-2015 §4 regression_weights() / sparse_synthetic_control() + # diagnostics can hold V FIXED (unlike leave_one_out / in_time_placebo, which + # re-search V). For nested / custom / inverse_variance the inner solve ran on X1s/X0s + # (standardized or raw per the branch above), so those locals are exactly right. For + # cv the final W is solved on the VALIDATION-window standardized predictors inside + # _outer_solve_V_cv; reuse the balance matrices (Xb1/Xb0, the identical + # validation-window re-aggregation) standardized the same way (== the solver's + # X1s_va/X0s_va, SC _outer_solve_V_cv). The J==1 short-circuit never reaches the cv + # solve, so guard on J > 1 (the diagnostics fail closed on J < 2 anyway). + if self.v_method == "cv" and J > 1: + fit_X1s, fit_X0s, _ = _standardize(Xb1, Xb0, self.standardize) + else: + fit_X1s, fit_X0s = X1s, X0s results._fit_snapshot = _SyntheticControlFitSnapshot( pivots=pivots, specs=list(specs), @@ -701,6 +717,11 @@ def fit( inner_max_iter=self.inner_max_iter, inner_min_decrease=self.inner_min_decrease, v_cv_t0=self.v_cv_t0, + # The V-fixed (X1s, X0s, V) triple captured just above (copied so post-fit + # mutation cannot change the diagnostics' output). + fit_X1s=np.array(fit_X1s, dtype=float, copy=True), + fit_X0s=np.array(fit_X0s, dtype=float, copy=True), + fit_v=np.array(v, dtype=float, copy=True), ) # Persist whether the treated unit's own fit reached a valid optimum — both # the inner Frank-Wolfe weight solve AND (on the nested path) the outer V @@ -2379,5 +2400,127 @@ def _truncate_snapshot_in_time( post_periods=new_post, custom_v=new_custom_v, v_cv_t0=new_v_cv_t0, + # Drop the baseline fit's (X1s, X0s, V) triple: it was built on the FULL pre-period + # predictors, so it must never be paired with these truncated specs. No in-time + # consumer reads it today, but nulling closes the latent trap up front. + fit_X1s=None, + fit_X0s=None, + fit_v=None, ) return snap_mod, dropped + + +# ============================================================================= +# ADH-2015 §4 tail diagnostics (opt-in) — used by SyntheticControlResults +# .regression_weights() / .sparse_synthetic_control() via function-level import +# ============================================================================= + + +def _regression_weights(X1s: np.ndarray, X0s: np.ndarray) -> Tuple[np.ndarray, bool, float]: + """Regression-weight extrapolation diagnostic (ADH 2015 §4, journal pp. 498-499). + + Returns the implied donor weights ``W^reg = X0a'(X0a X0a')^{-1} X1a`` of the + regression counterfactual ``B̂'X_1`` (``B̂ = (X0 X0')^{-1} X0 Y0'``), where ``X0a`` / + ``X1a`` prepend an intercept row of ones to the predictor matrices. With the constant, + ``ι'W^reg = 1`` (under full row rank), so the regression is ALSO a weighting estimator + summing to one — but with UNRESTRICTED weights (can be < 0 or > 1), i.e. it + extrapolates outside the donors' convex hull (the simplex-constrained synthetic control + cannot). Detectable by flagging ``W^reg`` entries outside ``[0, 1]``. + + Computed as the min-norm least-squares solution ``argmin_w ||X0a w - X1a||`` = + ``pinv(X0a) @ X1a`` (equal to the gram-inverse form under full row rank, but avoids + squaring the condition number; ``lstsq`` returns the rank directly). **At full row rank** + the system is exactly solvable, so ``W^reg`` is invariant to per-predictor row scaling + (standardization only divides rows by SD): the result is identical across the + standardized and raw predictor spaces. In the **rank-deficient** case the fit is inexact, + and row scaling reweights the least-squares residual metric, so the min-norm ``W^reg`` is + computed in the captured fit space and can differ across predictor spaces (still a valid + extrapolation witness there, just space-dependent). + + Returns ``(w_reg, rank_deficient, weight_sum)``. ``rank_deficient`` is True when the + augmented predictor matrix is not full ROW rank (``k+1 > J``, or collinear predictors); + the min-norm ``W^reg`` is then still an informative extrapolation witness, but need not + sum to 1 (so ``weight_sum`` may deviate from 1 exactly in that case). + """ + J = X0s.shape[1] + X0a = np.vstack([np.ones((1, J), dtype=float), np.asarray(X0s, dtype=float)]) # (k+1, J) + X1a = np.concatenate([[1.0], np.asarray(X1s, dtype=float).ravel()]) # (k+1,) + w_reg, _, rank, _ = np.linalg.lstsq(X0a, X1a, rcond=None) + rank_deficient = bool(rank < X0a.shape[0]) # not full ROW rank (rows = k+1) + w_reg = np.asarray(w_reg, dtype=float) + return w_reg, rank_deficient, float(np.sum(w_reg)) + + +def _sparse_search_size(snap: "_SyntheticControlFitSnapshot", size: int) -> Dict[str, Any]: + """Exhaustive best size-``size`` donor subset, holding V FIXED at the baseline fit. + + Iterates every ``C(J, size)`` donor subset in deterministic ``itertools.combinations`` + order, solves the inner simplex-constrained weight problem on that subset with the + baseline fit's FIXED ``V`` (``snap.fit_v`` — NOT re-searched, per ADH 2015 footnote 20), + and returns the subset minimizing the pre-period OUTCOME MSPE (strict ``<`` tie-break, so + the first subset in combination order wins ties). Non-converged subset solves are + excluded from the argmin and counted. Only the WINNER's full gap path / ATT is + materialized (the search itself computes only the cheap pre-MSPE per subset). + + Returns a dict always carrying ``n_subsets_evaluated`` / ``n_failed`` / ``all_failed``; + when ``all_failed`` is False it also carries the winner's ``donor_ids`` (tuple), + ``weights`` (dict), ``gap_path``, ``att``, ``pre_rmspe``, ``post_rmspe``, ``rmspe_ratio``. + """ + X1s = snap.fit_X1s + X0s = snap.fit_X0s + v = snap.fit_v + assert X1s is not None and X0s is not None and v is not None + donor_ids = snap.donor_ids + Y = snap.pivots[snap.outcome] + pre_periods = snap.pre_periods + post_periods = snap.post_periods + + # Pre-period treated + donor OUTCOMES for ranking (uniform-outcome MSPE — the SC pre-fit + # metric that the reported pre_rmspe uses; the pre-period gaps equal this residual). + Z1_pre = Y.loc[pre_periods, snap.treated_id].to_numpy(dtype=float) # (T0,) + Z0_pre_all = Y.loc[pre_periods, donor_ids].to_numpy(dtype=float) # (T0, J) + + n_eval = 0 + n_failed = 0 + best_mspe = np.inf + best: Optional[Tuple[List[int], np.ndarray]] = None + # Suppress the per-solve non-convergence / poor-fit warnings across the (up to + # max_subsets) inner solves; the aggregate is surfaced via n_failed instead. + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + for cols in itertools.combinations(range(len(donor_ids)), size): + n_eval += 1 + cols_arr = list(cols) + w_sub, conv = _inner_solve_W( + X1s, X0s[:, cols_arr], v, snap.inner_max_iter, snap.inner_min_decrease + ) + if not conv: + n_failed += 1 + continue + resid = Z1_pre - Z0_pre_all[:, cols_arr] @ w_sub + mspe = float(np.mean(resid**2)) + if mspe < best_mspe: # strict < → first subset in combination order wins ties + best_mspe = mspe + best = (cols_arr, w_sub) + + if best is None: + return {"n_subsets_evaluated": n_eval, "n_failed": n_failed, "all_failed": True} + + cols_arr, w_sub = best + sub_ids = [donor_ids[c] for c in cols_arr] + gap_path = _compute_gap_path(Y, w_sub, snap.treated_id, sub_ids, snap.all_periods) + pre_gaps = np.array([gap_path[p] for p in pre_periods], dtype=float) + post_gaps = np.array([gap_path[p] for p in post_periods], dtype=float) + scale = float(np.max(np.abs(Z1_pre))) if Z1_pre.size else 0.0 + return { + "donor_ids": tuple(sub_ids), + "weights": {sub_ids[i]: float(w_sub[i]) for i in range(len(sub_ids))}, + "gap_path": gap_path, + "att": float(np.mean(post_gaps)), + "pre_rmspe": float(np.sqrt(_mspe(gap_path, pre_periods))), + "post_rmspe": float(np.sqrt(_mspe(gap_path, post_periods))), + "rmspe_ratio": _rmspe_ratio(pre_gaps, post_gaps, scale), + "n_subsets_evaluated": n_eval, + "n_failed": n_failed, + "all_failed": False, + } diff --git a/diff_diff/synthetic_control_results.py b/diff_diff/synthetic_control_results.py index 0b8ba829..18f53129 100644 --- a/diff_diff/synthetic_control_results.py +++ b/diff_diff/synthetic_control_results.py @@ -15,6 +15,7 @@ import warnings from dataclasses import dataclass, field +from math import comb from typing import Any, Dict, List, Optional, Tuple import numpy as np @@ -67,6 +68,16 @@ class _SyntheticControlFitSnapshot: # None → len(pre_periods)//2 default. Carried so in-space/LOO/in-time placebo refits # reproduce the same CV split as the treated fit. v_cv_t0: Optional[int] + # The exact predictor matrices + diagonal V that produced the treated unit's FINAL + # donor weights, in the space the inner solve used (standardized for nested/custom, + # raw for inverse_variance, validation-window standardized for cv). Held FIXED (no + # re-search) by the ADH-2015 §4 regression_weights() / sparse_synthetic_control() + # diagnostics. Optional + nulled in _truncate_snapshot_in_time (a backdated in-time + # snapshot must never pair these full-pre matrices with truncated specs). fit_X1s is + # (k,), fit_X0s is (k, J), fit_v is (k,). + fit_X1s: Optional[np.ndarray] = None + fit_X0s: Optional[np.ndarray] = None + fit_v: Optional[np.ndarray] = None def _validate_conformal_bounds( @@ -359,6 +370,35 @@ def __post_init__(self) -> None: self._conformal_ci_df: Optional[pd.DataFrame] = None self._conformal_grid_df: Optional[pd.DataFrame] = None + # --- ADH 2015 §4 "tail" diagnostics (opt-in, populated by regression_weights() / --- + # --- sparse_synthetic_control()). Both read the fit snapshot's captured (X1s, X0s, V) --- + # --- triple; an unpickled result (snapshot nulled) fails closed in those methods. The --- + # --- small summary tables survive pickling; the per-size sparse gap paths --- + # --- (_sparse_gaps) are panel-derived and nulled by __getstate__. analytical --- + # --- se/t/p/ci stay NaN throughout. + # Regression-weight extrapolation diagnostic: + self._regw_df: Optional[pd.DataFrame] = None + # Status: None (not run), "ran", "treated_fit_nonconverged", "too_few_donors". + self._regw_status: Optional[str] = None + # True if the intercept-augmented predictor matrix was not full ROW rank (a min-norm + # W^reg is reported; the sum-to-1 property then need not hold — see _regw_weight_sum). + self._regw_rank_deficient: bool = False + # Number of donors whose implied regression weight falls outside [0, 1] — the + # extrapolation signal (regression weights are unrestricted, unlike the SC simplex). + self._regw_n_extrapolating: int = 0 + # Σ W^reg — a numerical self-check: ~1 under full row rank (the intercept forces + # ι'W^reg = 1), may deviate from 1 when _regw_rank_deficient. None until run. + self._regw_weight_sum: Optional[float] = None + # Sparse-SC subset search: + self._sparse_df: Optional[pd.DataFrame] = None + # Status: None (not run), "ran", "treated_fit_nonconverged", "too_few_donors". + self._sparse_status: Optional[str] = None + # Headline: max |att_sparse - baseline_att| over the searched sizes. None until run. + self._sparse_max_abs_delta_att: Optional[float] = None + # Per-size winning gap paths {size: {period: gap}} for the overlay plot; panel- + # derived, nulled by __getstate__. + self._sparse_gaps: Optional[Dict[int, Dict[Any, float]]] = None + def __getstate__(self) -> Dict[str, Any]: """Exclude panel-derived internal state from pickling. @@ -379,6 +419,10 @@ def __getstate__(self) -> Dict[str, Any]: # overlay gap accessors raise (re-fit to recompute). state["_loo_gaps"] = None state["_in_time_gaps"] = None + # Sparse-SC per-size winning gap paths are panel-derived (same hazard); the small + # _sparse_df / _regw_df summary tables survive so a round-tripped result still + # reports the diagnostic, but get_sparse_synthetic_control_gaps() raises (re-fit). + state["_sparse_gaps"] = None return state def __setstate__(self, state: Dict[str, Any]) -> None: @@ -391,7 +435,21 @@ def __setstate__(self, state: Dict[str, Any]) -> None: (the "no infeasible refits recorded" state) so a legacy result reports cleanly. """ self.__dict__.update(state) - for _attr, _default in (("n_infeasible", 0), ("_loo_n_infeasible", 0)): + for _attr, _default in ( + ("n_infeasible", 0), + ("_loo_n_infeasible", 0), + # ADH-2015 §4 tail diagnostics (added later): default so a pre-feature pickle + # does not AttributeError in the accessors / DiagnosticReport. + ("_regw_df", None), + ("_regw_status", None), + ("_regw_rank_deficient", False), + ("_regw_n_extrapolating", 0), + ("_regw_weight_sum", None), + ("_sparse_df", None), + ("_sparse_status", None), + ("_sparse_max_abs_delta_att", None), + ("_sparse_gaps", None), + ): if not hasattr(self, _attr): setattr(self, _attr, _default) @@ -1689,6 +1747,450 @@ def get_in_time_placebo_gaps(self) -> pd.DataFrame: ) return pd.DataFrame(rows, columns=["placebo_period", "period", "gap", "phase"]) + # ===================================================================== + # ADH-2015 §4 "tail" diagnostics: regression-weight extrapolation + + # sparse-SC subset search (opt-in; analytical inference unchanged) + # ===================================================================== + + _REGW_COLS = [ + "donor_id", + "w_reg", + "w_sc", + "extrapolates", + "abs_extrapolation", + ] + + def regression_weights(self) -> pd.DataFrame: + """ + Regression-weight extrapolation diagnostic (ADH 2015 §4, journal pp. 498-499). + + Computes the implied donor weights ``W^reg = X0a'(X0a X0a')^{-1} X1a`` of the + REGRESSION counterfactual ``B̂'X_1`` — the same predictor matrices the synthetic + control matched on, augmented with an intercept row of ones. Because a constant is + included, ``ι'W^reg = 1`` (under full row rank), so regression is ALSO a weighting + estimator summing to one — but with UNRESTRICTED weights (can be negative or exceed + 1), i.e. it extrapolates outside the donors' convex hull. The simplex-constrained + synthetic control cannot; comparing the two quantifies how much a regression + counterfactual would have to extrapolate. (In ADH's application regression assigned + negative weights to Greece/Italy/Portugal/Spain.) + + Pure linear algebra — NO solver re-fit — leaving the analytical inference contract + unchanged: ``se``/``t_stat``/``p_value``/``conf_int``/``is_significant`` stay bound + to the NaN analytical ``p_value``. + + Returns + ------- + pandas.DataFrame + One row per donor (all ``J`` donors), sorted by ``abs_extrapolation`` + descending. Columns: ``donor_id``, ``w_reg`` (implied regression weight), + ``w_sc`` (the synthetic-control weight, 0 if below the reporting floor), + ``extrapolates`` (bool: ``w_reg < 0`` or ``w_reg > 1``), ``abs_extrapolation`` + (``max(0, -w_reg, w_reg - 1)`` — the distance outside ``[0, 1]``). + + Raises + ------ + ValueError + If the fit snapshot is unavailable (e.g. this result was unpickled). + + Notes + ----- + When the intercept-augmented predictor matrix is not full ROW rank (``k+1 > J`` — + realistic with the default per-period outcome lags when ``T0 > J`` — or collinear + predictors), the reported ``W^reg`` is the MIN-NORM least-squares solution, a + ``UserWarning`` is emitted, and ``self._regw_rank_deficient`` is set True; it is + still an informative extrapolation witness, but ``Σ W^reg`` + (``self._regw_weight_sum``) need not equal 1 in that case. + """ + if self._fit_snapshot is None: + raise ValueError( + "regression_weights() requires the fit snapshot on the results object. " + "This result appears to have been loaded from serialization (which " + "excludes the snapshot) or produced by an older estimator version. " + "Re-fit to enable the regression-weight extrapolation diagnostic." + ) + from diff_diff.synthetic_control import _regression_weights + + snap = self._fit_snapshot + # Fail closed on a non-converged treated fit for CONSISTENCY with the other ADH-2015 + # diagnostics. (W^reg itself is well-defined regardless — pure linear algebra on the + # captured predictor matrices — so this is a uniform-behaviour POLICY, not a + # correctness necessity: a non-converged treated fit is untrustworthy overall.) + if not self._fit_converged: + warnings.warn( + "regression_weights() skipped: the treated unit's own SCM fit did not " + "converge at fit time, so the synthetic control it is compared against is " + "not a valid optimum. Re-fit with a larger inner_max_iter / more n_starts.", + UserWarning, + stacklevel=2, + ) + self._regw_status = "treated_fit_nonconverged" + self._regw_df = pd.DataFrame([], columns=self._REGW_COLS) + return self._regw_df.copy() + + donor_ids = snap.donor_ids + if len(donor_ids) < 2: + warnings.warn( + "regression_weights() requires at least 2 donors to be informative (with a " + f"single donor W^reg is trivially [1]); only {len(donor_ids)} available. " + "Returning an empty table.", + UserWarning, + stacklevel=2, + ) + self._regw_status = "too_few_donors" + self._regw_df = pd.DataFrame([], columns=self._REGW_COLS) + return self._regw_df.copy() + + # The fit snapshot exists (guarded above) and a converged J>=2 fit always captured the + # matrices, so these are non-None here (narrow for the type checker). + assert snap.fit_X1s is not None and snap.fit_X0s is not None + w_reg, rank_deficient, weight_sum = _regression_weights(snap.fit_X1s, snap.fit_X0s) + if rank_deficient: + warnings.warn( + "regression_weights(): the intercept-augmented predictor matrix is not full " + "row rank (more predictors+intercept than donors, or collinear predictors), so " + "the ADH Gram-inverse form is unavailable; W^reg is the MIN-NORM least-squares " + "solution and need not sum to 1 (and, being an inexact fit, can differ across " + "predictor spaces). It still witnesses extrapolation (weights outside [0, 1]); " + "uniqueness of the least-squares solution depends on the predictor COLUMN rank. " + "Reduce the predictor set or enlarge the donor pool for a full-row-rank W^reg.", + UserWarning, + stacklevel=2, + ) + rows: List[Dict[str, Any]] = [] + for j, d in enumerate(donor_ids): + wj = float(w_reg[j]) + extra = max(0.0, -wj, wj - 1.0) + rows.append( + { + "donor_id": d, + "w_reg": wj, + "w_sc": float(self.donor_weights.get(d, 0.0)), + "extrapolates": bool(wj < 0.0 or wj > 1.0), + "abs_extrapolation": float(extra), + } + ) + # Most-extrapolating donor first (the flagged donors surface at the top). + rows.sort(key=lambda r: r["abs_extrapolation"], reverse=True) + self._regw_rank_deficient = bool(rank_deficient) + self._regw_n_extrapolating = int(sum(1 for r in rows if r["extrapolates"])) + self._regw_weight_sum = float(weight_sum) + self._regw_status = "ran" + self._regw_df = pd.DataFrame(rows, columns=self._REGW_COLS) + return self._regw_df.copy() + + def get_regression_weights_df(self) -> pd.DataFrame: + """ + Get the regression-weight extrapolation table (see :meth:`regression_weights`). + + Survives pickling. Raises if :meth:`regression_weights` has not been run. + + Returns + ------- + pandas.DataFrame + """ + if self._regw_df is None: + raise ValueError("No regression-weight results yet; call regression_weights() first.") + return self._regw_df.copy() + + _SPARSE_COLS = [ + "size", + "donor_ids", + "weights", + "pre_rmspe", + "post_rmspe", + "rmspe_ratio", + "att", + "delta_att", + "n_subsets_evaluated", + "n_failed", + "status", + ] + + def sparse_synthetic_control( + self, sizes: Optional[Any] = None, max_subsets: int = 50000 + ) -> pd.DataFrame: + """ + Sparse synthetic-control subset search (ADH 2015 §4, journal pp. 506-507). + + For each target size ``l < J`` (the donor count), exhaustively searches ALL + ``C(J, l)`` donor subsets — HOLDING ``V`` FIXED at the baseline fit's V (ADH hold V + fixed to make the combinatorial search tractable, footnote 20) — refits the inner + simplex weight solve on each subset, and reports the best-fitting size-``l`` + synthetic (lowest pre-period outcome MSPE). This shows how the fit degrades and the + ATT moves as the synthetic is forced to be sparse (ADH: reducing to ``l = 4, 3, 2`` + degrades fit "moderately", ``l = 1`` much worse — a single-match design close to + DiD). A thin re-run of the validated inner solver: the analytical inference contract + is unchanged (``se``/``t_stat``/``p_value``/``conf_int``/``is_significant`` stay NaN). + + Parameters + ---------- + sizes : int or sequence of int, optional + Target sparsity size(s) ``l``. Default None sweeps ``[1, 2, 3]`` (clipped to + ``l < J``). A DEFAULTED size whose ``C(J, l)`` exceeds ``max_subsets`` is SKIPPED + with a warning (a defaulted call never raises); an EXPLICITLY requested ``l`` with + ``C(J, l) > max_subsets`` raises ValueError instead. Each explicit ``l`` must + satisfy ``1 <= l <= J - 1``. + max_subsets : int, default 50000 + Guard on the exhaustive search. An explicitly requested size exceeding it raises + ValueError with guidance (lower ``l``, curate the donor pool, or raise this cap). + + Returns + ------- + pandas.DataFrame + A ``status="baseline"`` row first (the full fit; ``size`` = the baseline support + count, ``delta_att = 0``), then one ``status="ran"`` row per searched size (or a + ``status="all_subsets_failed"`` row with NaN metrics if every subset of that size + failed to converge). Columns: ``size``, ``donor_ids`` (winning subset, a tuple), + ``weights`` (dict), ``pre_rmspe``, ``post_rmspe``, ``rmspe_ratio``, ``att``, + ``delta_att`` (``att_sparse - full_att``), ``n_subsets_evaluated``, ``n_failed``, + ``status``. + + Raises + ------ + ValueError + If the fit snapshot is unavailable (unpickled result); if ``max_subsets`` is not + a positive integer; if ``sizes`` is an empty sequence; or if an explicitly + requested size is out of range or exceeds ``max_subsets``. + + Notes + ----- + Pre-fit typically degrades as ``l`` shrinks, but strict monotonicity is NOT + guaranteed: subsets are ranked by the uniform-outcome pre-period MSPE while each + subset's weights are V-optimal on the *predictor* objective. The diagnostic's signal + is the degradation of fit and the movement of the ATT as you sparsify. + """ + if self._fit_snapshot is None: + raise ValueError( + "sparse_synthetic_control() requires the fit snapshot on the results " + "object. This result appears to have been loaded from serialization (which " + "excludes the snapshot) or produced by an older estimator version. Re-fit " + "to enable the sparse-SC subset search." + ) + from diff_diff.synthetic_control import _mspe, _sparse_search_size + + # Validate the search budget up front (before any work): a non-positive or non-integer + # max_subsets would otherwise silently mis-behave — e.g. NaN makes every `comb > NaN` + # comparison False, bypassing the cap entirely; <= 0 skips every default size. + if not isinstance(max_subsets, (int, np.integer)) or bool(max_subsets < 1): + raise ValueError(f"max_subsets must be a positive integer, got {max_subsets!r}.") + + snap = self._fit_snapshot + J = len(snap.donor_ids) + + # Baseline row: read DIRECTLY from the full fit (do NOT re-fit) so delta_att=0 is exact. + baseline_row = { + "size": len(snap.weighted_donor_ids), + "donor_ids": tuple(snap.weighted_donor_ids), + "weights": dict(self.donor_weights), + "pre_rmspe": float(self.pre_rmspe), + "post_rmspe": float(np.sqrt(_mspe(self.gap_path, snap.post_periods))), + "rmspe_ratio": float(self.rmspe_ratio), + "att": float(self.att), + "delta_att": 0.0, + "n_subsets_evaluated": 0, + "n_failed": 0, + "status": "baseline", + } + + # Fail closed on a non-converged treated fit: an under-optimized baseline ATT makes + # every sparse delta_att meaningless (mirrors leave_one_out()). + if not self._fit_converged: + warnings.warn( + "sparse_synthetic_control() skipped: the treated unit's own SCM fit did not " + "converge at fit time, so the baseline ATT is not a valid optimum to compare " + "sparse refits against. Re-fit with a larger inner_max_iter / more n_starts.", + UserWarning, + stacklevel=2, + ) + self._sparse_status = "treated_fit_nonconverged" + self._sparse_gaps = {} + self._sparse_df = pd.DataFrame([baseline_row], columns=self._SPARSE_COLS) + return self._sparse_df.copy() + + if J < 2: + warnings.warn( + "sparse_synthetic_control requires at least 2 donors (a sparse subset must " + f"be smaller than the pool); only {J} available. Returning the baseline " + "fit only.", + UserWarning, + stacklevel=2, + ) + self._sparse_status = "too_few_donors" + self._sparse_gaps = {} + self._sparse_df = pd.DataFrame([baseline_row], columns=self._SPARSE_COLS) + return self._sparse_df.copy() + + # Resolve the requested sizes: default sweep [1,2,3] (skip over-cap), or explicit. + explicit = sizes is not None + if sizes is None: + requested = [size for size in (1, 2, 3) if size < J] + else: + # Normalize a scalar to a 1-list; anything else must be a sequence. + raw = [sizes] if isinstance(sizes, (int, np.integer, float)) else None + if raw is None: + try: + raw = list(sizes) + except TypeError: + raise ValueError(f"sizes must be an int or a sequence of ints, got {sizes!r}.") + if not raw: + raise ValueError( + "sizes must be a non-empty int or sequence of ints (got an empty " + "sequence); pass e.g. sizes=[1, 2, 3] or leave sizes=None for the default." + ) + requested = [] + for s in raw: + # Reject bool (an int subclass) and non-integral values: int(2.9) would + # silently truncate to a DIFFERENT requested size than the caller intended. + if isinstance(s, bool) or not isinstance(s, (int, np.integer)): + raise ValueError( + f"sparse_synthetic_control sizes must be integer(s); got {s!r} " + f"(type {type(s).__name__})." + ) + requested.append(int(s)) + + search_sizes: List[int] = [] + for size in requested: + if not (1 <= size <= J - 1): + if explicit: + raise ValueError( + f"sparse_synthetic_control size l={size} is out of range; each size " + f"must satisfy 1 <= l <= J-1 = {J - 1}." + ) + continue # defaulted sizes are pre-clipped; belt-and-suspenders + n_sub = comb(J, size) + if n_sub > max_subsets: + if explicit: + raise ValueError( + f"sparse_synthetic_control size l={size} requires " + f"C({J},{size})={n_sub} inner solves, exceeding " + f"max_subsets={max_subsets}. Lower l, curate the donor pool, or raise " + "max_subsets (the search is exhaustive by design)." + ) + warnings.warn( + f"sparse_synthetic_control: skipping default size l={size} — " + f"C({J},{size})={n_sub} exceeds max_subsets={max_subsets}. Pass " + f"sizes=[{size}] with a larger max_subsets to force it.", + UserWarning, + stacklevel=2, + ) + continue + search_sizes.append(size) + + sparse_rows: List[Dict[str, Any]] = [baseline_row] + sparse_gaps: Dict[int, Dict[Any, float]] = {} + deltas: List[float] = [] + for size in search_sizes: + res = _sparse_search_size(snap, size) + if res["all_failed"]: + warnings.warn( + f"sparse_synthetic_control: all C({J},{size})=" + f"{res['n_subsets_evaluated']} size-{size} subsets failed to converge; " + "row reported with status='all_subsets_failed' and NaN metrics. Re-fit " + "with a larger inner_max_iter / looser inner_min_decrease.", + UserWarning, + stacklevel=2, + ) + sparse_rows.append( + { + "size": size, + "donor_ids": None, + "weights": None, + "pre_rmspe": np.nan, + "post_rmspe": np.nan, + "rmspe_ratio": np.nan, + "att": np.nan, + "delta_att": np.nan, + "n_subsets_evaluated": res["n_subsets_evaluated"], + "n_failed": res["n_failed"], + "status": "all_subsets_failed", + } + ) + continue + delta = float(res["att"]) - float(self.att) + deltas.append(delta) + sparse_gaps[size] = res["gap_path"] + sparse_rows.append( + { + "size": size, + "donor_ids": res["donor_ids"], + "weights": res["weights"], + "pre_rmspe": res["pre_rmspe"], + "post_rmspe": res["post_rmspe"], + "rmspe_ratio": res["rmspe_ratio"], + "att": float(res["att"]), + "delta_att": delta, + "n_subsets_evaluated": res["n_subsets_evaluated"], + "n_failed": res["n_failed"], + "status": "ran", + } + ) + + self._sparse_gaps = sparse_gaps + self._sparse_max_abs_delta_att = max((abs(d) for d in deltas), default=None) + self._sparse_status = "ran" + self._sparse_df = pd.DataFrame(sparse_rows, columns=self._SPARSE_COLS) + return self._sparse_df.copy() + + def get_sparse_synthetic_control_df(self) -> pd.DataFrame: + """ + Get the sparse synthetic-control table (see :meth:`sparse_synthetic_control`). + + Survives pickling. Raises if :meth:`sparse_synthetic_control` has not been run. + + Returns + ------- + pandas.DataFrame + """ + if self._sparse_df is None: + raise ValueError( + "No sparse synthetic-control results yet; call sparse_synthetic_control() first." + ) + return self._sparse_df.copy() + + def get_sparse_synthetic_control_gaps(self) -> pd.DataFrame: + """ + Long-form per-size sparse gap paths, for the overlay ("spaghetti") plot. + + One row per (size, period) for every searched size's winning subset. Columns: + ``size``, ``period``, ``gap``, ``phase`` (``"pre"``/``"post"``) — mirroring + :meth:`get_gap_df`. These per-period paths are panel-derived and are NOT retained + after pickling. + + Returns + ------- + pandas.DataFrame + + Raises + ------ + ValueError + If :meth:`sparse_synthetic_control` has not been run, or if the gap paths were + dropped on pickling (re-fit and re-run to recompute them). + """ + if self._sparse_df is None: + raise ValueError( + "No sparse synthetic-control results yet; call sparse_synthetic_control() first." + ) + if self._sparse_gaps is None: + raise ValueError( + "Sparse synthetic-control gap paths are not retained after pickling " + "(panel-derived); re-run sparse_synthetic_control() on a freshly fitted " + "result to recompute them." + ) + rows: List[Dict[str, Any]] = [] + for size, gap_path in self._sparse_gaps.items(): + for period in list(self.pre_periods) + list(self.post_periods): + if period in gap_path: + phase = "post" if period in self.post_periods else "pre" + rows.append( + { + "size": size, + "period": period, + "gap": gap_path[period], + "phase": phase, + } + ) + return pd.DataFrame(rows, columns=["size", "period", "gap", "phase"]) + # ===================================================================== # Confidence sets by test inversion (Firpo & Possebom 2018, §4) # ===================================================================== diff --git a/docs/api/synthetic_control.rst b/docs/api/synthetic_control.rst index fb76a0ef..420b3fa0 100644 --- a/docs/api/synthetic_control.rst +++ b/docs/api/synthetic_control.rst @@ -30,11 +30,17 @@ from the (NaN) ``p_value``; ``is_significant`` stays bound to ``p_value``. **Robustness diagnostics (ADH 2015 §4, opt-in):** :meth:`~diff_diff.SyntheticControlResults.leave_one_out` drops each reportably-weighted (weight > 1e-6) donor and re-fits (per-drop ATT / ``delta_att`` table — a large ``delta_att`` flags -single-donor dependence), and +single-donor dependence); :meth:`~diff_diff.SyntheticControlResults.in_time_placebo` reassigns the intervention to an earlier pre-date and checks for a spurious gap before the true treatment date (the -backdating placebo; ``placebo_att`` should be ~0). Both re-run the validated solver and -leave the analytical inference fields NaN. +backdating placebo; ``placebo_att`` should be ~0); +:meth:`~diff_diff.SyntheticControlResults.regression_weights` computes the implied +regression-counterfactual donor weights ``W^reg`` (intercept-augmented) and flags those +outside ``[0, 1]`` — the extrapolation an OLS counterfactual would incur but the +simplex-constrained synthetic control cannot (pure linear algebra, no refit); and +:meth:`~diff_diff.SyntheticControlResults.sparse_synthetic_control` exhaustively searches +size-``l`` donor subsets holding ``V`` fixed at the baseline, showing how fit / ATT degrade +as the synthetic is forced sparse. All leave the analytical inference fields NaN. **Confidence sets by test inversion (Firpo & Possebom 2018 §4, opt-in):** :meth:`~diff_diff.SyntheticControlResults.test_sharp_null` tests a sharp null @@ -120,6 +126,11 @@ Results container for synthetic control estimation. ~SyntheticControlResults.in_time_placebo ~SyntheticControlResults.get_in_time_placebo_df ~SyntheticControlResults.get_in_time_placebo_gaps + ~SyntheticControlResults.regression_weights + ~SyntheticControlResults.get_regression_weights_df + ~SyntheticControlResults.sparse_synthetic_control + ~SyntheticControlResults.get_sparse_synthetic_control_df + ~SyntheticControlResults.get_sparse_synthetic_control_gaps ~SyntheticControlResults.test_sharp_null ~SyntheticControlResults.confidence_set ~SyntheticControlResults.get_confidence_set_df diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index 94bcb3ec..76853066 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -2189,9 +2189,11 @@ Classic synthetic control (donor/unit weights only) for a single treated unit, d **Inference:** **No analytical standard error** (Section 2.4) — `se`/`t_stat`/`p_value`/`conf_int` are always NaN. Significance comes from **in-space placebo permutation inference** via `SyntheticControlResults.in_space_placebo()`: reassign treatment to each donor, refit a synthetic control for it, and rank the treated unit's post/pre RMSPE ratio (`rmspe_ratio` = `RMSPE_post / RMSPE_pre` = `sqrt(MSPE_post / MSPE_pre)`) among all units; `placebo_p_value = rank / (n_placebos + 1)`, where `rank = 1 + #{placebos with ratio ≥ treated ratio}` — an **upper-tail rank test on the (unsigned) RMSPE-ratio statistic**, ties counted conservatively via `≥`. Because the ratio squares the gaps it is direction-agnostic: a large ratio signals an effect of *either* sign, so this is NOT a signed/one-directional ("one-sided") hypothesis test on the treatment-effect direction. ADH 2010 §3.4 reports the *MSPE* ratio (the square of `rmspe_ratio`); the two are monotone-equivalent, so the rank and p-value are identical — only the reported statistic's scale differs. `rmspe_ratio` (the treated statistic) is computed at fit time; `placebo_p_value` / `n_placebos` / `n_failed` / `n_infeasible` are populated by the opt-in `in_space_placebo()` call (`n_failed` = solver non-convergences, `n_infeasible` = structural `cv` exclusions — see the split Note below). `get_placebo_df()` returns the per-unit ratio table used for the rank (the per-period placebo gap paths for a "spaghetti" plot are retained internally, not in this summary). -**ADH-2015 §4 robustness diagnostics (opt-in):** Two further diagnostics from Abadie-Diamond-Hainmueller (2015, §4), each a thin re-run of the validated solver — populated only when called and surfaced under `estimator_native_diagnostics` (the analytical inference contract is unchanged; `se`/`t_stat`/`p_value`/`conf_int`/`is_significant` stay bound to the NaN analytical `p_value`): +**ADH-2015 §4 robustness diagnostics (opt-in):** Four further diagnostics from Abadie-Diamond-Hainmueller (2015, §4), each a thin re-run of the validated solver (the regression-weight diagnostic is pure linear algebra, no re-fit) — populated only when called and surfaced under `estimator_native_diagnostics` (the analytical inference contract is unchanged; `se`/`t_stat`/`p_value`/`conf_int`/`is_significant` stay bound to the NaN analytical `p_value`): - **Leave-one-out donor robustness** (`leave_one_out()`): drops each **reportably-weighted** donor and re-fits the treated unit against the reduced pool, returning a per-drop ATT / `delta_att` table (a `status="baseline"` row first, then one row per dropped donor sorted by `|delta_att|`). A large `delta_att` flags single-donor dependence (a single *dominant* donor is still dropped — the others absorb its mass — and its large `delta_att` is the intended signal, not a failure). The reporting stack's headline donor-sensitivity number is `max_abs_delta_att` = `max |delta_att|` over the drops (baseline-relative, so a uniform shift of every drop away from the full-fit ATT is not masked the way a raw ATT range would be). `get_leave_one_out_gaps()` returns the per-drop trajectories for the overlay plot. Fails closed on a non-converged treated fit or `< 2` donors. - **In-time (backdating) placebo** (`in_time_placebo()`): reassigns the intervention to an earlier pre-date `t_f`, re-fits using ONLY pre-`t_f` information (TRUNCATE convention — see Note), and reports the placebo "effect" over the held-out window `[t_f, T0)` — ~0 if there is no real pre-period effect (ADH 2015 Fig. 4, German reunification backdated to 1975). Sweeps every feasible interior pre-date by default (≥2 pre-fake + ≥1 post-fake); an explicit post-period / non-pre date raises, a valid-but-dimensionally-infeasible date yields a `status="infeasible"` row (no raise). +- **Regression-weight extrapolation diagnostic** (`regression_weights()`, journal pp. 498-499): computes the implied donor weights `W^reg = X0a'(X0a X0a')^{-1} X1a` of the regression counterfactual `B̂'X_1` on the baseline fit's predictor matrices with an **intercept row prepended** (so `ι'W^reg = 1` under full row rank — matching ADH's "if a constant is included"). Regression is then also a weighting estimator summing to one, but with UNRESTRICTED weights: `w_reg < 0` or `> 1` flags donors it extrapolates on, outside the donors' convex hull the simplex-constrained SC never leaves (in ADH's application regression assigned negative weights to Greece/Italy/Portugal/Spain). Returns a per-donor `w_reg` / `w_sc` / `extrapolates` / `abs_extrapolation` table sorted by extrapolation magnitude; the reporting headline is `n_extrapolating`. Pure linear algebra (no re-fit), computed as the min-norm least-squares solution. **At full row rank** the system is exactly solvable, so `W^reg` is invariant to per-predictor row scaling — identical across the standardized (nested/custom) and raw (inverse_variance) predictor spaces (differing only for `cv`, which matches on validation-window predictors). Fails closed (for consistency with the other diagnostics) on a non-converged treated fit or `< 2` donors. **Note:** the sum-to-one property AND the row-scaling invariance both hold only at full ROW rank; when the intercept-augmented predictor matrix is rank-deficient (`k+1 > J` — e.g. the default per-period outcome lags when `T0 > J` — or collinear predictors) the reported `W^reg` is the MIN-NORM least-squares solution (a `UserWarning` fires; `_regw_rank_deficient` is set), `Σ W^reg` need not equal 1, and — because the inexact fit's residual metric is reweighted by row scaling — the min-norm `W^reg` is computed in the captured fit space and can itself differ across predictor spaces (still an informative extrapolation witness; within a fixed space, uniqueness of the least-squares solution depends on the predictor COLUMN rank — an overdetermined full-column-rank case has a unique solution that merely need not sum to 1). +- **Sparse-SC subset search** (`sparse_synthetic_control(sizes=None, max_subsets=50000)`, journal pp. 506-507): for each target size `l < J`, EXHAUSTIVELY searches all `C(J, l)` donor subsets **holding V FIXED at the baseline fit's V** (ADH hold V fixed to make the combinatorial search tractable, footnote 20 — this reuses the captured fit-time predictor matrices + V, unlike `leave_one_out()`/`in_time_placebo()` which re-search V), refits the inner simplex solve per subset, and reports the best size-`l` synthetic (lowest pre-period outcome MSPE). Shows how the fit degrades and the ATT moves as the synthetic is forced sparse (`l = 4, 3, 2` degrade "moderately", `l = 1` much worse — a single-match design close to DiD, footnote 23). A `status="baseline"` row (the full fit) first, then one row per size with the winning subset's ids / weights / `pre_rmspe` / `att` / `delta_att`; `get_sparse_synthetic_control_gaps()` returns the per-size winner paths for the overlay. `sizes=None` sweeps `[1, 2, 3]` (clipped to `l < J`); a DEFAULTED size whose `C(J, l)` exceeds `max_subsets` is **skipped with a warning** (a defaulted call never raises), while an **explicitly** requested `l` with `C(J, l) > max_subsets` raises `ValueError` (the user-facing "exhaustive or raise" contract — no silent approximation). Non-converged subset solves are excluded and counted (`n_failed`); a size with no converged subset gets `status="all_subsets_failed"`. Fails closed (warn + baseline-only) on a non-converged treated fit or `< 2` donors. **Note:** pre-fit typically degrades as `l` shrinks but strict monotonicity is NOT guaranteed — subsets are ranked by the uniform-outcome pre-period MSPE while each subset's weights are V-optimal on the predictor objective. **Confidence sets by test inversion (Firpo & Possebom 2018, §4, opt-in):** A confidence set for the treatment-effect path, built ON TOP of the in-space placebo and surfaced under `estimator_native_diagnostics` (the analytical inference contract is unchanged — see the non-analytical Note below). Under a common-effect sharp null `H_0^f: α_1t = f(t)` (Eq 11) the donor synthetic controls and the pre-period MSPE denominators do not change — only the post-period gaps shift by `f(t)` — so the test is a pure **re-ranking** of the gap paths `in_space_placebo()` already computed (no synthetic-control refits): - **`test_sharp_null(effect, gamma=0.1)`** forms the modified RMSPE ratio `RMSPE^f_j = sqrt(mean_post((α̂_jt − f(t))²) / pre_denom_j)` for every unit (Eq 12) and the permutation p-value `p^f = (1 + #{converged placebos j : RMSPE^f_j ≥ RMSPE^f_1}) / (n_ref + 1)` (Eq 13 at `φ=0`, `v=(1,…,1)` — the equal-weights benchmark). `effect` is a scalar (a constant-in-time effect) or a length-`n_post` array (an arbitrary post-period path `f(t)` — e.g. an intervention cost path or a theory-predicted shape). At `f≡0` this is identically the in-space `placebo_p_value` (Eq 5 = Eq 13), held bit-for-bit by reusing each unit's floored pre-period denominator persisted from the placebo run (the pre window is `f`-free, so the denominator is grid-invariant; the floor scale is **per-unit** `max|Z1_j|`). @@ -2230,7 +2232,7 @@ Classic synthetic control (donor/unit weights only) for a single treated unit, d - **Note (in-time placebo windowing — TRUNCATE):** ADH 2015 §4 says to re-estimate the in-time placebo "with the same predictors lagged accordingly." Because `diff_diff`'s predictor specs reference **absolute** periods, the in-time placebo re-cuts them by TRUNCATION: pre-period-outcome predictors become the pre-`t_f` outcomes, and covariate / special-predictor windows are intersected with the pre-`t_f` window; a window lying ENTIRELY in the held-out region `[t_f, T0)` is **dropped** (surfaced in the `n_dropped_specs` column + an aggregated warning), and `custom_v` is subset in lockstep with the surviving specs. For an outcome-predictor fit (the R-anchorable case) TRUNCATE is identical to ADH's "lag" — both equal a manual `Synth::synth` re-run with `time.optimize.ssr` cut at `t_f`. The held-out window never enters the fit (the placebo's `all_periods` is the pre-fake + post-fake span; the true post-treatment periods are excluded entirely), so there is no "peeking." This concrete convention is NOT spelled out in ADH 2015 (which gives only the qualitative "lag accordingly"). - **Note (in-time placebo requires ≥2 pre-fake periods):** the in-time placebo treats a date with fewer than 2 pre-fake periods as `status="infeasible"` (the default sweep starts at the 3rd pre-period). This is DELIBERATELY stricter than the base estimator's `T0 ≥ 1` allowance (which permits a single-pre-period fit but warns that nested-`V` selection is unreliable): an auto-swept placebo date with a single pre-fake period is a trivially-matchable, non-credible pre-fit, so it is dropped rather than surfaced as a `ran` placebo (mirrors `SyntheticDiD.in_time_placebo`'s `i ≥ 2` rule). A date whose surviving `custom_v` has zero mass after truncation is likewise infeasible (not a convergence failure). - **Note (leave-one-out weight floor):** ADH 2015 §4 leave-one-out omits "each donor that received positive weight." This implementation drops each donor with **reportable** weight — above the `1e-6` interpretability floor (`synthetic_control._MIN_REPORT_WEIGHT`), i.e. exactly the donors in `donor_weights` — rather than every strictly-positive weight. A donor with `0 < w ≤ 1e-6` is numerical dust whose removal moves the ATT by ~its weight (its `delta_att` would be ~0, an uninformative row), and the floor keeps the LOO table aligned with the reported donor support. The drop-set is **frozen at fit time** on the fit snapshot (`weighted_donor_ids`), so `leave_one_out()` is immune to post-fit mutation of the presentation-level `donor_weights` dict. -- **Note (ADH-2015 diagnostics validation):** R `Synth` has **no** in-time-placebo or leave-one-out function (verified against its full CRAN function index; `SCtools` adds only the *in-space* placebo battery, `scpi` only prediction-interval uncertainty), so there is no canonical R *output* to match for these diagnostics — in R they are hand-rolled by re-running `dataprep()`+`synth()`. They are validated instead by (a) the solver's existing Basque R parity (above), and (b) deterministic **self-consistency** tests proving each diagnostic equals a from-scratch `synthetic_control()` fit on the equivalent sub-problem — `leave_one_out()` drop-`d` == a fit on the donor pool minus `d`; `in_time_placebo([t_f])` == a fit on the backdated/truncated panel — both via a fixed `custom_v` (match to 1e-7). The remaining deferred ADH-2015 items (regression-weight `W^reg` extrapolation diagnostic, sparse-SC subset search) are tracked in `TODO.md`. +- **Note (ADH-2015 diagnostics validation):** R `Synth` has **no** in-time-placebo or leave-one-out function (verified against its full CRAN function index; `SCtools` adds only the *in-space* placebo battery, `scpi` only prediction-interval uncertainty), so there is no canonical R *output* to match for these diagnostics — in R they are hand-rolled by re-running `dataprep()`+`synth()`. They are validated instead by (a) the solver's existing Basque R parity (above), and (b) deterministic **self-consistency** tests proving each diagnostic equals a from-scratch `synthetic_control()` fit on the equivalent sub-problem — `leave_one_out()` drop-`d` == a fit on the donor pool minus `d`; `in_time_placebo([t_f])` == a fit on the backdated/truncated panel — both via a fixed `custom_v` (match to 1e-7). The two §4-tail diagnostics are likewise R-anchor-free (R `Synth` has neither): `regression_weights()` is validated by a **numpy oracle** re-implementing `W^reg = X0a'(X0a X0a')^{-1}X1a` on hand-built matrices (incl. the full-rank sum-to-one property, the rank-deficient min-norm branch, and row-scaling invariance across `v_method` spaces); `sparse_synthetic_control()` by **self-consistency** — its exhaustive size-`l` winner (and the winner's weights) match an independent brute-force enumeration using the SAME fixed baseline `V`, which also confirms `V` is held fixed rather than re-searched. - **Note (conformal proxy ≠ ADH V-matrix — deliberate):** the CWZ conformal layer fits its OWN counterfactual proxy — the canonical constrained-LS synthetic control on **raw outcomes over all periods** (eqs 3–4, no V-matrix) — NOT the headline ADH V-matrix weights (which match on pre-period predictors only). This is required, not incidental: CWZ's exactness theory (Lemma 1; Appendix D exchangeability under the null) holds for a **time-permutation-invariant** estimator, which the ADH pre-period V-fit is not (it treats pre and post asymmetrically). So the conformal counterfactual / `point_estimate` can differ from the headline `att`, and is reported as a separate object. - **Note (conformal permutation floor — `1/|Π|`, distinct from Firpo's `1/(J+1)`):** the conformal p-value is `(1/|Π|)·#{π : S(û_π) ≥ S(û)}` (eq 2). The permutation set `Π` **includes the identity** (`S(û_id)=S(û)`, counted under `≥`), so `p̂ ≥ 1/|Π|` automatically — there is **no extra `+1`** (unlike the cross-unit placebo / Firpo `(1+n)/(n+1)`, where the treated unit is not a member of the placebo reference set). `|Π| = T` (moving-block, joint), `T0+1` (pointwise sub-series), or `T/T*` (average-effect blocks), capped at `n_iid` for the i.i.d. scheme. - **Note (conformal per-period CI drops the other post-periods — paper-sourced):** for a pointwise CI of period `t`, CWZ §2.2 defines the data under the null as `Z = (Z_1,…,Z_{T0}, Z_t)` — the pre-periods plus ONLY period `t`. The other post-periods are dropped (not plugged in, not zeroed), making each per-period CI a clean `T*=1` conformal test on the `(T0+1)`-length sub-series. (Confirmed against arXiv:1712.09089v10 §2.2 + Algorithm 1.) The grid is centred on a **pre-only** proxy fit (predict the post slot from the pre slots) — the unconditional all-slot fit would soak the effect into the weights and bias the naïve residual toward 0, so it is not used as the centre. @@ -2253,7 +2255,8 @@ Classic synthetic control (donor/unit weights only) for a single treated unit, d - [x] In-time (backdating) placebo (`in_time_placebo()`, ADH 2015 §4): TRUNCATE windowing (drop held-out-window predictors + lockstep `custom_v` subset), feasible-date sweep, fail-closed. - [x] Confidence sets by test inversion (`test_sharp_null()` + `confidence_set()`, Firpo & Possebom 2018 §4): sharp-null `RMSPE^f` re-ranking of the in-space placebo gaps (Eqs 12–13) + constant/linear one-parameter sets (Eqs 14/16/18) with the strict `p^f > γ` boundary, EXACT piecewise-constant breakpoint inversion (no shape assumption; isolated/disjoint/unbounded sets handled), and fail-closed unbounded/empty/non-contiguous handling. *Deferred:* sensitivity weights (φ≠0), the general-θ menu (Eq 19), one-sided (§7), multiple-outcome/treated (§6). - [x] Conformal inference (`conformal_test()` + `conformal_confidence_intervals()` + `conformal_average_effect()`, Chernozhukov-Wüthrich-Zhu 2021): own constrained-LS proxy under the null on all periods (eqs 3–4, no V-matrix) + `S_q` statistic (`q=1,2,∞`) + permutation p-value (eq 2, `1/|Π|` floor) over moving-block (`Π_→`) / i.i.d. (`Π_all`) schemes; joint sharp-null test (eqs 1–2), pointwise per-period CIs (Algorithm 1, `Z=(pre,t)`), and the average-effect block-collapse CI (Appendix A.1); fail-closed grid_limited/empty/indeterminate handling; analytical `conf_int` stays NaN. *Deferred:* one-sided (§7), covariates folded into the proxy, AR/innovation-permutation path (Lemmas 5–7) — see `TODO.md`. -- [ ] *Deferred (ADH 2015):* regression-weight `W^reg` extrapolation diagnostic, sparse-SC subset search (see `TODO.md`). +- [x] Regression-weight `W^reg` extrapolation diagnostic (`regression_weights()`, ADH 2015 §4): intercept-augmented `W^reg = X0a'(X0a X0a')^{-1}X1a`, flag donors outside `[0,1]`; min-norm + rank-deficient handling; pure linear algebra, analytical inference unchanged. +- [x] Sparse-SC subset search (`sparse_synthetic_control()`, ADH 2015 §4): exhaustive `C(J,l)` subset search holding `V` fixed at the baseline; default-skip vs explicit-raise `max_subsets` guard; per-size winner table + overlay gaps. - [x] Predictor-leakage, absorbing-suffix/no-anticipation, empty-window, duplicate-label, and inner-non-convergence validation gates. --- diff --git a/docs/methodology/REPORTING.md b/docs/methodology/REPORTING.md index 67629c58..e4599ad9 100644 --- a/docs/methodology/REPORTING.md +++ b/docs/methodology/REPORTING.md @@ -268,10 +268,12 @@ a library setting. routes parallel-trends to the `scm_fit` analogue (`pre_rmspe`, verdict `design_enforced_pt`) and surfaces `pre_rmspe`, donor-weight concentration, the in-space placebo permutation p-value, the - ADH-2015 leave-one-out (`leave_one_out`) and in-time placebo - (`in_time_placebo`) blocks, the Firpo-Possebom (2018) test-inversion - `confidence_set`, and the Chernozhukov-Wüthrich-Zhu (2021) - `conformal_inference` block (joint / pointwise / average) under + ADH-2015 §4 leave-one-out (`leave_one_out`), in-time placebo + (`in_time_placebo`), regression-weight extrapolation + (`regression_weights`) and sparse-SC subset-search + (`sparse_synthetic_control`) blocks, the Firpo-Possebom (2018) + test-inversion `confidence_set`, and the Chernozhukov-Wüthrich-Zhu + (2021) `conformal_inference` block (joint / pointwise / average) under `estimator_native_diagnostics` — each is populated only when the caller has already run the corresponding opt-in method (DR never triggers a refit loop implicitly; otherwise a diff --git a/docs/methodology/papers/abadie-diamond-hainmueller-2015-review.md b/docs/methodology/papers/abadie-diamond-hainmueller-2015-review.md index 12127999..517e17c3 100644 --- a/docs/methodology/papers/abadie-diamond-hainmueller-2015-review.md +++ b/docs/methodology/papers/abadie-diamond-hainmueller-2015-review.md @@ -85,7 +85,8 @@ If a constant is included, `ι'W^{reg}=1` — i.e., regression is *also* a weigh - [ ] In-time placebo (date reassignment with predictor lagging). - [ ] Leave-one-out donor robustness (drop each positively-weighted donor). - [ ] Post/pre-RMSPE-ratio permutation p-value `(#{r_j ≥ r_1})/(J+1)`. -- [ ] Regression-weight (`W^{reg}`) extrapolation diagnostic (flag weights outside `[0,1]`). +- [x] Regression-weight (`W^{reg}`) extrapolation diagnostic (flag weights outside `[0,1]`) — implemented as `SyntheticControlResults.regression_weights()`. +- [x] Sparse-SC subset search (`l k+1 = 4 <= J so W^reg is full row rank. + df, _years, _T0 = _make_panel(n_donors=6, T=8, T0=3, seed=1) + res = _fit_iv(df) + assert res._fit_converged + snap = res._fit_snapshot + X1s, X0s = snap.fit_X1s, snap.fit_X0s + _k, J = X0s.shape + X0a = np.vstack([np.ones((1, J)), X0s]) + X1a = np.concatenate([[1.0], X1s.ravel()]) + # ADH's exact formula W^reg = X0a'(X0a X0a')^{-1} X1a (independent of the impl's lstsq). + w_expected = X0a.T @ np.linalg.solve(X0a @ X0a.T, X1a) + tab = res.regression_weights().set_index("donor_id") + got = np.array([tab.loc[d, "w_reg"] for d in snap.donor_ids]) + np.testing.assert_allclose(got, w_expected, atol=1e-10) + # Full row rank -> intercept forces the weights to sum to 1. + assert res._regw_rank_deficient is False + assert res._regw_status == "ran" + assert abs(res._regw_weight_sum - 1.0) < 1e-8 + # Flag columns are internally consistent with w_reg. + full = res.get_regression_weights_df() + for _, row in full.iterrows(): + w = row["w_reg"] + assert row["extrapolates"] == bool(w < 0.0 or w > 1.0) + assert row["abs_extrapolation"] == pytest.approx(max(0.0, -w, w - 1.0)) + assert res._regw_n_extrapolating == int(full["extrapolates"].sum()) + # The FULL analytical inference contract is untouched by the diagnostic (all NaN; the + # permutation-only significance stays off the analytical fields). + res.sparse_synthetic_control(sizes=[2]) + assert np.isnan(res.se) and np.isnan(res.p_value) and np.isnan(res.t_stat) + assert np.isnan(res.conf_int[0]) and np.isnan(res.conf_int[1]) + assert not res.is_significant + + +def test_regression_weights_invariant_to_predictor_row_scaling(): + # At FULL ROW RANK (k=T0=3 predictors < J=6 donors, exact fit) W^reg is invariant to + # per-predictor row scaling: a custom (standardized) fit and an inverse_variance (raw) fit + # on the same data give identical implied regression weights. (The invariance holds only + # under full row rank — see test_regression_weights_rank_deficient_warns_and_min_norm and + # the REGISTRY note; in the rank-deficient min-norm case row scaling can change W^reg.) + df, _years, T0 = _make_panel(n_donors=6, T=8, T0=3, seed=2) + r_std = SyntheticControl(v_method="custom", custom_v=np.ones(T0), seed=0, **_FAST).fit( + df, "y", "treated", "unit", "year" + ) + r_raw = _fit_iv(df) + t_std = r_std.regression_weights().set_index("donor_id")["w_reg"] + t_raw = r_raw.regression_weights().set_index("donor_id")["w_reg"] + assert not r_std._regw_rank_deficient and not r_raw._regw_rank_deficient # full-rank regime + order = sorted(t_std.index) + np.testing.assert_allclose( + t_std.reindex(order).to_numpy(), t_raw.reindex(order).to_numpy(), atol=1e-10 + ) + + +def test_regression_weights_rank_deficient_warns_and_min_norm(): + # k = T0 = 6 predictors, J = 3 donors -> k+1 = 7 > J: not full row rank. + df, _years, _T0 = _make_panel(n_donors=3, T=10, T0=6, seed=1) + res = _fit_iv(df) + with pytest.warns(UserWarning, match="not full row rank"): + tab = res.regression_weights() + assert res._regw_rank_deficient is True + assert res._regw_status == "ran" + assert len(tab) == 3 # all donors still reported + + +def test_regression_weights_fail_closed_on_unpickled(): + df, _years, _T0 = _make_panel(n_donors=5, T=8, T0=4, seed=1) + res = _fit_iv(df) + res2 = pickle.loads(pickle.dumps(res)) + with pytest.raises(ValueError, match="fit snapshot"): + res2.regression_weights() + + +def test_regression_weights_too_few_donors(): + df, _years, _T0 = _make_panel(n_donors=1, T=8, T0=4, seed=1) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = SyntheticControl(seed=0, **_FAST).fit(df, "y", "treated", "unit", "year") + with pytest.warns(UserWarning, match="at least 2 donors"): + tab = res.regression_weights() + assert res._regw_status == "too_few_donors" + assert tab.empty + + +def test_sparse_self_consistency_and_holds_v_fixed(): + from diff_diff.synthetic_control import _inner_solve_W + + df, _years, _T0 = _make_panel(n_donors=6, T=8, T0=4, seed=1) + res = _fit_iv(df) + snap = res._fit_snapshot + tab = res.sparse_synthetic_control(sizes=[2]) + row = tab[tab["status"] == "ran"].iloc[0] + won_ids = row["donor_ids"] + won_weights = row["weights"] + + # Independent brute-force best size-2 subset, using the FIXED baseline V (snap.fit_v) — + # if the diagnostic re-searched V instead of holding it fixed, this would diverge. + import itertools + + Y = snap.pivots[snap.outcome] + Z1_pre = Y.loc[snap.pre_periods, snap.treated_id].to_numpy(float) + Z0_pre = Y.loc[snap.pre_periods, snap.donor_ids].to_numpy(float) + best, best_mspe = None, np.inf + for cols in itertools.combinations(range(len(snap.donor_ids)), 2): + w, conv = _inner_solve_W( + snap.fit_X1s, + snap.fit_X0s[:, list(cols)], + snap.fit_v, + snap.inner_max_iter, + snap.inner_min_decrease, + ) + if not conv: + continue + m = float(np.mean((Z1_pre - Z0_pre[:, list(cols)] @ w) ** 2)) + if m < best_mspe: + best_mspe, best = m, (cols, w) + exp_ids = tuple(snap.donor_ids[c] for c in best[0]) + assert won_ids == exp_ids + got_w = np.array([won_weights[i] for i in won_ids]) + np.testing.assert_allclose(got_w, best[1], atol=1e-10) + # pre_rmspe reported equals sqrt of the winning subset's pre-MSPE. + assert row["pre_rmspe"] == pytest.approx(np.sqrt(best_mspe), abs=1e-10) + + +def test_sparse_l1_picks_best_single_donor(): + df, _years, _T0 = _make_panel(n_donors=5, T=8, T0=4, seed=3) + res = _fit_iv(df) + snap = res._fit_snapshot + tab = res.sparse_synthetic_control(sizes=[1]) + won = tab[tab["status"] == "ran"].iloc[0]["donor_ids"] + assert len(won) == 1 + # l=1 forces w=[1], so the synthetic IS that donor's series; the winner is the donor + # whose own pre-period outcomes best match the treated unit's. + Y = snap.pivots[snap.outcome] + Z1 = Y.loc[snap.pre_periods, snap.treated_id].to_numpy(float) + mspes = { + d: float(np.mean((Z1 - Y.loc[snap.pre_periods, d].to_numpy(float)) ** 2)) + for d in snap.donor_ids + } + assert won[0] == min(mspes, key=mspes.get) + + +def test_sparse_explicit_oversize_raises(): + df, _years, _T0 = _make_panel(n_donors=6, T=8, T0=4, seed=1) + res = _fit_iv(df) + with pytest.raises(ValueError, match="exceeding max_subsets"): + res.sparse_synthetic_control(sizes=3, max_subsets=5) # C(6,3)=20 > 5 + + +def test_sparse_default_skips_over_cap_without_raising(): + # J=8: C(8,1)=8 (ok), C(8,2)=28 and C(8,3)=56 (> cap 10) -> the two large defaults skip. + df, _years, _T0 = _make_panel(n_donors=8, T=8, T0=4, seed=1) + res = _fit_iv(df) + with pytest.warns(UserWarning, match="skipping default size"): + tab = res.sparse_synthetic_control(max_subsets=10) + assert set(tab[tab["status"] == "ran"]["size"]) == {1} + assert res._sparse_status == "ran" # skipped, NOT raised + + +def test_sparse_baseline_row_is_exact(): + df, _years, _T0 = _make_panel(n_donors=5, T=8, T0=4, seed=1) + res = _fit_iv(df) + tab = res.sparse_synthetic_control(sizes=[2]) + base = tab[tab["status"] == "baseline"].iloc[0] + assert base["delta_att"] == 0.0 + assert base["att"] == pytest.approx(res.att, abs=1e-12) + assert base["size"] == len(res._fit_snapshot.weighted_donor_ids) + + +def test_sparse_fail_closed_on_unpickled(): + df, _years, _T0 = _make_panel(n_donors=5, T=8, T0=4, seed=1) + res = _fit_iv(df) + res2 = pickle.loads(pickle.dumps(res)) + with pytest.raises(ValueError, match="fit snapshot"): + res2.sparse_synthetic_control() + # The summary table survives, but the panel-derived gap accessor fails closed. + res.sparse_synthetic_control(sizes=[2]) + res3 = pickle.loads(pickle.dumps(res)) + assert res3.get_sparse_synthetic_control_df() is not None # small table survives + with pytest.raises(ValueError, match="not retained after pickling"): + res3.get_sparse_synthetic_control_gaps() + + +def test_sparse_too_few_donors(): + df, _years, _T0 = _make_panel(n_donors=1, T=8, T0=4, seed=1) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + res = SyntheticControl(seed=0, **_FAST).fit(df, "y", "treated", "unit", "year") + with pytest.warns(UserWarning, match="at least 2 donors"): + tab = res.sparse_synthetic_control() + assert res._sparse_status == "too_few_donors" + assert list(tab["status"]) == ["baseline"] + + +def test_adh_tail_diagnostics_surface_in_diagnostic_report(): + df, _years, T = _make_panel(n_donors=5, T=8, T0=4, seed=1) + res = _fit_iv(df) + # not_run stubs before the opt-in methods are called. + nat0 = DiagnosticReport(res).to_dict()["estimator_native_diagnostics"] + assert nat0["regression_weights"]["status"] == "not_run" + assert nat0["sparse_synthetic_control"]["status"] == "not_run" + res.regression_weights() + res.sparse_synthetic_control(sizes=[2, 3]) + gaps = res.get_sparse_synthetic_control_gaps() + assert set(gaps["size"]) == {2, 3} + assert len(gaps) == 2 * 8 # 2 sizes x T periods + nat = DiagnosticReport(res).to_dict()["estimator_native_diagnostics"] + assert nat["regression_weights"]["status"] == "ran" + assert "n_extrapolating" in nat["regression_weights"] + assert nat["sparse_synthetic_control"]["status"] == "ran" + assert len(nat["sparse_synthetic_control"]["sizes"]) == 2 + + +def test_sparse_max_subsets_invalid_raises(): + df, _years, _T0 = _make_panel(n_donors=5, T=8, T0=4, seed=1) + res = _fit_iv(df) + for bad in (0, -1, 2.5, np.nan): + with pytest.raises(ValueError, match="max_subsets must be a positive integer"): + res.sparse_synthetic_control(sizes=[2], max_subsets=bad) + + +def test_sparse_empty_sizes_raises(): + df, _years, _T0 = _make_panel(n_donors=5, T=8, T0=4, seed=1) + res = _fit_iv(df) + with pytest.raises(ValueError, match="non-empty"): + res.sparse_synthetic_control(sizes=[]) + + +def test_cv_tail_diagnostics_use_validation_window_capture(): + # Exercise the special v_method="cv" snapshot capture: the fixed (X1s, X0s, V) triple is + # the VALIDATION-window standardized predictor matrices (re-aggregated per window), NOT the + # full-pre matrices. Both tail diagnostics must operate on that captured cv space. + from diff_diff.synthetic_control import _inner_solve_W + + df, _years, _T0 = _make_panel(n_donors=6, T=8, T0=6, seed=1) + res = _fit_cv(df) # v_method="cv", spanning special predictors (_CV_SPANNING, k=2) + assert res._fit_converged + snap = res._fit_snapshot + # Captured matrices are the 2 spanning specs x 6 donors, in the cv validation-window space. + assert snap.fit_X0s.shape == (2, 6) + assert snap.fit_X1s.shape == (2,) + assert snap.fit_v.shape == (2,) + + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + rw = res.regression_weights() + assert res._regw_status == "ran" + assert len(rw) == 6 + + # Sparse winner (size 2) matches an independent brute-force over the CAPTURED cv-space + # matrices with the fixed cv V — proving the cv path holds V fixed in the right space. + import itertools + + tab = res.sparse_synthetic_control(sizes=[2]) + won = tab[tab["status"] == "ran"].iloc[0]["donor_ids"] + Y = snap.pivots[snap.outcome] + Z1_pre = Y.loc[snap.pre_periods, snap.treated_id].to_numpy(float) + Z0_pre = Y.loc[snap.pre_periods, snap.donor_ids].to_numpy(float) + best, best_mspe = None, np.inf + for cols in itertools.combinations(range(len(snap.donor_ids)), 2): + w, conv = _inner_solve_W( + snap.fit_X1s, + snap.fit_X0s[:, list(cols)], + snap.fit_v, + snap.inner_max_iter, + snap.inner_min_decrease, + ) + if not conv: + continue + m = float(np.mean((Z1_pre - Z0_pre[:, list(cols)] @ w) ** 2)) + if m < best_mspe: + best_mspe, best = m, cols + assert won == tuple(snap.donor_ids[c] for c in best) + + +def test_sparse_non_integer_sizes_raise(): + df, _years, _T0 = _make_panel(n_donors=6, T=8, T0=4, seed=1) + res = _fit_iv(df) + # int(2.9) would silently truncate to size 2 -> reject non-integral / bool sizes. + for bad in ([2.9], [2, 3.0], True, [True], ["2"]): + with pytest.raises(ValueError, match="must be integer|int or a sequence"): + res.sparse_synthetic_control(sizes=bad) + # A valid numpy-int size still works. + tab = res.sparse_synthetic_control(sizes=[np.int64(2)]) + assert set(tab[tab["status"] == "ran"]["size"]) == {2}