Partially Pooled SCM (PPSCM)

When to Use This Estimator

PPSCM is a faithful port of augsynth::multisynth – the partially pooled synthetic control of Ben-Michael, Feller and Rothstein [PPSCM] for staggered adoption. Use it when several units are treated but at different times, with a pool of never-treated (or late-treated) comparison units, and you want a single estimate of the average treatment effect on the treated (ATT) over relative time (time-since-treatment), pooling information across cohorts.

The central idea is a pooling dial \(\nu\). Fitting a separate synthetic control for each treated unit gives the best per-unit pre-treatment fit but high variance; a fully pooled control (one synthetic match for the average treated unit) is stable but may fit any individual unit poorly. PPSCM interpolates between the two, choosing \(\nu\) to balance overall and unit-level imbalance. time_cohort=True collapses units sharing an adoption time into a single fully-pooled cohort (one synthetic control per cohort).

The problem PPSCM solves is that the two reflexive extensions of SCM to staggered adoption are each flawed. Separate SCM (fit a synthetic control per treated unit, then average – common practice) requires a good synthetic control for every treated unit, which often fails, and its strong per-unit fits can still leave the average poorly matched, biasing the ATT. Pooled SCM (match the average treated unit) nails the average fit but can fit individual units badly, biasing unit-level effects and the average when the data-generating process drifts over time. Ben-Michael, Feller and Rothstein bound the estimation error by both the average imbalance and the per-unit imbalances, and partially pooled SCM minimises a weighted combination of the two – the regime where neither extreme is trustworthy.

Reach for PPSCM when

• Several units are treated at different adoption times, with a pool of never-treated (or not-yet-treated) comparison units.

• You want an ATT over relative time (an event-study path), pooling information across cohorts instead of fitting each cohort alone.

• No single donor mix matches every treated unit, so separate SCM leaves you with unreliable per-unit fits – the partial-pooling dial lets the average fit borrow strength without abandoning unit-level fit.

• You want an estimator that nests the familiar special cases (separate and fully pooled SCM) and a principled way to choose between them.

Do not use PPSCM when

• All treated units adopt at the same time (a single cohort). The staggered machinery is unnecessary; use classic SC (Two-Step Synthetic Control, Forward Difference-in-Differences (FDID)) or, for many treated units at one time, Synthetic Difference-in-Differences (SDID).

• You are willing to assume parallel trends after weighting and want the DiD-flavoured double weighting / time weights. Synthetic Difference-in-Differences (SDID) (and, for efficiency under interactive fixed effects, Sequential Synthetic Difference-in-Differences (Sequential SDiD)) is the more natural home; PPSCM is a synthetic-control estimator, not a difference-in-differences one.

• Spillovers violate SUTVA across the donor pool – use Spatial Synthetic Difference-in-Differences (SpSyDiD).

• The treated paths lie outside the donor convex hull / the donor pool is large and noisy. Partial pooling cannot manufacture a hull that does not contain the treated units; a factor-model (Factor Model Approach (FMA)) or low-rank (Cluster Synthetic Controls (CLUSTERSC), Matrix Completion with Nuclear Norm Minimization (MCNNM)) approach is better.

• The pre-period is short relative to the donor pool. Synthetic control needs enough pre-periods to pin the weights down; with more donors than periods to balance, the simplex is wide enough that near-exact balance is available for free and says nothing about how well the donors track the treated units. This is the regime where difference in differences is the justified estimator – it asks for many units and few periods, which is the opposite trade – so reach for Synthetic Difference-in-Differences (SDID), or for PPSCM’s own method="callaway_santanna". design.underdetermined reports when a fit is in this regime.

• Distributional effects (quantiles, tails) – use Distributional Synthetic Control (DSC).

Notation

All units \(\mathcal{N} \coloneqq \{1, \ldots, N\}\) are observed over periods \(t \in \mathcal{T} \coloneqq \{1, \ldots, T\}\). A treated unit (or cohort) \(j\) adopts at period \(T_j\); never-treated units have \(T_j = \infty\) and form the donor pool \(\mathcal{N}_0 \coloneqq \{j \in \mathcal{N} : T_j = \infty\}\) of cardinality \(N_0\). The panel is split at the last adoption time, the canonical point \(T_0\), into a pre-period \(\mathcal{T}_1 \coloneqq \{t \in \mathcal{T} : t \le T_0\}\) of length \(T_0\) and a post-period \(\mathcal{T}_2 \coloneqq \{t \in \mathcal{T} : t > T_0\}\). For cohort \(j\), donor weights \(\mathbf{w}_j\) live on the simplex \(\Delta^{N_0} \coloneqq \{\mathbf{w} \in \mathbb{R}_{\ge 0}^{N_0} : \|\mathbf{w}\|_1 = 1\}\); the synthetic control matches the cohort’s pre-treatment residuals. The per-period effect is \(\tau_t\) and the average treatment effect on the treated is \(\widehat{\tau}\).

Method

PPSCM follows multisynth in three stages.

1. Two-way fixed effects (fixedeff=True, the default). A time effect is the never-treated units’ per-period mean; a unit effect is each unit’s mean over its own pre-adoption window. Both are removed and the synthetic control balances the residuals – the “intercept-shifted” estimator of the paper.

2. Partially pooled QP. With per-cohort pre-treatment imbalance \(\mathbf{q}_j \coloneqq \mathbf{x}_j - \mathbf{X}_{0,j}\mathbf{w}_j\) (residuals; the pooled imbalance aligned by relative time), the weights solve

\[\min_{\{\mathbf{w}_j \in \Delta^{N_0}\}} \; \frac{\nu}{\text{norm}_{\text{pool}}\,J^2} \Bigl\|\textstyle\sum_j \mathbf{q}_j\Bigr\|^2 + \frac{1-\nu}{\text{norm}_{\text{sep}}\,J} \sum_j \frac{\|\mathbf{q}_j\|^2}{\text{ndim}_j} + \lambda \sum_j \|\mathbf{w}_j\|^2 ,\]

where \(\text{norm}_{\text{pool}}\) and \(\text{norm}_{\text{sep}}\) are the separate-fit (nu=0) global and individual imbalance norms. Small \(\nu\) approaches a separate SCM per cohort; large \(\nu\) a fully pooled SCM.

3. Choosing \(\nu\). With nu="auto" (default) PPSCM uses augsynth’s triangle-inequality ratio \(\nu = \text{global\_l2}\cdot\sqrt{T_0}/\text{avg\_l2}\) from the separate fit; a float fixes it.

The program is posed on residuals divided by a power of two that brings them to unit magnitude, and the answer is read back in the caller’s units. Synthetic control is scale-equivariant, so this leaves the estimand alone: multiplying every series by a constant leaves the weights where they were and scales the effect with them. What it changes is what the solver is asked for. The separate fit normalizes by one, because it is the fit that produces \(\text{norm}_{\text{pool}}\) and \(\text{norm}_{\text{sep}}\), so its objective carries the square of whatever units the outcome happens to be in, while the solver’s convergence test is an absolute tolerance. On a panel running at \(10^5\) that fit takes twenty times as long as the identical problem at unit scale, and past about \(10^3\) it stops converging at all and a fallback solver answers to a looser standard. The divisor is the median absolute residual: stage 1 has already removed the level, and a panel of markets spans an order of magnitude in size, so the typical residual and not the largest sets the scale. Residuals already within \(2^3\) of unit magnitude are left alone.

Assumptions / Remarks.

Assumption 1 (no anticipation, parallel residual trends). After removing the two-way fixed effects, the treated cohorts’ residual paths would have matched a convex combination of donor residual paths absent treatment. Remark. This is the staggered-adoption analogue of the SCM identifying assumption; the fixed effects absorb level and common-time shifts so the weights only need to match the residual dynamics.

Assumption 2 (overlap / donor availability). Each cohort has eligible donors – never-treated units, or units treated more than n_leads periods later. Remark. Late-treated units can serve as “clean” controls for earlier cohorts until they themselves are treated, which the donor-eligibility rule enforces.

Remark (pooling). \(\nu\) is a bias–variance dial, not an identification parameter: the estimand (the wATET over the treated cohorts) is the same; \(\nu\) only trades per-cohort fit against stability of the pooled average.

Auxiliary covariates

By default PPSCM matches on the pre-treatment outcome path alone. Passing covariates=[...] also balances a set of auxiliary covariates, following the paper’s Section 5.2. Each covariate is z-scored against the never-treated controls and rescaled to the outcome scale, so covariate and outcome imbalance share a footing; the covariate imbalance is then stacked into both the pooled and the separate terms of the partially-pooled objective. Time-varying covariates are aggregated to their mean over the periods before the first adoption. Balancing covariates typically improves covariate balance at a small cost to the pre-treatment outcome fit – the usual bias/variance trade of matching on more.

res = PPSCM({"df": df, "outcome": "y", "treat": "d",
             "unitid": "unit", "time": "period",
             "covariates": ["income_1959", "student_teacher_ratio_1959"]}).fit()

This reproduces augsynth::multisynth’s covariate mode (y ~ d | income + ratio); see PPSCM — augsynth multisynth (Paglayan collective bargaining) and the ppscm_paglayan_covs benchmark for the cell-by-cell cross-check against a live augsynth 0.2.0 run.

Reaching Callaway-Sant’Anna and Sun-Abraham

Ben-Michael, Feller and Rothstein (2022, p.369) observe that with uniform donor weights their intercept-shifted estimator “is equivalent to recent proposals for DiD estimators that allow for treatment effect heterogeneity with a fixed donor set per treatment time cohort (see Callaway & Sant’Anna, 2020; Sun & Abraham, 2020)”. Measured, that is not an approximation: three independent implementations agree to 1e-14 once three conventions are aligned.

The three are separate settings, because each is independently meaningful:

donor_weights

"scm" (default) solves the partially-pooled QP; "uniform" puts equal weight on every admissible donor, which is the comparison Callaway-Sant’Anna and Sun-Abraham make. It is the \(\lambda \to \infty\) limit of the same program, written in closed form; the derivation is below.

base_period

"all_pre" (default) is augsynth’s: each unit’s mean over its whole pre-adoption window. "pre_treatment" is the single period \(g-1\) that Callaway-Sant’Anna normalise against. On its own the choice shifts each cohort’s level without moving the event-study shape.

