Autobounds in practice

Setup, solver mechanics, worked identification examples, and a small scaling probe

1 Overview

duarteguilherme/autobounds is a tool for asking a precise partial-identification question in discrete causal models: among all latent response-function distributions that agree with the observed data and the stated assumptions, how small or large can the target causal estimand be?

This note works through the package in that spirit. It starts with a reproducible local setup, then explains the response-function and optimization machinery, then solves a set of demo-style causal designs to show which estimands are point identified and which remain partially identified. The final section runs a small scaling probe to see what happens as discrete covariate support grows.

The core reference is Duarte, Finkelstein, Knox, Mummolo, and Shpitser, “An Automated Approach to Causal Inference in Discrete Settings”. The practical version is:

causal graph + observed distribution + assumptions + estimand
  -> polynomial program
  -> lower and upper causal bounds

The examples below are not a Monte Carlo study. They are a checked, executable tour of the workflow and the solver behavior.

2 Reproducible setup

The package can be cloned directly with the GitHub CLI:

gh repo clone duarteguilherme/autobounds autobounds

There is one macOS-specific wrinkle. The repository currently contains both:

autobounds/Bounder.py
autobounds/bounder.py

On a case-insensitive filesystem these paths collide. The lowercase file is a compatibility wrapper and the uppercase file contains the implementation, but the cleanest local setup is still to preserve both. The run below uses a small case-sensitive APFS disk image and clones into that mount:

base=/Users/alal/tmp/autobounds_lalten_20260617
img="$base/autobounds_case.sparseimage"
mount="$base/case_mount"

hdiutil create -size 4g -fs 'Case-sensitive APFS' \
  -volname autobounds_case_20260617 "$img"
hdiutil attach "$img" -mountpoint "$mount" -nobrowse

gh repo clone duarteguilherme/autobounds "$mount/autobounds"

The package declares requires-python >=3.10, while /usr/bin/python3 on this machine is 3.9. The run uses uv to create a clean Python 3.11 environment:

cd /Users/alal/tmp/autobounds_lalten_20260617/case_mount/autobounds
uv venv --python 3.11 .venv
uv pip install --python .venv/bin/python -e .

For this Quarto note, the minimal Jupyter stack is installed into the same environment:

uv pip install --python .venv/bin/python ipykernel nbformat nbclient jupyter

The resulting environment was:

autobounds_version python platform maxtime_per_solve_sec script
0 0.0.2 3.11.11 macOS-26.4.1-arm64-arm-64bit 3.0000 /Users/alal/tmp/autobounds_lalten_20260617/cas...

The first sanity check is the smallest point-identified example: a binary treatment and binary outcome with no confounding in the graph.

from autobounds import DAG, causalProblem
import pandas as pd

d = DAG("D -> Y")
problem = causalProblem(d)
df = pd.DataFrame({
    "D": [0,0,1,0,1,0,1,1,1,0,1,0,0,1],
    "Y": [0,1,1,0,0,0,1,0,1,1,0,0,1,1],
})
problem.load_data(raw=df)
problem.set_ate("D", "Y")
problem.solve(verbose_result=False, verbose_optimizer=False, maxtime=10)

That returns a point-identified ATE of about 0.142857 for this toy dataset. The later examples deliberately move away from this clean case into settings where the graph and observed distribution do not pin down a unique causal effect.

3 Autobounds methodology

Autobounds is easiest to understand as an automated principal-strata or response-function engine for discrete causal models.

Given a DAG, it builds a canonical nonparametric structural model. Each observed variable is represented as a deterministic response function of its parents and latent background variables. In a binary IV model, for example, treatment response types include never-takers, always-takers, compliers, and defiers. In larger graphs the same idea scales to many more response-function types.

The unknowns are probabilities over latent response-function types. The observed data impose constraints on those probabilities. Causal assumptions add more constraints: graphical restrictions, monotonicity, exclusion restrictions, measurement assumptions, missingness assumptions, or user-specified equalities and inequalities. The causal estimand is also expressible as a polynomial in the same unknown probabilities.

So the inferential problem becomes two optimization problems:

lower bound = minimize estimand over all feasible latent response-type distributions
upper bound = maximize estimand over all feasible latent response-type distributions

