panelkit

Fast, from-scratch causal-inference estimators for panel / geo experiments — written in Rust, exposed to Python.

panelkit reimplements the standard panel causal-inference toolbox on top of its own dependency-free numerical core (no BLAS/LAPACK, no ndarray, no rand), so the whole stack is self-contained, deterministic, and fast. The numerical core (panelkit-linalg) is a standalone crate intended to also back sibling projects (e.g. a future time-series library).

• Fast: ~45–70× a NumPy+SciPy synthetic control on a single fit, ~1400× on a full placebo test (multithreaded). Details.
• Self-contained: the numerical core is hand-written — matmul, Cholesky, QR, a one-sided Jacobi SVD, simplex solvers, and a PRNG, with zero numeric deps.
• Reproducible: all resampling inference is bit-identical regardless of thread count (deterministic per-replicate seed substreams).
• Modern: correct under staggered adoption (Callaway-Sant'Anna, Sun-Abraham) with a Goodman-Bacon diagnostic, plus a novel conformal-pooled SC family.

Install

pip install panelkit

# Or from a clone (needs a Rust toolchain — https://rustup.rs — and maturin):
pip install maturin numpy
maturin develop --release --manifest-path crates/pypanelkit/Cargo.toml

Data model — what your data should look like

Every estimator takes a single N × T NumPy array Y: one row per unit (store, region, DMA, country…), one column per time period (week, month…), in chronological order. Y[i, t] is the outcome for unit i at period t. The panel must be balanced (no missing cells) — fill or aggregate gaps before fitting. Treatment timing is passed separately, not encoded in Y.

import numpy as np

#            t=0    t=1    t=2    t=3    t=4    t=5   ...   (periods →)
Y = np.array([
    [102.0, 104.0, 101.0, 130.0, 133.0, 131.0],   # unit 0  ← treated
    [ 98.0,  99.0,  97.0,  98.0, 100.0,  99.0],    # unit 1  (control/donor)
    [201.0, 205.0, 199.0, 204.0, 206.0, 203.0],    # unit 2  (control/donor)
    # ... one row per unit ...
])
# Here T=6 periods; say treatment starts at period 3 (so periods 0–2 are
# "pre", periods 3–5 are "post"). Unit 0 jumps ~+30 post-treatment.

Two ways to declare treatment:

• Block treatment (SC family, MC-NNM, CP-ASC) — a list of treated unit indices and the first post-treatment column:

  SyntheticControl().fit(Y, treated=[0], treat_time=3)
  

• Staggered adoption (DiD family) — a per-unit treat_start array giving each unit's first-treated period; -1 (or None) means never-treated. A start at/after the last column means the unit never becomes treated within the sample, so it is used as a never-treated control (the R did package convention):

  # unit 0 treated at t=3, unit 2 at t=4, unit 1 never treated
  CallawaySantAnna().fit(Y, treat_start=[3, -1, 4])
  

Coming from a long/tidy DataFrame (unit, period, outcome rows)? Pivot to the matrix first:

Y = df.pivot(index="unit", columns="period", values="outcome").to_numpy()

Common gotchas: rows must be units and columns periods (not transposed); periods must be in time order; the array is float64 (panelkit copies/converts as needed); and at least a couple of pre-treatment periods are required (more is better — SC/SDID lean on pre-treatment fit).

Estimators at a glance

Examples

Synthetic Control (+ placebo inference)

import numpy as np
from panelkit import SyntheticControl

# Y: 50 units × 60 periods; unit 0 treated from period 45.
res = SyntheticControl(inference="placebo").fit(Y, treated=[0], treat_time=45)

res.att                 # average post-treatment effect
res.att_path            # per-period effects (length T_post)
res.counterfactual      # synthetic control's predicted path
res.weights             # donor weights (on the simplex)
res.donor_ids           # which units those weights correspond to
res.p_value             # in-space placebo p-value
print(res.summary())

Synthetic DiD — the robust default

from panelkit import SyntheticDiD

res = SyntheticDiD().fit(Y, treated=[0], treat_time=45)
print(res.att)          # unit + time weighted 2×2 DiD

Augmented SC, demeaned SC, robust SC and MC-NNM

from panelkit import AugmentedSC, DemeanedSC, MCNNM, RobustSC

AugmentedSC().fit(Y, treated=[0], treat_time=45).att          # ridge-corrected SC
DemeanedSC().fit(Y, treated=[0], treat_time=45).att           # Ferman-Pinto: SC on deviations
RobustSC().fit(Y, treated=[0], treat_time=45).att             # spectral de-noising, extrapolates
MCNNM().fit(Y, treated=[0], treat_time=45).att                # low-rank completion, λ by CV

DemeanedSC is the one to reach for when the treated unit's level sits outside the donor hull — plain SC spends its weights trying to close a gap it cannot close, and reports the remainder as effect.

Bounding donor weights

# no single donor may carry more than 20% of the synthetic unit …
SyntheticControl(max_weight=0.20).fit(Y, treated=[0], treat_time=45)
# … and every donor carries at least 1% (forced diversification)
SyntheticDiD(min_weight=0.01, max_weight=0.20).fit(Y, treated=[0], treat_time=45)

# Abadie & L'Hour (2021): prefer donors that individually resemble the treated unit
SyntheticControl(penalty=0.5).fit(Y, treated=[0], treat_time=45)

Weights stay on the simplex; the box is solved exactly (accelerated projected gradient with an active-set polish) rather than by clipping after the fact. max_weight is the guard against a counterfactual that is really just one donor market. The same options are on GeoDesign.power/evaluate/select_markets.

CP-ASC — conformal pooled SC (multiple treated units)

from panelkit import CPASC

treated = [0, 1, 2, 3, 4, 5]
res = CPASC(mode="mspe").fit(Y, treated, treat_time=22)   # CP-ASC
res.att                 # empirical-Bayes pooled ATT
res.p_value             # conformal block-permutation p-value
res.unit_att            # per-unit effects
res.unit_weight         # inverse-MSPE pooling weights

CPASC(mode="stratified").fit(Y, treated, 22)   # Strat-CP-ASC (size-robust)
CPASC(mode="cumulative").fit(Y, treated, 22)   # C-AS-CP-ASC (total-dollar target)

Difference-in-differences with staggered adoption

from panelkit import TWFE, CallawaySantAnna, SunAbraham, GoodmanBacon

# treat_start[i] = first treated period for unit i, or -1 if never treated.
cs = CallawaySantAnna().fit(Y, treat_start)
cs.att                  # overall ATT (cell-weighted "simple" aggregation)
cs.overall_group_att    # cohort-weighted "group" ATT (C&S's headline)
cs.event_time           # relative event times, e.g. [-5,...,-1, 0, 1,...]
cs.event_att            # event-study coefficients (clean pre-trends + dynamics)
cs.event_se             # influence-function standard errors
print(cs.summary())

# the C&S-recommended joint inference + a publication-style figure:
cs = CallawaySantAnna(inference="bootstrap", anticipation=0).fit(Y, treat_start)
cs.event_bands          # simultaneous sup-t bands (multiplier bootstrap)
cs.plot("event_study.png")

sa = SunAbraham().fit(Y, treat_start)           # interaction-weighted event study
twfe = TWFE().fit(Y, treat_start)               # baseline; biased under heterogeneity

# Goodman-Bacon: why TWFE is biased — decompose it into 2×2 comparisons.
bacon = GoodmanBacon().fit(Y, treat_start)
bacon.twfe              # == TWFE coefficient (Σ weightᵢ · estimateᵢ)
bacon.forbidden_weight  # weight on "already-treated as control" comparisons
print(bacon.summary())

Runnable scripts live in examples/: sc_demo.py, did_demo.py, cpasc_demo.py. See GUIDE.md for the estimand, assumptions, and valid inference for each estimator.

Geo test design (power analysis & market selection)

panelkit.design is the planning layer in front of a geo experiment — multi-method and robustness-first, with the heavy simulation in Rust. It answers: which markets should I treat, how big a lift can I detect, can I trust this design — and, once it's run, how big was the effect?

from panelkit.design import GeoDesign

# from a long/tidy DataFrame (location, date, outcome) — or GeoDesign(Y, names=…)
design = GeoDesign.from_long(df, location="dma", time="week", outcome="sales")

rep = design.power(treated=["chicago", "denver"], test_len=8, alpha=0.10)
print(rep.summary())          # plain-English report: MDE, confidence, warnings
rep.plot("design.png")        # the figure below

# guardrails: is this design trustworthy? (pre-fit, seasonality, holdout, warnings)
guard = design.diagnose(treated=["chicago", "denver"], test_len=8)
print(guard.summary())
guard.plot("guardrails.png")  # the guardrails figure below

# let it pick the markets for you:
ranked = design.select_markets(test_len=8, target_lift=0.05, max_treated=3)

# run several disjoint treatment cells at once (each vs. a shared donor pool):
mc = design.multi_cell(cells={"west": ["los_angeles", "san_diego"],
                              "east": ["boston", "philadelphia"]}, test_len=8)
print(mc.summary())           # per-cell MDE / confidence / holdout
mc.plot("multicell.png")      # the multi-cell figure below

# pin in must-have markets, drop ones you don't trust:
ranked = design.select_markets(test_len=8, target_lift=0.05, max_treated=3,
                               include=["chicago"], exclude=["miami"])

# already ran the test? measure it (SC/ASC/SDID/FP/RSC + a weighted-average ensemble):
ev = design.evaluate(treated=["chicago", "denver"], treat_start=52)
print(ev.summary())           # per-method + ensemble lift, CI, cumulative
ev.plot("evaluate.png")       # observed vs counterfactual + lift-by-method
ev.plot_effect_over_time("effect.png")   # pointwise + cumulative over time, w/ CIs

# or sweep specifications (length × #geos × significance) and recommend one:
grid = design.recommend(test_lengths=[4, 6, 8, 12], n_geos_options=[3, 5, 10, 20],
                        target_lift=0.05, alphas=[0.05, 0.10])
print(grid.summary())
grid.plot("tradeoffs.png")    # the tradeoffs figure below

Guardrails — can you trust the design? diagnose(...) visualizes the pre-period fit (treated vs synthetic control), seasonality, holdout share, and surfaces plain-language warnings when the design is risky:

Recommendations across specifications. recommend(...) sweeps test length × number of geos × significance level (alpha) and points you at the cheapest design that still detects your target lift — with a readable figure of the tradeoffs (MDE vs length per #geos, an intuitive heatmap, and alpha sensitivity):

Multi-cell tests. multi_cell(...) runs several disjoint treatment cells simultaneously — each measured against a shared donor pool that excludes every cell's treated markets, so cells never borrow each other as controls. You get a per-cell MDE/confidence/holdout report and a combined figure:

Evaluate a test that ran. evaluate(...) is the measurement counterpart to the power analysis: fit the SC family on a test that already happened, blend them into a weighted-average ensemble estimate, and report each one's lift, confidence interval (in-space placebo), and cumulative incremental — with an in-space placebo p-value:

And the effect over time — the pointwise effect across the full timeline (pre-period included, so you can see it sit flat in the noise band before the test and break out after) plus the running cumulative incremental, each as a point estimate with a confidence band:

Pin in / drop markets. select_markets/recommend take include=[…] (force must-treat markets into every candidate) and exclude=[…] (drop markets entirely — never treated, never a control). exclude is also accepted by power, diagnose, and evaluate to keep a market out of the donor pool.

Messy DataFrame? No problem. from_long coerces real-world data: outcome strings → numeric (with a clear error on genuinely non-numeric values), dates (string or unsorted) → chronological columns, locations → market names, duplicate rows aggregated with a warning, and a clear error (with a count) if the panel is gappy. You don't pre-clean dtypes.

What you get out of the box:

• Real-data power — historical placebo with injected lift on your actual panel (not an assumed variance), across SC, ASC, SDID, demeaned SC (Ferman-Pinto) and robust SC with a recommended method, plus a naive-DiD baseline for improvement-over-naive.
• MDE three ways — minimum detectable effect as a % lift, an absolute per-period change, and the cumulative incremental over the whole window, each with confidence intervals.
• Real-world guardrails — seasonality detection, a pre-period stability score, holdout share, and plain-language warnings when the design is risky (weak fit, too few donors, short history, volatile markets).
• A 0–100 confidence score and a one-line verdict, so non-experts get a clear go/no-go.
• Market selection that searches candidate treatment sets and ranks them by power, MDE, fit, holdout, and confidence.
• Multi-cell tests — several disjoint treatment cells powered at once against a shared donor pool, with a per-cell MDE/confidence report.
• A weighted-average ensemble of the base methods (combined per placebo window, with auto inverse-variance weights) for a steadier estimate than any one method — configurable via ensemble_members, and justified by the simulation study in benchmarks/sim_methods.py.
• Donor weight bounds — max_weight caps how much of the synthetic market any single donor may carry, min_weight forces diversification; both work on power(), evaluate(), select_markets() and the estimator classes.
• Post-test evaluation — evaluate() measures a test that already ran: per-method + ensemble lift, in-space placebo CIs, cumulative incremental, and a p-value.

See examples/geo_demo.py.

Inference

All bootstrap/permutation engines produce bit-identical results regardless of RAYON_NUM_THREADS, because replicate b always draws from Xoshiro256pp::substream(seed, b). Verified in CI at 1 and 8 threads.

Performance

From-scratch Rust + multithreaded inference vs the textbook NumPy + SciPy-SLSQP synthetic control. Numbers below are measured live by benchmarks/make_plots.py (median over 5 panels per size, Apple M4 Pro; absolute times vary by hardware — re-run to regenerate).

A single synthetic-control fit, across donor-pool sizes (log scale):

The per-fit win compounds under inference — a full in-space placebo test runs one fit per donor (here 200), multithreaded:

The constrained-weight estimators (SC, ASC, SDID) all beat the SLSQP reference — SDID most of all, since it solves two constrained weight problems per fit:

Estimates are identical (ATT |Δ| ≈ 1e-11; same placebo p-value) — this is pure implementation speed, not an approximation. panelkit is also far steadier: SciPy SLSQP has occasional convergence cliffs on near-collinear donor panels (one took 9.5 s in testing), which is why the table reports the median typical case.

The honest exception — MC-NNM. Matrix completion's bottleneck is the SVD, and there our hand-written one-sided Jacobi SVD is ~20× slower than NumPy's np.linalg.svd (≈112 ms vs ≈5 ms per fit at N=200). NumPy calls LAPACK — decades of hand-tuned, vendor-optimized assembly that a from-scratch SVD won't beat. We keep MC-NNM self-contained rather than fast; if you need it at scale, a LAPACK-backed SVD is the lever. See BENCHMARKS.md.

How the speed happens — Frank–Wolfe vs SLSQP, in plain English

Synthetic control picks donor weights that (a) are non-negative and (b) sum to 1 — a simplex-constrained least-squares problem. Two ways to solve it:

• SLSQP (what the SciPy reference uses) is a general-purpose constrained optimizer: it handles arbitrary nonlinear objectives and constraints by repeatedly solving quadratic sub-problems and doing line searches. Powerful and general — but it pays for that generality with overhead, and it can stall (those convergence cliffs) when the donor columns are nearly collinear.

• Frank–Wolfe (what panelkit uses) is purpose-built for exactly this shape. Each step it asks one cheap question — "which single donor most improves the fit right now?" — moves toward that donor, and (in our away-step variant) can also pull weight off a donor that's hurting. No matrix inverses, no general constraint machinery; every step stays on the simplex by construction.

  • Pros: very fast per step, naturally produces sparse interpretable weights, no convergence cliffs, and the same routine powers SC, ASC, and both of SDID's weight problems.
    • Cons: it's specialized — it only solves this convex constrained shape, not arbitrary optimization. (For the one problem that isn't this shape — MC-NNM's SVD — we don't have a tuned-assembly equivalent, hence the exception above.)

That trade — a specialized solver where the problem is regular, honest about the one place a general tuned library wins — is the whole performance story.

Monte Carlo, power analysis & robustness sweeps

Running an estimator across thousands of simulated panels is the common heavy workload (power curves, robustness checks). fit_many does the whole loop in Rust across all cores — far faster than a Python for loop calling .fit():

import numpy as np
from panelkit import SyntheticControl

# stack of R simulated panels, shape (R, N, T)
stack = np.stack([simulate_panel(tau=1.0) for _ in range(2000)])
atts = SyntheticControl().fit_many(stack, treated=[0], treat_time=45)  # (2000,) ATTs
power = np.mean(np.abs(atts) > crit)        # ← one cell of a power curve

~2000 SC fits in ~0.05 s here; a full 5-point power curve (12k fits) in ~3.6 s. See examples/power_demo.py. More efficiency levers for large runs are listed in BENCHMARKS.md.

Architecture

A four-crate Cargo workspace with a strict dependency DAG:

linalg  ←  estimators  ←  inference  ←  pypanelkit (PyO3 bindings)

Only pypanelkit touches Python; linalg depends on nothing but std.

Development

cargo test --workspace                       # Rust tests (incl. SVD cross-oracle)
cargo clippy --workspace --all-targets -- -D warnings
cargo fmt --all -- --check
maturin develop --release --manifest-path crates/pypanelkit/Cargo.toml
pytest python/tests
cargo bench -p panelkit-estimators           # criterion micro-benchmarks

License

MIT OR Apache-2.0.