donor_pool

"window" (default) admits any unit untreated through the cohort’s whole estimation window, \(g_i > g + H\). "never_treated" and "not_yet_treated" are the Callaway-Sant’Anna comparison groups. The first two coincide exactly when every other cohort adopts inside the window.

method="callaway_santanna" sets all three at once (and selects their standard error, below), leaving any convention the caller set explicitly alone:

res = PPSCM({"df": df, "outcome": "y", "treat": "d",
             "unitid": "unit", "time": "period",
             "method": "callaway_santanna"}).fit()

The three estimators diverge in exactly one regime, and it is a documented difference and not a defect: when a later cohort outlives an earlier cohort’s estimation window, augsynth admits it as a donor and Callaway-Sant’Anna do not. On four cohorts spread over a long window the gap is about 1.2e-02. A panel with adoptions spread widely lands there, so a difference of that size between donor_pool="window" and donor_pool="never_treated" is the conventions disagreeing, not a bug.

The equivalence is not a coincidence of two formulas. Callaway-Sant’Anna sits at the end of a path the partially-pooled program already contains, and the conventions above are the coordinates of that endpoint. Ben-Michael, Feller and Rothstein introduce the \(\lambda\) term (their Section 3) as

a term that penalizes the weights towards uniformity, with hyperparameter \(\lambda\). While we penalize the sum of the squared weights, there are many options, for example, an entropy or elastic net penalty

so uniformity is what the penalty is for. Following it to its limit is what produces the other estimator.

Step 1: the barycenter is what any such penalty selects. Let \(\Omega : \Delta^{N_0} \to \mathbb{R}\) be strictly convex and permutation-symmetric, so \(\Omega(\mathbf{P}\mathbf{w}) = \Omega(\mathbf{w})\) for every permutation matrix \(\mathbf{P}\). Strict convexity gives a unique minimiser \(\mathbf{w}^\star\); symmetry makes \(\mathbf{P}\mathbf{w}^\star\) a minimiser too, so \(\mathbf{P}\mathbf{w}^\star = \mathbf{w}^\star\) for all \(\mathbf{P}\), and the only permutation-invariant point of the simplex is its barycenter:

\[\operatorname*{arg\,min}_{\mathbf{w} \in \Delta^{N_0}} \Omega(\mathbf{w}) = \bar{\mathbf{w}} \coloneqq \tfrac{1}{N_0}\mathbf{1}.\]

The squared norm \(\sum_i w_i^2\) and the negative entropy \(\sum_i w_i \log w_i\) are both of this form, so the alternatives BFR list have the same endpoint. Geometrically \(\bar{\mathbf{w}}\) is the point of the simplex nearest the origin, which is why the ridge sends the weights there.

Step 2: the program converges to it, and forgets \(\nu\). Write the partially-pooled objective as \(f_\nu(\mathbf{w}) + \lambda \Omega(\mathbf{w})\), equivalently \(\lambda^{-1} f_\nu(\mathbf{w}) + \Omega(\mathbf{w})\). Since \(\Delta^{N_0}\) is compact and \(f_\nu\) continuous, \(\lambda^{-1} f_\nu \to 0\) uniformly, so \(\mathbf{w}_\lambda \to \bar{\mathbf{w}}\) for every \(\nu\). The pooling dial is inert in the limit: it weights a term that has been scaled away.