If the estimand is point identified, these two optima coincide. If it is partially identified, the two optima define the identified set. Internally, the package writes polynomial programs and solves them with SCIP/PySCIPOpt, using primal and dual bounds from spatial branch-and-bound.

3.1 Optimizer backend

The default optimizer path in this checkout is open-source SCIP through PySCIPOpt. No commercial solver such as Gurobi or Mosek was used in this run.

The relevant implementation is Program.run_scip() in autobounds/Program.py. In outline, it:

  1. imports Model from pyscipopt;
  2. builds two separate SCIP models, one for the lower bound and one for the upper bound;
  3. creates one variable per latent response-type probability, bounded between 0 and 1;
  4. creates an objvar variable for the causal estimand;
  5. adds the polynomial equality/inequality constraints generated by Autobounds;
  6. sets the objectives as minimize objvar and maximize objvar;
  7. launches the lower and upper solves in separate worker processes;
  8. returns SCIP’s primal and dual bounds.

In simplified form:

from pyscipopt import Model

M_lower = Model()
M_upper = Model()

for p in parameters:
    if p != "objvar":
        var_lower[p] = M_lower.addVar(p, lb=0.0, ub=1.0)
        var_upper[p] = M_upper.addVar(p, lb=0.0, ub=1.0)
    else:
        var_lower[p] = M_lower.addVar(p, lb=-1, ub=1)
        var_upper[p] = M_upper.addVar(p, lb=-1, ub=1)

M_lower.setObjective(var_lower["objvar"], sense="minimize")
M_upper.setObjective(var_upper["objvar"], sense="maximize")

In this local run the installed stack was:

Autobounds: 0.0.2
PySCIPOpt: 6.2.1
SCIP: 10.0.x

There are older or alternate methods in the repository, including run_couenne() through Pyomo/Couenne and a generic run_pyomo(solver_name="ipopt"). Those are not the solve path used here. causalProblem.solve() routes to the SCIP/PySCIPOpt path.

3.2 What branch-and-bound means here

It is tempting to describe this as “integer programming” because SCIP is famous as a mixed-integer optimizer. The Autobounds models in this note are not ordinary MILPs, though. The main unknowns are continuous probabilities over latent response-function types. The nonconvexity comes from polynomial products of these probabilities, especially when constraints or estimands multiply response-type terms.

So the relevant optimization problem is better described as a nonconvex polynomial program over a probability simplex-like feasible set. SCIP attacks this using spatial branch-and-bound / branch-and-cut:

  • primal bound: the best feasible response-type distribution found so far, hence an actual feasible data-generating process and an attainable estimand value;
  • dual bound: a certified relaxation bound showing how low or high the estimand could still be over an unexplored or relaxed feasible region;
  • branching: split the variable domain into smaller boxes;
  • bounding: solve relaxations on each box to rule out regions that cannot improve the current bound;
  • convergence: continue until the primal-dual gap is small, the requested threshold is reached, or the time limit stops the run.

This is why the result dictionary distinguishes point lb dual, point ub dual, point lb primal, and point ub primal. For a lower-bound problem, the dual bound is the conservative certified lower endpoint, while the primal bound is the best feasible lower-side witness found so far. For an upper-bound problem, the dual bound is the conservative certified upper endpoint, while the primal bound is the best feasible upper-side witness.

When SCIP proves optimality, primal and dual agree. When it times out, the dual bounds are the guaranteed outer interval and the primal bounds are feasible inner witnesses. The tables below report the dual endpoints because those are the conservative Autobounds-style partial-identification bounds.

For covariates, the current read_data(..., covariates=...) path is important. With discrete covariates it solves the same causal-bounds problem within each observed covariate stratum, then aggregates the stratum-specific bounds by empirical stratum weights. That means the runtime here scales mostly with the number of observed covariate strata, not with sample size directly.

4 A running synthetic design

The scaling probe uses a small IV-like design:

Code 
flowchart LR
  Z[Instrument Z] --> D[Treatment D]
  D --> Y[Outcome Y]
  U((Unobserved U)) -.-> D
  U -.-> Y
  classDef latent fill:#f7f7f7,stroke:#666,stroke-dasharray: 4 3;
  class U latent;

