From 1081d784da1e62491ec942ea9e4be9791ad63b21 Mon Sep 17 00:00:00 2001 From: igerber Date: Fri, 10 Jul 2026 07:58:09 -0400 Subject: [PATCH 1/2] feat(conformal): one-sided alternatives + covariates in the CWZ proxy MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Extends the SyntheticControl conformal layer (CWZ 2021) with both user-locked extensions: (a) alternative={"two-sided","greater","less"} on conformal_test / conformal_confidence_intervals / conformal_average_effect. One-sided tests use the SIGNED average-effect statistic S(u)=T_*^{-1/2} sum u_t, grounded in Remark 1's statistic freedom + eq 2 (the paper has NO section 7 - the prior TODO/REGISTRY "(§7 signed-t)" wording was a Firpo-Possebom cross-contamination, fixed here). q=1 enforced for one-sided. CI inversion reports half-line HULLS: the infinite side attaches because p -> 1 in the accepted direction; contiguous is RECOMPUTED against the final reported ray (Algorithm 1 has no monotonicity guarantee after per-candidate proxy refits); an all-rejected scanned grid reports the infinite side + an uncertified NaN finite endpoint as grid_limited (never a false empty set). (b) covariates=[...] stacks pivoted-variable matching rows RAW into the constrained-LS proxy objective per the paper's note after eq 6; residuals/statistic stay outcome-only; exchangeability preserved. EXOGENEITY ASSUMPTION warned on every covariate run: Z(theta0) null-imputes only the outcome, so covariates enter as observed (no-intervention) values - a policy-affected covariate leaks treated post-period information into the proxy. Block-collapsed average-effect covariates collapse with the same T*-block structure; pointwise CIs precompute the covariate panel once and slice per (pre, t) sub-series. DiagnosticReport's conformal block carries alternative + covariates for all three kinds. Scope excludes (c) the AR/innovation-permutation path (TODO row rewritten to track only that remainder). Tests (5 local review rounds): hand-rolled signed-statistic permutation oracle; directional rejection; half-line endpoints bracketing the two-sided interval; missed-ray grid_limited (9-block fixture - the single-post-slot tie structure floors p at 1/n_blocks); deterministic non-contiguous patterns incl. rejected-toward-infinity; proxy-weight mechanism; exogenous covariate fixture + assumption-warning assertion; report propagation; permutation floor with covariates. All pre-existing conformal + SC + diagnostic-report suites pass unchanged. Reference scoping: the authors' scinference R package anchors only the base machinery; extensions anchor to paper formulas via the oracle tests. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01Lbd6nqWmg4snvvBmegwqiw --- CHANGELOG.md | 13 + TODO.md | 2 +- diff_diff/conformal.py | 182 ++++++++++-- diff_diff/diagnostic_report.py | 5 + diff_diff/guides/llms-full.txt | 4 +- diff_diff/synthetic_control_results.py | 204 ++++++++++++- docs/methodology/REGISTRY.md | 12 +- tests/test_methodology_synthetic_control.py | 301 ++++++++++++++++++++ 8 files changed, 692 insertions(+), 31 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 407f59e0..40ead62d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -223,6 +223,19 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 assert between the two (defense in depth). ### Added +- **`SyntheticControl` conformal extensions: one-sided alternatives + covariates in the + proxy (Chernozhukov-Wüthrich-Zhu 2021).** `conformal_test` / `conformal_confidence_intervals` + / `conformal_average_effect` gain `alternative={"two-sided","greater","less"}` — one-sided + tests use the SIGNED average-effect statistic `S(û)=T_*^{-1/2}·Σ û_t` (grounded in Remark 1's + statistic freedom; **the paper has no §7** — this PR also fixes the prior "(§7 signed-t)" + citation, a cross-contamination from the Firpo-Possebom scope note), `q=1` enforced, and CI + inversion yields half-lines with the infinite side genuinely accepted. `covariates=[...]` + stacks pivoted-variable matching rows into the constrained-LS proxy objective per the paper's + note after eq 6; rows stack raw (pre-scale as needed), residuals/statistic stay outcome-only, + exchangeability preserved; block-collapsed average-effect covariates collapse with the same + T*-block structure. Locked by a hand-rolled signed-statistic permutation oracle, directional + rejection/half-line/composition tests, a proxy-weight mechanism check, and a permutation-floor + invariant; all 31 pre-existing conformal tests pass unchanged. - **Opt-in `df_convention="cluster"` inference-df knob (DiD / TWFE / MultiPeriodDiD + `LinearRegression`).** Clustered analytical fits historically compute t-statistics, p-values, and CIs at the fitted **residual df** (`n − K_full`), while `fixest`/Stata use diff --git a/TODO.md b/TODO.md index bf3ba813..831b9792 100644 --- a/TODO.md +++ b/TODO.md @@ -28,7 +28,7 @@ The `Origin` column (Actionable tables) and the `PR` column (Deferred tables) bo | Issue | Location | Origin | Effort | Priority | |-------|----------|--------|--------|----------| -| `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 | +| `SyntheticControl` conformal (CWZ 2021) AR / innovation-permutation path (Lemmas 5-7) for time-series proxies — the residual-permutation shortcut is only valid for time-permutation-invariant proxies (SC/Lasso/DiD); an AR proxy needs innovation permutation. One-sided alternatives (Remark 1 signed statistic) and proxy covariates (eq 4/6 note) SHIPPED 2026-07. | `diff_diff/conformal.py`, `diff_diff/synthetic_control_results.py` | CWZ-2021 | Heavy | Low | | `ContinuousDiD` CGBS-2024 extensions. (a) `covariates=` kwarg — **DONE (reg/dr)**; (b) discrete-treatment saturated regression (`treatment_type="discrete"`) — **DONE**; (c) lowest-dose-as-control per Remark 3.1 when `P(D=0)=0` (`control_group="lowest_dose"`) — **DONE** (discrete + continuous mass-point, single-cohort; estimand `ATT(d)−ATT(d_L)`; see REGISTRY Note #7). Remaining (all deferred `NotImplementedError`, documented): `estimation_method="ipw"` on the dose curve (scalar-adjustment / degenerate); `covariates=` × `survey_design=` (weighted OR + weighted nuisance IF); multi-cohort **heterogeneous-support** discrete aggregation (support-aware: average each dose only over the cohorts that observe it); **multi-cohort `lowest_dose`** (within-cohort `d_L` reference + support-aware cross-cohort aggregation); and **`covariates=` × `lowest_dose`** (conditional-PT-relative-to-`d_L` estimand). Single-cohort / 2-period / shared-support multi-cohort are supported. | `continuous_did.py` | CGBS-2024 | Heavy | Low | | TWFE's HC2/HC2-BM inline full-dummy build (`twfe.py:280-315`) duplicates the dummy-construction logic in `DifferenceInDifferences(fixed_effects=...)` (`estimators.py:478-486`). Extract a shared helper, or delegate TWFE's HC2/HC2-BM path to DiD's `fixed_effects=` branch (with TWFE-specific cluster-default threading), to reduce drift risk on FE naming / survey behavior / result-surface conventions. Substantive refactor — touches both estimators. | `twfe.py::fit`, `estimators.py::DifferenceInDifferences.fit` | follow-up | Heavy | Low | diff --git a/diff_diff/conformal.py b/diff_diff/conformal.py index 3cb4f76f..febe4c6b 100644 --- a/diff_diff/conformal.py +++ b/diff_diff/conformal.py @@ -21,8 +21,10 @@ **Proxy.** The proxy is the canonical CWZ constrained-LS synthetic control (eqs 3–4): simplex weights minimising ``Σ_{t}(Y^N_{1t} − Σ_j w_j Y^N_{jt})²`` over **all** periods under the null, ``w ≥ 0, Σ w = 1``, **no V-matrix, no -intercept, outcomes-only** (footnote 9: "we estimate w under the null based on -all the data"). This is DISTINCT from the headline ADH V-matrix weights — CWZ's +intercept, outcomes by default** (footnote 9: "we estimate w under the null +based on all the data"); optional RAW covariate-matching rows stack into the +same objective when ``covariates=`` is supplied (the paper's note after eq 6), +with outcome-only residuals. This is DISTINCT from the headline ADH V-matrix weights — CWZ's exactness theory (Lemma 1, Appendix D exchangeability) requires a time-symmetric proxy, which the ADH pre-period V-fit is not. Reuses the Frank-Wolfe simplex solver :func:`diff_diff.utils._sc_weight_fw`. @@ -53,6 +55,8 @@ def _cwz_proxy_fit( max_iter: int, min_decrease: float, init_weights: Optional[np.ndarray] = None, + x1_rows: Optional[np.ndarray] = None, + X0_rows: Optional[np.ndarray] = None, ) -> Tuple[np.ndarray, np.ndarray, bool]: """Fit the canonical CWZ constrained-LS SC proxy over ALL given periods (eq 4). @@ -67,16 +71,38 @@ def _cwz_proxy_fit( caller must pass a θ0-invariant scale (e.g. the pre-window outcome norm) so CI membership does not drift with the grid value via the tolerance. + **Covariates (CWZ §2.3.2, after eq 6):** "it is straightforward to + incorporate (transformations of) covariates X_jt into the estimation + problems (4) and (6)." ``x1_rows`` ``(R,)`` / ``X0_rows`` ``(R, J)`` are + OPTIONAL covariate-matching rows stacked under the outcome rows in the + same simplex LS objective — each row enters exactly like an outcome + period, so the stacked objective remains invariant under time + permutations of the outcome rows (exchangeability preserved; the + covariate rows are fixed features of every permuted dataset). The + RESIDUALS returned are OUTCOME rows only — the statistic and the + permutation p-value are unchanged in form. Covariates are stacked RAW + (no internal standardization): the paper's "(transformations of)" + delegates scaling to the caller, and a covariate on a large scale + dominates the fit exactly as it would in eq 4 — pre-scale accordingly. + Returns ``(w (J,), resid = y1 − Y0 @ w (T,), converged)``. ``J == 1`` is the degenerate single-donor case ``w = [1]`` (no optimisation). """ y1 = np.asarray(y1, dtype=float) Y0 = np.asarray(Y0, dtype=float) + n_out = y1.shape[0] + if x1_rows is not None: + x1_rows = np.asarray(x1_rows, dtype=float) + X0_rows = np.asarray(X0_rows, dtype=float) + y1_fit = np.concatenate([y1, x1_rows]) + Y0_fit = np.vstack([Y0, X0_rows]) + else: + y1_fit, Y0_fit = y1, Y0 _, J = Y0.shape if J == 1: w = np.array([1.0], dtype=float) return w, y1 - Y0 @ w, True - packed = np.column_stack([Y0, y1]) # (T, J+1); last column is the target + packed = np.column_stack([Y0_fit, y1_fit]) # (T+R, J+1); last column is the target with warnings.catch_warnings(): warnings.filterwarnings("ignore", message=r".*did not converge.*", category=UserWarning) w, converged = _sc_weight_fw( @@ -89,7 +115,7 @@ def _cwz_proxy_fit( return_convergence=True, ) w = np.asarray(w, dtype=float) - return w, y1 - Y0 @ w, bool(converged) + return w, y1[:n_out] - Y0[:n_out] @ w, bool(converged) # ============================================================================= @@ -160,7 +186,11 @@ def _make_perms(m: int, scheme: str, n_iid: int, rng: np.random.Generator) -> np def _cwz_pvalue( - u: np.ndarray, post_mask: np.ndarray, perms: np.ndarray, q: Any + u: np.ndarray, + post_mask: np.ndarray, + perms: np.ndarray, + q: Any, + alternative: str = "two-sided", ) -> Tuple[float, float, int]: """CWZ eq (2) permutation p-value for residuals ``u``. @@ -169,15 +199,34 @@ def _cwz_pvalue( identity ``π`` is in ``perms`` so ``p̂ ≥ 1/|Π|`` automatically (no extra ``+1``). Ties (``≥``) are counted conservatively. ``perms`` is the ``(|Π|, len(u))`` integer index array from :func:`_make_perms`. Returns ``(p, S_observed, |Π|)``. + + **One-sided alternatives (CWZ Remark 1).** The permutation argument is + statistic-agnostic ("other test statistics can be used as well"), so + ``alternative="greater"`` uses the SIGNED average-effect statistic + ``S(û) = T_*^{-1/2}·Σ_{t>T0} û_t`` — large positive residuals (treated + outcome above the counterfactual proxy, i.e. ``θ > θ0``) reject; + ``"less"`` negates the residuals and applies the same rule. NOTE the + paper has NO §7 and no named one-sided section — the variants ground in + Remark 1's statistic freedom + eq 2. ``q`` applies only to the two-sided + ``S_q`` family and must be 1 for one-sided alternatives (enforced by the + callers). """ - s_obs = _cwz_statistic(u, post_mask, q) - post_perm = np.abs(u[perms][:, post_mask]) # (|Π|, T*) - n_star = post_perm.shape[1] - if q == _INF: - s_perm = post_perm.max(axis=1) + if alternative == "less": + u = -u + if alternative in ("greater", "less"): + post = u[post_mask] + n_star = post.size + s_obs = float(np.sum(post) / np.sqrt(n_star)) + s_perm = u[perms][:, post_mask].sum(axis=1) / np.sqrt(n_star) else: - s_perm = (post_perm**q).sum(axis=1) / np.sqrt(n_star) - s_perm = s_perm ** (1.0 / q) + s_obs = _cwz_statistic(u, post_mask, q) + post_perm = np.abs(u[perms][:, post_mask]) # (|Π|, T*) + n_star = post_perm.shape[1] + if q == _INF: + s_perm = post_perm.max(axis=1) + else: + s_perm = (post_perm**q).sum(axis=1) / np.sqrt(n_star) + s_perm = s_perm ** (1.0 / q) n = perms.shape[0] tol = 1e-12 * max(abs(s_obs), 1.0) n_ge = int(np.sum(s_perm >= s_obs - tol)) @@ -200,6 +249,9 @@ def _single_null_pvalue( max_iter: int, min_decrease: float, init_weights: Optional[np.ndarray] = None, + alternative: str = "two-sided", + x1_rows: Optional[np.ndarray] = None, + X0_rows: Optional[np.ndarray] = None, ) -> Dict[str, Any]: """Build ``Z(θ0)``, fit the proxy under the null, and return the p-value. @@ -211,9 +263,15 @@ def _single_null_pvalue( y1n = np.asarray(y1_obs, dtype=float).copy() y1n[post_mask] = y1n[post_mask] - np.asarray(effect_on_post, dtype=float) w, resid, conv = _cwz_proxy_fit( - y1n, Y0, max_iter=max_iter, min_decrease=min_decrease, init_weights=init_weights + y1n, + Y0, + max_iter=max_iter, + min_decrease=min_decrease, + init_weights=init_weights, + x1_rows=x1_rows, + X0_rows=X0_rows, ) - p, s_obs, n = _cwz_pvalue(resid, post_mask, perms, q) + p, s_obs, n = _cwz_pvalue(resid, post_mask, perms, q, alternative=alternative) return {"p_value": p, "s_observed": s_obs, "n_perms": n, "converged": conv, "weights": w} @@ -292,6 +350,9 @@ def _invert_single_post( min_decrease: float, grid: Optional[np.ndarray] = None, n_grid: int = 100, + alternative: str = "two-sided", + x1_rows: Optional[np.ndarray] = None, + X0_rows: Optional[np.ndarray] = None, ) -> Dict[str, Any]: """Invert a single-post-slot conformal test over a grid (CWZ Algorithm 1). @@ -301,8 +362,13 @@ def _invert_single_post( (``y1[post_idx] − θ``), the proxy refit (warm-started across the grid — FW is convex so warm-starting only affects speed, not the optimum), and the permutation p-value recomputed. With one post slot ``S_q`` reduces to - ``|û_post|`` for every ``q``, so ``q`` is inert here (fixed to 1). Returns the CI - summary plus the ``grid`` table of ``(θ, p, in_set, converged)`` rows. + ``|û_post|`` for every ``q``, so ``q`` is inert here (fixed to 1); under + ``alternative="greater"``/``"less"`` the single-slot statistic is the SIGNED + ``û_post`` (resp. ``−û_post``), and the inversion yields a ONE-SIDED + confidence set — a half-line whose finite endpoint is the reported bound + (the infinite side is genuinely accepted, not grid-limited; the caller + overrides the edge-touch status accordingly). Returns the CI summary plus + the ``grid`` table of ``(θ, p, in_set, converged)`` rows. When ``alpha < 1/|Π|`` every candidate has ``p ≥ 1/|Π| > alpha`` (the identity is in ``Π``), so NO value is ever rejected and the confidence set is the whole line: @@ -323,7 +389,12 @@ def _invert_single_post( # relying on that naive residual). ``spread`` is the pre-fit residual noise scale. pre_sel = ~post_mask w_pre, resid_pre, _ = _cwz_proxy_fit( - y1_obs[pre_sel], Y0[pre_sel], max_iter=max_iter, min_decrease=min_decrease + y1_obs[pre_sel], + Y0[pre_sel], + max_iter=max_iter, + min_decrease=min_decrease, + x1_rows=x1_rows, + X0_rows=X0_rows, ) point_est = float(y1_obs[post_idx] - Y0[post_idx] @ w_pre) spread = float(np.std(resid_pre)) if resid_pre.size else 1.0 @@ -354,10 +425,16 @@ def _invert_single_post( y1n = y1_obs.copy() y1n[post_idx] = y1n[post_idx] - float(theta) w, resid, conv = _cwz_proxy_fit( - y1n, Y0, max_iter=max_iter, min_decrease=min_decrease, init_weights=w_prev + y1n, + Y0, + max_iter=max_iter, + min_decrease=min_decrease, + init_weights=w_prev, + x1_rows=x1_rows, + X0_rows=X0_rows, ) w_prev = w - p, _, _ = _cwz_pvalue(resid, post_mask, perms, 1) + p, _, _ = _cwz_pvalue(resid, post_mask, perms, 1, alternative=alternative) pvals[i] = p converged[i] = conv @@ -368,6 +445,73 @@ def _invert_single_post( return out +def _apply_one_sided_endpoints(res: Dict[str, Any], alternative: str) -> Dict[str, Any]: + """Convert a one-sided inversion's accepted half-line into ±inf endpoints. + + Under ``alternative="greater"`` (H1: θ > θ0) small candidates are rejected + and the accepted region extends to the grid's upper edge, so the upper + endpoint becomes ``+inf`` and only a LOWER-edge touch keeps + ``status="grid_limited"``. ``"less"`` is the mirror image. No-op for + ``"two-sided"`` or for unbounded results (everything accepted). + + **Hull convention (no monotonicity guarantee).** Algorithm 1 accepts any + candidate with ``p > α`` after REFITTING the proxy at that candidate; a + monotone p-curve is typical but NOT guaranteed, so an + accepted/rejected/accepted pattern is possible. Exactly as in the + two-sided convention, the reported endpoints are then the HULL of the + accepted set (with the infinite side attached) and ``contiguous`` is + RECOMPUTED against the final reported ray — any rejected scanned point + strictly inside it (including between the finite hull edge and the + attached infinity) flips ``contiguous=False``, so a rejected pocket is + disclosed, not hidden (callers warn on it). + An all-rejected grid (the scan missed the ray — narrow user ``bounds``) + reports the infinite side + an UNCERTIFIED (NaN) finite endpoint with + ``grid_limited`` status: the ray always exists beyond the grid on the + accepted side (``p → 1`` as ``θ0`` walks into the accepted direction). + """ + if alternative == "two-sided" or res["status"] == "unbounded": + return res + if res["status"] == "empty": + # Every scanned grid point rejected. Under a one-sided alternative + # the accepted ray ALWAYS exists beyond the grid in the accepted + # direction (as θ0 → ±∞ on that side, the signed statistic walks + # into the permutation distribution's far tail and p → 1), so an + # all-rejected grid means the scan missed it — typically narrow + # user-supplied bounds. Report the infinite side with the finite + # endpoint UNCERTIFIED (NaN) and grid_limited status rather than a + # (false) empty confidence set. + if alternative == "greater": + res["upper"] = float("inf") + else: + res["lower"] = -float("inf") + res["status"] = "grid_limited" + return res + grid = res.get("grid") or [] + if not grid: + return res + in_first = bool(grid[0][2]) + in_last = bool(grid[-1][2]) + if alternative == "greater": + res["upper"] = float("inf") + res["status"] = "grid_limited" if in_first else "ran" + # Recompute contiguity against the FINAL reported ray: the + # underlying flag only covered the finite accepted hull, so a + # rejected scanned point ABOVE the hull (inside [lower, +inf) after + # the rewrite) must flip it (REGISTRY hull-disclosure contract). + lower = res["lower"] + res["contiguous"] = bool(res["contiguous"]) and all( + bool(row[2]) for row in grid if row[0] > lower + ) + else: # "less" + res["lower"] = -float("inf") + res["status"] = "grid_limited" if in_last else "ran" + upper = res["upper"] + res["contiguous"] = bool(res["contiguous"]) and all( + bool(row[2]) for row in grid if row[0] < upper + ) + return res + + # ============================================================================= # Block collapse for the average-effect test (CWZ Appendix A.1) # ============================================================================= diff --git a/diff_diff/diagnostic_report.py b/diff_diff/diagnostic_report.py index b91e7936..7b8cf2f6 100644 --- a/diff_diff/diagnostic_report.py +++ b/diff_diff/diagnostic_report.py @@ -2690,6 +2690,11 @@ def _scm_native(self, r: Any) -> Dict[str, Any]: "scheme": ci_summary.get("scheme"), "n_perms": _to_python_scalar(ci_summary.get("n_perms")), "n_post": _to_python_scalar(ci_summary.get("n_post")), + # One-sided / covariate-proxy runs are uninterpretable without + # these two — a one-sided p-value or half-line CI reads wrong + # under the default two-sided assumption. + "alternative": ci_summary.get("alternative", "two-sided"), + "covariates": ci_summary.get("covariates"), } if kind == "joint": block["q"] = ci_summary.get("q") diff --git a/diff_diff/guides/llms-full.txt b/diff_diff/guides/llms-full.txt index 25e037a5..3039d113 100644 --- a/diff_diff/guides/llms-full.txt +++ b/diff_diff/guides/llms-full.txt @@ -621,7 +621,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); `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). +**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). All three accept `alternative={"two-sided","greater","less"}` (one-sided uses the SIGNED average-effect statistic per CWZ Remark 1, q fixed at 1; CIs become half-lines `[lower,+inf)`/`(-inf,upper]`) and `covariates=[...]` (pivoted-variable matching rows stacked RAW into the constrained-LS proxy per the note after eq 6; residuals/p-value stay outcome-only). `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:** @@ -1400,7 +1400,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()`, `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()` +**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, alternative="two-sided", covariates=None, 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,∞}`; one-sided `alternative` = signed statistic per Remark 1 with q=1; `covariates=` stacks pivoted-variable matching rows into the proxy), `conformal_confidence_intervals(alpha=0.1, alternative="two-sided", covariates=None, scheme="moving_block", bounds=None, n_grid=100, seed=None)` (pointwise per-period CIs, Algorithm 1 — each period uses Z=(pre, t); one-sided → half-lines), `conformal_average_effect(alpha=0.1, alternative="two-sided", covariates=None, scheme="moving_block", bounds=None, n_grid=200, seed=None)` (average-effect CI by T*-block collapse, Appendix A.1; one-sided → half-line; covariate rows block-collapse identically) + `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_results.py b/diff_diff/synthetic_control_results.py index 18f53129..f728b5f2 100644 --- a/diff_diff/synthetic_control_results.py +++ b/diff_diff/synthetic_control_results.py @@ -2607,11 +2607,101 @@ def _conformal_panel(self) -> Tuple[np.ndarray, np.ndarray, int, int, float, Lis pre_scale = max(float(np.linalg.norm(y1[:n_pre])), 1e-12) return y1, Y0, n_pre, n_post, pre_scale, donors + @staticmethod + def _coerce_alternative(alternative: str, q: Any = 1) -> str: + """Validate ``alternative`` and its interaction with ``q``. + + One-sided alternatives use the SIGNED average-effect statistic (CWZ + Remark 1 — the permutation framework is statistic-agnostic; the paper + has no dedicated one-sided section), for which the ``S_q`` norm order + does not apply: ``q`` must be left at its default 1. + """ + if alternative not in ("two-sided", "greater", "less"): + raise ValueError( + "alternative must be 'two-sided', 'greater', or 'less', " f"got {alternative!r}" + ) + if alternative != "two-sided" and q != 1: + raise ValueError( + "q applies only to the two-sided S_q statistic; one-sided " + "alternatives use the signed average-effect statistic " + f"(CWZ Remark 1). Got q={q!r} with alternative={alternative!r}." + ) + return alternative + + def _conformal_covariate_rows( + self, covariates: Optional[List[str]], cal_periods: List[Any] + ) -> Tuple[Optional[np.ndarray], Optional[np.ndarray]]: + """Build stacked covariate-matching rows for the CWZ proxy (eq 4 note). + + ``covariates`` names variables retained in the fit snapshot's pivots + (any variable the estimator pivoted — predictors or additional + columns). For each variable the treated column and donor columns over + ``cal_periods`` are stacked as extra matching rows. Rows are stacked + RAW — the paper's "(transformations of) covariates" delegates scaling + to the caller, so pre-scale variables of very different magnitudes. + Returns ``(None, None)`` when ``covariates`` is None/empty. + """ + if not covariates: + return None, None + # CWZ builds Z(θ0) by null-imputing ONLY the outcome; covariates + # enter as their NO-INTERVENTION values. Pulling observed pivot + # values is therefore only valid for covariates unaffected by the + # policy — a policy-affected covariate leaks treated post-period + # intervention information into the proxy weights and the p-value + # silently conditions on it. Warn on every covariate run so the + # assumption is on the record (no-silent-failures contract). + warnings.warn( + "conformal covariate matching assumes the covariates are " + "UNAFFECTED by the intervention (their observed post-period " + "values are used as no-intervention values in Z(theta0); only " + "the outcome is null-imputed). A policy-affected covariate " + "leaks treatment information into the proxy. Pass only " + "exogenous covariates.", + UserWarning, + stacklevel=3, + ) + snap = self._fit_snapshot + assert snap is not None + donors = list(snap.donor_ids) + x1_parts: List[np.ndarray] = [] + X0_parts: List[np.ndarray] = [] + for var in covariates: + if var == snap.outcome: + raise ValueError( + f"covariates must not include the outcome ({var!r}); the " + "outcome rows are always in the proxy objective." + ) + if var not in snap.pivots: + raise ValueError( + f"covariate {var!r} is not in the fit snapshot's pivoted " + f"variables {sorted(snap.pivots.keys())!r}. Pass a variable " + "the SyntheticControl fit pivoted (e.g. a predictor)." + ) + piv = snap.pivots[var] + missing = [c for c in cal_periods if c not in piv.index] + if missing: + raise ValueError( + f"covariate {var!r} is missing periods {missing[:5]!r} in the " "fitted panel." + ) + x1 = piv.loc[cal_periods, snap.treated_id].to_numpy(dtype=float) + X0 = piv.loc[cal_periods, donors].to_numpy(dtype=float) + if not (np.all(np.isfinite(x1)) and np.all(np.isfinite(X0))): + raise ValueError( + f"covariate {var!r} has non-finite cells over the requested " + "periods; conformal covariate matching requires a complete " + "finite panel." + ) + x1_parts.append(x1) + X0_parts.append(X0) + return np.concatenate(x1_parts), np.vstack(X0_parts) + def conformal_test( self, effect: Any, *, q: Any = 1, + alternative: str = "two-sided", + covariates: Optional[List[str]] = None, scheme: str = "moving_block", n_iid: int = 10000, seed: Optional[int] = None, @@ -2639,7 +2729,33 @@ def conformal_test( q : {1, 2, inf}, default 1 The ``S_q`` norm order. ``1`` (robust to heavy tails — the paper's application default), ``2`` (permanent effects), ``inf`` (= ``max|û_t|``, - large temporary effects). + large temporary effects). Two-sided only; one-sided alternatives + fix ``q=1`` (signed statistic). + alternative : {"two-sided", "greater", "less"}, default "two-sided" + ``"two-sided"`` uses the paper's ``S_q`` magnitude statistic. + One-sided alternatives use the SIGNED average-effect statistic + ``S(û) = T_*^{-1/2}·Σ_{t>T0} û_t`` (CWZ Remark 1 — the + permutation framework is valid for any statistic; the paper has + no dedicated one-sided section): ``"greater"`` rejects when the + treated outcomes sit ABOVE the counterfactual proxy (``θ > θ0``), + ``"less"`` mirrors it. Requires ``q=1``. On the CI surfaces the + inversion yields half-lines: ``[lower, +inf)`` for + ``"greater"``, ``(-inf, upper]`` for ``"less"``. + covariates : list of str, optional + Variable names (pivoted by the fit — e.g. ``predictors``) whose + treated/donor series are stacked as covariate-matching rows in + the CWZ proxy objective (the paper after eq 6: "straightforward + to incorporate (transformations of) covariates X_jt into the + estimation problems (4) and (6)"). Rows are stacked RAW — no + internal standardization; pre-scale covariates of very different + magnitudes. Residuals and the permutation p-value remain + outcome-only, so exactness is preserved (covariate rows are + fixed features of every permuted dataset). **Exogeneity + assumption (warned on every run):** ``Z(θ0)`` null-imputes ONLY + the outcome — covariate values enter as observed, i.e. as + no-intervention paths, so covariates must be unaffected by the + policy; a policy-affected covariate leaks treated post-period + information into the proxy. scheme : {"moving_block", "iid"}, default "moving_block" The permutation set. ``"moving_block"`` (``Π_→``, ``T`` cyclic shifts) is valid under serially-dependent / stationary weakly-dependent errors @@ -2668,11 +2784,16 @@ def conformal_test( from diff_diff.conformal import _INF, _make_perms, _single_null_pvalue q = self._coerce_q(q) + alternative = self._coerce_alternative(alternative, q) if scheme not in ("moving_block", "iid"): raise ValueError(f"scheme must be 'moving_block' or 'iid', got {scheme!r}") if not isinstance(n_iid, (int, np.integer)) or n_iid < 1: raise ValueError(f"n_iid must be a positive integer, got {n_iid!r}") y1, Y0, n_pre, n_post, pre_scale, _ = self._conformal_panel() + snap0 = self._fit_snapshot + assert snap0 is not None + cal = list(snap0.pre_periods) + list(snap0.post_periods) + x1_rows, X0_rows = self._conformal_covariate_rows(covariates, cal) f_post = self._coerce_effect_path(effect, n_post) if n_post >= n_pre: warnings.warn( @@ -2698,6 +2819,9 @@ def conformal_test( q, max_iter=snap.inner_max_iter, min_decrease=snap.inner_min_decrease * pre_scale, + alternative=alternative, + x1_rows=x1_rows, + X0_rows=X0_rows, ) if not res["converged"]: warnings.warn( @@ -2712,6 +2836,8 @@ def conformal_test( "kind": "joint", "scheme": scheme, "q": q_out, + "alternative": alternative, + "covariates": list(covariates) if covariates else None, "alpha": None, "n_perms": int(res["n_perms"]), "n_post": int(n_post), @@ -2724,6 +2850,7 @@ def conformal_test( "p_value": float(res["p_value"]), "S_observed": float(res["s_observed"]), "q": q_out, + "alternative": alternative, "scheme": scheme, "n_perms": int(res["n_perms"]), "n_post": int(n_post), @@ -2735,6 +2862,8 @@ def conformal_average_effect( self, *, alpha: float = 0.1, + alternative: str = "two-sided", + covariates: Optional[List[str]] = None, scheme: str = "moving_block", n_iid: int = 10000, bounds: Optional[Tuple[float, float]] = None, @@ -2760,6 +2889,14 @@ def conformal_average_effect( ---------- alpha : float, default 0.1 The confidence level is ``1 − alpha``; membership is ``p^θ̄0 > alpha``. + alternative : {"two-sided", "greater", "less"}, default "two-sided" + One-sided alternatives use the signed statistic (CWZ Remark 1) + and invert to a HALF-LINE: ``[lower, +inf)`` for ``"greater"``, + ``(-inf, upper]`` for ``"less"`` (the infinite side is genuinely + accepted, not grid-limited). See :meth:`conformal_test`. + covariates : list of str, optional + Covariate-matching rows stacked into the CWZ proxy objective + (raw — pre-scale as needed). See :meth:`conformal_test`. scheme : {"moving_block", "iid"}, default "moving_block" Permutation set over the collapsed blocks. n_iid : int, default 10000 @@ -2787,7 +2924,12 @@ def conformal_average_effect( If ``alpha`` / ``scheme`` / ``n_iid`` / ``n_grid`` / ``bounds`` are invalid, ``T0 < T*`` (no full pre-block), or the fit snapshot is unavailable. """ - from diff_diff.conformal import _block_collapse, _invert_single_post, _make_perms + from diff_diff.conformal import ( + _apply_one_sided_endpoints, + _block_collapse, + _invert_single_post, + _make_perms, + ) if scheme not in ("moving_block", "iid"): raise ValueError(f"scheme must be 'moving_block' or 'iid', got {scheme!r}") @@ -2797,6 +2939,7 @@ def conformal_average_effect( raise ValueError(f"alpha must be in (0, 1), got {alpha!r}") if not isinstance(n_grid, (int, np.integer)) or n_grid < 2: raise ValueError(f"n_grid must be an integer >= 2, got {n_grid!r}") + alternative = self._coerce_alternative(alternative) grid = _validate_conformal_bounds(bounds, int(n_grid)) y1, Y0, n_pre, n_post, _, _ = self._conformal_panel() if n_pre < n_post: @@ -2813,6 +2956,27 @@ def conformal_average_effect( stacklevel=2, ) y1b, Y0b, n_dropped = _block_collapse(y1, Y0, n_pre, n_post) + x1b_rows: Optional[np.ndarray] = None + X0b_rows: Optional[np.ndarray] = None + if covariates: + snap_cov = self._fit_snapshot + assert snap_cov is not None + cal_cov = list(snap_cov.pre_periods) + list(snap_cov.post_periods) + x1_rows, X0_rows = self._conformal_covariate_rows(covariates, cal_cov) + # Covariate rows collapse with the SAME T*-block structure so the + # collapsed panel remains a coherent Z (each block-averaged + # covariate row enters the proxy like a block-averaged outcome). + n_vars = x1_rows.shape[0] // len(cal_cov) + xb_parts, Xb_parts = [], [] + T_cal = len(cal_cov) + for v in range(n_vars): + xv = x1_rows[v * T_cal : (v + 1) * T_cal] + Xv = X0_rows[v * T_cal : (v + 1) * T_cal] + xvb, Xvb, _ = _block_collapse(xv, Xv, n_pre, n_post) + xb_parts.append(xvb) + Xb_parts.append(Xvb) + x1b_rows = np.concatenate(xb_parts) + X0b_rows = np.vstack(Xb_parts) if n_dropped: warnings.warn( f"conformal_average_effect: T0={n_pre} is not a multiple of T*={n_post}; " @@ -2847,12 +3011,18 @@ def conformal_average_effect( min_decrease=snap.inner_min_decrease * block_scale, grid=grid, n_grid=int(n_grid), + alternative=alternative, + x1_rows=x1b_rows, + X0_rows=X0b_rows, ) + res = _apply_one_sided_endpoints(res, alternative) _warn_conformal_ci_status(res, "conformal_average_effect") self.conformal_inference = { "kind": "average", "scheme": scheme, "alpha": float(alpha), + "alternative": alternative, + "covariates": list(covariates) if covariates else None, "n_perms": n_perms, "n_post": int(n_post), "n_blocks": n_blocks, @@ -2907,6 +3077,8 @@ def conformal_confidence_intervals( self, *, alpha: float = 0.1, + alternative: str = "two-sided", + covariates: Optional[List[str]] = None, scheme: str = "moving_block", n_iid: int = 10000, bounds: Optional[Tuple[float, float]] = None, @@ -2930,6 +3102,14 @@ def conformal_confidence_intervals( ---------- alpha : float, default 0.1 The confidence level is ``1 − alpha``; membership is ``p^c > alpha``. + alternative : {"two-sided", "greater", "less"}, default "two-sided" + One-sided alternatives use the signed statistic (CWZ Remark 1) + and invert to a HALF-LINE: ``[lower, +inf)`` for ``"greater"``, + ``(-inf, upper]`` for ``"less"`` (the infinite side is genuinely + accepted, not grid-limited). See :meth:`conformal_test`. + covariates : list of str, optional + Covariate-matching rows stacked into the CWZ proxy objective + (raw — pre-scale as needed). See :meth:`conformal_test`. scheme : {"moving_block", "iid"}, default "moving_block" Permutation set over the ``(T0+1)``-length sub-series. n_iid : int, default 10000 @@ -2957,7 +3137,11 @@ def conformal_confidence_intervals( If ``alpha`` / ``scheme`` / ``n_iid`` / ``n_grid`` / ``bounds`` are invalid or the fit snapshot is unavailable. """ - from diff_diff.conformal import _invert_single_post, _make_perms + from diff_diff.conformal import ( + _apply_one_sided_endpoints, + _invert_single_post, + _make_perms, + ) if scheme not in ("moving_block", "iid"): raise ValueError(f"scheme must be 'moving_block' or 'iid', got {scheme!r}") @@ -2967,6 +3151,7 @@ def conformal_confidence_intervals( raise ValueError(f"alpha must be in (0, 1), got {alpha!r}") if not isinstance(n_grid, (int, np.integer)) or n_grid < 2: raise ValueError(f"n_grid must be an integer >= 2, got {n_grid!r}") + alternative = self._coerce_alternative(alternative) grid_template = _validate_conformal_bounds(bounds, int(n_grid)) y1, Y0, n_pre, n_post, pre_scale, _ = self._conformal_panel() if n_pre <= 1: @@ -2997,8 +3182,15 @@ def conformal_confidence_intervals( grid_rows: List[Dict[str, Any]] = [] statuses: List[str] = [] any_noncontig = False + snap_cov = self._fit_snapshot + assert snap_cov is not None + pre_list = list(snap_cov.pre_periods) + post_list = list(snap_cov.post_periods) for k, period in enumerate(post_periods): sub_idx = list(range(n_pre)) + [n_pre + k] + # Covariate rows subset to the SAME sub-series periods (Z for the + # pointwise test is (Z_1..Z_T0, Z_t) — covariates ride along). + x1_rows, X0_rows = self._conformal_covariate_rows(covariates, pre_list + [post_list[k]]) res = _invert_single_post( y1[sub_idx], Y0[sub_idx], @@ -3009,7 +3201,11 @@ def conformal_confidence_intervals( min_decrease=md, grid=grid_template, n_grid=int(n_grid), + alternative=alternative, + x1_rows=x1_rows, + X0_rows=X0_rows, ) + res = _apply_one_sided_endpoints(res, alternative) ci_rows.append( { "period": period, @@ -3066,6 +3262,8 @@ def conformal_confidence_intervals( self.conformal_inference = { "kind": "pointwise", "scheme": scheme, + "alternative": alternative, + "covariates": list(covariates) if covariates else None, "alpha": float(alpha), "n_perms": n_perms, "n_post": int(n_post), diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index e6a77e7c..ddfc0f28 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -2435,12 +2435,12 @@ Classic synthetic control (donor/unit weights only) for a single treated unit, d - **`confidence_set(family="constant"|"linear", gamma=0.1, bounds=None, n_grid=200)`** inverts that test over a one-parameter family: `"constant"` → `f(t)=c` (Eq 15), a confidence **interval** for a constant effect (Eq 16); `"linear"` → `f(t)=c̃·(t−T0)` with the 1-based post-period index (Eq 17), a confidence **set** over the slope `c̃` (Eq 18). Membership is the paper's **strict** `p^c > γ` (Eq 14 — see the boundary Note). With `bounds=None` the set is recovered **EXACTLY**: `p^c` is a piecewise-constant step function (each placebo's indicator flips only at the real roots of `A_j(c)·D_1 = A_1(c)·D_j`, a quadratic in `c`), so the placebo breakpoints partition the line and `p` is evaluated once per induced open interval AND at each breakpoint — where a tie under `≥` can lift `p` above γ, so an isolated accepted singleton (a tangent / co-located root) is captured. The accepted set is the union of accepted intervals/points, with **no centering or monotonicity assumption** (a poor-pre-fit treated unit can have its accepted region in the tails, not around the point estimate). `bounds=(lo,hi)` instead scans a fixed grid (grid-limited membership). The summary is on `effect_confidence_set` (`status ∈ {"ran","empty","unbounded"}`) and the full grid on `get_confidence_set_df()`. **Fail-closed:** `γ < 1/(J+1)` ⇒ `p^c > γ` for every `c` ⇒ `"unbounded"` (`±inf` endpoints + warning — the discrete-granularity point, fn 8); a treated unit lacking the best pre-fit can give a one-sided unbounded edge; if no interval or breakpoint is accepted the set is `"empty"` (NaN endpoints); a non-contiguous accepted region (disjoint components / an isolated singleton) reports the `[lower, upper]` hull with `contiguous=False` + a warning; and `< 2` donors / a non-converged treated fit / an unpickled result (no placebo reference set) raise `ValueError`. **Scope:** sensitivity weights (`φ≠0`, Eqs 7–9), the general test-statistic menu `θ¹–θ⁵` (Eq 19), one-sided (§7's signed-`t` statistic), and the multiple-outcome / multiple-treated extensions (§6) are **deferred** (flagged in the paper review checklist). **Conformal inference (Chernozhukov-Wüthrich-Zhu 2021, opt-in):** Valid p-values for a hypothesized post-period effect trajectory and pointwise confidence intervals — what the placebo / test-inversion paths cannot give. Surfaced under `estimator_native_diagnostics`; the analytical `se`/`t_stat`/`p_value`/`conf_int` stay NaN (separate permutation object — see the non-analytical Note). Under a sharp null `H0: θ = θ0` the counterfactual treated outcomes are imputed (`Y^N_{1t} = Y_{1t} − θ0_t`, `t>T0`; pre unchanged), a **time-permutation-invariant proxy** is fit UNDER THE NULL on **all** periods, and the statistic `S_q(û) = ((1/√T*)·Σ_{t>T0}|û_t|^q)^{1/q}` (CWZ §2.2) is referred to its permutation distribution (eq 2) by reshuffling residuals over time: -- **Proxy (CWZ §2.3, eqs 3–4):** the canonical constrained-LS synthetic control — simplex weights `argmin_w Σ_{t=1}^T (Y^N_{1t} − Σ_j w_j Y^N_{jt})²` s.t. `w≥0, Σw=1`, **NO V-matrix, no intercept, outcomes-only**, fit over ALL periods (footnote 9: "we estimate w under the null based on all the data"). Reuses the Frank-Wolfe simplex solver `utils._sc_weight_fw` (packed `(T,J+1)`, `zeta=0, intercept=False`). -- **`conformal_test(effect, q=1, scheme="moving_block", n_iid=10000, seed=None)`** — joint sharp-null test (eqs 1–2). `effect` is a scalar (constant) or length-`n_post` path; `q ∈ {1,2,∞}` (`q=∞` is `max|û_t|`). The proxy is fit ONCE under the null and only residuals are permuted (footnote 7 — permuting residuals ≡ permuting data for a time-permutation-invariant proxy). Returns `p_value`, `S_observed`, `n_perms`, `proxy_converged`. -- **`conformal_confidence_intervals(alpha=0.1, scheme="moving_block", n_iid=10000, bounds=None, n_grid=100, seed=None)`** — pointwise per-period CIs by test inversion (Algorithm 1). Per the paper (§2.2), each per-period CI for period `t` uses `Z = (Z_1,…,Z_{T0}, Z_t)` — the pre-periods PLUS only period `t`, the **other post-periods dropped** — a clean `T*=1` conformal test (`q` is therefore inert: `S_q = |û_t|`). Returns one row per post period (`lower`/`upper`/`point_estimate`/`status`/`contiguous`); the full grid is on `get_conformal_grid_df()`. -- **`conformal_average_effect(alpha=0.1, scheme="moving_block", n_iid=10000, bounds=None, n_grid=200, seed=None)`** — CI for the average effect `T*^{-1}Σ_{t>T0}θ_t` by test inversion (Appendix A.1): **collapse** the panel into non-overlapping `T*`-blocks (per-unit block averages), fit the proxy on the collapsed `Z̄`, permute the `T/T*` **block** residuals. The earliest `T0 mod T*` pre-periods are dropped to make the block count integral. Requires `T0 ≥ T*`. +- **Proxy (CWZ §2.3, eqs 3–4):** the canonical constrained-LS synthetic control (outcomes by default; optional RAW covariate-matching rows stack into the same objective when `covariates=` is supplied — see the scope note) — simplex weights `argmin_w Σ_{t=1}^T (Y^N_{1t} − Σ_j w_j Y^N_{jt})²` s.t. `w≥0, Σw=1`, **NO V-matrix, no intercept, outcomes-only**, fit over ALL periods (footnote 9: "we estimate w under the null based on all the data"). Reuses the Frank-Wolfe simplex solver `utils._sc_weight_fw` (packed `(T,J+1)`, `zeta=0, intercept=False`). +- **`conformal_test(effect, q=1, alternative="two-sided", covariates=None, scheme="moving_block", n_iid=10000, seed=None)`** — joint sharp-null test (eqs 1–2). `effect` is a scalar (constant) or length-`n_post` path; `q ∈ {1,2,∞}` (`q=∞` is `max|û_t|`). The proxy is fit ONCE under the null and only residuals are permuted (footnote 7 — permuting residuals ≡ permuting data for a time-permutation-invariant proxy). Returns `p_value`, `S_observed`, `n_perms`, `proxy_converged`. +- **`conformal_confidence_intervals(alpha=0.1, alternative="two-sided", covariates=None, scheme="moving_block", n_iid=10000, bounds=None, n_grid=100, seed=None)`** — pointwise per-period CIs by test inversion (Algorithm 1). Per the paper (§2.2), each per-period CI for period `t` uses `Z = (Z_1,…,Z_{T0}, Z_t)` — the pre-periods PLUS only period `t`, the **other post-periods dropped** — a clean `T*=1` conformal test (`q` is therefore inert: `S_q = |û_t|`). Returns one row per post period (`lower`/`upper`/`point_estimate`/`status`/`contiguous`); the full grid is on `get_conformal_grid_df()`. +- **`conformal_average_effect(alpha=0.1, alternative="two-sided", covariates=None, scheme="moving_block", n_iid=10000, bounds=None, n_grid=200, seed=None)`** — CI for the average effect `T*^{-1}Σ_{t>T0}θ_t` by test inversion (Appendix A.1): **collapse** the panel into non-overlapping `T*`-blocks (per-unit block averages), fit the proxy on the collapsed `Z̄`, permute the `T/T*` **block** residuals. The earliest `T0 mod T*` pre-periods are dropped to make the block count integral. Requires `T0 ≥ T*`. - **Permutation schemes:** `"moving_block"` (`Π_→`, `m` cyclic shifts — valid under serially-dependent / stationary weakly-dependent errors, Assumption 2.2; the robust default) and `"iid"` (`Π_all`, sampled — valid under i.i.d. errors, Assumption 2.1, finer p-values). Both include the identity, so the permutation p-value is `≥ 1/|Π|` by construction. -- **Fail-closed:** `<1` donor or an unpickled result (no fit snapshot) or non-finite panel cells / `<2` periods → `ValueError`; a single donor warns (degenerate proxy `w=[1]`); `T*≥T0` warns (the validity caveat — large `T0` drives exactness); a non-converged grid point is **indeterminate**, not rejected — only a converged `p≤α` rejects, so the set is the complement of the rejections (excluding non-converged points would understate the CI width); an accepted set touching a grid edge is `status="grid_limited"`; an empty accepted set is `status="empty"`; `α < 1/|Π|` ⇒ every value accepted, short-circuited to `status="unbounded"` with `(−inf, +inf)` endpoints (warned). **Scope:** one-sided / signed variants (§7), covariates folded into the proxy, and the AR/innovation-permutation path (Lemmas 5–7, for time-series proxies) are **deferred** (flagged in the paper review checklist). +- **Fail-closed:** `<1` donor or an unpickled result (no fit snapshot) or non-finite panel cells / `<2` periods → `ValueError`; a single donor warns (degenerate proxy `w=[1]`); `T*≥T0` warns (the validity caveat — large `T0` drives exactness); a non-converged grid point is **indeterminate**, not rejected — only a converged `p≤α` rejects, so the set is the complement of the rejections (excluding non-converged points would understate the CI width); an accepted set touching a grid edge is `status="grid_limited"`; an empty accepted set is `status="empty"`; `α < 1/|Π|` ⇒ every value accepted, short-circuited to `status="unbounded"` with `(−inf, +inf)` endpoints (warned). **Scope:** one-sided alternatives and proxy covariates are **implemented (2026-07)** — `alternative={"two-sided","greater","less"}` on all three conformal methods uses the SIGNED average-effect statistic `S(û)=T_*^{-1/2}·Σ_{t>T0}û_t` for one-sided tests (grounded in **Remark 1**'s statistic freedom + eq 2; **the paper has no §7** — the prior "(§7 signed-t)" wording here was a cross-contamination from the Firpo-Possebom scope sentence above), with `q=1` enforced (the signed statistic has no norm order) and CI inversion reporting half-line HULLS (`[lower,+inf)` / `(-inf,upper]`; the infinite side is attached because `p → 1` in the accepted direction, and the finite endpoint is the accepted-set hull — Algorithm 1 has no monotonicity guarantee after per-candidate proxy refits, so `contiguous` is recomputed against the final reported ray (a rejected scanned point anywhere inside it — including above the finite hull edge — flips `contiguous=False`); an all-rejected grid reports the infinite side + an uncertified NaN finite endpoint as `grid_limited`); `covariates=[...]` stacks pivoted-variable matching rows into the constrained-LS proxy objective per the paper's note after eq 6 ("straightforward to incorporate (transformations of) covariates X_jt into the estimation problems (4) and (6)") — rows stack RAW (the paper's "(transformations of)" delegates scaling to the caller; documented in the docstring), residuals/statistic stay outcome-only, and exchangeability is preserved (covariate rows are fixed features of every permuted dataset; the stacked objective is invariant to time permutations of Z). **Exogeneity assumption (UserWarning on every covariate run):** `Z(θ0)` null-imputes only the OUTCOME — covariates enter as their observed values, which the framework treats as no-intervention paths (`X^N`), so covariates must be unaffected by the policy; a policy-affected covariate silently conditions the proxy and p-value on treated post-period intervention information. On the block-collapsed average-effect surface, covariate rows collapse with the same T*-block structure. Locked by `TestConformalOneSided` (hand-rolled signed-statistic oracle, directional rejection, half-line endpoints) + `TestConformalCovariates` (weight-mechanism, outcome-only residuals, permutation floor). The AR/innovation-permutation path (Lemmas 5–7, for time-series proxies) remains **deferred** (tracked in TODO.md). **Notes / deviations:** - **Note:** The standardization divisor `divisor = sqrt(apply(cbind(X0,X1), 1, var))` (per-predictor SD over donors+treated, ddof=1) and the inner/outer optimizer are **not specified in ADH 2010** (which defers these numerics to Abadie & Gardeazabal 2003 App. B / the `Synth` software). The divisor is pinned from the R `Synth::synth` source; `solution.v` lives in this scaled predictor space, so the deterministic R-parity test feeds `custom_v` in the same scaled space. @@ -2489,7 +2489,7 @@ Classic synthetic control (donor/unit weights only) for a single treated unit, d - [x] Leave-one-out donor robustness (`leave_one_out()`, ADH 2015 §4): per-drop ATT / `delta_att` table + overlay gaps; fail-closed. - [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`. +- [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. *Implemented (2026-07):* one-sided alternatives (Remark 1 signed statistic — the paper has no §7) + proxy covariates (eq 4/6 note). *Deferred:* AR/innovation-permutation path (Lemmas 5–7) — 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/tests/test_methodology_synthetic_control.py b/tests/test_methodology_synthetic_control.py index 98090375..8024ec9e 100644 --- a/tests/test_methodology_synthetic_control.py +++ b/tests/test_methodology_synthetic_control.py @@ -4563,3 +4563,304 @@ def test_sparse_non_integer_sizes_raise(): # A valid numpy-int size still works. tab = res.sparse_synthetic_control(sizes=[np.int64(2)]) assert set(tab[tab["status"] == "ran"]["size"]) == {2} + + +# =========================================================================== +# CWZ conformal extensions — one-sided alternatives (Remark 1) + covariates +# in the proxy (eq 4 note: "straightforward to incorporate (transformations +# of) covariates X_jt into the estimation problems (4) and (6)") +# =========================================================================== + + +class TestConformalOneSided: + """One-sided alternatives use the SIGNED average-effect statistic + S(u) = T*^{-1/2} sum_{t>T0} u_t — grounded in CWZ Remark 1's statistic + freedom (the paper has NO section 7; the prior TODO/REGISTRY "(§7)" + citation was a Firpo-Possebom cross-contamination, fixed in this PR).""" + + def test_signed_pvalue_matches_hand_rolled_oracle(self): + """Paper-formula oracle: eq 2 with the signed statistic, full + moving-block enumeration, computed by hand.""" + from diff_diff.conformal import _cwz_pvalue, _make_perms + + rng = np.random.default_rng(7) + T, T_star = 12, 3 + u = rng.normal(size=T) + post_mask = np.zeros(T, dtype=bool) + post_mask[-T_star:] = True + perms = _make_perms(T, "moving_block", 0, rng) + p, s_obs, n = _cwz_pvalue(u, post_mask, perms, 1, alternative="greater") + # hand-rolled: all T cyclic shifts, signed sum over the FIXED post slots + stats = [] + for j in range(T): + u_pi = u[(np.arange(T) + j) % T] + stats.append(u_pi[post_mask].sum() / np.sqrt(T_star)) + s0 = u[post_mask].sum() / np.sqrt(T_star) + p_oracle = np.mean(np.array(stats) >= s0 - 1e-12 * max(abs(s0), 1)) + assert n == T + np.testing.assert_allclose(s_obs, s0, rtol=1e-12) + np.testing.assert_allclose(p, p_oracle, rtol=1e-12) + # "less" is the mirrored statistic + p_less, s_less, _ = _cwz_pvalue(u, post_mask, perms, 1, alternative="less") + np.testing.assert_allclose(s_less, -s0, rtol=1e-12) + + def test_directional_rejection_and_q_guard(self): + """A strongly positive true effect: H1 'greater' rejects at the + two-sided level or better; H1 'less' never rejects. q != 1 with a + one-sided alternative raises (the signed statistic has no norm + order).""" + res = _fit_for_conformal(effect=4.0) + p2 = float(res.conformal_test(0.0)["p_value"]) + pg = float(res.conformal_test(0.0, alternative="greater")["p_value"]) + pl = float(res.conformal_test(0.0, alternative="less")["p_value"]) + assert pg <= p2 + 1e-12 + assert pl > 0.5 + with pytest.raises(ValueError, match="one-sided"): + res.conformal_test(0.0, q=2, alternative="greater") + with pytest.raises(ValueError, match="alternative"): + res.conformal_test(0.0, alternative="bogus") + + def test_one_sided_average_ci_is_half_line(self): + """'greater' inverts to [lower, +inf) and 'less' to (-inf, upper]; + the finite endpoints bracket the two-sided interval.""" + res = _fit_for_conformal(effect=4.0) + kw = dict(alpha=0.2, scheme="iid", n_iid=400, seed=3) + two = res.conformal_average_effect(**kw) + hi = res.conformal_average_effect(alternative="greater", **kw) + lo = res.conformal_average_effect(alternative="less", **kw) + assert np.isinf(hi["upper"]) and hi["upper"] > 0 + assert np.isinf(lo["lower"]) and lo["lower"] < 0 + assert hi["lower"] >= two["lower"] - 1e-9 + assert lo["upper"] <= two["upper"] + 1e-9 + assert hi["status"] in ("ran", "grid_limited") + + def test_one_sided_missed_ray_is_not_empty(self): + """Round-1 P1: an all-rejected one-sided grid (narrow bounds= that + miss the accepted ray) must NOT report a false empty set — the ray + always exists beyond the grid on the accepted side, so the result is + the infinite side + an UNCERTIFIED (NaN) finite endpoint with + grid_limited status. Both directions, both CI surfaces.""" + # T*=2 with T0=16 -> 9 collapsed blocks: the single-post-slot tie + # structure floors p at 1/n_blocks, so a 4-block fixture can never + # reject at alpha=0.2 — this one can (floor 1/9). + res = _fit_for_conformal(T=18, T0=16, effect=4.0) + kw = dict(alpha=0.2, scheme="iid", n_iid=400, seed=3) + # true average effect ~4: mild bounds below the accepted ray (still + # within FW convergence range — non-converged points are + # indeterminate-KEPT, which would mask the empty-grid path) + hi = res.conformal_average_effect(alternative="greater", bounds=(0.5, 1.5), **kw) + assert np.isinf(hi["upper"]) and hi["upper"] > 0 + assert np.isnan(hi["lower"]) + assert hi["status"] == "grid_limited" + lo = res.conformal_average_effect(alternative="less", bounds=(6.5, 7.5), **kw) + assert np.isinf(lo["lower"]) and lo["lower"] < 0 + assert np.isnan(lo["upper"]) + assert lo["status"] == "grid_limited" + cis = res.conformal_confidence_intervals(alternative="greater", bounds=(0.5, 1.5), **kw) + assert bool(np.isinf(cis["upper"]).all()) + assert bool(np.isnan(cis["lower"]).all()) + assert set(cis["status"]) == {"grid_limited"} + # pointwise metadata echoes the new params (round-1 P2) + assert res.conformal_inference["alternative"] == "greater" + assert res.conformal_inference["covariates"] is None + + def test_one_sided_noncontiguous_pattern_keeps_flag(self): + """Round-2 P1: Algorithm 1 has no monotonicity guarantee, so an + accepted/rejected/accepted one-sided grid must keep the hull + convention WITH contiguous=False preserved (a rejected pocket inside + the reported ray is disclosed, not hidden). Deterministic via a + crafted p-sequence.""" + import unittest.mock + + import diff_diff.conformal as cf + + res = _fit_for_conformal(T=18, T0=16, effect=4.0) + # p pattern over the 40-point grid: accepted / rejected pocket / accepted + pattern = [0.5] * 10 + [0.01] * 5 + [0.5] * 25 + calls = {"i": 0} + real = cf._cwz_pvalue + + def _scripted(u, post_mask, perms, q, alternative="two-sided"): + _, s_obs, n = real(u, post_mask, perms, q, alternative=alternative) + p = pattern[min(calls["i"], len(pattern) - 1)] + calls["i"] += 1 + return p, s_obs, n + + with unittest.mock.patch.object(cf, "_cwz_pvalue", _scripted): + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + out = res.conformal_average_effect( + alpha=0.2, + alternative="greater", + scheme="iid", + n_iid=400, + seed=3, + bounds=(0.0, 8.0), + n_grid=40, + ) + assert np.isinf(out["upper"]) + assert np.isfinite(out["lower"]) + assert out["contiguous"] is False or out["contiguous"] == False # noqa: E712 + + def test_one_sided_rejected_toward_infinity_flips_contiguous(self): + """Round-3 P1: rejected scanned points BETWEEN the finite hull edge + and the attached infinity are inside the final reported ray — the + contiguous flag must be recomputed against that ray, not the finite + hull. Accepted-pocket-then-rejected-to-the-edge, both directions.""" + import unittest.mock + + import diff_diff.conformal as cf + + res = _fit_for_conformal(T=18, T0=16, effect=4.0) + real = cf._cwz_pvalue + + def run(pattern, alternative): + calls = {"i": 0} + + def _scripted(u, post_mask, perms, q, alternative="two-sided"): + _, s_obs, n = real(u, post_mask, perms, q, alternative=alternative) + pv = pattern[min(calls["i"], len(pattern) - 1)] + calls["i"] += 1 + return pv, s_obs, n + + with unittest.mock.patch.object(cf, "_cwz_pvalue", _scripted): + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + return res.conformal_average_effect( + alpha=0.2, + alternative=alternative, + scheme="iid", + n_iid=400, + seed=3, + bounds=(0.0, 8.0), + n_grid=40, + ) + + # "greater": accepted pocket at the low end, rejected all the way to + # the upper edge -> ray [lower, +inf) contains rejected points + out = run([0.5] * 10 + [0.01] * 30, "greater") + assert np.isinf(out["upper"]) and np.isfinite(out["lower"]) + assert not out["contiguous"] + # "less": mirrored + out = run([0.01] * 30 + [0.5] * 10, "less") + assert np.isinf(out["lower"]) and np.isfinite(out["upper"]) + assert not out["contiguous"] + + def test_report_propagates_alternative_and_covariates(self): + """Round-2 P1: DiagnosticReport's conformal block must carry + alternative + covariates — a one-sided p-value is uninterpretable + under the default two-sided assumption.""" + res = _fit_for_conformal(effect=4.0) + res.conformal_test(0.0, alternative="greater") + block = res.conformal_inference + assert block["alternative"] == "greater" + from diff_diff.diagnostic_report import DiagnosticReport + + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + native = DiagnosticReport(res).to_dict()["estimator_native_diagnostics"] + ci = native["conformal_inference"] + assert ci["alternative"] == "greater" + assert ci["covariates"] is None + + def test_one_sided_pointwise_cis(self): + res = _fit_for_conformal(effect=4.0) + cis = res.conformal_confidence_intervals( + alpha=0.2, scheme="iid", n_iid=400, seed=3, alternative="greater" + ) + assert bool(np.isinf(cis["upper"]).all()) + assert np.isfinite(cis["lower"]).all() + + +class TestConformalCovariates: + """Covariate-matching rows stack into the CWZ proxy objective, raw + (the paper's "(transformations of)" delegates scaling to the caller); + residuals and the permutation p-value stay outcome-only.""" + + @staticmethod + def _fit_with_covariate(effect=4.0, seed=2, cov_signal=True): + df, years, T0 = _make_panel(n_donors=5, T=18, T0=14, effect=effect, seed=seed) + rng = np.random.default_rng(seed + 100) + # EXOGENOUS covariate: unit-level structure + common trend, NOT + # derived from the (effect-carrying) outcome — the CWZ covariate + # assumption is that X is unaffected by the intervention (round-4 + # review P1: a y-derived fixture would bless post-treatment leakage). + if cov_signal: + unit_level = { + u: h for u, h in zip(df["unit"].unique(), rng.normal(0, 1.0, df["unit"].nunique())) + } + df["z"] = ( + df["unit"].map(unit_level) + + 0.1 * df["year"].astype(float) + + rng.normal(0, 0.05, len(df)) + ) + else: + df["z"] = rng.normal(0, 1.0, len(df)) + res = SyntheticControl(**_FAST).fit( + df, + "y", + "treated", + "unit", + "year", + post_periods=years[T0:], + treated_unit="treated", + predictors=["z"], + ) + return res + + def test_covariate_rows_change_the_proxy_weights(self): + """Mechanism check: a covariate with independent signal must be able + to move the proxy weights (the stacked rows enter the objective).""" + from diff_diff.conformal import _cwz_proxy_fit + + rng = np.random.default_rng(11) + T, J = 20, 4 + Y0 = rng.normal(size=(T, J)) + y1 = Y0[:, 0] * 0.7 + Y0[:, 1] * 0.3 + rng.normal(0, 0.05, T) + # covariate says the target matches donor 2 instead + X0 = rng.normal(size=(T, J)) + x1 = X0[:, 2] + rng.normal(0, 0.01, T) + w_plain, resid_plain, _ = _cwz_proxy_fit(y1, Y0, max_iter=2000, min_decrease=1e-10) + w_cov, resid_cov, _ = _cwz_proxy_fit( + y1, Y0, max_iter=2000, min_decrease=1e-10, x1_rows=x1 * 10, X0_rows=X0 * 10 + ) + assert np.max(np.abs(w_cov - w_plain)) > 0.05 + assert w_cov[2] > w_plain[2] + # residuals are OUTCOME rows only (length T, not T + covariate rows) + assert resid_cov.shape == (T,) + np.testing.assert_allclose(resid_cov, y1 - Y0 @ w_cov, rtol=1e-12) + + def test_covariate_surface_and_validation(self): + res = self._fit_with_covariate() + with pytest.warns(UserWarning, match="UNAFFECTED by the intervention"): + s = res.conformal_test(0.0, covariates=["z"]) + assert 0.0 < float(s["p_value"]) <= 1.0 + assert res.conformal_inference["covariates"] == ["z"] + with pytest.raises(ValueError, match="not in the fit snapshot"): + res.conformal_test(0.0, covariates=["nope"]) + with pytest.raises(ValueError, match="outcome"): + res.conformal_test(0.0, covariates=["y"]) + + @pytest.mark.filterwarnings("ignore::UserWarning") + def test_covariates_compose_with_one_sided_and_cis(self): + res = self._fit_with_covariate(effect=4.0) + pg = float(res.conformal_test(0.0, alternative="greater", covariates=["z"])["p_value"]) + assert 0.0 < pg <= 1.0 + avg = res.conformal_average_effect( + alpha=0.2, scheme="iid", n_iid=300, seed=3, alternative="greater", covariates=["z"] + ) + assert np.isinf(avg["upper"]) + cis = res.conformal_confidence_intervals( + alpha=0.2, scheme="iid", n_iid=300, seed=3, covariates=["z"] + ) + assert np.isfinite(cis["lower"]).all() and np.isfinite(cis["upper"]).all() + + def test_permutation_floor_still_holds_with_covariates(self): + """The p-value floor 1/|Pi| (identity in Pi) is statistic- and + proxy-independent — exchangeability is preserved because covariate + rows are fixed features of every permuted dataset.""" + res = self._fit_with_covariate(effect=0.0) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + s = res.conformal_test(0.0, covariates=["z"], scheme="moving_block") + assert float(s["p_value"]) >= 1.0 / float(s["n_perms"]) - 1e-12 From 0d7ebe5a4714efd333514657d7ca1bc37aa7162a Mon Sep 17 00:00:00 2001 From: igerber Date: Fri, 10 Jul 2026 08:06:02 -0400 Subject: [PATCH 2/2] docs(registry): drop the contradictory outcomes-only phrase from the proxy bullet Review P3: the bullet both allowed covariate rows in the objective and said outcomes-only; the invariant is that residuals/statistic are outcome-only regardless of covariate rows. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01Lbd6nqWmg4snvvBmegwqiw --- docs/methodology/REGISTRY.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index ddfc0f28..f9ec0121 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -2435,7 +2435,7 @@ Classic synthetic control (donor/unit weights only) for a single treated unit, d - **`confidence_set(family="constant"|"linear", gamma=0.1, bounds=None, n_grid=200)`** inverts that test over a one-parameter family: `"constant"` → `f(t)=c` (Eq 15), a confidence **interval** for a constant effect (Eq 16); `"linear"` → `f(t)=c̃·(t−T0)` with the 1-based post-period index (Eq 17), a confidence **set** over the slope `c̃` (Eq 18). Membership is the paper's **strict** `p^c > γ` (Eq 14 — see the boundary Note). With `bounds=None` the set is recovered **EXACTLY**: `p^c` is a piecewise-constant step function (each placebo's indicator flips only at the real roots of `A_j(c)·D_1 = A_1(c)·D_j`, a quadratic in `c`), so the placebo breakpoints partition the line and `p` is evaluated once per induced open interval AND at each breakpoint — where a tie under `≥` can lift `p` above γ, so an isolated accepted singleton (a tangent / co-located root) is captured. The accepted set is the union of accepted intervals/points, with **no centering or monotonicity assumption** (a poor-pre-fit treated unit can have its accepted region in the tails, not around the point estimate). `bounds=(lo,hi)` instead scans a fixed grid (grid-limited membership). The summary is on `effect_confidence_set` (`status ∈ {"ran","empty","unbounded"}`) and the full grid on `get_confidence_set_df()`. **Fail-closed:** `γ < 1/(J+1)` ⇒ `p^c > γ` for every `c` ⇒ `"unbounded"` (`±inf` endpoints + warning — the discrete-granularity point, fn 8); a treated unit lacking the best pre-fit can give a one-sided unbounded edge; if no interval or breakpoint is accepted the set is `"empty"` (NaN endpoints); a non-contiguous accepted region (disjoint components / an isolated singleton) reports the `[lower, upper]` hull with `contiguous=False` + a warning; and `< 2` donors / a non-converged treated fit / an unpickled result (no placebo reference set) raise `ValueError`. **Scope:** sensitivity weights (`φ≠0`, Eqs 7–9), the general test-statistic menu `θ¹–θ⁵` (Eq 19), one-sided (§7's signed-`t` statistic), and the multiple-outcome / multiple-treated extensions (§6) are **deferred** (flagged in the paper review checklist). **Conformal inference (Chernozhukov-Wüthrich-Zhu 2021, opt-in):** Valid p-values for a hypothesized post-period effect trajectory and pointwise confidence intervals — what the placebo / test-inversion paths cannot give. Surfaced under `estimator_native_diagnostics`; the analytical `se`/`t_stat`/`p_value`/`conf_int` stay NaN (separate permutation object — see the non-analytical Note). Under a sharp null `H0: θ = θ0` the counterfactual treated outcomes are imputed (`Y^N_{1t} = Y_{1t} − θ0_t`, `t>T0`; pre unchanged), a **time-permutation-invariant proxy** is fit UNDER THE NULL on **all** periods, and the statistic `S_q(û) = ((1/√T*)·Σ_{t>T0}|û_t|^q)^{1/q}` (CWZ §2.2) is referred to its permutation distribution (eq 2) by reshuffling residuals over time: -- **Proxy (CWZ §2.3, eqs 3–4):** the canonical constrained-LS synthetic control (outcomes by default; optional RAW covariate-matching rows stack into the same objective when `covariates=` is supplied — see the scope note) — simplex weights `argmin_w Σ_{t=1}^T (Y^N_{1t} − Σ_j w_j Y^N_{jt})²` s.t. `w≥0, Σw=1`, **NO V-matrix, no intercept, outcomes-only**, fit over ALL periods (footnote 9: "we estimate w under the null based on all the data"). Reuses the Frank-Wolfe simplex solver `utils._sc_weight_fw` (packed `(T,J+1)`, `zeta=0, intercept=False`). +- **Proxy (CWZ §2.3, eqs 3–4):** the canonical constrained-LS synthetic control (outcomes by default; optional RAW covariate-matching rows stack into the same objective when `covariates=` is supplied — see the scope note) — simplex weights `argmin_w Σ_{t=1}^T (Y^N_{1t} − Σ_j w_j Y^N_{jt})²` s.t. `w≥0, Σw=1`, **NO V-matrix, no intercept** (residuals/statistic are outcome-only regardless of covariate rows), fit over ALL periods (footnote 9: "we estimate w under the null based on all the data"). Reuses the Frank-Wolfe simplex solver `utils._sc_weight_fw` (packed `(T,J+1)`, `zeta=0, intercept=False`). - **`conformal_test(effect, q=1, alternative="two-sided", covariates=None, scheme="moving_block", n_iid=10000, seed=None)`** — joint sharp-null test (eqs 1–2). `effect` is a scalar (constant) or length-`n_post` path; `q ∈ {1,2,∞}` (`q=∞` is `max|û_t|`). The proxy is fit ONCE under the null and only residuals are permuted (footnote 7 — permuting residuals ≡ permuting data for a time-permutation-invariant proxy). Returns `p_value`, `S_observed`, `n_perms`, `proxy_converged`. - **`conformal_confidence_intervals(alpha=0.1, alternative="two-sided", covariates=None, scheme="moving_block", n_iid=10000, bounds=None, n_grid=100, seed=None)`** — pointwise per-period CIs by test inversion (Algorithm 1). Per the paper (§2.2), each per-period CI for period `t` uses `Z = (Z_1,…,Z_{T0}, Z_t)` — the pre-periods PLUS only period `t`, the **other post-periods dropped** — a clean `T*=1` conformal test (`q` is therefore inert: `S_q = |û_t|`). Returns one row per post period (`lower`/`upper`/`point_estimate`/`status`/`contiguous`); the full grid is on `get_conformal_grid_df()`. - **`conformal_average_effect(alpha=0.1, alternative="two-sided", covariates=None, scheme="moving_block", n_iid=10000, bounds=None, n_grid=200, seed=None)`** — CI for the average effect `T*^{-1}Σ_{t>T0}θ_t` by test inversion (Appendix A.1): **collapse** the panel into non-overlapping `T*`-blocks (per-unit block averages), fit the proxy on the collapsed `Z̄`, permute the `T/T*` **block** residuals. The earliest `T0 mod T*` pre-periods are dropped to make the block count integral. Requires `T0 ≥ T*`.