Step 3: the rate, and what \(\nu\) does instead. The barycenter has every coordinate \(1/N_0 > 0\), so it lies in the relative interior and no non-negativity constraint is active near it. For large \(\lambda\) the program is therefore smooth on the affine hull \(\{\mathbf{1}'\mathbf{w} = 1\}\), and stationarity \(\nabla f_\nu(\mathbf{w}_\lambda) + 2\lambda \mathbf{w}_\lambda + \mu\mathbf{1} = \mathbf{0}\) linearised at \(\bar{\mathbf{w}}\) gives

\[\mathbf{w}_\lambda = \bar{\mathbf{w}} - \frac{1}{2\lambda}\,\mathbf{P}\,\nabla f_\nu(\bar{\mathbf{w}}) + O(\lambda^{-2}), \qquad \mathbf{P} \coloneqq \mathbf{I} - \tfrac{1}{N_0}\mathbf{1}\mathbf{1}' ,\]

with \(\mathbf{P}\) the projection onto the simplex’s tangent space. So the approach is first order in \(\lambda^{-1}\) – not the \(\lambda^{-1/2}\) a boundary solution would give – along a fixed direction that is tangent to the simplex. That direction is where \(\nu\) survives: it sets the direction in which partially pooled SCM departs from Callaway-Sant’Anna, having no say in where the path ends.

Measured on a three-cohort panel with 41 never-treated donors, against the Callaway-Sant’Anna estimate:

\(\lambda\)	\(|\widehat{\tau} - \widehat{\tau}_{\mathrm{CS}}|\)	\(\max_i |w_i - 1/N_0|\)
0	1.0e-01	2.6e-01
1e6	7.7e-07	1.4e-07
1e12	7.7e-13	1.4e-13

A decade of \(\lambda\) buys a decade of accuracy, in the weights and in the estimate alike, down to machine precision. The scaled departure \(\lambda(\mathbf{w}_\lambda - \bar{\mathbf{w}})\) settles to a fixed vector (norms 0.410880, 0.411064, 0.411066 at \(\lambda = 10^4, 10^6, 10^8\); direction cosine 1.000000 to eight decimals), and its size moves with \(\nu\) alone – 0.2846, 0.6152, 0.9541 at \(\nu = 0, 0.5, 1\) – while the limit does not move at all.

Step 4: the endpoint is the estimator. At \(\bar{\mathbf{w}}\), with base_period="pre_treatment" subtracting \(Y_{i,g-1}\) and donor_pool="never_treated" supplying the comparison group \(\mathcal{C}\), cohort \(g\)’s horizon-\(k\) effect is

\[\widehat{\tau}_{g,g+k} = \frac{1}{n_g}\sum_{i \in \mathcal{G}_g}\bigl(Y_{i,g+k} - Y_{i,g-1}\bigr) - \frac{1}{n_{\mathcal{C}}}\sum_{i \in \mathcal{C}} \bigl(Y_{i,g+k} - Y_{i,g-1}\bigr) = \widehat{ATT}(g, g+k),\]

the two-period, two-group difference in differences Callaway and Sant’Anna identify under their Assumptions 1-4 with a never-treated comparison group. The common time effect cancels between the two group means, so removing it changes nothing. This is BFR’s own reading of their Equation (9): with uniform weights it “is the simple average over all two-period, two-group DiD estimates”, which they call “equivalent to recent proposals … (see Callaway & Sant’Anna, 2020; Sun & Abraham, 2020)”.

One difference hides in that sentence. BFR average over all pre-treatment lags, where Callaway-Sant’Anna normalise on \(g-1\) alone – which is exactly the base_period setting, and why the equivalence needs all three conventions and not just uniform weights.

Aggregation closes it. PPSCM averages cohorts by size at each horizon and then averages horizons,

\[\widehat{\theta} = \frac{1}{H}\sum_{h} \widehat{\theta}_h , \qquad \widehat{\theta}_h = \sum_{k \in \mathcal{K}_h} \frac{p_k}{S_h}\,\widehat{\tau}_{g_k, g_k + h} , \quad S_h = \sum_{k \in \mathcal{K}_h} p_k ,\]

which is Callaway-Sant’Anna’s dynamic aggregation followed by an average over event time. It coincides with their simple aggregation when every cohort is the same size and reaches every horizon, and not otherwise – so equal-sized cohorts hide the distinction instead of establishing it.

The event-study window decides which cells enter that sum. n_leads defaults to the last cohort’s post window, which is the shortest, so every cohort reaches every horizon and \(\mathcal{K}_h\) is the full set of cohorts at each \(h\). Raising it adds horizons that only the early cohorts observe: the late cohorts report NaN there and \(\mathcal{K}_h\) thins as \(h\) grows. The ceiling is the longest post window, since no cohort observes anything past the end of the panel, and a larger request is cut to it.

That longer window is the one Callaway-Sant’Anna’s simple aggregation runs over, and the per-unit paths carry it. The mean over their finite entries weights every unit-post-period cell equally, which is what did::aggte(type = "simple") reports, while res.effects.att stays the dynamic aggregation above:

import numpy as np

res = PPSCM({"df": df, "outcome": "y", "treat": "d", "unitid": "unit",
             "time": "period", "method": "callaway_santanna",
             "n_leads": 6}).fit()

dynamic = res.effects.att                 # aggte(type="dynamic")$overall.att
paths = np.vstack([u.tau for u in res.per_unit.values()])
simple = np.nanmean(paths)                # aggte(type="simple")$overall.att

Raising the window is not free with donor_pool="window", where a unit is a donor to a cohort only if it stays untreated through that cohort’s estimation window: a longer window is a stricter pool. The Callaway-Sant’Anna comparison groups, never_treated and not_yet_treated, do not depend on it.

What this does not close is the donor pool. Uniform weights and the \(g-1\) baseline are choices inside the program; \(\mathcal{C}\) is the set the program ranges over. When a later cohort outlives an earlier cohort’s estimation window the two sets genuinely differ, and no \(\lambda\) reconciles them – which is the regime described above.

In practice donor_weights="uniform" is the limit written down instead of approached: exact, with no quadratic program to solve and no \(\lambda\) for the caller to guess. The path matters because it explains why the two estimators are the same object, and it is verified in test_ppscm_cs_ridge_limit.py, which pins every number quoted above and checks the limit against diff-diff itself where it is installed.

Inference

PPSCM reports the paper’s delete-one jackknife: drop each unit, refit the full estimator (holding \(\nu\) fixed), and form \(\widehat{\text{se}}^2 = \tfrac{N-1}{N}\sum_{j \in \mathcal{N}}(\widehat{\tau}_j - \bar{\tau})^2\) for the overall ATT and each relative-time horizon, with Wald intervals. inference_method="bootstrap" swaps in augsynth’s default Mammen wild bootstrap, which reweights the single fit instead of refitting.

Dropping a control and dropping a treated unit are not the same experiment. Removing a control moves the synthetic counterfactual a little. Removing a treated unit removes one of the few draws the effect is averaged over, and that is the sampling variability of the pooled estimand. A jackknife ensemble made entirely of control deletions measures donor substitution.

That distinction bites with a single treated unit. Deleting it leaves no treated unit, so the estimator cannot be refit and the replicate is skipped: every survivor is a control deletion, and the reported standard error cannot be moved by the treated unit’s own outcomes at all. Raising a planted effect from 2.0 to 52.0 on such a panel moves the ATT by 50 and leaves the standard error bit-identical.

PPSCM therefore refuses instead of reporting that number. The refusal names how many replicates were admitted and how many removed a treated unit, so the distinction is visible in the message. Two treated units give one treated deletion, which is a thin ensemble but not a vacuous one, and it is allowed.

With one treated unit, reach for an inference method that does not rest on deleting treated units — inference_method="bootstrap" reweights a single fit — or report the point estimate without an interval.

method="callaway_santanna" reaches their point estimate exactly, and an equal point estimate does not make an equal interval. The preset therefore also selects inference_method="influence_function", which is the standard error that goes with that estimate.

Under those conventions each cell is a two-period, two-group difference in differences and its influence function is available in closed form:

\[\widehat{\phi}_i(g,t) = \frac{\Delta_i - \bar{\Delta}_{\mathcal{G}_g}}{n_g}\ \ (i \in \mathcal{G}_g), \qquad \widehat{\phi}_i(g,t) = -\frac{\Delta_i - \bar{\Delta}_{\mathcal{C}}}{n_{\mathcal{C}}}\ \ (i \in \mathcal{C}),\]

with \(\Delta_i = Y_{it} - Y_{i,g-1}\) and \(\widehat{\text{se}}(g,t) = \sqrt{\sum_i \widehat{\phi}_i(g,t)^2}\). A standard error then costs one pass over the panel, where the jackknife costs one refit per unit.

PPSCM averages cohorts by size at each horizon and then averages horizons, \(\widehat{\theta} = \tfrac1H \sum_h \widehat{\theta}_h\) with \(\widehat{\theta}_h = \sum_{k \in \mathcal{K}_h} p_k \widehat{\tau}_{kh} / S_h\) and \(S_h = \sum_{k \in \mathcal{K}_h} p_k\). That is Callaway and Sant’Anna’s dynamic aggregation followed by an average over event time, and it coincides with their simple aggregation when every cohort is the same size and reaches every horizon. The cohort shares \(p_k\) are estimated from the same panel, so the aggregate carries their sampling error through \(\partial \theta_h / \partial p_j = (\tau_{jh} - \theta_h)/S_h\) – the term R’s did::aggte calls wif. Deleting it leaves a standard error that is finite, plausible and too small, so the test suite computes the aggregate both ways and pins the difference.

results.inference_detail then carries group_time_att and group_time_se keyed by the public (adoption time, time) labels, the per-unit influence matrix every reported standard error is a functional of, and two bands on the event-time path. cband=True tabulates the second one:

\[c_{1-\alpha} = \text{quantile}_{1-\alpha}\ \max_h \Bigl| \sum_i v_i \widehat{\psi}_{h,i} \Bigr| \big/ \widehat{\text{se}}_h ,\]

with \(v_i\) Mammen (1993) multipliers. One critical value covers every horizon at once, which is the level a reader assumes when they read the path as a path (“positive by horizon three and never back”); the pointwise band read that way covers less. No refit is involved – the multipliers act on the influence functions the point estimate already produced.

The derivation assumes the conventions that produce the Callaway-Sant’Anna estimate, so inference_method="influence_function" is available only with donor_weights="uniform", base_period="pre_treatment", donor_pool="never_treated" and fixedeff=True. Solved SCM weights are estimated too and contribute a term of their own, and a not-yet-treated pool changes the comparison group’s composition over time; outside the four the formula is wrong and not approximate, and the configuration raises MlsynthConfigError naming the convention that broke it. Setting a convention explicitly alongside the preset is a coherent question whose answer is the jackknife, so the preset stands down instead of raising.

Verification: the analytical standard errors are pinned cell by cell and in aggregate against benchmarks/reference/ppscm_cs/reference.py, a transcription of diff-diff 3.9.0 at commit d9cd475 kept in-tree so the check runs without a runtime dependency. Measured agreement is 0.0 per cell and 5.6e-17 on the aggregate.

Reading the balance diagnostics

design.global_l2 and design.ind_l2 are the pre-treatment imbalance the fitted weights leave, and pct_improve_global / pct_improve_ind express them against the uniform-weight baseline. A residual near zero means a good fit only if the program could have done worse, so the design also reports the shape of the problem that produced it:

max_donors

the widest admissible donor pool across cohorts.

balance_periods

the pre-periods actually balanced, capped by the cohort with the least history, since that cohort binds.

underdetermined

whether max_donors - 1 > balance_periods. Weights on the simplex carry max_donors - 1 degrees of freedom against balance_periods equations, so past that point exact balance is generically attainable.

Both panels diff-diff ships illustrate the difference. On mpdta – 500 counties over five years, the shape difference in differences is built for – the binding cohort adopts in the second period, leaving one pre-period against a pool of 480: global_l2 comes to 3.4e-06 and 100 percent better than uniform, which describes the geometry. On castle_doctrine, 32 donors against five periods cannot reach zero, and its 2.1e-02 is a fit.

underdetermined is a statement about the program’s shape and not a verdict. The simplex is bounded, so a wide pool whose convex hull misses the treated path still leaves a large residual, and that residual is informative. What the flag rules out is reading a small one as evidence.

Per-unit fits alongside the pooled report

Because partially pooled SCM fits a separate synthetic control per treated unit (or per cohort with time_cohort=True) and averages them into the ATT, the unit-level estimates are the components of the pooled one – so both are read off a single fit. results.per_unit is a dict keyed the same as donor_weights_by_cohort; each value is a PPSCMUnitFit carrying the unit’s att, its relative-time tau path, its donor_weights, its adoption time and member units, and its in-sample fit prefit_rmspe – the root-mean-square pre-treatment imbalance \(q_j\) of that unit’s synthetic control.

Each PPSCMUnitFit additionally carries a per-unit prediction interval on its time-averaged effect – ci_lower / ci_upper with a band-implied p_value – populated when run_inference is on. It is built by the CFPT/SCPI out-of-sample interval engine (mlsynth.utils.scpi_helpers, the same machinery behind MSQRT’s bands), applied to each unit’s own pre-period residuals and post-period gap with the synthetic-control weights held fixed. This is the per-unit analogue of the pooled inference above: the delete-one jackknife (or bootstrap) quantifies uncertainty across units for the aggregate, whereas the per-unit SCPI band quantifies each unit’s own effect. A naive permutation over the QP-optimised pre-period residuals would over-reject – the fit makes those residuals small, so they are not exchangeable with the post-period gaps – which is why the per-unit band uses the SCPI construction, not a residual permutation.

The two levels reconcile exactly, so the unit-level and pooled reports never disagree: the reported separate imbalance design.ind_l2 equals \(\sqrt{\tfrac1J\sum_j q_j^2}\), and the n_units-weighted per-horizon average of the unit tau paths reproduces event_study.tau and hence the aggregate effects.att. This makes it a one-line switch to serve either request – pooled error via design.ind_l2 / global_l2 and the aggregate ATT, or per-unit estimates and their in-sample error via results.per_unit.

A caveat for whoever reads the unit-level numbers: at a high \(\nu\) (heavily pooled), the per-unit synthetic controls fit poorly, so a unit’s att is only as trustworthy as its prefit_rmspe – read the two together, and prefer a lower \(\nu\) (toward separate SCM) when unit-level estimates are the deliverable.

res = PPSCM(config).fit()
res.design.ind_l2                       # pooled/separate in-sample error
res.effects.att                         # aggregate ATT
for label, uf in res.per_unit.items():  # per-unit estimates + in-sample error
    print(label, uf.att, uf.prefit_rmspe)

The cumulative effect per unit

The per-unit bands above cover a unit’s effect in a single period, or its average across the post-period. A different question is what a unit gained in total: the sum of its effects over the periods since adoption. Setting conformal_horizon adds a band for that total to every per_unit entry:

res = PPSCM({..., "conformal_horizon": 8}).fit()
for label, uf in res.per_unit.items():
    print(label, uf.cumulative_effect,
          (uf.cumulative_lower, uf.cumulative_upper), uf.cumulative_windows)

An interval for a running total is not the running total of the per-period intervals. Adding endpoints up treats the period errors as moving in lockstep, so the width grows with the number of periods, not with its square root; rescaling a single period’s interval by the horizon assumes the opposite. Which is right depends on how the errors accumulate, and neither assumption measures it.

So the band measures it. An origin slides across the pre-period, and at each one every treated unit is treated as if it had adopted there: partially-pooled SCM fits them all in a single solve, so one pass yields each unit’s summed error over the following window at the cost of one solve per origin, not one per unit per origin. Those sums are conformity scores for exactly the quantity being reported, and the half-width is the \(\lceil (m+1)(1-\alpha) \rceil\)-th order statistic of the centred scores (mlsynth.utils.conformal.cumulative_conformal_interval(), shared with VanillaSC’s inference="conformal_cumulative"). Each fit sees only data before the window it scores, so the scores carry no in-sample optimism, and origins step by a whole horizon, so the windows do not overlap.

Two things to note. The band is additional, not a mode: inference_method still chooses the bootstrap or jackknife behind the ATT, and leaving conformal_horizon unset changes nothing. And it costs pre-period. Non-overlapping windows of length \(L\) are scarce, so a \(1-\alpha\) band needs at least \(\lceil 1/\alpha \rceil - 1\) of them: roughly \(T_0 \gtrsim L/(0.7\,\alpha)\) counting the training block held back at the start. When they run out, cumulative_lower/cumulative_upper are infinite and cumulative_windows says how many were available, instead of a narrow band that does not cover.

Where the roll starts is conformal_min_train_frac, a fraction of the pre-period: the first origin sits at \(\max(10,\ \text{frac} \times T_0)\), so every calibration fit has that many periods to train on. The default 0.3 suits most panels, and two situations call for moving it. Periods spent on the warm-up are periods not available for calibration, so lowering it buys windows when the level is out of reach; against that, a fit trained on fewer periods than there are donors can interpolate its training window, and raising the fraction past the donor count removes those origins. The two pull opposite ways, so the choice belongs to whoever knows the panel. With \(T_0 = 120\), \(L = 7\) and 60 donors, the default starts at period 36 and yields twelve windows, four of them trained on fewer periods than there are donors; 0.5 starts at 60 and removes all four, leaving eight windows, which no longer supports a 90% band. Read cumulative_windows back off each unit to see what a given choice bought.

\(T_0\)	frac	windows	at \(L = 7\)
120	0.3	12	supports 90%, not 95%
120	0.4	10	supports 90%
120	0.5	8	below the 90% threshold

Which calibration set: conformal_method

The band above calibrates on non-overlapping windows of the pre-period, and that reference set is small. With \(m \approx 0.7\,T_0/L\) windows, a finite \(1-\alpha\) band needs \(\lceil (m+1)(1-\alpha) \rceil \le m\), which puts a floor of about \(12.8\,L\) pre-periods under it. Past the floor, at whatever level a given \(m\) just supports, the required rank equals \(m\) itself, so the order statistic never trims anything and the half-width is simply the largest calibration score:

\(T_0\)	L	windows	tightest level	rank
120	8	10	90%	10 of 10
120	4	21	95%	21 of 21
104	13	5	80%	5 of 5
156	8	13	90%	13 of 13

On a thirty-market weekly panel with 120 pre-weeks, an eight-week horizon and \(\alpha = 0.05\), that leaves ten windows against a rank of eleven, and every treated unit’s band is infinite. Marketing geo-lift panels sit here routinely: the horizon is a campaign flight and the pre-period is however much history the advertiser has.

conformal_method="cyclic" calibrates instead against the \(T\) cyclic shifts of the residual path, a reference set whose size is the length of the panel and not a count of windows. Neither the floor nor the rank-never-trims regime applies to it. On the same panel, 118 of 120 unit-fits return a finite band and 114 of those 118 cover, against 0 of 120 finite for the split band:

res = PPSCM({..., "conformal_horizon": 8,
                  "conformal_method": "cyclic"}).fit()
for label, uf in res.per_unit.items():
    print(label, uf.cumulative_effect,
          (uf.cumulative_lower, uf.cumulative_upper), uf.cumulative_p_value)

The price is a shape assumption. Inverting a test needs a null to subtract, and the null here is a constant per-period effect: a candidate \(\theta\) is subtracted from the treated unit’s post-period, the panel is refitted, and the adjusted residual path is compared against its own cyclic shifts by a moving-block statistic. The reported band is \(L\) times the range of the candidates the test accepts. An effect that ramps is outside that null family, and the honest outcome is an empty accepted set, reported as nan bounds – distinct from None, which means no band was asked for at all. The split band assumes only that the calibration windows are exchangeable with the post-period window, and reports the accumulated effect directly, so it makes no claim about the effect’s shape.

conformal_method="resample" attacks the same shortage from a third direction, and it is the cheapest of the three. The split band runs a rolling pass over the pre-period and reduces each window to its total, so \(m\) windows give \(m\) numbers. The same pass computes an \(L\)-period path on the way to each of those totals, and the resample band keeps them: its reference set is the \(m \times L\) per-period errors rather than the \(m\) totals. Each draw assembles a post-period path from circular blocks of a unit’s own errors, flipping each block’s sign with probability one half, and the band is the \(1-\alpha\) quantile of the absolute accumulated draw:

res = PPSCM({..., "conformal_horizon": 8,
                  "conformal_method": "resample"}).fit()

Because the reference set counts periods, the window floor does not apply, and a panel with seven windows against a required rank of eight – infinite under the split band – returns a finite one. It buys that without the cyclic band’s shape assumption: nothing is subtracted, no null family is posited, and an effect that ramps is as admissible as a flat one. It also refits nothing beyond the rolling pass the split band already pays for, where the cyclic band pays a refit per candidate in its grid.

What it does assume is that conformal_block is long enough to carry the serial correlation of the period errors. The variance of an \(H\)-period total is \(H\gamma_0 + 2\sum_k (H-k)\gamma_k\), so drawing periods independently – conformal_block=1, Wheeler’s original construction – keeps only the first term and reports a band too narrow whenever the errors are positively autocorrelated. The default, 0, means the whole horizon, the longest block the accumulated total is sensitive to. A block longer than the horizon is clamped to it.

The three report the same estimand, so conformal_method selects between them the way inference_method selects the bootstrap or the jackknife behind the ATT, and there is one set of bounds whichever is chosen. Their diagnostics differ and so occupy separate fields: cumulative_windows counts calibration windows and is filled by the split and resample bands, cumulative_p_value is the cyclic band’s permutation p-value of the no-effect null, and the unused one is None. cumulative_method says which produced the bounds. A window count means different things to the two bands that report it – the split band’s order statistic is taken over exactly those windows, while the resample band draws from the \(L\) periods inside each – which is what cumulative_method is there to disambiguate.

Their parameters differ too, and a parameter belonging to a method not chosen is refused by name at config time: conformal_n_nulls and conformal_grid_scale are the cyclic band’s, conformal_block and conformal_n_sim are the resample band’s, and conformal_min_train_frac is shared by the split and resample bands, which run the same pass and differ only in how they read it. The grid is an approximation and it errs in one direction: the band is the range of accepted candidates, so a coarse grid samples fewer of them and reports a band too narrow, converging upward as conformal_n_nulls rises – measured on a two-unit panel, a width of 3.14 at five candidates against 5.86 at thirty-one. Coarse is anti-conservative here, the opposite of the usual intuition about discretisation. An accepted set reaching an end of the grid is bounded by conformal_grid_scale and not by the data, and that end is reported as infinite.

Every candidate is a refit, so the cyclic band costs about fourteen times the split one – 9.6s against 0.7s per fit at the geometry above. Reach for it when the split band comes back infinite, or when the pre-period is too short for the floor.

The cumulative effect overall

The band above is per unit. The corresponding question about the pool is what the treated units gained in total over the first \(L\) periods, and cumulative_band=True answers it:

res = PPSCM({..., "cumulative_band": True}).fit()
band = res.inference_detail.cumulative
for L, point, lo, hi in zip(band.horizons, band.point, band.lower, band.upper):
    print(L, point, (lo, hi))

Both the jackknife and the wild bootstrap already fit the estimator many times and get a whole per-horizon path back from each fit. Those paths are kept on res.inference_detail.replicate_paths, and the band is built from them, so it costs no refits beyond the inference that was going to run anyway.

Keeping them is what makes the band possible. Collapsing each replicate to one standard error per horizon – which is all a per-period band needs – discards how the horizons move together, and that covariance is the entire content of a cumulative interval. A caller with only the collapsed standard errors has to choose an assumption instead: adding period interval endpoints treats the errors as moving in lockstep and grows the width like \(L\), while rescaling a single period’s interval assumes they are independent and grows it like \(\sqrt{L}\). Accumulating the replicates before taking the standard error measures which is true.

The band is simultaneous. A cumulative path is read as a path – “the total is positive by week six and stays there” is a claim about every horizon at once – and a pointwise band read that way covers at well below its nominal level, by more as the number of horizons grows. One shared critical value (mlsynth.utils.supt.supt_critical_value(), the sup-t construction of Montiel Olea and Plagborg-Moller) restores the level for the whole path.

That multiplier is the \(1 - \alpha\) quantile of \(\max_h |z_h| / s_h\), where \(z\) is the vector of horizon errors and \(s_h\) the standard error each is divided by. The construction assumes \(s_h\) is the true standard error, in which case the ratio is normal and its maximum has a known law. Here \(s_h\) is estimated from the same replicates that supply \(z\), so the ratio is a Student-t and not a \(z\), and its maximum is wider – by more as the ensemble shrinks. Reading the multiplier off the normal law would therefore hand back a band narrower than \(1 - \alpha\) at every finite ensemble size, and by a margin that grows exactly where the ensemble is small enough to make it matter.

supt_critical_value simulates the whole statistic instead. It draws \(z \sim N(0, R)\) and an independent \(Q \sim \text{Wishart}(n - 1, R)\) for the estimated variances, which is exact for Gaussian replicates because a sample mean and a sample covariance are independent, and takes the quantile of \(\max_h |z_h| / \sqrt{Q_{hh} / (n - 1)}\). At one horizon this reduces to a Student-t on \(n - 1\) degrees of freedom, which is what a single-horizon interval should have been all along. Passing reference="normal" restores the older behaviour for a caller reproducing a number from an earlier release.

How much it changes depends on the ensemble, and the two ensembles PPSCM offers sit at opposite ends. A wild bootstrap draws as many replicates as it is asked for – at 999 and thirteen horizons the correction is a tenth of a percent, so the two references agree to the digit. A delete-one jackknife has one replicate per treated unit, and a design with a handful of treated units is where the correction is the whole story: at \(n = 50\) it widens the multiplier by 4 percent, at \(n = 20\) by 13. Simulated against a known zero over 3000 draws at thirteen horizons, a nominal 95 percent band read off the normal law covers at 0.932 with fifty replicates and 0.910 with twenty; the studentized reference returns 0.949 and 0.957 on the same draws.

Which ensemble produced the band is recorded on band.method, because the two are not interchangeable. The wild bootstrap reweights each unit’s residual by an independent multiplier, which does not cancel the common factors the synthetic weights cancel in the point estimate, so its replicate variance is inflated where factor structure is strong. The delete-one jackknife refits the weights on each leave-one-out, so the factors re-cancel per replicate. The jackknife replicates also carry the delete-one inflation and the bootstrap draws do not, since the latter are already on the estimator’s scale; the band applies whichever matches.

Empirical Illustration: mandatory collective bargaining

The multisynth vignette studies the effect of state mandatory collective-bargaining laws on log per-pupil education expenditure (Paglayan 2018), a staggered design. basedata/Teachingaugsynth.scv ships the panel; the analysis restricts to 1959–1997, drops DC and Wisconsin, and treats a state from the year it required bargaining.

import numpy as np
import pandas as pd
from mlsynth import PPSCM

url = "https://raw.githubusercontent.com/jgreathouse9/mlsynth/refs/heads/main/basedata/Teachingaugsynth.scv"
df = pd.read_csv(url)
df = df[~df["State"].isin(["DC", "WI"])]
df = df[(df["year"] >= 1959) & (df["year"] <= 1997)].copy()
df["cbr"] = (df["year"] >= df["YearCBrequired"].fillna(np.inf)).astype(int)

res = PPSCM({"df": df, "outcome": "lnppexpend", "treat": "cbr",
             "unitid": "State", "time": "year", "display_graphs": True}).fit()

print(f"nu (auto)   : {res.design.nu_used:.4f}")
print(f"Average ATT : {res.att:.3f}  (SE {res.inference.se:.3f})")

This prints:

nu (auto)   : 0.2607
Average ATT : -0.011  (SE 0.020)

reproducing the augsynth vignette (nu = 0.2607, Average ATT -0.011). Setting time_cohort=True collapses to adoption-time cohorts and gives nu = 0.3939, Average ATT -0.017 (augsynth: -0.018).

Verification

Note

Exact replication of augsynth. On the Paglayan data PPSCM matches augsynth::multisynth to high precision: the auto-\(\nu\) agrees to four decimals (0.2607 default, 0.3939 time-cohort), the Average ATT matches (\(-0.011\) default; \(-0.017\) vs \(-0.018\) time-cohort), and the raw global/individual L2 imbalances agree (0.003 / 0.028). The full relative-time event study matches the vignette’s per-horizon averages to 3–4 decimals. The decisive fidelity detail is aligning the pooled imbalance by relative time on top of two-way fixed effects. The jackknife SE (0.020) is close to augsynth’s default wild-bootstrap SE (0.022); they differ only by inference procedure. This is locked in by test_matches_augsynth_vignette in mlsynth/tests/test_ppscm.py.

A cross-package case on real data covers both modes at once: the cannabis-alcohol panel of Ronczewski (2026), which runs augsynth and did side by side on one sample. PPSCM’s default reproduces the published multisynth ATT to 8.2e-08 and its callaway_santanna mode reproduces the published did::aggte simple aggregate to 4.9e-17, with all six dynamic event-study coefficients to 2.0e-16. Both need n_leads = 6, which the paper sets against a default of 2. See benchmarks/cases/ronczewski_cannabis.py.

Core API

Partially Pooled Synthetic Control (PPSCM) estimator.

A thin orchestration over mlsynth.utils.ppscm_helpers, faithfully porting augsynth::multisynth:

Ben-Michael, E., Feller, A., & Rothstein, J. (2022). “Synthetic Controls with Staggered Adoption.” JRSS-B 84(2):351-381.

PPSCM removes two-way fixed effects, balances the residuals with a partially-pooled QP (nu interpolating between separate and fully pooled SCM), and reports a relative-time event study and overall ATT with the paper’s delete-one jackknife. time_cohort=True collapses units sharing an adoption time into one fully-pooled cohort.

class mlsynth.estimators.ppscm.PPSCM(config: PPSCMConfig | dict)

Bases: object

Partially Pooled SCM estimator (augsynth::multisynth port).

Parameters:

config (PPSCMConfig or dict) – Validated configuration. Reads nu (pooling, or "auto"), fixedeff, n_leads, n_lags, time_cohort, lam, run_inference and alpha beyond the common panel fields.

Returns:

PPSCMResults – Design (pooling level + balance diagnostics), relative-time event study, overall ATT with jackknife inference, donor weights, and per_unit – the per-treated-unit (or per-cohort) fits (att, in-sample prefit_rmspe, tau path, donor weights) that are the components of the pooled estimate and reconstruct design.ind_l2 and the aggregate ATT.

fit() → PPSCMResults

Fit PPSCM and return the typed result container.

Configuration

class mlsynth.config_models.PPSCMConfig(*, df: ~pandas.DataFrame, outcome: str, treat: str, unitid: str, time: str, display_graphs: bool = True, save: bool | str = False, counterfactual_color: ~typing.List[str] = <factory>, treated_color: str = 'black', plot: ~mlsynth.config_models.PlotConfig = <factory>, nu: float | ~typing.Literal['auto'] = 'auto', fixedeff: bool = True, n_leads: ~typing.Annotated[int | None, ~annotated_types.Ge(ge=1)] = None, n_lags: ~typing.Annotated[int | None, ~annotated_types.Ge(ge=1)] = None, time_cohort: bool = False, donor_weights: ~typing.Literal['scm', 'uniform'] = 'scm', base_period: ~typing.Literal['all_pre', 'pre_treatment'] = 'all_pre', donor_pool: ~typing.Literal['window', 'never_treated', 'not_yet_treated'] = 'window', method: ~typing.Literal['callaway_santanna'] | None = None, lam: ~typing.Annotated[float, ~annotated_types.Ge(ge=0)] = 0.0, solver: ~typing.Any = None, run_inference: bool = True, inference_method: str = 'jackknife', cband: bool = False, n_boot: ~typing.Annotated[int, ~annotated_types.Ge(ge=1)] = 1000, seed: int = 0, alpha: ~typing.Annotated[float, ~annotated_types.Gt(gt=0.0), ~annotated_types.Lt(lt=1.0)] = 0.05, conformal_horizon: int | None = None, conformal_min_train_frac: float = 0.3, conformal_method: ~typing.Literal['split', 'cyclic', 'resample'] = 'split', conformal_block: int = 0, conformal_n_sim: int = 2000, conformal_n_nulls: int = 25, conformal_grid_scale: float = 3.0, covariates: ~typing.List[str] | None = None, cumulative_band: bool = False)

Configuration for the Partially Pooled SCM (PPSCM) estimator.

Implements Ben-Michael, Feller & Rothstein (2022, JRSS-B 84(2):351-381). Targets staggered-adoption designs by minimizing a weighted average of the per-treated-unit imbalance q_sep and the average-treated imbalance q_pool, with weighting hyper- parameter nu.

alpha: float

base_period: Literal['all_pre', 'pre_treatment']

cband: bool

conformal_block: int

conformal_grid_scale: float

conformal_horizon: int | None

conformal_method: Literal['split', 'cyclic', 'resample']

conformal_min_train_frac: float

conformal_n_nulls: int

conformal_n_sim: int

covariates: List[str] | None

cumulative_band: bool

donor_pool: Literal['window', 'never_treated', 'not_yet_treated']

donor_weights: Literal['scm', 'uniform']

fixedeff: bool

inference_method: str

lam: float

method: Literal['callaway_santanna'] | None

model_config: ClassVar[ConfigDict] = {'arbitrary_types_allowed': True, 'extra': 'forbid'}

Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

n_boot: int

n_lags: int | None

n_leads: int | None

nu: float | Literal['auto']

run_inference: bool

seed: int

solver: Any

time_cohort: bool

Result Containers

PPSCM.fit() returns a PPSCMResults: the PPSCMDesign (pooling level and balance diagnostics), the relative-time PPSCMEventStudy, the overall PPSCMInference, and the per-cohort donor weights.

Typed, NumPy-first result containers for Partially Pooled SCM (staggered).

PPSCM ports augsynth::multisynth (Ben-Michael, Feller & Rothstein 2022): a partially-pooled synthetic control for staggered adoption that interpolates, via nu, between a separate SCM per treated unit (nu small) and a fully pooled SCM (nu large), on top of two-way fixed effects.

class mlsynth.utils.ppscm_helpers.structures.PPSCMCumulativeBand(horizons: ndarray, point: ndarray, lower: ndarray, upper: ndarray, se: ndarray, critical_value: float, alpha: float, n_replicates: int, method: str)

Bases: object

Simultaneous band for the cumulative (running-total) effect path.

lower[L] and upper[L] bound the total effect over horizons 0..L, and the band covers every L at once with probability 1 - alpha – so a statement about the path (“positive by week six and never back”) is covered at the stated level, which a pointwise band read the same way is not.

se is the standard error of the running total at each horizon, taken from the replicate paths themselves. That is what makes the band’s growth an observation instead of an assumption: independent period errors grow it like sqrt(L), perfectly correlated ones like L, and the replicates carry whichever is true.

alpha: float

critical_value: float

horizons: ndarray

lower: ndarray

method: str

n_replicates: int

point: ndarray

se: ndarray

upper: ndarray

class mlsynth.utils.ppscm_helpers.structures.PPSCMDesign(nu_used: float, lam: float, fixedeff: bool, time_cohort: bool, n_leads: int, n_lags: int, global_l2: float, ind_l2: float, scaled_global_l2: float, scaled_ind_l2: float, conventions: Dict[str, Any] | None = None, max_donors: int | None = None, balance_periods: int | None = None)

Bases: object

The fitted design: pooling level and balance diagnostics.

balance_periods: int | None = None

conventions: Dict[str, Any] | None = None

The three conventions this fit ran under – donor weighting, unit-effect baseline and donor eligibility – plus which inference produced the interval. Together they decide whether the fit is augsynth’s multisynth or the Callaway-Sant’Anna estimator, which the numbers alone do not say.

fixedeff: bool

global_l2: float

ind_l2: float

lam: float

max_donors: int | None = None

The shape of the balance problem the imbalances above came out of. max_donors is the widest admissible donor pool across cohorts and balance_periods the pre-periods actually balanced, capped by the cohort with the least history since that cohort binds. Read them with global_l2: a residual near zero means a good fit when the program was constrained and describes the geometry when it was not.

n_lags: int

n_leads: int

nu_used: float

property pct_improve_global: float

property pct_improve_ind: float

scaled_global_l2: float

scaled_ind_l2: float

time_cohort: bool

property underdetermined: bool | None

Whether the balance system had more freedom than constraints.

Weights on the simplex carry max_donors - 1 degrees of freedom against balance_periods equations. When the first exceeds the second, exact balance is generically attainable and a near-zero global_l2 says nothing about how well the donors track the treated units. It is a statement about the program’s shape and not a verdict on the fit: the simplex is bounded, so a wide pool whose hull misses the treated path still leaves a large residual, and that residual is informative.

class mlsynth.utils.ppscm_helpers.structures.PPSCMEventStudy(horizons: ndarray, tau: ndarray, se: ndarray, ci: ndarray)

Bases: object

Relative-time (time-since-treatment) average ATT path.

ci: ndarray

horizons: ndarray

se: ndarray

tau: ndarray

class mlsynth.utils.ppscm_helpers.structures.PPSCMInference(att: float, se: float, ci: Tuple[float, float], method: str, replicate_paths: ndarray | None = None, cumulative: PPSCMCumulativeBand | None = None, group_time_att: Dict[Tuple[Any, Any], float] | None = None, group_time_se: Dict[Tuple[Any, Any], float] | None = None, pointwise_band: ndarray | None = None, uniform_band: ndarray | None = None, critical_value: float | None = None, influence: ndarray | None = None)

Bases: object

Overall (post-period average) ATT and its inference.

att: float

ci: Tuple[float, float]

critical_value: float | None = None

cumulative: PPSCMCumulativeBand | None = None

group_time_att: Dict[Tuple[Any, Any], float] | None = None

group_time_se: Dict[Tuple[Any, Any], float] | None = None

influence: ndarray | None = None

method: str

pointwise_band: ndarray | None = None

replicate_paths: ndarray | None = None

se: float

uniform_band: ndarray | None = None

class mlsynth.utils.ppscm_helpers.structures.PPSCMInputs(Xy: ndarray, trt: ndarray, n_pre: int, time_labels: ndarray, units: ndarray, outcome: str, intervention_time: Any, Z: ndarray | None = None, cov_names: tuple | None = None)

Bases: object

Preprocessed staggered panel (the only pandas touchpoint is setup).

Parameters:

• Xy (np.ndarray) – Full outcome matrix, shape (n, T) (units x all periods).

• trt (np.ndarray) – Adoption index per unit (position in time_labels); inf for never-treated controls.

• n_pre (int) – Number of pre-treatment periods (columns before the last adoption).

• time_labels (np.ndarray) – Sorted time labels, length T.

• units (np.ndarray) – Unit labels, length n.

• outcome (str) – Outcome column name.

• intervention_time (Any) – The last adoption time (pre/post split point).

• Z (np.ndarray or None) – Per-unit auxiliary-covariate matrix, shape (n, d_cov), aggregated to the pre-first-adoption mean. None when no covariates are given.

• cov_names (tuple of str or None) – Names of the covariate columns, length d_cov.

Xy: ndarray

Z: ndarray | None = None

property control_units: ndarray

cov_names: tuple | None = None

intervention_time: Any

property n: int

n_pre: int

outcome: str

time_labels: ndarray

property treated_units: ndarray

trt: ndarray

units: ndarray

class mlsynth.utils.ppscm_helpers.structures.PPSCMResults(*, effects: ~mlsynth.config_models.EffectsResults | None = None, fit_diagnostics: ~mlsynth.config_models.FitDiagnosticsResults | None = None, time_series: ~mlsynth.config_models.TimeSeriesResults | None = None, weights: ~mlsynth.config_models.WeightsResults | None = None, inference: ~mlsynth.config_models.InferenceResults | None = None, method_details: ~mlsynth.config_models.MethodDetailsResults | None = None, sub_method_results: ~typing.Dict[str, ~typing.Any] | None = None, additional_outputs: ~typing.Dict[str, ~typing.Any] | None = None, raw_results: ~typing.Dict[str, ~typing.Any] | None = None, execution_summary: ~typing.Dict[str, ~typing.Any] | None = None, plot_config: ~mlsynth.config_models.PlotConfig | None = None, inputs: ~mlsynth.utils.ppscm_helpers.structures.PPSCMInputs, design: ~mlsynth.utils.ppscm_helpers.structures.PPSCMDesign, event_study: ~mlsynth.utils.ppscm_helpers.structures.PPSCMEventStudy, inference_detail: ~mlsynth.utils.ppscm_helpers.structures.PPSCMInference, donor_weights_by_cohort: ~typing.Dict[~typing.Any, ~typing.Dict[~typing.Any, float]], per_unit: ~typing.Dict[~typing.Any, ~mlsynth.utils.ppscm_helpers.structures.PPSCMUnitFit] = <factory>, metadata: ~typing.Dict[str, ~typing.Any] = <factory>)

Bases: BaseEstimatorResults

Top-level container returned by mlsynth.PPSCM.fit().

An EffectResult. PPSCM is a staggered / partially-pooled estimator, so the standardized time_series carries the pooled event-time effect path (gap = horizon effect, counterfactual = no-effect baseline), and effects.att is the aggregate ATT – mirroring the SequentialSDID convention. The native objects are preserved: inference_detail (the PPSCMInference, formerly inference) and donor_weights_by_cohort (the nested per-cohort weights, formerly donor_weights); the contract names inference / donor_weights are taken by the base contract.

design: PPSCMDesign

donor_weights_by_cohort: Dict[Any, Dict[Any, float]]

event_study: PPSCMEventStudy

inference_detail: PPSCMInference

inputs: PPSCMInputs

metadata: Dict[str, Any]

model_config: ClassVar[ConfigDict] = {'arbitrary_types_allowed': True, 'extra': 'forbid', 'frozen': True, 'json_encoders': {<class 'numpy.ndarray'>: <function BaseEstimatorResults.Config.<lambda>>}}

Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

property nu: float

per_unit: Dict[Any, PPSCMUnitFit]

class mlsynth.utils.ppscm_helpers.structures.PPSCMUnitFit(label: str, adoption_time: Any, member_units: List[str], n_units: int, att: float, prefit_rmspe: float, tau: ndarray, pre_imbalance: ndarray, donor_weights: Dict[Any, float], ci_lower: float | None = None, ci_upper: float | None = None, p_value: float | None = None, tau_lower: ndarray | None = None, tau_upper: ndarray | None = None, cumulative_effect: float | None = None, cumulative_lower: float | None = None, cumulative_upper: float | None = None, cumulative_windows: int | None = None, cumulative_p_value: float | None = None, cumulative_method: str | None = None)

Bases: object

The synthetic-control fit for one treated unit (or cohort) in the pool.

Partially-pooled SCM fits a separate synthetic control per treated unit (or per adoption cohort with time_cohort=True) and averages them into the ATT, so these are the components of the pooled estimate at the chosen nu – not a separate re-run. The two aggregates reconstruct exactly: the reported design.ind_l2 equals sqrt(mean_j prefit_rmspe_j**2), and the n1-weighted per-horizon average of the tau paths reproduces the pooled event study.

label

The unit label (time_cohort=False) or adoption-time label (time_cohort=True); matches the key in donor_weights_by_cohort.

Type:

str

adoption_time

The (public) time label at which this unit / cohort adopts treatment.

Type:

Any

member_units

Treated unit label(s) in this group (one unless a cohort pools several).

Type:

list of str

n_units

Cohort size (len(member_units)); the aggregation weight n1.

Type:

int

att

This unit’s/cohort’s average post-treatment effect (mean of tau).

Type:

float

prefit_rmspe

Pre-treatment in-sample fit error q_j – the root-mean-square pre-period imbalance of this synthetic control (residual, fixed-effect- removed space, matching the estimator’s balance objective). A large value flags a poorly fit unit whose att should not be over-trusted (the nu-pooling caveat). Aggregates to design.ind_l2.

Type:

float

tau

Relative-time effect path (length n_leads); NaN past this unit’s observed horizon.

Type:

np.ndarray

tau_lower, tau_upper

Per-period (pointwise) CFPT/SCPI prediction band for the effect path, aligned with tau (length n_leads, NaN past the horizon); the per-horizon counterpart of ci_lower / ci_upper and wider than them (the time average shrinks by sqrt(L), a single period does not). None when inference is off.

Type:

np.ndarray or None

pre_imbalance

The pre-treatment imbalance vector (front-padded to the balance window) whose weighted RMS is prefit_rmspe; the per-period in-sample residual.

Type:

np.ndarray

donor_weights

{donor_label: weight} for this unit’s synthetic control (nonneg, sums to 1).

Type:

dict

adoption_time: Any

att: float

ci_lower: float | None = None

ci_upper: float | None = None

cumulative_effect: float | None = None

cumulative_lower: float | None = None

cumulative_method: str | None = None

cumulative_p_value: float | None = None

cumulative_upper: float | None = None

cumulative_windows: int | None = None

donor_weights: Dict[Any, float]

label: str

member_units: List[str]

n_units: int

p_value: float | None = None

pre_imbalance: ndarray

prefit_rmspe: float

tau: ndarray

tau_lower: ndarray | None = None

tau_upper: ndarray | None = None

Helper Modules

Staggered long-to-wide formatting (the only DataFrame touchpoint): derive adoption times, split pre/post at the last adoption.

Long-DataFrame -> NumPy boundary for PPSCM (staggered adoption).

Mirrors augsynth::format_data_stag: derive each unit’s first treated period, split the panel at the last adoption time into pre (X) and post (y), and index adoption by position in the sorted time vector (Inf for never-treated).

mlsynth.utils.ppscm_helpers.setup.prepare_ppscm_inputs(df: DataFrame, *, outcome: str, treat: str, unitid: str, time: str, covariates: List[str] | None = None) → PPSCMInputs

The engine: two-way fixed effects (fit_feff), the partially-pooled QP, auto-\(\nu\), and the relative-time event study / ATT.

Core staggered-adoption engine for PPSCM, ported faithfully from augsynth::multisynth (Ben-Michael, Feller & Rothstein 2022).

Pipeline (one call = one fit):

1. fit_feff removes fixed effects (force=3 two-way: time effect from never-treated column means + per-cohort unit pre-mean) and balances the residuals.

2. solve_cohort_qp solves the partially-pooled QP over donor weights, with the pooled imbalance aligned by relative time (front-padded) and the pooled/separate terms normalized by the separate fit’s norms.

3. run_multisynth chooses nu (triangle-inequality ratio when “auto”), refits, and produces the relative-time event study and ATT.

Validated to reproduce the multisynth vignette exactly (default nu=0.2607, ATT=-0.011; time_cohort nu=0.3939, ATT=-0.017).

class mlsynth.utils.ppscm_helpers.engine.Conventions(donor_weights: str = 'scm', base_period: str = 'all_pre', donor_pool: str = 'window')

The three choices that decide which estimator a fit is.

Carried as one object because they travel together and are forgotten separately: every place that refits the panel – the jackknife, the conformal calibration’s rolling origins – has to run the estimator the caller configured, and three keyword arguments with defaults let a refit site keep an older signature and silently pin itself to augsynth’s. A fourth convention added here reaches every refit without anyone editing them.

The defaults are augsynth’s multisynth, so a fit that says nothing is the port it has always been.

base_period: str = 'all_pre'

donor_pool: str = 'window'

donor_weights: str = 'scm'

mlsynth.utils.ppscm_helpers.engine.balance_shape(adopt_of, donors, groups, n_lags: int)

The balance problem’s shape: widest donor pool, and periods balanced.

balance_periods is capped by the cohort with the least pre-history, because that cohort’s block is what the padded design can actually constrain – asking for ten lags of a cohort adopting at t=4 balances three. See PPSCMDesign.underdetermined.

mlsynth.utils.ppscm_helpers.engine.eligible_donors(trt: ndarray, adopt: int, n_leads: int, donor_pool: str) → ndarray

Indices admissible as donors for a cohort adopting at adopt.

Never-treated units carry a non-finite adoption time, so they satisfy every rule below. The three differ only in which later-adopting units they also admit, and that choice is what separates this estimator from Callaway-Sant’Anna when a later cohort outlives an earlier one’s window:

• "window" – untreated through this cohort’s whole estimation window, trt > adopt + n_leads. augsynth’s rule, and the default.

• "never_treated" – never treated at all, which is what CS and Sun-Abraham use by default.

• "not_yet_treated" – untreated as of this cohort’s adoption.

"window" and "never_treated" coincide exactly when every other cohort adopts inside the window, which is why the three estimators agree to machine precision there and diverge otherwise.

mlsynth.utils.ppscm_helpers.engine.fit_feff(Xy: ndarray, trt: ndarray, adopt_indices, fixedeff: bool, base_period: str = 'all_pre') → Dict[int, ndarray]

Residualize Xy per cohort.

Returns {adoption_index: residual_matrix (n, T)}. With fixedeff the time effect is the never-treated column mean and the unit effect is a pre-adoption baseline; without it, only the time effect (control averages) is removed.

base_period selects the baseline, and the two are different estimators:

• "all_pre" – each unit’s mean over its whole pre-adoption window [:tj]. This is augsynth’s fit_feff (rowMeans(residuals[, 1:tj])) and the default, verified against it to 0.0 over 300 random panels.

• "pre_treatment" – the single period tj - 1, which is the base period Callaway-Sant’Anna and Sun-Abraham normalise against. Choosing it is one of the three conventions that make this estimator theirs (#465); on its own it shifts each cohort’s level without moving the event-study shape.

mlsynth.utils.ppscm_helpers.engine.predict_tau(res, groups, adopt_of, members, donors, W, n1, H, n, bs_weight=None)

Relative-time tau per cohort, plus the n1-weighted event study and ATT.

With bs_weight (per-unit multipliers, default all ones) this is augsynth’s predict.multisynth(bs_weight=...) written in residual space: the fixed-effect terms cancel between the treated mean and the synthetic, leaving the treated residuals (scaled by bs_weight, averaged over the cohort) minus the donor residuals weighted by W[g] * bs_weight. bs_weight = ones reproduces the point estimate exactly.

mlsynth.utils.ppscm_helpers.engine.run_multisynth(Xy: ndarray, trt: ndarray, d: int, n_leads: int, n_lags: int, *, fixedeff: bool = True, time_cohort: bool = False, nu: float | None = None, lam: float = 0.0, solver: Any = None, Z: ndarray | None = None, conventions: Conventions = Conventions(donor_weights='scm', base_period='all_pre', donor_pool='window')) → Dict[str, Any]

Run one multisynth fit; returns weights, event study, ATT, diagnostics.

mlsynth.utils.ppscm_helpers.engine.solve_cohort_qp(res, groups, adopt_of, members, donors, n1, d, n, n_lags, nu, norm_pool, norm_sep, lam, solver, zt=None, Zc=None) → Dict[Any, ndarray]

Partially-pooled QP: per-cohort simplex weights (summing to cohort size).

When zt/Zc (per-cohort scaled auxiliary-covariate target sums and donor blocks) are supplied, the covariate imbalance is stacked into the pooled and separate terms (normalized by the number of covariates), following augsynth::multisynth Sec 5.2.

mlsynth.utils.ppscm_helpers.engine.solve_scale(res: Dict[Any, ndarray]) → float

Power-of-two divisor bringing the residuals to unit magnitude.

The summary is the median absolute residual, not the panel’s median absolute level: the level is removed by fit_feff before anything reaches the program, so a panel of markets with a large mean and small variation would be divided by a factor its residuals never had, and the objective would land under eps_abs instead of on it. The median is the right summary of what is left because a panel of markets spans an order of magnitude in size and one very large market should not set the scale.

Returns 1.0 when the residuals are already within 2 ** 3 of unit magnitude, so the arithmetic of a fit that never had this problem is untouched, and when there is nothing finite and non-zero to measure.

mlsynth.utils.ppscm_helpers.engine.uniform_weights(donors, groups, n) → Dict[Any, ndarray]

Equal weight on every admissible donor – the CS/Sun-Abraham comparison.

The lam -> infinity limit of the partially-pooled program, in closed form. The barycenter minimises sum_i w_i^2 over the simplex, so the solved weights approach it as the ridge grows – at O(1/lam), for every nu, since the barycenter is interior and no non-negativity constraint binds there. Setting this reaches the limit exactly and skips the program; the QP would otherwise need a lam the caller has to guess, and would still return the answer to within solver tolerance instead of exactly.

The paper’s delete-one jackknife inference.

Delete-one jackknife inference for PPSCM (Ben-Michael et al. 2022).

The paper’s jackknife drops each unit i (treated or control), refits the full staggered estimator on the remaining n - 1 units (holding nu fixed), and forms

se^2 = (n - 1) / n * sum_i (theta_i - mean_i theta_i)^2

separately for the overall ATT and each relative-time horizon. Wald intervals are built from these SEs around the full-sample point estimates.

Every function here that refits the panel – the jackknife, and the conformal band’s rolling origins – takes the fit’s Conventions and passes it on. A replicate has to be a refit of the estimator that produced the point estimate, and #467 is what happens when it is not: the replicates ran augsynth’s donor weighting, baseline and donor pool while the estimate ran the caller’s, giving a standard error that was finite, plausible, and for a different estimator.

Two guards follow from the same incident. run_multisynth refuses a non-finite nu, since a uniform-weight fit poses no program and reports nu_used as NaN; and a jackknife that ends with fewer than two usable replicates raises instead of returning NaN, because a missing standard error read as a degenerate panel for a whole review.

mlsynth.utils.ppscm_helpers.inference.bootstrap_inference(fit: dict, *, alpha: float, n_boot: int, seed: int, per_time_full: ndarray, att_full: float, return_paths: bool = False)

augsynth’s default Mammen wild/multiplier bootstrap (weighted_bootstrap_multi).

Reweights the single fit by per-unit multipliers Z (no refit): for each draw, predict_tau(bs_weight=Z) - (sum(Z)/n_treated) * point_estimate; the bootstrap SE is the root-mean-square of the centered draws. Returns (att, se, ci, per_time_se, per_time_ci) matching jackknife_inference.

mlsynth.utils.ppscm_helpers.inference.cumulative_conformal_per_unit(Xy: ndarray, trt: ndarray, d: int, n_leads: int, n_lags: int, *, fixedeff: bool, time_cohort: bool, nu_used: float, lam: float, solver: Any, alpha: float, horizon: int, min_train_frac: float = 0.3, conventions: Conventions = Conventions(donor_weights='scm', base_period='all_pre', donor_pool='window')) → Tuple[ndarray, ndarray, ndarray, ndarray]

Per-unit conformal band for each treated unit’s cumulative effect.

The point estimates come from the full fit’s tau_rel; the calibration comes from rolling_pooled_block_sums() and mlsynth.utils.conformal.cumulative_conformal_interval() – the same combiner VanillaSC uses, so the order statistic has one definition.

Returns:

tuple of numpy.ndarray – (point, lower, upper, n_scores), each of shape (J,) in groups order. A unit with too few calibration windows for the requested level gets an infinite band rather than a narrow one that does not cover.

mlsynth.utils.ppscm_helpers.inference.cumulative_supt_band(per_time_full: ndarray, replicate_paths: ndarray, *, alpha: float, jackknife: bool = True, n_sims: int = 200000, seed: int | None = 0, method: str = 'jackknife')

Simultaneous band for the running total, from the replicate paths.

An interval for a cumulative effect is not the running total of the per-period intervals. Adding endpoints treats the period errors as moving in lockstep, so the width grows with the number of periods; rescaling one period’s interval assumes the opposite. Here the replicate paths are accumulated first and the standard error taken after, so whatever correlation the errors have is the correlation the band inherits.

The band is simultaneous over horizons (mlsynth.utils.supt.supt_critical_value()), because a cumulative path is read as a path.

Parameters:

• per_time_full (np.ndarray, shape (H,)) – The per-horizon effect path from the full fit.

• replicate_paths (np.ndarray, shape (n_replicates, H)) – One per-horizon path per replicate. Rows containing NaN are dropped, so a leave-one-out refit that failed is absent instead of counted as zero.

• alpha (float) – The band is simultaneous at 1 - alpha.

• jackknife (bool, default True) – Apply the delete-one inflation to the standard error. True for leave-one-out replicates, which differ from the full estimate by O(1/m); False for bootstrap draws, already on the estimator’s scale.

• n_sims, seed – Tabulation of the sup-t critical value.

• method (str) – Which ensemble produced the paths, recorded on the result – a jackknife band and a bootstrap band are not interchangeable numbers.

Returns:

PPSCMCumulativeBand

mlsynth.utils.ppscm_helpers.inference.cwz_cumulative_per_unit(Xy: ndarray, trt: ndarray, d: int, n_leads: int, n_lags: int, *, fixedeff: bool, time_cohort: bool, nu_used: float, lam: float, solver: Any, alpha: float, horizon: int, n_nulls: int = 25, grid_scale: float = 3.0, conventions: Conventions = Conventions(donor_weights='scm', base_period='all_pre', donor_pool='window')) → Tuple[ndarray, ndarray, ndarray, ndarray]

Per-unit cumulative band by inverting a moving-block conformal test.

The counterpart of cumulative_conformal_per_unit(), calibrated against the T cyclic shifts of the residual path instead of a disjoint split of the pre-period. The split version’s window count is roughly 0.7 * T0 / L and a finite 1 - alpha band needs ceil((m+1)(1-alpha)) <= m, so at the 90 percent level it needs nine windows – a floor of about 12.8 * L pre-periods before a band exists at all, and every feasible design then sits in the regime where the rank never trims and the half-width is simply the largest score. The cyclic reference set does not depend on the horizon, so neither the floor nor that regime applies.

The price is a shape assumption. Test inversion needs a null to subtract, so the null here is a constant per-period effect and the reported band is horizon times the accepted range of that effect. An effect that ramps is not in the null family, and the honest outcome then is an empty accepted set, which confidence_set_bounds() returns as (nan, nan).

Two construction details decide the answer, both established by measurement; the full account is in moving_block_pvalue().

The statistic is mean_abs, the reference implementation’s. The absolute block sum is the intuitive choice for a running total and is invalid here: the quadratic program leaves end-of-window residuals sign-coherent, so block sums ramp toward the end of the window while magnitudes stay flat, and the trailing block always occupies the most inflated position. It is not a compromise – displacing a mean-zero block by delta raises the mean of its absolute values, so the test has power against precisely the constant shift being inverted.

The null refit balances every period. Under the null the adjusted series is an untreated series, so all of it is fitting data, and leaving a period out would put the trailing block partly outside the fit its reference blocks come from. Balancing the whole window requires an explicit nu: the automatic rule sits exactly on its boundary when one unit is treated, and cvxpy refuses the program. That is why nu_used is a parameter, not a default.

The grid is an approximation and it errs in one direction. The band is the range of accepted candidates, so a coarse grid samples fewer of them and reports a band that is too narrow, converging upward as n_nulls rises – measured on a two-unit panel, a width of 3.14 at five candidates against 5.86 at thirty-one. Coarse is anti-conservative here, not conservative, which is the opposite of the usual intuition about discretisation.

conventions reaches both refits, the observed fit and each null one, so the band is calibrated on the estimator the caller configured. The default is augsynth’s, so a call that says nothing gets what it always got.

An accepted set that reaches an end of the grid is bounded by grid_scale and not by the data, and that end is returned as infinite. Reporting the endpoint instead would understate the band silently, since it looks like any other number.

Returns:

tuple of numpy.ndarray – (point, lower, upper, p_zero), each of shape (J,) in groups order. lower and upper are on the cumulative scale and may be infinite; p_zero is the permutation p-value of the no-effect null.

mlsynth.utils.ppscm_helpers.inference.jackknife_inference(Xy: ndarray, trt: ndarray, d: int, n_leads: int, n_lags: int, *, fixedeff: bool, time_cohort: bool, nu_used: float, lam: float, solver: Any, alpha: float, per_time_full: ndarray, att_full: float, conventions: Conventions = Conventions(donor_weights='scm', base_period='all_pre', donor_pool='window'), return_paths: bool = False) → Tuple[float, float, Tuple[float, float], ndarray, ndarray]

Return (att, se, ci, per_time_se, per_time_ci).

With return_paths the leave-one-out per-horizon paths are appended. They are computed either way; keeping them lets a caller build a cumulative band without refitting the whole jackknife a second time.

mlsynth.utils.ppscm_helpers.inference.per_unit_intervals(M: ndarray, tau_rel: ndarray, *, alpha: float, time_dependence: str = 'iid') → Tuple[ndarray, ndarray, ndarray, ndarray, ndarray]

Per-unit CFPT/SCPI prediction intervals for each unit’s effect path.

The pooled bootstrap / jackknife measures variability across units and so cannot give one treated unit its own interval. This builds a per-unit band from that unit’s own fit and reuses mlsynth’s out-of-sample interval engine (the same CFPT/SCPI machinery MSQRT uses), so PPSCM’s per-unit bands are methodologically consistent with MSQRT’s.

For unit k the bands come from its post-period effect path tau_rel[k, :] (the CFPT effects) and its pre-period residuals M[:, k] (the CFPT pre_residuals): the residual moments set the sub-Gaussian scale of the counterfactual prediction error, which correctly accounts for the in-sample fit – unlike a naive permutation over the QP-optimised pre-residuals, which are not exchangeable with the post gaps and over-reject. The engine is called per unit (one column at a time), so units with different post horizons (ragged NaN) are handled by trimming.

A single engine call returns the full CFPT family, so both the time-averaged band (TAUS) and the per-period pointwise bands (TSUS) come out of the same computation: the pointwise bands are the TAUS band’s per-horizon counterpart and are wider (TAUS shrinks by sqrt(L) under time_dependence="iid"; a single period does not).

Parameters:

• M (numpy.ndarray) – Pre-period residual columns, shape (d, J) (a 1-D array is a single unit). NaN entries are dropped per unit.

• tau_rel (numpy.ndarray) – Post-period relative-time effect paths, shape (J, H) (a 1-D array is a single unit). NaN (past a unit’s horizon) is dropped per unit.

• alpha (float) – Total miscoverage level; the interval is 100 * (1 - alpha) percent. Keyword-only.

• time_dependence ({“iid”, “general”}, default “iid”) – Time-averaging bound passed through to the CFPT engine (it affects only the time-averaged band, never the per-period bands). Keyword-only.

Returns:

tuple of numpy.ndarray – (ci_lower, ci_upper, p_value, tau_lower, tau_upper). The first three have shape (J,): the per-unit band bounds on the time-averaged ATT and a band-implied two-sided p-value (the house convention 2 * (alpha/2) ** ((point/half_width) ** 2), clamped to [0, 1]). The last two have shape (J, H) – the per-unit, per-period band bounds, aligned with tau_rel (NaN where tau_rel is NaN). A unit with no usable residuals yields NaN throughout.

mlsynth.utils.ppscm_helpers.inference.resample_cumulative_per_unit(Xy: ndarray, trt: ndarray, d: int, n_leads: int, n_lags: int, *, fixedeff: bool, time_cohort: bool, nu_used: float, lam: float, solver: Any, alpha: float, horizon: int, min_train_frac: float = 0.3, block: int = 0, n_sim: int = 2000, seed: int = 0, conventions: Conventions = Conventions(donor_weights='scm', base_period='all_pre', donor_pool='window')) → Tuple[ndarray, ndarray, ndarray, ndarray]

Per-unit cumulative band from block-resampled rolling-origin errors.

The third calibration set behind PPSCM’s per-unit cumulative band, alongside cumulative_conformal_per_unit() (a disjoint split of the pre-period) and cwz_cumulative_per_unit() (the cyclic shifts of the residual path). All three report the same estimand – the total a treated unit gained over horizon periods – so conformal_method chooses between them.

It runs the rolling pass the split band already pays for and reads the m * horizon per-period errors inside those windows instead of the m totals. Each draw assembles a post-period path from circular blocks of a unit’s own errors, flipping each block’s sign with probability one half, and the band is the 1 - alpha quantile of the absolute accumulated draw.

Two things follow, and they are the reason to reach for it. The split band needs ceil((m+1)(1-alpha)) <= m before its order statistic exists, a floor of roughly 12.8 * horizon pre-periods; drawing from periods rather than totals has no such floor, so a panel that leaves the split band infinite still gets a finite one. And unlike the cyclic band it assumes no shape for the effect and refits nothing beyond the pass itself, where cyclic pays a refit per candidate in its grid.

What it does assume is that a unit’s pre-period errors are exchangeable with its post-period ones, which is the same assumption the split band makes, and that block is long enough to carry their serial correlation. Drawing periods independently (block = 1) is Wheeler’s original construction and understates the spread of the total whenever the period errors are positively autocorrelated – the variance of an H-period total is H * gamma_0 + 2 * sum_k (H - k) * gamma_k, and independent draws keep only the first term.

Parameters:

• block (int, optional) – Block length in periods. 0 (default) means the whole horizon, the longest block the accumulated total is sensitive to. Clamped to the horizon; see mlsynth.utils.conformal.resolve_block().

• n_sim (int, optional) – Paths drawn per unit (default 2000). These cost no refits – the pass is the cost – so precision here is cheap.

• seed (int, optional) – Base seed. Each unit draws from seed + k so a unit’s band does not depend on how many other units the panel happens to carry.

Returns:

tuple of numpy.ndarray – (point, lower, upper, n_windows), each of shape (J,) in groups order. n_windows counts the calibration windows the pass realised, the same number cumulative_conformal_per_unit() reports; the draw itself uses the n_windows * horizon periods inside them. A unit whose pass realised no window gets an infinite band, not a narrow one that does not cover.

mlsynth.utils.ppscm_helpers.inference.rolling_pooled_block_sums(Xy: ndarray, trt: ndarray, d: int, n_leads: int, n_lags: int, *, fixedeff: bool, time_cohort: bool, nu_used: float, lam: float, solver: Any, horizon: int, min_train_frac: float = 0.3, conventions: Conventions = Conventions(donor_weights='scm', base_period='all_pre', donor_pool='window')) → List[ndarray]

Each calibration window reduced to its total – the split band’s scores.

One number per window, which is the estimand’s own scale, so the order statistic in mlsynth.utils.conformal.cumulative_conformal_interval() is taken over quantities directly comparable with the reported cumulative effect.

Returns:

list of numpy.ndarray – One array of m finite window totals per treated cohort, in groups order.

mlsynth.utils.ppscm_helpers.inference.rolling_pooled_period_errors(Xy: ndarray, trt: ndarray, d: int, n_leads: int, n_lags: int, *, fixedeff: bool, time_cohort: bool, nu_used: float, lam: float, solver: Any, horizon: int, min_train_frac: float = 0.3, conventions: Conventions = Conventions(donor_weights='scm', base_period='all_pre', donor_pool='window')) → List[ndarray]

The same windows left unreduced – the resample band’s calibration series.

rolling_pooled_block_sums() throws away the cross-period structure inside each window on its way to the total. A band for a running total needs it: the spread of an H-period sum depends on how the period errors move together, and the totals alone cannot say. So this returns the paths, and the caller block-resamples them.

The reference set is the m * horizon periods instead of the m totals, which is why this construction stays finite on a pre-period that leaves the split band’s order statistic undefined.

Returns:

list of numpy.ndarray – One (m, horizon) array per treated cohort, in groups order. Row i is the effect path over window i; arr.sum(axis=1) is exactly what rolling_pooled_block_sums() returns for the same panel.