flowchart LR
  Z[Instrument Z] --> D[Treatment D]
  D --> Y[Outcome Y]
  U((Unobserved U)) -.-> D
  U -.-> Y
  classDef latent fill:#f7f7f7,stroke:#666,stroke-dasharray: 4 3;
  class U latent;

U is unobserved, so the ATE of D on Y is not point identified from the observed joint distribution of (Z,D,Y). Optional discrete covariates X0, X1, ... affect treatment and outcome probabilities in the simulated data. Autobounds then treats those covariates as observed discrete strata.

The data-generating process used in the run:

u ~ Bernoulli(0.5)
z ~ Bernoulli(0.5)

p_d = clip(0.18 + 0.42*z + 0.22*u + 0.12*x_score, 0.01, 0.99)
d ~ Bernoulli(p_d)

p_y0 = clip(0.08 + 0.25*u + 0.10*x_score, 0.01, 0.90)
p_y1 = clip(p_y0 + 0.22, 0.01, 0.99)
y ~ Bernoulli(p_y1 if d == 1 else p_y0)

With no covariates, the solved interval is:

label rows observed_covariate_strata naive_difference lb_dual ub_dual width elapsed_sec status
0 simple_no_covariates 2000 1 0.2703 -0.0947 0.4684 0.5631 0.5336 ok

The naive difference in means sits inside the partial-identification interval, but the interval is much wider because the graph permits unobserved treatment-outcome confounding.

5 Worked designs from the demo notebooks

The repository’s autobounds_demo/ folder is useful because it repeats the same core workflow across very different causal designs:

DAG -> causalProblem -> load_data/read_data -> set_ate or set_estimand -> write_program/solve

The main variation is not the surface API. It is the graph, the target estimand, and the extra constraints imposed on the latent response types. Each numbered subsection below is self-contained: DAG, data construction or demo-data load, Autobounds setup, solve call, and solved interval.

For this note, a design is classified as numerically point identified when the SCIP lower and upper dual endpoints agree to tolerance:

point identified if abs(upper_dual - lower_dual) <= 1e-4
partially identified otherwise

This is a numerical criterion for these runs, not a theorem about every possible dataset under each named graph. It is still the right operational diagnostic: if lower and upper collapse, the assumptions and observed distribution pin down the estimand; if they do not, Autobounds is exposing the remaining compatible causal worlds.

label classification interval width elapsed_sec status
0 Unconfounded ATE point [0.279, 0.279] 0.0000 0.5300 ok
1 Confounded observational ATE partial [-0.364, 0.636] 1.0000 0.5400 ok
2 Standard IV ATE partial [0.015, 0.508] 0.4930 0.5400 ok
3 IV LATE, no defiers point [0.334, 0.334] 0.0000 0.5400 ok
4 Cautious IV ATE partial [-0.322, 0.678] 1.0000 0.5500 ok
5 Ternary-instrument IV ATE partial [-0.500, 0.500] 1.0000 0.5400 ok
6 Outcome selection ATE partial [-0.500, 0.640] 1.1410 0.5300 ok
7 Measurement-error ATE partial [-0.620, 1.000] 1.6200 0.5300 ok
8 Mediation CDE at Z=0 partial [-0.628, 0.372] 1.0000 0.5400 ok
9 KLM conditional ATE partial [-0.997, 0.382] 1.3790 10.5700 ok

Two rows in the summary table are calibration cases rather than demo-gallery subsections: the unconfounded ATE, which should collapse, and the confounded observational ATE, which should not.

classification count
0 partial 8
1 point 2

The unconfounded ATE collapses because the DAG says treatment is as-good-as-random in the toy setup. The confounded observational ATE does not collapse because hidden U can explain both treatment assignment and outcomes. The IV examples then show the more interesting split: LATE is point identified after adding no-defiers and targeting compliers, while the full-population ATE remains partially identified.

5.1 Standard IV

From autobounds_demo/iv_demo.ipynb, autobounds_demo/apsa2023/iv_demo.ipynb, and autobounds_demo/iv_more_categories_demo.ipynb.

