diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 0000000..82a1564 --- /dev/null +++ b/.gitattributes @@ -0,0 +1,2 @@ +# Preserve the authors' BOM and CRLF in the frozen state-list input. +scspill/data/files/cao_dowd/state_list.csv whitespace=cr-at-eol diff --git a/CHANGELOG.md b/CHANGELOG.md index 9393c03..300942a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,7 +2,29 @@ ## Unreleased -Authorship and citation metadata only; no functional change. +### New model: `cd` (Cao & Dowd) + +- **`method="cd"` implements Cao and Dowd (arXiv:1902.07343v4)**: an + intercept-adjusted simplex SCM for every unit, a fixed exposure structure + given by `affected_units` or a labeled `effect_matrix`, period-specific + direct and spillover effects, Andrews-type P tests for `C α_t = d`, and the + `κ_A` specification test. `load_california_cd(51)` and + `load_california_cd(41)` bundle the authors' 50-state-plus-DC panel. +- **The defaults follow the authors' Proposition 99 replication code**, so a + default fit reproduces the paper's Section 6 tables and Figure 4: + interpolated critical values (`quantile_method="linear"`) and intervals + that add the signed pre-treatment error quantiles to the estimate + (`ci_method="published"`, returned by `effect_ci`). The paper's Section 4 + definitions are available as `quantile_method="empirical"`, + `ci_method="inversion"` and `reference_mode="leave_one_period_out"`; the + differences are listed under Caveats in `docs/models/cd.qmd`. +- `scspill/tests/test_cd.py` checks the 51- and 41-unit fits against the + authors' frozen output (`jcao0/synthetic-control-spillover`, commit + `86004020`), stored in `benchmarks/reference/cao_dowd/`. + +### Authorship and citation + +No functional change. - **The software has three authors and one citation.** Shosei Sakaguchi and Hayato Tagawa have joined the project as co-authors, so `CITATION.cff`, the diff --git a/README.md b/README.md index 637c367..f5a98da 100644 --- a/README.md +++ b/README.md @@ -24,18 +24,19 @@ express at all. | Model | `method` | Reference | Class | Docs | |---|---|---|---|---| | Bayesian spatial-autoregressive spillover SCM | `"sar"` | [Sakaguchi & Tagawa (2026)](https://doi.org/10.1093/ectj/utag006), *The Econometrics Journal* | `SCSPILL` | [`sar`](https://quarcs-lab.github.io/scspill/models/sar.html) | +| SCM with spillover effects | `"cd"` | [Cao & Dowd, v4](https://arxiv.org/html/1902.07343v4) | `SCSPILL` | [`cd`](https://quarcs-lab.github.io/scspill/models/cd.html) | -One model ships today. `sar` routes spillovers through spatial weights you +`sar` routes spillovers through spatial weights you supply, scaled by a single intensity ρ that is estimated rather than assumed, fits unconstrained horseshoe-shrunk synthetic weights, and returns full Bayesian uncertainty for both estimands. At ρ = 0 it collapses exactly to the Bayesian horseshoe synthetic control. -Three further spillover-aware models — Cao & Dowd, the inclusive SCM of -Di Stefano & Mellace, and the partial-interference SCG of Grossi et al. — are -on the [roadmap](https://quarcs-lab.github.io/scspill/models/#planned). They -are **not implemented**, and `SCSPILLConfig` rejects their names rather than -falling back silently. +`cd` uses intercept-adjusted SCM equations for every unit and a specified +exposure structure to recover direct and spillover effects. It includes +Andrews P tests and a structure test. The inclusive SCM of Di Stefano & +Mellace and partial-interference SCG of Grossi et al. remain on the +[roadmap](https://quarcs-lab.github.io/scspill/models/#planned). The estimator architecture follows [mlsynth](https://github.com/jgreathouse9/mlsynth) (a pydantic config in, a @@ -76,10 +77,10 @@ result.plot(kind="panel") # counterfactual | effect | top spillovers | Subpackage | What it does | Docs | |---|---|---| -| `scspill` | The model layer. `SCSPILL(config).fit()` — today the `sar` model's two-step Bayesian sampler (horseshoe synthetic weights, SAR spillover block, adaptive Metropolis for the spillover intensity) and its identification formulas | [Get started](https://quarcs-lab.github.io/scspill/get-started.html) | +| `scspill` | The model layer. `SCSPILL(config).fit()` dispatches the Bayesian `sar` and Cao–Dowd `cd` estimators | [Models](https://quarcs-lab.github.io/scspill/models/) | | `scspill.validation` | `sar`'s sampler validation: the Geweke (2004) joint distribution test, prior-sensitivity grids, prior predictive checks | [Validation](https://quarcs-lab.github.io/scspill/articles/validation.html) | | `scspill.simulate` | `sar`'s Monte Carlo engine: rook-lattice SAR DGP, SCM/BSCM/SCSPILL comparison, the Tables 1–2 grid | [Simulation study](https://quarcs-lab.github.io/scspill/articles/simulation-study.html) | -| `scspill.data` | The bundled California Prop 99 and Sudan secession spillover panels — model-agnostic, usable by any model added later | [Datasets](https://quarcs-lab.github.io/scspill/articles/datasets.html) | +| `scspill.data` | California and Sudan panels, plus the authors' 51/41-unit Cao–Dowd Proposition 99 data | [Datasets](https://quarcs-lab.github.io/scspill/articles/datasets.html) | ## `sar` is validated, not just implemented diff --git a/benchmarks/reference/cao_dowd/README.md b/benchmarks/reference/cao_dowd/README.md new file mode 100644 index 0000000..4b4e231 --- /dev/null +++ b/benchmarks/reference/cao_dowd/README.md @@ -0,0 +1,41 @@ +# Cao–Dowd reference inputs + +Source: [jcao0/synthetic-control-spillover](https://github.com/jcao0/synthetic-control-spillover), +commit `86004020a248c7ee18f3c8a68b1e41db6abc0dda`, for +[Cao and Dowd, arXiv:1902.07343v4](https://arxiv.org/html/1902.07343v4). + +- `scspill/data/files/cao_dowd/cigs_consumption.csv` and `state_list.csv` are + copied from `replication files/application/prop99 spillover 51 units/data/`. + These are the CDC-based public-domain input data: 50 states plus DC, + 1970–2000, with treatment beginning in 1989. They are distinct from the + package's older 39-state, 1988-onward California panel. +- `results_51.json` and `results_41.json` are the authors' frozen application + outputs in the respective `prop99 spillover * units/output/` directories. +- `inference_replication_1.npz` contains the six first-replication outcome + panels generated by the authors' inference DGP and the corresponding + SP-A direct-effect estimates at effect 5 from + `inference/provenance/published_estimates.npz`. The settings are stationary + and integrated factors, each at pure-donor distance $r=0,0.5,1$. +- `loo_reference_stationary_r05_rep1.npy` holds all 50 squared SP-A + leave-one-pre-period-out reference statistics for the stationary, $r=0.5$ + first inference replication. They were independently generated by the + authors' `inference/src/procedures.py::andrews_loo_references` at the same + pinned commit. The slow estimator test compares all 50 statistics. +- `kappa_rep1.json` records the author's range and homogeneous-structure + $\kappa_A$ statistics and all 50 full-sample reference values on the + stationary, $r=0.5$ first inference panel. The statistic was generated with + `kappa_a/src/structures.py` and the fitted SCM equations with + `inference/src/scm.py` at the pinned commit. + +The bundled application CSV SHA-256 checksums are +`2b96994df9f28f7a28948fe765a20efa65270cb36662f44422899bbf98119ddc` +and `8e35c51b30e233bb8059be282d22b7f72637dc2950d7275158f1be26bd4bece8`, +respectively. The frozen output checksums are +`06357ece7e87a6b981e23bbe6338219849259b09303f917544d0f718a5eb55aa` +for 51 units, `0cf41c805905c6c989a66262a721800d674632ec958acc7f835be76f78aaa1f3` +for 41 units. The derived simulation fixture has SHA-256 +`7634351bfb1ea95242b74a89dfe8d10d865f48f72572ddbde50a4944cdec88de`. +The leave-one-out reference vector has SHA-256 +`b07e8cf0ed30fbdbc8d654ac8e2e11819239547cb8ff2c664f7ec4345a1d7e1a`. +The $\kappa_A$ reference JSON has SHA-256 +`267c7b9fd9e09e0c458ffcad78044468e99c36d5c8bec86d21af358ae15f9a45`. diff --git a/benchmarks/reference/cao_dowd/inference_replication_1.npz b/benchmarks/reference/cao_dowd/inference_replication_1.npz new file mode 100644 index 0000000..dbd1e5b Binary files /dev/null and b/benchmarks/reference/cao_dowd/inference_replication_1.npz differ diff --git a/benchmarks/reference/cao_dowd/kappa_rep1.json b/benchmarks/reference/cao_dowd/kappa_rep1.json new file mode 100644 index 0000000..956bce6 --- /dev/null +++ b/benchmarks/reference/cao_dowd/kappa_rep1.json @@ -0,0 +1,112 @@ +{ + "range": { + "statistic": 5.866044590640894, + "reference": [ + 6.73742584461722, + 5.32690209764272, + 5.598498441459483, + 5.623980426657359, + 5.192919548247222, + 4.682020135408817, + 5.364715300901857, + 6.183503725580241, + 5.826501334400295, + 5.881860095929147, + 5.439145864861086, + 5.969207242167086, + 5.18314959713628, + 5.600658356411405, + 4.500061965114871, + 5.6278601296189175, + 4.530813003116822, + 4.915967585201005, + 5.281668471926136, + 4.868597638775217, + 5.069550240358595, + 5.37807597568166, + 5.582020615270328, + 5.106798781334146, + 6.1078932913427995, + 5.579041201566973, + 5.061085601965185, + 5.754565605202067, + 5.185274584811025, + 5.7658920860929355, + 4.982686537801549, + 5.13090731067515, + 4.660774167963435, + 5.9948776299576965, + 5.599733263206939, + 5.425845828009287, + 7.505218418787916, + 5.233609737801514, + 5.427125269986897, + 5.379546413464169, + 5.005961183906842, + 5.893387001024662, + 4.951031661362545, + 5.567837253081451, + 4.997314126976024, + 5.983917896087184, + 5.5179743853345, + 4.410461998685994, + 5.111617104499174, + 5.293821498779172 + ] + }, + "homogeneous": { + "statistic": 6.450032327736446, + "reference": [ + 7.325305140324069, + 6.092936846234802, + 5.97787511650872, + 6.319762474386711, + 5.59939581465682, + 5.947490124653801, + 5.938469879304639, + 6.989714138313834, + 6.599462236925682, + 6.253629247435942, + 6.268952003306721, + 6.71238972278363, + 5.998773989074392, + 6.642238583748685, + 5.060370994545393, + 5.943647054228487, + 5.6774922105981735, + 6.087482382227443, + 6.457883341478321, + 6.5537880825688735, + 6.081276862424513, + 6.016756696360435, + 6.214165698736111, + 5.804321372017244, + 6.502287593492488, + 6.3767956537928265, + 5.907801616258618, + 6.431346743910252, + 6.053265550006875, + 6.416817655386247, + 6.524426051519425, + 6.282811751213115, + 5.46229311304185, + 6.665560005113521, + 6.445135203011679, + 6.073972320808133, + 8.065765837744179, + 6.248489395531836, + 6.616909732035294, + 6.491114065730435, + 5.916013841867676, + 7.006978715841566, + 6.2500968200449725, + 6.2878245525865655, + 5.924322617261755, + 7.092850667523054, + 6.828398469832565, + 6.314887883517232, + 5.8329156361119665, + 5.960573920968978 + ] + } +} diff --git a/benchmarks/reference/cao_dowd/loo_reference_stationary_r05_rep1.npy b/benchmarks/reference/cao_dowd/loo_reference_stationary_r05_rep1.npy new file mode 100644 index 0000000..0a32698 Binary files /dev/null and b/benchmarks/reference/cao_dowd/loo_reference_stationary_r05_rep1.npy differ diff --git a/benchmarks/reference/cao_dowd/results_41.json b/benchmarks/reference/cao_dowd/results_41.json new file mode 100644 index 0000000..3fc89db --- /dev/null +++ b/benchmarks/reference/cao_dowd/results_41.json @@ -0,0 +1,363 @@ +{ + "spec": 41, + "units": 41, + "k": 4, + "exposed": [ + "CA", + "AZ", + "NV", + "OR" + ], + "scm_pool_units": 39, + "years": [ + 1970, + 1971, + 1972, + 1973, + 1974, + 1975, + 1976, + 1977, + 1978, + 1979, + 1980, + 1981, + 1982, + 1983, + 1984, + 1985, + 1986, + 1987, + 1988, + 1989, + 1990, + 1991, + 1992, + 1993, + 1994, + 1995, + 1996, + 1997, + 1998, + 1999, + 2000 + ], + "T_pre": 19, + "obs": [ + 123.0, + 121.0, + 123.5, + 124.4, + 126.7, + 127.1, + 128.0, + 126.4, + 126.1, + 121.9, + 120.2, + 118.6, + 115.4, + 110.8, + 104.8, + 102.8, + 99.7, + 97.5, + 90.1, + 82.4, + 77.8, + 68.7, + 67.5, + 63.4, + 58.6, + 56.4, + 54.5, + 53.8, + 52.3, + 47.2, + 41.6 + ], + "scm": [ + 121.03559675240484, + 121.9435247894134, + 124.62564982136527, + 124.04070839569889, + 126.29778109217665, + 126.12802575303073, + 128.65051415323472, + 126.31397426102063, + 124.96673051012118, + 122.44839588146868, + 120.9371650160719, + 118.95092013779544, + 116.04728550944611, + 111.7811444065151, + 103.27498774397753, + 102.38763020072791, + 98.95958806608812, + 97.3878117968743, + 91.82256571256887, + 88.18328500618372, + 82.09672831252293, + 76.0301778241935, + 73.57338269601927, + 72.25277628640136, + 69.39971913133152, + 69.4133396221185, + 67.02714436932618, + 66.70761472316771, + 67.67314150905867, + 65.8674057593783, + 58.97908560234633 + ], + "sp": [ + 121.24025201985182, + 121.85164771954157, + 124.34375048250212, + 124.32559952884318, + 126.1149625330545, + 125.9488663972927, + 128.5750054120392, + 126.50253096052398, + 125.10683533268713, + 122.81241832324596, + 120.98739586607343, + 119.26192101749842, + 116.03209277628494, + 110.96481527315402, + 103.45003108534273, + 102.41104629479163, + 98.90710067337824, + 97.2006109945064, + 91.96311730938797, + 83.36690460908169, + 73.61966627261683, + 73.529456566541, + 72.01898304693351, + 71.95298465583045, + 67.53608153102081, + 70.36233068413517, + 68.66252141140399, + 69.73643399235016, + 68.27189945432146, + 67.443906563139, + 58.42401727770659 + ], + "te": [ + -0.9669046090816749, + 4.180333727383178, + -4.829456566541002, + -4.5189830469335135, + -8.552984655830452, + -8.936081531020804, + -13.96233068413518, + -14.162521411403993, + -15.936433992350157, + -15.971899454321463, + -20.243906563138985, + -16.824017277706588 + ], + "te_ci": [ + -4.301451484935192, + 3.4234562426899378 + ], + "te_pvalue": [ + 0.5789473684210527, + 0.05263157894736842, + 0.0, + 0.05263157894736842, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0 + ], + "spillover_test": [ + 1, + 1, + 1, + 1, + 1, + 1, + 0, + 0, + 0, + 0, + 0, + 0 + ], + "lambda_min": 0.5432517687088025, + "lambda_max": 2.3584104656262226, + "cond": 4.341284467847529, + "support_median": 6.0, + "support_max": 11, + "uAZ": [ + 3.3103681336869215, + -6.45325824639265, + 3.3754690968767562, + 4.752151936238064, + 6.201123949806444, + -4.196255414695088, + -5.020251445678099, + -3.659270385555118, + 0.026866500090392265, + 1.7948122473492987, + 4.8009740048636615, + -4.998047750576575, + -2.859002330333297, + -2.0229462425067055, + 2.150631271380945, + -0.45283409812044795, + 2.2082359378272542, + -0.4757340405273567, + 1.5169668762656556, + 8.265425316203078, + 4.78129556781991, + 0.13074069608386005, + -0.05658923251886558, + 0.347340626512099, + 0.938785768581426, + 4.628281856655049, + -1.6890816683903536, + -8.052261396753291, + -7.3654406794747445, + -8.022062225080722, + -2.0070304615289194 + ], + "ciAZ": [ + -7.096274347804453, + 6.19328592671738 + ], + "pAZ": [ + 0.0, + 0.3684210526315789, + 1.0, + 1.0, + 1.0, + 0.8421052631578947, + 0.3684210526315789, + 0.7894736842105263, + 0.0, + 0.0, + 0.0, + 0.7368421052631579 + ], + "uNV": [ + 8.346087256579665, + 2.938534673413301, + -0.5050726679093316, + 2.2714838398662494, + 5.434554301550122, + 7.259931115009579, + -0.47498006819421335, + -5.379828508176388, + -5.1709323951714055, + -1.8430041294749913, + -1.2963666921102117, + -2.027046043799544, + -3.6956423564202225, + 0.40888817048113424, + -14.092264759048703, + -2.0825094973670844, + -1.2560763615432506, + 7.081009407992468, + 4.083234714322913, + 22.4600579714719, + 38.20393610063767, + 17.88399353402303, + 13.21104535400626, + 8.40732800647911, + 10.593387490447382, + -1.2109781484234468, + -3.974768699056019, + -11.628995745580315, + -5.817554476847964, + 0.013252027763283457, + 4.319996095319082 + ], + "ciNV": [ + -13.865471892448491, + 18.001750858109517 + ], + "pNV": [ + 0.0, + 0.0, + 0.05263157894736842, + 0.15789473684210525, + 0.3684210526315789, + 0.2631578947368421, + 0.7894736842105263, + 0.631578947368421, + 0.15789473684210525, + 0.47368421052631576, + 1.0, + 0.5789473684210527 + ], + "uOR": [ + 4.812230565217504, + 2.6924433329947988, + -0.5550691833147141, + -1.0659292174167092, + -0.1762138560278288, + -1.3545985116879535, + -5.472820102695124, + -1.0006508474309328, + 3.4973808213142092, + 3.116760578636274, + 2.8602830301705637, + 2.3339963324607993, + -0.887560780063211, + -6.849005564256743, + -0.1408795123076345, + 2.0495466199495525, + -0.8302669611155977, + -4.877888997242877, + 1.8482422528155524, + 18.362811645086147, + 32.06086366680145, + 31.2613822548075, + 32.379474175997075, + 27.77217698600626, + 24.661771379314796, + 17.57889797879717, + 19.298490892215817, + 7.672237178342751, + 2.4385243959894547, + -1.5638299984691824, + 7.368234118625371 + ], + "ciOR": [ + -11.378110249050598, + 12.068001901431465 + ], + "pOR": [ + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.21052631578947367, + 0.5789473684210527, + 0.5789473684210527, + 0.21052631578947367 + ], + "ca_weights": { + "CT": 0.24098156758036102, + "NV": 0.16377009026101794, + "NE": 0.15511874473566314, + "OR": 0.12074738829099117, + "CO": 0.11994544277848386, + "IL": 0.07894471165396805, + "NC": 0.035914677722353586, + "NH": 0.03293865571177756, + "KS": 0.031129946785812313, + "MT": 0.02050877421432828 + }, + "ca_intercept": -20.032016652525527, + "ca_rmspe_pre": 0.9180912378829023 +} \ No newline at end of file diff --git a/benchmarks/reference/cao_dowd/results_51.json b/benchmarks/reference/cao_dowd/results_51.json new file mode 100644 index 0000000..9a0c298 --- /dev/null +++ b/benchmarks/reference/cao_dowd/results_51.json @@ -0,0 +1,374 @@ +{ + "spec": 51, + "units": 51, + "k": 14, + "exposed": [ + "CA", + "AK", + "AZ", + "DC", + "FL", + "HI", + "MA", + "MD", + "MI", + "NJ", + "NV", + "NY", + "OR", + "WA" + ], + "scm_pool_units": 39, + "years": [ + 1970, + 1971, + 1972, + 1973, + 1974, + 1975, + 1976, + 1977, + 1978, + 1979, + 1980, + 1981, + 1982, + 1983, + 1984, + 1985, + 1986, + 1987, + 1988, + 1989, + 1990, + 1991, + 1992, + 1993, + 1994, + 1995, + 1996, + 1997, + 1998, + 1999, + 2000 + ], + "T_pre": 19, + "obs": [ + 123.0, + 121.0, + 123.5, + 124.4, + 126.7, + 127.1, + 128.0, + 126.4, + 126.1, + 121.9, + 120.2, + 118.6, + 115.4, + 110.8, + 104.8, + 102.8, + 99.7, + 97.5, + 90.1, + 82.4, + 77.8, + 68.7, + 67.5, + 63.4, + 58.6, + 56.4, + 54.5, + 53.8, + 52.3, + 47.2, + 41.6 + ], + "scm": [ + 121.03559675240484, + 121.9435247894134, + 124.62564982136527, + 124.04070839569889, + 126.29778109217665, + 126.12802575303073, + 128.65051415323472, + 126.31397426102063, + 124.96673051012118, + 122.44839588146868, + 120.9371650160719, + 118.95092013779544, + 116.04728550944611, + 111.7811444065151, + 103.27498774397753, + 102.38763020072791, + 98.95958806608812, + 97.3878117968743, + 91.82256571256887, + 88.18328500618372, + 82.09672831252293, + 76.0301778241935, + 73.57338269601927, + 72.25277628640136, + 69.39971913133152, + 69.4133396221185, + 67.02714436932618, + 66.70761472316771, + 67.67314150905867, + 65.8674057593783, + 58.97908560234633 + ], + "sp": [ + 121.92477664581043, + 121.43813025084717, + 124.58294958790344, + 123.86259716809457, + 126.40515656113212, + 126.87895357348475, + 128.1980541000673, + 126.58778195224929, + 125.37729644994093, + 121.865364636675, + 120.56657160036374, + 118.8518894769351, + 115.80284471431872, + 110.89873511012036, + 103.73309752689457, + 103.51188423529727, + 99.84228459368529, + 96.82189439469875, + 90.84973742148075, + 82.31539627769534, + 74.08353058672053, + 72.45593330017549, + 70.9246531593901, + 71.01264847055259, + 69.51252951174925, + 69.23257831473674, + 67.58224452878808, + 68.71214038203074, + 68.37973974258463, + 66.15802036663672, + 57.08840913445824 + ], + "te": [ + 0.08460372230465855, + 3.716469413279463, + -3.755933300175487, + -3.4246531593901097, + -7.61264847055259, + -10.912529511749254, + -12.832578314736738, + -13.082244528788076, + -14.912140382030742, + -16.07973974258463, + -18.95802036663671, + -15.488409134458237 + ], + "te_ci": [ + -3.9593189900494346, + 3.1873882165643317 + ], + "te_pvalue": [ + 1.0, + 0.05263157894736842, + 0.05263157894736842, + 0.05263157894736842, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0 + ], + "spillover_test": [ + 0, + 1, + 1, + 1, + 0, + 1, + 1, + 1, + 1, + 1, + 1, + 1 + ], + "lambda_min": 0.221540129924041, + "lambda_max": 2.7658216904010406, + "cond": 12.484517777205204, + "support_median": 7.0, + "support_max": 14, + "uAZ": [ + 2.292997202226587, + -5.198484845445833, + 3.9976029636314223, + 0.8341107181678211, + 3.7519889306348375, + -3.8650601614093887, + -2.102713206885653, + -1.7527470889570793, + 2.7847254162654167, + 2.269796253899287, + 2.8731008704518217, + -3.9967741350591552, + -3.220927624281856, + -0.6132337181997656, + -0.389195846825217, + 0.2765672056979298, + 1.2433705851905046, + -0.4193289878144384, + 1.2342054687125454, + 4.990980276234831, + -11.241447345143545, + -15.365465721044684, + -16.54884994458616, + -15.136040055574123, + -16.051291848670065, + -2.958225438679426, + -7.195445817822032, + -10.47803389531696, + -9.691153778055874, + -10.433079588322817, + -7.474339703561894 + ], + "ciAZ": [ + -4.986875678947999, + 6.720650965301914 + ], + "pAZ": [ + 0.15789473684210525, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.42105263157894735, + 0.05263157894736842, + 0.0, + 0.0, + 0.0, + 0.05263157894736842 + ], + "uNV": [ + 5.632028894647021, + -0.6453364058157831, + -4.217962477787239, + -1.8724447311115568, + 3.5837051645347913, + 5.953151968555034, + 0.44528016036582585, + -1.6595125472648533, + -1.0221319703341578, + 1.4132613373523952, + 1.3516843609417322, + 0.1970844535527192, + -2.3333942815142734, + 1.6259094940288321, + -14.180041948890334, + -2.5905159633315122, + -1.9606618744867887, + 6.643019081300224, + 3.636877285257583, + 14.967365232650836, + 26.867691299812964, + 3.831292134152011, + -1.6090392380305385, + -5.11965642898089, + 2.670576730965645, + -9.684297383683875, + -12.39608608581913, + -13.869418503679515, + -8.657621299393275, + -1.4649268144588146, + -1.8933719030512286 + ], + "ciNV": [ + -12.442938604918666, + 12.685665747299975 + ], + "pNV": [ + 0.0, + 0.0, + 0.5263157894736842, + 0.7894736842105263, + 0.47368421052631576, + 0.5789473684210527, + 0.21052631578947367, + 0.15789473684210525, + 0.0, + 0.21052631578947367, + 0.8421052631578947, + 0.7368421052631579 + ], + "uOR": [ + 4.7468854507849585, + 2.4784980163187242, + -0.9839245054496928, + -1.5741195984281546, + -0.016740424443554502, + -1.8488602257302773, + -5.136193219851776, + 0.006410872304485338, + 4.387384210938649, + 3.112182121412502, + 3.0466562347176875, + 1.8804887465886395, + -1.4648783668055785, + -6.651272353711562, + -0.4951175607804714, + 1.8427428596517927, + -1.1568684088146028, + -4.129676005556291, + 1.9564021568544234, + 13.90300976863315, + 26.222668606786367, + 23.45568383956229, + 23.332388795383388, + 19.7609324074306, + 19.42945173942901, + 11.959813554933358, + 14.469738836173747, + 6.0048931174590425, + 0.9922225086838468, + -2.5216571468559836, + 4.710534443942803 + ], + "ciOR": [ + -8.565954982928405, + 8.960787151832687 + ], + "pOR": [ + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.21052631578947367, + 0.7894736842105263, + 0.5263157894736842, + 0.2631578947368421 + ], + "ca_weights": { + "OR": 0.27550305237933004, + "MA": 0.20628963656673088, + "AZ": 0.14803408405414048, + "AK": 0.10078165512882611, + "NV": 0.06899440483918341, + "CT": 0.06131172516757177, + "MN": 0.03565451939125902, + "HI": 0.0345533459713232, + "KS": 0.03322940251049153, + "NH": 0.0305515791416624, + "DC": 0.005096594615076314 + }, + "ca_intercept": -16.163852521865493, + "ca_rmspe_pre": 0.589402935565328 +} \ No newline at end of file diff --git a/docs/_quarto.yml b/docs/_quarto.yml index af86949..fe48b18 100644 --- a/docs/_quarto.yml +++ b/docs/_quarto.yml @@ -25,6 +25,8 @@ website: - text: "---" - text: "sar — Bayesian spatial-autoregressive spillover SCM" file: models/sar.qmd + - text: "cd — Cao–Dowd SCM with spillovers" + file: models/cd.qmd - text: "Sudan case study" file: sudan.qmd - text: "For AI / LLMs" @@ -83,18 +85,16 @@ quartodoc: - title: Estimator desc: >- Configure once, fit once. `SCSPILL(config).fit()` runs the spillover - model named by `method` — today `"sar"`, whose two-step sampler fits - the horseshoe synthetic weights and then the SAR spillover block — and - returns a `SCSPILLResults` with the treatment effect, its credible - interval, and the spillover received by every donor. The + model named by `method`: `"sar"` runs the Bayesian spatial sampler; + `"cd"` runs Cao–Dowd all-unit SCM and P-test inference. It returns + the model-specific result and a shared effect surface. The [model catalogue](../models/index.qmd) says which models exist. package: scspill.estimators.scspill contents: - SCSPILL - title: Configuration desc: >- - The pydantic input contract — `method` (the spillover model; `"sar"` - is the only accepted value today), the long panel and its column + The pydantic input contract for `"sar"` — the long panel and its column names, the spatial weights `spatial_w` / `spatial_W` (label-aligned), covariates, the MCMC budget, priors, and the paper-correct switches (`beta_prior`, `propagate_alpha`, `adapt_rho`) with their R-compatible @@ -102,6 +102,14 @@ quartodoc: package: scspill.utils.scspill_helpers.config contents: - SCSPILLConfig + - title: Cao–Dowd configuration and results + desc: >- + Fixed exposure sets or labeled A matrices, all-unit SCM equations, + pointwise P tests, and kappa_A structure tests. + package: scspill.utils.scspill_helpers.cd + contents: + - CDConfig + - CDResults - title: Results desc: >- The standardized posterior report — flat accessors (`att`, `att_ci`, @@ -164,14 +172,14 @@ quartodoc: - mcmc_summary - title: Datasets desc: >- - Two bundled spillover case studies, shipped inside the wheel — the - California Proposition 99 tobacco panel (39 states, 1970–2000, rook - contiguity weights) and the 2011 Sudan secession panel (34 African - countries, 2000–2015, bilateral-trade weights). Each loader returns a - `SpillPanel` whose `config_kwargs()` feeds any spillover model's - configuration directly. + Three bundled panels: the 39-state California and Sudan SAR case + studies, plus the 51/41-unit California data used by Cao–Dowd with + treatment beginning in 1989. Each loader supplies model-specific + configuration fields. package: scspill.data contents: - load_california - load_sudan - SpillPanel + - load_california_cd + - CDPanel diff --git a/docs/articles/datasets.qmd b/docs/articles/datasets.qmd index 4f59882..b452e08 100644 --- a/docs/articles/datasets.qmd +++ b/docs/articles/datasets.qmd @@ -1,9 +1,9 @@ --- title: "The bundled datasets" -subtitle: "Provenance and dictionaries for the two case studies" +subtitle: "Provenance and dictionaries for the SAR and Cao–Dowd case studies" --- -Both of the paper's empirical applications ship inside the wheel as plain +Both of the SAR paper's empirical applications ship inside the wheel as plain CSVs (about 150 kB total), copied verbatim from the replication package's nonproprietary export. Each loader returns a `SpillPanel` — the long panel with a 0/1 `treated` column, the raw spatial weights aligned by donor label, @@ -11,6 +11,14 @@ and a `config_kwargs()` helper that feeds `SCSPILLConfig` directly. The weights are stored *unnormalized*; the estimator row-normalizes `spatial_W` and scales `spatial_w` to sum to one, exactly as the R package does. +## Cao–Dowd Proposition 99 + +`load_california_cd(51)` and `load_california_cd(41)` load a separate +51-region panel (50 states and DC, 1970–2000, treatment from 1989) taken +from the authors' [replication repository](https://github.com/jcao0/synthetic-control-spillover). +It is not interchangeable with `load_california()`; see the +[`cd` model page](../models/cd.qmd). + ## California Proposition 99 ```{python} diff --git a/docs/index.qmd b/docs/index.qmd index 22571ad..a88c9ac 100644 --- a/docs/index.qmd +++ b/docs/index.qmd @@ -67,12 +67,11 @@ RMSE, and coverage across estimators. | Model | `method` | Reference | Status | |---|---|---|---| | [Bayesian spatial-autoregressive spillover SCM](models/sar.qmd) | `"sar"` | Sakaguchi & Tagawa (2026), *Econometrics Journal* | **Available** | -| SCM with spillover effects | `"cd"` | Cao & Dowd (2019, rev. 2026) | Planned | +| [SCM with spillover effects](models/cd.qmd) | `"cd"` | Cao & Dowd (2019, rev. 2026) | **Available** | | Inclusive synthetic control | `"iscm"` | Di Stefano & Mellace (2024) | Planned | | Partial-interference SCG | `"grossi"` | Grossi et al. (2025), *JRSS-A* | Planned | -The planned rows are a roadmap, not a feature list: they are not implemented -and their `method` names are rejected rather than silently substituted. The +The remaining planned rows are a roadmap. Their `method` names are rejected. The [model catalogue](models/index.qmd) explains what each would assume, and why having several matters — the assumptions are not nested, so two models bracket the answer. diff --git a/docs/llms.txt b/docs/llms.txt index 443d1fa..c73ac43 100644 --- a/docs/llms.txt +++ b/docs/llms.txt @@ -8,13 +8,14 @@ Install: `pip install scspill` (or `pip install "scspill[numba]"` for JIT-compiled samplers). -Models: exactly ONE is implemented -- `method="sar"` (Sakaguchi & Tagawa, -The Econometrics Journal 2026), class SCSPILL, and the default. It routes +Models: TWO are implemented -- `method="sar"` (Sakaguchi & Tagawa, +The Econometrics Journal 2026), the default, and `method="cd"` (Cao & Dowd, +arXiv:1902.07343v4). `sar` routes spillovers through user-supplied spatial weights scaled by one estimated intensity rho, via two-step MCMC (horseshoe synthetic weights, then a SAR -spillover block). The catalogue also lists `cd` (Cao & Dowd), `iscm` (Di -Stefano & Mellace) and `grossi` (Grossi et al.) as PLANNED: they are NOT -implemented and SCSPILLConfig rejects them. Do not emit code using them. +spillover block). `cd` fits intercept-adjusted all-unit SCM with a fixed +exposure matrix and P-test inference. `iscm` (Di Stefano & Mellace) and +`grossi` (Grossi et al.) remain planned and are not implemented. Catalogue: https://quarcs-lab.github.io/scspill/models/index.html Minimal contract (model `sar`): SCSPILL({df, outcome, treat, unitid, time, @@ -24,11 +25,19 @@ result with .att, .att_ci, .rho_hat, .rho_ci, .counterfactual, treated unit and date are inferred from the 0/1 `treat` column; spatial weights align by donor label. +Minimal contract (model `cd`): SCSPILL({method: "cd", df, outcome, treat, +unitid, time, affected_units? or effect_matrix?, reference_mode?, +quantile_method?, ci_method?}).fit() -> CDResults with .effects_panel, +.effect_ci (= .published_ci by default, or .pointwise_ci), +.pointwise_p_values, .kappa_test, and .test_hypothesis(C). The 51-unit +1989-onward replication panel is load_california_cd(51). + ## Docs - [Get started (California Prop 99)](https://quarcs-lab.github.io/scspill/get-started.html) - [Models: catalogue and roadmap](https://quarcs-lab.github.io/scspill/models/index.html) - [The `sar` model (Sakaguchi & Tagawa)](https://quarcs-lab.github.io/scspill/models/sar.html) +- [The `cd` model (Cao & Dowd)](https://quarcs-lab.github.io/scspill/models/cd.html) - [Sudan case study](https://quarcs-lab.github.io/scspill/sudan.html) - [Validating the sampler](https://quarcs-lab.github.io/scspill/articles/validation.html) - [The Monte Carlo simulation study](https://quarcs-lab.github.io/scspill/articles/simulation-study.html) @@ -42,7 +51,7 @@ weights align by donor label. ### Top level -`SCSPILL`, `SCSPILLConfig`, `SCSPILLResults` +`SCSPILL`, `SCSPILLConfig`, `SCSPILLResults`, `CDConfig`, `CDResults` ### scspill.validation @@ -54,7 +63,7 @@ weights align by donor label. ### scspill.data -`SpillPanel`, `load_california`, `load_sudan` +`CDPanel`, `SpillPanel`, `load_california`, `load_california_cd`, `load_sudan` ## Full corpus diff --git a/docs/models/cd.qmd b/docs/models/cd.qmd new file mode 100644 index 0000000..3bd6d76 --- /dev/null +++ b/docs/models/cd.qmd @@ -0,0 +1,152 @@ +--- +title: "`cd` — Cao–Dowd synthetic control with spillovers" +subtitle: "Single treated unit, fixed exposure structure, and end-of-sample inference" +--- + +`method="cd"` implements [Cao and Dowd, v4](https://arxiv.org/html/1902.07343v4). +It fits an intercept and simplex synthetic-control weights for **every** unit +using pre-treatment outcomes. The researcher supplies a set of potentially +exposed units or a labeled, fixed structure matrix $A$. The model permits a +different direct and spillover effect in each post-treatment period. + +## Fit + +```python +from scspill import SCSPILL +from scspill.data import load_california_cd + +panel = load_california_cd(51) # 50 states + DC; 1970–2000; post starts 1989 +result = SCSPILL({**panel.config_kwargs(), "display_graphs": False}).fit() +result.effects_panel.loc[1990, ["CA", "NV"]] +result.effect_ci["CA"].loc[1990] +``` + +The existing `load_california()` is a different 39-state, 1988-onward case +study and does not reproduce this paper. `load_california_cd(41)` supplies +the paper's reduced specification. + +For a custom panel, pass `df`, `outcome`, `treat`, `unitid`, `time`, and either +`affected_units` or a pandas `effect_matrix`. The treatment flag is 1 only +for the single treated unit from the first post period onward. The +`effect_matrix` rows must cover every unit by label and include one column +equal to the treated-unit indicator. Other columns may encode shared or +heterogeneous spillovers. All columns of $(I-\widehat B)A$ must be linearly +independent; the fit raises an identification error otherwise. + +Writing $R=I-\widehat B$, the period-specific effect is + +$$ +\widehat\alpha_t=A[(RA)'(RA)]^{-1}(RA)'(RY_t-\widehat a). +$$ + +## Inference + +The defaults follow the authors' Proposition 99 replication code, so a +default fit reproduces the paper's Section 6 tables and Figure 4. The +departures that convention makes from the paper's own Section 4 definitions +are listed under [Caveats](#caveats). + +`result.test_hypothesis(C, d=0, period=1990)` tests a general linear +restriction $C\alpha_t=d$ by the squared P statistic. `C` can be a labeled +DataFrame or a matrix in `result.unit_labels` order. Identity and inverse +projected second-moment weights are available. Set +`reference_mode="leave_one_period_out"` to refit all SCM equations for each +held-out pre-treatment period; `"full"` matches the application code and is +the default. + +`result.effect_ci[unit]` is the pointwise interval selected by `ci_method`. +Both versions are always stored: `published_ci` (the default, as in the +paper's tables) and `pointwise_ci` (the inverse of the squared P test). +`pointwise_p_values` holds the period-wise P values, and +`kappa_test[period]` checks the assumed structure $A$ period by period +(Section 5.1.2). `att` is the arithmetic mean of the post-period direct +estimates; no ATT confidence interval is implied by the pointwise intervals. + +## Caveats {#caveats} + +The paper's empirical results and its stated procedure differ in three +places. `cd` follows the empirical results by default; each difference has a +switch. + +| # | Paper, Section 4 | Authors' application code (default) | Switch | +|---|---|---|---| +| 1 | Critical value $\inf\{x:\widehat F(x)\ge 1-\tau\}$, the left inverse of the empirical CDF | NumPy's linearly interpolated quantile | `quantile_method="empirical"` | +| 2 | Intervals by inverting the squared P test, which makes them symmetric about $\widehat\alpha$ | $\widehat\alpha + [q_{\tau/2},\, q_{1-\tau/2}]$ of the signed pre-treatment errors $e_i'\widehat G\widehat u_t$ | `ci_method="inversion"` | +| 3 | Leave-one-period-out references recommended for finite samples | Full-sample references | `reference_mode="leave_one_period_out"` | + +What each one changes in practice: + +1. With only 19 pre-treatment periods the two rules pick different order + statistics. The paper's 95% critical value is + $\lceil 0.95 \times 19 \rceil = 19$th, the largest reference statistic, + so it rejects only when the post-treatment statistic exceeds all 19. The + interpolated quantile sits at position 18.1, just above the second + largest. In the 51-unit application the joint test for any spillover + rejects in 1990–1992 and 1994–2000 under the default, as in Figure 4. In + 1994 and 1995 the statistic (1745.8 and 1725.3) lies between the + interpolated cut-off (1482.8) and the largest reference value (1828.9), + so the paper's definition does not reject, and the paper's own P value, + $\#\{\widehat P_t \ge P\}/19 = 0.053$, is above 0.05. +2. Since $\widehat\alpha-\alpha\approx Gu_{T+1}$, a percentile interval for + $\alpha$ subtracts the quantiles, $[\widehat\alpha-q_{1-\tau/2},\, + \widehat\alpha-q_{\tau/2}]$. The published interval adds them, which + agrees with that only when the error distribution is symmetric; for + California in 1989 it gives $[-3.87, 3.27]$ where the subtracting form + gives $[-3.11, 4.04]$. In the authors' repository the Proposition 99 + code (`run_spec.py` computes the quantiles, `make_figures.py` plots + `te + lo` to `te + hi`) adds, while their continuous-treatment + application (`continuous_treatment.m`) subtracts. +3. Leave-one-out references widen the intervals substantially in the + California application; for California itself the range of the + reference errors more than doubles. + +The authors' own studies are not uniform either: their Monte Carlo size and +power studies and the $\kappa_A$ study use the left-inverse critical value +(`quantile_method="empirical"`), and the repository's README notes that its +confidence intervals were checked against the original MATLAB output only +for the point estimates. + +## Reproducing the paper + +The bundled panel is the input of the [authors' replication +package](https://github.com/jcao0/synthetic-control-spillover) at commit +`86004020`. A default fit reproduces Section 6: + +```python +import numpy as np +import pandas as pd +from scspill import SCSPILL +from scspill.data import load_california_cd + +fit = SCSPILL({**load_california_cd(51).config_kwargs(), "display_graphs": False}).fit() + +# Table S.6: estimate, 95% interval and P value for CA (direct) and NV, AZ, OR +table_s6 = pd.concat( + { + unit: fit.effect_ci[unit].assign( + estimate=fit.effects_panel[unit], p_value=fit.pointwise_p_values[unit] + ) + for unit in ["CA", "NV", "AZ", "OR"] + }, + axis=1, +).round(2) + +# Figure 4b shading: years in which "no spillover anywhere" is rejected +no_spill = np.eye(len(fit.unit_labels))[1:] +[t for t in fit.effects_panel.index if fit.test_hypothesis(no_spill, period=t).reject] + +# Table S.8: California's synthetic-control weights and intercept +fit.scm_matrix.loc["CA"].nlargest(10), fit.intercepts["CA"] +``` + +`load_california_cd(41)` gives the robustness specification of Table S.7. +`scspill/tests/test_cd.py` checks both fits against the authors' saved +output: P values match exactly, and effects agree to within 0.01 packs. +The remaining gap is the simplex solver. With the authors' settings, SciPy's +SLSQP stops short of the optimum for every unit of the 51-unit panel; `cd` +uses scaled objectives with analytic gradients and reaches an objective at +or below theirs for every unit, so a few published two-decimal values +differ in the last digit. + +Multiple treated units, covariates, time-varying $A$, and alternative +estimators are outside this implementation. diff --git a/docs/models/index.qmd b/docs/models/index.qmd index ddf5d1b..a30202a 100644 --- a/docs/models/index.qmd +++ b/docs/models/index.qmd @@ -27,17 +27,13 @@ worth having more than one — run two and they bracket the answer. ## Available now -::: {.callout-important} -## One model ships today -scspill implements **one** model. Everything under -[Planned](#planned) below is a roadmap: none of it is implemented, and none of -those `method` names is accepted. Passing one raises a validation error rather -than silently falling back. -::: +scspill implements two models. `sar` remains the default, and `cd` uses a +specified spillover structure and pre-treatment P-test reference distribution. | Model | `method` | Reference | Class | Page | |---|---|---|---|---| | Bayesian spatial-autoregressive spillover SCM | `"sar"` | Sakaguchi & Tagawa (2026), *The Econometrics Journal*, [doi:10.1093/ectj/utag006](https://doi.org/10.1093/ectj/utag006) | `SCSPILL` | [`sar` →](sar.qmd) | +| SCM with spillover effects | `"cd"` | Cao & Dowd, [arXiv:1902.07343v4](https://arxiv.org/html/1902.07343v4) | `SCSPILL` | [`cd` →](cd.qmd) | ### `sar` — spillovers as a spatial process @@ -68,18 +64,27 @@ Full statement of the model, the identification theorem, the samplers, and the six documented departures from the authors' R replication package are on [the `sar` page](sar.qmd). +### `cd` — spillovers under a specified structure + +`cd` estimates an intercept and simplex synthetic-control weights for each +unit, then uses a fixed exposure matrix $A$ to separate direct and spillover +effects by post-treatment period. Andrews P tests use pre-treatment variation; +the $\kappa_A$ test checks whether the specified structure fits. The +[model page](cd.qmd) documents the 51-state Proposition 99 replication, +inference options, and the distinction between the paper's P-test interval +and the published application's interval calculation. + ## Planned {#planned} ::: {.callout-warning} ## Not implemented -The three models below are candidates, not features. They are listed so the +The two models below are candidates, not features. They are listed so the package's intended scope is legible and so the `method` names stay coherent — **none of them exists in scspill today**, and there is no release date. ::: | Model | Reserved `method` | Reference | How it *would* identify the spillover | |---|---|---|---| -| SCM with spillover effects | `"cd"` | Cao & Dowd, [arXiv:1902.07343](https://arxiv.org/abs/1902.07343) (2019, rev. 2026) | *Keep* the exposed donors and impose a researcher-specified spillover-structure matrix, linear in unknown parameters, recovering direct and per-unit spillover effects jointly. Comes with a misspecification test for the structure you assumed. | | Inclusive synthetic control | `"iscm"` | Di Stefano & Mellace, [arXiv:2403.17624](https://arxiv.org/abs/2403.17624) (2024) | *Keep* the exposed donors with no parametric spillover structure at all: build a synthetic control for every affected unit, then invert the system of their own cross-weights to de-contaminate the gaps. | | Partial-interference synthetic control group | `"grossi"` | Grossi, Mariani, Mattei, Lattarulo & Öner (2025), *JRSS-A* 188(1):223–240, [arXiv:2004.05027](https://arxiv.org/abs/2004.05027) | *Drop* the treated unit's whole cluster from the pool and rebuild from the far, clean clusters only. Cleanest identification when partial interference is defensible, at the cost of pre-treatment fit. | diff --git a/docs/models/sar.qmd b/docs/models/sar.qmd index 5c82bd3..8a095b1 100644 --- a/docs/models/sar.qmd +++ b/docs/models/sar.qmd @@ -5,8 +5,8 @@ aliases: - /articles/method.html --- -`sar` is scspill's first model, and today its only one — the -[model catalogue](index.qmd) lists what else is planned. It implements +`sar` is scspill's default model. The [model catalogue](index.qmd) also lists +the implemented [`cd` model](cd.qmd) and future work. `sar` implements Sakaguchi & Tagawa's *Identification and Bayesian Inference for Synthetic Control Methods with Spillover Effects*. This page states the model, the identification result, and the samplers, and documents precisely where this diff --git a/docs/use-with-llms.qmd b/docs/use-with-llms.qmd index 9dd7364..055dd81 100644 --- a/docs/use-with-llms.qmd +++ b/docs/use-with-llms.qmd @@ -51,6 +51,11 @@ result.diagnostics() # DataFrame: mean/sd/quantiles/ESS/R-hat per chain result.plot(kind="panel") # counterfactual | effect | top spillovers ``` +For Cao–Dowd, pass `load_california_cd(51).config_kwargs()` to `SCSPILL`; +the `CDResults` it returns has `effects_panel`, `effect_ci`, +`pointwise_p_values`, `test_hypothesis` and `kappa_test`. See the +[`cd` model page](models/cd.qmd). + ## Rules of thumb for agents - The treated unit and the treatment date are inferred from the `treat` @@ -65,7 +70,7 @@ result.plot(kind="panel") # counterfactual | effect | top spillovers `ScspillPlottingError` — all subclass `ScspillError`. - Sampler validation lives in `scspill.validation` (`geweke_test`, `prior_sensitivity`, `prior_predictive`); Monte Carlo tooling in - `scspill.simulate`; the two case studies in `scspill.data`. + `scspill.simulate`; bundled panels in `scspill.data`. - The estimator's architecture mirrors [mlsynth](https://github.com/jgreathouse9/mlsynth) (pydantic config in, standardized results object out), so mlsynth patterns transfer. diff --git a/pyproject.toml b/pyproject.toml index c64ee26..ebc0720 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -87,7 +87,7 @@ exclude = ["scspill.tests*"] [tool.setuptools.package-data] scspill = ["py.typed"] -"scspill.data" = ["files/*.csv"] +"scspill.data" = ["files/*.csv", "files/cao_dowd/*.csv"] # --------------------------------------------------------------------------- # Tooling diff --git a/scspill/__init__.py b/scspill/__init__.py index 01a93e1..46df9c5 100644 --- a/scspill/__init__.py +++ b/scspill/__init__.py @@ -6,7 +6,7 @@ by every control unit. The package is organized around a model layer, selected by ``method``. -Exactly one model ships today: +Two models are implemented: ``"sar"`` The Bayesian spatial-autoregressive model of Sakaguchi & Tagawa (2026), @@ -19,15 +19,20 @@ block (latent AR(1) factors, covariates, and a random-walk Metropolis step for ``rho``) conditional on the posterior mean of ``alpha``. +``"cd"`` + Cao and Dowd's intercept-adjusted, all-unit SCM with a researcher-supplied + spillover structure. It reports period-specific direct and spillover + effects, Andrews P tests, confidence intervals, and specification tests. + Further spillover-aware models are planned behind the same interface; the catalogue at https://quarcs-lab.github.io/scspill/models/ tracks them and states plainly which are implemented. Public API:: - from scspill import SCSPILL, SCSPILLConfig + from scspill import SCSPILL, SCSPILLConfig, CDConfig - result = SCSPILL(config).fit() # -> SCSPILLResults (method="sar") + result = SCSPILL(config).fit() # -> SCSPILLResults or CDResults plus the companion subpackages :mod:`scspill.validation` (the ``sar`` model's Geweke joint distribution test, prior sensitivity, and prior predictive @@ -53,6 +58,8 @@ __all__ = [ "SCSPILL", + "CDConfig", + "CDResults", "SCSPILLConfig", "SCSPILLResults", "ScspillConfigError", @@ -69,6 +76,8 @@ "SCSPILL": ("scspill.estimators.scspill", "SCSPILL"), "SCSPILLConfig": ("scspill.utils.scspill_helpers.config", "SCSPILLConfig"), "SCSPILLResults": ("scspill.utils.scspill_helpers.structures", "SCSPILLResults"), + "CDConfig": ("scspill.utils.scspill_helpers.cd", "CDConfig"), + "CDResults": ("scspill.utils.scspill_helpers.cd", "CDResults"), } diff --git a/scspill/config_models.py b/scspill/config_models.py index e24084f..9e9c42f 100644 --- a/scspill/config_models.py +++ b/scspill/config_models.py @@ -582,6 +582,7 @@ def plot(self, kind: str = "auto", *, ax: Any = None, **overrides: Any) -> Any: # --------------------------------------------------------------------------- _RELOCATED_CONFIGS = { "SCSPILLConfig": "scspill.utils.scspill_helpers.config", + "CDConfig": "scspill.utils.scspill_helpers.cd", } diff --git a/scspill/data/__init__.py b/scspill/data/__init__.py index f72ef9a..003ab13 100644 --- a/scspill/data/__init__.py +++ b/scspill/data/__init__.py @@ -1,18 +1,20 @@ """Bundled case-study spillover panels for scspill. Ships the two empirical applications of the ``sar`` model's source article -(Sakaguchi & Tagawa 2026) inside the wheel, each as a long panel plus its -spatial-weight specification. Nothing about them is model-specific: they are -ordinary spillover panels and will serve any model added later. +(Sakaguchi & Tagawa 2026), plus Cao and Dowd's 51/41-unit California panel. +The latter has a different sample and intervention date from the SAR +California example. * :func:`load_california` -- the California Proposition 99 tobacco panel (Abadie, Diamond & Hainmueller 2010): 39 states, 1970-2000, treatment in 1988; rook-contiguity spatial weights. * :func:`load_sudan` -- the 2011 Sudan secession panel: 34 African countries, 2000-2015, treatment in 2011; average bilateral-trade spatial weights. +* :func:`load_california_cd` -- the Cao--Dowd Proposition 99 panel: 51 or + 41 regions, 1970-2000, treatment from 1989; labeled affected-unit set. -The CSV files are verbatim copies of the replication package's nonproprietary -data exports. Both weight objects are stored raw (unnormalized) -- +The SAR CSV files are verbatim copies of that replication package's +nonproprietary data exports. Both SAR weight objects are stored raw -- row-normalization of ``W`` and sum-normalization of ``w`` happen inside the estimator, exactly as in the R package. """ @@ -27,7 +29,53 @@ from ..exceptions import ScspillDataError -__all__ = ["SpillPanel", "load_california", "load_sudan"] +__all__ = ["CDPanel", "SpillPanel", "load_california", "load_california_cd", "load_sudan"] + + +@dataclass(frozen=True) +class CDPanel: + """Cao--Dowd Proposition 99 panel and published exposure structure.""" + + df: pd.DataFrame + affected_units: tuple[str, ...] + specification: int + + def config_kwargs(self) -> dict[str, Any]: + """Return inputs for ``SCSPILL({"method": "cd", ...})``.""" + return { + "method": "cd", + "df": self.df, + "outcome": "cigs", + "treat": "treated", + "unitid": "state", + "time": "year", + "affected_units": list(self.affected_units), + } + + +def load_california_cd(specification: int = 51) -> CDPanel: + """Load the authors' 1970--2000 Proposition 99 panel, treated from 1989. + + ``51`` uses every state and DC with the 12 policy-contaminated states and + three neighbors as exposed units. ``41`` removes ten non-neighbor policy + states, leaving Arizona, Nevada and Oregon as exposed controls. + """ + if specification not in (41, 51): + raise ScspillDataError("Cao--Dowd California specification must be 41 or 51.") + root = files("scspill.data").joinpath("files", "cao_dowd") + with root.joinpath("cigs_consumption.csv").open("r", encoding="utf-8") as fh: + df = pd.read_csv(fh) + with root.joinpath("state_list.csv").open("r", encoding="utf-8-sig") as fh: + flags = pd.read_csv(fh) + policy = set(flags.loc[flags["missing12"] == 1, "state"]) + neighbors = set(flags.loc[flags["neighbor"] == 1, "state"]) + if specification == 41: + df = df.loc[~df["state"].isin(policy - neighbors)].copy() + affected = tuple(sorted(neighbors)) + else: + affected = tuple(sorted(policy | neighbors)) + df["treated"] = ((df["state"] == "CA") & (df["year"] >= 1989)).astype(int) + return CDPanel(df=df, affected_units=affected, specification=specification) #: Mapping from the R-mangled WDI column headers in ``sudan_panel.csv`` to the diff --git a/scspill/data/files/cao_dowd/cigs_consumption.csv b/scspill/data/files/cao_dowd/cigs_consumption.csv new file mode 100644 index 0000000..cce1e59 --- /dev/null +++ b/scspill/data/files/cao_dowd/cigs_consumption.csv @@ -0,0 +1,1582 @@ +state,year,cigs +AK,1970,121.3 +AK,1971,123 +AK,1972,130 +AK,1973,125.8 +AK,1974,130.4 +AK,1975,150.4 +AK,1976,164.8 +AK,1977,159.3 +AK,1978,143.6 +AK,1979,138.1 +AK,1980,134 +AK,1981,144 +AK,1982,148.3 +AK,1983,146.9 +AK,1984,133.6 +AK,1985,128.9 +AK,1986,126.3 +AK,1987,110.4 +AK,1988,80 +AK,1989,94.4 +AK,1990,100.2 +AK,1991,101.8 +AK,1992,98.5 +AK,1993,95.2 +AK,1994,94.1 +AK,1995,89.5 +AK,1996,89.9 +AK,1997,81.7 +AK,1998,78.6 +AK,1999,70.5 +AK,2000,66 +AL,1970,89.8 +AL,1971,95.4 +AL,1972,101.1 +AL,1973,102.9 +AL,1974,108.2 +AL,1975,111.7 +AL,1976,116.2 +AL,1977,117.1 +AL,1978,123 +AL,1979,121.4 +AL,1980,123.2 +AL,1981,119.6 +AL,1982,119.1 +AL,1983,116.3 +AL,1984,113 +AL,1985,114.5 +AL,1986,116.3 +AL,1987,114 +AL,1988,112.1 +AL,1989,105.6 +AL,1990,108.6 +AL,1991,107.9 +AL,1992,109.1 +AL,1993,108.5 +AL,1994,107.1 +AL,1995,102.6 +AL,1996,101.4 +AL,1997,104.9 +AL,1998,106.2 +AL,1999,100.7 +AL,2000,96.2 +AR,1970,100.3 +AR,1971,104.1 +AR,1972,103.9 +AR,1973,108 +AR,1974,109.7 +AR,1975,114.8 +AR,1976,119.1 +AR,1977,122.6 +AR,1978,127.3 +AR,1979,126.5 +AR,1980,131.8 +AR,1981,128.7 +AR,1982,127.4 +AR,1983,128 +AR,1984,123.1 +AR,1985,125.8 +AR,1986,126 +AR,1987,122.3 +AR,1988,121.5 +AR,1989,118.3 +AR,1990,113.1 +AR,1991,116.8 +AR,1992,126 +AR,1993,113.8 +AR,1994,108.8 +AR,1995,113 +AR,1996,110.7 +AR,1997,108.7 +AR,1998,109.5 +AR,1999,104.8 +AR,2000,99.4 +AZ,1970,115.2 +AZ,1971,109.6 +AZ,1972,125 +AZ,1973,128.3 +AZ,1974,133.1 +AZ,1975,121.8 +AZ,1976,122.3 +AZ,1977,121.7 +AZ,1978,124.7 +AZ,1979,124.6 +AZ,1980,126.8 +AZ,1981,113.8 +AZ,1982,113.5 +AZ,1983,111.1 +AZ,1984,107.1 +AZ,1985,107.1 +AZ,1986,107.9 +AZ,1987,106.1 +AZ,1988,102.2 +AZ,1989,96.8 +AZ,1990,88.9 +AZ,1991,81.2 +AZ,1992,79 +AZ,1993,80.3 +AZ,1994,79 +AZ,1995,77.3 +AZ,1996,71.4 +AZ,1997,64.6 +AZ,1998,66 +AZ,1999,61.7 +AZ,2000,58.5 +CA,1970,123 +CA,1971,121 +CA,1972,123.5 +CA,1973,124.4 +CA,1974,126.7 +CA,1975,127.1 +CA,1976,128 +CA,1977,126.4 +CA,1978,126.1 +CA,1979,121.9 +CA,1980,120.2 +CA,1981,118.6 +CA,1982,115.4 +CA,1983,110.8 +CA,1984,104.8 +CA,1985,102.8 +CA,1986,99.7 +CA,1987,97.5 +CA,1988,90.1 +CA,1989,82.4 +CA,1990,77.8 +CA,1991,68.7 +CA,1992,67.5 +CA,1993,63.4 +CA,1994,58.6 +CA,1995,56.4 +CA,1996,54.5 +CA,1997,53.8 +CA,1998,52.3 +CA,1999,47.2 +CA,2000,41.6 +CO,1970,124.8 +CO,1971,125.5 +CO,1972,134.3 +CO,1973,137.9 +CO,1974,132.8 +CO,1975,131 +CO,1976,134.2 +CO,1977,132 +CO,1978,129.2 +CO,1979,131.5 +CO,1980,131 +CO,1981,133.8 +CO,1982,130.5 +CO,1983,125.3 +CO,1984,119.7 +CO,1985,112.4 +CO,1986,109.9 +CO,1987,102.4 +CO,1988,94.6 +CO,1989,88.8 +CO,1990,87.4 +CO,1991,90.2 +CO,1992,88.3 +CO,1993,88.6 +CO,1994,89.1 +CO,1995,85.4 +CO,1996,83.1 +CO,1997,81.3 +CO,1998,81.2 +CO,1999,79.6 +CO,2000,73 +CT,1970,120 +CT,1971,117.6 +CT,1972,110.8 +CT,1973,109.3 +CT,1974,112.4 +CT,1975,110.2 +CT,1976,113.4 +CT,1977,117.3 +CT,1978,117.5 +CT,1979,117.4 +CT,1980,118 +CT,1981,116.4 +CT,1982,114.7 +CT,1983,114.1 +CT,1984,112.5 +CT,1985,111 +CT,1986,108.5 +CT,1987,109 +CT,1988,104.8 +CT,1989,100.6 +CT,1990,91.5 +CT,1991,86.7 +CT,1992,83.6 +CT,1993,79.1 +CT,1994,76.6 +CT,1995,79.3 +CT,1996,76 +CT,1997,75.9 +CT,1998,75.5 +CT,1999,73.4 +CT,2000,71.4 +DC,1970,200.4 +DC,1971,213 +DC,1972,220.6 +DC,1973,209.4 +DC,1974,182.7 +DC,1975,176.5 +DC,1976,167.7 +DC,1977,136.7 +DC,1978,130.3 +DC,1979,128.5 +DC,1980,129.7 +DC,1981,131.9 +DC,1982,132 +DC,1983,124.5 +DC,1984,124.6 +DC,1985,122.1 +DC,1986,117.2 +DC,1987,113 +DC,1988,105.4 +DC,1989,101.1 +DC,1990,84.1 +DC,1991,82.2 +DC,1992,82.6 +DC,1993,68 +DC,1994,54.1 +DC,1995,57.9 +DC,1996,56 +DC,1997,54.3 +DC,1998,50.9 +DC,1999,50.8 +DC,2000,50.2 +DE,1970,155 +DE,1971,161.1 +DE,1972,156.3 +DE,1973,154.7 +DE,1974,151.3 +DE,1975,147.6 +DE,1976,153 +DE,1977,153.3 +DE,1978,155.5 +DE,1979,150.2 +DE,1980,150.5 +DE,1981,152.6 +DE,1982,154.1 +DE,1983,149.6 +DE,1984,144 +DE,1985,144.5 +DE,1986,142.4 +DE,1987,141 +DE,1988,137.1 +DE,1989,131.7 +DE,1990,127.2 +DE,1991,118.8 +DE,1992,120 +DE,1993,123.8 +DE,1994,126.1 +DE,1995,127.2 +DE,1996,128.3 +DE,1997,124.1 +DE,1998,132.8 +DE,1999,139.5 +DE,2000,140.7 +FL,1970,123.6 +FL,1971,123 +FL,1972,126 +FL,1973,132.6 +FL,1974,136 +FL,1975,131.9 +FL,1976,130.3 +FL,1977,133.1 +FL,1978,133.4 +FL,1979,135.4 +FL,1980,139 +FL,1981,132.8 +FL,1982,131.8 +FL,1983,127.6 +FL,1984,124.7 +FL,1985,124.2 +FL,1986,124.7 +FL,1987,125.1 +FL,1988,119.5 +FL,1989,113.8 +FL,1990,110.1 +FL,1991,98.2 +FL,1992,95.9 +FL,1993,96.8 +FL,1994,93.6 +FL,1995,95.4 +FL,1996,91.8 +FL,1997,93 +FL,1998,91.7 +FL,1999,88.2 +FL,2000,82.5 +GA,1970,109.9 +GA,1971,115.7 +GA,1972,117 +GA,1973,119.8 +GA,1974,123.7 +GA,1975,122.9 +GA,1976,125.9 +GA,1977,127.9 +GA,1978,130.6 +GA,1979,131 +GA,1980,134 +GA,1981,131.7 +GA,1982,131.2 +GA,1983,128.6 +GA,1984,126.3 +GA,1985,128.8 +GA,1986,129 +GA,1987,129.3 +GA,1988,124.1 +GA,1989,117.1 +GA,1990,113.8 +GA,1991,109.6 +GA,1992,109.2 +GA,1993,109.2 +GA,1994,107.8 +GA,1995,100.3 +GA,1996,102.7 +GA,1997,100.6 +GA,1998,100.5 +GA,1999,97.1 +GA,2000,88.1 +HI,1970,82.1 +HI,1971,85 +HI,1972,81.6 +HI,1973,79.2 +HI,1974,89.3 +HI,1975,92.4 +HI,1976,99.4 +HI,1977,95.1 +HI,1978,93.2 +HI,1979,88.8 +HI,1980,89.4 +HI,1981,85.6 +HI,1982,77.7 +HI,1983,80 +HI,1984,75.7 +HI,1985,69.8 +HI,1986,61.2 +HI,1987,54.7 +HI,1988,55.6 +HI,1989,57.4 +HI,1990,48.9 +HI,1991,54.7 +HI,1992,51.3 +HI,1993,52.4 +HI,1994,44.2 +HI,1995,45.7 +HI,1996,51.6 +HI,1997,49.1 +HI,1998,35.5 +HI,1999,32.5 +HI,2000,31.9 +IA,1970,108.5 +IA,1971,108.4 +IA,1972,109.4 +IA,1973,110.6 +IA,1974,116.1 +IA,1975,120.5 +IA,1976,124.4 +IA,1977,125.5 +IA,1978,127.1 +IA,1979,124.2 +IA,1980,124.6 +IA,1981,132.9 +IA,1982,116.2 +IA,1983,115.6 +IA,1984,111.2 +IA,1985,109.4 +IA,1986,104.1 +IA,1987,101.1 +IA,1988,100.2 +IA,1989,94.4 +IA,1990,95.4 +IA,1991,97.1 +IA,1992,95.2 +IA,1993,92.5 +IA,1994,93.4 +IA,1995,93 +IA,1996,94 +IA,1997,93.9 +IA,1998,94 +IA,1999,91.7 +IA,2000,88.9 +ID,1970,102.4 +ID,1971,108.5 +ID,1972,126.1 +ID,1973,121.8 +ID,1974,125.6 +ID,1975,123.3 +ID,1976,125.1 +ID,1977,125 +ID,1978,122.8 +ID,1979,117.5 +ID,1980,115.2 +ID,1981,114.1 +ID,1982,111.5 +ID,1983,111.3 +ID,1984,103.6 +ID,1985,100.7 +ID,1986,96.7 +ID,1987,95 +ID,1988,84.5 +ID,1989,78.4 +ID,1990,90.1 +ID,1991,85.4 +ID,1992,85.1 +ID,1993,86.7 +ID,1994,93 +ID,1995,78.2 +ID,1996,73.6 +ID,1997,75 +ID,1998,78.9 +ID,1999,75.1 +ID,2000,66.9 +IL,1970,124.8 +IL,1971,125.6 +IL,1972,126.6 +IL,1973,124.4 +IL,1974,131.9 +IL,1975,131.8 +IL,1976,134.4 +IL,1977,134 +IL,1978,136.7 +IL,1979,135.3 +IL,1980,135.2 +IL,1981,133 +IL,1982,130.7 +IL,1983,127.9 +IL,1984,124 +IL,1985,121.6 +IL,1986,118.2 +IL,1987,109.5 +IL,1988,107.6 +IL,1989,104.6 +IL,1990,94.1 +IL,1991,96.1 +IL,1992,94.8 +IL,1993,94.6 +IL,1994,85.7 +IL,1995,84.3 +IL,1996,81.8 +IL,1997,79.6 +IL,1998,80.3 +IL,1999,72.2 +IL,2000,70 +IN,1970,134.6 +IN,1971,139.3 +IN,1972,149.2 +IN,1973,156 +IN,1974,159.6 +IN,1975,162.4 +IN,1976,166.6 +IN,1977,173 +IN,1978,150.9 +IN,1979,148.9 +IN,1980,146.9 +IN,1981,148.5 +IN,1982,147.7 +IN,1983,143 +IN,1984,137.8 +IN,1985,135.3 +IN,1986,137.6 +IN,1987,134 +IN,1988,134 +IN,1989,132.5 +IN,1990,128.3 +IN,1991,127.2 +IN,1992,128.2 +IN,1993,126.8 +IN,1994,128.2 +IN,1995,135.4 +IN,1996,135.1 +IN,1997,135.3 +IN,1998,135.9 +IN,1999,133.3 +IN,2000,125.5 +KS,1970,114 +KS,1971,102.8 +KS,1972,111 +KS,1973,115.2 +KS,1974,118.6 +KS,1975,123.4 +KS,1976,127.7 +KS,1977,127.9 +KS,1978,127.1 +KS,1979,126.4 +KS,1980,127.1 +KS,1981,132 +KS,1982,130.9 +KS,1983,127.6 +KS,1984,121.7 +KS,1985,115.7 +KS,1986,109.4 +KS,1987,105.2 +KS,1988,103.2 +KS,1989,96.5 +KS,1990,94.3 +KS,1991,91.8 +KS,1992,90 +KS,1993,89.9 +KS,1994,89.1 +KS,1995,90.1 +KS,1996,88.7 +KS,1997,89.2 +KS,1998,87.6 +KS,1999,83.3 +KS,2000,79.8 +KY,1970,155.8 +KY,1971,163.5 +KY,1972,179.4 +KY,1973,201.9 +KY,1974,212.4 +KY,1975,223 +KY,1976,230.9 +KY,1977,229.4 +KY,1978,224.7 +KY,1979,214.9 +KY,1980,215.3 +KY,1981,209.7 +KY,1982,210.6 +KY,1983,201.1 +KY,1984,183.2 +KY,1985,182.4 +KY,1986,179.8 +KY,1987,171.2 +KY,1988,173.2 +KY,1989,171.6 +KY,1990,182.5 +KY,1991,170.8 +KY,1992,167.6 +KY,1993,167.6 +KY,1994,170.2 +KY,1995,175.3 +KY,1996,179 +KY,1997,186.8 +KY,1998,171.3 +KY,1999,165.3 +KY,2000,156.2 +LA,1970,115.9 +LA,1971,119.8 +LA,1972,125.3 +LA,1973,126.7 +LA,1974,129.9 +LA,1975,133.6 +LA,1976,139.6 +LA,1977,140 +LA,1978,142.7 +LA,1979,140.1 +LA,1980,143.8 +LA,1981,144 +LA,1982,143.9 +LA,1983,133.7 +LA,1984,128.9 +LA,1985,125 +LA,1986,121.2 +LA,1987,116.5 +LA,1988,110.9 +LA,1989,103.6 +LA,1990,101.5 +LA,1991,107.2 +LA,1992,108.5 +LA,1993,106.2 +LA,1994,105.3 +LA,1995,105.7 +LA,1996,106.8 +LA,1997,105.3 +LA,1998,103.2 +LA,1999,101 +LA,2000,104.3 +MA,1970,124.3 +MA,1971,121.4 +MA,1972,117.9 +MA,1973,121.2 +MA,1974,124.3 +MA,1975,126.1 +MA,1976,116.9 +MA,1977,118.9 +MA,1978,120.5 +MA,1979,118.2 +MA,1980,120.5 +MA,1981,123.2 +MA,1982,122.3 +MA,1983,119.8 +MA,1984,116 +MA,1985,117.2 +MA,1986,115.8 +MA,1987,113.8 +MA,1988,111.7 +MA,1989,104.4 +MA,1990,99.4 +MA,1991,93.5 +MA,1992,90.5 +MA,1993,84.1 +MA,1994,77 +MA,1995,76.9 +MA,1996,74.6 +MA,1997,66.7 +MA,1998,63.6 +MA,1999,60.4 +MA,2000,58.3 +MD,1970,123.5 +MD,1971,126.7 +MD,1972,133.2 +MD,1973,137.3 +MD,1974,144.8 +MD,1975,146.1 +MD,1976,137.1 +MD,1977,135.7 +MD,1978,134.1 +MD,1979,131.9 +MD,1980,131.5 +MD,1981,137.1 +MD,1982,131.3 +MD,1983,128.2 +MD,1984,122.5 +MD,1985,121.9 +MD,1986,121.8 +MD,1987,119.3 +MD,1988,115 +MD,1989,109 +MD,1990,102.8 +MD,1991,97.6 +MD,1992,91.7 +MD,1993,82.6 +MD,1994,80.8 +MD,1995,77.1 +MD,1996,75 +MD,1997,72.7 +MD,1998,72 +MD,1999,71.4 +MD,2000,57.7 +ME,1970,128.5 +ME,1971,133.2 +ME,1972,136.5 +ME,1973,138 +ME,1974,142.1 +ME,1975,140.7 +ME,1976,144.9 +ME,1977,145.6 +ME,1978,143.9 +ME,1979,138.5 +ME,1980,141.2 +ME,1981,138.9 +ME,1982,139.5 +ME,1983,135.4 +ME,1984,135.5 +ME,1985,127.9 +ME,1986,119 +ME,1987,125 +ME,1988,125 +ME,1989,122.4 +ME,1990,117.5 +ME,1991,116.1 +ME,1992,114.7 +ME,1993,108.5 +ME,1994,101.6 +ME,1995,102.3 +ME,1996,100 +ME,1997,101.1 +ME,1998,94.5 +ME,1999,85.5 +ME,2000,82.9 +MI,1970,128.6 +MI,1971,126.1 +MI,1972,131.4 +MI,1973,133.7 +MI,1974,138.8 +MI,1975,136.8 +MI,1976,138 +MI,1977,140.3 +MI,1978,141.8 +MI,1979,140.9 +MI,1980,140.7 +MI,1981,139.5 +MI,1982,137.3 +MI,1983,128.3 +MI,1984,128.2 +MI,1985,126.6 +MI,1986,126.5 +MI,1987,122.1 +MI,1988,121.1 +MI,1989,116.6 +MI,1990,114.4 +MI,1991,110.3 +MI,1992,107.6 +MI,1993,106.5 +MI,1994,112.9 +MI,1995,82.8 +MI,1996,82.9 +MI,1997,75.6 +MI,1998,73.7 +MI,1999,81.7 +MI,2000,83.7 +MN,1970,104.3 +MN,1971,116.4 +MN,1972,96.8 +MN,1973,106.8 +MN,1974,110.6 +MN,1975,111.5 +MN,1976,116.7 +MN,1977,117.2 +MN,1978,118.9 +MN,1979,118.3 +MN,1980,117.7 +MN,1981,120.8 +MN,1982,119.4 +MN,1983,113.2 +MN,1984,110.8 +MN,1985,113 +MN,1986,104.8 +MN,1987,108.8 +MN,1988,94.1 +MN,1989,92.3 +MN,1990,90.7 +MN,1991,86.2 +MN,1992,83.8 +MN,1993,81.6 +MN,1994,83.4 +MN,1995,84.1 +MN,1996,81.7 +MN,1997,84.1 +MN,1998,83.2 +MN,1999,80.7 +MN,2000,76 +MO,1970,121.3 +MO,1971,127.6 +MO,1972,130 +MO,1973,132.1 +MO,1974,135.4 +MO,1975,135.6 +MO,1976,139.5 +MO,1977,140.8 +MO,1978,141.8 +MO,1979,140.2 +MO,1980,142.1 +MO,1981,140.5 +MO,1982,139.7 +MO,1983,134.1 +MO,1984,130 +MO,1985,129.2 +MO,1986,128.8 +MO,1987,128.7 +MO,1988,127.4 +MO,1989,122.8 +MO,1990,119.1 +MO,1991,119.9 +MO,1992,122.3 +MO,1993,121.6 +MO,1994,119.4 +MO,1995,124 +MO,1996,124.1 +MO,1997,120.6 +MO,1998,120.1 +MO,1999,118 +MO,2000,113.8 +MS,1970,93.4 +MS,1971,105.4 +MS,1972,112.1 +MS,1973,115 +MS,1974,117.1 +MS,1975,116.8 +MS,1976,120.9 +MS,1977,122.1 +MS,1978,124.9 +MS,1979,123.9 +MS,1980,127 +MS,1981,125.3 +MS,1982,125.8 +MS,1983,122.3 +MS,1984,116.4 +MS,1985,115.3 +MS,1986,113.2 +MS,1987,110 +MS,1988,109 +MS,1989,108.3 +MS,1990,101.8 +MS,1991,105.6 +MS,1992,103.9 +MS,1993,105.4 +MS,1994,106 +MS,1995,107.5 +MS,1996,106.9 +MS,1997,106.3 +MS,1998,107 +MS,1999,103.9 +MS,2000,97.2 +MT,1970,111.2 +MT,1971,115.6 +MT,1972,122.2 +MT,1973,119.9 +MT,1974,121.9 +MT,1975,123.7 +MT,1976,124.9 +MT,1977,127 +MT,1978,127.2 +MT,1979,120.3 +MT,1980,122 +MT,1981,121.1 +MT,1982,122.4 +MT,1983,113.7 +MT,1984,110.1 +MT,1985,103.6 +MT,1986,97.8 +MT,1987,91.7 +MT,1988,87.1 +MT,1989,86.2 +MT,1990,84.7 +MT,1991,82.9 +MT,1992,86.6 +MT,1993,86 +MT,1994,88.2 +MT,1995,90.5 +MT,1996,87.3 +MT,1997,88.9 +MT,1998,89.1 +MT,1999,82.6 +MT,2000,75.5 +NC,1970,172.4 +NC,1971,187.6 +NC,1972,214.1 +NC,1973,226.5 +NC,1974,227.3 +NC,1975,226 +NC,1976,230.2 +NC,1977,217 +NC,1978,205.5 +NC,1979,197.3 +NC,1980,187.8 +NC,1981,179.3 +NC,1982,179 +NC,1983,169.7 +NC,1984,160.6 +NC,1985,156.3 +NC,1986,154.4 +NC,1987,150.5 +NC,1988,146 +NC,1989,139.3 +NC,1990,133.7 +NC,1991,132.7 +NC,1992,128.8 +NC,1993,129.7 +NC,1994,112.7 +NC,1995,124.9 +NC,1996,129.7 +NC,1997,125.6 +NC,1998,126 +NC,1999,113.1 +NC,2000,109 +ND,1970,93.8 +ND,1971,98.5 +ND,1972,103.8 +ND,1973,108.7 +ND,1974,110.5 +ND,1975,117.9 +ND,1976,125.4 +ND,1977,122.2 +ND,1978,121.9 +ND,1979,121.3 +ND,1980,123.7 +ND,1981,125.7 +ND,1982,126.8 +ND,1983,119.6 +ND,1984,109.4 +ND,1985,103.2 +ND,1986,99.8 +ND,1987,92.3 +ND,1988,87.1 +ND,1989,84.1 +ND,1990,77.1 +ND,1991,85.2 +ND,1992,74.4 +ND,1993,83 +ND,1994,81 +ND,1995,80.6 +ND,1996,80.8 +ND,1997,77.5 +ND,1998,79.1 +ND,1999,74.7 +ND,2000,72.5 +NE,1970,108.1 +NE,1971,108.6 +NE,1972,104.9 +NE,1973,106.6 +NE,1974,110.5 +NE,1975,114.1 +NE,1976,118.1 +NE,1977,117.7 +NE,1978,117.4 +NE,1979,116.1 +NE,1980,116.3 +NE,1981,117 +NE,1982,117.1 +NE,1983,110.8 +NE,1984,107.7 +NE,1985,105.1 +NE,1986,103.1 +NE,1987,101.3 +NE,1988,92.9 +NE,1989,93.8 +NE,1990,89.9 +NE,1991,92.4 +NE,1992,90.6 +NE,1993,91.1 +NE,1994,85.9 +NE,1995,88.5 +NE,1996,86.2 +NE,1997,85.5 +NE,1998,79.7 +NE,1999,86.8 +NE,2000,77.6 +NH,1970,265.7 +NH,1971,278 +NH,1972,296.2 +NH,1973,279 +NH,1974,269.8 +NH,1975,269.1 +NH,1976,290.5 +NH,1977,278.8 +NH,1978,269.6 +NH,1979,254.6 +NH,1980,247.8 +NH,1981,245.4 +NH,1982,239.8 +NH,1983,232.9 +NH,1984,215.1 +NH,1985,201.1 +NH,1986,195.9 +NH,1987,195.1 +NH,1988,180.4 +NH,1989,172.9 +NH,1990,152.4 +NH,1991,144.8 +NH,1992,143.7 +NH,1993,148.9 +NH,1994,153.8 +NH,1995,158.5 +NH,1996,158 +NH,1997,174.4 +NH,1998,173.8 +NH,1999,171.7 +NH,2000,147.3 +NJ,1970,120.7 +NJ,1971,125.1 +NJ,1972,128.4 +NJ,1973,120.2 +NJ,1974,121.4 +NJ,1975,122.3 +NJ,1976,122.4 +NJ,1977,122.8 +NJ,1978,123.9 +NJ,1979,123.8 +NJ,1980,124.3 +NJ,1981,127.9 +NJ,1982,126.3 +NJ,1983,121.3 +NJ,1984,117.8 +NJ,1985,116.9 +NJ,1986,115.5 +NJ,1987,111.6 +NJ,1988,107.8 +NJ,1989,103 +NJ,1990,98.6 +NJ,1991,86.7 +NJ,1992,86.6 +NJ,1993,82.9 +NJ,1994,81.6 +NJ,1995,81 +NJ,1996,79.1 +NJ,1997,77 +NJ,1998,72.6 +NJ,1999,63.6 +NJ,2000,60.1 +NM,1970,90 +NM,1971,92.6 +NM,1972,99.3 +NM,1973,98.9 +NM,1974,100.3 +NM,1975,103.1 +NM,1976,102.4 +NM,1977,102.4 +NM,1978,103.1 +NM,1979,101 +NM,1980,102.7 +NM,1981,103 +NM,1982,97.5 +NM,1983,96.3 +NM,1984,88.9 +NM,1985,88 +NM,1986,88.2 +NM,1987,82.3 +NM,1988,77.7 +NM,1989,74.4 +NM,1990,70.8 +NM,1991,69.9 +NM,1992,71.4 +NM,1993,69 +NM,1994,68.2 +NM,1995,67 +NM,1996,65.7 +NM,1997,61.8 +NM,1998,62.6 +NM,1999,59.7 +NM,2000,53.8 +NV,1970,189.5 +NV,1971,190.5 +NV,1972,198.6 +NV,1973,201.5 +NV,1974,204.7 +NV,1975,205.2 +NV,1976,201.4 +NV,1977,190.8 +NV,1978,187 +NV,1979,183.3 +NV,1980,177.7 +NV,1981,171.9 +NV,1982,165.1 +NV,1983,159.2 +NV,1984,136.6 +NV,1985,146.7 +NV,1986,142.6 +NV,1987,147.7 +NV,1988,141.9 +NV,1989,137.9 +NV,1990,137.3 +NV,1991,115.5 +NV,1992,110 +NV,1993,108.1 +NV,1994,105.2 +NV,1995,100.9 +NV,1996,99 +NV,1997,95.6 +NV,1998,102.4 +NV,1999,103.9 +NV,2000,93.2 +NY,1970,119 +NY,1971,123.3 +NY,1972,119.9 +NY,1973,118.7 +NY,1974,121.6 +NY,1975,123.9 +NY,1976,124.6 +NY,1977,125.5 +NY,1978,126.4 +NY,1979,124.4 +NY,1980,127.6 +NY,1981,130.1 +NY,1982,128.5 +NY,1983,124.6 +NY,1984,119.3 +NY,1985,115.9 +NY,1986,114 +NY,1987,108.9 +NY,1988,106.2 +NY,1989,100.1 +NY,1990,93.7 +NY,1991,86.4 +NY,1992,82.4 +NY,1993,78.2 +NY,1994,71.6 +NY,1995,70.5 +NY,1996,65.9 +NY,1997,64.5 +NY,1998,64.9 +NY,1999,62.9 +NY,2000,57.8 +OH,1970,121.6 +OH,1971,124.6 +OH,1972,124.4 +OH,1973,120.5 +OH,1974,122.1 +OH,1975,122.5 +OH,1976,124.6 +OH,1977,127.3 +OH,1978,131.3 +OH,1979,130.9 +OH,1980,133.5 +OH,1981,132.8 +OH,1982,134 +OH,1983,130 +OH,1984,127.1 +OH,1985,126.7 +OH,1986,126.3 +OH,1987,124.6 +OH,1988,122.4 +OH,1989,118.6 +OH,1990,115.5 +OH,1991,113.2 +OH,1992,112.3 +OH,1993,108.9 +OH,1994,108.6 +OH,1995,111.7 +OH,1996,107.6 +OH,1997,108.6 +OH,1998,106.4 +OH,1999,104 +OH,2000,99.9 +OK,1970,108.4 +OK,1971,115.4 +OK,1972,121.7 +OK,1973,124.1 +OK,1974,130.5 +OK,1975,132.9 +OK,1976,138.6 +OK,1977,140.4 +OK,1978,143.6 +OK,1979,141.6 +OK,1980,141.6 +OK,1981,143.7 +OK,1982,147 +OK,1983,140 +OK,1984,128.1 +OK,1985,124.2 +OK,1986,119.9 +OK,1987,113.1 +OK,1988,103.6 +OK,1989,97.5 +OK,1990,88.4 +OK,1991,87.8 +OK,1992,86.3 +OK,1993,86.2 +OK,1994,104.8 +OK,1995,109.5 +OK,1996,110.8 +OK,1997,111.8 +OK,1998,112.2 +OK,1999,111.4 +OK,2000,108.9 +OR,1970,157 +OR,1971,157.9 +OR,1972,160.4 +OR,1973,155.3 +OR,1974,155.8 +OR,1975,154.4 +OR,1976,156.5 +OR,1977,155.6 +OR,1978,156.3 +OR,1979,150.7 +OR,1980,146.5 +OR,1981,143.1 +OR,1982,135.3 +OR,1983,124.6 +OR,1984,118 +OR,1985,118.6 +OR,1986,112 +OR,1987,108.9 +OR,1988,108.2 +OR,1989,106.6 +OR,1990,102.5 +OR,1991,100.1 +OR,1992,101.6 +OR,1993,98.7 +OR,1994,96.6 +OR,1995,94.6 +OR,1996,94.3 +OR,1997,89.5 +OR,1998,84.2 +OR,1999,79.9 +OR,2000,72 +PA,1970,107.3 +PA,1971,106.3 +PA,1972,109 +PA,1973,110.7 +PA,1974,114.2 +PA,1975,114.6 +PA,1976,118.8 +PA,1977,120.1 +PA,1978,122.3 +PA,1979,122.6 +PA,1980,124 +PA,1981,125.2 +PA,1982,123.3 +PA,1983,125.3 +PA,1984,115.3 +PA,1985,115.8 +PA,1986,113.9 +PA,1987,110.6 +PA,1988,107.6 +PA,1989,107.1 +PA,1990,101.3 +PA,1991,102.5 +PA,1992,96.2 +PA,1993,94.7 +PA,1994,95.4 +PA,1995,95.4 +PA,1996,93.3 +PA,1997,92.9 +PA,1998,92.1 +PA,1999,91.1 +PA,2000,87.9 +RI,1970,123.9 +RI,1971,123.2 +RI,1972,134.4 +RI,1973,142 +RI,1974,146.1 +RI,1975,154.7 +RI,1976,150.2 +RI,1977,148.8 +RI,1978,146.8 +RI,1979,145.8 +RI,1980,149.3 +RI,1981,151.2 +RI,1982,146.3 +RI,1983,135.8 +RI,1984,136.9 +RI,1985,133.4 +RI,1986,136.3 +RI,1987,124.4 +RI,1988,138 +RI,1989,120.8 +RI,1990,101.4 +RI,1991,103.6 +RI,1992,100.1 +RI,1993,94.1 +RI,1994,91.9 +RI,1995,90.8 +RI,1996,87.5 +RI,1997,90 +RI,1998,88.7 +RI,1999,86.9 +RI,2000,83.1 +SC,1970,103.6 +SC,1971,115 +SC,1972,118.7 +SC,1973,125.5 +SC,1974,129.7 +SC,1975,130.5 +SC,1976,136.8 +SC,1977,137.2 +SC,1978,140.4 +SC,1979,135.7 +SC,1980,138.3 +SC,1981,136.1 +SC,1982,136 +SC,1983,131.1 +SC,1984,127 +SC,1985,125.4 +SC,1986,126.6 +SC,1987,126.6 +SC,1988,124.4 +SC,1989,122.4 +SC,1990,118.6 +SC,1991,121.5 +SC,1992,112.9 +SC,1993,115.2 +SC,1994,112.2 +SC,1995,109.2 +SC,1996,102.9 +SC,1997,124.5 +SC,1998,126.9 +SC,1999,109.4 +SC,2000,103.9 +SD,1970,92.7 +SD,1971,96.7 +SD,1972,103 +SD,1973,103.5 +SD,1974,108.4 +SD,1975,113.5 +SD,1976,116.7 +SD,1977,115.6 +SD,1978,116.9 +SD,1979,117.4 +SD,1980,114.7 +SD,1981,115.7 +SD,1982,113 +SD,1983,109.8 +SD,1984,105.7 +SD,1985,104.4 +SD,1986,97 +SD,1987,95.8 +SD,1988,91.9 +SD,1989,87.4 +SD,1990,88.3 +SD,1991,91.8 +SD,1992,93 +SD,1993,91.6 +SD,1994,94.8 +SD,1995,98.6 +SD,1996,92.3 +SD,1997,88.8 +SD,1998,88.3 +SD,1999,83.5 +SD,2000,75.1 +TN,1970,99.8 +TN,1971,106.3 +TN,1972,111.5 +TN,1973,109.7 +TN,1974,114.8 +TN,1975,117.4 +TN,1976,121.7 +TN,1977,124.6 +TN,1978,127.3 +TN,1979,127.2 +TN,1980,130.4 +TN,1981,129.1 +TN,1982,131.4 +TN,1983,129 +TN,1984,125.1 +TN,1985,128.7 +TN,1986,107.9 +TN,1987,130.6 +TN,1988,125.3 +TN,1989,124.7 +TN,1990,121.8 +TN,1991,120.6 +TN,1992,121 +TN,1993,120.8 +TN,1994,118.8 +TN,1995,125.4 +TN,1996,119.2 +TN,1997,118.9 +TN,1998,119.7 +TN,1999,115.8 +TN,2000,108.7 +TX,1970,106.4 +TX,1971,108.9 +TX,1972,108.6 +TX,1973,110.4 +TX,1974,114.7 +TX,1975,116 +TX,1976,121.4 +TX,1977,124.2 +TX,1978,126.6 +TX,1979,126.4 +TX,1980,129.7 +TX,1981,129 +TX,1982,131.2 +TX,1983,126.4 +TX,1984,117.2 +TX,1985,115.9 +TX,1986,113.7 +TX,1987,105.8 +TX,1988,96.5 +TX,1989,94.5 +TX,1990,85.6 +TX,1991,79.6 +TX,1992,77.2 +TX,1993,81.3 +TX,1994,78.8 +TX,1995,75.2 +TX,1996,74.6 +TX,1997,72.6 +TX,1998,73.2 +TX,1999,67.8 +TX,2000,69.3 +UT,1970,65.5 +UT,1971,67.7 +UT,1972,71.3 +UT,1973,72.7 +UT,1974,75.6 +UT,1975,75.8 +UT,1976,77.9 +UT,1977,78.9 +UT,1978,79.6 +UT,1979,79.1 +UT,1980,74.8 +UT,1981,77.6 +UT,1982,73.6 +UT,1983,69 +UT,1984,66.3 +UT,1985,66.5 +UT,1986,64.4 +UT,1987,67.7 +UT,1988,55 +UT,1989,57 +UT,1990,53.4 +UT,1991,53.5 +UT,1992,55 +UT,1993,56.2 +UT,1994,55.8 +UT,1995,52 +UT,1996,54 +UT,1997,57 +UT,1998,42.3 +UT,1999,43.9 +UT,2000,40.7 +VA,1970,124.3 +VA,1971,128.4 +VA,1972,137 +VA,1973,143.1 +VA,1974,149.6 +VA,1975,152.7 +VA,1976,158.1 +VA,1977,157.7 +VA,1978,155.9 +VA,1979,151.8 +VA,1980,148.9 +VA,1981,149.9 +VA,1982,147.4 +VA,1983,144.7 +VA,1984,136.8 +VA,1985,134.6 +VA,1986,135.8 +VA,1987,133 +VA,1988,129.5 +VA,1989,122.5 +VA,1990,118.9 +VA,1991,109.1 +VA,1992,108.2 +VA,1993,105.4 +VA,1994,106.2 +VA,1995,106.7 +VA,1996,104.6 +VA,1997,108 +VA,1998,105.6 +VA,1999,102.1 +VA,2000,96.7 +VT,1970,122.6 +VT,1971,124.4 +VT,1972,138 +VT,1973,146.8 +VT,1974,151.8 +VT,1975,155.5 +VT,1976,171.1 +VT,1977,169.4 +VT,1978,162.4 +VT,1979,160.9 +VT,1980,161.6 +VT,1981,163.8 +VT,1982,162.3 +VT,1983,153.8 +VT,1984,144.3 +VT,1985,144.5 +VT,1986,131.2 +VT,1987,128.3 +VT,1988,128.7 +VT,1989,120.9 +VT,1990,124.3 +VT,1991,120.9 +VT,1992,126.5 +VT,1993,117.2 +VT,1994,120.3 +VT,1995,123.2 +VT,1996,102.5 +VT,1997,97.7 +VT,1998,97 +VT,1999,94.1 +VT,2000,88.9 +WA,1970,96.7 +WA,1971,97 +WA,1972,88.5 +WA,1973,91 +WA,1974,98.6 +WA,1975,99.5 +WA,1976,100.3 +WA,1977,99.3 +WA,1978,101.3 +WA,1979,101.4 +WA,1980,101.4 +WA,1981,112.7 +WA,1982,106.6 +WA,1983,102.4 +WA,1984,96.2 +WA,1985,96.5 +WA,1986,94.2 +WA,1987,91.3 +WA,1988,88.2 +WA,1989,86.1 +WA,1990,83.4 +WA,1991,78.7 +WA,1992,81.1 +WA,1993,79.4 +WA,1994,72.6 +WA,1995,67.4 +WA,1996,58.1 +WA,1997,55.6 +WA,1998,56.4 +WA,1999,55.1 +WA,2000,52.8 +WI,1970,106.4 +WI,1971,105.4 +WI,1972,108.8 +WI,1973,109.5 +WI,1974,111.8 +WI,1975,113.5 +WI,1976,115.4 +WI,1977,117.2 +WI,1978,116.7 +WI,1979,117.1 +WI,1980,117.6 +WI,1981,119.9 +WI,1982,115.6 +WI,1983,106.3 +WI,1984,105.6 +WI,1985,107 +WI,1986,105.4 +WI,1987,106 +WI,1988,102.6 +WI,1989,100.3 +WI,1990,94.1 +WI,1991,95.5 +WI,1992,96.2 +WI,1993,91.2 +WI,1994,91.8 +WI,1995,93.5 +WI,1996,92.1 +WI,1997,91.9 +WI,1998,88.7 +WI,1999,84.4 +WI,2000,80.1 +WV,1970,114.5 +WV,1971,111.5 +WV,1972,117.5 +WV,1973,116.6 +WV,1974,119.9 +WV,1975,123.2 +WV,1976,129.7 +WV,1977,133.9 +WV,1978,131.6 +WV,1979,122.1 +WV,1980,122.3 +WV,1981,120.5 +WV,1982,119.8 +WV,1983,115.7 +WV,1984,111.9 +WV,1985,109.1 +WV,1986,112.1 +WV,1987,107.5 +WV,1988,109.1 +WV,1989,104 +WV,1990,104.1 +WV,1991,100.1 +WV,1992,98 +WV,1993,111 +WV,1994,104.2 +WV,1995,115.8 +WV,1996,112.7 +WV,1997,114.5 +WV,1998,114.6 +WV,1999,112.4 +WV,2000,107.9 +WY,1970,132.2 +WY,1971,131.7 +WY,1972,140 +WY,1973,141.2 +WY,1974,145.8 +WY,1975,160.7 +WY,1976,161.5 +WY,1977,160.4 +WY,1978,160.3 +WY,1979,168.6 +WY,1980,158.1 +WY,1981,163.1 +WY,1982,157.7 +WY,1983,141.2 +WY,1984,128.9 +WY,1985,125.7 +WY,1986,124.8 +WY,1987,110.4 +WY,1988,114.3 +WY,1989,111.4 +WY,1990,96.9 +WY,1991,109.1 +WY,1992,110.8 +WY,1993,108.4 +WY,1994,111.2 +WY,1995,115 +WY,1996,110.3 +WY,1997,108.8 +WY,1998,102.9 +WY,1999,104.8 +WY,2000,90.5 diff --git a/scspill/data/files/cao_dowd/state_list.csv b/scspill/data/files/cao_dowd/state_list.csv new file mode 100644 index 0000000..f1c84be --- /dev/null +++ b/scspill/data/files/cao_dowd/state_list.csv @@ -0,0 +1,52 @@ +state,treated,missing12,neighbor +CA,1,0,0 +AK,0,1,0 +AL,0,0,0 +AR,0,0,0 +AZ,0,1,1 +CO,0,0,0 +CT,0,0,0 +DC,0,1,0 +DE,0,0,0 +FL,0,1,0 +GA,0,0,0 +HI,0,1,0 +IA,0,0,0 +ID,0,0,0 +IL,0,0,0 +IN,0,0,0 +KS,0,0,0 +KY,0,0,0 +LA,0,0,0 +MA,0,1,0 +MD,0,1,0 +ME,0,0,0 +MI,0,1,0 +MN,0,0,0 +MO,0,0,0 +MS,0,0,0 +MT,0,0,0 +NC,0,0,0 +ND,0,0,0 +NE,0,0,0 +NH,0,0,0 +NJ,0,1,0 +NM,0,0,0 +NV,0,0,1 +NY,0,1,0 +OH,0,0,0 +OK,0,0,0 +OR,0,1,1 +PA,0,0,0 +RI,0,0,0 +SC,0,0,0 +SD,0,0,0 +TN,0,0,0 +TX,0,0,0 +UT,0,0,0 +VA,0,0,0 +VT,0,0,0 +WA,0,1,0 +WI,0,0,0 +WV,0,0,0 +WY,0,0,0 \ No newline at end of file diff --git a/scspill/estimators/scspill.py b/scspill/estimators/scspill.py index 6d88a1a..af04c13 100644 --- a/scspill/estimators/scspill.py +++ b/scspill/estimators/scspill.py @@ -1,6 +1,6 @@ -"""The ``sar`` model: Bayesian spatial-autoregressive spillover SCM. +"""Dispatch spillover SCM models selected by ``method``. -scspill's first model -- today its only one -- implementing: +The default ``sar`` model implements: Sakaguchi, S., & Tagawa, H. (2026). "Identification and Bayesian Inference for Synthetic Control Methods with Spillover Effects." The @@ -34,6 +34,7 @@ ) from ..utils.datautils import balance from ..utils.results_helpers import build_effect_submodels, make_weights_results +from ..utils.scspill_helpers.cd import CDConfig, CDResults, fit_cd from ..utils.scspill_helpers.config import SCSPILLConfig from ..utils.scspill_helpers.sar.inference import build_inference, prediction_interval_spec from ..utils.scspill_helpers.sar.pipeline import run_scspill @@ -42,17 +43,18 @@ class SCSPILL: - """Spillover-aware synthetic control estimator (model ``sar``). + """Spillover-aware synthetic control estimator (models ``sar`` and ``cd``). - Fits the Bayesian spatial-autoregressive spillover model of Sakaguchi & - Tagawa (2026) -- scspill's first and, today, only model. Select it with - ``method="sar"``; that is the default, so existing code needs no change. + The default ``sar`` fit is the Bayesian spatial-autoregressive model of + Sakaguchi and Tagawa (2026). ``method="cd"`` fits the Cao--Dowd (2019, + v4) all-unit SCM and end-of-sample inference. Parameters ---------- - config : SCSPILLConfig or dict + config : SCSPILLConfig, CDConfig, or dict The estimation configuration; a dict is coerced into - :class:`scspill.config_models.SCSPILLConfig`. + :class:`scspill.config_models.SCSPILLConfig` or ``CDConfig`` according + to ``method``. Examples -------- @@ -70,14 +72,20 @@ class SCSPILL: ``` """ - def __init__(self, config: SCSPILLConfig | dict) -> None: + def __init__(self, config: SCSPILLConfig | CDConfig | dict) -> None: if isinstance(config, dict): try: - config = SCSPILLConfig(**config) + config = ( + CDConfig(**config) + if config.get("method", "sar") == "cd" + else SCSPILLConfig(**config) + ) except ValidationError as exc: raise ScspillConfigError(f"Invalid SCSPILL configuration: {exc}") from exc - if not isinstance(config, SCSPILLConfig): - raise ScspillConfigError("config must be a SCSPILLConfig or a dict of its fields.") + if not isinstance(config, (SCSPILLConfig, CDConfig)): + raise ScspillConfigError( + "config must be a SCSPILLConfig, CDConfig, or a dict of fields." + ) self.config = config # Convenience attribute unpacking (mlsynth convention). self.df = config.df @@ -88,7 +96,7 @@ def __init__(self, config: SCSPILLConfig | dict) -> None: self.method = config.method self.display_graphs = config.display_graphs - def fit(self) -> SCSPILLResults: + def fit(self) -> SCSPILLResults | CDResults: """Estimate the model and return standardized results. Runs the balance check, panel preparation, the two-step sampler, the @@ -109,6 +117,12 @@ def fit(self) -> SCSPILLResults: """ cfg = self.config + if isinstance(cfg, CDConfig): + cd_result = fit_cd(cfg) + if cfg.display_graphs: + cd_result.plot(kind="auto") + return cd_result + try: balance(cfg.df, cfg.unitid, cfg.time) inputs = prepare_scspill_inputs( diff --git a/scspill/tests/test_cd.py b/scspill/tests/test_cd.py new file mode 100644 index 0000000..a308b1a --- /dev/null +++ b/scspill/tests/test_cd.py @@ -0,0 +1,263 @@ +"""Cao--Dowd estimator, inference, and Proposition 99 replication checks.""" + +import json +from pathlib import Path + +import numpy as np +import pandas as pd +import pytest + +from scspill import SCSPILL, CDConfig, CDResults, ScspillConfigError, ScspillEstimationError +from scspill.data import load_california_cd +from scspill.utils.scspill_helpers.cd import _scm_all + +REFERENCE = Path(__file__).resolve().parents[2] / "benchmarks/reference/cao_dowd" + + +def test_scm_is_optimal_with_tiny_maxiter_and_falls_back_to_fista(monkeypatch): + rng = np.random.default_rng(17) + y = rng.normal(size=(4, 12)) + y[0] = 0.85 * y[1] + 0.15 * y[2] + 0.02 * rng.normal(size=12) + short_intercepts, short_weights = _scm_all(y, maxiter=1, ftol=1e-12) + regular_intercepts, regular_weights = _scm_all(y, maxiter=3000, ftol=1e-12) + short_loss = np.sum((y[0] - short_intercepts[0] - short_weights[0] @ y) ** 2) + regular_loss = np.sum((y[0] - regular_intercepts[0] - regular_weights[0] @ y) ** 2) + assert short_loss <= regular_loss + 1e-5 + + # A failed SLSQP falls back to projected gradient with the same optimum. + from scipy.optimize import OptimizeResult + + def failing(fun, x0, **kwargs): + return OptimizeResult(x=np.full_like(x0, np.nan), success=False) + + monkeypatch.setattr("scspill.utils.scspill_helpers.cd.minimize", failing) + fb_intercepts, fb_weights = _scm_all(y, maxiter=3000, ftol=1e-12) + fb_loss = np.sum((y[0] - fb_intercepts[0] - fb_weights[0] @ y) ** 2) + assert fb_loss <= regular_loss + 1e-4 + assert fb_weights[0].sum() == pytest.approx(1) + + +@pytest.mark.parametrize("spec", [51, 41]) +def test_application_against_authors_frozen_output(spec): + panel = load_california_cd(spec) + expected = json.loads((REFERENCE / f"results_{spec}.json").read_text()) + result = SCSPILL({**panel.config_kwargs(), "display_graphs": False}).fit() + assert isinstance(result, CDResults) + assert len(result.unit_labels) == spec + assert result.t_pre == 19 + assert result.rank == expected["k"] + assert result.condition_number == pytest.approx(expected["cond"], abs=0.03) + np.testing.assert_allclose(result.effects_panel["CA"], expected["te"], atol=0.03) + np.testing.assert_allclose(result.counterfactual, expected["sp"], atol=0.03) + np.testing.assert_allclose( + result.published_ci["CA"].iloc[0].values - result.effects_panel["CA"].iloc[0], + expected["te_ci"], + atol=0.02, + ) + np.testing.assert_allclose(result.pointwise_p_values["CA"], expected["te_pvalue"]) + for state in ("AZ", "NV", "OR"): + np.testing.assert_allclose( + result.effects_panel[state], expected["u" + state][19:], atol=0.03 + ) + np.testing.assert_allclose( + result.pre_residuals[state], expected["u" + state][:19], atol=0.03 + ) + np.testing.assert_allclose(result.pointwise_p_values[state], expected["p" + state]) + + +def test_defaults_follow_the_authors_application_code(): + """Default critical values and intervals reproduce Figure 4 and Table S.6.""" + panel = load_california_cd(51) + result = SCSPILL({**panel.config_kwargs(), "display_graphs": False}).fit() + assert (result.quantile_method, result.ci_method) == ("linear", "published") + assert result.effect_ci is result.published_ci + spill = np.eye(len(result.unit_labels))[1:] + years = result.effects_panel.index + linear = [t for t in years if result.test_hypothesis(spill, period=t).reject] + paper = [ + t + for t in years + if result.test_hypothesis(spill, period=t, quantile_method="empirical").reject + ] + assert linear == [1990, 1991, 1992, *range(1994, 2001)] + assert 1994 not in paper and 1995 not in paper + # Prediction band is the published CI mapped to the counterfactual scale. + band = result.time_series.counterfactual_lower + observed = result.time_series.observed_outcome + ci = result.published_ci["CA"] + np.testing.assert_allclose(band[19:], observed[19:] - ci["upper"].to_numpy()) + + inverted = SCSPILL( + {**panel.config_kwargs(), "ci_method": "inversion", "display_graphs": False} + ).fit() + assert inverted.effect_ci is inverted.pointwise_ci + width = inverted.pointwise_ci["CA"].eval("upper - lower") + assert np.allclose(width, width.iloc[0]) # symmetric, constant half-width + + +def test_general_hypothesis_and_distinct_intervals(): + panel = load_california_cd(41) + result = SCSPILL({**panel.config_kwargs(), "display_graphs": False}).fit() + c = pd.DataFrame(0.0, index=["CA+NV"], columns=result.unit_labels) + c.loc["CA+NV", ["CA", "NV"]] = 1 + test = result.test_hypothesis(c, period=1990, d=0) + assert test.reference_statistics.shape == (19,) + assert 0 <= test.p_value <= 1 + assert test.statistic >= 0 + assert result.pointwise_ci["CA"].loc[1990, "lower"] != pytest.approx( + result.published_ci["CA"].loc[1990, "lower"] + ) + assert result.att_ci is None # a pointwise interval is not an ATT interval + assert result.time_series.has_prediction_interval + assert len(result.kappa_test) == 12 + + +def test_labeled_effect_matrix_and_rank_failure(): + panel = load_california_cd(41) + labels = list(pd.unique(panel.df["state"])) + a = pd.DataFrame(0.0, index=labels[::-1], columns=["direct", "shared"]) + a.loc["CA", "direct"] = 1 + a.loc[["AZ", "NV", "OR"], "shared"] = 1 + config = CDConfig( + **{ + **panel.config_kwargs(), + "affected_units": None, + "effect_matrix": a, + "display_graphs": False, + } + ) + result = SCSPILL(config).fit() + assert result.effect_matrix.loc["CA", "direct"] == 1 + assert result.effects_panel.loc[1990, "AZ"] == pytest.approx( + result.effects_panel.loc[1990, "NV"] + ) + bad = pd.DataFrame(np.eye(len(labels)), index=labels, columns=labels) + with pytest.raises(ScspillEstimationError, match="unidentified"): + SCSPILL( + { + **panel.config_kwargs(), + "affected_units": None, + "effect_matrix": bad, + "display_graphs": False, + } + ).fit() + + +def test_leave_one_period_out_and_input_validation(): + panel = load_california_cd(41) + small = panel.df.loc[ + panel.df["state"].isin(["CA", "AZ", "NV", "OR", "CO"]) + & panel.df["year"].between(1970, 1977) + ].copy() + small.loc[small["state"] == "CA", "treated"] = ( + small.loc[small["state"] == "CA", "year"] >= 1976 + ).astype(int) + cfg = { + **panel.config_kwargs(), + "df": small, + "affected_units": ["AZ", "NV"], + "display_graphs": False, + } + full = SCSPILL(cfg).fit() + loo = SCSPILL({**cfg, "reference_mode": "leave_one_period_out"}).fit() + assert not np.allclose(full.placebo_effects, loo.placebo_effects) + assert loo.kappa_test[1976].reference_statistics.size == 6 + with pytest.raises(ScspillConfigError): + SCSPILL({**cfg, "affected_units": ["unknown"]}).fit() + with pytest.raises(ScspillConfigError): + loo.test_hypothesis(np.ones((1, 2))) + with pytest.raises(ScspillConfigError): + loo.test_hypothesis(np.empty((0, len(loo.unit_labels)))) + with pytest.raises(ScspillConfigError): + loo.test_hypothesis(np.eye(len(loo.unit_labels))[0], d=[0, 1]) + + +def test_monte_carlo_first_replication_against_published_estimates(): + """Check all six Section S.3.2 settings against author DGP draws and estimates.""" + with np.load(REFERENCE / "inference_replication_1.npz") as fixture: + panels = fixture["panels"] + expected = fixture["sp_estimates"] + labels = [f"U{i:02d}" for i in range(51)] + for panel, expected_estimate in zip(panels, expected, strict=True): + df = pd.DataFrame(panel, index=labels, columns=range(51)).stack().rename("y").reset_index() + df.columns = ["unit", "time", "y"] + df["treated"] = ((df["unit"] == "U00") & (df["time"] == 50)).astype(int) + result = SCSPILL( + { + "method": "cd", + "df": df, + "outcome": "y", + "treat": "treated", + "unitid": "unit", + "time": "time", + "affected_units": labels[1:14], + "display_graphs": False, + } + ).fit() + assert result.effects_panel.iloc[0, 0] == pytest.approx(expected_estimate, abs=0.001) + + +def test_kappa_against_author_reference(): + """Compare range and homogeneous structure statistics with author code.""" + with np.load(REFERENCE / "inference_replication_1.npz") as fixture: + panel = fixture["panels"][1] + expected = json.loads((REFERENCE / "kappa_rep1.json").read_text()) + labels = [f"U{i:02d}" for i in range(51)] + df = pd.DataFrame(panel, index=labels, columns=range(51)).stack().rename("y").reset_index() + df.columns = ["unit", "time", "y"] + df["treated"] = ((df["unit"] == "U00") & (df["time"] == 50)).astype(int) + base = { + "method": "cd", + "df": df, + "outcome": "y", + "treat": "treated", + "unitid": "unit", + "time": "time", + "display_graphs": False, + } + range_fit = SCSPILL({**base, "affected_units": labels[1:14]}).fit() + homogeneous = pd.DataFrame(0.0, index=labels, columns=["direct", "shared"]) + homogeneous.loc["U00", "direct"] = 1 + homogeneous.loc[labels[1:14], "shared"] = 1 + homogeneous_fit = SCSPILL({**base, "effect_matrix": homogeneous}).fit() + for key, result in (("range", range_fit), ("homogeneous", homogeneous_fit)): + assert result.kappa_test[50].statistic == pytest.approx( + expected[key]["statistic"], abs=0.001 + ) + np.testing.assert_allclose( + result.kappa_test[50].reference_statistics, + expected[key]["reference"], + atol=0.001, + rtol=0, + ) + + +@pytest.mark.slow +def test_inference_leave_one_out_against_author_reference(): + """Check all 50 squared SP-A placebo statistics, not just the effect.""" + with np.load(REFERENCE / "inference_replication_1.npz") as fixture: + panel = fixture["panels"][1] # stationary factors, r=0.5, replication 1 + author = np.load(REFERENCE / "loo_reference_stationary_r05_rep1.npy") + labels = [f"U{i:02d}" for i in range(51)] + df = pd.DataFrame(panel, index=labels, columns=range(51)).stack().rename("y").reset_index() + df.columns = ["unit", "time", "y"] + df["treated"] = ((df["unit"] == "U00") & (df["time"] == 50)).astype(int) + result = SCSPILL( + { + "method": "cd", + "df": df, + "outcome": "y", + "treat": "treated", + "unitid": "unit", + "time": "time", + "affected_units": labels[1:14], + "reference_mode": "leave_one_period_out", + "display_graphs": False, + } + ).fit() + observed = result.placebo_effects["U00"].to_numpy() ** 2 + np.testing.assert_allclose(observed, author, atol=0.001, rtol=0) + # The authors' Monte Carlo code uses the paper's left-inverse critical value. + test = result.test_hypothesis(np.eye(51)[0], period=50, quantile_method="empirical") + assert test.critical_value == pytest.approx(np.sort(author)[47], abs=0.001) + assert test.p_value == float(np.mean(author >= test.statistic)) diff --git a/scspill/tests/test_config.py b/scspill/tests/test_config.py index ead48fe..e7dc93f 100644 --- a/scspill/tests/test_config.py +++ b/scspill/tests/test_config.py @@ -123,7 +123,7 @@ def test_r_argument_mapping_documented(): def test_method_defaults_to_sar(panel): - """Omitting `method` must keep selecting the only implemented model.""" + """Omitting `method` must keep selecting the SAR model.""" assert _cfg(panel).method == "sar" @@ -131,8 +131,8 @@ def test_method_accepts_sar_explicitly(panel): assert _cfg(panel, method="sar").method == "sar" -def test_planned_methods_are_rejected_not_silently_accepted(panel): - """The roadmap names must not validate: they are not implemented.""" - for planned in ("cd", "iscm", "grossi"): +def test_non_sar_methods_are_rejected_by_sar_config(panel): + """Model-specific fields must not validate against the SAR config.""" + for method in ("cd", "iscm", "grossi"): with pytest.raises(ValidationError): - _cfg(panel, method=planned) + _cfg(panel, method=method) diff --git a/scspill/utils/scspill_helpers/cd.py b/scspill/utils/scspill_helpers/cd.py new file mode 100644 index 0000000..18ba669 --- /dev/null +++ b/scspill/utils/scspill_helpers/cd.py @@ -0,0 +1,520 @@ +"""Cao--Dowd synthetic control with spillovers and end-of-sample inference. + +The defaults follow the authors' Proposition 99 replication code at commit +86004020 (Cao and Dowd, arXiv:1902.07343v4, Section 6), which is what the +paper's empirical tables report: + +* critical values use the linearly interpolated sample quantile + (``quantile_method="linear"``), where the paper's Section 4 defines the + left inverse of the empirical CDF (``"empirical"``); +* pointwise intervals add the signed pre-treatment projected-error quantiles + to the estimate (``ci_method="published"``), where inverting the squared + P test gives a symmetric interval (``"inversion"``). + +Both intervals are always stored; ``effect_ci`` is the one selected. +""" + +from __future__ import annotations + +from typing import Any, Literal + +import numpy as np +import pandas as pd +from pydantic import BaseModel, ConfigDict, Field +from scipy.optimize import minimize + +from ...config_models import ( + BaseEstimatorConfig, + BaseEstimatorResults, + InferenceResults, + MethodDetailsResults, +) +from ...exceptions import ScspillConfigError, ScspillDataError, ScspillEstimationError +from ..results_helpers import build_effect_submodels, make_weights_results + + +class CDConfig(BaseEstimatorConfig): + """Cao--Dowd configuration for one treated unit and a fixed spillover map. + + ``affected_units`` creates one effect parameter per named unit, including + the treated unit automatically. Alternatively, ``effect_matrix`` is a + DataFrame with one row per unit and one column per effect parameter; it + must contain a direct-effect column supported only on the treated unit. + The latter permits homogeneous or otherwise structured spillovers. + """ + + method: Literal["cd"] = "cd" + affected_units: list[Any] | None = None + effect_matrix: pd.DataFrame | None = None + reference_mode: Literal["full", "leave_one_period_out"] = "full" + quantile_method: Literal["linear", "empirical"] = Field( + default="linear", + description=( + "Critical-value rule for every P and kappa test. 'linear' is the " + "interpolated sample quantile of the authors' application code; " + "'empirical' is the paper's left inverse of the empirical CDF, also " + "used in the authors' Monte Carlo studies." + ), + ) + ci_method: Literal["published", "inversion"] = Field( + default="published", + description=( + "Interval reported as effect_ci. 'published' adds signed placebo " + "quantiles to the estimate (the paper's tables); 'inversion' inverts " + "the squared P test and is symmetric." + ), + ) + ci: float = Field(default=0.95, gt=0, lt=1) + scm_maxiter: int = Field(default=3000, ge=1) + scm_ftol: float = Field(default=1e-12, gt=0) + + +def _empirical_quantile(values: np.ndarray, probability: float) -> float: + """Left inverse of the empirical CDF, as defined in v4 Section 4.""" + x = np.sort(np.asarray(values, dtype=float).reshape(-1)) + return float(x[max(0, int(np.ceil(probability * len(x))) - 1)]) + + +def _critical_value(values: np.ndarray, probability: float, method: str) -> float: + """Critical value under ``quantile_method`` (see :class:`CDConfig`).""" + if method == "empirical": + return _empirical_quantile(values, probability) + if method == "linear": + return float(np.quantile(values, probability)) + raise ScspillConfigError(f"Unknown quantile_method {method!r}.") + + +class CDTestResult(BaseModel): + """A P test with its pre-intervention reference distribution.""" + + model_config = ConfigDict(arbitrary_types_allowed=True, extra="forbid") + statistic: float + critical_value: float + p_value: float + reject: bool + reference_statistics: np.ndarray + + +class CDResults(BaseEstimatorResults): + """Cao--Dowd fit, including effects, SCM equations, and P-test results.""" + + model_config = ConfigDict(arbitrary_types_allowed=True, extra="forbid") + method: Literal["cd"] = "cd" + treated_label: Any + unit_labels: tuple[Any, ...] + time_labels: np.ndarray + t_pre: int + intercepts: pd.Series + scm_matrix: pd.DataFrame + effect_matrix: pd.DataFrame + effects_panel: pd.DataFrame + spillover_panel: pd.DataFrame + pre_residuals: pd.DataFrame + placebo_effects: pd.DataFrame + pointwise_p_values: pd.DataFrame + pointwise_ci: dict[Any, pd.DataFrame] + published_ci: dict[Any, pd.DataFrame] + kappa_test: dict[Any, CDTestResult] + rank: int + condition_number: float + reference_mode: str + quantile_method: Literal["linear", "empirical"] + ci_method: Literal["published", "inversion"] + ci_level: float + + @property + def effect_ci(self) -> dict[Any, pd.DataFrame]: + """Per-unit pointwise intervals selected by ``ci_method``.""" + return self.published_ci if self.ci_method == "published" else self.pointwise_ci + + def _critical(self, refs: np.ndarray, quantile_method: str | None) -> float: + return _critical_value(refs, self.ci_level, quantile_method or self.quantile_method) + + def _parse_c_d( + self, C: np.ndarray | pd.DataFrame, d: np.ndarray | float | None + ) -> tuple[np.ndarray, np.ndarray]: + """Validate the contrast matrix ``C`` and null vector ``d``.""" + if isinstance(C, pd.DataFrame): + if set(C.columns) != set(self.unit_labels) or C.columns.has_duplicates: + raise ScspillConfigError("C columns must match the unit labels exactly.") + c = C.loc[:, list(self.unit_labels)].to_numpy(dtype=float) + else: + c = np.atleast_2d(np.asarray(C, dtype=float)) + if ( + c.ndim != 2 + or c.shape[0] == 0 + or c.shape[1] != len(self.unit_labels) + or not np.isfinite(c).all() + ): + raise ScspillConfigError("C must be a finite matrix with one column per unit.") + q = c.shape[0] + try: + null = np.zeros(q) if d is None else np.broadcast_to(np.asarray(d, dtype=float), (q,)) + except ValueError as exc: + raise ScspillConfigError("d must have one entry per row of C.") from exc + if not np.isfinite(null).all(): + raise ScspillConfigError("d must be finite.") + return c, null + + def test_hypothesis( + self, + C: np.ndarray | pd.DataFrame, + *, + d: np.ndarray | float | None = None, + period: Any | None = None, + weight: np.ndarray | Literal["inverse_covariance"] | None = None, + quantile_method: Literal["empirical", "linear"] | None = None, + ) -> CDTestResult: + """Test ``C @ effect_t = d`` by the pre-intervention P distribution. + + ``C`` has one column per unit in ``unit_labels`` order, or labeled + columns when passed as a DataFrame. ``period`` is a post-period label. + The default statistic is the squared Euclidean norm; an explicit + positive-definite weight or the inverse placebo second moment is allowed. + ``quantile_method`` overrides the fit's critical-value rule. + """ + if period is None: + period = self.effects_panel.index[0] + if period not in self.effects_panel.index: + raise ScspillConfigError(f"Unknown post-treatment period {period!r}.") + c, null = self._parse_c_d(C, d) + q = c.shape[0] + pre = c @ self.placebo_effects.to_numpy(dtype=float).T + if weight is None: + w = np.eye(q) + elif isinstance(weight, str): + if weight != "inverse_covariance" or pre.shape[1] <= q: + raise ScspillConfigError( + "Inverse covariance requires more reference periods than rows of C." + ) + second_moment = pre @ pre.T / pre.shape[1] + try: + w = np.linalg.inv(np.atleast_2d(second_moment)) + except np.linalg.LinAlgError as exc: + raise ScspillConfigError("Projected placebo second moment is singular.") from exc + else: + w = np.asarray(weight, dtype=float) + if w.shape != (q, q) or not np.isfinite(w).all() or not np.allclose(w, w.T): + raise ScspillConfigError("weight must be a finite symmetric square matrix.") + if np.linalg.eigvalsh(w).min() <= 0: + raise ScspillConfigError("weight must be positive definite.") + observed = c @ self.effects_panel.loc[period].to_numpy(dtype=float) - null + statistic = float(observed @ w @ observed) + refs = np.einsum("it,ij,jt->t", pre, w, pre) + critical = self._critical(refs, quantile_method) + return CDTestResult( + statistic=statistic, + critical_value=critical, + p_value=float(np.mean(refs >= statistic)), + reject=bool(statistic > critical), + reference_statistics=refs, + ) + + +def _prepare( + config: CDConfig, +) -> tuple[np.ndarray, tuple[Any, ...], np.ndarray, int, Any, pd.DataFrame]: + """Validate and arrange the balanced panel as units by periods.""" + df = config.df + if df[[config.unitid, config.time]].isna().any().any(): + raise ScspillDataError("Unit and time labels cannot be missing.") + if df.duplicated([config.unitid, config.time]).any(): + raise ScspillDataError("Panel has duplicate unit-time observations.") + if df[[config.outcome, config.treat]].isna().any().any(): + raise ScspillDataError("Outcomes and treatment flags must be complete.") + units = list(pd.unique(df[config.unitid])) + times = np.asarray(sorted(pd.unique(df[config.time]))) + y = df.pivot(index=config.unitid, columns=config.time, values=config.outcome).reindex( + index=units, columns=times + ) + treat = df.pivot(index=config.unitid, columns=config.time, values=config.treat).reindex( + index=units, columns=times + ) + if y.isna().any().any() or treat.isna().any().any(): + raise ScspillDataError("Cao--Dowd requires a balanced panel.") + if not np.isfinite(y.to_numpy(dtype=float)).all(): + raise ScspillDataError("Outcomes must be finite.") + flag = treat.to_numpy(dtype=float) + if not np.isin(flag, [0, 1]).all(): + raise ScspillDataError("Treatment flags must be binary.") + treated = np.flatnonzero(flag.sum(axis=1)) + if len(treated) != 1: + raise ScspillDataError("Exactly one treated unit is required.") + idx = int(treated[0]) + t0 = int(np.argmax(flag[idx] == 1)) + if t0 < 3 or t0 >= len(times) or not np.all(flag[idx, t0:] == 1): + raise ScspillDataError("Treatment must begin after at least three pre periods and persist.") + label = units.pop(idx) + units.insert(0, label) + yy = y.loc[units].to_numpy(dtype=float) + if len(units) < 3: + raise ScspillDataError("At least three units are required for all-unit SCM.") + if config.effect_matrix is not None and config.affected_units is not None: + raise ScspillConfigError("Choose effect_matrix or affected_units, not both.") + if config.effect_matrix is not None: + a = config.effect_matrix + if a.index.has_duplicates or a.columns.has_duplicates or set(a.index) != set(units): + raise ScspillConfigError("effect_matrix needs unique rows matching every unit label.") + a = a.loc[units].astype(float) + if not np.isfinite(a.to_numpy()).all() or a.shape[1] == 0: + raise ScspillConfigError("effect_matrix must be finite and nonempty.") + direct = np.zeros(len(units)) + direct[0] = 1 + if not np.any(np.all(np.isclose(a.to_numpy().T, direct), axis=1)): + raise ScspillConfigError( + "effect_matrix must contain a direct-effect column for the treated unit." + ) + else: + affected = [] if config.affected_units is None else list(config.affected_units) + if len(set(affected)) != len(affected) or any(x not in units for x in affected): + raise ScspillConfigError("affected_units must be distinct known unit labels.") + selected = [label] + [x for x in affected if x != label] + a = pd.DataFrame(0.0, index=units, columns=selected) + for x in selected: + a.loc[x, x] = 1.0 + return yy, tuple(units), times, t0, label, a + + +def _scm_all(y_pre: np.ndarray, *, maxiter: int, ftol: float) -> tuple[np.ndarray, np.ndarray]: + """Fit an intercept and simplex donor weights for every unit.""" + n = y_pre.shape[0] + intercepts = np.zeros(n) + b_matrix = np.zeros((n, n)) + for i in range(n): + target = y_pre[i] + donors = np.delete(y_pre, i, axis=0).T + centered_y = target - target.mean() + centered_x = donors - donors.mean(axis=0) + + # Scale-free objective with analytic derivatives. SLSQP is the primary + # solve; projected gradient is the fallback if it fails or is infeasible. + scale = max(float(np.std(centered_y)), 1.0) + scaled_y = centered_y / scale + scaled_x = centered_x / scale + + def objective( + candidate: np.ndarray, + target: np.ndarray = scaled_y, + donor_matrix: np.ndarray = scaled_x, + ) -> float: + return float(np.sum((target - donor_matrix @ candidate) ** 2)) + + def gradient( + candidate: np.ndarray, + target: np.ndarray = scaled_y, + donor_matrix: np.ndarray = scaled_x, + ) -> np.ndarray: + return 2 * donor_matrix.T @ (donor_matrix @ candidate - target) + + fit = minimize( + objective, + np.full(n - 1, 1 / (n - 1)), + jac=gradient, + method="SLSQP", + bounds=[(0, None)] * (n - 1), + constraints=[ + {"type": "eq", "fun": lambda b: b.sum() - 1, "jac": lambda b: np.ones(n - 1)} + ], + options={"maxiter": maxiter, "ftol": ftol}, + ) + weights = np.clip(fit.x, 0, None) + if ( + not fit.success + or not np.isfinite(weights).all() + or weights.sum() <= 0 + or abs(weights.sum() - 1) > 1e-4 + ): + weights = _simplex_fista(centered_y, centered_x) + weights /= weights.sum() + intercepts[i] = target.mean() - donors.mean(axis=0) @ weights + b_matrix[i] = np.insert(weights, i, 0.0) + return intercepts, b_matrix + + +def _simplex_project(values: np.ndarray) -> np.ndarray: + """Project one vector onto the nonnegative unit simplex.""" + ordered = np.sort(values)[::-1] + thresholds = (np.cumsum(ordered) - 1) / np.arange(1, len(values) + 1) + active = np.flatnonzero(ordered > thresholds) + if active.size == 0: + raise ScspillEstimationError("Simplex projection has no active coordinate.") + projected = np.maximum(values - thresholds[active[-1]], 0) + return projected / projected.sum() + + +def _simplex_fista(target: np.ndarray, donors: np.ndarray) -> np.ndarray: + """Solve the centered SCM problem by restarted projected gradient.""" + n = donors.shape[1] + current = np.full(n, 1 / n) + extrapolated = current.copy() + momentum = 1.0 + lipschitz = max(2 * float(np.linalg.norm(donors, ord=2)) ** 2, np.finfo(float).eps) + previous_loss = float("inf") + for _ in range(50_000): + gradient = 2 * donors.T @ (donors @ extrapolated - target) + candidate = _simplex_project(extrapolated - gradient / lipschitz) + loss = float(np.sum((target - donors @ candidate) ** 2)) + if loss > previous_loss + 1e-13: + extrapolated = current.copy() + momentum = 1.0 + gradient = 2 * donors.T @ (donors @ extrapolated - target) + candidate = _simplex_project(extrapolated - gradient / lipschitz) + loss = float(np.sum((target - donors @ candidate) ** 2)) + candidate_gradient = 2 * donors.T @ (donors @ candidate - target) + fixed_point = _simplex_project(candidate - candidate_gradient / lipschitz) + if np.max(np.abs(candidate - fixed_point)) <= 1e-9: + return candidate + next_momentum = (1 + np.sqrt(1 + 4 * momentum**2)) / 2 + extrapolated = candidate + (momentum - 1) / next_momentum * (candidate - current) + current = candidate + momentum = next_momentum + previous_loss = loss + raise ScspillEstimationError("Simplex SCM did not converge after 50,000 fallback iterations.") + + +def _projection(b: np.ndarray, a: np.ndarray) -> tuple[np.ndarray, np.ndarray, int, float]: + """Build the effect map G and residual-space projector after rank checks.""" + r = np.eye(len(b)) - b + ra = r @ a + singular = np.linalg.svd(ra, compute_uv=False) + tol = max(ra.shape) * np.finfo(float).eps * singular[0] + rank = int(np.sum(singular > tol)) + if rank != a.shape[1]: + raise ScspillEstimationError( + f"Cao--Dowd effects unidentified: rank((I-B)A)={rank} < {a.shape[1]}." + ) + pinv = np.linalg.pinv(ra) + g = a @ pinv + projector = ra @ pinv + # The paper reports the condition number of A'(I-B)'(I-B)A, whose + # eigenvalue ratio is the square of the singular-value ratio of (I-B)A. + return g, projector, rank, float((singular[0] / singular[-1]) ** 2) + + +def fit_cd(config: CDConfig) -> CDResults: + """Fit Cao--Dowd and compute pointwise and specification-test references.""" + y, labels, times, t0, treated, a_df = _prepare(config) + a_matrix = a_df.to_numpy(dtype=float) + intercepts, b = _scm_all(y[:, :t0], maxiter=config.scm_maxiter, ftol=config.scm_ftol) + r = np.eye(len(labels)) - b + g, projector, rank, condition = _projection(b, a_matrix) + residuals = r @ y[:, :t0] - intercepts[:, None] + post_target = r @ y[:, t0:] - intercepts[:, None] + alpha = g @ post_target + if config.reference_mode == "full": + placebo = g @ residuals + kappa_pre_vectors = (np.eye(len(labels)) - projector) @ residuals + else: + placebo = np.empty((len(labels), t0)) + kappa_pre_vectors = np.empty((len(labels), t0)) + for held in range(t0): + kept = np.arange(t0) != held + fold_a, fold_b = _scm_all( + y[:, :t0][:, kept], maxiter=config.scm_maxiter, ftol=config.scm_ftol + ) + fold_g, fold_projector, _, _ = _projection(fold_b, a_matrix) + fold_resid = (np.eye(len(labels)) - fold_b) @ y[:, held] - fold_a + placebo[:, held] = fold_g @ fold_resid + kappa_pre_vectors[:, held] = (np.eye(len(labels)) - fold_projector) @ fold_resid + kappa_post_vectors = (np.eye(len(labels)) - projector) @ post_target + kappa_ref = np.linalg.norm(kappa_pre_vectors, axis=0) + kappa_post = np.linalg.norm(kappa_post_vectors, axis=0) + post_times = times[t0:] + effects = pd.DataFrame(alpha.T, index=post_times, columns=labels) + placebo_df = pd.DataFrame(placebo.T, index=times[:t0], columns=labels) + residual_df = pd.DataFrame(residuals.T, index=times[:t0], columns=labels) + pointwise = {} + published = {} + pvalues = pd.DataFrame(index=post_times, columns=labels, dtype=float) + for j, label in enumerate(labels): + pre = placebo[j] + threshold = np.sqrt(_critical_value(pre**2, config.ci, config.quantile_method)) + pointwise[label] = pd.DataFrame( + {"lower": alpha[j] - threshold, "upper": alpha[j] + threshold}, index=post_times + ) + qlo, qhi = np.quantile(pre, [(1 - config.ci) / 2, (1 + config.ci) / 2]) + published[label] = pd.DataFrame( + {"lower": alpha[j] + qlo, "upper": alpha[j] + qhi}, index=post_times + ) + pvalues[label] = np.mean(pre[:, None] ** 2 >= alpha[j][None, :] ** 2, axis=0) + kappa = {} + critical_kappa = _critical_value(kappa_ref, config.ci, config.quantile_method) + for period, statistic in zip(post_times, kappa_post, strict=True): + kappa[period] = CDTestResult( + statistic=float(statistic), + critical_value=critical_kappa, + p_value=float(np.mean(kappa_ref >= statistic)), + reject=bool(statistic > critical_kappa), + reference_statistics=kappa_ref.copy(), + ) + selected_ci = published if config.ci_method == "published" else pointwise + cf = np.r_[intercepts[0] + b[0] @ y[:, :t0], y[0, t0:] - alpha[0]] + pre_spill = np.zeros((t0, len(labels) - 1)) + spill = pd.DataFrame(np.vstack([pre_spill, alpha[1:].T]), index=times, columns=list(labels[1:])) + weights = make_weights_results( + {labels[j]: b[0, j] for j in range(1, len(labels))}, + "intercept-adjusted simplex SCM", + ) + submodels = build_effect_submodels( + observed_outcome=y[0], + counterfactual_outcome=cf, + n_pre_periods=t0, + n_post_periods=len(times) - t0, + time_periods=times, + intervention_time=times[t0], + weights=weights, + inference=InferenceResults( + method=f"cao_dowd_P_{config.reference_mode}", + confidence_level=config.ci, + details="Pointwise and joint inference is in CDResults; no ATT interval is implied.", + ), + prediction_interval={ + "lower": y[0, t0:] - selected_ci[treated]["upper"].to_numpy(), + "upper": y[0, t0:] - selected_ci[treated]["lower"].to_numpy(), + "periods": post_times, + "level": config.ci, + "kind": ( + "published_signed_quantile" + if config.ci_method == "published" + else "pointwise_P_inversion" + ), + }, + ) + return CDResults( + **submodels, + method_details=MethodDetailsResults( + method_name="SCSPILL/cd", + is_recommended=True, + parameters_used={ + "method": "cd", + "reference_mode": config.reference_mode, + "quantile_method": config.quantile_method, + "ci_method": config.ci_method, + "ci": config.ci, + "scm_maxiter": config.scm_maxiter, + "scm_ftol": config.scm_ftol, + }, + ), + treated_label=treated, + unit_labels=labels, + time_labels=times, + t_pre=t0, + intercepts=pd.Series(intercepts, index=labels), + scm_matrix=pd.DataFrame(b, index=labels, columns=labels), + effect_matrix=a_df, + effects_panel=effects, + spillover_panel=spill, + pre_residuals=residual_df, + placebo_effects=placebo_df, + pointwise_p_values=pvalues, + pointwise_ci=pointwise, + published_ci=published, + kappa_test=kappa, + rank=rank, + condition_number=condition, + reference_mode=config.reference_mode, + quantile_method=config.quantile_method, + ci_method=config.ci_method, + ci_level=config.ci, + plot_config=config.resolved_plot(), + ) diff --git a/scspill/utils/scspill_helpers/config.py b/scspill/utils/scspill_helpers/config.py index 372c371..a407184 100644 --- a/scspill/utils/scspill_helpers/config.py +++ b/scspill/utils/scspill_helpers/config.py @@ -11,7 +11,7 @@ class SCSPILLConfig(BaseEstimatorConfig): - """Configuration for a spillover synthetic-control fit (today, model ``sar``). + """Configuration for the Bayesian ``sar`` spillover synthetic-control fit. Implements: @@ -72,9 +72,8 @@ class SCSPILLConfig(BaseEstimatorConfig): "spatial-autoregressive Bayesian SCM: donor outcomes follow a SAR " "process under the supplied weights (spatial_w, spatial_W), and a " "two-step sampler estimates the synthetic weights alpha and the " - "spillover intensity rho. It is the only implemented model, and the " - "default, so omitting this field is correct. Models added later " - "select here rather than through a new class; see " + "spillover intensity rho. It remains the default; Cao--Dowd uses " + "a separate CDConfig selected by method='cd'. See " "https://quarcs-lab.github.io/scspill/models/ for the catalogue." ), ) diff --git a/scspill/utils/scspill_helpers/structures.py b/scspill/utils/scspill_helpers/structures.py index 55d459c..874459b 100644 --- a/scspill/utils/scspill_helpers/structures.py +++ b/scspill/utils/scspill_helpers/structures.py @@ -289,7 +289,7 @@ class SCSPILLResults(BaseEstimatorResults): :attr:`method` records which spillover model produced the fit. It mirrors ``SCSPILLConfig.method`` and the suffix of ``method_details.method_name`` (``"SCSPILL/"``). Only ``"sar"`` - exists today; a model added later declares its own posterior blocks + uses this result; another model declares its own estimation blocks alongside :attr:`alpha_posterior` / :attr:`sar_posterior`, which are SAR-specific, while :attr:`inputs`, :attr:`effects_detail` and the flat effect surface are the cross-model contract. diff --git a/tools/build_llms_txt.py b/tools/build_llms_txt.py index 6f8a260..cd89fe9 100644 --- a/tools/build_llms_txt.py +++ b/tools/build_llms_txt.py @@ -27,6 +27,7 @@ ("Get started (California Prop 99)", "get-started.qmd", "get-started.html"), ("Models: catalogue and roadmap", "models/index.qmd", "models/index.html"), ("The `sar` model (Sakaguchi & Tagawa)", "models/sar.qmd", "models/sar.html"), + ("The `cd` model (Cao & Dowd)", "models/cd.qmd", "models/cd.html"), ("Sudan case study", "sudan.qmd", "sudan.html"), ("Validating the sampler", "articles/validation.qmd", "articles/validation.html"), ( @@ -40,7 +41,11 @@ ] API_GROUPS = [ - ("Top level", "scspill", ["SCSPILL", "SCSPILLConfig", "SCSPILLResults"]), + ( + "Top level", + "scspill", + ["SCSPILL", "SCSPILLConfig", "SCSPILLResults", "CDConfig", "CDResults"], + ), ("scspill.validation", "scspill.validation", None), ("scspill.simulate", "scspill.simulate", None), ("scspill.data", "scspill.data", None), @@ -60,13 +65,14 @@ def canonical() -> str: 'Install: `pip install scspill` (or `pip install "scspill[numba]"` for', "JIT-compiled samplers).", "", - 'Models: exactly ONE is implemented -- `method="sar"` (Sakaguchi & Tagawa,', - "The Econometrics Journal 2026), class SCSPILL, and the default. It routes", + 'Models: TWO are implemented -- `method="sar"` (Sakaguchi & Tagawa,', + 'The Econometrics Journal 2026), the default, and `method="cd"` (Cao & Dowd,', + "arXiv:1902.07343v4). `sar` routes", "spillovers through user-supplied spatial weights scaled by one estimated", "intensity rho, via two-step MCMC (horseshoe synthetic weights, then a SAR", - "spillover block). The catalogue also lists `cd` (Cao & Dowd), `iscm` (Di", - "Stefano & Mellace) and `grossi` (Grossi et al.) as PLANNED: they are NOT", - "implemented and SCSPILLConfig rejects them. Do not emit code using them.", + "spillover block). `cd` fits intercept-adjusted all-unit SCM with a fixed", + "exposure matrix and P-test inference. `iscm` (Di Stefano & Mellace) and", + "`grossi` (Grossi et al.) remain planned and are not implemented.", f"Catalogue: {BASE}/models/index.html", "", "Minimal contract (model `sar`): SCSPILL({df, outcome, treat, unitid, time,", @@ -76,6 +82,13 @@ def canonical() -> str: "treated unit and date are inferred from the 0/1 `treat` column; spatial", "weights align by donor label.", "", + 'Minimal contract (model `cd`): SCSPILL({method: "cd", df, outcome, treat,', + "unitid, time, affected_units? or effect_matrix?, reference_mode?,", + "quantile_method?, ci_method?}).fit() -> CDResults with .effects_panel,", + ".effect_ci (= .published_ci by default, or .pointwise_ci),", + ".pointwise_p_values, .kappa_test, and .test_hypothesis(C). The 51-unit", + "1989-onward replication panel is load_california_cd(51).", + "", "## Docs", "", ]