diff --git a/CHANGELOG.md b/CHANGELOG.md index 669aafd7..d046b829 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] ### Added +- **`StackedDiD` covariate balancing (CBWSDID; Ustyuzhanin 2026, arXiv:2604.02293).** New constructor parameter `balance="entropy"` plus `fit(..., covariates=[...])` add a within-sub-experiment design stage: entropy balancing (Hainmueller 2012) reweights the clean controls toward the treated cohort's covariate means (read at the last pre-treatment period), and the resulting design weights `b_sa` compose with the Wing et al. (2024) corrective weights via the effective control mass into the final stacked weights `W_sa`. This is **control-only reweighting**, so it estimates untreated trends under *conditional* parallel trends while preserving the trimmed-aggregate-ATT estimand (at `b_sa=1` it reduces to the paper's unit-count weighted stacked DID, equal to `StackedDiD(weighting="aggregate")` on balanced event windows). Inference reuses the existing conditional-on-weights cluster-robust path. Scope: requires `weighting="aggregate"` and **balanced event windows** (ragged windows raise — the unit-count vs observation-count convention is unresolved off balanced panels); `population`/`sample_share`/`survey_design=` and matching-based balancing / the repeated-treatment extension are not supported (raise `NotImplementedError`). Infeasible cohorts fail closed with a clear error. New `diff_diff/balancing.py` (entropy-balancing solver). Estimand validated end-to-end against the closed-form CBWSDID formula (`tests/test_methodology_stacked_did.py`). - **`SyntheticControl` conformal inference (Chernozhukov, Wüthrich & Zhu 2021, *JASA* 116(536)).** Three opt-in `SyntheticControlResults` methods give valid p-values for the post-period effect trajectory and pointwise confidence intervals — what the in-space placebo / Firpo-Possebom test-inversion paths cannot. Unlike the Firpo path (which re-ranks the cross-unit placebo gaps), the conformal layer fits its **own** time-permutation-invariant constrained-LS synthetic-control proxy (CWZ §2.3 eqs 3–4 — simplex weights on raw outcomes over **all** periods under the null, no `V`-matrix, no intercept) and permutes residuals **over time** for the single treated unit (CWZ's exactness theory requires a time-symmetric proxy, which the headline ADH `V`-matrix fit is not). **`conformal_test(effect, q=1, scheme="moving_block", n_iid=10000, seed=None)`** computes the joint sharp-null permutation p-value (eqs 1–2) of `S_q(û) = ((1/√T*)·Σ_{t>T0}|û_t|^q)^{1/q}` (`q ∈ {1, 2, ∞}`); the proxy is fit once and only residuals are permuted (footnote 7). **`conformal_confidence_intervals(alpha=0.1, scheme="moving_block", bounds=None, n_grid=100, seed=None)`** returns pointwise per-period CIs by test inversion (Algorithm 1 — each period `t` uses `Z = (pre-periods, t)` with the other post-periods dropped, a clean `T*=1` test). **`conformal_average_effect(alpha=0.1, scheme="moving_block", bounds=None, n_grid=200, seed=None)`** returns a CI for the average post-period effect by collapsing the panel into non-overlapping `T*`-blocks and permuting the block residuals (Appendix A.1). Permutation schemes: `"moving_block"` (`Π_→` cyclic shifts, valid under serial dependence — the default) and `"iid"` (`Π_all`, sampled, finer p-values); both include the identity so the p-value floor is `1/|Π|` (no extra `+1`). Fail-closed handling for `<1` donor / unpickled result / non-finite panel / non-converged grid points (treated as indeterminate, not rejected) / grid-limited / empty / unbounded sets; a single donor and `T*≥T0` warn. Surfaced under `conformal_inference` / `get_conformal_grid_df()` and `DiagnosticReport`'s `estimator_native_diagnostics`; the analytical `se`/`t_stat`/`p_value`/`conf_int`/`is_significant` stay NaN throughout. Core in the new `diff_diff/conformal.py` (reuses the Frank-Wolfe simplex solver). *Deferred:* one-sided variants (§7), covariates folded into the proxy, and the AR/innovation-permutation path (Lemmas 5–7). ### Changed diff --git a/README.md b/README.md index ae7b5cfb..3a3a26b1 100644 --- a/README.md +++ b/README.md @@ -112,7 +112,7 @@ Full guide: `diff_diff.get_llm_guide("practitioner")`. - [TripleDifference](https://diff-diff.readthedocs.io/en/stable/api/triple_diff.html) - triple difference (DDD) estimator for designs requiring two criteria for treatment eligibility - [ContinuousDiD](https://diff-diff.readthedocs.io/en/stable/api/continuous_did.html) - Callaway, Goodman-Bacon & Sant'Anna (2024) continuous treatment DiD with dose-response curves - [HeterogeneousAdoptionDiD](https://diff-diff.readthedocs.io/en/stable/api/had.html) - de Chaisemartin, Ciccia, D'Haultfœuille & Knau (2026) for designs where **no unit remains untreated**; local-linear estimator at the dose support boundary returning Weighted Average Slope (WAS) on Design 1' (`d̲ = 0` / QUG) or `WAS_{d̲}` on Design 1 (`d̲ > 0`, continuous-near-d̲ or mass-point), with a multi-period event-study extension (last-treatment cohort, pointwise CIs). **Panel-only** in this release - repeated cross-sections rejected by the validator. Alias `HAD`. -- [StackedDiD](https://diff-diff.readthedocs.io/en/stable/api/stacked_did.html) - Wing, Freedman & Hollingsworth (2024) stacked DiD with Q-weights and sub-experiments +- [StackedDiD](https://diff-diff.readthedocs.io/en/stable/api/stacked_did.html) - Wing, Freedman & Hollingsworth (2024) stacked DiD with Q-weights and sub-experiments; optional covariate balancing (Ustyuzhanin 2026) - [EfficientDiD](https://diff-diff.readthedocs.io/en/stable/api/efficient_did.html) - Chen, Sant'Anna & Xie (2025) efficient DiD with optimal weighting for tighter SEs - [TROP](https://diff-diff.readthedocs.io/en/stable/api/trop.html) - Triply Robust Panel estimator (Athey et al. 2025) with nuclear norm factor adjustment - [StaggeredTripleDifference](https://diff-diff.readthedocs.io/en/stable/api/staggered.html#staggeredtripledifference) - Ortiz-Villavicencio & Sant'Anna (2025) staggered DDD with group-time ATT diff --git a/TODO.md b/TODO.md index d7903c69..69497c39 100644 --- a/TODO.md +++ b/TODO.md @@ -74,6 +74,7 @@ Deferred items from PR reviews that were not addressed before merge. | Issue | Location | PR | Priority | |-------|----------|----|----------| +| CBWSDID covariate balancing (`StackedDiD(balance="entropy")`) v1 supports only balanced event windows + `weighting="aggregate"`; unbalanced/ragged panels fail closed (the unit-count vs observation-count corrector convention is unresolved off balanced panels). Matching-based balancing and the repeated `0→1`/`1→0` episode extension are also deferred (out-of-scope guards raise). Documented in REGISTRY.md StackedDiD "Covariate balancing (CBWSDID)" Notes. | `stacked_did.py`, `balancing.py`, `docs/methodology/REGISTRY.md` | follow-up | Low | | `SyntheticControl` cv: `in_space_placebo()` / `leave_one_out()` report a cv refit excluded for STRUCTURAL infeasibility (donor-indistinguishable re-aggregated window) with the generic `status="failed"` — same machine-readable status as a genuine inner-solver non-convergence. The failure warnings now distinguish the two causes (and the correct remediation) under cv, and `in_time_placebo()` already splits structural→`"infeasible"` vs `"failed"`, but in-space/LOO do not yet emit a separate machine-readable status/reason-code. Thread a reason code from `_outer_solve_V_cv()`/`_placebo_fit_unit()` and add an `"infeasible"` status + count to the in-space/LOO outputs (mirror the in-time split). | `synthetic_control.py`, `synthetic_control_results.py` | follow-up | Low | | dCDH: Phase 1 per-period placebo DID_M^pl has NaN SE (no IF derivation for the per-period aggregation path). Multi-horizon placebos (L_max >= 1) have valid SE. | `chaisemartin_dhaultfoeuille.py` | #294 | Low | | dCDH: Survey cell-period allocator's post-period attribution is a library convention, not derived from the observation-level survey linearization. MC coverage is empirically close to nominal on the test DGP; a formal derivation (or a covariance-aware two-cell alternative) is deferred. Documented in REGISTRY.md survey IF expansion Note. | `chaisemartin_dhaultfoeuille.py`, `docs/methodology/REGISTRY.md` | #408 | Medium | diff --git a/benchmarks/R/generate_cbwsdid_golden.R b/benchmarks/R/generate_cbwsdid_golden.R new file mode 100644 index 00000000..49a6c008 --- /dev/null +++ b/benchmarks/R/generate_cbwsdid_golden.R @@ -0,0 +1,52 @@ +#!/usr/bin/env Rscript +# Generate the cross-language golden fixture for StackedDiD's covariate-balancing +# (CBWSDID) path against the reference R package `cbwsdid` (Ustyuzhanin 2026). +# +# Unlike generate_stacked_did_golden.R (which operates on a PRE-stacked CSV so the +# R side is independent of Python stacking logic), `cbwsdid` does its OWN stacking +# + balancing, so this harness hands it the raw panel and dumps the dynamic +# event-study ATTs. The Python side (StackedDiD(balance="entropy", ...)) reproduces +# them via its independent entropy-balancing solver + effective-mass W_sa. +# +# Refinement: refinement.method="weightit", method="ebal" = entropy balancing +# (Hainmueller 2012) on covs.formula=~x, matching StackedDiD(balance="entropy", +# covariates=["x"]). Install: remotes::install_github("vadvu/cbwsdid"). +# +# Usage: Rscript benchmarks/R/generate_cbwsdid_golden.R + +suppressMessages({ + library(cbwsdid) + library(jsonlite) +}) + +# Run from the repository root: Rscript benchmarks/R/generate_cbwsdid_golden.R +panel_csv <- "benchmarks/data/cbwsdid_balance_panel.csv" +out_json <- "benchmarks/data/cbwsdid_golden.json" + +df <- read.csv(panel_csv) + +m <- cbwsdid( + data = df, y = "y", d = "d", id = c("unit", "time"), + kappa = c(-2, 2), design = "absorbing", post_path = "stable", + refinement.method = "weightit", covs.formula = ~x, + refinement.args = list(method = "ebal"), pooled = TRUE +) +qoi <- cbwsdid_qoi(m, type = "dynamic") + +golden <- list( + meta = list( + package = "cbwsdid", + R_version = R.version.string, + panel = "benchmarks/data/cbwsdid_balance_panel.csv", + estimator = "cbwsdid(design='absorbing', refinement.method='weightit', method='ebal', covs.formula=~x)", + kappa = c(-2L, 2L) + ), + dynamic = list( + event_time = as.integer(qoi$et), + estimate = as.numeric(qoi$estimate), + std_error = as.numeric(qoi$std.error) + ) +) +write_json(golden, out_json, auto_unbox = TRUE, digits = 15, pretty = TRUE) +cat("wrote", out_json, "\n") +print(data.frame(et = qoi$et, estimate = qoi$estimate, se = qoi$std.error)) diff --git a/benchmarks/data/cbwsdid_balance_panel.csv b/benchmarks/data/cbwsdid_balance_panel.csv new file mode 100644 index 00000000..ccbb0fb8 --- /dev/null +++ b/benchmarks/data/cbwsdid_balance_panel.csv @@ -0,0 +1,1597 @@ +unit,time,first_treat,x,y,d +1,1,5,0.0619784,-1.59798253,0 +1,2,5,0.0619784,-1.52627612,0 +1,3,5,0.0619784,-1.47772199,0 +1,4,5,0.0619784,-1.43535074,0 +1,5,5,0.0619784,0.5649883,1 +1,6,5,0.0619784,0.60370807,1 +1,7,5,0.0619784,0.65067989,1 +1,8,5,0.0619784,0.72669429,1 +1,9,5,0.0619784,0.7980041,1 +1,10,5,0.0619784,0.82746107,1 +1,11,5,0.0619784,0.92154498,1 +1,12,5,0.0619784,0.92395958,1 +1,13,5,0.0619784,0.93183999,1 +1,14,5,0.0619784,1.00056904,1 +2,1,5,0.65736913,0.38047965,0 +2,2,5,0.65736913,0.55596233,0 +2,3,5,0.65736913,0.83815501,0 +2,4,5,0.65736913,1.02382106,0 +2,5,5,0.65736913,3.26572812,1 +2,6,5,0.65736913,3.52920877,1 +2,7,5,0.65736913,3.76577671,1 +2,8,5,0.65736913,3.9713327,1 +2,9,5,0.65736913,4.18295648,1 +2,10,5,0.65736913,4.45709044,1 +2,11,5,0.65736913,4.69669344,1 +2,12,5,0.65736913,4.90749324,1 +2,13,5,0.65736913,5.16012842,1 +2,14,5,0.65736913,5.39255219,1 +3,1,5,1.21564633,0.33251432,0 +3,2,5,1.21564633,0.65499706,0 +3,3,5,1.21564633,1.05079385,0 +3,4,5,1.21564633,1.3506669,0 +3,5,5,1.21564633,3.65835611,1 +3,6,5,1.21564633,4.02131486,1 +3,7,5,1.21564633,4.33357149,1 +3,8,5,1.21564633,4.72262335,1 +3,9,5,1.21564633,5.03430878,1 +3,10,5,1.21564633,5.37009223,1 +3,11,5,1.21564633,5.67865626,1 +3,12,5,1.21564633,6.03268207,1 +3,13,5,1.21564633,6.39183666,1 +3,14,5,1.21564633,6.66193614,1 +4,1,5,2.6461679,0.14122941,0 +4,2,5,2.6461679,0.93350853,0 +4,3,5,2.6461679,1.69030628,0 +4,4,5,2.6461679,2.5977661,0 +4,5,5,2.6461679,5.35034265,1 +4,6,5,2.6461679,6.18673203,1 +4,7,5,2.6461679,6.9652389,1 +4,8,5,2.6461679,7.75381914,1 +4,9,5,2.6461679,8.56808641,1 +4,10,5,2.6461679,9.45542808,1 +4,11,5,2.6461679,10.28120835,1 +4,12,5,2.6461679,11.064938,1 +4,13,5,2.6461679,11.81215469,1 +4,14,5,2.6461679,12.65756657,1 +5,1,5,2.09149197,0.83846805,0 +5,2,5,2.09149197,1.38464557,0 +5,3,5,2.09149197,2.0126814,0 +5,4,5,2.09149197,2.66333951,0 +5,5,5,2.09149197,5.24974434,1 +5,6,5,2.09149197,5.86244331,1 +5,7,5,2.09149197,6.47885123,1 +5,8,5,2.09149197,7.04855436,1 +5,9,5,2.09149197,7.69673647,1 +5,10,5,2.09149197,8.28678262,1 +5,11,5,2.09149197,8.8741409,1 +5,12,5,2.09149197,9.46215496,1 +5,13,5,2.09149197,10.140564,1 +5,14,5,2.09149197,10.68088998,1 +6,1,5,1.82314769,-0.07559539,0 +6,2,5,1.82314769,0.46698753,0 +6,3,5,1.82314769,1.00120045,0 +6,4,5,1.82314769,1.61740068,0 +6,5,5,1.82314769,4.13300594,1 +6,6,5,1.82314769,4.73614101,1 +6,7,5,1.82314769,5.22815521,1 +6,8,5,1.82314769,5.82035717,1 +6,9,5,1.82314769,6.31762669,1 +6,10,5,1.82314769,6.96464463,1 +6,11,5,1.82314769,7.47195508,1 +6,12,5,1.82314769,8.00460021,1 +6,13,5,1.82314769,8.57598633,1 +6,14,5,1.82314769,9.09139995,1 +7,1,5,1.11430533,-0.67679771,0 +7,2,5,1.11430533,-0.29650171,0 +7,3,5,1.11430533,0.0334142,0 +7,4,5,1.11430533,0.37557146,0 +7,5,5,1.11430533,2.68045595,1 +7,6,5,1.11430533,3.12140448,1 +7,7,5,1.11430533,3.45869282,1 +7,8,5,1.11430533,3.79318919,1 +7,9,5,1.11430533,4.18961833,1 +7,10,5,1.11430533,4.51383963,1 +7,11,5,1.11430533,4.94423734,1 +7,12,5,1.11430533,5.2172892,1 +7,13,5,1.11430533,5.57686051,1 +7,14,5,1.11430533,5.94371125,1 +8,1,5,0.18125616,0.10172873,0 +8,2,5,0.18125616,0.11965237,0 +8,3,5,0.18125616,0.16518538,0 +8,4,5,0.18125616,0.26254506,0 +8,5,5,0.18125616,2.39479531,1 +8,6,5,0.18125616,2.49548195,1 +8,7,5,0.18125616,2.58514595,1 +8,8,5,0.18125616,2.6487408,1 +8,9,5,0.18125616,2.69642823,1 +8,10,5,0.18125616,2.80027309,1 +8,11,5,0.18125616,2.87442813,1 +8,12,5,0.18125616,2.91441279,1 +8,13,5,0.18125616,2.97055854,1 +8,14,5,0.18125616,3.10453609,1 +9,1,5,0.08424425,-0.3517667,0 +9,2,5,0.08424425,-0.29629018,0 +9,3,5,0.08424425,-0.35228664,0 +9,4,5,0.08424425,-0.28739137,0 +9,5,5,0.08424425,1.74548241,1 +9,6,5,0.08424425,1.73209131,1 +9,7,5,0.08424425,1.78792241,1 +9,8,5,0.08424425,1.83998869,1 +9,9,5,0.08424425,1.82736338,1 +9,10,5,0.08424425,1.87922359,1 +9,11,5,0.08424425,1.85391412,1 +9,12,5,0.08424425,1.96037724,1 +9,13,5,0.08424425,1.97217235,1 +9,14,5,0.08424425,1.96498664,1 +10,1,5,0.02663859,0.40738306,0 +10,2,5,0.02663859,0.4623977,0 +10,3,5,0.02663859,0.47033454,0 +10,4,5,0.02663859,0.4412097,0 +10,5,5,0.02663859,2.56857834,1 +10,6,5,0.02663859,2.5825816,1 +10,7,5,0.02663859,2.57376204,1 +10,8,5,0.02663859,2.63423355,1 +10,9,5,0.02663859,2.68622431,1 +10,10,5,0.02663859,2.68704927,1 +10,11,5,0.02663859,2.74094222,1 +10,12,5,0.02663859,2.79544999,1 +10,13,5,0.02663859,2.83966723,1 +10,14,5,0.02663859,2.80862978,1 +11,1,5,1.54096723,0.20395887,0 +11,2,5,1.54096723,0.60806005,0 +11,3,5,1.54096723,1.04498675,0 +11,4,5,1.54096723,1.5505272,0 +11,5,5,1.54096723,3.89704921,1 +11,6,5,1.54096723,4.35118279,1 +11,7,5,1.54096723,4.74459683,1 +11,8,5,1.54096723,5.22891687,1 +11,9,5,1.54096723,5.6520534,1 +11,10,5,1.54096723,6.08143056,1 +11,11,5,1.54096723,6.52405468,1 +11,12,5,1.54096723,6.93956961,1 +11,13,5,1.54096723,7.38123571,1 +11,14,5,1.54096723,7.82127972,1 +12,1,5,1.21737111,-0.64425019,0 +12,2,5,1.21737111,-0.27520081,0 +12,3,5,1.21737111,0.10637152,0 +12,4,5,1.21737111,0.43171935,0 +12,5,5,1.21737111,2.80187323,1 +12,6,5,1.21737111,3.23343585,1 +12,7,5,1.21737111,3.60106621,1 +12,8,5,1.21737111,3.96057988,1 +12,9,5,1.21737111,4.36560138,1 +12,10,5,1.21737111,4.69971872,1 +12,11,5,1.21737111,5.07748057,1 +12,12,5,1.21737111,5.46235658,1 +12,13,5,1.21737111,5.89042934,1 +12,14,5,1.21737111,6.22506572,1 +13,1,5,-0.29590772,0.81201162,0 +13,2,5,-0.29590772,0.72378966,0 +13,3,5,-0.29590772,0.64283886,0 +13,4,5,-0.29590772,0.43426942,0 +13,5,5,-0.29590772,2.33371759,1 +13,6,5,-0.29590772,2.30386932,1 +13,7,5,-0.29590772,2.09549987,1 +13,8,5,-0.29590772,1.97952887,1 +13,9,5,-0.29590772,1.85934101,1 +13,10,5,-0.29590772,1.71407616,1 +13,11,5,-0.29590772,1.68601792,1 +13,12,5,-0.29590772,1.56327924,1 +13,13,5,-0.29590772,1.43954586,1 +13,14,5,-0.29590772,1.27732065,1 +14,1,5,1.14617944,-0.27313147,0 +14,2,5,1.14617944,0.10287029,0 +14,3,5,1.14617944,0.42133981,0 +14,4,5,1.14617944,0.79289005,0 +14,5,5,1.14617944,3.14325199,1 +14,6,5,1.14617944,3.46825441,1 +14,7,5,1.14617944,3.8307704,1 +14,8,5,1.14617944,4.22061798,1 +14,9,5,1.14617944,4.52842087,1 +14,10,5,1.14617944,4.89711292,1 +14,11,5,1.14617944,5.25413536,1 +14,12,5,1.14617944,5.59029585,1 +14,13,5,1.14617944,5.95328292,1 +14,14,5,1.14617944,6.33763715,1 +15,1,5,1.24908067,0.90252178,0 +15,2,5,1.24908067,1.28892054,0 +15,3,5,1.24908067,1.64173335,0 +15,4,5,1.24908067,2.00607916,0 +15,5,5,1.24908067,4.37109301,1 +15,6,5,1.24908067,4.68709485,1 +15,7,5,1.24908067,5.09144331,1 +15,8,5,1.24908067,5.48134854,1 +15,9,5,1.24908067,5.82635493,1 +15,10,5,1.24908067,6.14064878,1 +15,11,5,1.24908067,6.47909026,1 +15,12,5,1.24908067,6.85417136,1 +15,13,5,1.24908067,7.22920883,1 +15,14,5,1.24908067,7.54099748,1 +16,1,5,1.19172283,1.4971106,0 +16,2,5,1.19172283,1.90292522,0 +16,3,5,1.19172283,2.18463403,0 +16,4,5,1.19172283,2.57604714,0 +16,5,5,1.19172283,4.93348906,1 +16,6,5,1.19172283,5.31454143,1 +16,7,5,1.19172283,5.61854664,1 +16,8,5,1.19172283,5.94915643,1 +16,9,5,1.19172283,6.32185647,1 +16,10,5,1.19172283,6.68870333,1 +16,11,5,1.19172283,6.99862151,1 +16,12,5,1.19172283,7.34532352,1 +16,13,5,1.19172283,7.73179207,1 +16,14,5,1.19172283,7.99747907,1 +17,1,5,-0.18890405,1.3213682,0 +17,2,5,-0.18890405,1.25468992,0 +17,3,5,-0.18890405,1.11539778,0 +17,4,5,-0.18890405,1.02036195,0 +17,5,5,-0.18890405,2.85781751,1 +17,6,5,-0.18890405,2.80228872,1 +17,7,5,-0.18890405,2.69666259,1 +17,8,5,-0.18890405,2.63148312,1 +17,9,5,-0.18890405,2.41774276,1 +17,10,5,-0.18890405,2.41293205,1 +17,11,5,-0.18890405,2.27294063,1 +17,12,5,-0.18890405,2.1983224,1 +17,13,5,-0.18890405,2.06246496,1 +17,14,5,-0.18890405,2.0360752,1 +18,1,5,1.97559081,2.5860472,0 +18,2,5,1.97559081,3.20213157,0 +18,3,5,1.97559081,3.80670198,0 +18,4,5,1.97559081,4.28730135,0 +18,5,5,1.97559081,6.88137861,1 +18,6,5,1.97559081,7.56215858,1 +18,7,5,1.97559081,8.08686109,1 +18,8,5,1.97559081,8.64908739,1 +18,9,5,1.97559081,9.30321083,1 +18,10,5,1.97559081,9.84038631,1 +18,11,5,1.97559081,10.38992486,1 +18,12,5,1.97559081,11.02526776,1 +18,13,5,1.97559081,11.56581042,1 +18,14,5,1.97559081,12.14888286,1 +19,1,7,-1.769458,-0.93089824,0 +19,2,7,-1.769458,-1.48986803,0 +19,3,7,-1.769458,-1.85729208,0 +19,4,7,-1.769458,-2.46232227,0 +19,5,7,-1.769458,-2.8830552,0 +19,6,7,-1.769458,-3.36454725,0 +19,7,7,-1.769458,-1.87285446,1 +19,8,7,-1.769458,-2.34856075,1 +19,9,7,-1.769458,-2.87342149,1 +19,10,7,-1.769458,-3.34596513,1 +19,11,7,-1.769458,-3.87138329,1 +19,12,7,-1.769458,-4.37701319,1 +19,13,7,-1.769458,-4.84568211,1 +19,14,7,-1.769458,-5.35545896,1 +20,1,7,1.17920186,0.6962865,0 +20,2,7,1.17920186,1.06080703,0 +20,3,7,1.17920186,1.36579228,0 +20,4,7,1.17920186,1.70319799,0 +20,5,7,1.17920186,2.00781164,0 +20,6,7,1.17920186,2.36451747,0 +20,7,7,1.17920186,4.67251961,1 +20,8,7,1.17920186,5.10201071,1 +20,9,7,1.17920186,5.39791445,1 +20,10,7,1.17920186,5.76742859,1 +20,11,7,1.17920186,6.06707457,1 +20,12,7,1.17920186,6.50193806,1 +20,13,7,1.17920186,6.80653435,1 +20,14,7,1.17920186,7.10740768,1 +21,1,7,0.66879027,-0.71638375,0 +21,2,7,0.66879027,-0.52127945,0 +21,3,7,0.66879027,-0.34949208,0 +21,4,7,0.66879027,-0.15457919,0 +21,5,7,0.66879027,0.06966585,0 +21,6,7,0.66879027,0.2502128,0 +21,7,7,0.66879027,2.49806362,1 +21,8,7,0.66879027,2.65262735,1 +21,9,7,0.66879027,2.7729355,1 +21,10,7,0.66879027,3.04782087,1 +21,11,7,0.66879027,3.2272434,1 +21,12,7,0.66879027,3.43199623,1 +21,13,7,0.66879027,3.5809193,1 +21,14,7,0.66879027,3.77352629,1 +22,1,7,0.69480967,3.91333397,0 +22,2,7,0.69480967,4.09072096,0 +22,3,7,0.69480967,4.30570233,0 +22,4,7,0.69480967,4.522243,0 +22,5,7,0.69480967,4.74119048,0 +22,6,7,0.69480967,5.02439111,0 +22,7,7,0.69480967,7.20185814,1 +22,8,7,0.69480967,7.38133013,1 +22,9,7,0.69480967,7.65934752,1 +22,10,7,0.69480967,7.7971287,1 +22,11,7,0.69480967,7.99285497,1 +22,12,7,0.69480967,8.23263107,1 +22,13,7,0.69480967,8.46983265,1 +22,14,7,0.69480967,8.73286745,1 +23,1,7,1.63741745,1.63039168,0 +23,2,7,1.63741745,2.09752562,0 +23,3,7,1.63741745,2.62818551,0 +23,4,7,1.63741745,3.18544222,0 +23,5,7,1.63741745,3.65712293,0 +23,6,7,1.63741745,4.02960026,0 +23,7,7,1.63741745,6.56398957,1 +23,8,7,1.63741745,7.01914942,1 +23,9,7,1.63741745,7.49708466,1 +23,10,7,1.63741745,7.97035094,1 +23,11,7,1.63741745,8.46558247,1 +23,12,7,1.63741745,8.99876645,1 +23,13,7,1.63741745,9.44318373,1 +23,14,7,1.63741745,9.91803391,1 +24,1,7,0.61461586,-0.1324507,0 +24,2,7,0.61461586,0.02216676,0 +24,3,7,0.61461586,0.19251282,0 +24,4,7,0.61461586,0.40592971,0 +24,5,7,0.61461586,0.64900823,0 +24,6,7,0.61461586,0.7336341,0 +24,7,7,0.61461586,2.92852307,1 +24,8,7,0.61461586,3.05209748,1 +24,9,7,0.61461586,3.24299833,1 +24,10,7,0.61461586,3.45435016,1 +24,11,7,0.61461586,3.57951638,1 +24,12,7,0.61461586,3.76025993,1 +24,13,7,0.61461586,3.97824383,1 +24,14,7,0.61461586,4.12292921,1 +25,1,7,0.17826182,-0.76388521,0 +25,2,7,0.17826182,-0.71777462,0 +25,3,7,0.17826182,-0.63716607,0 +25,4,7,0.17826182,-0.52383659,0 +25,5,7,0.17826182,-0.42204807,0 +25,6,7,0.17826182,-0.33995984,0 +25,7,7,0.17826182,1.75994314,1 +25,8,7,0.17826182,1.8638537,1 +25,9,7,0.17826182,1.90854368,1 +25,10,7,0.17826182,2.07417732,1 +25,11,7,0.17826182,2.20292778,1 +25,12,7,0.17826182,2.25954773,1 +25,13,7,0.17826182,2.34909522,1 +25,14,7,0.17826182,2.46732926,1 +26,1,7,2.6952652,-0.30571824,0 +26,2,7,2.6952652,0.49839644,0 +26,3,7,2.6952652,1.29330915,0 +26,4,7,2.6952652,2.1185548,0 +26,5,7,2.6952652,2.90319557,0 +26,6,7,2.6952652,3.69345488,0 +26,7,7,2.6952652,6.52088365,1 +26,8,7,2.6952652,7.35239264,1 +26,9,7,2.6952652,8.15005858,1 +26,10,7,2.6952652,8.91168421,1 +26,11,7,2.6952652,9.72435336,1 +26,12,7,2.6952652,10.49714276,1 +26,13,7,2.6952652,11.31137528,1 +26,14,7,2.6952652,12.14258372,1 +27,1,7,0.26133846,-0.06473003,0 +27,2,7,0.26133846,-0.03433362,0 +27,3,7,0.26133846,-0.02725789,0 +27,4,7,0.26133846,-0.00073601,0 +27,5,7,0.26133846,0.12916186,0 +27,6,7,0.26133846,0.11938033,0 +27,7,7,0.26133846,2.28909492,1 +27,8,7,0.26133846,2.27229641,1 +27,9,7,0.26133846,2.29552851,1 +27,10,7,0.26133846,2.37233349,1 +27,11,7,0.26133846,2.43198261,1 +27,12,7,0.26133846,2.50383358,1 +27,13,7,0.26133846,2.56178384,1 +27,14,7,0.26133846,2.63910176,1 +28,1,7,2.54221534,1.60224025,0 +28,2,7,2.54221534,2.4354198,0 +28,3,7,2.54221534,3.0965727,0 +28,4,7,2.54221534,3.83062195,0 +28,5,7,2.54221534,4.56478621,0 +28,6,7,2.54221534,5.30581796,0 +28,7,7,2.54221534,8.132134,1 +28,8,7,2.54221534,8.83096189,1 +28,9,7,2.54221534,9.51050418,1 +28,10,7,2.54221534,10.28565743,1 +28,11,7,2.54221534,11.00525689,1 +28,12,7,2.54221534,11.74694847,1 +28,13,7,2.54221534,12.52993607,1 +28,14,7,2.54221534,13.29234606,1 +29,1,7,0.70806507,0.9495518,0 +29,2,7,0.70806507,1.11055892,0 +29,3,7,0.70806507,1.39149501,0 +29,4,7,0.70806507,1.6016901,0 +29,5,7,0.70806507,1.82249104,0 +29,6,7,0.70806507,2.0864023,0 +29,7,7,0.70806507,4.29761179,1 +29,8,7,0.70806507,4.55005105,1 +29,9,7,0.70806507,4.71405014,1 +29,10,7,0.70806507,4.99377233,1 +29,11,7,0.70806507,5.20917,1 +29,12,7,0.70806507,5.4599991,1 +29,13,7,0.70806507,5.68931673,1 +29,14,7,0.70806507,5.90473119,1 +30,1,7,0.89886939,1.69914488,0 +30,2,7,0.89886939,1.9938903,0 +30,3,7,0.89886939,2.24857234,0 +30,4,7,0.89886939,2.49252053,0 +30,5,7,0.89886939,2.76028642,0 +30,6,7,0.89886939,2.93956269,0 +30,7,7,0.89886939,5.26668457,1 +30,8,7,0.89886939,5.51908161,1 +30,9,7,0.89886939,5.75535197,1 +30,10,7,0.89886939,5.99787732,1 +30,11,7,0.89886939,6.2284592,1 +30,12,7,0.89886939,6.55013979,1 +30,13,7,0.89886939,6.77826223,1 +30,14,7,0.89886939,7.07909447,1 +31,1,7,0.07033098,-1.25438788,0 +31,2,7,0.07033098,-1.22800484,0 +31,3,7,0.07033098,-1.23057788,0 +31,4,7,0.07033098,-1.16268669,0 +31,5,7,0.07033098,-1.18692565,0 +31,6,7,0.07033098,-1.18222467,0 +31,7,7,0.07033098,0.81234809,1 +31,8,7,0.07033098,0.79349822,1 +31,9,7,0.07033098,0.86797082,1 +31,10,7,0.07033098,0.76398998,1 +31,11,7,0.07033098,0.81208143,1 +31,12,7,0.07033098,0.7842721,1 +31,13,7,0.07033098,0.8018374,1 +31,14,7,0.07033098,0.78868903,1 +32,1,7,1.68282384,-0.27572013,0 +32,2,7,1.68282384,0.21470232,0 +32,3,7,1.68282384,0.7080657,0 +32,4,7,1.68282384,1.2203077,0 +32,5,7,1.68282384,1.71622125,0 +32,6,7,1.68282384,2.18346623,0 +32,7,7,1.68282384,4.70739965,1 +32,8,7,1.68282384,5.20112356,1 +32,9,7,1.68282384,5.66420854,1 +32,10,7,1.68282384,6.26012466,1 +32,11,7,1.68282384,6.7076865,1 +32,12,7,1.68282384,7.14231833,1 +32,13,7,1.68282384,7.62364944,1 +32,14,7,1.68282384,8.17638787,1 +33,1,7,-0.07140757,-0.69746676,0 +33,2,7,-0.07140757,-0.7791546,0 +33,3,7,-0.07140757,-0.82005093,0 +33,4,7,-0.07140757,-0.85962899,0 +33,5,7,-0.07140757,-0.91489615,0 +33,6,7,-0.07140757,-0.90458685,0 +33,7,7,-0.07140757,1.03228419,1 +33,8,7,-0.07140757,1.00506771,1 +33,9,7,-0.07140757,0.95142468,1 +33,10,7,-0.07140757,0.9807123,1 +33,11,7,-0.07140757,0.86223688,1 +33,12,7,-0.07140757,0.85785161,1 +33,13,7,-0.07140757,0.79776473,1 +33,14,7,-0.07140757,0.71440504,1 +34,1,7,1.18092228,0.78168934,0 +34,2,7,1.18092228,1.24323852,0 +34,3,7,1.18092228,1.52502097,0 +34,4,7,1.18092228,1.91740774,0 +34,5,7,1.18092228,2.26991905,0 +34,6,7,1.18092228,2.67302863,0 +34,7,7,1.18092228,5.12666421,1 +34,8,7,1.18092228,5.41869442,1 +34,9,7,1.18092228,5.77454804,1 +34,10,7,1.18092228,6.15045723,1 +34,11,7,1.18092228,6.52645234,1 +34,12,7,1.18092228,6.89445317,1 +34,13,7,1.18092228,7.26531642,1 +34,14,7,1.18092228,7.6904745,1 +35,1,7,2.37121725,0.67516056,0 +35,2,7,2.37121725,1.39279659,0 +35,3,7,2.37121725,2.13662046,0 +35,4,7,2.37121725,2.82963196,0 +35,5,7,2.37121725,3.54041975,0 +35,6,7,2.37121725,4.26268816,0 +35,7,7,2.37121725,6.96495799,1 +35,8,7,2.37121725,7.67703742,1 +35,9,7,2.37121725,8.38665343,1 +35,10,7,2.37121725,9.12493824,1 +35,11,7,2.37121725,9.78784686,1 +35,12,7,2.37121725,10.54408408,1 +35,13,7,2.37121725,11.24733989,1 +35,14,7,2.37121725,11.93841693,1 +36,1,7,1.09868075,1.01313034,0 +36,2,7,1.09868075,1.34812128,0 +36,3,7,1.09868075,1.56634909,0 +36,4,7,1.09868075,1.97941854,0 +36,5,7,1.09868075,2.23032436,0 +36,6,7,1.09868075,2.57208432,0 +36,7,7,1.09868075,4.83209779,1 +36,8,7,1.09868075,5.11159147,1 +36,9,7,1.09868075,5.42873436,1 +36,10,7,1.09868075,5.80525514,1 +36,11,7,1.09868075,6.06222693,1 +36,12,7,1.09868075,6.3486017,1 +36,13,7,1.09868075,6.7356342,1 +36,14,7,1.09868075,6.95100273,1 +37,1,9,-0.02335637,-2.30314919,0 +37,2,9,-0.02335637,-2.28958647,0 +37,3,9,-0.02335637,-2.21449073,0 +37,4,9,-0.02335637,-2.23960048,0 +37,5,9,-0.02335637,-2.23494718,0 +37,6,9,-0.02335637,-2.18984267,0 +37,7,9,-0.02335637,-2.18365708,0 +37,8,9,-0.02335637,-2.1625995,0 +37,9,9,-0.02335637,-0.15088298,1 +37,10,9,-0.02335637,-0.09798289,1 +37,11,9,-0.02335637,-0.10932186,1 +37,12,9,-0.02335637,-0.12468418,1 +37,13,9,-0.02335637,-0.08286975,1 +37,14,9,-0.02335637,-0.07749066,1 +38,1,9,0.57591315,0.24134147,0 +38,2,9,0.57591315,0.42238335,0 +38,3,9,0.57591315,0.67358584,0 +38,4,9,0.57591315,0.89328705,0 +38,5,9,0.57591315,1.07903282,0 +38,6,9,0.57591315,1.26187869,0 +38,7,9,0.57591315,1.35771596,0 +38,8,9,0.57591315,1.57771265,0 +38,9,9,0.57591315,3.77130402,1 +38,10,9,0.57591315,3.99405597,1 +38,11,9,0.57591315,4.13782155,1 +38,12,9,0.57591315,4.34320594,1 +38,13,9,0.57591315,4.5244316,1 +38,14,9,0.57591315,4.73434146,1 +39,1,9,2.04803794,0.95354806,0 +39,2,9,2.04803794,1.63845368,0 +39,3,9,2.04803794,2.23417888,0 +39,4,9,2.04803794,2.87880593,0 +39,5,9,2.04803794,3.48001107,0 +39,6,9,2.04803794,4.17691971,0 +39,7,9,2.04803794,4.80393383,0 +39,8,9,2.04803794,5.42986331,0 +39,9,9,2.04803794,8.04206355,1 +39,10,9,2.04803794,8.67206005,1 +39,11,9,2.04803794,9.36766648,1 +39,12,9,2.04803794,10.01512265,1 +39,13,9,2.04803794,10.65016029,1 +39,14,9,2.04803794,11.33940438,1 +40,1,9,0.27439432,-0.84243157,0 +40,2,9,0.27439432,-0.78031295,0 +40,3,9,0.27439432,-0.78528402,0 +40,4,9,0.27439432,-0.64609682,0 +40,5,9,0.27439432,-0.67463095,0 +40,6,9,0.27439432,-0.59209698,0 +40,7,9,0.27439432,-0.60297847,0 +40,8,9,0.27439432,-0.52061873,0 +40,9,9,0.27439432,1.47819515,1 +40,10,9,0.27439432,1.52152392,1 +40,11,9,0.27439432,1.58354591,1 +40,12,9,0.27439432,1.61019761,1 +40,13,9,0.27439432,1.64933016,1 +40,14,9,0.27439432,1.6593388,1 +41,1,9,1.33278514,0.06124999,0 +41,2,9,1.33278514,0.49514705,0 +41,3,9,1.33278514,0.90646541,0 +41,4,9,1.33278514,1.28547681,0 +41,5,9,1.33278514,1.64366273,0 +41,6,9,1.33278514,2.04350652,0 +41,7,9,1.33278514,2.46015371,0 +41,8,9,1.33278514,2.83456642,0 +41,9,9,1.33278514,5.28466523,1 +41,10,9,1.33278514,5.65580079,1 +41,11,9,1.33278514,6.05562003,1 +41,12,9,1.33278514,6.42811755,1 +41,13,9,1.33278514,6.7479916,1 +41,14,9,1.33278514,7.2063016,1 +42,1,9,0.93772643,1.06794723,0 +42,2,9,0.93772643,1.38324858,0 +42,3,9,0.93772643,1.75192118,0 +42,4,9,0.93772643,2.11097241,0 +42,5,9,0.93772643,2.41849495,0 +42,6,9,0.93772643,2.70987295,0 +42,7,9,0.93772643,3.04849956,0 +42,8,9,0.93772643,3.35202373,0 +42,9,9,0.93772643,5.68275127,1 +42,10,9,0.93772643,5.98726103,1 +42,11,9,0.93772643,6.32756382,1 +42,12,9,0.93772643,6.6873886,1 +42,13,9,0.93772643,7.0071445,1 +42,14,9,0.93772643,7.30461842,1 +43,1,9,0.29287216,-0.00900984,0 +43,2,9,0.29287216,0.07513345,0 +43,3,9,0.29287216,0.10010729,0 +43,4,9,0.29287216,0.17531445,0 +43,5,9,0.29287216,0.2050564,0 +43,6,9,0.29287216,0.30446309,0 +43,7,9,0.29287216,0.35077423,0 +43,8,9,0.29287216,0.42674868,0 +43,9,9,0.29287216,2.47378165,1 +43,10,9,0.29287216,2.53794965,1 +43,11,9,0.29287216,2.60048783,1 +43,12,9,0.29287216,2.66668126,1 +43,13,9,0.29287216,2.67844988,1 +43,14,9,0.29287216,2.78048552,1 +44,1,9,-0.48356888,-0.8833158,0 +44,2,9,-0.48356888,-1.02187099,0 +44,3,9,-0.48356888,-1.13736524,0 +44,4,9,-0.48356888,-1.32908451,0 +44,5,9,-0.48356888,-1.40633121,0 +44,6,9,-0.48356888,-1.49660096,0 +44,7,9,-0.48356888,-1.58962295,0 +44,8,9,-0.48356888,-1.7678052,0 +44,9,9,-0.48356888,0.13692757,1 +44,10,9,-0.48356888,-0.0640963,1 +44,11,9,-0.48356888,-0.08598583,1 +44,12,9,-0.48356888,-0.29642704,1 +44,13,9,-0.48356888,-0.40815878,1 +44,14,9,-0.48356888,-0.54350954,1 +45,1,9,0.90279289,2.01686643,0 +45,2,9,0.90279289,2.31279676,0 +45,3,9,0.90279289,2.56718212,0 +45,4,9,0.90279289,2.91476725,0 +45,5,9,0.90279289,3.1361303,0 +45,6,9,0.90279289,3.45760828,0 +45,7,9,0.90279289,3.74284828,0 +45,8,9,0.90279289,3.91504431,0 +45,9,9,0.90279289,6.24725799,1 +45,10,9,0.90279289,6.49599316,1 +45,11,9,0.90279289,6.7539757,1 +45,12,9,0.90279289,7.04688353,1 +45,13,9,0.90279289,7.27467229,1 +45,14,9,0.90279289,7.50534982,1 +46,1,9,1.36512074,0.54290496,0 +46,2,9,1.36512074,1.02205292,0 +46,3,9,1.36512074,1.46944665,0 +46,4,9,1.36512074,1.79001862,0 +46,5,9,1.36512074,2.23770915,0 +46,6,9,1.36512074,2.71377247,0 +46,7,9,1.36512074,3.12917286,0 +46,8,9,1.36512074,3.55304259,0 +46,9,9,1.36512074,5.99348265,1 +46,10,9,1.36512074,6.44532473,1 +46,11,9,1.36512074,6.82572389,1 +46,12,9,1.36512074,7.26966028,1 +46,13,9,1.36512074,7.74293036,1 +46,14,9,1.36512074,8.09985678,1 +47,1,9,0.81599385,-0.38623635,0 +47,2,9,0.81599385,-0.12001672,0 +47,3,9,0.81599385,0.01560652,0 +47,4,9,0.81599385,0.30032658,0 +47,5,9,0.81599385,0.52837487,0 +47,6,9,0.81599385,0.70685538,0 +47,7,9,0.81599385,0.9643763,0 +47,8,9,0.81599385,1.15148682,0 +47,9,9,0.81599385,3.43752336,1 +47,10,9,0.81599385,3.63460113,1 +47,11,9,0.81599385,3.80328712,1 +47,12,9,0.81599385,4.05866868,1 +47,13,9,0.81599385,4.28232068,1 +47,14,9,0.81599385,4.53130627,1 +48,1,9,0.69987487,0.35070934,0 +48,2,9,0.69987487,0.58010672,0 +48,3,9,0.69987487,0.79978655,0 +48,4,9,0.69987487,0.97092667,0 +48,5,9,0.69987487,1.24307203,0 +48,6,9,0.69987487,1.46388243,0 +48,7,9,0.69987487,1.64029418,0 +48,8,9,0.69987487,1.82644605,0 +48,9,9,0.69987487,4.06070566,1 +48,10,9,0.69987487,4.31874436,1 +48,11,9,0.69987487,4.49514478,1 +48,12,9,0.69987487,4.74405908,1 +48,13,9,0.69987487,4.93242636,1 +48,14,9,0.69987487,5.1750701,1 +49,1,9,-1.14509609,0.72279287,0 +49,2,9,-1.14509609,0.39139985,0 +49,3,9,-1.14509609,0.02764879,0 +49,4,9,-1.14509609,-0.36261003,0 +49,5,9,-1.14509609,-0.85272459,0 +49,6,9,-1.14509609,-1.15608097,0 +49,7,9,-1.14509609,-1.53666477,0 +49,8,9,-1.14509609,-1.89827943,0 +49,9,9,-1.14509609,-0.27033541,1 +49,10,9,-1.14509609,-0.69843953,1 +49,11,9,-1.14509609,-1.07711703,1 +49,12,9,-1.14509609,-1.45271676,1 +49,13,9,-1.14509609,-1.85232125,1 +49,14,9,-1.14509609,-2.2050032,1 +50,1,9,0.68508563,-1.18322663,0 +50,2,9,0.68508563,-1.05406361,0 +50,3,9,0.68508563,-0.89426441,0 +50,4,9,0.68508563,-0.69559207,0 +50,5,9,0.68508563,-0.53378429,0 +50,6,9,0.68508563,-0.39655585,0 +50,7,9,0.68508563,-0.20238909,0 +50,8,9,0.68508563,-0.00891307,0 +50,9,9,0.68508563,2.14020593,1 +50,10,9,0.68508563,2.2923739,1 +50,11,9,0.68508563,2.54860175,1 +50,12,9,0.68508563,2.64890539,1 +50,13,9,0.68508563,2.8241236,1 +50,14,9,0.68508563,2.99258688,1 +51,1,9,2.26778856,1.17153315,0 +51,2,9,2.26778856,1.86157334,0 +51,3,9,2.26778856,2.5985609,0 +51,4,9,2.26778856,3.31196215,0 +51,5,9,2.26778856,3.98024026,0 +51,6,9,2.26778856,4.66264972,0 +51,7,9,2.26778856,5.47126058,0 +51,8,9,2.26778856,6.17232399,0 +51,9,9,2.26778856,8.87608126,1 +51,10,9,2.26778856,9.54857631,1 +51,11,9,2.26778856,10.26172321,1 +51,12,9,2.26778856,10.96753546,1 +51,13,9,2.26778856,11.7259937,1 +51,14,9,2.26778856,12.46842801,1 +52,1,9,1.08584385,-0.58112316,0 +52,2,9,1.08584385,-0.3139937,0 +52,3,9,1.08584385,0.03254066,0 +52,4,9,1.08584385,0.31567298,0 +52,5,9,1.08584385,0.58101912,0 +52,6,9,1.08584385,0.97414089,0 +52,7,9,1.08584385,1.23142814,0 +52,8,9,1.08584385,1.5189514,0 +52,9,9,1.08584385,3.82034457,1 +52,10,9,1.08584385,4.19181682,1 +52,11,9,1.08584385,4.46867313,1 +52,12,9,1.08584385,4.77398021,1 +52,13,9,1.08584385,5.15864682,1 +52,14,9,1.08584385,5.35333218,1 +53,1,9,1.21890949,-2.55247163,0 +53,2,9,1.21890949,-2.26112615,0 +53,3,9,1.21890949,-1.85584795,0 +53,4,9,1.21890949,-1.41730445,0 +53,5,9,1.21890949,-1.05417326,0 +53,6,9,1.21890949,-0.61628982,0 +53,7,9,1.21890949,-0.23449127,0 +53,8,9,1.21890949,0.19427479,0 +53,9,9,1.21890949,2.57600124,1 +53,10,9,1.21890949,2.98534451,1 +53,11,9,1.21890949,3.32542224,1 +53,12,9,1.21890949,3.7846396,1 +53,13,9,1.21890949,4.19183213,1 +53,14,9,1.21890949,4.58826869,1 +54,1,9,-0.16317057,2.12356684,0 +54,2,9,-0.16317057,2.03072472,0 +54,3,9,-0.16317057,2.02355491,0 +54,4,9,-0.16317057,1.94948253,0 +54,5,9,-0.16317057,1.89889557,0 +54,6,9,-0.16317057,1.8556241,0 +54,7,9,-0.16317057,1.7577265,0 +54,8,9,-0.16317057,1.74622285,0 +54,9,9,-0.16317057,3.69537021,1 +54,10,9,-0.16317057,3.63229466,1 +54,11,9,-0.16317057,3.57273417,1 +54,12,9,-0.16317057,3.50381799,1 +54,13,9,-0.16317057,3.43874407,1 +54,14,9,-0.16317057,3.40159688,1 +55,1,0,1.58054535,0.90156353,0 +55,2,0,1.58054535,1.38568412,0 +55,3,0,1.58054535,1.8555898,0 +55,4,0,1.58054535,2.31677611,0 +55,5,0,1.58054535,2.83323135,0 +55,6,0,1.58054535,3.26583592,0 +55,7,0,1.58054535,3.79785286,0 +55,8,0,1.58054535,4.26080618,0 +55,9,0,1.58054535,4.73157597,0 +55,10,0,1.58054535,5.17521734,0 +55,11,0,1.58054535,5.68051611,0 +55,12,0,1.58054535,6.19862717,0 +55,13,0,1.58054535,6.65732067,0 +55,14,0,1.58054535,7.12862139,0 +56,1,0,-0.21703623,0.25606051,0 +56,2,0,-0.21703623,0.15303417,0 +56,3,0,-0.21703623,0.14219721,0 +56,4,0,-0.21703623,0.12227549,0 +56,5,0,-0.21703623,0.09231433,0 +56,6,0,-0.21703623,0.07298236,0 +56,7,0,-0.21703623,0.00860778,0 +56,8,0,-0.21703623,-0.03358761,0 +56,9,0,-0.21703623,-0.02852995,0 +56,10,0,-0.21703623,-0.03860497,0 +56,11,0,-0.21703623,-0.06811465,0 +56,12,0,-0.21703623,-0.09640895,0 +56,13,0,-0.21703623,-0.1669684,0 +56,14,0,-0.21703623,-0.15253021,0 +57,1,0,-2.37697981,0.42032095,0 +57,2,0,-2.37697981,-0.24405915,0 +57,3,0,-2.37697981,-0.99397751,0 +57,4,0,-2.37697981,-1.69567313,0 +57,5,0,-2.37697981,-2.42827878,0 +57,6,0,-2.37697981,-3.16011489,0 +57,7,0,-2.37697981,-3.90856306,0 +57,8,0,-2.37697981,-4.60862923,0 +57,9,0,-2.37697981,-5.27776029,0 +57,10,0,-2.37697981,-6.02046028,0 +57,11,0,-2.37697981,-6.72360688,0 +57,12,0,-2.37697981,-7.48253578,0 +57,13,0,-2.37697981,-8.17322054,0 +57,14,0,-2.37697981,-8.90000979,0 +58,1,0,-0.35184034,-0.71169753,0 +58,2,0,-0.35184034,-0.74483188,0 +58,3,0,-0.35184034,-0.89226584,0 +58,4,0,-0.35184034,-0.93062579,0 +58,5,0,-0.35184034,-1.013612,0 +58,6,0,-0.35184034,-1.13106368,0 +58,7,0,-0.35184034,-1.20369731,0 +58,8,0,-0.35184034,-1.26111968,0 +58,9,0,-0.35184034,-1.34238375,0 +58,10,0,-0.35184034,-1.39812319,0 +58,11,0,-0.35184034,-1.49123833,0 +58,12,0,-0.35184034,-1.62610318,0 +58,13,0,-0.35184034,-1.64802658,0 +58,14,0,-0.35184034,-1.74664289,0 +59,1,0,-0.72388845,-0.25692459,0 +59,2,0,-0.72388845,-0.42909954,0 +59,3,0,-0.72388845,-0.64228268,0 +59,4,0,-0.72388845,-0.8668892,0 +59,5,0,-0.72388845,-1.09399806,0 +59,6,0,-0.72388845,-1.26516184,0 +59,7,0,-0.72388845,-1.47082585,0 +59,8,0,-0.72388845,-1.68661123,0 +59,9,0,-0.72388845,-1.91228068,0 +59,10,0,-0.72388845,-2.13629991,0 +59,11,0,-0.72388845,-2.33667407,0 +59,12,0,-0.72388845,-2.46256099,0 +59,13,0,-0.72388845,-2.79844738,0 +59,14,0,-0.72388845,-2.90260278,0 +60,1,0,-0.48363748,-0.29845815,0 +60,2,0,-0.48363748,-0.50089599,0 +60,3,0,-0.48363748,-0.61305137,0 +60,4,0,-0.48363748,-0.79432557,0 +60,5,0,-0.48363748,-0.9706834,0 +60,6,0,-0.48363748,-1.15706078,0 +60,7,0,-0.48363748,-1.3390131,0 +60,8,0,-0.48363748,-1.46702496,0 +60,9,0,-0.48363748,-1.6841186,0 +60,10,0,-0.48363748,-1.79442394,0 +60,11,0,-0.48363748,-1.98872829,0 +60,12,0,-0.48363748,-2.21268868,0 +60,13,0,-0.48363748,-2.34568845,0 +60,14,0,-0.48363748,-2.49690867,0 +61,1,0,-0.42738475,0.15224757,0 +61,2,0,-0.42738475,-0.01504092,0 +61,3,0,-0.42738475,-0.14713962,0 +61,4,0,-0.42738475,-0.19431336,0 +61,5,0,-0.42738475,-0.39565808,0 +61,6,0,-0.42738475,-0.52270025,0 +61,7,0,-0.42738475,-0.64912181,0 +61,8,0,-0.42738475,-0.75385476,0 +61,9,0,-0.42738475,-0.90572636,0 +61,10,0,-0.42738475,-1.01281817,0 +61,11,0,-0.42738475,-1.14318888,0 +61,12,0,-0.42738475,-1.27542952,0 +61,13,0,-0.42738475,-1.42896872,0 +61,14,0,-0.42738475,-1.54822176,0 +62,1,0,0.30562977,-1.10915863,0 +62,2,0,0.30562977,-1.05854258,0 +62,3,0,0.30562977,-0.89722768,0 +62,4,0,0.30562977,-0.81499401,0 +62,5,0,0.30562977,-0.74007191,0 +62,6,0,0.30562977,-0.64370385,0 +62,7,0,0.30562977,-0.52917124,0 +62,8,0,0.30562977,-0.45149234,0 +62,9,0,0.30562977,-0.32346851,0 +62,10,0,0.30562977,-0.18826982,0 +62,11,0,0.30562977,-0.20131184,0 +62,12,0,0.30562977,-0.00991409,0 +62,13,0,0.30562977,0.09491213,0 +62,14,0,0.30562977,0.21390497,0 +63,1,0,1.18315075,1.51706677,0 +63,2,0,1.18315075,1.90864142,0 +63,3,0,1.18315075,2.33034607,0 +63,4,0,1.18315075,2.61108269,0 +63,5,0,1.18315075,3.06231714,0 +63,6,0,1.18315075,3.36504135,0 +63,7,0,1.18315075,3.72735239,0 +63,8,0,1.18315075,4.00465805,0 +63,9,0,1.18315075,4.49716339,0 +63,10,0,1.18315075,4.84042599,0 +63,11,0,1.18315075,5.13871637,0 +63,12,0,1.18315075,5.52228039,0 +63,13,0,1.18315075,5.86254323,0 +63,14,0,1.18315075,6.23537076,0 +64,1,0,0.69269082,-0.21621487,0 +64,2,0,0.69269082,-0.10941651,0 +64,3,0,0.69269082,0.1285544,0 +64,4,0,0.69269082,0.31962655,0 +64,5,0,0.69269082,0.42923612,0 +64,6,0,0.69269082,0.68284635,0 +64,7,0,0.69269082,0.85627462,0 +64,8,0,0.69269082,1.01844415,0 +64,9,0,0.69269082,1.18345361,0 +64,10,0,0.69269082,1.3736362,0 +64,11,0,0.69269082,1.53037365,0 +64,12,0,0.69269082,1.66932924,0 +64,13,0,0.69269082,1.90148165,0 +64,14,0,0.69269082,2.03497568,0 +65,1,0,0.15346359,1.32518599,0 +65,2,0,0.15346359,1.42449867,0 +65,3,0,0.15346359,1.47474371,0 +65,4,0,0.15346359,1.55352904,0 +65,5,0,0.15346359,1.69424676,0 +65,6,0,0.15346359,1.73709368,0 +65,7,0,0.15346359,1.80400389,0 +65,8,0,0.15346359,1.91393132,0 +65,9,0,0.15346359,1.97891153,0 +65,10,0,0.15346359,2.03088831,0 +65,11,0,0.15346359,2.09701578,0 +65,12,0,0.15346359,2.20216223,0 +65,13,0,0.15346359,2.3470251,0 +65,14,0,0.15346359,2.36133988,0 +66,1,0,1.43242481,-0.7552191,0 +66,2,0,1.43242481,-0.3054247,0 +66,3,0,1.43242481,0.11880881,0 +66,4,0,1.43242481,0.53360177,0 +66,5,0,1.43242481,0.99395459,0 +66,6,0,1.43242481,1.4502063,0 +66,7,0,1.43242481,1.89923881,0 +66,8,0,1.43242481,2.35404999,0 +66,9,0,1.43242481,2.66343312,0 +66,10,0,1.43242481,3.1312325,0 +66,11,0,1.43242481,3.53535049,0 +66,12,0,1.43242481,3.91489722,0 +66,13,0,1.43242481,4.40471446,0 +66,14,0,1.43242481,4.79422176,0 +67,1,0,-2.34298948,-3.77867871,0 +67,2,0,-2.34298948,-4.39264927,0 +67,3,0,-2.34298948,-5.13539899,0 +67,4,0,-2.34298948,-5.89163276,0 +67,5,0,-2.34298948,-6.52655033,0 +67,6,0,-2.34298948,-7.32314508,0 +67,7,0,-2.34298948,-8.04151879,0 +67,8,0,-2.34298948,-8.75232539,0 +67,9,0,-2.34298948,-9.54996455,0 +67,10,0,-2.34298948,-10.19723407,0 +67,11,0,-2.34298948,-10.95437425,0 +67,12,0,-2.34298948,-11.68934375,0 +67,13,0,-2.34298948,-12.40885584,0 +67,14,0,-2.34298948,-13.04754267,0 +68,1,0,-1.50171659,-2.27421908,0 +68,2,0,-1.50171659,-2.73714776,0 +68,3,0,-1.50171659,-3.22685606,0 +68,4,0,-1.50171659,-3.72853747,0 +68,5,0,-1.50171659,-4.21775976,0 +68,6,0,-1.50171659,-4.6886591,0 +68,7,0,-1.50171659,-5.06446408,0 +68,8,0,-1.50171659,-5.56843192,0 +68,9,0,-1.50171659,-6.09476345,0 +68,10,0,-1.50171659,-6.57447951,0 +68,11,0,-1.50171659,-6.9932541,0 +68,12,0,-1.50171659,-7.5054013,0 +68,13,0,-1.50171659,-7.95657258,0 +68,14,0,-1.50171659,-8.44424555,0 +69,1,0,-0.38916756,1.42210452,0 +69,2,0,-0.38916756,1.28170529,0 +69,3,0,-0.38916756,1.14642473,0 +69,4,0,-0.38916756,0.99245409,0 +69,5,0,-0.38916756,0.85078438,0 +69,6,0,-0.38916756,0.72747438,0 +69,7,0,-0.38916756,0.57209241,0 +69,8,0,-0.38916756,0.47968538,0 +69,9,0,-0.38916756,0.34197325,0 +69,10,0,-0.38916756,0.18376222,0 +69,11,0,-0.38916756,0.04599472,0 +69,12,0,-0.38916756,-0.05058671,0 +69,13,0,-0.38916756,-0.22180723,0 +69,14,0,-0.38916756,-0.33525017,0 +70,1,0,0.94795794,-0.29181369,0 +70,2,0,0.94795794,-0.01326701,0 +70,3,0,0.94795794,0.25261343,0 +70,4,0,0.94795794,0.45369005,0 +70,5,0,0.94795794,0.73897714,0 +70,6,0,0.94795794,0.96019885,0 +70,7,0,0.94795794,1.26762747,0 +70,8,0,0.94795794,1.48632499,0 +70,9,0,0.94795794,1.69413587,0 +70,10,0,0.94795794,2.00270345,0 +70,11,0,0.94795794,2.19214852,0 +70,12,0,0.94795794,2.51112518,0 +70,13,0,0.94795794,2.74898078,0 +70,14,0,0.94795794,2.99164804,0 +71,1,0,0.06958891,-0.16403003,0 +71,2,0,0.06958891,-0.1668906,0 +71,3,0,0.06958891,-0.12395097,0 +71,4,0,0.06958891,-0.08910075,0 +71,5,0,0.06958891,-0.03277997,0 +71,6,0,0.06958891,-0.02121622,0 +71,7,0,0.06958891,-0.07579358,0 +71,8,0,0.06958891,0.04755849,0 +71,9,0,0.06958891,0.04584969,0 +71,10,0,0.06958891,0.03365915,0 +71,11,0,0.06958891,0.08572376,0 +71,12,0,0.06958891,0.14084486,0 +71,13,0,0.06958891,0.15325613,0 +71,14,0,0.06958891,0.13953346,0 +72,1,0,-0.68703783,1.4457824,0 +72,2,0,-0.68703783,1.20639313,0 +72,3,0,-0.68703783,1.03542087,0 +72,4,0,-0.68703783,0.82570815,0 +72,5,0,-0.68703783,0.56418543,0 +72,6,0,-0.68703783,0.38916052,0 +72,7,0,-0.68703783,0.14461866,0 +72,8,0,-0.68703783,-0.04562296,0 +72,9,0,-0.68703783,-0.23559265,0 +72,10,0,-0.68703783,-0.51116879,0 +72,11,0,-0.68703783,-0.69497472,0 +72,12,0,-0.68703783,-0.84880285,0 +72,13,0,-0.68703783,-1.01257675,0 +72,14,0,-0.68703783,-1.23747264,0 +73,1,0,-0.73634771,0.26303912,0 +73,2,0,-0.73634771,0.00899379,0 +73,3,0,-0.73634771,-0.20084727,0 +73,4,0,-0.73634771,-0.39289394,0 +73,5,0,-0.73634771,-0.64379523,0 +73,6,0,-0.73634771,-0.91657284,0 +73,7,0,-0.73634771,-1.0140854,0 +73,8,0,-0.73634771,-1.25920465,0 +73,9,0,-0.73634771,-1.51249711,0 +73,10,0,-0.73634771,-1.74953414,0 +73,11,0,-0.73634771,-1.97923743,0 +73,12,0,-0.73634771,-2.12153538,0 +73,13,0,-0.73634771,-2.39631277,0 +73,14,0,-0.73634771,-2.59108915,0 +74,1,0,0.46092395,-0.30814815,0 +74,2,0,0.46092395,-0.18300324,0 +74,3,0,0.46092395,0.01808658,0 +74,4,0,0.46092395,0.13606293,0 +74,5,0,0.46092395,0.36835717,0 +74,6,0,0.46092395,0.49773705,0 +74,7,0,0.46092395,0.6397167,0 +74,8,0,0.46092395,0.77762642,0 +74,9,0,0.46092395,0.89409556,0 +74,10,0,0.46092395,1.02870297,0 +74,11,0,0.46092395,1.20855823,0 +74,12,0,0.46092395,1.33119345,0 +74,13,0,0.46092395,1.5256443,0 +74,14,0,0.46092395,1.73792452,0 +75,1,0,0.75815412,1.72588827,0 +75,2,0,0.75815412,1.9783213,0 +75,3,0,0.75815412,2.18574581,0 +75,4,0,0.75815412,2.43847456,0 +75,5,0,0.75815412,2.65286803,0 +75,6,0,0.75815412,2.85493662,0 +75,7,0,0.75815412,3.11310069,0 +75,8,0,0.75815412,3.31129808,0 +75,9,0,0.75815412,3.60078527,0 +75,10,0,0.75815412,3.800725,0 +75,11,0,0.75815412,4.03286584,0 +75,12,0,0.75815412,4.22495143,0 +75,13,0,0.75815412,4.42672001,0 +75,14,0,0.75815412,4.67592213,0 +76,1,0,-2.18213482,1.6985665,0 +76,2,0,-2.18213482,1.12181202,0 +76,3,0,-2.18213482,0.48887991,0 +76,4,0,-2.18213482,-0.18475336,0 +76,5,0,-2.18213482,-0.85355981,0 +76,6,0,-2.18213482,-1.35885966,0 +76,7,0,-2.18213482,-2.05642252,0 +76,8,0,-2.18213482,-2.59954108,0 +76,9,0,-2.18213482,-3.31007255,0 +76,10,0,-2.18213482,-3.87472255,0 +76,11,0,-2.18213482,-4.50991183,0 +76,12,0,-2.18213482,-5.11809168,0 +76,13,0,-2.18213482,-5.75786437,0 +76,14,0,-2.18213482,-6.38569115,0 +77,1,0,-0.24164157,-1.36182126,0 +77,2,0,-0.24164157,-1.45378131,0 +77,3,0,-0.24164157,-1.50835269,0 +77,4,0,-0.24164157,-1.62872003,0 +77,5,0,-0.24164157,-1.72527612,0 +77,6,0,-0.24164157,-1.8140916,0 +77,7,0,-0.24164157,-1.92752541,0 +77,8,0,-0.24164157,-1.96054708,0 +77,9,0,-0.24164157,-2.03753448,0 +77,10,0,-0.24164157,-2.09572254,0 +77,11,0,-0.24164157,-2.20311603,0 +77,12,0,-0.24164157,-2.35612658,0 +77,13,0,-0.24164157,-2.36459314,0 +77,14,0,-0.24164157,-2.52493452,0 +78,1,0,0.43575038,-0.31659758,0 +78,2,0,0.43575038,-0.10465918,0 +78,3,0,0.43575038,0.03316598,0 +78,4,0,0.43575038,0.20954722,0 +78,5,0,0.43575038,0.278118,0 +78,6,0,0.43575038,0.50146619,0 +78,7,0,0.43575038,0.66718639,0 +78,8,0,0.43575038,0.78986818,0 +78,9,0,0.43575038,0.96962,0 +78,10,0,0.43575038,1.20875249,0 +78,11,0,0.43575038,1.27819446,0 +78,12,0,0.43575038,1.49165589,0 +78,13,0,0.43575038,1.62648467,0 +78,14,0,0.43575038,1.79571772,0 +79,1,0,0.88890768,1.95573442,0 +79,2,0,0.88890768,2.17538437,0 +79,3,0,0.88890768,2.46867024,0 +79,4,0,0.88890768,2.71325283,0 +79,5,0,0.88890768,3.03666697,0 +79,6,0,0.88890768,3.29002935,0 +79,7,0,0.88890768,3.51261382,0 +79,8,0,0.88890768,3.78422261,0 +79,9,0,0.88890768,4.02386271,0 +79,10,0,0.88890768,4.31675767,0 +79,11,0,0.88890768,4.57866794,0 +79,12,0,0.88890768,4.84820163,0 +79,13,0,0.88890768,5.12300133,0 +79,14,0,0.88890768,5.34864138,0 +80,1,0,0.61236265,-0.83223841,0 +80,2,0,0.61236265,-0.70884572,0 +80,3,0,0.61236265,-0.4286798,0 +80,4,0,0.61236265,-0.1798602,0 +80,5,0,0.61236265,-0.06173553,0 +80,6,0,0.61236265,0.20174259,0 +80,7,0,0.61236265,0.35944202,0 +80,8,0,0.61236265,0.55905141,0 +80,9,0,0.61236265,0.80160015,0 +80,10,0,0.61236265,0.94302646,0 +80,11,0,0.61236265,1.18045401,0 +80,12,0,0.61236265,1.41001957,0 +80,13,0,0.61236265,1.65835701,0 +80,14,0,0.61236265,1.82542175,0 +81,1,0,-0.84653551,-0.30186321,0 +81,2,0,-0.84653551,-0.49228253,0 +81,3,0,-0.84653551,-0.7362166,0 +81,4,0,-0.84653551,-1.04009901,0 +81,5,0,-0.84653551,-1.22679441,0 +81,6,0,-0.84653551,-1.52961285,0 +81,7,0,-0.84653551,-1.78900594,0 +81,8,0,-0.84653551,-2.04500873,0 +81,9,0,-0.84653551,-2.26425151,0 +81,10,0,-0.84653551,-2.53902769,0 +81,11,0,-0.84653551,-2.7960582,0 +81,12,0,-0.84653551,-3.01766862,0 +81,13,0,-0.84653551,-3.31918104,0 +81,14,0,-0.84653551,-3.55123605,0 +82,1,0,-0.30962931,-1.36039469,0 +82,2,0,-0.30962931,-1.48011854,0 +82,3,0,-0.30962931,-1.60714468,0 +82,4,0,-0.30962931,-1.57205277,0 +82,5,0,-0.30962931,-1.70405212,0 +82,6,0,-0.30962931,-1.78911423,0 +82,7,0,-0.30962931,-1.94755467,0 +82,8,0,-0.30962931,-1.99790364,0 +82,9,0,-0.30962931,-2.0779828,0 +82,10,0,-0.30962931,-2.18794295,0 +82,11,0,-0.30962931,-2.30987006,0 +82,12,0,-0.30962931,-2.35676427,0 +82,13,0,-0.30962931,-2.42624622,0 +82,14,0,-0.30962931,-2.56800011,0 +83,1,0,2.01760342,-0.1542964,0 +83,2,0,2.01760342,0.44198177,0 +83,3,0,2.01760342,1.14560801,0 +83,4,0,2.01760342,1.74549111,0 +83,5,0,2.01760342,2.32459169,0 +83,6,0,2.01760342,3.02908144,0 +83,7,0,2.01760342,3.59514071,0 +83,8,0,2.01760342,4.18083812,0 +83,9,0,2.01760342,4.84188397,0 +83,10,0,2.01760342,5.33456232,0 +83,11,0,2.01760342,6.03712135,0 +83,12,0,2.01760342,6.67078861,0 +83,13,0,2.01760342,7.3370149,0 +83,14,0,2.01760342,7.81081035,0 +84,1,0,-0.04102571,0.14979333,0 +84,2,0,-0.04102571,0.06173188,0 +84,3,0,-0.04102571,0.06618106,0 +84,4,0,-0.04102571,0.09309177,0 +84,5,0,-0.04102571,0.07965821,0 +84,6,0,-0.04102571,0.09629449,0 +84,7,0,-0.04102571,0.03062446,0 +84,8,0,-0.04102571,0.0326756,0 +84,9,0,-0.04102571,-0.01277945,0 +84,10,0,-0.04102571,-0.08252767,0 +84,11,0,-0.04102571,-0.05875817,0 +84,12,0,-0.04102571,-0.07820592,0 +84,13,0,-0.04102571,-0.10514677,0 +84,14,0,-0.04102571,-0.09049598,0 +85,1,0,1.51886789,1.39781771,0 +85,2,0,1.51886789,1.85890013,0 +85,3,0,1.51886789,2.32084782,0 +85,4,0,1.51886789,2.80717466,0 +85,5,0,1.51886789,3.34330134,0 +85,6,0,1.51886789,3.76842679,0 +85,7,0,1.51886789,4.32313485,0 +85,8,0,1.51886789,4.71802677,0 +85,9,0,1.51886789,5.18387047,0 +85,10,0,1.51886789,5.66059034,0 +85,11,0,1.51886789,6.12602434,0 +85,12,0,1.51886789,6.60413424,0 +85,13,0,1.51886789,7.09013758,0 +85,14,0,1.51886789,7.55543954,0 +86,1,0,0.09675537,1.2482013,0 +86,2,0,0.09675537,1.29904246,0 +86,3,0,0.09675537,1.33662018,0 +86,4,0,0.09675537,1.38577264,0 +86,5,0,0.09675537,1.43359688,0 +86,6,0,0.09675537,1.41577091,0 +86,7,0,0.09675537,1.46297865,0 +86,8,0,0.09675537,1.57248394,0 +86,9,0,0.09675537,1.55582926,0 +86,10,0,0.09675537,1.57245625,0 +86,11,0,0.09675537,1.60997446,0 +86,12,0,0.09675537,1.64186738,0 +86,13,0,0.09675537,1.65357958,0 +86,14,0,0.09675537,1.77370839,0 +87,1,0,0.35151612,0.77896482,0 +87,2,0,0.35151612,0.89835531,0 +87,3,0,0.35151612,0.99929643,0 +87,4,0,0.35151612,0.97459127,0 +87,5,0,0.35151612,1.12190902,0 +87,6,0,0.35151612,1.19046073,0 +87,7,0,0.35151612,1.34197261,0 +87,8,0,0.35151612,1.40291105,0 +87,9,0,0.35151612,1.49945419,0 +87,10,0,0.35151612,1.58438214,0 +87,11,0,0.35151612,1.75729215,0 +87,12,0,0.35151612,1.82430529,0 +87,13,0,0.35151612,1.91395051,0 +87,14,0,0.35151612,2.02344833,0 +88,1,0,-0.76876415,-0.27382279,0 +88,2,0,-0.76876415,-0.47931058,0 +88,3,0,-0.76876415,-0.78501531,0 +88,4,0,-0.76876415,-0.94311415,0 +88,5,0,-0.76876415,-1.08641579,0 +88,6,0,-0.76876415,-1.39820096,0 +88,7,0,-0.76876415,-1.56406219,0 +88,8,0,-0.76876415,-1.77590955,0 +88,9,0,-0.76876415,-2.0028704,0 +88,10,0,-0.76876415,-2.20241111,0 +88,11,0,-0.76876415,-2.36212405,0 +88,12,0,-0.76876415,-2.55007676,0 +88,13,0,-0.76876415,-2.78043329,0 +88,14,0,-0.76876415,-3.02269475,0 +89,1,0,-0.84141692,-1.3190862,0 +89,2,0,-0.84141692,-1.52859453,0 +89,3,0,-0.84141692,-1.79411443,0 +89,4,0,-0.84141692,-1.98875227,0 +89,5,0,-0.84141692,-2.2452544,0 +89,6,0,-0.84141692,-2.45771555,0 +89,7,0,-0.84141692,-2.69975584,0 +89,8,0,-0.84141692,-2.94157538,0 +89,9,0,-0.84141692,-3.19121696,0 +89,10,0,-0.84141692,-3.37735773,0 +89,11,0,-0.84141692,-3.63911715,0 +89,12,0,-0.84141692,-3.86336309,0 +89,13,0,-0.84141692,-4.13423596,0 +89,14,0,-0.84141692,-4.3739612,0 +90,1,0,-1.94007321,0.71766548,0 +90,2,0,-1.94007321,0.14726656,0 +90,3,0,-1.94007321,-0.45697469,0 +90,4,0,-1.94007321,-1.04026272,0 +90,5,0,-1.94007321,-1.67658612,0 +90,6,0,-1.94007321,-2.1822746,0 +90,7,0,-1.94007321,-2.80209795,0 +90,8,0,-1.94007321,-3.37384103,0 +90,9,0,-1.94007321,-3.97501072,0 +90,10,0,-1.94007321,-4.58948767,0 +90,11,0,-1.94007321,-5.12846636,0 +90,12,0,-1.94007321,-5.73024914,0 +90,13,0,-1.94007321,-6.37754422,0 +90,14,0,-1.94007321,-6.88660035,0 +91,1,0,-0.28135865,0.63446253,0 +91,2,0,-0.28135865,0.51741378,0 +91,3,0,-0.28135865,0.41016952,0 +91,4,0,-0.28135865,0.34282025,0 +91,5,0,-0.28135865,0.26501162,0 +91,6,0,-0.28135865,0.21785899,0 +91,7,0,-0.28135865,0.14839987,0 +91,8,0,-0.28135865,0.08458606,0 +91,9,0,-0.28135865,-0.01824613,0 +91,10,0,-0.28135865,-0.12732634,0 +91,11,0,-0.28135865,-0.17607901,0 +91,12,0,-0.28135865,-0.24153786,0 +91,13,0,-0.28135865,-0.37472745,0 +91,14,0,-0.28135865,-0.4304573,0 +92,1,0,-0.03358079,-1.50651613,0 +92,2,0,-0.03358079,-1.46556536,0 +92,3,0,-0.03358079,-1.44978833,0 +92,4,0,-0.03358079,-1.49875405,0 +92,5,0,-0.03358079,-1.5106182,0 +92,6,0,-0.03358079,-1.47787629,0 +92,7,0,-0.03358079,-1.52467346,0 +92,8,0,-0.03358079,-1.55862045,0 +92,9,0,-0.03358079,-1.54525806,0 +92,10,0,-0.03358079,-1.62297645,0 +92,11,0,-0.03358079,-1.66605909,0 +92,12,0,-0.03358079,-1.64543643,0 +92,13,0,-0.03358079,-1.65991167,0 +92,14,0,-0.03358079,-1.66466395,0 +93,1,0,-1.11115755,-0.71552524,0 +93,2,0,-1.11115755,-0.95693157,0 +93,3,0,-1.11115755,-1.33594085,0 +93,4,0,-1.11115755,-1.62252967,0 +93,5,0,-1.11115755,-1.95868148,0 +93,6,0,-1.11115755,-2.27834084,0 +93,7,0,-1.11115755,-2.56472384,0 +93,8,0,-1.11115755,-2.89816426,0 +93,9,0,-1.11115755,-3.21122316,0 +93,10,0,-1.11115755,-3.47738722,0 +93,11,0,-1.11115755,-3.78859132,0 +93,12,0,-1.11115755,-4.10999847,0 +93,13,0,-1.11115755,-4.4350847,0 +93,14,0,-1.11115755,-4.76945217,0 +94,1,0,-1.04568632,1.31321029,0 +94,2,0,-1.04568632,1.07623678,0 +94,3,0,-1.04568632,0.68178024,0 +94,4,0,-1.04568632,0.37009652,0 +94,5,0,-1.04568632,0.02157703,0 +94,6,0,-1.04568632,-0.33483955,0 +94,7,0,-1.04568632,-0.63419515,0 +94,8,0,-1.04568632,-0.99341621,0 +94,9,0,-1.04568632,-1.35824697,0 +94,10,0,-1.04568632,-1.66961045,0 +94,11,0,-1.04568632,-2.02685801,0 +94,12,0,-1.04568632,-2.36430934,0 +94,13,0,-1.04568632,-2.70076647,0 +94,14,0,-1.04568632,-3.01136844,0 +95,1,0,-0.23544848,-0.14245627,0 +95,2,0,-0.23544848,-0.16405499,0 +95,3,0,-0.23544848,-0.27024385,0 +95,4,0,-0.23544848,-0.34327877,0 +95,5,0,-0.23544848,-0.37373744,0 +95,6,0,-0.23544848,-0.46206315,0 +95,7,0,-0.23544848,-0.52224244,0 +95,8,0,-0.23544848,-0.58821305,0 +95,9,0,-0.23544848,-0.60467527,0 +95,10,0,-0.23544848,-0.72611125,0 +95,11,0,-0.23544848,-0.79362248,0 +95,12,0,-0.23544848,-0.90606381,0 +95,13,0,-0.23544848,-0.97409729,0 +95,14,0,-0.23544848,-1.03058039,0 +96,1,0,0.84833973,-0.03759139,0 +96,2,0,0.84833973,0.19742433,0 +96,3,0,0.84833973,0.43831584,0 +96,4,0,0.84833973,0.68322093,0 +96,5,0,0.84833973,0.91054192,0 +96,6,0,0.84833973,1.18265408,0 +96,7,0,0.84833973,1.37873188,0 +96,8,0,0.84833973,1.60625765,0 +96,9,0,0.84833973,1.82820988,0 +96,10,0,0.84833973,2.10336725,0 +96,11,0,0.84833973,2.37001919,0 +96,12,0,0.84833973,2.55428253,0 +96,13,0,0.84833973,2.86226233,0 +96,14,0,0.84833973,3.03314674,0 +97,1,0,0.86828264,1.31528369,0 +97,2,0,0.86828264,1.58308227,0 +97,3,0,0.86828264,1.91763087,0 +97,4,0,0.86828264,2.19828222,0 +97,5,0,0.86828264,2.57674233,0 +97,6,0,0.86828264,2.80196005,0 +97,7,0,0.86828264,3.1868346,0 +97,8,0,0.86828264,3.48590138,0 +97,9,0,0.86828264,3.82693004,0 +97,10,0,0.86828264,4.10559545,0 +97,11,0,0.86828264,4.44056949,0 +97,12,0,0.86828264,4.73369915,0 +97,13,0,0.86828264,5.05172798,0 +97,14,0,0.86828264,5.36246126,0 +98,1,0,-0.3132031,-1.00305528,0 +98,2,0,-0.3132031,-1.10866669,0 +98,3,0,-0.3132031,-1.24390296,0 +98,4,0,-0.3132031,-1.28862389,0 +98,5,0,-0.3132031,-1.38995454,0 +98,6,0,-0.3132031,-1.4741515,0 +98,7,0,-0.3132031,-1.57986628,0 +98,8,0,-0.3132031,-1.70793509,0 +98,9,0,-0.3132031,-1.87471202,0 +98,10,0,-0.3132031,-1.91404961,0 +98,11,0,-0.3132031,-2.01264995,0 +98,12,0,-0.3132031,-2.13801949,0 +98,13,0,-0.3132031,-2.19120009,0 +98,14,0,-0.3132031,-2.29580255,0 +99,1,0,-1.38183282,0.20100353,0 +99,2,0,-1.38183282,-0.15928647,0 +99,3,0,-1.38183282,-0.59948375,0 +99,4,0,-1.38183282,-0.94483473,0 +99,5,0,-1.38183282,-1.3274795,0 +99,6,0,-1.38183282,-1.80052299,0 +99,7,0,-1.38183282,-2.180972,0 +99,8,0,-1.38183282,-2.55168549,0 +99,9,0,-1.38183282,-2.96644999,0 +99,10,0,-1.38183282,-3.30142449,0 +99,11,0,-1.38183282,-3.70949243,0 +99,12,0,-1.38183282,-4.1601258,0 +99,13,0,-1.38183282,-4.51845748,0 +99,14,0,-1.38183282,-4.9992723,0 +100,1,0,-1.83193478,0.28495733,0 +100,2,0,-1.83193478,-0.24209999,0 +100,3,0,-1.83193478,-0.79688703,0 +100,4,0,-1.83193478,-1.30894223,0 +100,5,0,-1.83193478,-1.74715788,0 +100,6,0,-1.83193478,-2.32104877,0 +100,7,0,-1.83193478,-2.79077689,0 +100,8,0,-1.83193478,-3.31317856,0 +100,9,0,-1.83193478,-3.8603181,0 +100,10,0,-1.83193478,-4.40524866,0 +100,11,0,-1.83193478,-4.91206247,0 +100,12,0,-1.83193478,-5.39172395,0 +100,13,0,-1.83193478,-5.8875056,0 +100,14,0,-1.83193478,-6.44062223,0 +101,1,0,0.16135302,-0.18067827,0 +101,2,0,0.16135302,-0.15014324,0 +101,3,0,0.16135302,-0.14613974,0 +101,4,0,0.16135302,-0.01177725,0 +101,5,0,0.16135302,-0.02263935,0 +101,6,0,0.16135302,0.04518623,0 +101,7,0,0.16135302,0.05108814,0 +101,8,0,0.16135302,0.10550541,0 +101,9,0,0.16135302,0.21507987,0 +101,10,0,0.16135302,0.30351479,0 +101,11,0,0.16135302,0.36572969,0 +101,12,0,0.16135302,0.37167169,0 +101,13,0,0.16135302,0.43037421,0 +101,14,0,0.16135302,0.45916981,0 +102,1,0,-0.03290985,1.47994119,0 +102,2,0,-0.03290985,1.46383587,0 +102,3,0,-0.03290985,1.47092409,0 +102,4,0,-0.03290985,1.45290577,0 +102,5,0,-0.03290985,1.42982469,0 +102,6,0,-0.03290985,1.47221549,0 +102,7,0,-0.03290985,1.35122109,0 +102,8,0,-0.03290985,1.36601765,0 +102,9,0,-0.03290985,1.43207577,0 +102,10,0,-0.03290985,1.37907767,0 +102,11,0,-0.03290985,1.34627594,0 +102,12,0,-0.03290985,1.37016735,0 +102,13,0,-0.03290985,1.34288116,0 +102,14,0,-0.03290985,1.29231074,0 +103,1,0,-0.22756915,-1.85127946,0 +103,2,0,-0.22756915,-1.90754139,0 +103,3,0,-0.22756915,-2.04761498,0 +103,4,0,-0.22756915,-2.07879877,0 +103,5,0,-0.22756915,-2.16619119,0 +103,6,0,-0.22756915,-2.18158611,0 +103,7,0,-0.22756915,-2.19976851,0 +103,8,0,-0.22756915,-2.26807405,0 +103,9,0,-0.22756915,-2.38306084,0 +103,10,0,-0.22756915,-2.44529512,0 +103,11,0,-0.22756915,-2.51409659,0 +103,12,0,-0.22756915,-2.59393641,0 +103,13,0,-0.22756915,-2.64371077,0 +103,14,0,-0.22756915,-2.67284393,0 +104,1,0,-0.66367978,-0.65504824,0 +104,2,0,-0.66367978,-0.84354219,0 +104,3,0,-0.66367978,-0.99927207,0 +104,4,0,-0.66367978,-1.22890495,0 +104,5,0,-0.66367978,-1.4434219,0 +104,6,0,-0.66367978,-1.60798034,0 +104,7,0,-0.66367978,-1.83044867,0 +104,8,0,-0.66367978,-2.04198263,0 +104,9,0,-0.66367978,-2.23967766,0 +104,10,0,-0.66367978,-2.47448329,0 +104,11,0,-0.66367978,-2.65801481,0 +104,12,0,-0.66367978,-2.87016089,0 +104,13,0,-0.66367978,-3.09490116,0 +104,14,0,-0.66367978,-3.30208896,0 +105,1,0,-0.12998564,0.47655418,0 +105,2,0,-0.12998564,0.42417084,0 +105,3,0,-0.12998564,0.47672522,0 +105,4,0,-0.12998564,0.41010284,0 +105,5,0,-0.12998564,0.39858764,0 +105,6,0,-0.12998564,0.41650573,0 +105,7,0,-0.12998564,0.2881141,0 +105,8,0,-0.12998564,0.35622521,0 +105,9,0,-0.12998564,0.32179954,0 +105,10,0,-0.12998564,0.31047824,0 +105,11,0,-0.12998564,0.36952042,0 +105,12,0,-0.12998564,0.27352985,0 +105,13,0,-0.12998564,0.29564857,0 +105,14,0,-0.12998564,0.29627155,0 +106,1,0,0.17425923,0.79132595,0 +106,2,0,0.17425923,0.85020735,0 +106,3,0,0.17425923,0.87423276,0 +106,4,0,0.17425923,0.93790857,0 +106,5,0,0.17425923,0.97185118,0 +106,6,0,0.17425923,1.03353323,0 +106,7,0,0.17425923,1.04953844,0 +106,8,0,0.17425923,1.07849014,0 +106,9,0,0.17425923,1.06152904,0 +106,10,0,0.17425923,1.12616937,0 +106,11,0,0.17425923,1.17137951,0 +106,12,0,0.17425923,1.20235388,0 +106,13,0,0.17425923,1.24250771,0 +106,14,0,0.17425923,1.26389897,0 +107,1,0,1.01070341,-1.17815734,0 +107,2,0,1.01070341,-0.79938146,0 +107,3,0,1.01070341,-0.52821709,0 +107,4,0,1.01070341,-0.17201651,0 +107,5,0,1.01070341,0.16836286,0 +107,6,0,1.01070341,0.47573517,0 +107,7,0,1.01070341,0.75867997,0 +107,8,0,1.01070341,1.12274261,0 +107,9,0,1.01070341,1.47012755,0 +107,10,0,1.01070341,1.85180945,0 +107,11,0,1.01070341,2.11605264,0 +107,12,0,1.01070341,2.50782584,0 +107,13,0,1.01070341,2.82722193,0 +107,14,0,1.01070341,3.14960956,0 +108,1,0,-1.22001302,0.48146911,0 +108,2,0,-1.22001302,0.07547274,0 +108,3,0,-1.22001302,-0.30025403,0 +108,4,0,-1.22001302,-0.71948493,0 +108,5,0,-1.22001302,-1.15089627,0 +108,6,0,-1.22001302,-1.5751796,0 +108,7,0,-1.22001302,-2.00246622,0 +108,8,0,-1.22001302,-2.34777328,0 +108,9,0,-1.22001302,-2.83661294,0 +108,10,0,-1.22001302,-3.19276087,0 +108,11,0,-1.22001302,-3.63662429,0 +108,12,0,-1.22001302,-4.04402714,0 +108,13,0,-1.22001302,-4.44230456,0 +108,14,0,-1.22001302,-4.89627854,0 +109,1,0,-1.12724787,-1.2626749,0 +109,2,0,-1.12724787,-1.679597,0 +109,3,0,-1.12724787,-1.96277421,0 +109,4,0,-1.12724787,-2.33207189,0 +109,5,0,-1.12724787,-2.62711136,0 +109,6,0,-1.12724787,-3.01028434,0 +109,7,0,-1.12724787,-3.35914544,0 +109,8,0,-1.12724787,-3.70140283,0 +109,9,0,-1.12724787,-4.06324301,0 +109,10,0,-1.12724787,-4.40966793,0 +109,11,0,-1.12724787,-4.87116532,0 +109,12,0,-1.12724787,-5.12397026,0 +109,13,0,-1.12724787,-5.45632416,0 +109,14,0,-1.12724787,-5.83354097,0 +110,1,0,0.17000237,1.54567662,0 +110,2,0,0.17000237,1.58744864,0 +110,3,0,0.17000237,1.60512901,0 +110,4,0,0.17000237,1.62867298,0 +110,5,0,0.17000237,1.60846368,0 +110,6,0,0.17000237,1.63549388,0 +110,7,0,0.17000237,1.73782842,0 +110,8,0,0.17000237,1.75537707,0 +110,9,0,0.17000237,1.7742746,0 +110,10,0,0.17000237,1.80088507,0 +110,11,0,0.17000237,1.76723413,0 +110,12,0,0.17000237,1.8782747,0 +110,13,0,0.17000237,1.85492038,0 +110,14,0,0.17000237,1.90954896,0 +111,1,0,1.13189732,-0.40196854,0 +111,2,0,1.13189732,-0.03834484,0 +111,3,0,1.13189732,0.32443231,0 +111,4,0,1.13189732,0.69948509,0 +111,5,0,1.13189732,1.07419676,0 +111,6,0,1.13189732,1.40940451,0 +111,7,0,1.13189732,1.82680992,0 +111,8,0,1.13189732,2.11649698,0 +111,9,0,1.13189732,2.47146142,0 +111,10,0,1.13189732,2.79929226,0 +111,11,0,1.13189732,3.18965164,0 +111,12,0,1.13189732,3.48418655,0 +111,13,0,1.13189732,3.88613776,0 +111,14,0,1.13189732,4.25161611,0 +112,1,0,-0.91959659,-0.33898614,0 +112,2,0,-0.91959659,-0.63664352,0 +112,3,0,-0.91959659,-0.95293423,0 +112,4,0,-0.91959659,-1.23042649,0 +112,5,0,-0.91959659,-1.54660786,0 +112,6,0,-0.91959659,-1.87473563,0 +112,7,0,-0.91959659,-2.09209356,0 +112,8,0,-0.91959659,-2.42648691,0 +112,9,0,-0.91959659,-2.80821113,0 +112,10,0,-0.91959659,-3.02623721,0 +112,11,0,-0.91959659,-3.30499623,0 +112,12,0,-0.91959659,-3.68790929,0 +112,13,0,-0.91959659,-3.96045286,0 +112,14,0,-0.91959659,-4.25681575,0 +113,1,0,0.21203234,0.37971994,0 +113,2,0,0.21203234,0.46616862,0 +113,3,0,0.21203234,0.41248281,0 +113,4,0,0.21203234,0.51576116,0 +113,5,0,0.21203234,0.52126912,0 +113,6,0,0.21203234,0.63808167,0 +113,7,0,0.21203234,0.67799792,0 +113,8,0,0.21203234,0.76313164,0 +113,9,0,0.21203234,0.78242812,0 +113,10,0,0.21203234,0.85692304,0 +113,11,0,0.21203234,0.90513849,0 +113,12,0,0.21203234,0.9357872,0 +113,13,0,0.21203234,1.00272189,0 +113,14,0,0.21203234,1.07363208,0 +114,1,0,-0.51235137,-0.44576762,0 +114,2,0,-0.51235137,-0.50462483,0 +114,3,0,-0.51235137,-0.73095194,0 +114,4,0,-0.51235137,-0.84582122,0 +114,5,0,-0.51235137,-0.97891404,0 +114,6,0,-0.51235137,-1.14029351,0 +114,7,0,-0.51235137,-1.28060655,0 +114,8,0,-0.51235137,-1.43727429,0 +114,9,0,-0.51235137,-1.62245179,0 +114,10,0,-0.51235137,-1.65128442,0 +114,11,0,-0.51235137,-1.80134406,0 +114,12,0,-0.51235137,-2.00996569,0 +114,13,0,-0.51235137,-2.1310453,0 +114,14,0,-0.51235137,-2.25309008,0 diff --git a/benchmarks/data/cbwsdid_golden.json b/benchmarks/data/cbwsdid_golden.json new file mode 100644 index 00000000..d26dd41b --- /dev/null +++ b/benchmarks/data/cbwsdid_golden.json @@ -0,0 +1,14 @@ +{ + "meta": { + "package": "cbwsdid", + "R_version": "R version 4.5.2 (2025-10-31)", + "panel": "benchmarks/data/cbwsdid_balance_panel.csv", + "estimator": "cbwsdid(design='absorbing', refinement.method='weightit', method='ebal', covs.formula=~x)", + "kappa": [-2, 2] + }, + "dynamic": { + "event_time": [-2, 0, 1, 2], + "estimate": [0.008614905402558678, 1.998507474812239, 2.000435324046538, 1.977837212056277], + "std_error": [0.05870423491347857, 0.05704878175034945, 0.1113623933607992, 0.1700804801053713] + } +} diff --git a/diff_diff/balancing.py b/diff_diff/balancing.py new file mode 100644 index 00000000..98cabd7c --- /dev/null +++ b/diff_diff/balancing.py @@ -0,0 +1,209 @@ +"""Covariate balancing weights (entropy balancing). + +Implements **entropy balancing** (Hainmueller, J. (2012). "Entropy Balancing for +Causal Effects: A Multivariate Reweighting Method to Produce Balanced Samples in +Observational Studies." *Political Analysis*, 20(1), 25-46. +https://doi.org/10.1093/pan/mpr025). + +Entropy balancing finds nonnegative control weights ``w_i`` that exactly match a set +of target covariate moments (here: the treated-group covariate means) while staying as +close as possible — in the Kullback-Leibler sense — to a set of base weights (uniform by +default). The solution is obtained from the convex dual + + minimize over λ: L(λ) = log( Σ_i q_i exp(Z_iᵀ λ) ), Z_i = X_i − target, + +whose stationary point ``∇L = Σ_i w_i Z_i = 0`` is exactly first-moment balance, with +``w_i = q_i exp(Z_iᵀ λ) / Σ_j q_j exp(Z_jᵀ λ)``. ``L`` is convex (log-sum-exp), so a +damped Newton iteration (gradient = weighted mean of the centered moments, Hessian = +weighted covariance) converges to the balancing weights whenever the target lies in the +interior of the control covariate convex hull; otherwise no finite λ balances the +moments and the problem is **infeasible**. + +This module is dependency-light (numpy, with an optional scipy L-BFGS fallback) and is +used by ``StackedDiD`` to construct the within-sub-experiment design weights ``b_{sa}`` +for Covariate-Balanced Weighted Stacked DID (Ustyuzhanin 2026). +""" + +from __future__ import annotations + +from typing import Any, Dict, Optional, Tuple + +import numpy as np + +__all__ = ["entropy_balance", "BalanceError"] + + +class BalanceError(ValueError): + """Raised when entropy balancing fails to achieve first-moment balance. + + Carries the achieved ``max_residual`` and the per-covariate residual vector so + callers (e.g. ``StackedDiD``) can attach cohort context and report the worst- + balanced covariate. + """ + + def __init__(self, message: str, *, max_residual: float, residuals: np.ndarray): + super().__init__(message) + self.max_residual = max_residual + self.residuals = residuals + + +def entropy_balance( + X: np.ndarray, + target_means: np.ndarray, + base_weights: Optional[np.ndarray] = None, + *, + max_iter: int = 200, + tol: float = 1e-8, +) -> Tuple[np.ndarray, Dict[str, Any]]: + """Solve entropy balancing for control weights matching ``target_means``. + + Parameters + ---------- + X : np.ndarray, shape (n, k) + Covariate matrix for the ``n`` control units (``k`` covariates). + target_means : np.ndarray, shape (k,) + Target first moments to match — for CBWSDID these are the treated-group + covariate means. + base_weights : np.ndarray, shape (n,), optional + Nonnegative base weights ``q_i`` (the KL reference). Defaults to uniform. + Internally renormalized to sum to one. + max_iter : int + Maximum damped-Newton iterations (then a scipy L-BFGS fallback is attempted). + tol : float + Convergence tolerance on the maximum absolute (raw-scale) moment residual + ``max_r |Σ_i w_i X_{i,r} − target_r|``. + + Returns + ------- + weights : np.ndarray, shape (n,) + Nonnegative weights summing to one with ``Σ_i w_i X_i ≈ target_means``. + info : dict + ``converged`` (bool), ``max_residual`` (float), ``n_iter`` (int), + ``ess`` (effective sample size ``1 / Σ_i w_i²``), ``solver`` (str). + + Raises + ------ + BalanceError + If neither the damped-Newton nor the L-BFGS fallback drives the maximum moment + residual below ``tol`` (the target is outside the control covariate hull, i.e. + infeasible). + ValueError + On malformed inputs (shape mismatch, non-finite, negative base weights). + """ + X = np.asarray(X, dtype=np.float64) + if X.ndim != 2: + raise ValueError(f"X must be 2-D (n, k); got shape {X.shape}") + n, k = X.shape + target = np.asarray(target_means, dtype=np.float64).reshape(-1) + if target.shape[0] != k: + raise ValueError(f"target_means length {target.shape[0]} != n_covariates {k}") + if not np.all(np.isfinite(X)) or not np.all(np.isfinite(target)): + raise ValueError("X and target_means must be finite") + if n == 0: + raise ValueError("X has no control rows") + + if base_weights is None: + q = np.full(n, 1.0 / n) + else: + q = np.asarray(base_weights, dtype=np.float64).reshape(-1) + if q.shape[0] != n: + raise ValueError(f"base_weights length {q.shape[0]} != n_control {n}") + if np.any(q < 0) or not np.all(np.isfinite(q)): + raise ValueError("base_weights must be nonnegative and finite") + s = q.sum() + if s <= 0: + raise ValueError("base_weights sum to zero") + q = q / s + + # Centered moments; standardize columns for conditioning (balance set is invariant + # to the linear rescaling — it is absorbed into the dual variable λ). + Z = X - target + scale = Z.std(axis=0) + scale[scale < 1e-12] = 1.0 + Zs = Z / scale + + def weights_at(lam: np.ndarray) -> np.ndarray: + logits = Zs @ lam + logits -= logits.max() + ew = q * np.exp(logits) + return ew / ew.sum() + + def dual_loss(lam: np.ndarray) -> float: + logits = Zs @ lam + m = logits.max() + return float(m + np.log(np.sum(q * np.exp(logits - m)))) + + def raw_residual(w: np.ndarray) -> np.ndarray: + return w @ X - target + + lam = np.zeros(k) + solver = "newton" + n_iter = 0 + for n_iter in range(1, max_iter + 1): + w = weights_at(lam) + if np.max(np.abs(raw_residual(w))) < tol: + break + g = w @ Zs # gradient of the dual loss (standardized scale) + Zc = Zs - g + H = (w[:, None] * Zc).T @ Zc # weighted covariance (PSD) + ridge = 1e-10 * (np.trace(H) / k + 1e-12) + try: + direction = -np.linalg.solve(H + ridge * np.eye(k), g) + except np.linalg.LinAlgError: + direction = -np.linalg.lstsq(H, g, rcond=None)[0] + # Backtracking (Armijo) line search on the convex dual loss. + base = dual_loss(lam) + slope = float(g @ direction) # < 0 (descent) + step = 1.0 + for _ in range(40): + if dual_loss(lam + step * direction) <= base + 1e-4 * step * slope: + break + step *= 0.5 + lam = lam + step * direction + else: + w = weights_at(lam) + + w = weights_at(lam) + resid = raw_residual(w) + max_resid = float(np.max(np.abs(resid))) + + if max_resid >= tol: + # Fallback: scipy L-BFGS-B on the convex dual (robust to poor Newton scaling). + try: + from scipy.optimize import minimize + + res = minimize( + dual_loss, + lam, + jac=lambda L: weights_at(L) @ Zs, + method="L-BFGS-B", + options={"maxiter": 500, "gtol": 1e-12}, + ) + w_lbfgs = weights_at(res.x) + resid_lbfgs = raw_residual(w_lbfgs) + if np.max(np.abs(resid_lbfgs)) < max_resid: + lam, w, resid = res.x, w_lbfgs, resid_lbfgs + max_resid = float(np.max(np.abs(resid))) + solver = "lbfgs" + except Exception: # pragma: no cover - scipy always present, defensive + pass + + converged = max_resid < tol + info: Dict[str, Any] = { + "converged": converged, + "max_residual": max_resid, + "n_iter": n_iter, + "ess": float(1.0 / np.sum(w**2)), + "solver": solver, + } + if not converged: + worst = int(np.argmax(np.abs(resid))) + raise BalanceError( + "entropy balancing did not converge to first-moment balance " + f"(max moment residual {max_resid:.3e} >= tol {tol:.1e}; worst covariate " + f"index {worst}). The target mean is likely outside the convex hull of the " + "control covariates (infeasible).", + max_residual=max_resid, + residuals=resid, + ) + return w, info diff --git a/diff_diff/guides/llms-full.txt b/diff_diff/guides/llms-full.txt index 14b342d3..8a35b773 100644 --- a/diff_diff/guides/llms-full.txt +++ b/diff_diff/guides/llms-full.txt @@ -816,6 +816,7 @@ StackedDiD( anticipation: int = 0, rank_deficient_action: str = "warn", vcov_type: str = "hc1", # {"hc1","hc2_bm"}; classical/hc2 rejected (intrinsically clustered), conley deferred. survey_design=... requires hc1 + balance: str = "none", # {"none","entropy"}; "entropy" = CBWSDID covariate balancing (Ustyuzhanin 2026), requires fit(covariates=[...]) + weighting="aggregate", no survey_design ) ``` @@ -832,6 +833,7 @@ stacked.fit( first_treat: str, aggregate: str = None, # None, "simple", or "event_study" population: str = None, # Required when weighting="population" + covariates: list[str] = None, # Columns to balance (requires balance="entropy"); values read at t=a-1-anticipation; balanced windows only ) -> StackedDiDResults ``` diff --git a/diff_diff/guides/llms.txt b/diff_diff/guides/llms.txt index 44a218c4..54ea1593 100644 --- a/diff_diff/guides/llms.txt +++ b/diff_diff/guides/llms.txt @@ -64,7 +64,7 @@ Full practitioner guide: call `diff_diff.get_llm_guide("practitioner")` - [TripleDifference](https://diff-diff.readthedocs.io/en/stable/api/triple_diff.html): Triple difference (DDD) estimator for designs requiring two criteria for treatment eligibility - [ContinuousDiD](https://diff-diff.readthedocs.io/en/stable/api/continuous_did.html): Callaway, Goodman-Bacon & Sant'Anna (2024) continuous treatment DiD with dose-response curves - [HeterogeneousAdoptionDiD](https://diff-diff.readthedocs.io/en/stable/api/had.html): de Chaisemartin, Ciccia, D'Haultfœuille & Knau (2026) for designs where **no unit remains untreated**; local-linear estimator at the dose support boundary returning Weighted Average Slope (WAS) on Design 1' (`d̲=0` / QUG) or `WAS_{d̲}` on Design 1 (`d̲>0`, continuous-near-d̲ or mass-point), with multi-period event-study extension (last-treatment cohort, pointwise CIs). **Panel-only** in this release (repeated cross-sections rejected by the validator). Alias `HAD`. -- [StackedDiD](https://diff-diff.readthedocs.io/en/stable/api/stacked_did.html): Wing, Freedman & Hollingsworth (2024) stacked DiD with Q-weights and sub-experiments +- [StackedDiD](https://diff-diff.readthedocs.io/en/stable/api/stacked_did.html): Wing, Freedman & Hollingsworth (2024) stacked DiD with Q-weights and sub-experiments; optional covariate balancing (`balance="entropy"`, Ustyuzhanin 2026) - [EfficientDiD](https://diff-diff.readthedocs.io/en/stable/api/efficient_did.html): Chen, Sant'Anna & Xie (2025) efficient DiD with optimal weighting for tighter SEs - [TROP](https://diff-diff.readthedocs.io/en/stable/api/trop.html): Triply Robust Panel estimator (Athey et al. 2025) with nuclear norm factor adjustment - [StaggeredTripleDifference](https://diff-diff.readthedocs.io/en/stable/api/staggered.html#staggeredtripledifference): Ortiz-Villavicencio & Sant'Anna (2025) staggered DDD with group-time ATT diff --git a/diff_diff/stacked_did.py b/diff_diff/stacked_did.py index 698be759..20dd2f47 100644 --- a/diff_diff/stacked_did.py +++ b/diff_diff/stacked_did.py @@ -24,6 +24,7 @@ import numpy as np import pandas as pd +from diff_diff.balancing import BalanceError, entropy_balance from diff_diff.linalg import solve_ols from diff_diff.stacked_did_results import StackedDiDResults # noqa: F401 (re-export) from diff_diff.utils import safe_inference @@ -114,6 +115,23 @@ class StackedDiD: raised — the survey Taylor-series linearization (or replicate-weight refit) variance overrides the analytical sandwich. Use the default ``vcov_type="hc1"`` for survey designs. + balance : {"none", "entropy"}, default="none" + Within-sub-experiment covariate balancing (Covariate-Balanced Weighted + Stacked DID; Ustyuzhanin 2026). With ``"entropy"`` and a ``fit(..., + covariates=[...])`` list, each clean-control group is reweighted by + entropy balancing (Hainmueller 2012) so its covariate means match the + treated cohort's (measured at the last pre-treatment period), and the + resulting design weights ``b_sa`` are composed with the Wing corrective + weights via the effective control mass into the final stacked weights + ``W_sa``. This is **control-only reweighting**, so it preserves the + trimmed-aggregate-ATT estimand (it changes only how untreated trends are + estimated, not the treated-cohort weights); at ``b_sa=1`` it reduces to + the paper's unit-count weighted stacked DID, equal to + ``weighting="aggregate"`` on balanced event windows. v1 requires + ``weighting="aggregate"`` and **balanced event windows** (ragged windows + raise a ``ValueError``), and does not support ``survey_design=``; + matching-based balancing and the repeated-treatment extension are out of + scope. Default ``"none"`` reproduces plain weighted stacked DID. Attributes ---------- @@ -165,6 +183,7 @@ def __init__( anticipation: int = 0, rank_deficient_action: str = "warn", vcov_type: str = "hc1", + balance: str = "none", ): if weighting not in ("aggregate", "population", "sample_share"): raise ValueError( @@ -186,6 +205,7 @@ def __init__( # vcov_type validation (Phase 1b 2/8: thread through StackedDiD). # Factored into _validate_vcov_type so set_params() can re-validate. self._validate_vcov_type(vcov_type) + self._validate_balance(balance) self.kappa_pre = kappa_pre self.kappa_post = kappa_post @@ -196,6 +216,7 @@ def __init__( self.anticipation = anticipation self.rank_deficient_action = rank_deficient_action self.vcov_type = vcov_type + self.balance = balance self.is_fitted_ = False self.results_: Optional[StackedDiDResults] = None @@ -235,6 +256,20 @@ def _validate_vcov_type(vcov_type: str) -> None: "'hc2_bm' (CR2 Bell-McCaffrey)." ) + @staticmethod + def _validate_balance(balance: str) -> None: + """Validate the covariate-balancing method (CBWSDID, Ustyuzhanin 2026). + + Called from __init__ AND set_params (mirrors _validate_vcov_type) so + sklearn-style mutation hits the estimator-level guard. v1 supports only + entropy balancing; matching-based balancing and IPW are out of scope + (see docs/methodology/REGISTRY.md StackedDiD).""" + if balance not in ("none", "entropy"): + raise ValueError( + f"balance must be 'none' or 'entropy', got '{balance}'. " + "Matching-based balancing and IPW are out of scope for v1." + ) + def fit( self, data: pd.DataFrame, @@ -245,6 +280,7 @@ def fit( aggregate: Optional[str] = None, population: Optional[str] = None, survey_design=None, + covariates: Optional[List[str]] = None, ) -> StackedDiDResults: """ Fit the stacked DiD estimator. @@ -275,6 +311,14 @@ def fit( Survey design specification for design-based inference. When provided, uses Taylor Series Linearization for variance estimation and applies sampling weights to the regression. + covariates : list of str, optional + Covariate column names to balance the clean controls toward the + treated cohort (requires ``balance="entropy"``; see the constructor + ``balance`` parameter). Values are read at the last pre-treatment + period ``t = a-1-anticipation`` per sub-experiment, so balancing uses + only pre-treatment information (Assumption 4). Raises ``ValueError`` + if ``balance="none"`` (or vice versa), if a name is absent from + ``data``, or if a cohort cannot be balanced (infeasible). Returns ------- @@ -309,6 +353,44 @@ def fit( if self.weighting == "population" and population is None: raise ValueError("population column must be specified when weighting='population'") + # ---- Covariate balancing (CBWSDID, Ustyuzhanin 2026) validation + guards ---- + if isinstance(covariates, str): + raise TypeError( + "covariates must be a list of column names, not a string (got " + f"{covariates!r}). Use covariates=[{covariates!r}]." + ) + balancing = self.balance != "none" + if balancing and not covariates: + raise ValueError( + "balance='entropy' requires a non-empty covariates= list (the " + "columns to balance the clean controls toward the treated cohort). " + "Use balance='none' for unrefined weighted stacked DID." + ) + if covariates and not balancing: + raise ValueError( + "covariates= was provided but balance='none'. Set balance='entropy' " + "to enable covariate balancing, or drop covariates=." + ) + if balancing: + assert covariates is not None # guaranteed by the cross-validation above + # Deduplicate (repeated columns are redundant moments) while preserving order. + covariates = list(dict.fromkeys(covariates)) + if self.weighting != "aggregate": + raise NotImplementedError( + f"balance='entropy' is only supported with weighting='aggregate' " + f"(got weighting='{self.weighting}'); the CBWSDID corrective weight " + "uses the Wing aggregate (treated-share) form. v1 scope." + ) + if survey_design is not None: + raise NotImplementedError( + "balance='entropy' with survey_design= is not supported in v1 " + "(design-weight + survey-weight composition is out of scope). " + "Drop survey_design= or set balance='none'." + ) + missing_cov = [c for c in covariates if c not in data.columns] + if missing_cov: + raise ValueError(f"covariates not found in data columns: {missing_cov}") + # ---- Resolve survey design ---- from diff_diff.survey import ( SurveyDesign, @@ -432,6 +514,17 @@ def fit( # ---- Compute Q-weights ---- stacked_df = self._compute_q_weights(stacked_df, unit, population) + # ---- Covariate balancing: design weights b_sa -> effective-mass W_sa ---- + # When balancing, this OVERWRITES `_Q_weight` with the CBWSDID final weights + # W_sa (paper §3.1) so the existing WLS path downstream consumes them + # transparently; the raw design weights are preserved in `_b_sa`. + balance_diagnostics: Optional[Dict[Any, Dict[str, Any]]] = None + if balancing: + assert covariates is not None # narrowed by the cross-validation above + stacked_df, balance_diagnostics = self._compute_balancing_weights( + stacked_df, df, unit, time, first_treat, covariates + ) + # ---- Count units ---- treated_units = stacked_df.loc[stacked_df["_D_sa"] == 1, unit].unique() control_units = stacked_df.loc[stacked_df["_D_sa"] == 0, unit].unique() @@ -845,6 +938,9 @@ def _refit_stacked(w_r): cluster_name=self.cluster, n_clusters=int(np.unique(cluster_ids).size), survey_metadata=survey_metadata, + balance=self.balance, + covariates=list(covariates) if balancing else None, + balance_diagnostics=balance_diagnostics, ) self.is_fitted_ = True @@ -1187,6 +1283,174 @@ def _compute_q_weights_aggregate(self, stacked_df: pd.DataFrame) -> pd.DataFrame stacked_df["_Q_weight"] = weights return stacked_df + # ========================================================================= + # Covariate balancing (CBWSDID, Ustyuzhanin 2026) + # ========================================================================= + + def _compute_balancing_weights( + self, + stacked_df: pd.DataFrame, + df: pd.DataFrame, + unit: str, + time: str, + first_treat: str, + covariates: List[str], + ) -> Tuple[pd.DataFrame, Dict[Any, Dict[str, Any]]]: + """Compute CBWSDID covariate-balancing weights (Ustyuzhanin 2026, §3.1). + + For each sub-experiment ``a``, balance the clean controls' covariate means — + measured at the last pre-treatment period ``t = a-1-anticipation`` from the + SOURCE data, so the design weights use only pre-treatment information + (Assumption 4) — toward the treated-cohort means via entropy balancing, + yielding nonnegative design weights ``b_sa`` (treated keep ``b=1``). Then + compose the final stacked weights with the EFFECTIVE control mass + ``Ñ^C_a = Σ_{s∈C_a} b_sa``:: + + W_sa = b_sa · (N^D_a / N^D_Ω) / (Ñ^C_a / Ñ^C_Ω) for s ∈ C_a + W_sa = 1 for s ∈ D_a + + A naive ``b_sa · Q_aggregate`` multiply is **not** equivalent: it aggregates the + cohort control means with weights ∝ (N^D_a/N^D_Ω)·(Ñ^C_a/N^C_a) instead of the + required ∝ (N^D_a/N^D_Ω), biasing the estimate unless ``b_sa`` is uniform. + + Overwrites ``_Q_weight`` with W_sa (so the downstream WLS consumes it + transparently) and records the raw design weights in ``_b_sa``. Fail-closed: + raises a cohort-named ``ValueError`` if a cohort's balance is infeasible — + dropping the cohort would silently shift the estimand to an overlap-trimmed ATT. + """ + balance_tol = 1e-8 + sub_exps = list(pd.unique(stacked_df["_sub_exp"])) + + b_lookup: Dict[Tuple[Any, Any], float] = {} + N_D: Dict[Any, float] = {} + Nt_C: Dict[Any, float] = {} + diagnostics: Dict[Any, Dict[str, Any]] = {} + + # Balanced event windows are required: the paper's unit-count corrective and + # diff-diff's observation-count "aggregate" Q-weights coincide only when every + # unit is observed at every event time. On ragged windows they diverge (the + # count-convention is unresolved — out of scope for v1), so fail closed rather + # than silently producing unit-count estimates that differ from balance="none". + # The check validates exact (unit x event_time) coverage, not just row counts: + # it catches (a) eligible units with zero rows in the window (silently dropped + # by _build_sub_experiment, so invisible to a count), (b) wrong row counts, and + # (c) duplicate (unit, event_time) rows that a count-only check would let pass + # alongside a compensating missing row. + expected_events = list(range(-self.kappa_pre - self.anticipation, self.kappa_post + 1)) + n_expected = len(expected_events) + ft_by_unit = df.drop_duplicates(subset=[unit]).set_index(unit)[first_treat] + + def _expected_units(a_val: Any) -> set: + treated = set(ft_by_unit[ft_by_unit == a_val].index) + if self.clean_control == "not_yet_treated": + controls = set(ft_by_unit[ft_by_unit > a_val + self.kappa_post].index) + elif self.clean_control == "strict": + controls = set( + ft_by_unit[ft_by_unit > a_val + self.kappa_post + self.kappa_pre].index + ) + else: # never_treated + controls = set(ft_by_unit[np.isinf(ft_by_unit)].index) + return treated | controls + + for a in sub_exps: + sub = stacked_df[stacked_df["_sub_exp"] == a] + counts = sub.groupby(unit).size() + present_units = set(counts.index) + missing_eligible = _expected_units(a) - present_units # zero-row eligible units + wrong_count = set(counts[counts != n_expected].index) + dup_mask = sub.duplicated(subset=[unit, "_event_time"], keep=False) + dup_units = set(sub.loc[dup_mask, unit].unique()) + bad = missing_eligible | wrong_count | dup_units + if bad: + raise ValueError( + f"balance='entropy' requires balanced event windows, but cohort a={a} " + f"has {len(bad)} treated/clean-control unit(s) without exact coverage of " + f"the {n_expected} event times {expected_events} (zero-row, missing/extra, " + f"or duplicated (unit, event_time) rows; e.g. {list(bad)[:3]}). Covariate " + "balancing on unbalanced/ragged panels is out of scope for v1 because the " + "paper's unit-count corrective and diff-diff's observation-count " + "'aggregate' Q-weights diverge off balanced panels. Use balance='none', or " + "restrict to a balanced window." + ) + treated_units = list(pd.unique(sub.loc[sub["_D_sa"] == 1, unit])) + control_units = list(pd.unique(sub.loc[sub["_D_sa"] == 0, unit])) + + ref_time = a - 1 - self.anticipation + pre = df.loc[df[time] == ref_time].drop_duplicates(subset=[unit]).set_index(unit) + missing_units = [u for u in treated_units + control_units if u not in pre.index] + if missing_units: + raise ValueError( + f"Covariate balancing for cohort a={a}: {len(missing_units)} unit(s) " + f"have no observation at the pre-treatment reference period " + f"t={ref_time} (e.g. {missing_units[:3]}). Balancing requires the " + "covariate values at t=a-1-anticipation for every treated and " + "clean-control unit." + ) + Xt = pre.loc[treated_units, covariates].to_numpy(dtype=np.float64) + Xc = pre.loc[control_units, covariates].to_numpy(dtype=np.float64) + if not (np.all(np.isfinite(Xt)) and np.all(np.isfinite(Xc))): + raise ValueError( + f"Covariate balancing for cohort a={a}: covariates contain NaN/inf at " + f"t={ref_time}. Balancing requires finite covariate values." + ) + + target = Xt.mean(axis=0) + pre_imbalance = float(np.max(np.abs(Xc.mean(axis=0) - target))) + try: + b_control, info = entropy_balance(Xc, target, tol=balance_tol) + except BalanceError as exc: + worst = covariates[int(np.argmax(np.abs(exc.residuals)))] + raise ValueError( + f"Covariate balancing failed for cohort a={a}: could not match the " + f"treated covariate means (worst covariate '{worst}', residual " + f"{exc.max_residual:.3e}). The treated cohort's covariate profile lies " + "outside the clean-control support (infeasible). Remove this cohort, " + "reduce the covariate set, or use balance='none'." + ) from exc + + for u in treated_units: + b_lookup[(a, u)] = 1.0 + for u, b in zip(control_units, b_control): + b_lookup[(a, u)] = float(b) + N_D[a] = float(len(treated_units)) + Nt_C[a] = float(np.sum(b_control)) + + ess = float(info["ess"]) + if ess < max(2.0, 0.05 * len(control_units)): + warnings.warn( + f"Covariate balancing for cohort a={a} produced highly concentrated " + f"control weights (effective sample size {ess:.1f} of " + f"{len(control_units)} controls); estimates for this cohort may be " + "unstable.", + UserWarning, + stacklevel=3, + ) + diagnostics[a] = { + "n_treated": int(len(treated_units)), + "n_control": int(len(control_units)), + "effective_control_mass": Nt_C[a], + "ess": ess, + "max_imbalance_pre": pre_imbalance, + "max_imbalance_post": float(info["max_residual"]), + "balance_solver": info["solver"], + } + + N_D_Omega = sum(N_D.values()) + Nt_C_Omega = sum(Nt_C.values()) + corr = {a: (N_D[a] / N_D_Omega) / (Nt_C[a] / Nt_C_Omega) for a in sub_exps} + + sub_vals = stacked_df["_sub_exp"].to_numpy() + unit_vals = stacked_df[unit].to_numpy() + d_vals = stacked_df["_D_sa"].to_numpy() + b_vals = np.array([b_lookup[(sub_vals[i], unit_vals[i])] for i in range(len(stacked_df))]) + corr_vals = np.array([corr[a] for a in sub_vals]) + W = np.where(d_vals == 0, b_vals * corr_vals, 1.0) + + stacked_df = stacked_df.copy() + stacked_df["_b_sa"] = b_vals + stacked_df["_Q_weight"] = W + return stacked_df, diagnostics + # ========================================================================= # sklearn-compatible interface # ========================================================================= @@ -1203,6 +1467,7 @@ def get_params(self) -> Dict[str, Any]: "anticipation": self.anticipation, "rank_deficient_action": self.rank_deficient_action, "vcov_type": self.vcov_type, + "balance": self.balance, } def set_params(self, **params: Any) -> "StackedDiD": @@ -1217,6 +1482,8 @@ def set_params(self, **params: Any) -> "StackedDiD": # error surface as __init__ applies. if "vcov_type" in params: self._validate_vcov_type(params["vcov_type"]) + if "balance" in params: + self._validate_balance(params["balance"]) for key, value in params.items(): if hasattr(self, key): setattr(self, key, value) @@ -1252,6 +1519,7 @@ def stacked_did( aggregate: Optional[str] = None, population: Optional[str] = None, survey_design=None, + covariates: Optional[List[str]] = None, **kwargs: Any, ) -> StackedDiDResults: """ @@ -1281,8 +1549,13 @@ def stacked_did( Population column for weighting="population". survey_design : SurveyDesign, optional Survey design specification for design-based inference. + covariates : list of str, optional + Covariate columns to balance the clean controls toward the treated + cohort (pass ``balance="entropy"`` via ``**kwargs`` to enable). See + ``StackedDiD.fit``. **kwargs - Additional keyword arguments passed to StackedDiD constructor. + Additional keyword arguments passed to StackedDiD constructor + (e.g. ``balance="entropy"``, ``weighting``, ``cluster``, ``vcov_type``). Returns ------- @@ -1308,4 +1581,5 @@ def stacked_did( aggregate=aggregate, population=population, survey_design=survey_design, + covariates=covariates, ) diff --git a/diff_diff/stacked_did_results.py b/diff_diff/stacked_did_results.py index 5a5cbeed..52a60ddf 100644 --- a/diff_diff/stacked_did_results.py +++ b/diff_diff/stacked_did_results.py @@ -107,6 +107,16 @@ class StackedDiDResults: n_clusters: Optional[int] = None # Survey design metadata (SurveyMetadata instance from diff_diff.survey) survey_metadata: Optional[Any] = field(default=None) + # --- Covariate balancing (CBWSDID, Ustyuzhanin 2026) --- + # balance: "none" (default, plain weighted stacked DID) or "entropy". When + # "entropy", `covariates` lists the balanced columns and `balance_diagnostics` + # maps each sub-experiment a to {n_treated, n_control, effective_control_mass + # (Ñ^C_a), ess, max_imbalance_pre, max_imbalance_post, balance_solver}. When + # balancing, `stacked_data` carries `_b_sa` (raw design weights) and the + # `_Q_weight` column holds the composed final weights W_sa. + balance: str = "none" + covariates: Optional[List[str]] = None + balance_diagnostics: Optional[Dict[Any, Dict[str, Any]]] = field(default=None) # --- Inference-field aliases (balance/external-adapter compatibility) --- @property diff --git a/docs/api/stacked_did.rst b/docs/api/stacked_did.rst index 20c2da8c..5c572e46 100644 --- a/docs/api/stacked_did.rst +++ b/docs/api/stacked_did.rst @@ -125,5 +125,5 @@ Comparison with Other Staggered Estimators - Full stacked dataset accessible - Group-time effects accessible * - Covariates - - Not yet supported + - Entropy balancing via ``balance="entropy"`` + ``fit(covariates=...)`` (CBWSDID, Ustyuzhanin 2026); requires ``weighting="aggregate"`` + balanced windows, no ``survey_design`` - Supported (OR, IPW, DR) diff --git a/docs/doc-deps.yaml b/docs/doc-deps.yaml index f7720541..ae3616a1 100644 --- a/docs/doc-deps.yaml +++ b/docs/doc-deps.yaml @@ -529,6 +529,15 @@ sources: # ── StackedDiD (stacked_did group) ───────────────────────────────── + diff_diff/balancing.py: + drift_risk: medium + docs: + - path: docs/methodology/REGISTRY.md + section: "StackedDiD (Covariate balancing / CBWSDID)" + type: methodology + - path: docs/references.rst + type: user_guide + diff_diff/stacked_did.py: drift_risk: medium docs: diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index dfb0431b..9a9c65d2 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -1538,6 +1538,30 @@ The paper text states a stricter bound (T_min + 1) but the R code by the co-auth - [x] Survey design support (Phase 3): Q-weights compose multiplicatively with survey weights; TSL vcov on composed weights; survey design columns propagated through sub-experiments. Replicate weights supported via estimator-level refit with Q-weight composition (see Replicate Weight Variance section). - **Note:** Survey weights compose multiplicatively with Q-weights for StackedDiD; only `weight_type="pweight"` (default) is supported — `fweight` and `aweight` are rejected because Q-weight composition changes weight semantics (non-integer for fweight, non-inverse-variance for aweight) +### Covariate balancing (CBWSDID) + +**Primary source:** Ustyuzhanin, V. (2026). Covariate-Balanced Weighted Stacked Difference-in-Differences. arXiv:2604.02293v1. https://arxiv.org/abs/2604.02293 (in-repo paper review: `docs/methodology/papers/ustyuzhanin-2026-review.md`). Within-sub-experiment refinement via **entropy balancing** (Hainmueller 2012). + +Optional `balance="entropy"` (constructor) + `fit(..., covariates=[...])` adds a within-sub-experiment design stage that reweights the clean controls toward the treated cohort under **conditional** parallel trends, before the Wing et al. (2024) corrective aggregation. + +*Design weights (Ustyuzhanin 2026 §3.1):* +For each sub-experiment `a`, entropy balancing produces nonnegative control weights `b_{sa}` matching the treated covariate means (covariates read at the last pre-treatment period `t = a-1-anticipation`; treated keep `b=1`). The final stacked weights compose `b_{sa}` with the corrective factor via the **effective control mass** `Ñ^C_a = Σ_{s∈C_a} b_{sa}`: + + W_{sa} = b_{sa} · (N^D_a / N^D_Ω) / (Ñ^C_a / Ñ^C_Ω) for s ∈ C_a + W_{sa} = 1 for s ∈ D_a + +The pooled estimator is `DID^{CBWSDID}_e = Σ_a (N^D_a/N^D_Ω)(Δ̄^D_{a,e} − Δ̄^{C,b}_{a,e})`, recovered by the existing Q-weighted WLS when `W_{sa}` is injected. **Estimand preservation:** because only controls are reweighted (treated cohorts and their shares `N^D_a/N^D_Ω` are unchanged), the target remains the trimmed aggregate ATT `θ^e_κ` — the refinement changes only how untreated trends are estimated. At `b_{sa}=1` this reduces to the paper's **unit-count** weighted stacked DID, which equals `StackedDiD(weighting="aggregate")` on **balanced** event windows (where unit and observation counts coincide). Validated end-to-end by `tests/test_methodology_stacked_did.py::TestCBWSDIDCovariateBalance` (closed-form `DID^{CBWSDID}_e` anchor at 1e-8), `::TestCBWSDIDEffectiveMass` (effective-mass corrector is load-bearing vs a naive `b·Q` multiply), and `::TestCBWSDIDRParity` (cross-language parity vs the R `cbwsdid` package, `refinement.method="weightit"` / `method="ebal"`, at 1e-5 — golden in `benchmarks/data/cbwsdid_golden.json`, regenerate via `benchmarks/R/generate_cbwsdid_golden.R`). + +- **Note:** The effective-mass `W_{sa}` is computed directly from cohort unit-counts + `Ñ^C_a` (a naive `b_{sa}·Q_aggregate` multiply is NOT equivalent — it aggregates control means with weights ∝ `(N^D_a/N^D_Ω)(Ñ^C_a/N^C_a)` instead of the required `∝ (N^D_a/N^D_Ω)`, biased unless `b` is uniform). +- **Note:** Inference is conditional-on-the-estimated-weights cluster-robust (the existing `hc1`/`hc2_bm` path with `W_{sa}` as the WLS weights) — the paper's default. The paper's weight-re-estimating bootstrap is NOT implemented in v1 (deliberate scope; entropy balancing is smooth so the Abadie–Imbens (2008) nonsmooth-matching bootstrap caveat does not apply). `cluster` is orthogonal to `b_{sa}` (weights conditioned-on); default `unit` matches the paper. +- **Note:** v1 scope — only `balance="entropy"` with `weighting="aggregate"`. `balance` + `population`/`sample_share` and `balance` + `survey_design=` raise `NotImplementedError`; matching-based balancing and the repeated `0→1/1→0` episode extension are out of scope. + +*Covariate-balancing edge cases:* +- Infeasible cohort (treated covariate mean outside the clean-control hull → entropy balancing cannot match the moments): **fail-closed** `ValueError` naming the cohort and worst covariate — NOT silently dropped (dropping a cohort would shift the estimand to an overlap-trimmed ATT, Ustyuzhanin 2026 §3.1). +- Degenerate design weights (`Ñ^C_a → 0` / highly concentrated weights): low effective sample size → `UserWarning` with the per-cohort diagnostic. +- Missing pre-treatment row, or covariate absent / `balance`↔`covariates` mismatch: `ValueError` at `fit()`. +- **Ragged / unbalanced event windows** (a unit not observed at every event time in a sub-experiment): **fail-closed `ValueError`** — `balance="entropy"` requires balanced windows. The paper assumes balanced event windows; off them the unit-count corrector and the observation-count `aggregate` Q diverge (the count-convention is unresolved, deferred). `balance="none"` continues to support unbalanced panels via observation-count Q. + --- ## WooldridgeDiD (ETWFE) diff --git a/docs/references.rst b/docs/references.rst index 4ef9a432..847ee7d5 100644 --- a/docs/references.rst +++ b/docs/references.rst @@ -220,6 +220,14 @@ Multi-Period and Staggered Adoption - **Wing, C., Freedman, S. M., & Hollingsworth, A. (2024).** "Stacked Difference-in-Differences." *NBER Working Paper* 32054. https://www.nber.org/papers/w32054 +- **Ustyuzhanin, V. (2026).** "Covariate-Balanced Weighted Stacked Difference-in-Differences." *arXiv preprint* arXiv:2604.02293v1. https://arxiv.org/abs/2604.02293 + + Primary source for ``StackedDiD``'s covariate-balancing path (the ``balance="entropy"`` parameter). Adds a within-sub-experiment design stage — nonnegative control design weights ``b_{sa}`` — that composes with the Wing et al. (2024) corrective weights via the effective control mass into the final stacked weights ``W_{sa}``, so untreated counterfactual trends are estimated under *conditional* parallel trends while the treated-cohort weights (and hence the trimmed-aggregate-ATT estimand) are unchanged. + +- **Hainmueller, J. (2012).** "Entropy Balancing for Causal Effects: A Multivariate Reweighting Method to Produce Balanced Samples in Observational Studies." *Political Analysis*, 20(1), 25-46. https://doi.org/10.1093/pan/mpr025 + + Foundational source for entropy balancing. ``StackedDiD(balance="entropy")`` implements this within-sub-experiment reweighting: nonnegative control weights that exactly match the treated covariate means while minimizing the Kullback-Leibler divergence from uniform. + - **Chen, X., Sant'Anna, P. H. C., & Xie, H. (2025).** "Efficient Difference-in-Differences and Event Study Estimators." *arXiv preprint* arXiv:2506.17729v1. https://arxiv.org/abs/2506.17729v1 (also Cowles Foundation Discussion Paper No. 2470). Primary source for the optimal-weighting / PT-All / PT-Post efficient DiD implemented in our ``EfficientDiD`` class. diff --git a/tests/test_balancing.py b/tests/test_balancing.py new file mode 100644 index 00000000..a65c8a57 --- /dev/null +++ b/tests/test_balancing.py @@ -0,0 +1,83 @@ +"""Unit tests for the entropy-balancing solver (diff_diff/balancing.py). + +Entropy balancing (Hainmueller 2012) must achieve exact first-moment balance to the +target whenever the target is in the interior of the control covariate hull, and must +fail loudly (BalanceError) when it is not. +""" + +import numpy as np +import pytest + +from diff_diff.balancing import BalanceError, entropy_balance + + +def _treated_control(seed, n_control=200, n_treated=80, k=3, shift=0.6): + """Two groups whose covariate means differ (treated shifted), so balancing is + non-trivial but feasible (treated mean stays inside the control spread).""" + rng = np.random.default_rng(seed) + Xc = rng.normal(0.0, 1.0, size=(n_control, k)) + Xt = rng.normal(shift, 1.0, size=(n_treated, k)) + return Xc, Xt + + +class TestEntropyBalance: + def test_exact_first_moment_balance(self): + Xc, Xt = _treated_control(seed=1) + target = Xt.mean(axis=0) + # before balancing the control means differ from the target + assert np.max(np.abs(Xc.mean(axis=0) - target)) > 0.1 + w, info = entropy_balance(Xc, target, tol=1e-10) + assert info["converged"] + # weighted control means match the treated means to tolerance + np.testing.assert_allclose(w @ Xc, target, atol=1e-9) + # weights are a valid nonnegative distribution + assert np.all(w >= 0) + np.testing.assert_allclose(w.sum(), 1.0, atol=1e-12) + assert info["max_residual"] < 1e-10 + assert 0 < info["ess"] <= Xc.shape[0] + + def test_trivial_target_is_near_uniform(self): + Xc, _ = _treated_control(seed=2) + target = Xc.mean(axis=0) # already balanced + w, info = entropy_balance(Xc, target, tol=1e-12) + assert info["converged"] + np.testing.assert_allclose(w, np.full(Xc.shape[0], 1.0 / Xc.shape[0]), atol=1e-8) + + def test_base_weights_respected(self): + Xc, Xt = _treated_control(seed=3) + target = Xt.mean(axis=0) + rng = np.random.default_rng(99) + q = rng.uniform(0.5, 2.0, size=Xc.shape[0]) + w, info = entropy_balance(Xc, target, base_weights=q, tol=1e-10) + assert info["converged"] + np.testing.assert_allclose(w @ Xc, target, atol=1e-9) + np.testing.assert_allclose(w.sum(), 1.0, atol=1e-12) + + def test_collinear_covariates_ridge(self): + # duplicate a column -> singular Hessian; ridge/lstsq must still balance the + # identified moments without blowing up. + Xc, Xt = _treated_control(seed=4, k=2) + Xc = np.column_stack([Xc, Xc[:, 0]]) # 3rd col == 1st + Xt = np.column_stack([Xt, Xt[:, 0]]) + target = Xt.mean(axis=0) + w, info = entropy_balance(Xc, target, tol=1e-8) + assert info["converged"] + np.testing.assert_allclose(w @ Xc, target, atol=1e-7) + + def test_infeasible_target_raises(self): + Xc, _ = _treated_control(seed=5, k=2) + # target far outside the control hull on covariate 0 -> no finite balancing λ + target = np.array([Xc[:, 0].max() + 5.0, 0.0]) + with pytest.raises(BalanceError) as exc: + entropy_balance(Xc, target, tol=1e-8, max_iter=100) + assert exc.value.max_residual >= 1e-8 + assert exc.value.residuals.shape == (2,) + + def test_input_validation(self): + Xc, _ = _treated_control(seed=6, k=3) + with pytest.raises(ValueError): + entropy_balance(Xc, np.zeros(2)) # wrong target length + with pytest.raises(ValueError): + entropy_balance(Xc, np.zeros(3), base_weights=np.ones(5)) # wrong q length + with pytest.raises(ValueError): + entropy_balance(Xc[:, 0], np.zeros(3)) # 1-D X diff --git a/tests/test_guides.py b/tests/test_guides.py index b2ab8f68..fadc85fc 100644 --- a/tests/test_guides.py +++ b/tests/test_guides.py @@ -719,3 +719,54 @@ def test_llms_full_had_pretests_assumption_labels_correct(self): f"Assumption 7 is pre-trends (step 2, only " f"covered on the event-study path). Line: {line!r}" ) + + +class TestLLMsFullStackedDiDCoverage: + """Pin the StackedDiD section of llms-full.txt to the real API. + + Adding a public parameter (here: balance= on __init__, covariates= on fit()) + requires updating diff_diff/guides/llms-full.txt — these tests catch drift. + """ + + def _stacked_section(self): + text = get_llm_guide("full") + start = text.index("### StackedDiD") + nxt = text.index("\n### ", start + 1) + return text[start:nxt] + + def test_llms_full_has_stacked_section(self): + assert "### StackedDiD" in get_llm_guide("full") + + def test_llms_full_stacked_constructor_signature_matches_real_api(self): + import inspect + + from diff_diff import StackedDiD + + sig_params = set(inspect.signature(StackedDiD.__init__).parameters) + sig_params.discard("self") + section = self._stacked_section() + block_start = section.index("StackedDiD(") + block_end = section.index("\n)", block_start) + ctor_block = section[block_start:block_end] + for param in sig_params: + assert f"{param}:" in ctor_block or f"{param} " in ctor_block, ( + f"StackedDiD constructor block in llms-full.txt is missing the real " + f"public parameter {param!r} (adding a public param requires updating " + f"the guide)." + ) + + def test_llms_full_stacked_fit_documents_covariates(self): + import inspect + + from diff_diff import StackedDiD + + assert "covariates" in inspect.signature(StackedDiD.fit).parameters + section = self._stacked_section() + fit_start = section.index("stacked.fit(") + fit_block = section[fit_start : section.index("\n)", fit_start)] + assert "covariates" in fit_block, ( + "StackedDiD.fit() exposes covariates= but the llms-full.txt fit() block " + "does not document it." + ) + # balance= must be documented somewhere in the section (constructor param) + assert "balance" in section diff --git a/tests/test_methodology_stacked_did.py b/tests/test_methodology_stacked_did.py index 641d14a3..3257ef1a 100644 --- a/tests/test_methodology_stacked_did.py +++ b/tests/test_methodology_stacked_did.py @@ -27,11 +27,13 @@ import os import numpy as np +import pandas as pd import pytest from scipy.optimize import brentq from scipy.stats import t as t_dist from diff_diff import StackedDiD +from diff_diff.balancing import entropy_balance from diff_diff.prep_dgp import generate_staggered_data _GOLDEN_PATH = os.path.join( @@ -667,3 +669,348 @@ def test_sample_share_weighting_hc2_bm_matches_clubsandwich_cr2(self, goldens, p rtol=1e-6, err_msg="sample_share hc2_bm overall BM DOF must match R Wald_test(HTZ)", ) + + +# ============================================================================= +# CBWSDID covariate balancing (Ustyuzhanin 2026) — estimand validation. +# Self-contained (no R golden needed): the closed-form paper formula is the +# primary, R-independent estimand anchor (productionized validation spike). +# ============================================================================= + +_CB_KP = 2 +_CB_KPOST = 2 + + +def _cbwsdid_panel( + seed, + *, + hetero=False, + n_per_cohort=25, + n_never=100, + cohorts=(4, 5, 6), + n_periods=10, + trend_gamma=0.5, + te0=2.0, + te_gamma=1.0, + noise=0.03, +): + """Staggered panel where untreated trends depend (linearly) on covariate x and + treated cohorts are selected on x (so UNCONDITIONAL parallel trends fails). With + hetero=True the effect is heterogeneous in x (tau_i = te0 + te_gamma*x_i).""" + rng = np.random.default_rng(seed) + rows = [] + uid = 0 + specs = [(c, n_per_cohort) for c in cohorts] + [(np.inf, n_never)] + for ft, n in specs: + x_mean = 0.8 if np.isfinite(ft) else 0.0 # selection on x + for _ in range(n): + uid += 1 + x = rng.normal(x_mean, 1.0) + alpha = rng.normal() + slope = trend_gamma * x + rng.normal(0, 0.02) # x-driven untreated trend + for t in range(1, n_periods + 1): + y0 = alpha + slope * t + rng.normal(0, noise) + te = ( + (te0 + (te_gamma * x if hetero else 0.0)) + if (np.isfinite(ft) and t >= ft) + else 0.0 + ) + rows.append((uid, t, ft if np.isfinite(ft) else 0, x, y0 + te)) + return pd.DataFrame(rows, columns=["unit", "time", "first_treat", "x", "y"]) + + +def _closed_form_cbwsdid(df, omega, kp, kpost): + """DID^CBWSDID_e = Σ_a (N^D_a/N^D_Ω)(Δ̄^D_{a,e} − Δ̄^{C,b}_{a,e}), computed + independently with never-treated controls and entropy-balanced b_sa.""" + ftmap = df.drop_duplicates("unit").set_index("unit")["first_treat"].replace(0, np.inf) + ylook = df.set_index(["unit", "time"])["y"].to_dict() + xlook = df.drop_duplicates("unit").set_index("unit")["x"].to_dict() + never = ftmap[~np.isfinite(ftmap)].index.tolist() + event_times = [h for h in range(-kp, kpost + 1) if h != -1] + omega = sorted(omega) + ND = {a: int((ftmap == a).sum()) for a in omega} + ND_Om = sum(ND.values()) + out = {} + for h in event_times: + acc = 0.0 + for a in omega: + treated = ftmap[ftmap == a].index.tolist() + Xc = np.array([[xlook[u]] for u in never]) + target = np.array([np.mean([xlook[u] for u in treated])]) + # Use the SAME tol the estimator uses internally so the entropy weights + # are identical (the solver stops at the same iterate) — then the closed + # form matches the fitted regression to WLS-algebra precision. + b, _ = entropy_balance(Xc, target, tol=1e-8) + mt = float(np.mean([ylook[(u, a + h)] - ylook[(u, a - 1)] for u in treated])) + dyc = np.array([ylook[(u, a + h)] - ylook[(u, a - 1)] for u in never]) + acc += (ND[a] / ND_Om) * (mt - float(b @ dyc)) + out[h] = acc + return out + + +class TestCBWSDIDCovariateBalance: + """CBWSDID estimand validation (Ustyuzhanin 2026).""" + + def _fit(self, df, balance="entropy"): + est = StackedDiD( + kappa_pre=_CB_KP, + kappa_post=_CB_KPOST, + weighting="aggregate", + clean_control="never_treated", + cluster="unit", + balance=balance, + ) + return est.fit( + df, + outcome="y", + unit="unit", + time="time", + first_treat="first_treat", + aggregate="event_study", + covariates=["x"] if balance != "none" else None, + ) + + def test_matches_closed_form_paper_formula(self): + # PRIMARY anchor: fitted CBWSDID event-study coefs == the closed-form paper + # formula (Ustyuzhanin 2026 §3.1) to machine precision. Pins the effective- + # mass W_sa composition end-to-end without R. + df = _cbwsdid_panel(seed=11) + res = self._fit(df) + fitted = { + h: res.event_study_effects[h]["effect"] + for h in res.event_study_effects + if res.event_study_effects[h]["n_obs"] > 0 + } + closed = _closed_form_cbwsdid(df, res.groups, _CB_KP, _CB_KPOST) + assert set(fitted) == {h for h in closed} + for h in fitted: + assert abs(fitted[h] - closed[h]) < 1e-8, (h, fitted[h], closed[h]) + + def test_estimand_integrity_under_heterogeneity(self): + # The failure mode of the abandoned regression-adjustment attempt: with + # heterogeneous effects (tau_i = te0 + te_gamma*x_i) AND x-confounded untreated + # trends, control REWEIGHTING (CBWSDID) recovers the treated-average ATT while + # plain StackedDiD is biased. Assert BOTH directions so the test proves the + # mechanism, not just no-exception. + df = _cbwsdid_panel(seed=7, hetero=True, te0=2.0, te_gamma=1.0) + units = df.drop_duplicates("unit").copy() + units["ft"] = units["first_treat"].replace(0, np.inf) + treated_x_mean = units.loc[np.isfinite(units["ft"]), "x"].mean() + true_att = 2.0 + 1.0 * treated_x_mean # treated-average ATT + + cb = self._fit(df, balance="entropy").overall_att + plain = self._fit(df, balance="none").overall_att + + assert abs(cb - true_att) < 0.15, f"CBWSDID should recover {true_att:.3f}, got {cb:.3f}" + assert abs(plain - true_att) > 0.3, f"plain StackedDiD should be biased, got {plain:.3f}" + assert abs(cb - true_att) < abs(plain - true_att) + + def test_balanced_covariate_reduces_to_plain_stacked(self): + # When the covariate is independent of treatment (already balanced), entropy + # weights are ~uniform and CBWSDID ≈ plain weighted stacked DID. + rng = np.random.default_rng(3) + rows = [] + uid = 0 + for ft, n in [(4, 25), (5, 25), (6, 25), (np.inf, 100)]: + for _ in range(n): + uid += 1 + x = rng.normal() # SAME distribution for treated and control + alpha = rng.normal() + for t in range(1, 11): + y = ( + alpha + + 0.1 * t + + rng.normal(0, 0.05) + + (2.0 if (np.isfinite(ft) and t >= ft) else 0.0) + ) + rows.append((uid, t, ft if np.isfinite(ft) else 0, x, y)) + df = pd.DataFrame(rows, columns=["unit", "time", "first_treat", "x", "y"]) + cb = self._fit(df, balance="entropy").overall_att + plain = self._fit(df, balance="none").overall_att + assert abs(cb - plain) < 0.05 + assert abs(cb - 2.0) < 0.1 + + +def _spread_panel(seed=21, n_per=15, n_never=50): + """Cohorts spread out so not-yet-treated control SETS (and counts N^C_a) DIFFER + across sub-experiments — required to distinguish the effective-mass corrective + from a naive b*Q multiply (they coincide when every cohort shares one control set).""" + rng = np.random.default_rng(seed) + rows = [] + uid = 0 + for ft, n in [(4, n_per), (8, n_per), (12, n_per), (np.inf, n_never)]: + treated = np.isfinite(ft) + xm = 0.6 if treated else 0.0 + for _ in range(n): + uid += 1 + x = rng.normal(xm, 1.0) + alpha = rng.normal() + slope = 0.3 * x + rng.normal(0, 0.02) + for t in range(1, 17): + y = ( + alpha + + slope * t + + rng.normal(0, 0.03) + + (2.0 if (treated and t >= ft) else 0.0) + ) + rows.append((uid, t, ft if treated else 0, x, y)) + return pd.DataFrame(rows, columns=["unit", "time", "first_treat", "x", "y"]) + + +def _closed_form_from_result(df, res, kp, kpost, *, naive=False): + """Closed-form CBWSDID using the estimator's OWN per-sub-experiment control sets + (robust to clean_control mode) and independently-solved entropy weights. With + naive=True, aggregate the control means with the WRONG raw-count corrective + weights ∝ (N^D_a/N^D_Ω)(Ñ^C_a/N^C_a) instead of the effective-mass ∝ (N^D_a/N^D_Ω).""" + sd = res.stacked_data + ylook = df.set_index(["unit", "time"])["y"].to_dict() + xlook = df.drop_duplicates("unit").set_index("unit")["x"].to_dict() + event_times = [h for h in range(-kp, kpost + 1) if h != -1] + omega = sorted(res.groups) + members = {} + ND = {} + NC = {} + NtC = {} + bmap = {} + for a in omega: + sub = sd[sd["_sub_exp"] == a] + treated = list(pd.unique(sub.loc[sub["_D_sa"] == 1, "unit"])) + controls = list(pd.unique(sub.loc[sub["_D_sa"] == 0, "unit"])) + members[a] = (treated, controls) + ND[a] = len(treated) + NC[a] = len(controls) + Xc = np.array([[xlook[u]] for u in controls]) + target = np.array([np.mean([xlook[u] for u in treated])]) + b, _ = entropy_balance(Xc, target, tol=1e-8) + bmap[a] = b + NtC[a] = float(np.sum(b)) + ND_Om = sum(ND.values()) + out = {} + for h in event_times: + treated_side = 0.0 + # control side aggregated with effective-mass (correct) or raw-count (naive) weights + cnum = 0.0 + cden = 0.0 + for a in omega: + treated, controls = members[a] + mt = float(np.mean([ylook[(u, a + h)] - ylook[(u, a - 1)] for u in treated])) + dyc = np.array([ylook[(u, a + h)] - ylook[(u, a - 1)] for u in controls]) + mcb = float(bmap[a] @ dyc) + share = ND[a] / ND_Om + treated_side += share * mt + w_a = share * (NtC[a] / NC[a]) if naive else share + cnum += w_a * mcb + cden += w_a + out[h] = treated_side - cnum / cden + return out + + +class TestCBWSDIDEffectiveMass: + """The effective-mass corrective is load-bearing (vs a naive b*Q multiply).""" + + KP = KPOST = 2 + + def _fit(self, df): + est = StackedDiD( + kappa_pre=self.KP, + kappa_post=self.KPOST, + weighting="aggregate", + clean_control="not_yet_treated", + cluster="unit", + balance="entropy", + ) + return est.fit( + df, + outcome="y", + unit="unit", + time="time", + first_treat="first_treat", + aggregate="event_study", + covariates=["x"], + ) + + def test_anchor_and_effective_mass_is_load_bearing(self): + df = _spread_panel() + res = self._fit(df) + fitted = { + h: res.event_study_effects[h]["effect"] + for h in res.event_study_effects + if res.event_study_effects[h]["n_obs"] > 0 + } + # control counts genuinely vary across cohorts (else the test is vacuous) + ncs = {a: d["n_control"] for a, d in res.balance_diagnostics.items()} + assert len(set(ncs.values())) > 1, f"need varying N^C_a, got {ncs}" + + effective = _closed_form_from_result(df, res, self.KP, self.KPOST, naive=False) + naive = _closed_form_from_result(df, res, self.KP, self.KPOST, naive=True) + for h in fitted: + # (1) fitted == effective-mass closed form (discriminating anchor) + assert abs(fitted[h] - effective[h]) < 1e-8, (h, fitted[h], effective[h]) + # (2) effective-mass differs MATERIALLY from the naive b*Q aggregation + max_gap = max(abs(effective[h] - naive[h]) for h in fitted) + assert max_gap > 1e-3, f"effective-mass vs naive gap too small ({max_gap:.2e})" + # (3) the estimator did NOT implement the naive form + for h in fitted: + assert abs(fitted[h] - naive[h]) > 1e-4 + + +_CBWSDID_GOLDEN = os.path.join( + os.path.dirname(__file__), "..", "benchmarks", "data", "cbwsdid_golden.json" +) +_CBWSDID_PANEL = os.path.join( + os.path.dirname(__file__), "..", "benchmarks", "data", "cbwsdid_balance_panel.csv" +) + + +@pytest.mark.skipif( + not os.path.exists(_CBWSDID_GOLDEN), + reason="cbwsdid_golden.json not generated; run Rscript benchmarks/R/generate_cbwsdid_golden.R", +) +class TestCBWSDIDRParity: + """Cross-language parity against the reference R package `cbwsdid` (Ustyuzhanin 2026). + + The committed golden holds dynamic event-study ATTs from + ``cbwsdid(design='absorbing', refinement.method='weightit', method='ebal', + covs.formula=~x)`` on ``benchmarks/data/cbwsdid_balance_panel.csv``. + ``StackedDiD(balance='entropy', covariates=['x'])`` — with an INDEPENDENT + entropy-balancing solver (custom Newton vs WeightIt's ebal) and the effective-mass + ``W_sa`` composition — reproduces them to ~3e-7 (the residual is the cross-solver + balancing tolerance); asserted at 1e-5. Regenerate the golden with + ``Rscript benchmarks/R/generate_cbwsdid_golden.R`` (requires + ``remotes::install_github('vadvu/cbwsdid')``). + """ + + def test_dynamic_atts_match_r_cbwsdid(self): + with open(_CBWSDID_GOLDEN) as f: + g = json.load(f) + df = pd.read_csv(_CBWSDID_PANEL) + res = StackedDiD( + kappa_pre=2, + kappa_post=2, + weighting="aggregate", + clean_control="not_yet_treated", + cluster="unit", + balance="entropy", + ).fit( + df, + outcome="y", + unit="unit", + time="time", + first_treat="first_treat", + aggregate="event_study", + covariates=["x"], + ) + for et, r_est, r_se in zip( + g["dynamic"]["event_time"], g["dynamic"]["estimate"], g["dynamic"]["std_error"] + ): + py = res.event_study_effects[et]["effect"] + assert abs(py - r_est) < 1e-5, f"event {et}: python {py} vs R cbwsdid {r_est}" + # SEs also agree (conditional-on-weights cluster-robust); the ~0.2-0.5% + # gap is the small-sample correction (diff-diff CR1S vs cbwsdid's vcov). + py_se = res.event_study_effects[et]["se"] + np.testing.assert_allclose( + py_se, + r_se, + rtol=0.02, + err_msg=f"event {et}: python SE {py_se} vs R cbwsdid SE {r_se}", + ) diff --git a/tests/test_stacked_did.py b/tests/test_stacked_did.py index 7fdf9925..348df20f 100644 --- a/tests/test_stacked_did.py +++ b/tests/test_stacked_did.py @@ -1715,3 +1715,224 @@ def test_summary_renders_clustered_variance_label(self, staggered_data): assert any( "at unit_subexp" in line for line in variance_lines ), f"cluster='unit_subexp' should render in label: {variance_lines}" + + +# ============================================================================= +# Covariate balancing (CBWSDID, Ustyuzhanin 2026) +# ============================================================================= + + +def _balance_panel( + seed=5, + *, + hetero_x=True, + n_per=20, + n_never=60, + const_x=False, + treated_x_mean=0.7, + infeasible=False, +): + """Staggered panel with a covariate x; treated cohorts optionally selected on x.""" + rng = np.random.default_rng(seed) + rows = [] + uid = 0 + for ft, n in [(4, n_per), (5, n_per), (6, n_per), (np.inf, n_never)]: + treated = np.isfinite(ft) + if const_x: + xm, xs = 1.0, 0.0 + elif infeasible and treated: + xm, xs = 12.0, 0.1 # treated x far outside the control hull + else: + xm, xs = (treated_x_mean if (treated and hetero_x) else 0.0), 1.0 + for _ in range(n): + uid += 1 + x = xm + xs * rng.normal() + alpha = rng.normal() + slope = 0.3 * x + rng.normal(0, 0.02) + for t in range(1, 11): + y = ( + alpha + + slope * t + + rng.normal(0, 0.04) + + (2.0 if (treated and t >= ft) else 0.0) + ) + rows.append((uid, t, ft if treated else 0, x, y)) + return pd.DataFrame(rows, columns=["unit", "time", "first_treat", "x", "y"]) + + +def _cb_fit(df, **kw): + params = dict( + kappa_pre=2, + kappa_post=2, + weighting="aggregate", + clean_control="never_treated", + cluster="unit", + ) + params.update( + {k: v for k, v in kw.items() if k in ("weighting", "balance", "cluster", "vcov_type")} + ) + est = StackedDiD(**params) + return est.fit( + df, + outcome="y", + unit="unit", + time="time", + first_treat="first_treat", + aggregate=kw.get("aggregate", "event_study"), + covariates=kw.get("covariates"), + survey_design=kw.get("survey_design"), + population=kw.get("population"), + ) + + +class TestStackedDiDCovariateBalance: + """CBWSDID covariate-balancing path on StackedDiD.""" + + # ---- validation + scope guards ---- + def test_invalid_balance_value_raises(self): + with pytest.raises(ValueError, match="balance must be"): + StackedDiD(balance="match") + + def test_set_params_revalidates_balance(self): + est = StackedDiD() + with pytest.raises(ValueError, match="balance must be"): + est.set_params(balance="ipw") + + def test_balance_without_covariates_raises(self): + df = _balance_panel() + with pytest.raises(ValueError, match="requires a non-empty covariates"): + _cb_fit(df, balance="entropy", covariates=None) + + def test_covariates_without_balance_raises(self): + df = _balance_panel() + with pytest.raises(ValueError, match="balance='none'"): + _cb_fit(df, balance="none", covariates=["x"]) + + def test_balance_with_population_weighting_raises(self): + df = _balance_panel() + df["pop"] = 1.0 + with pytest.raises(NotImplementedError, match="weighting='aggregate'"): + _cb_fit( + df, balance="entropy", covariates=["x"], weighting="population", population="pop" + ) + + def test_missing_covariate_column_raises(self): + df = _balance_panel() + with pytest.raises(ValueError, match="not found in data"): + _cb_fit(df, balance="entropy", covariates=["nonexistent"]) + + def test_infeasible_cohort_raises_named(self): + df = _balance_panel(infeasible=True) + with pytest.raises(ValueError, match="balancing failed for cohort"): + _cb_fit(df, balance="entropy", covariates=["x"]) + + # ---- correctness / properties ---- + def test_balance_constraint_satisfied(self): + df = _balance_panel(seed=2) + res = _cb_fit(df, balance="entropy", covariates=["x"]) + assert res.balance == "entropy" + assert res.covariates == ["x"] + assert res.balance_diagnostics is not None + for a, d in res.balance_diagnostics.items(): + # entropy balancing achieves exact first-moment balance + assert d["max_imbalance_post"] < 1e-7 + assert d["max_imbalance_pre"] > d["max_imbalance_post"] + # stacked_data carries the raw design weights + assert "_b_sa" in res.stacked_data.columns + + def test_constant_covariate_reduces_to_plain_stacked(self): + # constant covariate -> uniform entropy weights -> identical to plain StackedDiD + df = _balance_panel(seed=4, const_x=True) + cb = _cb_fit(df, balance="entropy", covariates=["x"]) + plain = _cb_fit(df, balance="none") + for h in plain.event_study_effects: + if plain.event_study_effects[h]["n_obs"] > 0: + assert ( + abs( + cb.event_study_effects[h]["effect"] - plain.event_study_effects[h]["effect"] + ) + < 1e-6 + ) + + def test_covariate_scale_invariance(self): + # entropy balancing is invariant to covariate rescaling: x vs 2*x give the + # same weights and hence the same estimate. + df = _balance_panel(seed=6) + df["x2"] = 2.0 * df["x"] + a = _cb_fit(df, balance="entropy", covariates=["x"]) + b = _cb_fit(df, balance="entropy", covariates=["x2"]) + assert abs(a.overall_att - b.overall_att) < 1e-8 + + def test_balance_on_default_aggregate_mode(self): + # balancing must work on the overall-ATT (aggregate=None) path, not only + # event_study (the BM-DOF machinery branches on aggregate). + df = _balance_panel(seed=8) + res = _cb_fit(df, balance="entropy", covariates=["x"], aggregate=None) + assert np.isfinite(res.overall_att) + assert res.balance == "entropy" + + # ---- sklearn surface ---- + def test_get_params_includes_balance_and_clone(self): + est = StackedDiD(balance="entropy") + assert est.get_params()["balance"] == "entropy" + clone = StackedDiD(**est.get_params()) # round-trip + assert clone.balance == "entropy" + + def test_convenience_function_threads_covariates(self): + df = _balance_panel(seed=9) + res = stacked_did( + df, + "y", + "unit", + "time", + "first_treat", + kappa_pre=2, + kappa_post=2, + aggregate="event_study", + clean_control="never_treated", + balance="entropy", + covariates=["x"], + ) + assert res.balance == "entropy" + assert res.covariates == ["x"] + + def test_ragged_event_window_raises(self): + # balance="entropy" requires balanced event windows: on a ragged panel (a unit + # missing a non-reference event-time row) the unit-count corrector and the + # observation-count aggregate Q diverge, so fail closed. balance="none" still + # supports the unbalanced panel. + df = _balance_panel(seed=10, const_x=True) + u = int(df.loc[df["first_treat"] == 4, "unit"].iloc[0]) + drop = (df["unit"] == u) & (df["time"] == 5) # a post-period row (keeps t=a-1 ref) + dfr = df[~drop].reset_index(drop=True) + with pytest.raises(ValueError, match="balanced event windows"): + _cb_fit(dfr, balance="entropy", covariates=["x"]) + assert np.isfinite(_cb_fit(dfr, balance="none").overall_att) + + def test_zero_row_eligible_control_raises(self): + # An eligible (never-treated) control with ZERO rows in a cohort's window is + # silently dropped by _build_sub_experiment — the count-only guard misses it, + # the exact-coverage guard catches it (missing_units). + df = _balance_panel(seed=11, const_x=True) + u = int(df.loc[df["first_treat"] == 0, "unit"].iloc[0]) # a never-treated unit + drop = (df["unit"] == u) & (df["time"].between(2, 6)) # cohort-4 window [2,6] + dfz = df[~drop].reset_index(drop=True) + with pytest.raises(ValueError, match="balanced event windows"): + _cb_fit(dfz, balance="entropy", covariates=["x"]) + + def test_duplicate_and_missing_event_time_raises(self): + # A unit with a duplicated event time AND a missing one nets the expected row + # count, bypassing a count-only guard; the exact-coverage guard catches the + # duplicate (unit, event_time). + df = _balance_panel(seed=12, const_x=True) + u = int(df.loc[df["first_treat"] == 4, "unit"].iloc[0]) + dupe = df[(df["unit"] == u) & (df["time"] == 3)] # event_time -1 for cohort 4 + df2 = df[~((df["unit"] == u) & (df["time"] == 5))] # drop event_time +1 + dfd = pd.concat([df2, dupe], ignore_index=True) # same row count, ragged coverage + with pytest.raises(ValueError, match="balanced event windows"): + _cb_fit(dfd, balance="entropy", covariates=["x"]) + + def test_string_covariate_raises_typeerror(self): + df = _balance_panel(seed=13, const_x=True) + with pytest.raises(TypeError, match="must be a list"): + _cb_fit(df, balance="entropy", covariates="x")