Code 
flowchart LR
  Uz((Uz)) -.-> Z[Instrument Z]
  Z --> X[Treatment X]
  X --> Y[Outcome Y]
  Uxy((Uxy)) -.-> X
  Uxy -.-> Y
  classDef latent fill:#f7f7f7,stroke:#666,stroke-dasharray: 4 3;
  class Uz,Uxy latent;

flowchart LR
  Uz((Uz)) -.-> Z[Instrument Z]
  Z --> X[Treatment X]
  X --> Y[Outcome Y]
  Uxy((Uxy)) -.-> X
  Uxy -.-> Y
  classDef latent fill:#f7f7f7,stroke:#666,stroke-dasharray: 4 3;
  class Uz,Uxy latent;

Data and solution code:

from autobounds import DAG, causalProblem
import numpy as np
import pandas as pd

rng = np.random.default_rng(101)
n = 3_000
z = rng.binomial(1, 0.5, n)

# Response type: 0 never-taker, 1 complier, 2 always-taker.
typ = rng.choice([0, 1, 2], size=n, p=[0.25, 0.50, 0.25])
x0 = np.where(typ == 2, 1, 0)
x1 = np.where(typ == 0, 0, 1)
x = np.where(z == 1, x1, x0)

p_y0 = np.clip(0.08 + 0.08 * typ, 0.01, 0.90)
p_y1 = np.clip(p_y0 + 0.28, 0.01, 0.99)
y = rng.binomial(1, np.where(x == 1, p_y1, p_y0))

iv_data = pd.DataFrame({"Z": z, "X": x, "Y": y})
data_summary = iv_data.value_counts(sort=False).reset_index(name="n")
data_summary["prob"] = data_summary["n"] / data_summary["n"].sum()
data_summary = data_summary.drop(columns="n")

dag = DAG()
dag.from_structure(
    edges="Uz -> Z, Z -> X, X -> Y, Uxy -> X, Uxy -> Y",
    unob="Uz, Uxy",
)

ate = causalProblem(dag)
ate.load_data(data_summary)
ate.set_ate(ind="X", dep="Y")
ate_result = ate.solve(
    ci=False,
    maxtime=8,
    theta=0.005,
    limits=[-1, 1],
    verbose_optimizer=False,
    verbose_result=False,
)

late = causalProblem(dag)
late.load_data(data_summary)
late.add_constraint(late.p("X(Z=0)=1&X(Z=1)=0"))  # no defiers
late.set_ate(
    ind="X",
    dep="Y",
    cond="X(Z=0)=0&X(Z=1)=1",
)
late_result = late.solve(
    ci=False,
    maxtime=8,
    theta=0.005,
    limits=[-1, 1],
    verbose_optimizer=False,
    verbose_result=False,
)

Solved intervals:

label classification interval width elapsed_sec status
2 Standard IV ATE partial [0.015, 0.508] 0.4930 0.5400 ok
3 IV LATE, no defiers point [0.334, 0.334] 0.0000 0.5400 ok

The full-population ATE remains partially identified. The complier LATE is point identified after monotonicity/no-defiers because the estimand matches what the IV assumptions pin down.

5.2 Cautious IV with a direct instrument-outcome path

The IV demo also includes a more cautious design that relaxes the exclusion restriction by allowing Z -> Y.

Code 
flowchart LR
  Uz((Uz)) -.-> Z[Instrument Z]
  Z --> X[Treatment X]
  Z --> Y[Outcome Y]
  X --> Y
  Uxy((Uxy)) -.-> X
  Uxy -.-> Y
  classDef latent fill:#f7f7f7,stroke:#666,stroke-dasharray: 4 3;
  class Uz,Uxy latent;

flowchart LR
  Uz((Uz)) -.-> Z[Instrument Z]
  Z --> X[Treatment X]
  Z --> Y[Outcome Y]
  X --> Y
  Uxy((Uxy)) -.-> X
  Uxy -.-> Y
  classDef latent fill:#f7f7f7,stroke:#666,stroke-dasharray: 4 3;
  class Uz,Uxy latent;

Data and solution code:

from autobounds import DAG, causalProblem
import numpy as np
import pandas as pd

