Synthetic Controls for Experimental Design (MAREX)

Contents

Synthetic Controls for Experimental Design (MAREX)#

When to Use This Estimator#

The estimators elsewhere in mlsynth are retrospective: a treatment has already happened and you reweight donors to reconstruct the treated unit’s counterfactual. MAREX, due to Abadie and Zhao (2026) [ABADIE2024], is prospective — it designs an experiment. Before any treatment is assigned, and using only pre-experimental data, it chooses which aggregate units to treat and which to hold out as controls, so that the experiment you are about to run yields a credible estimate.

The motivating setting is a firm (say a ride-sharing company) that wants to test a new policy but can only deploy it in a few whole markets. A within-market A/B test is contaminated by interference (treated and control drivers compete); randomizing whole markets to treatment is unbiased ex ante but, with a handful of large units, routinely produces treated and control groups with very different baselines, so any single realisation is badly off. MAREX instead picks the treated and control markets so their pre-experiment predictors match the population — a non-randomized design that, the paper shows, substantially reduces estimation bias relative to randomization.

Reach for MAREX when:

  • Units are large aggregates (markets, regions, stores) and only one or a few can be treated.

  • You control the assignment and want to choose it well, not estimate after the fact.

  • Interference or equity rules out within-unit randomization, forcing whole-unit treatment.

Notation#

There are \(N\) units \(\mathcal{N} \coloneqq \{1, \dots, N\}\) and \(T\) periods \(t \in \mathcal{T} \coloneqq \{1, \dots, T\}\), 1-indexed; the experiment takes effect after period \(T_0\), splitting \(\mathcal{T}\) into the pre-experiment window \(\mathcal{T}_1 \coloneqq \{t \in \mathcal{T} : t \le T_0\}\) (of length \(T_0\)) and the experimental window \(\mathcal{T}_2 \coloneqq \{t \in \mathcal{T} : t > T_0\}\). Because MAREX designs an experiment instead of reweighting around one already-treated unit, units are indexed generically by \(i, j \in \mathcal{N}\) with no forced treated unit. Each unit has a pre-intervention predictor vector \(\mathbf{x}_j\) (pre-period outcomes and optional covariates); \(\bar{\mathbf{x}} = \sum_j f_j \mathbf{x}_j\) is the population predictor mean for known weights \(f_j\) (e.g. market shares, or \(1/N\)). The experimenter chooses treated weights \(\mathbf{w}\) and control weights \(\mathbf{v}\), both on the simplex, and disjoint:

\[\sum_j w_j = 1,\quad \sum_j v_j = 1,\quad w_j, v_j \ge 0,\quad w_j v_j = 0 \;\;\forall j.\]

Units with \(w_j > 0\) are treated; among the rest, units with \(v_j > 0\) form the synthetic control. Writing \(y_{jt}\) for the observed outcome (treated units realise \(y_{jt}^I\) over \(\mathcal{T}_2\), everyone else \(y_{jt}^N\)), the design estimator of the average effect is

\[\widehat{\tau}_t(\mathbf{w}, \mathbf{v}) = \sum_j w_j y_{jt} - \sum_j v_j y_{jt}, \qquad t \in \mathcal{T}_2.\]

Assumptions#

Assumption 1 (linear factor model). Potential outcomes follow

\[y_{jt}^N = \delta_t + \boldsymbol{\theta}_t^\top \mathbf{z}_j + \boldsymbol{\lambda}_t^\top \boldsymbol{\mu}_j + \varepsilon_{jt}, \qquad y_{jt}^I = \upsilon_t + \boldsymbol{\gamma}_t^\top \mathbf{z}_j + \boldsymbol{\eta}_t^\top \boldsymbol{\mu}_j + \xi_{jt},\]

with observed covariates \(\mathbf{z}_j\), unobserved factors \(\boldsymbol{\mu}_j\), and mean-zero idiosyncratic noise.

Remark. This is the interactive-fixed-effects model of the SC literature (Abadie-Diamond-Hainmueller 2010), extended with a separate factor structure for the treated potential outcome — necessary because a design must choose a treatment group, not just a comparison group.

Assumption 2 (regularity). The factor loadings are non-degenerate (\(F \le T_E\), smallest eigenvalue bounded below) and the noise is i.i.d. sub-Gaussian with common variance, independent across the two potential outcomes; dependence across units is allowed.

Remark. These conditions keep the factor structure recoverable from the pre-experiment window and the noise well-behaved, so the population predictor mean \(\bar{\mathbf{x}}\) is a meaningful matching target, not an artifact of a degenerate loading matrix.

Assumption 3 / 4 (fit quality). A weight vector reproducing the population predictor means exists exactly (Assumption 3), or approximately within a tolerance \(d\) (Assumption 4). This is the design-time analogue of “the treated unit lies in the convex hull of the donors.”

Remark. Under these conditions Abadie & Zhao bound the bias of \(\widehat{\tau}_t(\mathbf{w}, \mathbf{v})\) and develop the permutation test below; the better the pre-experiment match, the smaller the bias.

Assumption 5 (exchangeable periods). The factors \(\boldsymbol{\lambda}_t\) are exchangeable across the blank and post periods, and the noise is continuously distributed.