rng = np.random.default_rng(101)
n = 3_000
z = rng.binomial(1, 0.5, n)
typ = rng.choice([0, 1, 2], size=n, p=[0.25, 0.50, 0.25])
x0 = np.where(typ == 2, 1, 0)
x1 = np.where(typ == 0, 0, 1)
x = np.where(z == 1, x1, x0)
p_y0 = np.clip(0.08 + 0.08 * typ, 0.01, 0.90)
p_y1 = np.clip(p_y0 + 0.28, 0.01, 0.99)
y = rng.binomial(1, np.where(x == 1, p_y1, p_y0))

iv_data = pd.DataFrame({"Z": z, "X": x, "Y": y})
data_summary = iv_data.value_counts(sort=False).reset_index(name="n")
data_summary["prob"] = data_summary["n"] / data_summary["n"].sum()
data_summary = data_summary.drop(columns="n")

dag = DAG()
dag.from_structure(
    edges="Uz -> Z, Z -> X, Z -> Y, X -> Y, Uxy -> X, Uxy -> Y",
    unob="Uz, Uxy",
)

problem = causalProblem(dag)
problem.load_data(data_summary)
problem.set_ate(ind="X", dep="Y")
result = problem.solve(
    ci=False,
    maxtime=8,
    theta=0.005,
    limits=[-1, 1],
    verbose_optimizer=False,
    verbose_result=False,
)

Solved interval:

label classification interval width elapsed_sec status
4 Cautious IV ATE partial [-0.322, 0.678] 1.0000 0.5500 ok

The important point is methodological: changing the graph changes the feasible response-function distributions. Allowing a direct instrument-outcome effect destroys the standard exclusion restriction, so the full ATE interval widens relative to the cleaner IV graph.

5.3 Multi-category IV

iv_more_categories_demo.ipynb keeps the same IV graph but changes the support size of Z, X, and/or Y.

Code 
flowchart LR
  Z[Z: 2 or 3 levels] --> X[X: 2 or 3 levels]
  X --> Y[Y: 2 or 3 levels]
  Uxy((Uxy)) -.-> X
  Uxy -.-> Y
  classDef latent fill:#f7f7f7,stroke:#666,stroke-dasharray: 4 3;
  class Uxy latent;

flowchart LR
  Z[Z: 2 or 3 levels] --> X[X: 2 or 3 levels]
  X --> Y[Y: 2 or 3 levels]
  Uxy((Uxy)) -.-> X
  Uxy -.-> Y
  classDef latent fill:#f7f7f7,stroke:#666,stroke-dasharray: 4 3;
  class Uxy latent;

Data and solution code:

from autobounds import DAG, causalProblem
import pandas as pd

rows = []
for z in [0, 1, 2]:
    for x in [0, 1]:
        for y in [0, 1]:
            rows.append({"Z": z, "X": x, "Y": y, "prob": 1 / 12})
data_summary = pd.DataFrame(rows)

dag = DAG()
dag.from_structure(
    edges="Z -> X, X -> Y, Uxy -> X, Uxy -> Y",
    unob="Uxy",
)

problem = causalProblem(dag, {"Z": 3})
problem.load_data(data_summary)
problem.set_ate(ind="X", dep="Y")
result = problem.solve(
    ci=False,
    maxtime=8,
    theta=0.005,
    limits=[-1, 1],
    verbose_optimizer=False,
    verbose_result=False,
)

Solved interval:

label classification interval width elapsed_sec status
5 Ternary-instrument IV ATE partial [-0.500, 0.500] 1.0000 0.5400 ok

This is the design analogue of the scaling probe above. Cardinality increases the number of response-function types inside the polynomial program, while discrete covariate support increases the number of stratum-level problems.

5.4 Outcome selection

From outcome_selection_demo.ipynb and apsa2023/outcome_selection_demo.ipynb. The simple demo uses Y -> S, so selection depends on the outcome. The APSA demo adds a separate selection disturbance.

Code 
flowchart LR
  X[Treatment X] --> Y[Outcome Y]
  Y --> S[Selection / observed S]
  Us((Us)) -.-> S
  Ux((Ux)) -.-> X
  Ux -.-> Y
  classDef latent fill:#f7f7f7,stroke:#666,stroke-dasharray: 4 3;
  class Us,Ux latent;

flowchart LR
  X[Treatment X] --> Y[Outcome Y]
  Y --> S[Selection / observed S]
  Us((Us)) -.-> S
  Ux((Ux)) -.-> X
  Ux -.-> Y
  classDef latent fill:#f7f7f7,stroke:#666,stroke-dasharray: 4 3;
  class Us,Ux latent;

Data and solution code:

from autobounds import DAG, causalProblem
import pandas as pd

raw = pd.read_csv("replication_files/data/selection_obsqty.csv")
if "prob" in raw.columns:
    data_summary = raw
else:
    data_summary = raw.value_counts(sort=False).reset_index(name="n")
    data_summary["prob"] = data_summary["n"] / data_summary["n"].sum()
    data_summary = data_summary.drop(columns="n")

dag = DAG()
dag.from_structure(
    edges="Us -> S, Y -> S, X -> Y, Ux -> X, Ux -> Y",
    unob="Us, Ux",
)
problem = causalProblem(dag)
problem.load_data(data_summary)
problem.set_ate(ind="X", dep="Y")
result = problem.solve(
    ci=False,
    maxtime=8,
    theta=0.005,
    limits=[-1, 1],
    verbose_optimizer=False,
    verbose_result=False,
)

Solved interval:

label classification interval width elapsed_sec status
6 Outcome selection ATE partial [-0.500, 0.640] 1.1410 0.5300 ok

This design asks for the treatment effect on the latent/target outcome distribution when the observed data are filtered by a selection process.

5.5 Measurement error / outcome proxy

From measurement_error_demo.ipynb. Here S is a measured proxy or report that depends on true Y, and the measurement process shares unobserved variation with Y.

Code 
flowchart LR
  X[Treatment X] --> Y[True outcome Y]
  Y --> S[Observed proxy S]
  Uys((Uys)) -.-> Y
  Uys -.-> S
  classDef latent fill:#f7f7f7,stroke:#666,stroke-dasharray: 4 3;
  class Uys latent;

flowchart LR
  X[Treatment X] --> Y[True outcome Y]
  Y --> S[Observed proxy S]
  Uys((Uys)) -.-> Y
  Uys -.-> S
  classDef latent fill:#f7f7f7,stroke:#666,stroke-dasharray: 4 3;
  class Uys latent;

Data and solution code:

from autobounds import DAG, causalProblem
import pandas as pd

raw = pd.read_csv("replication_files/data/measurement_error.csv")
if "prob" in raw.columns:
    data_summary = raw
else:
    data_summary = raw.value_counts(sort=False).reset_index(name="n")
    data_summary["prob"] = data_summary["n"] / data_summary["n"].sum()
    data_summary = data_summary.drop(columns="n")

dag = DAG()
dag.from_structure(
    "X -> Y, Y -> S, Uys -> S, Uys -> Y",
    unob="Uys",
)
problem = causalProblem(dag)
problem.load_data(data_summary)
problem.add_constraint(problem.p("S(Y=0)=1&S(Y=1)=0"))
problem.set_ate(ind="X", dep="Y")
result = problem.solve(
    ci=False,
    maxtime=8,
    theta=0.005,
    limits=[-1, 1],
    verbose_optimizer=False,
    verbose_result=False,
)

Solved interval:

label classification interval width elapsed_sec status
7 Measurement-error ATE partial [-0.620, 1.000] 1.6200 0.5300 ok

The estimand is about the true outcome Y, not merely the observed proxy S.

5.6 Mediation and decomposed effects

From mediation_demo.ipynb. The notebook names the upstream treatment Z, the mediator X, and the outcome Y.

Code 
flowchart LR
  Z[Treatment Z] --> X[Mediator X]
  Z --> Y[Outcome Y]
  X --> Y
  Uxy((Uxy)) -.-> X
  Uxy -.-> Y
  classDef latent fill:#f7f7f7,stroke:#666,stroke-dasharray: 4 3;
  class Uxy latent;

flowchart LR
  Z[Treatment Z] --> X[Mediator X]
  Z --> Y[Outcome Y]
  X --> Y
  Uxy((Uxy)) -.-> X
  Uxy -.-> Y
  classDef latent fill:#f7f7f7,stroke:#666,stroke-dasharray: 4 3;
  class Uxy latent;