Remark. This is the condition of Abadie and Zhao’s Theorem 2, under which the test is exact: \(\alpha - 1/|\Pi| \le \Pr(\hat p \le \alpha) \le \alpha\). The potential outcomes may still trend through \(\delta_t\) and \(\boldsymbol{\theta}_t\). For fixed factors their Theorem OA.1 makes the p-value approximately valid as the fitting window grows. Persistent or trending factors are where it is least reliable; a unit root in the factors breaks it whatever the window lengths (#744).

Mathematical Formulation#

The Design Optimization#

MAREX chooses \(\mathbf{w}, \mathbf{v}\) (and a binary selection mask \(\mathbf{z}\), with \(w_j \le z_j\), \(v_j \le 1 - z_j\), so a unit is treated or a control, never both) to match the population predictor mean. The base design minimises

\[\min_{\mathbf{w}, \mathbf{v}, \mathbf{z}}\; \Bigl\| \bar{\mathbf{x}} - \sum_j w_j \mathbf{x}_j \Bigr\|_2^2 + \Bigl\| \bar{\mathbf{x}} - \sum_j v_j \mathbf{x}_j \Bigr\|_2^2 \quad \text{s.t. the simplex / disjointness / cardinality constraints,}\]

with the number of treated units pinned by m_eq (exactly) or bounded by m_min/m_max. This is a mixed-integer quadratic program (the binary z); mlsynth solves it with SCIP by default, or — via relaxed=True — relaxes z to \([0, 1]\), solves the QP, and discretises post hoc.

A treated unit is one the synthetic treated unit puts weight on. Abadie and Zhao bound \(\|\mathbf{w}\|_0\), the number of nonzero treated weights, and the selection \(\mathbf{z}\) alone does not: \(w_j \le z_j\) keeps an unselected market’s weight at zero but lets a selected one carry none, and that market is then neither treated nor a control. With \(\|\mathbf{w}\|_0 = m\) read literally the program can have no minimum: where the best fit of a treated set leaves one of its markets at zero, giving that market weight \(\varepsilon\) costs less as \(\varepsilon\) falls and never reaches the zero-weight design. MAREX therefore counts a treated set as a design of size \(m\) only when its own fit gives every one of its \(m\) markets a weight above \(10^{-8}\), the threshold its results use to list treated markets. The exact backend refuses any other set during its search, and the MIQP cuts such a set and solves again. Under m_min/m_max this changes no optimum, since the same set without its unweighted market fits the treated side as well and leaves one more control. Under m_eq it can: on small simulated panels about one design in eight had a market at zero weight, and a binding budget on the treated markets’ total cost made it common.

From SCIP’s solution MAREX keeps only the selection. SCIP counts a binary within its integrality tolerance of an integer as integral, so a selection variable left near \(10^{-6}\) rounds to “not treated” while still allowing a treated weight of that size, and a market can then carry weight on both sides of the design. On one five-market test panel SCIP’s point did exactly that: it used a market as a treated unit with weight \(10^{-6}\) and as a control with weight 0.146, and reached an objective below that of any feasible design. MAREX therefore fixes the selection at its rounding and refits the weights with the selection held constant, which is a convex QP, so the weights and objective it reports belong to the design the selection defines.

Most of the MIQP’s cost is SCIP proving a design optimal, which grows steeply with the number of markets and treated units. Two options manage that. A warm_start (a list of treated unit labels) seeds the search with a known good design — for instance LEXSCM’s top candidate, which solves a near-identical problem by lexicographic search, not by proof: MAREX(..., warm_start=lexscm_warm_start(lex.fit())) (the helper lexscm_warm_start lives in mlsynth.utils.marex_helpers.warmstart). The seed is only a hint, so a solve that runs to completion returns the identical proven optimum; it just starts the branch-and-bound from a strong incumbent instead of hunting for one. Paired with time_limit (seconds), MAREX then returns the best design found within the budget — typically the seed, confirmed or refined — instead of paying the full optimality proof. Both default to off and apply to the exact MIQP only. For every design, backend="exact" skips the MIQP and finds the same optimum by searching treated sets directly (Choosing the treated markets: overview).

mlsynth exposes four objective variants through design (clear names that map to the paper’s formulations):

  • "standard" — match each predictor mean with both synthetic units (formulation 5);

  • "weakly_targeted" — match the treated synthetic to the mean and softly tie the control synthetic to it (weight beta);

  • "penalized" — standard plus a distance penalty that down-weights units far from the population mean (lambda1 / lambda2);

  • "unit_penalized" — standard plus unit-level penalties (xi / lambda1_unit).

The xi term is formulation (10): it charges each candidate treated unit for how well the remaining controls reproduce that unit, weighted by how much of the synthetic treated it carries. The charge for a given unit is a minimum over the control weights fitted to it alone, so it does not depend on the design being chosen; mlsynth solves those per-unit problems once and the penalty is then linear in the treated weights. Writing it as the treated weight times the joint design’s own fit to that unit gives a different objective, and one that is not convex in the two weight vectors together.

The paper’s control-side counterpart lambda2_unit is accepted only as 0. It pairs the treated weights against the control weights through a distance matrix, which is bilinear in two decision variables; lambda1_unit charges the treated side against the population target, and the penalized design charges both sides against it separately.

Control penalty#

The weakly_targeted design takes one more term, control_penalty (\(\lambda\)), from Vives-i-Bastida (2022, program 1):

\[\Bigl\| \bar{\mathbf{x}} - \mathbf{u} \Bigr\|_2^2 + \beta \Bigl\| \mathbf{u} - \sum_j v_j \mathbf{x}_j \Bigr\|_2^2 + \beta \lambda \sum_j v_j \bigl\| \mathbf{u} - \mathbf{x}_j \bigr\|_2^2, \qquad \mathbf{u} = \sum_j w_j \mathbf{x}_j.\]

The last term is Abadie and L’Hour’s (2021) penalty on the controls: each control market is charged its distance from the synthetic treated unit, so the synthetic control is built from a few nearby markets instead of many distant ones that average out. Vives-i-Bastida’s design for the Catalan UBI pilot uses \(\lambda = 0.1\). At control_penalty=0, the default, the design is the weakly-targeted design as before.

As written, the term multiplies a control weight by a function of the treated weights. Because the control weights sum to one it equals \(\|\mathbf{u}\|^2 - 2\,\mathbf{u}^\top \mathbf{s} + \sum_j v_j \|\mathbf{x}_j\|^2\) with \(\mathbf{s} = \sum_j v_j \mathbf{x}_j\), a quadratic in the two synthetic units plus a linear term in \(\mathbf{v}\). The whole objective is then convex in \((\mathbf{w}, \mathbf{v})\) exactly when \(\beta \lambda (1 + \lambda) \le 1\), and MAREX refuses a configuration past that bound when it is built, since neither backend could certify a minimum there. At \(\beta = 1\) the bound is \(\lambda \le 0.618\); as \(\beta \to 0\), the lexicographic limit of Vives-i-Bastida’s Algorithm 1, it places no limit on \(\lambda\). The MIQP receives the objective as a sum of squares in the centred synthetic units; the exact backend handles it as described under Choosing the treated markets: overview.

Covariates#

By default the design matches on pre-period outcomes. Passing covariates (column names) appends them to the predictor vector, \(\mathbf{x}_j = [\mathbf{y}^E_j ; \mathbf{z}_j]\), exactly as in the paper — the synthetic treated and control are then balanced on both pre-period outcomes and covariates (with an optional covariate_weight scale). With standardize=True every predictor column is first divided by its standard deviation across units, and the covariate columns are then multiplied by covariate_weight, so a weight of 0.1 makes each covariate count one tenth as much as one period of the outcome. This is the order of Vives-i-Bastida (2022, Algorithm 1), whose design puts the covariates at 1/10 of the outcomes’ weight. covariate_weight also takes a mapping from covariate name to weight, so each covariate can carry its own (Vives-i-Bastida’s diagonal \(\mathbf{H}\)): for example covariate_weight={"age": 0.1, "income": 0.3}. The mapping must name exactly the listed covariates, every weight must be non-negative, and a weight of 0 leaves that covariate out of the match. A column that does not vary across units is left as it is: the weights on each side sum to one, so such a column adds nothing to the objective. When the pre-period is long, the outcomes already encode the covariates’ contribution, so covariates matter most when few pre-periods are available.

Each covariate enters as its per-unit mean over the fitting periods, the first T0 - blank_periods pre-periods, the same window the outcomes are matched on. A covariate that does not change over time is its value. A covariate that does change is averaged over the fitting window only: the blank periods are kept out of the choice of the design so that their gaps can serve as the placebo distribution for inference, and their covariate values are kept out with them. Vives-i-Bastida (2022) averages the covariates over the fitting years in the same way. The balance diagnostics below use the same means.

Clustering, Costs, and Budgets#

Passing a cluster column solves the design within each cluster (one or a few treated units per cluster), which better approximates the population predictor distribution and limits interpolation bias (paper OA.1). Per-unit costs and a budget (scalar or per-cluster) add the constraint \(\sum_j c_j w_j \le B\) on the treated weights. Because the weights sum to one, that bounds the treated group’s weighted average cost per market, not the total the test spends. A scalar budget is split evenly across clusters. Both backends take it; the exact backend’s handling is below. costs is listed in MAREX’s unit order, which sorts the unit labels, so with labels such as u0 to u11 the order is u0, u1, u10, u11, u2, ....

A budget on what the test spends is a rule on which markets are treated. market_cost_col names a column holding each market’s cost of being in the test (the media held back, the fees to run it), and max_total_cost caps the sum over the treated markets,

\[\sum_j c_j z_j \le B,\]

with \(z_j = 1\) when market \(j\) is treated. Six markets that cost one million each spend six million, whatever their weights. Both backends take it. The MIQP adds the row; the exact backend drops every treated set over budget before fitting it, the way it drops a set that breaks a border conflict, so its bounds and its proof of optimality stand as they are.

With several clusters the budget ties them together. Each cluster is searched on its own, and the best affordable design can take a cluster’s fiftieth-best treated set because it is cheap. The exact backend first drops, in each cluster, every set that would overrun the budget even if every other cluster took its cheapest admissible markets. It then keeps each cluster’s best designs in order, with their costs, and picks one per cluster to minimise the summed objective within the budget, a multiple-choice knapsack solved by a branch and bound over the clusters. A combination that takes a design past a cluster’s list scores at least that list’s last value plus every other cluster’s best. When that floor does not clear the answer, the lists it could still beat are doubled and the step repeats, so the design is the optimum over every admissible combination, the same one the MIQP returns. The lists start at 16 designs and the step refuses past 4,096, naming backend='mip'. When no admissible combination fits the budget, the fit raises MlsynthEstimationError.

res = MAREX(MAREXConfig(df=df, outcome="orders", unitid="market", time="week",
                        T0=150, m_eq=6, backend="exact",
                        market_cost_col="test_cost", max_total_cost=2_000_000)).fit()

The tests are in mlsynth/tests/test_marex_spend_budget.py, against an oracle that enumerates every treated set within budget.

Geographic design restrictions#

MAREX already carries most geographic structure natively. Its cluster field is the paper’s “region as a distinct experimental design” (p.13: weather patterns dependent across cities in the same region motivate treating each region separately); m_min / m_max per cluster is a stratum quota (“at least / at most this many treated per region”); costs / budget is the cost bound; and because each cluster’s control synthetic \(v_{\cdot,k}\) is built only from cluster-k members, donors are automatically same-region. The same restriction vocabulary SYNDES and LEXSCM expose then adds the rest, as constraints on the MIP:

  • to_be_treated / not_to_be_treated – force a market in (one leadership already committed to) or out (a regulated market, or one already under another test); a forbidden market stays a donor.

  • adjacency + spillover_threshold – no two treated markets may share a border (\(\sum_k z_{ik} + \sum_k z_{jk} \le 1\)), for when one test market’s campaign would bleed into a neighbouring one and bias the lift.

  • size_col + min_size / max_size – a treated-unit size band (the floor is a power minimum, the ceiling is synthesizability); out-of-band markets stay donors.

  • cluster_col – at most one treated market per cluster value (\(\sum_k z_{ik} + \sum_k z_{jk} \le 1\) for same-cluster pairs), a no-two-from-one-cluster rule for markets that share a media buy or retail footprint. This is distinct from the cluster design grouping above: it adds conflict constraints to the treated set instead of splitting the objective.

  • stratum_col + min_per_stratum / max_per_stratum – a coverage quota on the treated set of one design (at least / at most this many treated per stratum). Use this when the design is a single cluster but the read still has to span regions; when each region is its own design, cluster with per-cluster m_min / m_max expresses the same quota.

  • exclude_bordering_donors – drop a treated market’s bordering neighbours from its (within-cluster) control pool, so a market partly treated by spillover never sits in its counterfactual. Requires adjacency.

These are supported in the exact MIQP and in the search backend (Choosing the treated markets: overview), not with relaxed=True, whose post-hoc rounding cannot guarantee them. An over-constrained design – e.g. forbidding every market of a region, leaving its required treated synthetic with no member – raises a translated MlsynthEstimationError naming the restrictions. The example below runs on the bundled real US DMA contiguity map and Census regions, with a region-grouped linear factor sales model (after Liao, Shi & Zheng [RelaxSC]); the geography is real, the sales reproducible.

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

base = ("https://raw.githubusercontent.com/jgreathouse9/mlsynth/"
        "refs/heads/main/basedata/markets")
adj = pd.read_csv(f"{base}/dma_adjacency.csv", index_col=0)
meta = pd.read_csv(f"{base}/dma_metadata.csv")
CDC = {**{s: "Midwest" for s in ["OH", "IN", "MI", "KY"]},
       **{s: "South" for s in ["WV", "VA", "TN", "GA", "FL", "SC"]}}
meta = meta[meta["state"].isin(CDC)].copy()
meta["region"] = meta["state"].map(CDC)
names = [n for n in meta["dma_name"] if n in adj.index]
meta = meta[meta["dma_name"].isin(names)].reset_index(drop=True)
adj = adj.loc[names, names]

rng = np.random.default_rng(7)
reg = meta["region"].to_numpy(); n, r, T, T0 = len(names), 3, 30, 24
Lam = (np.array([{g: rng.normal(size=r) for g in sorted(set(reg))}[g] for g in reg])
       + 0.15 * rng.normal(size=(n, r)))
F = np.cumsum(rng.normal(size=(T, r)), axis=0)
pop = np.round(rng.lognormal(12.5, 0.7, n)).astype(int)
Y = 100 + rng.normal(0, 5, n) + F @ Lam.T + rng.normal(0, 1, (T, n))
df = pd.DataFrame([{"dma": names[j], "week": t, "sales": float(Y[t, j]),
                    "region": reg[j], "population": int(pop[j])}
                   for j in range(n) for t in range(T)])

B = dict(df=df, outcome="sales", unitid="dma", time="week",
         cluster="region", T0=T0)                    # cluster = Census region

MAREX({**B, "m_max": 2}).fit()                       # native quota: <=2 per region
MAREX({**B, "m_max": 2, "to_be_treated": ["Atlanta, GA"],
       "not_to_be_treated": ["Miami-Fort Lauderdale, FL"]}).fit()
MAREX({**B, "m_max": 2, "adjacency": adj,
       "spillover_threshold": 0.5}).fit()            # no two treated border
MAREX({**B, "m_max": 2, "size_col": "population",
       "min_size": 200_000, "max_size": 3_000_000}).fit()
MAREX({**B, "m_max": 1, "adjacency": adj, "spillover_threshold": 0.5,
       "exclude_bordering_donors": True}).fit()      # non-bordering donors

Across groups: cluster vs. quotas vs. restrictions#

Three mechanisms touch on grouping/geography; they are distinct, not interchangeable, and cluster is not made obsolete by the restrictions:

  • A separate experiment per group – each group’s own representativeness target and donor pool: cluster. This is the design’s objective, not a constraint: each cluster k reconstructs its own predictor mean \(\overline{X}_k\). It is the SYNDES analog of that estimator’s arm (Synthetic Design (SYNDES)), but baked into the objective, not run as separate solves – which is why, in MAREX, the restrictions compose with cluster (they apply within each cluster) and the per-cluster cardinality m_min / m_max is the stratum quota.

  • One design with geographic / forcing limits within those clusters: the restriction suite above (force in/out, border conflict, size band, donor exclusion).

So a stratum quota is not an alternative to cluster; in MAREX it is a property of the clustering. Drop cluster (a single global cluster) and you recover one design against the whole-population mean – a different estimand, not a constrained version of the clustered one.

Choosing the treated markets: overview#

This section is for the analyst who will run MAREX and explain its choice. The next one, “Choosing the treated markets: theory”, gives the optimization behind it.

Within each cluster, MAREX looks for the treated markets whose weighted average looks like the cluster as a whole before the experiment, while the remaining markets, weighted, look like it too. Both composites then stand in for the same thing, the cluster’s average market, so a gap that opens between them after treatment is the effect. The exact backend (backend="exact") finds that treated set and proves no other set is better, in three layers.

  • Pricing a set. Once the treated markets are named, the best weights come from two small fits, one per side, solved exactly in microseconds. The set’s score is the sum of the two fits, \(f(\mathcal{S})\) below.

  • Trying every set. Up to twenty million candidate sets in a cluster, the search considers them all, with a shortcut: a formula gives a floor under each set’s treated fit without solving it, and a set whose floor is already above the best score found so far is skipped. Every set is either fitted or skipped on a floor that proves it cannot win.

  • Branch and bound. Above twenty million sets the search reasons about whole families at once: every set that treats some given markets and avoids others. It computes a floor under every set in a family from a relaxation, the problem with its rules loosened so that every set in the family still fits, and discards the family when that floor is above the best set found so far. A family that survives is split, one more market decided at a time, and a small one is priced set by set (“Large clusters: a certified branch and bound”).

What the result lets you say: the design has the smallest pre-period imbalance among all treated sets the cardinality and restrictions allow, and with certify=True the result carries the proof, a floor under every alternative beside the design’s own score (“Certifying the search”).

Choosing the treated markets: theory#

This section is for readers who want the optimization behind the search: the split of the design into two fits, the bounds that rule out treated sets, the certificate, and the branch and bound above the candidate limit. None of it is needed to run MAREX; “Choosing the treated markets: overview” above is enough for that.

What the search reports#

With backend="exact" every result also says how its search reached the design. res.clusters[label].search summarises one cluster’s search and res.search merges them; both are None with backend="mip". The summary has the same size whatever the number of candidates, about 850 bytes. Filling it adds about 0.3 percent to the function calls a search makes on 40 GeoLift markets at four treated, and no measurable time.

Each admissible treated set gets exactly one fate, and res.search.counts counts them:

  • cleared_by_bound: its closed-form bound already cleared the best design found so far, so it was never fitted;

  • cleared_by_coupling: the weakly_targeted design’s hull-distance bound cleared it;

  • cleared_by_treated_fit: its treated-side fit alone cleared the incumbent;

  • without_control: donor exclusions left it no control market;

  • cleared_by_fit_cutoff: the coupled fit stopped once its certified bound cleared the incumbent;

  • fitted: fitted to completion.

The six add up to candidates. Above the candidate limit the branch and bound reports the bounds it computed (bnb_bounds), the sets it priced one at a time (bnb_priced), the designs it fitted, its work in the cost model’s seconds (bnb_work) and the nodes its work limit left open (bnb_open, zero when it proved its design optimal).

res.search.clearance is a histogram of how far each discarded set’s bound sat above the incumbent, relative to it, in log-spaced bins; quantile(q) reads it. Sets that cleared by a wide margin say the bound did its work; a mass of small clearances says other designs came close to the one chosen. res.search.numerics counts the bounds that fell back to zero and the coupled fits’ steps and largest stopping gap; res.search.timing splits the seconds between pricing, fitting and the branch and bound.

The candidate count, the kind of search and the fallbacks are facts about the problem. The other counts depend on the order the search visits candidates in, and the times on the machine, so they are reported, not pinned. res.search.order_independent() returns the first group. Two conditions warn with a MAREXSearchWarning: a bound that fell back to zero, which leaves those sets searched on no bound (markets that duplicate one another cause it), and, with certify=True, a certificate that left sets unresolved, which marks a possible tie.

res = MAREX(MAREXConfig(..., backend="exact")).fit()
c = res.search.counts
c.candidates, c.cleared_by_bound, c.fitted
res.search.clearance.quantile(0.5)

Large clusters: a certified branch and bound#

Eight treated among forty markets is 77 million treated sets, ten is 848 million, and the number grows as \(\binom{N}{m}\). Above twenty million sets in a cluster, backend="exact" stops pricing every set and searches the standard, penalized and unit_penalized designs by branch and bound. Nothing changes at or under the limit, and nothing needs to be set: the search returns the optimum the enumeration would, and with certify=True the certificate of the previous section, which it produces as it goes.

A branch and bound organises the treated sets as a tree. A node fixes some markets into the treated set (\(\mathcal{I}\)) and some out of it (\(\mathcal{O}\)), and leaves the rest free (\(\mathcal{F}\)); below it lie the admissible treated sets of size \(m\) that contain \(\mathcal{I}\) and avoid \(\mathcal{O}\). At each node the search does one of two things. It prices the node’s sets one at a time, exactly as the enumeration prices a cluster, when there are few enough of them. Otherwise it computes a bound, a number no design below the node can go under: if the bound is above the best design found so far, the whole node is discarded unseen, and if not, the node is split on one free market into a node that treats it and one that does not. It bounds a node only when pricing the node’s sets would take longer than computing the bound, at about a microsecond per set and a few seconds per bound at forty markets.

The bound has to see the products of the weights. Relaxing the choice of treated markets to fractions between zero and one is not enough: spreading the weights evenly over every market is an average of designs, and it scores exactly zero, because the centred predictors sum to zero. The bound used here, persp+N, works with the matrix of weight products. For each side of the design (treated and control) it keeps the weights \(\mathbf{x}\) on the markets that side can use and a matrix \(\mathbf{X}\) standing for \(\mathbf{x}\mathbf{x}^\top\), with \(z_j\) the fraction of market \(j\) that is treated, and requires

\[\begin{split}\begin{pmatrix} 1 & \mathbf{x}^\top \\ \mathbf{x} & \mathbf{X} \end{pmatrix} \succeq 0, \qquad \mathbf{X} \ge 0, \qquad \mathbf{X}\mathbf{1} = \mathbf{x}, \qquad \mathbf{1}^\top \mathbf{x} = 1, \qquad \begin{pmatrix} z_j & x_j \\ x_j & X_{jj} \end{pmatrix} \succeq 0, \qquad x_j \le z_j,\end{split}\]

for every free market (with \(1 - z_j\) in place of \(z_j\) on the control side), and \(\sum_{j \in \mathcal{F}} z_j = m - |\mathcal{I}|\). The objective is \(\langle \mathbf{G}, \mathbf{X} \rangle + \mathbf{a}^\top \mathbf{x}\) per side. Every design meets every constraint with \(\mathbf{X} = \mathbf{x}\mathbf{x}^\top\), so the minimum is a lower bound on every design below the node. The 2x2 blocks are the perspective constraint \(X_{jj} z_j \ge x_j^2\) of Frangioni and Gentile (2006) and Günlük and Linderoth (2010): a market carrying weight must be treated, and its weight concentrates the second moments on few markets. With the semidefinite block they form the relaxation of Dong, Chen and Linderoth as given in Wei, Gómez and Küçükyavuz (2022, eq. 24), the primal of the best-diagonal program of Zheng, Sun and Li (2014). \(\mathbf{X} \ge 0\) is the entrywise half of the doubly nonnegative relaxation of Burer’s (2009) completely positive form, and it carries most of the bound: dropping it from the full doubly nonnegative relaxation took that bound from 99 percent of the optimum to 29.5 percent at the root of GeoLift’s tree with four treated.

A solver returns the bound in floating point, and the search does not take it on trust. It reads the solver’s multipliers, clips each to the sign its constraint requires, and assembles the Lagrangian they define in exact arithmetic from the exact Gram matrix of the data. Every design satisfies \(\operatorname{tr}\begin{pmatrix} 1 & \mathbf{x}^\top \\ \mathbf{x} & \mathbf{x}\mathbf{x}^\top \end{pmatrix} \le 2\) on each side, so the bound is the Lagrangian’s constant, plus twice the smallest eigenvalue of its matrix when that is negative, plus the \(m - |\mathcal{I}|\) smallest coefficients of the treated fractions. The smallest eigenvalue is bounded below exactly, from a floating-point Cholesky factor of the shifted matrix and its exact residual. This is the rigorous bound of Jansson, Chaykin and Keil (2007) for semidefinite programs, with exact arithmetic in place of their interval arithmetic: any multipliers give a valid bound, and accurate ones a tight one. A node whose solve fails keeps the bound it had, and is split or priced like any other.

Every set the search discards is discarded on a certified lower bound above a certified upper bound on a design already found, so it returns the enumeration’s optimum, and its certificate lists the sets that might tie it, as the certificate above does. The restrictions carry over: forced and forbidden units start the tree, a market treated at a node sends its conflict partners out of the treated set, stratum quotas and forced units decide how many sets a node holds, and every set priced passes the cluster’s rules. Donor exclusions can only raise a control fit, so every bound holds without them and every fitted design honours them.

On one core, with GeoLift’s 40 markets as one cluster over 90 days and Walmart’s 45 stores as the paper’s design sees them, the search returns the enumeration’s design every time, certified with a relative gap of at most \(1.5 \times 10^{-12}\) and no unresolved set:

The branch and bound against enumeration with its limit raised#

Panel

Treated

Treated sets

Branch and bound

Bounds computed

Enumeration

GeoLift

8

76,904,685

41.9 s

7

105.1 s

GeoLift

10

847,660,528

36.9 s

8

1,866.7 s

Walmart

6

8,145,060

14.0 s

2

7.5 s

Walmart

8

215,553,195

211.3 s

22

305.2 s

Below about ten million sets, enumeration is the faster of the two, which is why it keeps every cluster under the limit. Through MAREX, GeoLift with ten treated, which the exact backend refused before, takes 35.5 seconds with certify=True. The runs are benchmarks/studies/marex_bnb/library_search.py.

The search counts its work in the same cost model that decides where it spends time: a microsecond per set priced and \(1 + 3 (k / 40)^3\) seconds per bound with \(k\) markets on the larger side (on the treated side when the bound has one block, below). It stops at 600 such seconds, and does not start a pricing or a bound that would take it past them. Every search in the table above closes well inside the limit (Walmart with eight treated does about 150 of them). A search stopped there returns the best design it found and warns with a MAREXSearchWarning giving its relative gap. Its certificate stays valid: the lower bound is the least bound among the nodes it left open, so the certificate brackets the optimum, and res.search.counts.bnb_open counts those nodes. Only a search stopped before it found any design raises.

Memory binds before time. An interior-point solver carries, for a semidefinite block of order \(s\), a dense scaling block with \((s (s + 1) / 2)^2\) entries in the system it factors each step (Goulart and Chen, Clarabel, sec. 2.5; the Schur complement of Todd 2001, Theorem 5.3), so a bound’s memory grows as the fourth power of the markets on a side. One root bound measured 0.11 GB at 40 markets and 0.86 GB at 70; the fitted law puts 150 markets at about 15 GB and 210 at about 56. Before each bound the search estimates its memory against three quarters of the memory the process may use (the tighter of its control group’s limit and the machine’s memory).

When the two blocks will not fit, the search bounds the node on the treated block alone. The two fits share nothing once the treated set is named, and a design’s controls lie among the markets the node has not treated, where a fit is no worse than over any subset of them, so the treated block’s bound plus the control fit over all those markets, certified like a fitted design, is still a bound below every design in the node. It is weaker, by one to two points of the optimum at the root, and needs about half the memory: 0.06 GB at 40 markets and 0.64 GB at 80, about 7.4 GB at 150. The search returns the same design either way. Measured with each form at every node, the one block took twice as long on GeoLift with ten of forty markets treated, tied at forty to fifty markets elsewhere, and was 3.8 times faster at eighty. The search raises a MlsynthConfigError naming the estimate only when even the one block will not fit. A child uses fewer markets than its parent, so in practice only the root can be refused, before the search has done any work. Clusters or restrictions on the treated markets bring each search back within it; the treated count does not, since every market can still be treated. The weakly-targeted design couples the treated and control fits, which the bound splits, so it is refused above the limit. The tests are in mlsynth/tests/test_marex_branch.py, and the study that chose the bound, with the papers behind each piece, is in benchmarks/studies/marex_bnb.

Inference#

When inference=True with blank_periods > 0, the last few pre-experiment periods are held out as blanks: the design never sees them, so there the synthetic treated minus synthetic control is an out-of-sample error, and its distribution calibrates inference for the post-period effect. MAREX reports a permutation p-value for the global null of no effect, per-period p-values, and a pointwise confidence band from split-conformal prediction (Lei et al. 2018; Vovk, Gammerman, and Shafer 2005), all on MAREXInference. Every one of them rearranges only the blank and post periods, not all periods as in Chernozhukov, Wüthrich and Zhu (2021). alpha sets the level of the bands and of the power curve, and T_post the number of post periods the inference reads, all of them by default.

The global p-value is Abadie and Zhao’s eq. (17). Its statistic is the mean absolute post-period effect, compared with the same statistic on every way of choosing \(T_1\) of the \(T_b + T_1\) blank and post periods. When those arrangements number at most max_combinations MAREX uses all of them, the observed one included, and the p-value is exactly the paper’s. Otherwise it draws max_combinations of them, \(B\), and counts the observed arrangement as one more, reporting \((1 + k)/(1 + B)\) with \(k\) the draws at least as large. That keeps the test valid at any \(B\) (Phipson and Smyth 2010), and the p-value never reaches zero.

Under Assumption 5 the test is exact (Abadie and Zhao’s Theorem 2).

Exchangeability fails when the contrast is serially correlated, which weekly sales usually are. A random subset scatters the post periods across the window, while the observed post block is consecutive, and a run of correlated errors moves a consecutive block further than a scattered one. The test then rejects a true null too often. global_test="cyclic" replaces the subsets with Chernozhukov, Wüthrich and Zhu’s (2021) moving blocks: the \(T_b + T_1\) cyclic shifts of the blank and post contrast, each contributing the block of \(T_1\) consecutive periods that ends the shifted series. A shift keeps neighbouring periods together, so the test needs the contrast to be stationary and weakly dependent, not exchangeable. The statistic is their Remark 1’s, the magnitude of the mean post contrast \(|\bar u_{\text{post}}|\), in place of eq. (16)’s mean magnitude. Against an effect of one sign it is the more powerful of the two, and post contrasts that alternate in sign and cancel give it no evidence. The p-value is a multiple of \(1/(T_b + T_1)\), at least that, and draws nothing, so max_combinations and random_state do not apply. The option changes the global p-value of every contrast MAREX tests (each cluster, the aggregate, and with unit_inference=True each treated unit and eq. (11)’s aggregate) and nothing else; the joint test over treated units keeps its subsets.

res = MAREX(dict(df=df, outcome="y", unitid="unit", time="week",
                 T0=104, m_eq=3, inference=True,
                 global_test="cyclic")).fit()
res.globres.inference.global_p_value

The marex_global_test_dependence study measures both tests on Abadie and Zhao’s Section 5 panel with serially correlated factors, 2,000 draws at 104 and 156 pre-experiment weeks and 12 post weeks. Rejection rates of a true null at 5%:

Factors

"subsets", eq. (17)

"cyclic"

independent

5.2 / 5.2

4.7 / 2.7

AR(1) 0.5

7.8 / 7.8

6.5 / 4.0

AR(1) 0.8

16.9 / 16.5

8.6 / 6.7

AR(1) 0.95

31.9 / 27.2

16.6 / 12.5

random walk

33.6 / 31.5

24.8 / 18.2

With independent factors the cyclic test detects an effect of one error standard deviation 39.7% and 38.6% of the time, against 23.8% and 25.0% for eq. (17). Neither test holds its level near a unit root: the moving blocks need the contrast to be stationary, and a drift that the synthetic control does not cancel is not.

The pointwise band lays a constant half-width \(\hat q\) around each post period’s effect. \(\hat q\) is the \(\lceil (T_b + 1)(1 - \alpha) \rceil\)-th smallest absolute blank effect, the split-conformal quantile. Abadie and Zhao’s eq. (18) uses the empirical quantile, the \(\lceil T_b (1 - \alpha) \rceil\)-th, and their Theorem 3 shows pointwise coverage tending to \(1 - \alpha\) as the fitting and blank windows grow, provided treatment does not change the noise distribution. The split-conformal rank is never smaller. It gives coverage of at least \(1 - \alpha\) in finite samples when the blank and post errors are exchangeable, and the two ranks agree as \(T_b\) grows. On untreated panels in #744 (40 units, three factors autocorrelated at 0.8, a nominal 0.05), the empirical rank gave false-positive rates of 0.097, 0.074 and 0.066 at 19, 24 and 48 blank periods, and the split-conformal rank 0.055, 0.032 and 0.035.

The per-period p-value is mlsynth’s own, the one the band inverts: \((1 + k)/(T_b + 1)\), with \(k\) the blank effects at least as large in magnitude as that period’s. The band excludes zero exactly when it is at most \(\alpha\).

The split-conformal rank exists only when it is at most \(T_b\), which needs \(T_b \ge \lceil 1/\alpha \rceil - 1\) blank periods: 19 at \(\alpha = 0.05\), 9 at 0.10. With fewer, MAREX reports the band as \((-\infty, \infty)\) and warns; the p-values do not depend on it. Left to its default, blank_periods aims for 30% of the pre-period and is held between two floors: at least \(\lceil 1/\alpha \rceil - 1\) blank periods, and at least \(\max(2r, 10)\) fitting periods, where \(r\) is the effective rank of the pre-period panel (Ahn and Horenstein’s 2013 eigenvalue ratio, after removing unit and period means). The second floor is Assumption 2’s \(F \le T_E\) with a margin. A design fitted on too few periods interpolates, and its blank effects then describe the interpolation: in #744, five fitting periods fitted the pre-period exactly and left the cumulative false-positive rate at 0.38. A pre-period too short for both floors, under 29 periods at \(\alpha = 0.05\) and a low-rank panel, is refused with both floors named. A larger alpha or an explicit blank_periods gets past the refusal; an explicit window is taken as given.

The total effect#

The band above bounds one period at a time. How much the intervention added over the whole experiment is the total \(\sum_{t > T_0} \tau_t\) over the \(H\) post periods, and its band is not arithmetic on the pointwise endpoints: summing them assumes the period errors move together, and dividing one by \(H\) assumes they are independent. cumulative_band=True attaches a band for the total to inference.cumulative_band, on the aggregate design and on every cluster, built in one of two ways chosen by cumulative_method. Both are mlsynth’s own constructions: Abadie and Zhao’s Theorem 3 covers one period at a time, and their paper gives no band for the total.

res = MAREX(MAREXConfig(df=df, outcome="y", unitid="unit", time="time",
                        T0=T0, m_eq=3, inference=True, blank_periods=26,
                        cumulative_band=True,
                        cumulative_method="inversion")).fit()
band = res.globres.inference.cumulative_band
band.point, band.lower, band.upper

point is MAREX’s own estimate of the total, the sum of the post-period effects, under both methods.

cumulative_method="resample", the default, splits the blank window into horizon-length sub-windows and resamples their effects into whole paths (cumulative_block sets the block length, cumulative_n_sim the number of paths); the half-width is a quantile of the accumulated totals. It uses only errors observed before the post period, and every band built that way meets a floor set by how many of them there are. With \(n\) blank errors, suppose one error, placed at a period chosen at random, is far larger than the rest. With probability \(H/(n+H)\) it lands in the post window, where no blank error reveals it. A band that covers the total for every exchangeable error law (one whose distribution does not change when the periods are reordered) must then be unbounded with probability at least \(1 - \alpha(n/H + 1)\). A band that stays finite needs \(n \ge H(1/\alpha - 1)\), which is 247 blank periods for a 95% band over 13 post periods.

cumulative_method="inversion" is not subject to that floor because it reads the post outcomes, which an experiment observes: they are the counterfactual errors plus the effect. It inverts the test of Chernozhukov, Wüthrich and Zhu (2021). For a candidate per-period effect \(\theta\), MAREX subtracts \(\theta\) from the synthetic treated series in each post period, refits a synthetic control for that series on every period from the units the design left untreated, and takes the residuals \(\hat u_t\). The statistic is the size of their post sum,

\[S(\hat{\mathbf{u}}) = \frac{\bigl|\sum_{t > T_0} \hat u_t\bigr|}{\sqrt{H}},\]

and the p-value is the share of cyclic shifts of the residual path, taken over the blank and post periods, whose post block gives a statistic at least as large. A cyclic shift moves every residual \(k\) periods along and wraps the end around to the start, so it keeps the order in which errors follow one another. The band is \(H\) times the candidates with \(p > \alpha\), and its endpoints need not be symmetric about point.

Assumptions#

  1. The errors are stationary and weakly dependent over the blank and post periods: their joint law does not drift, and periods far apart are close to independent.

    Remark. This is the condition under which cyclic shifts give an approximately valid test (Chernozhukov, Wüthrich and Zhu 2021, Assumption 2.2). Shifts preserve the serial dependence that random reorderings would destroy, at the price of a coarser p-value: there are only \(T_b + H\) of them, so no p-value falls below \(1/(T_b + H)\). With fewer than \(1/\alpha\) blank and post periods nothing is rejected, the band is \((-\infty, \infty)\), and MAREX warns. A 95% band needs at least 20.

  2. The null is a constant effect, and the band is read as a band for the total.

    Remark. A one-dimensional band needs a one-dimensional null, so the test imposes \(\tau_t = \theta\) in every post period. The statistic depends on the residuals only through their sum, so its observed value depends on the effect path only through the total, and the band covers the total when the effect grows or decays as well. With the default statistic of the general routine, the mean absolute residual, that fails: on simulated MAREX designs with independent factors and a ramp-shaped effect the 95% band covered the total 11% of the time, against 99.7% for the signed sum. What breaks the signed sum is a single error so large that it dominates the post sum; heavy tails alone did not, in a search over Cauchy errors.

  3. The fitting periods stay out of the reference set.

    Remark. MAREX chooses the treated set on the fitting periods, and the refit residuals there are smaller, by about 30% in mean square on the simulated designs, than on the blank and post periods, which behave alike. A reference that includes them understates the spread of a 13-period post sum, by a factor of 1.17 with independent factors and 1.59 with factors autocorrelated at 0.8, and the band under-covers. The refit still matches on every period, which lets it track a trend the controls share; only the reference starts at the blank window. Coverage of the total at 95% / 80%, with the median width of the 95% band, over 300 MAREX designs with 40 fitting, 26 blank and 13 post periods:

    factors

    every period referenced

    blank and post referenced

    independent

    0.950 / 0.773 (25.0)

    0.967 / 0.827 (37.5)

    autocorrelated at 0.5

    0.927 / 0.767 (30.0)

    0.973 / 0.787 (45.0)

    autocorrelated at 0.8

    0.910 / 0.730 (42.5)

    0.967 / 0.777 (60.0)

    trending

    0.983 / 0.833 (32.5)

    0.983 / 0.807 (38.8)

    Whitening the residuals by an autoregressive model before shifting them, which Chernozhukov, Wüthrich and Zhu propose for serially dependent errors, did not close the gap: the difference is between two stretches of the series, not a stationary dependence a filter can remove.

For a cluster, the refit draws on the cluster’s own unselected members, the units its control was built from; for the aggregate design, on every untreated unit. Each candidate costs one simplex fit, and the search brackets and bisects each endpoint, so the band costs some tens of fits per design and no extra solve. cumulative_block, cumulative_n_sim and cumulative_seed apply to "resample" alone. The band is a CumulativeInversionBand; its reference_start is the first blank period and its n_reference is \(T_b + H\). The study behind these choices, with the floor above and its proof, is in benchmarks/studies/marex_cumulative_conformal, and the tests are in mlsynth/tests/test_marex_cumulative_inversion.py.

Each treated unit#

The inference above is for the aggregate contrast. unit_inference=True adds it for each treated unit, on res.globres.unit_level.

res = MAREX(MAREXConfig(df=df, outcome="y", unitid="unit", time="time",
                        T0=T0, m_eq=2, cluster="region", inference=True,
                        design="unit_penalized", xi=1.0, backend="exact",
                        unit_inference=True)).fit()
lvl = res.globres.unit_level
lvl.units["u07"].inference.ci        # the unit's pointwise band
lvl.units["u07"].p_bonferroni        # its p-value, adjusted for the family
lvl.units["u07"].blank_check         # bias of its control on the blank periods
lvl.joint_p_value                    # no effect on any treated unit
lvl.eq11                             # eq. (11)'s aggregate

Each treated unit \(j\) gets its own synthetic control, Abadie and Zhao’s eq. (10): weights \(\mathbf{v}_{\cdot j}\) on the simplex over the untreated units of \(j\)’s cluster, fitted on the fitting periods after the design is chosen. The unit’s contrast,

\[u_{jt} = y_{jt} - \sum_i v_{ij} y_{it},\]

has a blank stretch and a post stretch, so every output of the aggregate’s inference is computed for it: the Theorem 2 p-value, the per-period p-values, the pointwise band and, with cumulative_band=True, the band for the unit’s total. Theorem 2’s argument needs only weights chosen on the fitting periods and Assumption 5, so it covers \(u_j\), with Assumption 3(i) read as an exact fit of unit \(j\)’s own predictors. Abadie and Zhao define the Unit-level design without clusters; with clusters, each unit’s control comes from its own cluster, the pool its cluster’s controls are drawn from.

A unit’s estimated effect is unbiased when its control reproduces it. That needs the unit inside its donors’ convex hull, the set of outcome paths a weighted average of the donors with nonnegative weights summing to one can reach; weights chosen before the post periods; donors the treatment does not reach; and a factor model that holds over the whole timeline. The blank periods test the first and the last of these. They are untreated and the fit never saw them, so a control that reproduces the unit leaves a contrast there that averages zero. blank_check reports that average, its standard error allowing for AR(1) dependence over periods, their ratio t with a two-sided p-value against Student’s t on the effective number of periods, and the average as a share of the unit’s mean absolute level over the blank periods. A large t says the control misses the unit’s level, for instance because no blend of its cluster’s untreated units reaches it. The same bias then enters the unit’s post-period effects, and its band is centred on the wrong number. The check cannot see a bias that begins with the treatment, such as a spillover onto the donors.

Several treated units are several tests: rejecting whenever any unit’s p-value is below \(\alpha\) falsely rejects at \(1 - (1 - \alpha)^m\) with \(m\) units. p_bonferroni is \(\min(1, m\,p_j)\). The joint test asks whether any treated unit has an effect: its statistic is the largest, over units, of the mean absolute post-period contrast divided by the unit’s mean absolute contrast over the blank and post periods, and one arrangement of periods is applied to every unit, so the units’ common factors stay aligned. The divisor is a mean no arrangement changes, so the statistic ignores each unit’s scale and the test remains a permutation test.

Eq. (11) writes the aggregate as the design’s treated weights applied to the per-unit contrasts, \(\sum_j w_j u_{jt}\); eq11 carries its inference beside the shared contrast’s. With cumulative_method="inversion" the band for a total refits the treated series on the untreated units and never reads a control, so eq. (11)’s band for the total is the shared contrast’s, and only its point differs.

On simulated MAREX designs the per-unit inference held where the aggregate’s does. On Abadie and Zhao’s Section 5 panel with exchangeable factors, the per-unit Theorem 2 test rejected 4.7-6.0% of the time at 5%, the pointwise band covered 96-97% at 95%, and the band for a unit’s total 97-99%. A unit’s band was about twice as wide as the aggregate’s, because a single unit’s contrast does not average noise across the treated units. Over a long pre-period with a common cycle and persistent factors, the pointwise band and the band for a unit’s total still held, while the per-unit Theorem 2 test rejected up to three times as often as nominal at 5%; there the band for the total (cumulative_band=True, cumulative_method="inversion") is the per-unit test to read. The blank length mattered little once it cleared the floor.

unit_penalized designs, whose penalty prices each unit’s own fit, gave tighter per-unit bands. On a standard design eq. (11) was four to five times noisier than the shared contrast, since that design never chose its treated units for their own fits; on a unit_penalized design the two agreed. The study is in benchmarks/studies/marex_unit_conformal, and the tests in mlsynth/tests/test_marex_unit_inference.py.

Standardized Post-Fit and Power Analysis#

Every call to MAREX.fit() attaches a SyntheticControlPostFit to res.post_fit. This is the single, estimator-agnostic surface for the diagnostic numbers a consumer of the design typically needs: effects, fit RMSEs, conformal / permutation inference, covariate balance (when covariates were used), and power analysis. It is computed by compute_post_fit() from MAREX’s own synthetic_treated / synthetic_control trajectories and weight vectors, so by construction it agrees with what the underlying optimization produced.

pf = res.post_fit                          # SyntheticControlPostFit
pf.ate, pf.ate_percent, pf.total_effect    # treatment-effect scalars
pf.rmse_fit, pf.rmse_blank, pf.rmse_post   # fit quality, per phase
pf.p_value, pf.ci_lower, pf.ci_upper       # inference (when computed)
pf.covariate_smd                           # treated-vs-control SMD dict
pf.covariate_smd_treated_vs_pop            # treated-vs-population
pf.covariate_smd_control_vs_pop            # control-vs-population
pf.power                                   # PowerAnalysis (see below)

Three Standardized Mean Differences#

When covariates=[...] is set, the post-fit reports the three covariate balance diagnostics that match the structure of Abadie & Zhao’s objective. Each is a per-covariate signed dict (covariate_smd_*) plus two summary scalars (max absolute SMD, sum of squared SMDs). With \(\bar{\mathbf{x}}\) the population covariate aggregate, \(\mathbf{x}_w \coloneqq \sum_j w_j \mathbf{x}_j\), \(\mathbf{x}_v \coloneqq \sum_j v_j \mathbf{x}_j\), and \(s_m\) the cross-unit standard deviation of covariate \(m\), each comparison is the unit-free vector

\[\mathrm{SMD}_m^{(a,b)} = \frac{x_a[m] - x_b[m]}{s_m}.\]

The three pairs (a, b) reported are:

  • covariate_smd — \((\mathbf{x}_w, \mathbf{x}_v)\): synthetic treated vs synthetic control. The internal-validity check (“is the experiment apples-to-apples?”).

  • covariate_smd_treated_vs_pop — \((\mathbf{x}_w, \bar{\mathbf{x}})\): synthetic treated vs population aggregate. Tracks the first term of MAREX’s objective, \(\|\bar{\mathbf{x}} - \sum_j w_j \mathbf{x}_j\|^2\). Tells you whether the chosen treated group represents the population.

  • covariate_smd_control_vs_pop — \((\mathbf{x}_v, \bar{\mathbf{x}})\): synthetic control vs population aggregate. Tracks the second term of the objective. Tells you whether the control set represents the population.

A rule-of-thumb threshold of \(|\mathrm{SMD}| < 0.1\) is conventionally “well balanced”; below \(0.25\) is acceptable; above is a red flag.

Power Analysis and Minimum Detectable Effect#

Power analysis answers the pre-experiment planning question: given the design I’ve chosen, how large a treatment effect can I detect with high probability? This is the dual of inference: inference asks “is the observed effect distinguishable from noise?”, power asks “what effect sizes would be?”

The paper’s inference is built to stay valid under serial dependence and non-stationarity; staying squarely within that framework, mlsynth also reports a planning-oriented minimum detectable effect – an analytical, AR(1)-inflated Gaussian MDE computed from the same residual series the permutation test draws on. Set inference=True and the result auto-populates res.post_fit.power — blank_periods defaults to max(1, floor(0.3 * T0)) so you do not need to pick a scalar yourself (matching the LEXSCM / SYNDES / PANGEO convention).

Where the noise standard deviation comes from#

Under the linear factor model of Assumption 1, the per-period contrast \(g_t \coloneqq \sum_j w_j y_{jt} - \sum_j v_j y_{jt}\) has expectation zero under the no-effect null. Its sample SD on the blank window \(\mathcal{B}\) (the held-out tail of the pre-period) is the natural estimator of the noise scale:

\[\widehat{\sigma}_{\text{placebo}} = \sqrt{\frac{1}{|\mathcal{B}| - 1} \sum_{t \in \mathcal{B}} \bigl(g_t - \bar g\bigr)^2}.\]

When no blank window is carved out (inference=False) the pre-period gap serves as the placebo proxy. The blank-window estimator is preferred because it uses periods that played no role in fitting the weights — it is honest in exactly the same sense Chernozhukov-Wuthrich-Zhu’s conformal residuals are.

Serial correlation matters#

Synthetic-control gap residuals are virtually always serially correlated: the donor weighting absorbs the level but the persistent components of the factor structure (business cycles, seasonality, slow trends) leak through. Ignoring this systematically under-states the SE at long horizons. We model it as an AR(1) process with lag-1 autocorrelation

\[\widehat{\rho} = \frac{\sum_t g_t g_{t-1}}{\sum_t g_t^2},\]

clipped to \((-0.99, 0.99)\) for numerical safety. The variance of the mean of \(T\) consecutive AR(1) periods, expressed as a multiple of \(\sigma^2\), is the variance inflation factor

\[\mathrm{VIF}(T, \rho) = \frac{1}{T}\!\left(1 + 2 \sum_{k=1}^{T-1}\!\Bigl(1 - \frac{k}{T}\Bigr)\rho^k\right),\]

which collapses to the textbook \(1/T\) when \(\rho = 0\) and grows substantially for \(\rho > 0.3\). The same formula is used by PANGEO’s power module.

The MDE formula#

Combining: the standard error of the mean of \(T\) post-period contrasts under \(H_0\) is \(\mathrm{SE}(T) = \widehat{\sigma}_{\text{placebo}} \, \sqrt{\mathrm{VIF}(T, \widehat{\rho})}\). For a two-sided test at level \(\alpha\) with target power \(1 - \beta\), the minimum detectable effect is

\[\mathrm{MDE}(T) = \bigl(z_{1-\alpha/2} + z_{1-\beta}\bigr) \cdot \widehat{\sigma}_{\text{placebo}} \cdot \sqrt{\mathrm{VIF}(T, \widehat{\rho})}.\]

The corresponding power to detect a given true effect \(\tau\) at horizon \(T\) is

\[\pi(\tau, T) = \Phi\!\Bigl(\frac{|\tau|}{\mathrm{SE}(T)} - z_{1-\alpha/2}\Bigr) + \Phi\!\Bigl(-\frac{|\tau|}{\mathrm{SE}(T)} - z_{1-\alpha/2}\Bigr),\]

which is reported as power_at_observed for each horizon point using the realised \(\widehat{\tau}\).

What the surface looks like#

p = res.post_fit.power                  # PowerAnalysis dataclass

p.headline.mde_absolute                 # MDE at the realised T_post
p.headline.mde_pct                      # ... as % of post-period baseline
p.headline.se                           # implied SE of mean(g_t) over T_post
p.headline.power_at_observed            # power to detect res.post_fit.ate

p.curve                                 # tuple of MDEPoint, one per horizon
for pt in p.curve:
    print(pt.post_periods, pt.mde_absolute, pt.mde_pct, pt.power_at_observed)

p.sigma_placebo                         # σ̂ used (from blank or pre window)
p.serial_correlation                    # ρ̂ AR(1) of the placebo gaps
p.placebo_bias                          # mean placebo gap; 0 for an unbiased design
p.placebo_bias_pvalue                   # is that offset distinguishable from 0?
p.n_placebo, p.sigma_ci                 # periods σ̂ rests on, and its interval
p.baseline                              # mean(synthetic_control) on post window
p.alpha, p.power_target                 # 0.05 / 0.80 by default
p.method                                # "analytical_ar1_mean_gap"

pt = p.headline
pt.critical_value                       # t quantile the MDE is built at
pt.mde_ci                               # MDE at the ends of p.sigma_ci

Test an observed effect against pt.critical_value, not against 1.96. The noise scale is estimated from the placebo window, so the pivot follows a Student-t on that window’s effective degrees of freedom; the Gaussian quantile gives a nominal 5% test a true size above 10% on a short, persistent window.

The default horizon grid covers \(T \in \{1, 2, 4, 6, 8, 12\}\) plus the realised n_post, so the table also doubles as a “how long do I need to run?” answer — pick the smallest \(T\) whose MDE drops below your target effect size.

How far to trust the analytical MDE#

The analytical MDE prices a normal test of a mean from the blank window’s noise. Two things make it optimistic for the test MAREX reports. A persistent factor accumulates over a test faster than the AR(1) inflation allows, and the cumulative band is a conformal test inversion, which runs conservative. On 32 simulated panels (50 markets in four regions, 150 pre-periods, AR(0.9) regional factors, six test markets, a 15-period test), post_fit.power reported a 0.95% MDE and the pool’s Newey-West MDE 0.68%, while the cumulative band by inversion had a realized 80%-power MDE of 1.55% (#947). Use the analytical MDE to compare designs, and rehearse the test (below) for the size of effect it can find.

Rehearsing the test#

Before launch nothing is treated, so the planned test can be run on the panel itself. With power_placebos=B, MAREX draws \(B\) random subsets of power_share of the markets (85% by default). On each it solves the same design with the launch moved power_horizon periods earlier, and runs the configured cumulative band on the power_horizon pseudo-post periods. Those periods precede the real launch, so the placebo’s estimate is pure error.

Divide placebo \(b\)’s band \([L_b, U_b]\) by its test composite’s total over the pseudo-post periods, \(S_b\), giving \([\ell_b, u_b] = [L_b / S_b, U_b / S_b]\). A relative effect \(e\) on the test markets (\(e = -0.02\) is a 2% drop) moves the total by \(e S_b\), and both cumulative bands move with an additive change in the post period, so one band per placebo prices every effect:

\[\pi(e) = \frac{1}{B} \sum_{b=1}^{B} \mathbf{1}\bigl\{ 0 \notin [\ell_b + e,\; u_b + e] \bigr\}.\]

The MDE for a drop is the smallest \(m \ge 0\) past which \(\pi(-m') \ge 1 - \beta\) for every \(m' > m\); the MDE for a lift reads \(\pi(m')\) the same way.

res = MAREX(dict(df=df, outcome="sales", unitid="market", time="week",
                 T0=150, m_eq=6, design="penalized", lambda2=1e-3,
                 backend="exact", power_placebos=30,
                 power_horizon=15)).fit()
r = res.power_rehearsal                 # MAREXPowerRehearsal
r.mde_drop                              # e.g. 0.018: a 1.8% drop at 80% power
r.power_at(-0.02)                       # power against a 2% drop
r.estimates.std()                       # spread of the placebo errors
r.lower, r.upper                        # each placebo's band, relative

On the 32 panels above, with 30 placebos of 85% of the markets each, the rehearsal predicted a 1.80% MDE against the realized 1.55%, and the realized test rejected at each panel’s predicted MDE 84% of the time against a target of 80%. The 15% conservatism that remains comes from the pseudo-launch, which shortens the fitting window by power_horizon periods. With placebos of two thirds of the controls the prediction was 2.05%: a placebo’s inversion refit draws on its own untreated markets, so smaller placebos run wider bands.

Three limits. The rehearsal is a planning number for a design setting: its per-panel predictions did not rank panels by their realized error (rank correlation -0.07), so it says how sensitive this kind of test is, not how precise one launch will be. The placebos keep the design, its penalties, the number of test markets, the covariates and the blank window, and drop the clusters, the geographic restrictions and the cost budgets, which choose markets the placebos do not have. It costs power_placebos design solves and as many cumulative bands, so it runs only when asked.

Opting out#

The power computation is wrapped in a try/except in solve_marex() — a power analysis failure (e.g. degenerate residual variance) never breaks the fit, res.post_fit.power is just left as None. To compute power on a non-default horizon grid or significance level, call the free function directly:

from mlsynth.utils.post_fit import compute_power_analysis

alt = compute_power_analysis(
    res.post_fit, alpha=0.10, power_target=0.90,
    post_grid=[2, 4, 8, 16, 32, 52],     # weekly horizons out to a year
)

Monte Carlo: Recovering the Treatment Effect#

The block below replicates the qualitative finding of the paper’s simulation study (Section 5) using mlsynth’s own reimplementation of the linear-factor DGP. A sample is drawn, the design is fit on the pre-period, the treated units realise \(y^I\) in the experiment, and the estimate is compared to the true average effect.

import numpy as np
import pandas as pd
from mlsynth import MAREX
from mlsynth.utils.marex_helpers.simulation import generate_marex_sample

rng = np.random.default_rng(0)

def design_mae(sample, **card):
    J, T = sample.Y_N.shape
    T0 = sample.T0
    df = pd.DataFrame(
        [{"unit": f"u{j}", "time": t, "y": float(sample.Y_N[j, t])}
         for j in range(J) for t in range(T)]
    )
    res = MAREX({"df": df, "outcome": "y", "unitid": "unit",
                 "time": "time", "T0": T0, **card}).fit()
    w = res.globres.treated_weights_agg
    v = res.globres.control_weights_agg
    treated = np.where(w > 1e-8)[0]
    Y_obs = sample.Y_N.copy()
    Y_obs[treated, T0:] = sample.Y_I[treated, T0:]   # experiment realises Y^I
    tau_hat = w @ Y_obs[:, T0:] - v @ Y_obs[:, T0:]
    return np.mean(np.abs(tau_hat - sample.tau[T0:]))

maes, scales = [], []
for _ in range(5):
    s = generate_marex_sample(J=12, T=30, T0=25, rng=rng)
    maes.append(design_mae(s, m_min=1, m_max=11))     # Unconstrained
    scales.append(np.mean(np.abs(s.tau[s.T0:])))
print(f"MAE {np.mean(maes):.2f}  vs  effect scale {np.mean(scales):.2f}")

The synthetic-control design recovers the average treatment effect with mean absolute error far below the effect’s own scale (≈ 4.4 vs. ≈ 14, i.e. under a third), the central message of the paper’s Table 2. Over the paper’s full 1000 simulations the error also decreases as more units are allowed into the treated group (the Unconstrained design is best), with the largest gains moving from one to two or three treated units.

Note

This is a Path-B replication: it reproduces the simulation study’s conclusions from public DGPs and mlsynth code, with no dependency on the authors’ replication package. It is locked in as mlsynth.tests.test_marex_replication.

Empirical Application: Walmart (Placebo Experiment)#

We replicate the paper’s empirical illustration (Section 4) on the Walmart store-sales panel (basedata/walmart_weekly_sales.csv): weekly sales for 45 stores over 143 weeks (Feb 2010 – Oct 2012). Following the paper, we design a placebo experiment with a fictitious intervention at week 129: \(T_0 = 128\) pre-experiment weeks, of which the first \(T_E = 100\) are the fitting period and the last 28 are blank, leaving 15 experimental weeks. The design uses the constrained formulation with \(m = 2\) treated stores, uniform weights, and predictors normalised to unit variance (standardize).

import pandas as pd
from mlsynth import MAREX

df = pd.read_csv(
    "https://raw.githubusercontent.com/jgreathouse9/mlsynth/"
    "refs/heads/main/basedata/walmart_weekly_sales.csv"
)

res = MAREX({
    "df": df, "outcome": "sales", "unitid": "store", "time": "week",
    "T0": 128, "blank_periods": 28, "T_post": 15,   # TE=100, 28 blank, 15 post
    "m_eq": 2,                  # constrained design, two treated stores
    "design": "standard",
    "standardize": True,        # unit-variance predictors (paper's normalisation)
    "inference": True,
    "display_graph": True,
}).fit()

print("treated stores:", res.treated_units)              # [1, 15]
print("placebo p-value:", round(res.globres.inference.global_p_value, 3))

Because the intervention is a placebo (no real effect), a correct design should produce synthetic treated and control units that track closely and an estimated effect near zero. mlsynth reproduces exactly that — and the paper’s headline number:

Walmart placebo design (m = 2)#

Quantity

mlsynth

Paper (Section 4)

Pre-fit RMSE / mean sales

2.2%

small (close tracking)

Experimental ATT / mean sales

-1.0%

near zero

Placebo permutation p-value

0.937

0.933

Confidence band covers zero

yes (all post weeks)

yes

The synthetic treated and control units track to within ~2% of mean sales over the fitting and blank periods, the estimated placebo effect is ~1% of sales, and the permutation test fails to reject the null of no effect (\(p = 0.937\), matching the paper’s \(0.933\)) — exactly the “no spurious effect” result a good design should deliver on a placebo.

Note

This uses the exact MIQP (relaxed=False, the default) with standardize=True; the unit-variance normalisation is essential here because Walmart stores differ enormously in sales level, and without it the level differences dominate the match. The solve takes roughly a minute with the open-source SCIP solver (the paper used commercial Gurobi). With backend="exact" the whole fit takes about 0.03 seconds and returns the same two stores (Choosing the treated markets: overview).

Note

Cross-validated against the authors’ code, without Gurobi. The durable benchmark benchmarks/cases/marex_walmart.py runs this design on the full 45-store panel, additionally matching on the four store-level covariates the reference data carries (temperature, fuel price, CPI, unemployment – the R code’s “few covariates” configuration), and compares it cell-by-cell against a live run of SCDesign’s own design code. The authors’ published routine is a Gurobi non-convex MIQP, but their constrained (cardinality-\(K\)) design – Synthetic_Experiment_Cardinality_Constraint, which enumerates the partitions of size \(\le K\) and solves each treated/control synthetic control through the open quadprog path – is exactly the design m_eq solves and needs no commercial solver. MAREX (SCIP) and SCDesign (quadprog) select the same two treated stores with treated weights agreeing to \(2\times10^{-4}\), the same pre-period fit, and the same placebo effect to \(10^{-4}\) of mean sales; see benchmarks/reference/marex_walmart/ for the captured reference run.

Correspondence with the Authors’ Code#

The authors’ R replication code (Random_Data_Generator.R, Synthetic_Experiments.R, Different_optimization_methods.R) maps directly onto mlsynth’s implementation, which was checked against it:

Authors’ R ↔ mlsynth MAREX#

Authors’ R

mlsynth

DGP (Random_Data_Generator.R)

generate_marex_sample()

Formulation (5), Gurobi non-convex QCQP

design="standard" (MIQP with binary z; same optimum)

Penalization formulation

design="penalized" (identical \(\lambda\) distance penalty)

Cardinality formulation

m_eq / m_min / m_max

Cardinality routine (every partition, quadprog on each side)

backend="exact" (the same enumeration, with a bound that skips most of the fits)

Predictors \(\mathbf{x} = [\mathbf{y}^E ; \mathbf{z}]\)

covariates=[...] (matched on pre-outcomes + covariates)

“treated = smaller set” swap

applied where either labelling is optimal; not where m_eq, a cost bound or an asymmetric objective already settles it (see below)

Exact permutation test (sum statistic)

permutation inference (mlsynth defaults to a mean statistic / sampled permutations)

The reference’s labelling rule is kept, and only where it applies. The authors’ base program in SCdesign_LazyRun.R is two simplex rows and the disjointness \(\mathbf{w}'\mathbf{v} = 0\), solved as a continuous non-convex QCQP (vtype = 'C') with no selection variable at all. Exchanging the two weight vectors is feasible there at the same cost, so the labelling is free and “treatment regions are the smaller set of regions” is the tie-break that program needs. Their cardinality routine enumerates treated sets of every size from 1 up to \(K\), defaulting to \(\lfloor N/2 \rfloor\), so the treated group is never the majority there either.

mlsynth’s program can be either case, so the rule is applied on a test of which. Swapping \(\mathbf{w}\) and \(\mathbf{v}\) is feasible at the same cost when the objective is symmetric in them (the standard design), no exact m_eq is set, the complement’s size also satisfies m_min and m_max, and there is no cost bound or geographic restriction – those price and restrict \(\mathbf{w}\) and \(\mathbf{z}\) alone. That case is ordinary: m_min=1, m_max=N-1 is exactly it, and it is the cell the Section 5 Monte Carlo compares against the authors’ numbers. There the convention decides, as it does for them.

Where the program does settle the labelling it is read off \(\mathbf{w}\), the variable the program constrains. Applying the convention there instead returned three treated markets for m_eq = 6 on twelve, and they were the control synthetic’s markets – the group \(\sum_j z_j = m\) does not constrain and \(\sum_j c_j w_j \le B\) does not price. It also negated the estimated effect, which is the treated synthetic minus the control one.

Driving mlsynth’s solve_design on the authors’ exact DGP and predictor matrix recovers the average treatment effect to within its scale, and the effect estimate degrades gracefully as the noise SD rises from 1 to 5 to 10 (the figures 2-7 settings) — matching the paper’s qualitative findings.

Note

Two faithfulness details from their R code: rnorm(N, 0, noise.variance) passes the value as a standard deviation, so the figures’ “variance” 1/5/10 are SDs; and the R code uses random population weights \(f_j\) whereas the 2026 paper (and mlsynth) use \(f_j = 1/J\). Also, the unconstrained standard design (formulation 5) is degenerate — many disjoint splits match \(\bar X\) equally well, so the realised design (and hence a single ATE estimate) is solver-dependent; the cardinality-constrained design is the stable, recommended choice.

Example#

import pandas as pd
from mlsynth import MAREX

# long panel: one row per (market, period)
res = MAREX({
    "df": df, "outcome": "revenue", "unitid": "market", "time": "week",
    "T0": 40,            # 40 pre-experiment weeks
    # Equivalently: pass a 0/1 column marking the experiment window
    # "post_col": "in_experiment",
    "m_eq": 2,           # treat exactly two markets
    "design": "standard",
    "inference": True,   # blank_periods defaults to floor(0.3 * T0) = 12
    "display_graph": True,
}).fit()

print("treated markets:", res.treated_units)
print("global p-value:", res.globres.inference.global_p_value)
for label, c in res.clusters.items():
    print(label, c.unit_weight_map["Treated"])

Verification#

Validated against Abadie & Zhao’s Section 4 Walmart application (their reference code is jinglongzhao2/SCDesign): on a 10-store subset of walmart_weekly_sales.csv MAREX’s exact MIQP designs a placebo experiment that tracks closely pre-period (pre-fit RMSE ~2.7% of mean sales, matching LEXSCM) and yields a placebo effect indistinguishable from zero (~1% of mean, CI covering zero) – the paper’s “no spurious effect” result. This is an independent commit-stamped check (MAREX’s own optimizer) complementing the LEXSCM Walmart benchmark. See MAREX — Abadie & Zhao (2026) Walmart design (live R cross-validation); run it with python benchmarks/run_benchmarks.py marex_walmart.

The exact backend (backend="exact") is checked separately against a brute-force enumeration of every treated set in R, across every design and objective option it solves: benchmarks/cases/marex_exact_bruteforce.py, described in MAREX — Abadie & Zhao (2026) Walmart design (live R cross-validation).

On Abadie and Zhao’s own 1000 Section 5 simulations, regenerated in Python from their seeds, the base design is unbiased: the mean signed error is within one standard error of zero at every cardinality, and the blank-period gap post_fit.power.placebo_bias tracks the bias each single design carries with a slope near one. benchmarks/cases/marex_section5_bias.py, described in MAREX — Abadie & Zhao (2026) Walmart design (live R cross-validation).

The branch and bound above twenty million sets in a cluster has a regression gate: benchmarks/cases/marex_bnb_work.py pins the bounds it computes and the sets it prices on Walmart and two GeoLift subsets, beside the design and its certificate. See Branch and bound work — SYNDES and MAREX regression gates.

Note

The benchmark uses the exact MIQP (free SCIP), not the relaxed continuous-z mode: the relaxation shares the design objective but drops the integrality that defines the selection, so its top-m rounding is degenerate and non-deterministic for small treated counts. The authors’ full 45-store MIQP uses Gurobi, so the validator is Path A on a subset, not a live R cross-validation.

Core API#

MAREX: Synthetic Controls for Experimental Design (Abadie & Zhao 2026).

MAREX designs an experiment on aggregate units (e.g. markets): using only pre-experimental data it selects which units to treat (treated weights w) and which untreated units form the synthetic control (control weights v), on the simplex and disjoint (a unit is treated or a control, never both). The synthetic treated and synthetic control units are built to reproduce population predictor means, so their post-period difference estimates the average treatment effect. Optional clustering treats one (or a few) units per cluster; optional blank-period placebo inference yields p-values and confidence bands.

class mlsynth.estimators.scexp.MAREX(config: MAREXConfig | dict)#

Bases: object

Synthetic-control experimental design estimator (Abadie & Zhao 2026).

Parameters:

config (MAREXConfig or dict) – Configuration object. See mlsynth.config_models.MAREXConfig.

Returns:

MAREXResults – Per-cluster and aggregate treated/control weights, synthetic series, the selected treated units, and (optionally) placebo inference.

fit() → MAREXResults#

Run the MAREX design and return MAREXResults.

Configuration#

class mlsynth.config_models.MAREXConfig(*, df: DataFrame, outcome: str, unitid: str, time: str, T0: int | None = None, post_col: str | None = None, cluster: str | None = None, design: str = 'standard', covariates: List[str] | None = None, covariate_weight: float | Dict[str, float] = 1.0, standardize: bool = False, program_type: str = 'MIQP', display_graph: bool = False, beta: float = 1e-06, control_penalty: float = 0.0, lambda1: float = 0.0, lambda2: float = 0.0, xi: float = 0.0, lambda1_unit: float = 0.0, lambda2_unit: float = 0.0, costs: List[float] | None = None, budget: float | Dict[Any, float] | None = None, blank_periods: int | None = None, m_eq: int | None = None, m_min: int | None = None, m_max: int | None = None, exclusive: bool = True, relaxed: bool = False, backend: Literal['mip', 'exact'] = 'mip', certify: bool = False, solver: Any = None, verbose: bool = False, warm_start: List | None = None, time_limit: Annotated[float | None, Gt(gt=0.0)] = None, inference: bool = False, unit_inference: bool = False, T_post: int | None = None, global_test: Literal['subsets', 'cyclic'] = 'subsets', cumulative_band: bool = False, cumulative_method: Literal['resample', 'inversion'] = 'resample', cumulative_block: int = 0, cumulative_n_sim: int = 2000, cumulative_seed: int = 0, to_be_treated: List | None = None, not_to_be_treated: List | None = None, adjacency: DataFrame | None = None, spillover_threshold: float = 0.0, exclude_bordering_donors: bool = False, size_col: str | None = None, min_size: float | None = None, max_size: float | None = None, market_cost_col: str | None = None, max_total_cost: float | None = None, cluster_col: str | None = None, stratum_col: str | None = None, min_per_stratum: Annotated[int | None, Ge(ge=1)] = None, max_per_stratum: Annotated[int | None, Ge(ge=1)] = None, top_K: Annotated[int, Ge(ge=1)] = 1, power_weight: Annotated[float, Gt(gt=0)] = 0.51, fit_weight: Annotated[float, Gt(gt=0)] = 0.49, max_shortlist: Annotated[int, Ge(ge=1)] = 5, alpha: Annotated[float, Gt(gt=0.0), Lt(lt=1.0)] = 0.05, power_target: Annotated[float, Gt(gt=0.0), Lt(lt=1.0)] = 0.8, power_placebos: Annotated[int, Ge(ge=0)] = 0, power_horizon: int | None = None, power_share: float = 0.85, power_seed: int = 0)#

Configuration for the Synthetic Experiment Design estimator (MAREX) in mlsynth.

T0: int | None#
T_post: int | None#
adjacency: pd.DataFrame | None#
alpha: float#
backend: Literal['mip', 'exact']#
beta: float#
blank_periods: int | None#
budget: float | Dict[Any, float] | None#
certify: bool#
cluster: str | None#
cluster_col: str | None#
control_penalty: float#
costs: List[float] | None#
covariate_weight: float | Dict[str, float]#
covariates: List[str] | None#
cumulative_band: bool#
cumulative_block: int#
cumulative_method: Literal['resample', 'inversion']#
cumulative_n_sim: int#
cumulative_seed: int#
design: str#
display_graph: bool#
exclude_bordering_donors: bool#
exclusive: bool#
fit_weight: float#
global_test: Literal['subsets', 'cyclic']#
inference: bool#
lambda1: float#
lambda1_unit: float#
lambda2: float#
lambda2_unit: float#
m_eq: int | None#
m_max: int | None#
m_min: int | None#
market_cost_col: str | None#
max_per_stratum: int | None#
max_shortlist: int#
max_size: float | None#
max_total_cost: float | None#
min_per_stratum: int | None#
min_size: float | 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].

not_to_be_treated: List | None#
post_col: str | None#
power_horizon: int | None#
power_placebos: int#
power_seed: int#
power_share: float#
power_target: float#
power_weight: float#
program_type: str#
relaxed: bool#
size_col: str | None#
solver: Any#
spillover_threshold: float#
standardize: bool#
stratum_col: str | None#
time_limit: float | None#
to_be_treated: List | None#
top_K: int#
unit_inference: bool#
classmethod validate_design_params(values: Any) → Any#
verbose: bool#
warm_start: List | None#
xi: float#

Result Containers#

MAREX.fit() returns a MAREXResults: a dict of per-cluster MAREXClusterDesign objects, the aggregate MAREXGlobalDesign, the MAREXStudy hyperparameters, and (optionally) MAREXInference.

Note

MAREX.fit() returns a DesignResult (the experimental-design family, not an EffectResult): MAREX designs an experiment, so it exposes the standardized design surface – res.report (the realized effect as an EffectResult, the single source for ATT / CI / pre-fit; res.report.att / res.report.counterfactual / …), res.selected_units / res.assignment (treated vs control), res.design_weights, and res.power. The MAREX-specific design detail stays on res.clusters / res.study / res.globres / res.post_fit (the same SyntheticControlPostFit that res.report is built from).

The result also carries what a step after the fit reads, as typed models:

  • res.windows, a DesignWindows: the fitting window res.windows.fit, the blank window res.windows.blank and the post window res.windows.post as slices of the timeline, with T0, fit_end, blank_periods, post_periods and the period labels;

  • res.composite, a DesignComposite: the aggregate treated and control weights by unit label, the units the design drew from, and pool, the untreated units a refit of the treated composite may use as donors;

  • res.cluster_composites: the same for each cluster, keyed like res.clusters;

  • res.inference_settings, a MAREXInferenceSettings: the level, T_post, the global test, the subset draws and seed, and the cumulative band’s options. None when no inference ran.

With these and the panel, every contrast’s inference is reproduced exactly:

from mlsynth.utils.marex_helpers.inference import compute_inference

w, c, s = res.windows, res.composite, res.inference_settings
Y = df.pivot(index="unit", columns="time", values="y").loc[c.members].to_numpy()
wt = np.array([c.treated.get(u, 0.0) for u in c.members])
wc = np.array([c.control.get(u, 0.0) for u in c.members])
pool = Y[[c.members.index(u) for u in c.pool]]
inf = compute_inference(wt @ Y, wc @ Y, w.T0, w.fit_end, w.blank_periods,
                        control_outcomes=pool, **s.inference_kwargs())
inf.global_p_value == res.globres.inference.global_p_value   # True

res.metadata keeps T0 and blank_periods as mirrors of res.windows.

Frozen dataclass containers for the MAREX (synthetic experimental design) pipeline.

Implements the containers for:

Abadie, A., & Zhao, J. (2026). “Synthetic Controls for Experimental Design.”

MAREX designs an experiment on aggregate units: it chooses treated weights w and control weights v (on the simplex, disjoint via w_j v_j = 0) so the synthetic treated and synthetic control units reproduce population predictor means. All containers are frozen (immutable) per the repository convention; inference, when requested, is computed up front and embedded.

class mlsynth.utils.marex_helpers.structures.MAREXBlankCheck(mean: float, se: float, t: float, p_value: float, rho: float, n: int, relative: float)#

Bases: object

Whether a treated unit’s contrast averages zero over the blank periods.

The blank periods are untreated and unused by the fit, so a unit whose synthetic control reproduces it has a contrast there that averages zero. A mean far from zero says the control misses the unit’s level: the unit lies outside its donors’ hull, or its fit on the fitting periods did not carry over. The same bias then enters the unit’s post-period effects.

Parameters:

  • mean (float) – Mean of the contrast over the blank periods.

  • se (float) – Its standard error, sd * sqrt(f(n, rho)) with sd the sample standard deviation and f the AR(1) variance of a mean of n periods relative to one period’s variance, (1 + 2 sum_{k<n} (1 - k/n) rho^k) / n. NaN with one blank period.

  • t (float) – mean / se; zero for a contrast identically zero, and infinite for a nonzero constant one.

  • p_value (float) – Two-sided, against Student’s t on n_eff - 1 degrees of freedom, n_eff = n (1 - rho) / (1 + rho) capped at n and floored at 2. An approximation: the contrast is serially dependent and the AR(1) model for it is fitted to the same n periods.

  • rho (float) – Lag-1 autocorrelation of the contrast, demeaned, corrected for its small-sample bias, and clipped to (-0.99, 0.99).

  • n (int) – Number of blank periods.

  • relative (float) – mean over the unit’s mean absolute outcome on the blank periods: the bias as a share of the unit’s level. NaN at a zero level.

mean: float#
n: int#
p_value: float#
relative: float#
rho: float#
se: float#
t: float#
class mlsynth.utils.marex_helpers.structures.MAREXCertificate(lower: float, upper: float, clusters: Dict[str, ClusterCertificate])#

Bases: object

A certified optimality gap for an exact-backend design.

lower <= f* <= f(design) <= upper, with f* the optimal objective and f(design) the returned design’s, both exact for the stored data (mlsynth.utils.marex_helpers.certify). Clusters are searched apart, so the bounds are the sums of the clusters’ bounds, rounded outward.

Parameters:

  • lower, upper (float) – The bounds.

  • clusters (dict of {str: ClusterCertificate}) – Each cluster’s certificate, its unresolved treated sets as unit labels.

clusters: Dict[str, ClusterCertificate]#
property gap: float#

the design is this close to optimal.

Type:

upper - lower, rounded up

lower: float#
classmethod of(clusters: Dict[str, ClusterCertificate]) → MAREXCertificate#

The certificate of the whole design, from its clusters’.

upper: float#
class mlsynth.utils.marex_helpers.structures.MAREXClusterDesign(label: str, members: List[Any], cardinality: int, treated_weights: ndarray, control_weights: ndarray, selection_indicators: ndarray, synthetic_treated: ndarray, synthetic_control: ndarray, pre_treatment_means: ndarray, rmse: float, unit_weight_map: Dict[str, Dict[Any, float]], inference: MAREXInference | None = None, search: ExactSearchDiagnostics | None = None)#

Bases: object

Design and synthetics for a single cluster.

Parameters:

  • label (str) – Cluster label.

  • members (list) – Unit labels in this cluster.

  • cardinality (int) – Number of units in the cluster.

  • treated_weights (np.ndarray) – Treated weights w for this cluster’s column, shape (N,).

  • control_weights (np.ndarray) – Control weights v for this cluster’s column, shape (N,).

  • selection_indicators (np.ndarray) – Binary selection mask z over the cluster’s members.

  • synthetic_treated (np.ndarray) – Synthetic treated outcome over the full timeline, shape (T,).

  • synthetic_control (np.ndarray) – Synthetic control outcome over the full timeline, shape (T,).

  • pre_treatment_means (np.ndarray) – Cluster predictor means used as the matching target.

  • rmse (float) – Pre-treatment fit RMSE (synthetic treated vs control).

  • unit_weight_map (dict) – {"Treated": {unit: w}, "Control": {unit: v}} for non-zero weights.

  • inference (MAREXInference, optional) – Inference for this cluster (None unless requested).

  • search (ExactSearchDiagnostics, optional) – With backend='exact', what this cluster’s search did: each candidate’s fate, the clearances of the discarded ones, the bound fallbacks and the time by phase. None with backend='mip'.

cardinality: int#
control_weights: ndarray#
inference: MAREXInference | None = None#
label: str#
members: List[Any]#
pre_treatment_means: ndarray#
rmse: float#
search: ExactSearchDiagnostics | None = None#
selection_indicators: ndarray#
synthetic_control: ndarray#
synthetic_treated: ndarray#
treated_weights: ndarray#
unit_weight_map: Dict[str, Dict[Any, float]]#
class mlsynth.utils.marex_helpers.structures.MAREXGlobalDesign(Y_full: ndarray, Y_fit: ndarray, Y_blank: ndarray | None, treated_weights_agg: ndarray, control_weights_agg: ndarray, synthetic_treated: ndarray, synthetic_control: ndarray, inference: MAREXInference | None = None, objective: float | None = None, unit_level: MAREXUnitLevel | None = None)#

Bases: object

Aggregated (population-level) design and synthetics.

Parameters:

  • Y_full (np.ndarray) – Observed outcome matrix, shape (N, T).

  • Y_fit (np.ndarray) – Fitting slice, shape (N, T_fit).

  • Y_blank (np.ndarray, optional) – Held-out blank pre-periods, shape (N, Tb) (None if none).

  • treated_weights_agg (np.ndarray) – Cluster-size-weighted aggregate treated weights, shape (N,).

  • control_weights_agg (np.ndarray) – Cluster-size-weighted aggregate control weights, shape (N,).

  • synthetic_treated (np.ndarray) – Aggregate synthetic treated outcome, shape (T,).

  • synthetic_control (np.ndarray) – Aggregate synthetic control outcome, shape (T,).

  • inference (MAREXInference, optional) – Aggregate inference (None unless requested).

  • objective (float, optional) – The design objective at the solution: the fit terms of every cluster, plus any penalty the design adds. Comparable across backends, so it is how a search and the MIQP are checked against each other; SCIP reports it to its own tolerance.

  • unit_level (MAREXUnitLevel, optional) – Inference for each treated unit and eq. (11)’s aggregate, present when unit_inference=True was requested with a blank window.

Y_blank: ndarray | None#
Y_fit: ndarray#
Y_full: ndarray#
control_weights_agg: ndarray#
inference: MAREXInference | None = None#
objective: float | None = None#
synthetic_control: ndarray#
synthetic_treated: ndarray#
treated_weights_agg: ndarray#
unit_level: MAREXUnitLevel | None = None#
class mlsynth.utils.marex_helpers.structures.MAREXInference(treated_effects: ndarray, placebo_effects: ndarray, fulltreated_effects: ndarray, s_obs: float, global_p_value: float, per_period_pvals: ndarray, ci: ndarray, alpha: float = 0.05, cumulative_band: CumulativeConformalBand | CumulativeInversionBand | None = None)#

Bases: object

Placebo/permutation inference for one synthetic treated-vs-control pair.

Parameters:

  • treated_effects (np.ndarray) – Post-period synthetic treated minus synthetic control, shape (T1,).

  • placebo_effects (np.ndarray) – The same contrast on the blank (held-out pre) periods, shape (Tb,).

  • fulltreated_effects (np.ndarray) – The contrast over the whole timeline, shape (T,).

  • s_obs (float) – Observed test statistic: the mean absolute post-period effect under global_test="subsets", the magnitude of the mean post-period effect under "cyclic".

  • global_p_value (float) – Permutation p-value for the global null of no effect. Under global_test="subsets" it is Abadie and Zhao’s eq. (17): it uses every T1-subset of the blank and post periods when there are at most max_combinations of them, and otherwise max_combinations drawn subsets, as (1 + k) / (1 + B). Under "cyclic" the reference is the Tb + T1 cyclic shifts of the blank and post contrast, and the p-value is a multiple of 1 / (Tb + T1). Either way it is never zero.

  • per_period_pvals (np.ndarray) – Per-post-period p-values, shape (T1,): (1 + k) / (Tb + 1), with k the blank effects at least as large in magnitude.

  • ci (np.ndarray) – Split-conformal confidence band over the full timeline, shape (T, 2). The half-width is the ceil((Tb + 1)(1 - alpha))-th order statistic of the absolute blank effects; with fewer than ceil(1 / alpha) - 1 blank periods it does not exist and the rows are (-inf, inf). Rows outside the T1 post periods are NaN.

  • alpha (float) – Two-sided significance level.

  • cumulative_band (CumulativeConformalBand or CumulativeInversionBand, optional) – Band for the TOTAL post-period effect, present only when cumulative_band=True was requested. ci above is pointwise: a constant half-width laid around each period’s effect, which bounds one period at a time. Stacking those endpoints would bound the total only if the period errors moved in lockstep, and dividing one by the horizon only if they were independent. Under cumulative_method="resample" it is a CumulativeConformalBand, which block-resamples the blank-period effects into whole paths and takes the quantile of their totals; its n_scores is the number of blank periods drawn from. Under "inversion" it is a CumulativeInversionBand, the constant effects a test against the post outcomes does not reject; its n_reference is the number of blank and post periods, and its endpoints need not be symmetric about point. Either way point is the sum of treated_effects.

alpha: float = 0.05#
ci: ndarray#
cumulative_band: CumulativeConformalBand | CumulativeInversionBand | None = None#
fulltreated_effects: ndarray#
global_p_value: float#
per_period_pvals: ndarray#
placebo_effects: ndarray#
s_obs: float#
treated_effects: ndarray#
class mlsynth.utils.marex_helpers.structures.MAREXInferenceSettings(*, alpha: float, T_post: int | None = None, global_test: str = 'subsets', max_combinations: int = 1000, random_state: int | None = None, cumulative_band: bool = False, cumulative_method: str = 'resample', cumulative_block: int = 0, cumulative_n_sim: int = 2000, cumulative_seed: int = 0)#

Bases: DesignInferenceSettings

The settings MAREX’s inference ran with, as compute_inference() takes them.

inference_kwargs with the design’s windows and composites reproduces every contrast’s inference exactly.

T_post: int | None#
cumulative_band: bool#
cumulative_block: int#
cumulative_method: str#
cumulative_n_sim: int#
cumulative_seed: int#
global_test: str#
inference_kwargs() → Dict[str, Any]#

The settings as keyword arguments of the estimator’s inference.

max_combinations: int#
model_config: ClassVar[ConfigDict] = {'extra': 'forbid', 'frozen': True}#

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

random_state: int | None#
class mlsynth.utils.marex_helpers.structures.MAREXPowerRehearsal(horizon: int, n_placebos: int, n_units: int, share: float, seed: int, alpha: float, power_target: float, cumulative_method: str, lower: ndarray, upper: ndarray, estimates: ndarray, mde_drop: float, mde_lift: float)#

Bases: object

The planned test, rehearsed where nothing happened.

Each placebo is MAREX’s design re-solved on a random n_units of the markets and launched horizon periods before the real launch, with the configured cumulative band run on its horizon pseudo-post periods. Its band and point are divided by its test composite’s pseudo-post total, so every number below is a share of what the test markets would have done: an effect of -0.02 is a 2% drop.

A relative effect e moves a placebo’s band by e, so the placebo rejects “no effect” at e when 0 is outside [lower + e, upper + e]. power_at(e) is the share of placebos that do.

Parameters:

  • horizon (int) – Periods the planned test runs.

  • n_placebos, n_units (int) – Placebo designs rehearsed, and the markets each draws.

  • share, seed (float, int) – n_units as a share of the panel’s markets, and the draw’s seed.

  • alpha, power_target (float) – The test’s level and the power the MDEs are read at.

  • cumulative_method (str) – The cumulative band the rehearsal ran: the one the real test reports.

  • lower, upper, estimates (np.ndarray) – Each placebo’s band and point estimate, shape (n_placebos,), as shares of its composite’s pseudo-post total. The estimates are errors: nothing happened.

  • mde_drop, mde_lift (float) – The smallest drop, and the smallest lift, as a positive share, past which power stays at or above power_target. inf when unbounded bands keep it out of reach.

alpha: float#
cumulative_method: str#
estimates: ndarray#
horizon: int#
lower: ndarray#
mde_drop: float#
mde_lift: float#
n_placebos: int#
n_units: int#
power_at(effect: float) → float#

Share of placebos that reject “no effect” at the relative effect.

power_target: float#
seed: int#
share: float#
upper: ndarray#
class mlsynth.utils.marex_helpers.structures.MAREXRecommendation(winner: Dict[str, Any], shortlist: List[Dict[str, Any]], pareto: List[int], weights: Dict[str, float], status: str)#

Bases: object

Composite power-vs-fit recommendation over a MAREX solution pool.

Parameters:

  • winner (dict) – The recommended pool entry (lowest composite score among power-feasible designs).

  • shortlist (list of dict) – Pool entries ordered by composite score, truncated to max_shortlist.

  • pareto (list of int) – Indices (into the pool) of the designs on the fit-vs-power Pareto front.

  • weights (dict) – Normalised {"power": pw, "fit": fw} (sum to one).

  • status (str) – "OK" when at least one design has a finite MDE, else "POWER_NOT_ESTABLISHED".

pareto: List[int]#
shortlist: List[Dict[str, Any]]#
status: str#
weights: Dict[str, float]#
winner: Dict[str, Any]#
class mlsynth.utils.marex_helpers.structures.MAREXResults(*, report: BaseEstimatorResults | None = None, assignment: Any | None = None, selected_units: Any | None = None, design_weights: WeightsResults | None = None, power: Any | None = None, metadata: Dict[str, Any] | None = None, windows: DesignWindows | None = None, composite: DesignComposite | None = None, cluster_composites: Dict[str, DesignComposite] | None = None, inference_settings: DesignInferenceSettings | None = None, clusters: Dict[str, MAREXClusterDesign], study: MAREXStudy, globres: MAREXGlobalDesign, post_fit: Any | None = None, pool: List[Dict[str, Any]] | None = None, recommendation: Any | None = None, certificate: MAREXCertificate | None = None, search: ExactSearchDiagnostics | None = None, power_rehearsal: MAREXPowerRehearsal | None = None, **extra_data: Any)#

Bases: DesignResult

User-facing output of the MAREX estimator.

A DesignResult (the experimental-design family): MAREX designs an experiment (it chooses which units to treat), so it populates the standardized design surface – report (the realized effect as an EffectResult, the single source for ATT / CI / pre-fit), selected_units / assignment / design_weights, and power – while the MAREX-specific design detail stays in the typed fields below.

Parameters:

  • clusters (dict of {str: MAREXClusterDesign}) – Per-cluster design (a single "0" entry when no cluster column).

  • study (MAREXStudy) – Design hyperparameters.

  • globres (MAREXGlobalDesign) – Aggregate design and synthetics.

  • post_fit (SyntheticControlPostFit, optional) – Standardized post-fit diagnostics (ATE / total effect / percentage lift / fit RMSEs / inference / covariate SMDs). Computed via mlsynth.utils.post_fit.compute_post_fit(). None only when an estimator failure leaves the result partially constructed. Also mirrored into report (the standardized effect view).

  • certificate (MAREXCertificate, optional) – With certify=True and backend='exact', certified bounds on the optimal objective and the design’s; None otherwise.

  • search (ExactSearchDiagnostics, optional) – With backend='exact', the clusters’ search summaries merged: the candidates each stage decided, the clearances, the bound fallbacks and the time by phase (diagnostics). None with backend='mip'.

  • power_rehearsal (MAREXPowerRehearsal, optional) – With power_placebos > 0, the planned test’s power and MDE by rehearsal; None otherwise.

certificate: MAREXCertificate | None#
clusters: Dict[str, MAREXClusterDesign]#
property control_units: List[Any]#

Units assigned to control (non-zero aggregate control weight).

globres: MAREXGlobalDesign#
property mode: str#

Solver mode reported to downstream consumers.

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

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

pool: List[Dict[str, Any]] | None#
post_fit: Any | None#
power_rehearsal: MAREXPowerRehearsal | None#
recommendation: Any | None#
search: ExactSearchDiagnostics | None#
study: MAREXStudy#
property synthetic_control: ndarray#

Aggregate synthetic control outcome, shape (T,).

property synthetic_treated: ndarray#

Aggregate synthetic treated outcome, shape (T,).

property treated_units: List[Any]#

Units assigned to treatment (non-zero aggregate treated weight).

class mlsynth.utils.marex_helpers.structures.MAREXStudy(design: str, T0: int, blank_periods: int, beta: float = 1e-06, lambda1: float = 0.0, lambda2: float = 0.0, xi: float = 0.0, backend: str = 'mip')#

Bases: object

Design hyperparameters of a MAREX study (was StudyConfig).

backend records which solver produced the design: "mip" for SCIP’s mixed-integer program, "exact" for the search over treated sets in mlsynth.utils.marex_helpers.exact.

T0: int#
backend: str = 'mip'#
beta: float = 1e-06#
blank_periods: int#
design: str#
lambda1: float = 0.0#
lambda2: float = 0.0#
xi: float = 0.0#
class mlsynth.utils.marex_helpers.structures.MAREXUnitInference(unit: str, cluster: str, control_weights: Dict[str, float], synthetic_control: ndarray, inference: MAREXInference, p_bonferroni: float, blank_check: MAREXBlankCheck)#

Bases: object

Inference for one treated unit against its own synthetic control.

Parameters:

  • unit (str) – The treated unit’s label.

  • cluster (str) – Its cluster’s label.

  • control_weights (dict) – Its synthetic control, Abadie and Zhao’s v_.j: simplex weights on its cluster’s untreated units, fitted on the fitting periods, keyed by unit label. Every untreated unit of the cluster appears, at zero where the fit puts nothing.

  • synthetic_control (np.ndarray) – sum_i v_ij y_i over the whole timeline, shape (T,).

  • inference (MAREXInference) – MAREX’s inference on the contrast y_j - synthetic_control: the Theorem 2 p-value, per-period p-values, the pointwise band and, when requested, the band for the total.

  • p_bonferroni (float) – min(1, m p), with p the unit’s Theorem 2 p-value and m the number of treated units in the design.

  • blank_check (MAREXBlankCheck) – The mean of the contrast over the blank periods, with its standard error: a check for bias in the unit’s synthetic control.

blank_check: MAREXBlankCheck#
cluster: str#
control_weights: Dict[str, float]#
inference: MAREXInference#
p_bonferroni: float#
synthetic_control: ndarray#
unit: str#
class mlsynth.utils.marex_helpers.structures.MAREXUnitLevel(units: Dict[str, MAREXUnitInference], joint_p_value: float, eq11: MAREXInference, eq11_control: ndarray)#

Bases: object

Inference for every treated unit, and Abadie and Zhao’s eq. (11).

Parameters:

  • units (dict) – MAREXUnitInference by treated unit label.

  • joint_p_value (float) – A permutation test of no effect on any treated unit. Its statistic is the largest, over units, of the mean absolute post-period contrast divided by the unit’s mean absolute contrast over the blank and post periods; one arrangement of periods is applied to every unit, so the units’ common factors stay aligned, and the scale is one no arrangement changes. Exact when the arrangements number at most max_combinations, (1 + k) / (1 + B) otherwise.

  • eq11 (MAREXInference) – Inference on eq. (11)’s aggregate, the synthetic treated minus eq11_control. Under cumulative_method="inversion" its band for the total is the shared contrast’s: the inversion refits the synthetic treated series on the untreated units and never reads a control, so only the point, eq. (11)’s own total, differs.

  • eq11_control (np.ndarray) – sum_j w_j sum_i v_ij y_i with w the aggregate treated weights, shape (T,): each treated unit’s own control, weighted as the design weights the unit.

eq11: MAREXInference#
eq11_control: ndarray#
joint_p_value: float#
units: Dict[str, MAREXUnitInference]#

In addition, MAREX.fit() attaches a SyntheticControlPostFit as results.post_fit: the standardized diagnostics container shared across the MAREX family (LEXSCM, MAREX, SYNDES, PANGEO). It carries the ATE / total / lift / per-period / cumulative effect summaries, the inference triple (\(p\), CI), the pre-/blank-/post-period RMSEs, the three standardized-mean-difference blocks (treated-vs-control, treated-vs-population, control-vs-population), and — when a valid noise window exists — a PowerAnalysis block with the headline MDE and the MDE-versus-horizon curve.

class mlsynth.utils.post_fit.SyntheticControlPostFit(treated_series: ndarray, control_series: ndarray, gap_series: ndarray, n_fit: int, n_blank: int, n_post: int, ate: float | None = None, total_effect: float | None = None, ate_percent: float | None = None, ate_per_period: ndarray | None = None, cumulative_effect: ndarray | None = None, p_value: float | None = None, ci_lower: float | None = None, ci_upper: float | None = None, inference_method: str | None = None, rmse_fit: float | None = None, rmse_blank: float | None = None, rmse_post: float | None = None, covariate_names: Tuple[str, ...] = (), covariate_smd: Dict[str, float] | None = None, covariate_smd_abs_max: float | None = None, covariate_smd_squared_sum: float | None = None, covariate_smd_treated_vs_pop: Dict[str, float] | None = None, covariate_smd_treated_vs_pop_abs_max: float | None = None, covariate_smd_treated_vs_pop_squared_sum: float | None = None, covariate_smd_control_vs_pop: Dict[str, float] | None = None, covariate_smd_control_vs_pop_abs_max: float | None = None, covariate_smd_control_vs_pop_squared_sum: float | None = None, power: PowerAnalysis | None = None)#

Bases: object

Standardized post-fit diagnostics for a single synthetic control design.

Field semantics are estimator-agnostic; every MAREX-family adapter populates the same shape. Any field that isn’t naturally computable for the producing estimator is left None.

ate: float | None = None#
ate_per_period: ndarray | None = None#
ate_percent: float | None = None#
ci_lower: float | None = None#
ci_upper: float | None = None#
control_series: ndarray#
covariate_names: Tuple[str, ...] = ()#
covariate_smd: Dict[str, float] | None = None#
covariate_smd_abs_max: float | None = None#
covariate_smd_control_vs_pop: Dict[str, float] | None = None#
covariate_smd_control_vs_pop_abs_max: float | None = None#
covariate_smd_control_vs_pop_squared_sum: float | None = None#
covariate_smd_squared_sum: float | None = None#
covariate_smd_treated_vs_pop: Dict[str, float] | None = None#
covariate_smd_treated_vs_pop_abs_max: float | None = None#
covariate_smd_treated_vs_pop_squared_sum: float | None = None#
cumulative_effect: ndarray | None = None#
gap_series: ndarray#
inference_method: str | None = None#
n_blank: int#
n_fit: int#
n_post: int#
p_value: float | None = None#
power: PowerAnalysis | None = None#
rmse_blank: float | None = None#
rmse_fit: float | None = None#
rmse_post: float | None = None#
total_effect: float | None = None#
treated_series: ndarray#
class mlsynth.utils.post_fit.PowerAnalysis(headline: MDEPoint, curve: Tuple[MDEPoint, ...], alpha: float, power_target: float, sigma_placebo: float, serial_correlation: float, baseline: float, method: str = 'analytical_ar1_mean_gap', placebo_bias: float = 0.0, placebo_bias_pvalue: float = nan, n_placebo: int = 0, sigma_ci: Tuple[float, float] = (nan, nan))#

Bases: object

Standardized power-analysis output attached to SyntheticControlPostFit.

Built from the placebo / blank-period gap variance and an analytical Gaussian approximation, with AR(1) variance inflation to handle serial correlation in the gap residuals. The intent matches the per-estimator power modules already in the library (PangeoPower, SPCDPowerAnalysis, SYNDESPower) but consumes the same SyntheticControlPostFit shape so every covariate-aware SCM-family estimator gets the surface for free.

headline#

MDE for the actual n_post horizon of the realised design.

Type:

MDEPoint

curve#

MDE / power values across the requested post_grid horizons (so callers can read a detectability curve).

Type:

list of MDEPoint

alpha#

Two-sided significance level assumed.

Type:

float

power_target#

Target power the MDEs are computed at (default 0.80).

Type:

float

sigma_placebo#

Standard deviation of the placebo gap series, about that window’s own mean, used as the noise scale. The offset itself is reported separately as placebo_bias and enters the standard error there.

Type:

float

serial_correlation#

Lag-1 (AR(1)) autocorrelation of the placebo gap residuals, estimated on the demeaned series and small-sample corrected, used to inflate the variance for serial dependence.

Type:

float

placebo_bias#

Mean of the placebo gap. Zero for a design that is unbiased on the window it was not fitted to; otherwise the offset the MDE must clear alongside the noise.

Type:

float

placebo_bias_pvalue#

Two-sided p-value for placebo_bias == 0.

Type:

float

n_placebo#

Periods sigma_placebo is estimated from.

Type:

int

sigma_ci#

Chi-square confidence interval for sigma_placebo at alpha, on the placebo window’s effective degrees of freedom.

Type:

tuple of float

baseline#

Mean of the control trajectory on the post window (denominator for mde_pct). NaN when no post window exists.

Type:

float

method#

"analytical_ar1_mean_gap": the closed-form MDE for a two-sided test on the mean post-period gap, with AR(1) variance inflation and a Student-t critical value. This is a different test from the moving-block permutation on mean|e| that LEXSCM’s design search ranks candidates by (Vives-i-Bastida 2022, Section 2), so the two MDEs need not agree. Reserved for future "monte_carlo" extensions.

Type:

str

alpha: float#
baseline: float#
curve: Tuple[MDEPoint, ...]#
headline: MDEPoint#
mde_by_horizon() → Dict[int, float]#

{post_periods: mde_pct} for quick lookup.

method: str = 'analytical_ar1_mean_gap'#
n_placebo: int = 0#
placebo_bias: float = 0.0#
placebo_bias_pvalue: float = nan#
power_target: float#
serial_correlation: float#
sigma_ci: Tuple[float, float] = (nan, nan)#
sigma_placebo: float#
class mlsynth.utils.post_fit.MDEPoint(post_periods: int, mde_absolute: float, mde_pct: float, se: float, power_at_observed: float | None = None, critical_value: float = nan, mde_ci: Tuple[float, float] = (nan, nan))#

Bases: object

Minimum detectable effect at a single post-treatment horizon.

critical_value: float = nan#
mde_absolute: float#
mde_ci: Tuple[float, float] = (nan, nan)#
mde_pct: float#
post_periods: int#
power_at_observed: float | None = None#
se: float#

Helper Modules#

Input preparation for MAREX: long panel -> design-ready arrays.

class mlsynth.utils.marex_helpers.setup.MAREXPanel(Y_full: DataFrame, clusters: ndarray, T0: int, blank_periods: int, unit_index: IndexSet, time_index: IndexSet, covariates: ndarray | None = None, covariate_names: Tuple[str, ...] = ())#

Prepared MAREX inputs.

unit_index / time_index are the single source of truth for unit and time identity (label <-> integer index). Y_full rows align to unit_index.labels and its columns to time_index.labels; the clusters vector is in unit_index order. Every downstream consumer (optimizer, orchestrator, restriction builder) indexes through these IndexSets rather than re-deriving identity from a DataFrame.

T0: int#
Y_full: DataFrame#
blank_periods: int#
clusters: ndarray#
covariate_names: Tuple[str, ...] = ()#
covariates: ndarray | None = None#
time_index: IndexSet#
unit_index: IndexSet#
mlsynth.utils.marex_helpers.setup.prepare_marex_panel(df: DataFrame, outcome: str, unitid: str, time: str, cluster: str | None, T0: int | None, inference: bool, blank_periods: int, T_post: int | None, covariates: List[str] | None = None) → MAREXPanel#

Pivot the long panel to units x time and forward the resolved T0 / blank_periods.

The MAREX config validator is the single source of truth for resolving T0 (from either an explicit scalar or a post_col 0/1 column) and the default 30%-of-pre-tail blank window — by the time this helper runs both are concrete integers. covariates columns are aggregated to a per-unit mean over the fitting window, the first T0 - blank_periods periods, via mlsynth.utils.datautils.build_covariate_matrix() and returned as an (N, R) matrix aligned to the unit order. The matrix is left un-normalised here so MAREX’s existing standardize=True flag (applied to the combined [Y_fit; covariate_weight * Z] predictor matrix in marex_helpers.optimization) keeps its previous behaviour.

Design-formulation primitives for MAREX (Abadie & Zhao 2026).

The experimenter chooses treated weights w and control weights v per cluster on the simplex, with a binary selection mask z linking them (w_j <= z_j, v_j <= 1 - z_j) so a unit is either treated or a control, never both (the disjointness w_j v_j = 0). These helpers build the cvxpy variables, constraints, and the design-specific objective; the objective form is selected by design:

  • "standard" – match each cluster mean with both synthetic units;

  • "weakly_targeted" – match the treated synthetic to the mean and softly tie the control synthetic to it (weight beta);

  • "penalized" – standard plus cluster-level distance penalties (lambda1 / lambda2);

  • "unit_penalized" – standard plus unit-level penalties (xi / lambda1_unit).

Two further penalties appear in the paper’s text but have no convex form in these variables, and are refused with a reason: lambda2_unit pairs w against v through a distance matrix (bilinear in two decision variables), and zeta minimises z(1 - z) (concave). The relaxed path reaches integrality by discretising after the solve, not by penalising during it.

xi is the eq. (10) unit-level term. It weights each candidate’s own reproducibility by that candidate’s treated weight, so the inner quantity is a minimum over the per-unit control weights and is therefore constant in w: per_unit_reproducibility() evaluates those minima once and the penalty enters linearly. Multiplying w by a convex function of the shared joint v instead is not the same objective and is not convex in (w, v).

mlsynth.utils.marex_helpers.formulation.WEIGHT_FLOOR = 1e-08#

A treated market carries weight above this or it is not treated. Abadie and Zhao bound ||w||_0, the count of nonzero treated weights, and MAREX’s results call a market treated when its weight is above 1e-8.

mlsynth.utils.marex_helpers.formulation.build_constraints(w, v, z, M, cluster_members, cluster_labels, m_eq, m_min, m_max, costs, budget_dict, exclusive, restrictions=None, selection_budget=None)#

Simplex, disjointness, cardinality, cost, exclusivity, and (optional) geographic-restriction constraints.

selection_budget is (c, B), one cost per row and a budget: the treated markets’ costs, each counted once across clusters, sum to at most B (sum_j c_j sum_k z[j,k] <= B).

mlsynth.utils.marex_helpers.formulation.build_membership_mask(clusters, label_to_k, N, K)#

Boolean (N, K) mask of unit-to-cluster membership.

mlsynth.utils.marex_helpers.formulation.build_objective(Y_fit, Xbar_clusters, cluster_members, w, v, z, design, beta=1e-06, lambda1=0.0, lambda2=0.0, xi=0.0, lambda1_unit=0.0, lambda2_unit=0.0, D1=None, D2_list=None, zeta=0.0, control_penalty=0.0)#

Design-specific cvxpy objective (see module docstring).

lambda2_unit and a non-zero zeta raise MlsynthConfigError: neither has a convex form in these variables. The relaxed (continuous-z) solve discretises after the fact via post_hoc_discretize.

mlsynth.utils.marex_helpers.formulation.carries_weight(w) → bool#

Whether every market of a treated set’s fit carries weight.

mlsynth.utils.marex_helpers.formulation.compute_cluster_means_members(Y_fit, M, cluster_labels)#

Per-cluster predictor means and member index arrays.

mlsynth.utils.marex_helpers.formulation.get_per_cluster_param(param, klabel, default=None)#

Resolve a possibly-per-cluster parameter to its value for klabel.

mlsynth.utils.marex_helpers.formulation.init_cvxpy_variables(N, K, boolean=True)#

Treated (w), control (v) weights and selection (z).

z is binary for the exact MIQP (boolean=True) or continuous in [0, 1] for the relaxed QP (boolean=False).

mlsynth.utils.marex_helpers.formulation.per_unit_reproducibility(Y_fit, members)#

min_v ||x_j - sum_i v_i x_i||^2 for each unit, over the simplex.

Equation (10)’s penalty is xi * sum_j w_j * min_{v_.j} ||x_j - sum_i v_ij x_i||^2. The inner weights appear in exactly one term of the objective, multiplied by the non-negative scalar w_j, and a positive scalar cannot move an argmin – so v*_.j does not depend on w and the minimum is a constant. That is what makes the term linear in w, and therefore something a convex program can carry.

Written instead as w_j times an evaluated squared norm of a joint variable, the term is a product of two variables: the Hessian of w v^2 is [[0, 2v], [2v, 2w]] with determinant -4 v^2, so it is genuinely non-convex and no solver will take it.

Each unit is fitted against every other unit in the panel. Equation (10) excludes the treated set from one another’s donor pools, which cannot be imposed here because the assignment is itself a decision variable; the difference is the m - 1 donors a chosen tuple would remove, measured elsewhere at a median 0.0 to 0.7 per cent in the penalty with the same design selected in 11 of 12 cells.

Only members are fitted; every other entry of the returned length-N array is nan. A unit’s charge depends on the panel, not on the cluster asking for it, and each cluster asks for its own members, so a design with K clusters fits each unit once on the exact backend and twice on the MIQP path, not K and 2K times.

mlsynth.utils.marex_helpers.formulation.precompute_distances(Y_fit, Xbar_clusters, cluster_members)#

Unit-to-cluster-mean distances D1 and within-cluster pairwise D2.

mlsynth.utils.marex_helpers.formulation.prepare_clusters(Y_full, clusters)#

Coerce Y_full/clusters to arrays and return cluster bookkeeping.

mlsynth.utils.marex_helpers.formulation.prepare_fit_slices(Y_full_np, T0, blank_periods)#

Split the pre-period into a fitting slice and a held-out blank slice.

mlsynth.utils.marex_helpers.formulation.validate_costs_budget(costs, budget, N, cluster_labels, K)#

Validate cost/budget inputs; return (costs_np, budget_dict).

mlsynth.utils.marex_helpers.formulation.validate_scm_inputs(Y_full, T0, blank_periods, design, beta=1e-06, lambda1=0.0, lambda2=0.0, xi=0.0, lambda1_unit=0.0, lambda2_unit=0.0, control_penalty=0.0)#

Validate shapes and design/parameter compatibility (raises ValueError).

MAREX design optimizers (Abadie & Zhao 2026).

solve_design solves the exact mixed-integer design (binary selection z); solve_design_relaxed relaxes z to [0, 1], solves the QP, then discretizes post hoc. Both return a raw result dict consumed by the orchestrator.

SCIP accepts a binary within its integrality tolerance of an integer, so a selection left near 1e-6 rounds to “not selected” while still admitting a treated weight of that size, and a market can carry weight on both sides of the design. solve_design therefore takes only the selection from SCIP: it fixes z at its rounding and refits the weights by refit_selection() with the selection held constant, so the weights and objective reported are those of the design the selection defines.

mlsynth.utils.marex_helpers.optimization.MAX_WEIGHT_CUTS = 100#

Re-solves solve_design() makes to cut treated sets with an unweighted market before it reports the design infeasible.

mlsynth.utils.marex_helpers.optimization.finish_design(Y_full, Y_full_np, w_opt, v_opt, z_opt, Xbar_clusters, cluster_labels, cluster_members, N, K, T0, blank_periods, T_fit, Y_fit, Y_blank, design, clusters, objective)#

The design dict both backends return, from their solved weights.

solve_design (the MIQP) and mlsynth.utils.marex_helpers.exact.solve_design_exact() (the search) reach w, v and z differently and hand them here, so everything downstream – inference, post-fit diagnostics, the result object and the solution pool – reads one structure whichever backend produced it.

mlsynth.utils.marex_helpers.optimization.post_hoc_discretize(w_opt, v_opt, cluster_members, cluster_labels, m_eq=None, m_min=None, m_max=None, trim_threshold=0.01, Y_fit=None, Y_blank=None)#

Round relaxed weights to a feasible integer design (was internal).

mlsynth.utils.marex_helpers.optimization.refit_selection(X_fit, Xbar_clusters, cluster_members, cluster_labels, M, z_sel, *, m_eq=None, m_min=None, m_max=None, costs=None, budget_dict=None, exclusive=True, restrictions=None, design='standard', beta=1e-06, lambda1=0.0, lambda2=0.0, xi=0.0, lambda1_unit=0.0, lambda2_unit=0.0, D1=None, D2_list=None, solver='SCIP', control_penalty=0.0)#

The weights of a fixed selection z_sel, and their objective.

z_sel is rounded at 0.5 before anything else reads it.

For the standard, penalized and unit_penalized designs with no cost budget and exclusive arms, fixing the selection splits each cluster into a treated fit and a control fit, and both are solved exactly by the active-set simplex solver: the computation the exact backend’s search runs on each candidate, so a weight the optimum puts at zero is exactly zero. It reads the exact backend’s own penalty and fit helpers, so comparing the two backends cannot check those helpers; an oracle built from the definitions does (tests/test_marex_exact_oracle.py, #783). Otherwise the MIQP’s own constraints and objective are rebuilt with z continuous and pinned to z_sel, a convex QP, and handed to solver, the solver the MIQP used.

Raises:

MlsynthEstimationError – When the selection breaks the cardinality or a restriction, leaves a cluster with no control market, or the QP is not solved to optimality.

mlsynth.utils.marex_helpers.optimization.resolve_covariate_weight(covariate_weight, covariate_names)#

The covariate weight as _augment_fit takes it.

A scalar weights every covariate alike and is returned as given. A mapping from covariate name to weight, the diagonal H of Vives-i-Bastida (2022, program 1), becomes a vector in the order of covariate_names, so it is matched by name and the order the covariates were listed in does not matter.

Raises:

MlsynthConfigError – When a mapping does not name exactly the covariates.

mlsynth.utils.marex_helpers.optimization.solve_design(Y_full, T0, clusters, *, forbidden=None, **kwargs)#

Exact mixed-integer MAREX design (was SCMEXP), every treated market weighted.

The MIQP links a market’s weight to its selection one way, w_j <= z_j: an unselected market carries no weight, but a selected one may carry none either, so sum_j z_j = m counts markets, not the nonzero treated weights Abadie and Zhao bound. Their program, ||w||_0 = m, has no minimum where the best fit wants a market at zero weight: the value of giving it weight eps falls toward the zero-weight design’s as eps falls, and never reaches it. A treated set is therefore a design only when its own fit weights every market in it, the rule the exact backend applies.

After each solve the selection is refit (refit_selection()); where a cluster’s treated set has a market at weight <= WEIGHT_FLOOR, that set is cut from that cluster, sum z[S] <= |S| - 1, and the program solved again. Arguments are _solve_design_once()’s.

mlsynth.utils.marex_helpers.optimization.solve_design_pool(Y_full, T0, clusters, *, top_K=1, **kwargs)#

Enumerate up to top_K distinct exact designs via no-good cuts.

Calls solve_design() repeatedly, forbidding each chosen assignment so the next solve returns the next-best distinct design. Stops early when the feasible region is exhausted (the solver becomes infeasible).

mlsynth.utils.marex_helpers.optimization.solve_design_relaxed(Y_full, T0, clusters, blank_periods=0, m_eq=None, m_min=None, m_max=None, exclusive=True, design='standard', beta=1e-06, lambda1=0.0, lambda2=0.0, xi=0.0, lambda1_unit=0.0, lambda2_unit=0.0, costs=None, budget=None, covariates=None, covariate_weight=1.0, standardize=False, solver=None, verbose=False, zeta=0.0, trim_threshold=0.01, control_penalty=0.0)#

Relaxed (continuous-z) design with post-hoc discretization (was SCMEXP_REL).

Placebo/permutation inference for a MAREX synthetic treated-vs-control pair.

The held-out blank pre-periods act as placebos: the synthetic treated minus synthetic control there should be noise, so its distribution calibrates a permutation p-value and a split-conformal confidence band for the post-period effect (Abadie & Zhao 2026, OA; Chernozhukov-Wuthrich-Zhu 2021).

The global test has two reference sets. "subsets" is Abadie and Zhao’s eq. (17): the mean absolute post contrast against the T_post-subsets of the blank and post periods, exact when those periods are exchangeable. "cyclic" is Chernozhukov, Wuthrich and Zhu’s moving-block test with their Remark 1 statistic: the magnitude of the post block’s signed sum against the same statistic on every cyclic shift of the blank and post contrast. The shifts keep neighbouring periods together, so it holds its level under stationary, weakly dependent contrasts, where eq. (17) over-rejects (16.5% at 5% with factors autocorrelated at 0.8; benchmarks/studies/marex_global_test_dependence/).

The band for the total effect has two constructions. "resample" draws paths from the blank-period effects alone. "inversion" also reads the post outcomes: it inverts CWZ’s test of a constant effect with their Remark 1 statistic, the magnitude of the signed post sum of the residuals, against the cyclic shifts of the residual path over the blank and post periods. The fitting periods stay out of that reference because the design was chosen on them. On simulated MAREX designs their residuals are about 30% smaller in mean square than the blank and post residuals, which behave alike, and a reference that includes them understates the spread of a 13-period post sum by a factor of 1.17 with independent factors and 1.59 with factors autocorrelated at 0.8. The study is in benchmarks/studies/marex_cumulative_conformal/.

mlsynth.utils.marex_helpers.inference.CUMULATIVE_METHODS = ('resample', 'inversion')#

How the band for the total effect is built; see the module docstring.

mlsynth.utils.marex_helpers.inference.GLOBAL_TESTS = ('subsets', 'cyclic')#

The global test’s reference set; see the module docstring.

mlsynth.utils.marex_helpers.inference.compute_inference(Y_treated: ndarray, Y_control: ndarray, T0: int, TcE: int, Tb: int, alpha: float = 0.05, max_combinations: int = 1000, random_state: int | None = None, cumulative_band: bool = False, cumulative_block: int = 0, cumulative_n_sim: int = 2000, cumulative_seed: int = 0, cumulative_method: str = 'resample', control_outcomes: ndarray | None = None, T_post: int | None = None, global_test: str = 'subsets') → MAREXInference#

Compute permutation inference for one synthetic contrast.

Parameters:

  • Y_treated, Y_control (np.ndarray) – Synthetic treated / control outcomes over the full timeline, shape (T,).

  • T0 (int) – Number of pre-treatment periods.

  • TcE (int) – Start index of the blank (held-out) window.

  • Tb (int) – Number of blank periods.

  • alpha (float) – Two-sided significance level.

  • max_combinations (int) – Number of permutation draws for the global test.

  • random_state (int, optional) – Seed for reproducibility.

  • cumulative_band (bool) – Whether to add a band for the total post-period effect.

  • cumulative_block (int) – "resample" only: block length for the draw; 0 means the whole horizon.

  • cumulative_n_sim (int) – "resample" only: number of accumulated paths to draw.

  • cumulative_seed (int) – "resample" only: seed for the draw, held separate from random_state so the permutation test and the band are reproducible independently.

  • cumulative_method ({“resample”, “inversion”}) – How the band for the total is built; see the module docstring.

  • control_outcomes (np.ndarray, optional) – "inversion" only: outcomes of the units the refit may draw on, shape (J, T), aligned to Y_treated. MAREX passes the units the design left untreated.

  • T_post (int, optional) – Post periods the inference reads; every one after T0 by default.

  • global_test ({“subsets”, “cyclic”}) – The global test’s reference set and statistic; see the module docstring. max_combinations and random_state apply to "subsets" alone.

Returns:

MAREXInference

Raises:

MlsynthConfigError – On an unknown cumulative_method or global_test, or a band by "inversion" without control_outcomes.

Warns:

UserWarning – When the inversion’s reference is too short for any p-value to reach alpha, so the band is unbounded on both sides.

Top-level MAREX solve: run the design optimizer and assemble frozen results.

mlsynth.utils.marex_helpers.orchestration.solve_marex(Y_full, T0, clusters, design='standard', blank_periods=0, m_eq=None, m_min=None, m_max=None, exclusive=True, beta=1e-06, lambda1=0.0, lambda2=0.0, xi=0.0, lambda1_unit=0.0, lambda2_unit=0.0, costs=None, budget=None, covariates=None, covariate_names=(), covariate_weight=1.0, standardize=False, solver=None, verbose=False, relaxed=False, inference=False, alpha=0.05, T_post=None, unit_inference=False, max_combinations=1000, random_state=42, unit_index=None, time_index=None, restrictions=None, warm_start=None, time_limit=None, cumulative_band=False, cumulative_method='resample', cumulative_block=0, cumulative_n_sim=2000, cumulative_seed=0, global_test='subsets', backend='mip', control_penalty=0.0, certify=False, market_costs=None, max_total_cost=None) → MAREXResults#

Solve the MAREX design and return a frozen MAREXResults.

With relaxed=True the continuous-z QP with post-hoc discretization is used; otherwise the exact MIQP. With inference=True, blank-period placebo inference is computed for every cluster and the aggregate, at level alpha and over the first T_post periods after T0 (all of them when None). unit_inference=True adds it for each treated unit against its own synthetic control, with eq. (11)’s aggregate (mlsynth.utils.marex_helpers.unit_inference).

backend="exact" solves the design by searching treated sets (mlsynth.utils.marex_helpers.exact.solve_design_exact()) instead of handing the MIQP to SCIP; it returns the same design dict, so everything after the solve is shared. With certify=True it also certifies the design’s optimality gap (mlsynth.utils.marex_helpers.certify).

market_costs (one per row of Y_full) with max_total_cost caps the total cost of the treated markets, sum_j c_j z_j <= B, on either backend.

cumulative_band=True adds a band for the total post-period effect to that inference, block-resampling the blank-period effects with block length cumulative_block (0 means the whole horizon) into cumulative_n_sim accumulated paths, seeded by cumulative_seed. It needs no extra solve. global_test chooses the global test of every contrast: "subsets", the paper’s eq. (17), or "cyclic", the cyclic shifts of the blank and post contrast (mlsynth.utils.marex_helpers.inference).

cumulative_method="inversion" builds that band by test inversion instead, refitting each synthetic treated series on the units its design left untreated: a cluster’s own unselected members, or every untreated unit for the aggregate.

Linear-factor DGP for the MAREX simulation study (Abadie & Zhao 2026, Sec. 5).

Reimplements the paper’s baseline data-generating process (Assumption 1, equations 12a/12b) so the simulation study can be replicated without the authors’ code. Potential outcomes are

Y^N_jt = delta_t + theta_t’ Z_j + lambda_t’ mu_j + eps_jt Y^I_jt = upsilon_t + gamma_t’ Z_j + eta_t’ mu_j + xi_jt

with Z_j (R observed) and mu_j (F unobserved) covariates, sorted time effects, and i.i.d. Normal(0, sigma^2) noise.

generate_marex_sample() draws from a numpy generator. generate_section5_sample() reproduces the authors’ own simulations draw for draw: simulation i of their 1000 is set.seed(i) followed by their Generate_Model_Primitives (SCDesign, “2. Many Simulations”, Main_LazyRun.R), and the function returns the same panel from R’s stream, ported here.

class mlsynth.utils.marex_helpers.simulation.MAREXSample(Y_N: ndarray, Y_I: ndarray, tau: ndarray, T0: int, Z: ndarray, mu: ndarray | None = None)#

One simulated sample.

Y_N#

Potential outcomes under no treatment, shape (J, T).

Type:

np.ndarray

Y_I#

Potential outcomes under treatment, shape (J, T).

Type:

np.ndarray

tau#

True average treatment effect tau_t per period, shape (T,) (zero in the pre-period).

Type:

np.ndarray

T0#

Number of pre-treatment periods.

Type:

int

Z#

Observed (time-invariant) covariates that generate the outcome, shape (J, R).

Type:

np.ndarray

mu#

Unobserved (time-invariant) covariates, shape (J, F). A design never sees them; a simulation that knows them can compute the bias a design carries.

Type:

np.ndarray, optional

T0: int#
Y_I: ndarray#
Y_N: ndarray#
Z: ndarray#
mu: ndarray | None = None#
tau: ndarray#
mlsynth.utils.marex_helpers.simulation.generate_marex_sample(J: int = 15, R: int = 7, F: int = 11, T: int = 30, T0: int = 25, sigma: float = 1.0, rng: Generator | None = None) → MAREXSample#

Draw one sample from the paper’s baseline linear-factor DGP (Sec. 5).

Returns:

MAREXSample

mlsynth.utils.marex_helpers.simulation.generate_section5_sample(seed: int) → MAREXSample#

Simulation seed of Abadie and Zhao’s Section 5, as their R code draws it.

The panel is what set.seed(seed) followed by the authors’ Generate_Model_Primitives produces: J = 15 units, R = 7 observed and F = 11 unobserved covariates on Uniform(0, 1), T = 30 periods with T0 = 25, coefficients on Uniform(0, 10), intercepts on Uniform(0, 20) sorted small to large (delta over all 30 periods together, upsilon over the five experimental ones), and N(0, 1) noise. Every coefficient vector is drawn afresh each period. The treated outcomes exist only on the experimental periods, which is where the authors define them; Y_I is NaN before T0.

Their driver runs simulation i of 1000 as set.seed(i), so seeds 1 to 1000 give the panels behind their Table 2.

Parameters:

seed (int) – From 0 to 2^31 - 1.

Returns:

MAREXSample – With mu, the unobserved covariates, filled.

Raises:

MlsynthConfigError – If seed is not an integer in R’s range.

Plotting for MAREX: synthetic treated vs control (or the treatment effect).

mlsynth.utils.marex_helpers.plotter.plot_marex(results: MAREXResults, clusters: List[str] | None = None, plot_type: str = 'treatment', global_result: bool = True, figsize: tuple = (12, 6), donor_cloud: bool = False) → None#

Plot MAREX treatment effects (or predictions), one panel per cluster + global.

Parameters:

  • results (MAREXResults) – Output of mlsynth.estimators.MAREX.

  • clusters (list of str, optional) – Cluster labels to plot (default: all; a lone "0" cluster is skipped).

  • plot_type ({“treatment”, “prediction”}) – Plot the treated-minus-control effect, or both synthetic series.

  • global_result (bool) – Include the aggregate (global) panel.

  • donor_cloud (bool) – On the global "prediction" panel, overlay one faint line per unit (the rows of results.globres.Y_full) behind the series – the “observed data” cloud of Abadie & Zhao’s Figure 4. Ignored on the "treatment" (effect) plot and on cluster panels, which carry no per-unit outcome matrix. The population mean (a thick solid line) is drawn whenever the per-unit matrix is available, with or without the cloud; the synthetic treated and control are thinner dashed red / blue lines so they stand out against it.

The shared post-fit module — compute_smd(), compute_post_fit(), and compute_power_analysis() — lives outside the marex_helpers package so the other MAREX-family estimators (LEXSCM, SYNDES, PANGEO) can call into the same one-source-of-truth diagnostics:

Standardized post-fit diagnostics for synthetic control designs and the matching power-analysis surface that consumes them.

After any MAREX-family estimator (LEXSCM, MAREX, SYNDES, PANGEO, …) solves its design problem, downstream consumers (the SAGE dashboard, paper-style reports, comparison tables) all need the same numbers:

  • the post-treatment ATT, total effect, percentage lift, per-period gap;

  • pre / blank / post root-mean-squared-error of the synthetic gap;

  • inference scalars (p-value, CI bounds) when computed;

  • covariate-balance standardized mean differences (SMDs) when covariates were used in the design.

This module exposes one frozen dataclass (SyntheticControlPostFit) and three free functions:

  • compute_smd() – standalone, panel-independent SMD

    from any (cov_matrix, treated_w, control_w);

  • compute_post_fit() – the full diagnostic bundle from

    trajectories + boundaries + (optional) covariate matrix + (optional) inference;

  • compute_post_fit_marex() – adapter that builds the bundle from a

    MAREXResults + MAREXPanel pair.

The free-function entry points are deliberately small and reusable, so the LEXSCM / SYNDES / PANGEO equivalents can be added one-at-a-time without touching this module: they just compose the same primitives.

class mlsynth.utils.post_fit.MDEPoint(post_periods: int, mde_absolute: float, mde_pct: float, se: float, power_at_observed: float | None = None, critical_value: float = nan, mde_ci: Tuple[float, float] = (nan, nan))#

Bases: object

Minimum detectable effect at a single post-treatment horizon.

critical_value: float = nan#
mde_absolute: float#
mde_ci: Tuple[float, float] = (nan, nan)#
mde_pct: float#
post_periods: int#
power_at_observed: float | None = None#
se: float#
class mlsynth.utils.post_fit.PowerAnalysis(headline: MDEPoint, curve: Tuple[MDEPoint, ...], alpha: float, power_target: float, sigma_placebo: float, serial_correlation: float, baseline: float, method: str = 'analytical_ar1_mean_gap', placebo_bias: float = 0.0, placebo_bias_pvalue: float = nan, n_placebo: int = 0, sigma_ci: Tuple[float, float] = (nan, nan))#

Bases: object

Standardized power-analysis output attached to SyntheticControlPostFit.

Built from the placebo / blank-period gap variance and an analytical Gaussian approximation, with AR(1) variance inflation to handle serial correlation in the gap residuals. The intent matches the per-estimator power modules already in the library (PangeoPower, SPCDPowerAnalysis, SYNDESPower) but consumes the same SyntheticControlPostFit shape so every covariate-aware SCM-family estimator gets the surface for free.

headline#

MDE for the actual n_post horizon of the realised design.

Type:

MDEPoint

curve#

MDE / power values across the requested post_grid horizons (so callers can read a detectability curve).

Type:

list of MDEPoint

alpha#

Two-sided significance level assumed.

Type:

float

power_target#

Target power the MDEs are computed at (default 0.80).

Type:

float

sigma_placebo#

Standard deviation of the placebo gap series, about that window’s own mean, used as the noise scale. The offset itself is reported separately as placebo_bias and enters the standard error there.

Type:

float

serial_correlation#

Lag-1 (AR(1)) autocorrelation of the placebo gap residuals, estimated on the demeaned series and small-sample corrected, used to inflate the variance for serial dependence.

Type:

float

placebo_bias#

Mean of the placebo gap. Zero for a design that is unbiased on the window it was not fitted to; otherwise the offset the MDE must clear alongside the noise.

Type:

float

placebo_bias_pvalue#

Two-sided p-value for placebo_bias == 0.

Type:

float

n_placebo#

Periods sigma_placebo is estimated from.

Type:

int

sigma_ci#

Chi-square confidence interval for sigma_placebo at alpha, on the placebo window’s effective degrees of freedom.

Type:

tuple of float

baseline#

Mean of the control trajectory on the post window (denominator for mde_pct). NaN when no post window exists.

Type:

float

method#

"analytical_ar1_mean_gap": the closed-form MDE for a two-sided test on the mean post-period gap, with AR(1) variance inflation and a Student-t critical value. This is a different test from the moving-block permutation on mean|e| that LEXSCM’s design search ranks candidates by (Vives-i-Bastida 2022, Section 2), so the two MDEs need not agree. Reserved for future "monte_carlo" extensions.

Type:

str

alpha: float#
baseline: float#
curve: Tuple[MDEPoint, ...]#
headline: MDEPoint#
mde_by_horizon() → Dict[int, float]#

{post_periods: mde_pct} for quick lookup.

method: str = 'analytical_ar1_mean_gap'#
n_placebo: int = 0#
placebo_bias: float = 0.0#
placebo_bias_pvalue: float = nan#
power_target: float#
serial_correlation: float#
sigma_ci: Tuple[float, float] = (nan, nan)#
sigma_placebo: float#
class mlsynth.utils.post_fit.SyntheticControlPostFit(treated_series: ndarray, control_series: ndarray, gap_series: ndarray, n_fit: int, n_blank: int, n_post: int, ate: float | None = None, total_effect: float | None = None, ate_percent: float | None = None, ate_per_period: ndarray | None = None, cumulative_effect: ndarray | None = None, p_value: float | None = None, ci_lower: float | None = None, ci_upper: float | None = None, inference_method: str | None = None, rmse_fit: float | None = None, rmse_blank: float | None = None, rmse_post: float | None = None, covariate_names: Tuple[str, ...] = (), covariate_smd: Dict[str, float] | None = None, covariate_smd_abs_max: float | None = None, covariate_smd_squared_sum: float | None = None, covariate_smd_treated_vs_pop: Dict[str, float] | None = None, covariate_smd_treated_vs_pop_abs_max: float | None = None, covariate_smd_treated_vs_pop_squared_sum: float | None = None, covariate_smd_control_vs_pop: Dict[str, float] | None = None, covariate_smd_control_vs_pop_abs_max: float | None = None, covariate_smd_control_vs_pop_squared_sum: float | None = None, power: PowerAnalysis | None = None)#

Bases: object

Standardized post-fit diagnostics for a single synthetic control design.

Field semantics are estimator-agnostic; every MAREX-family adapter populates the same shape. Any field that isn’t naturally computable for the producing estimator is left None.

ate: float | None = None#
ate_per_period: ndarray | None = None#
ate_percent: float | None = None#
ci_lower: float | None = None#
ci_upper: float | None = None#
control_series: ndarray#
covariate_names: Tuple[str, ...] = ()#
covariate_smd: Dict[str, float] | None = None#
covariate_smd_abs_max: float | None = None#
covariate_smd_control_vs_pop: Dict[str, float] | None = None#
covariate_smd_control_vs_pop_abs_max: float | None = None#
covariate_smd_control_vs_pop_squared_sum: float | None = None#
covariate_smd_squared_sum: float | None = None#
covariate_smd_treated_vs_pop: Dict[str, float] | None = None#
covariate_smd_treated_vs_pop_abs_max: float | None = None#
covariate_smd_treated_vs_pop_squared_sum: float | None = None#
cumulative_effect: ndarray | None = None#
gap_series: ndarray#
inference_method: str | None = None#
n_blank: int#
n_fit: int#
n_post: int#
p_value: float | None = None#
power: PowerAnalysis | None = None#
rmse_blank: float | None = None#
rmse_fit: float | None = None#
rmse_post: float | None = None#
total_effect: float | None = None#
treated_series: ndarray#
mlsynth.utils.post_fit.compute_post_fit(treated_series: ndarray, control_series: ndarray, *, n_fit: int, n_blank: int = 0, n_post: int | None = None, cov_matrix: ndarray | None = None, cov_names: Sequence[str] | None = None, cov_scales: ndarray | None = None, treated_weights: ndarray | None = None, control_weights: ndarray | None = None, population_weights: ndarray | None = None, inference: Any | None = None, n_treated_units: int | None = None) → SyntheticControlPostFit#

Compute a SyntheticControlPostFit from trajectories + boundaries.

The trajectories treated_series and control_series are the estimator’s own synthetic constructs (Σⱼ wⱼ Yⱼ and Σⱼ vⱼ Yⱼ in Abadie-Zhou notation). n_post defaults to len(treated_series) - n_fit - n_blank.

Covariate balance fields are populated when cov_matrix + treated_weights + control_weights are all supplied (the natural inputs for any MAREX-family design). The compute_smd() helper does the work, so the SMD numbers are exactly consistent with a standalone call to compute_smd().

Inference scalars are pulled from the estimator’s inference object via _extract_inference(), which knows about the four common shapes (LEXSCM Inference, MAREX MAREXInference, SYNDES SYNDESInference, or a plain dict). All inference fields are optional.

mlsynth.utils.post_fit.compute_post_fit_marex(raw, panel, *, cov_scales: ndarray | None = None) → SyntheticControlPostFit#

Adapt a MAREXResults + MAREXPanel pair into a SyntheticControlPostFit.

Pulls the aggregate synthetic-treated / synthetic-control trajectories from raw.globres, the (T0, blank_periods) split from panel.T0 and panel.blank_periods, the inference object from raw.globres.inference, and the covariate matrix from panel.covariates (when present).

mlsynth.utils.post_fit.compute_post_fit_tbrmm(treated_series, control_series, *, n_fit: int, n_post: int, n_treated_units: int, ci: Tuple[float, float] | None = None, inference_method: str | None = None) → SyntheticControlPostFit#

Adapt a measured TBRMM design into a SyntheticControlPostFit.

A sibling of compute_post_fit_marex(), for the same reason: the generic builder takes two trajectories, and TBR has no weight vector to build them from. Its counterfactual comes from the augmented DiD regression instead, so control_series here is the fitted projection \(\hat\delta_1 + \hat\delta_2 \bar{y}_{co,t}\) averaged over the treated geos, not a weighted donor combination.

There is no blank window: the search holds pretest periods out inside its own A/A gate, and those periods are part of the fit here.

Parameters:

  • treated_series, control_series (np.ndarray) – The treated geos’ average observed path and the average of their counterfactuals, over the whole panel.

  • n_fit (int) – Pretest periods the regressions were fitted on.

  • n_post (int) – Realized periods measured.

  • n_treated_units (int) – Treated geos behind the average.

  • ci (tuple of float, optional) – Lower and upper bound for the realized effect, in the units ate is reported in. The design supplies TBR’s posterior interval rescaled to the mean per-period effect, so report.inference and report.effects.att sit on one scale.

  • inference_method (str, optional) – What produced ci, recorded on the result.

Returns:

SyntheticControlPostFit – The standard bundle, so report is built by the same adapter every other design estimator uses.

mlsynth.utils.post_fit.compute_power_analysis(post_fit: SyntheticControlPostFit, *, alpha: float = 0.05, power_target: float = 0.8, post_grid: Sequence[int] | None = None) → PowerAnalysis#

Analytical MDE + power curve for a design’s SyntheticControlPostFit.

Uses the placebo / blank-period gap residuals (or the pre-period gap when no blank window was carved out) to estimate the noise standard deviation sigma_placebo and the AR(1) autocorrelation rho, then computes the minimum detectable effect for each horizon T in post_grid via the Gaussian formula

MDE(T) = (z_{1-alpha/2} + z_{power}) * sigma_placebo * sqrt(VIF(T, rho)),

where VIF(T, rho) = Var(mean of T AR(1) periods) / sigma_placebo^2. The headline MDE uses T = post_fit.n_post (the realised post window).

Parameters:

  • post_fit (SyntheticControlPostFit) – The standardized post-fit from any MAREX-family estimator.

  • alpha (float, default 0.05) – Two-sided significance level.

  • power_target (float, default 0.80) – Target power for the MDE.

  • post_grid (sequence of int, optional) – Post-treatment horizons at which to compute MDE. Defaults to a small geometric grid centered on post_fit.n_post so users see the detectability tradeoff vs. running the experiment longer.

Returns:

PowerAnalysis – Headline MDE + a curve over the requested horizons.

mlsynth.utils.post_fit.compute_smd(cov_matrix: ndarray, treated_weights: ndarray, control_weights: ndarray, *, cov_names: Sequence[str] | None = None, cov_scales: ndarray | None = None) → Dict[str, Any]#

Standardized mean differences between weighted treated and control means.

Parameters:

  • cov_matrix (ndarray, shape (N, M)) – Per-unit covariate values; rows align to treated_weights and control_weights.

  • treated_weights, control_weights (ndarray, shape (N,)) – Non-negative weights with disjoint supports. They are renormalised to sum to 1 internally (so callers may pass raw sums-to-K weights).

  • cov_names (sequence of str, optional) – Names for the M covariates. Defaults to ("cov_0", "cov_1", ...).

  • cov_scales (ndarray, shape (M,), optional) – Pre-computed per-covariate standardization scales (cross-unit std). Defaults to the std of cov_matrix columns. Passing the value already cached by build_covariate_matrix is the right move.

Returns:

  • dict with keys smd (the per-covariate dict), smd_abs_max,

  • and smd_squared_sum. Returns empty / NaN summaries if either weight

  • vector is all-zero.

mlsynth.utils.post_fit.to_effect_result(pf: SyntheticControlPostFit, *, time_periods: ndarray | None = None, intervention_time: Any | None = None, method_name: str | None = None, donor_weights: Dict[str, float] | None = None) → Any#

Convert a SyntheticControlPostFit into a standardized EffectResult.

The single, family-wide adapter from the rich post-fit bundle to the contract’s EffectResult view, so every MAREX-family estimator (LEXSCM, MAREX, SYNDES, PANGEO) gets report for free instead of hand-copying fields. The realized effect’s standard scalars populate the standard sub-models; everything the contract has no slot for (per-period effects, cumulative effect, covariate SMDs, and the full post_fit object itself) is carried in additional_outputs so it remains discoverable.

References#

Abadie, A., & Zhao, J. (2026). “Synthetic Controls for Experimental Design.” See [ABADIE2024].

Abadie, A., Diamond, A., & Hainmueller, J. (2010). “Synthetic Control Methods for Comparative Case Studies.” Journal of the American Statistical Association 105(490):493-505.

Chernozhukov, V., Wuthrich, K., & Zhu, Y. (2021). “An Exact and Robust Conformal Inference Method for Counterfactual and Synthetic Controls.” Journal of the American Statistical Association 116(536):1849-1864.

Vives-i-Bastida, J. (2022). “Synthetic Experimental Design for a UBI Pilot Study.” Report, MIT.

Abadie, A., & L’Hour, J. (2021). “A Penalized Synthetic Control Estimator for Disaggregated Data.” Journal of the American Statistical Association 116(536):1817-1834.

Wolfe, P. (1976). “Finding the Nearest Point in a Polytope.” Mathematical Programming 11(1):128-149.

Frank, M., & Wolfe, P. (1956). “An Algorithm for Quadratic Programming.” Naval Research Logistics Quarterly 3(1-2):95-110.

Neumaier, A., & Shcherbina, O. (2004). “Safe Bounds in Linear and Mixed-Integer Linear Programming.” Mathematical Programming 99:283-296.

Burer, S. (2009). “On the Copositive Representation of Binary and Continuous Nonconvex Quadratic Programs.” Mathematical Programming 120:479-495.

Frangioni, A., & Gentile, C. (2006). “Perspective Cuts for a Class of Convex 0-1 Mixed Integer Programs.” Mathematical Programming 106:225-236.

Günlük, O., & Linderoth, J. (2010). “Perspective Reformulations of Mixed Integer Nonlinear Programs with Indicator Variables.” Mathematical Programming 124:183-205.

Jansson, C., Chaykin, D., & Keil, C. (2007). “Rigorous Error Bounds for the Optimal Value in Semidefinite Programming.” SIAM Journal on Numerical Analysis 46(1):180-200.

Wei, L., Gómez, A., & Küçükyavuz, S. (2022). “Ideal Formulations for Constrained Convex Optimization Problems with Indicator Variables.” Mathematical Programming 192:57-88.

Zheng, X., Sun, X., & Li, D. (2014). “Improving the Performance of MIQP Solvers for Quadratic Programs with Cardinality and Minimum Threshold Constraints: A Semidefinite Program Approach.” INFORMS Journal on Computing 26(4):690-703.

Higham, N. J. (2002). Accuracy and Stability of Numerical Algorithms, 2nd ed. SIAM.

Lei, J., G’Sell, M., Rinaldo, A., Tibshirani, R. J., & Wasserman, L. (2018). “Distribution-Free Predictive Inference for Regression.” Journal of the American Statistical Association 113(523):1094-1111.

Vovk, V., Gammerman, A., & Shafer, G. (2005). Algorithmic Learning in a Random World. Springer.

Phipson, B., & Smyth, G. K. (2010). “Permutation P-values Should Never Be Zero: Calculating Exact P-values When Permutations Are Randomly Drawn.” Statistical Applications in Genetics and Molecular Biology 9(1), Article 39.

Ahn, S. C., & Horenstein, A. R. (2013). “Eigenvalue Ratio Test for the Number of Factors.” Econometrica 81(3):1203-1227.

Contents