Data and solution code for the controlled direct effect at Z=0:

from autobounds import DAG, causalProblem
import pandas as pd

raw = pd.read_csv("replication_files/data/iv.csv")
if "prob" in raw.columns:
    data_summary = raw
else:
    data_summary = raw.value_counts(sort=False).reset_index(name="n")
    data_summary["prob"] = data_summary["n"] / data_summary["n"].sum()
    data_summary = data_summary.drop(columns="n")

dag = DAG()
dag.from_structure(
    edges="Z -> X, X -> Y, Z -> Y, Uxy -> X, Uxy -> Y",
    unob="Uxy",
)
problem = causalProblem(dag)
problem.load_data(data_summary)
problem.set_estimand(
    problem.p("Y(X=1,Z=0)=1") -
    problem.p("Y(X=0,Z=0)=1")
)
result = problem.solve(
    ci=False,
    maxtime=8,
    theta=0.005,
    limits=[-1, 1],
    verbose_optimizer=False,
    verbose_result=False,
)

Solved interval:

label classification interval width elapsed_sec status
8 Mediation CDE at Z=0 partial [-0.628, 0.372] 1.0000 0.5400 ok

This is a good example of Autobounds as a symbolic query system: once potential-outcome events are parseable, the estimand can be any polynomial expression over those events.

5.7 KLM stop/use-of-force design

From apsa2023/klm_demo.ipynb and klm_demo_relaxassumptions.ipynb. The variables are race/treatment D, stop/search mediator M, and force outcome Y.

Code 
flowchart LR
  Ud((Ud)) -.-> D[Race / treatment D]
  D --> M[Stop M]
  D --> Y[Force Y]
  M --> Y
  Umy((Umy)) -.-> M
  Umy -.-> Y
  classDef latent fill:#f7f7f7,stroke:#666,stroke-dasharray: 4 3;
  class Ud,Umy latent;

flowchart LR
  Ud((Ud)) -.-> D[Race / treatment D]
  D --> M[Stop M]
  D --> Y[Force Y]
  M --> Y
  Umy((Umy)) -.-> M
  Umy -.-> Y
  classDef latent fill:#f7f7f7,stroke:#666,stroke-dasharray: 4 3;
  class Ud,Umy latent;

Data and solution code:

from autobounds import DAG, Q, causalProblem
import pandas as pd

raw = pd.read_csv("autobounds_demo/apsa2023/data/klm_demo_bw_anyforce.csv")
if "prob" in raw.columns:
    data_summary = raw
else:
    data_summary = raw.value_counts(sort=False).reset_index(name="n")
    data_summary["prob"] = data_summary["n"] / data_summary["n"].sum()
    data_summary = data_summary.drop(columns="n")

dag = DAG()
dag.from_structure(
    edges="Ud -> D, D -> M, D -> Y, M -> Y, Umy -> M, Umy -> Y",
    unob="Ud, Umy",
)
problem = causalProblem(dag)
problem.load_data(data_summary, cond=["M"])
problem.add_constraint(problem.p("M=1") - Q(0.01), ">=")
problem.set_ate(ind="D", dep="Y", cond="M=1")
result = problem.solve(
    ci=False,
    maxtime=8,
    theta=0.005,
    limits=[-1, 1],
    verbose_optimizer=False,
    verbose_result=False,
)

Solved interval:

label classification interval width elapsed_sec status
9 KLM conditional ATE partial [-0.997, 0.382] 1.3790 10.5700 ok

The demo then layers substantive restrictions as polynomial constraints, for example no force without a stop, mediator monotonicity variants, and side information about stop disparities. This is the pattern to remember: graphs set the response-function universe; assumptions and outside information carve away parts of that universe.

6 Scaling experiment

The worked designs above show how identification changes with the graph and estimand. A separate practical question is how the computation scales when the graph is fixed but the observed covariate support grows.

Every scaling case uses the same synthetic graph from the running example and a 3-second SCIP time cap per solve. The number of raw rows is not the main driver here; rows are increased for some high-support cases only to avoid extremely sparse strata. The key variable is observed_covariate_strata, because Autobounds solves the same causal-bounds problem inside each observed discrete stratum and then aggregates by empirical stratum weights.

6.1 Increasing cardinality with fixed covariates

First, fix the number of covariates and increase the number of levels. This grows the support through cardinality.

label rows num_covariates levels_per_covariate observed_covariate_strata lb_dual ub_dual width elapsed_sec
1 cardinality_p1_L2 2000 1 2 2 -0.1144 0.4591 0.5735 1.0958
2 cardinality_p1_L4 2000 1 4 4 -0.1333 0.4934 0.6267 2.1787
3 cardinality_p1_L8 2000 1 8 8 -0.0897 0.4810 0.5707 4.3770
4 cardinality_p1_L16 2000 1 16 16 -0.1206 0.4477 0.5683 8.7505
5 cardinality_p1_L32 3840 1 32 32 -0.1086 0.4710 0.5796 17.4858
6 cardinality_p2_L2 2000 2 2 4 -0.0768 0.4769 0.5538 2.1795
7 cardinality_p2_L3 2000 2 3 9 -0.1117 0.4489 0.5606 4.9205
8 cardinality_p2_L4 2000 2 4 16 -0.0988 0.5027 0.6015 8.6895
9 cardinality_p2_L5 3000 2 5 25 -0.1101 0.4548 0.5649 13.3331

6.2 Increasing the number of covariates

Second, keep each covariate binary and increase the number of covariates. This makes support grow as \(2^p\).

label rows num_covariates levels_per_covariate observed_covariate_strata lb_dual ub_dual width elapsed_sec
10 covariate_count_p1_L2 2000 1 2 2 -0.0692 0.4862 0.5554 1.0905
11 covariate_count_p2_L2 2000 2 2 4 -0.0584 0.4906 0.5490 2.1784
12 covariate_count_p3_L2 2000 3 2 8 -0.0737 0.5037 0.5774 4.2753
13 covariate_count_p4_L2 2000 4 2 16 -0.1287 0.4569 0.5856 8.7780
14 covariate_count_p5_L2 3840 5 2 32 -0.0993 0.4790 0.5783 17.4466

6.3 Timing plot

Elapsed wall time against observed covariate support. Points from both scaling sweeps mostly fall on the same line because this experiment uses Autobounds’ discrete-covariate stratification path.

For this particular stratified path, the timing is close to linear in observed support:

sec_per_stratum
mean 0.5435
min 0.5333
max 0.5486

7 Takeaways

The package installed and solved the examples successfully once the case-sensitive checkout issue was avoided. The dependency most likely to fail, pyscipopt, installed cleanly as a wheel in the Python 3.11 venv.

Methodologically, Autobounds is doing exactly what one wants from a discrete partial-identification engine: it turns a DAG, observed distribution, assumptions, and estimand into lower and upper optimization problems over latent response-type probabilities.

Substantively, the solved design gallery is the useful identification sanity check. Clean unconfounded ATE and monotone-IV LATE are point identified in the constructed examples. Full-population ATE under IV, confounded observational ATE, cautious IV, selection, measurement-error/proxy, mediation, and KLM conditional-effect examples remain partially identified.

Computationally, the small scaling probe suggests a simple first-order rule for categorical covariates: support size matters. Increasing one covariate from 2 to 32 levels took the run from about 1.1s to 17.5s. Increasing binary covariates from p=1 to p=5 produced the same support sequence, 2 to 32 strata, and nearly the same runtime sequence.

That does not mean all Autobounds problems are easy. The harder cases are those where the causal graph and variable cardinalities create a large response-function space inside each stratum. But in this discrete-covariate path, many binary covariates and a few high-cardinality covariates are computationally similar if they induce the same observed support.

8 Reproducibility files

The note is backed by two small scripts and their outputs:

analysis/autobounds_scaling_experiment.py
analysis_outputs/autobounds_scaling_results.csv

analysis/autobounds_solved_designs.py
analysis_outputs/autobounds_solved_designs.csv

The cloned repo commit was:

git rev-parse HEAD
# 4e6d75886a732d46e46a0e549b987670beac2e4b