Dynamic Treatment Effects

Regression blips, parallel-trends SNMMs, and dynamic covariate balance

Dynamic causal estimators must respect the order in which treatment, outcomes, and time-varying confounders are observed. Two methods below parameterize the effect of a treatment blip at one time on later outcomes; the third estimates a mean final outcome under an entire treatment path. All avoid interpreting an ordinary autoregressive distributed-lag regression as causal when it conditions on variables caused by earlier treatment.

This page documents three deliberately distinct estimators:

Two timing conventions

The methods use different observation orderings, so their array contracts are intentionally different.

Contemporaneous regression blips

RegressionBlip treats \(A_t\) as capable of affecting \(Y_t\). Both arrays therefore have shape (n_units, n_periods). Let \(\bar A_t=(A_0,\ldots,A_t)\). For an additive, history-invariant impulse response,

\[ b_t(\bar a_t,j;\gamma) = E\left[ Y_t(\bar a_{t-j},\bar 0_j)-Y_t(\bar a_{t-j-1},\bar 0_{j+1}) \mid \bar A_{t-j}=\bar a_{t-j} \right] =\gamma_j a_{t-j}. \]

After estimating lags \(0,\ldots,j-1\), define the partially blipped-down outcome

\[ \widetilde Y_t^{\,j} =Y_t-\sum_{s=0}^{j-1}\widehat\gamma_s A_{t-s}. \]

At stage \(j\), the implementation regresses \(\widetilde Y_t^{\,j}\) on \(A_{t-j}\) and information available when \(A_{t-j}\) was assigned: \(A_{t-j-1}\), \(Y_{t-j-1}\), optional user-supplied history \(L_{t-j}\), and optional outcome-time indicators. The coefficient on \(A_{t-j}\) is \(\widehat\gamma_j\). The recursion avoids conditioning on \(L_t\) when estimating the effect of \(A_{t-j}\) if \(L_t\) could have been caused by that treatment.

The identifying restriction is sequential ignorability. For every treatment history and assignment time \(s\),

\[ \{Y_r(\bar a_r):r\ge s\}\perp A_s\mid V_s, \]

where \(V_s\) contains the treatment and covariate history sufficient to remove confounding of \(A_s\). The assumption permits \(A_s\) to affect future covariates and future treatment.

Orthogonal estimating equation

Let \(\pi_m(L_m)=E(A_m\mid L_m,\bar A_{m-1})\) and

\[ \nu_{mk}(L_m;\psi)=E\{H^\dagger_{mk}(\psi)\mid L_m,\bar A_{m-1}\}. \]

Shahn et al.’s doubly robust moment is

\[ U_{mk}(\psi) =s_{mk}(L_m,\bar A_{m-1}) \{A_m-\pi_m(L_m)\} \{H^\dagger_{mk}(\psi)-\nu_{mk}(L_m;\psi)\}. \]

Its expectation is zero if either the conditional treatment mean or the conditional blipped-trend mean is correctly specified. For a linear blip, write

\[ H^\dagger_{mk}(\psi)=\Delta Y_k-D_{mk}'\psi. \]

The implementation separately residualizes \(\Delta Y_k\) and \(D_{mk}\) on the history features. This produces a linear orthogonal equation,

\[ \sum_{i,m,k} Z_{imk}\left[ \{\Delta Y_{ik}-\widehat q_{mk}(L_{im})\} -\{D_{imk}-\widehat r_{mk}(L_{im})\}'\psi \right]=0, \]

where \(Z_{imk}\) places \(A_{im}-\widehat\pi_m(L_{im})\) in the coordinate for horizon \(k-m\). Treatment, trend, and blip-design conditional means are fitted with ridge regressions on training folds and evaluated on held-out units. The covariance is an influence-function sandwich over unit-level sums of the orthogonal scores.

Dynamic covariate balance

Viviano and Bradic target the mean final potential outcome under a complete treatment path,

\[ \mu_T(d_{1:T})=E\{Y_T(d_{1:T})\}. \]

Let \(H_{it}\) contain information observed before \(A_{it}\), including baseline covariates, prior treatment, and any relevant lagged outcomes or time-varying covariates. Their potential local-projection model is

\[ E\{Y_{iT}(d_{1:T})\mid H_{it}(d_{1:t-1})\} =H_{it}(d_{1:t-1})'\beta_d^{(t)}. \]

The fitted projections are obtained backwards. At \(t=T\), regress \(Y_T\) on \(H_T\) among units with \(A_{1:T}=d_{1:T}\). At every earlier \(t\), regress the next fitted conditional mean \(H_{t+1}'\widehat\beta_d^{(t+1)}\) on \(H_t\) among units with \(A_{1:t}=d_{1:t}\). The implementation uses ridge regularization with an unpenalized intercept; the paper’s high-dimensional theory instead permits appropriately convergent lasso projections.

Starting from \(\widehat\gamma_{i0}=1/n\), period-\(t\) weights solve the exact-balance special case

\[ \begin{aligned} \widehat\gamma_t =\arg\min_{\gamma\in\mathbb R^n}\;&\sum_{i=1}^n\gamma_i^2\\ \text{subject to }\;& \sum_i\gamma_iH_{it} =\sum_i\widehat\gamma_{i,t-1}H_{it},\\ &\sum_i\gamma_i=1,\qquad 0\leq\gamma_i\leq w_{\max},\\ &\gamma_i=0\quad\text{if }A_{i,1:t}\neq d_{1:t}. \end{aligned} \]

Thus the first-period path group is balanced to the whole sample. At later periods, the longer path group is balanced to the preceding path group after applying the previous period’s weights. DynamicCovariateBalance calls the same autoscaling, bounded quadratic dual map, solver sequence, normalization, and diagnostics used by BalancingWeights(objective="quadratic"); this is one shared calibration engine rather than a second weighting implementation.

Given the recursive predictions \(\widehat q_{it}=H_{it}'\widehat\beta_d^{(t)}\), the augmented path-mean estimator is

\[ \widehat\mu_T(d) =\sum_i\widehat\gamma_{iT}Y_{iT} -\sum_{t=2}^T\sum_i (\widehat\gamma_{it}-\widehat\gamma_{i,t-1})\widehat q_{it} -\sum_i(\widehat\gamma_{i1}-1/n)\widehat q_{i1}. \]

The paper permits coordinatewise \(\ell_\infty\) balance slack and explicit high-dimensional tuning. The first package implementation uses exact balance, ridge projections, and the existing scalar weight cap. It reports calibration diagnostics but deliberately does not report an analytic standard error yet; Monte Carlo dispersion or a unit bootstrap is required for uncertainty until the paper’s longitudinal influence-function variance is implemented.

API

import crabbymetrics as cm

# Blackwell--Glynn recursive regression blips.
reg = cm.RegressionBlip(max_lag=2, time_effects=True)
reg.fit(y_same_period, treatment_same_period, history)
reg.coef
reg.stage_standard_errors
reg.blip_down(y_same_period, treatment_same_period)
reg.summary()

# Shahn et al. parallel-trends SNMM.
pt = cm.ParallelTrendsSNMM(
    max_horizon=3,
    treatment_mode="initiation",
    n_folds=3,
    nuisance_penalty=1e-6,
    propensity_clip=0.01,
    seed=42,
)
pt.fit(y_forward, absorbing_treatment_status, history)
pt.coef
pt.standard_errors
pt.summary()

# Viviano--Bradic dynamic covariate balance for one treatment path.
dcb = cm.DynamicCovariateBalance(
    nuisance_penalty=1e-6,
    autoscale=True,
    max_weight=1.0,
)
dcb.fit(final_outcome, treatment, history, target_path=[0, 1, 1, 1])
dcb.potential_outcome
dcb.get_weights()  # n_units by n_periods
dcb.summary()

treatment_mode="blip" uses the supplied period-specific treatment directly. treatment_mode="initiation" requires a binary absorbing treatment-status panel, constructs a first-treatment pulse internally, and restricts each assignment-time equation to units not treated before that time. In initiation mode, \(\psi_h\) is the effect at horizon \(h\) of initiating at \(m\) rather than remaining untreated from \(m\) onward, among the relevant initiation history.

DynamicCovariateBalance.fit(...) takes only the final outcome because its estimand is a path-specific final mean. The treatment panel and target path have the same number of periods; history[:, t, :] must contain only pre-treatment information for period \(t\). A path contrast requires two fitted objects, one for each path. summary() reports prefix support, effective sample size, maximum original-scale imbalance, per-period solver status, and the fitted path mean.

Numerical experiment 1: post-treatment covariate feedback

The first design makes \(Z_t\) a mediator of \(A_{t-1}\) and a confounder for \(A_t\):

\[ Z_t=0.8A_{t-1}+U_t, \qquad \Pr(A_t=1\mid Z_t,A_{t-1})=\operatorname{logit}^{-1}(-0.2+0.8Z_t+0.2A_{t-1}), \]

\[ Y_t=2A_t+0.5A_{t-1}+1.5Z_t+\varepsilon_t. \]

The contemporaneous blip is \(\gamma_0=2\). The first-lag blip includes both the direct path and the mediated path through \(Z_t\):

\[ \gamma_1=0.5+1.5(0.8)=1.7. \]

Code 
import numpy as np
import pandas as pd
import crabbymetrics as cm

rng = np.random.default_rng(1234)
n, periods = 2_000, 7
a = np.zeros((n, periods))
y = np.zeros((n, periods))
z = np.zeros((n, periods, 1))

for t in range(periods):
    a_lag = a[:, t - 1] if t else 0.0
    z[:, t, 0] = 0.8 * a_lag + rng.normal(size=n)
    p = 1 / (1 + np.exp(-(-0.2 + 0.8 * z[:, t, 0] + 0.2 * a_lag)))
    a[:, t] = rng.binomial(1, p)
    y[:, t] = (
        2.0 * a[:, t]
        + 0.5 * a_lag
        + 1.5 * z[:, t, 0]
        + rng.normal(scale=0.7, size=n)
    )

reg = cm.RegressionBlip(max_lag=1)
reg.fit(y, a, z)
reg_result = reg.summary()
pd.DataFrame({
    "lag": reg_result["lags"].astype(int),
    "truth": [2.0, 1.7],
    "estimate": reg_result["coef"],
    "conditional_stage_se": reg_result["stage_se"],
})
lag truth estimate conditional_stage_se
0 0 2.0 1.981614 0.013754
1 1 1.7 1.637764 0.035863

An ADL regression that includes \(Z_t\) while interpreting the coefficient on \(A_{t-1}\) as its total lagged effect blocks the mediated path and targets \(0.5\), not \(1.7\). Omitting \(Z_t\) instead confounds the contemporaneous effect. The recursive blip regressions change the adjustment time with the causal lag.

Numerical experiment 2: staggered initiation with level confounding

The second design contains an unobserved \(U_i\) that affects both treatment initiation and the untreated outcome level. It does not affect untreated increments. The observed history \(L_{im}\) predicts both initiation and untreated increments:

\[ Y_{i0}(0)=2U_i+\varepsilon_{i0}, \qquad Y_{i,m+1}(0)-Y_{im}(0)=0.35L_{im}+\varepsilon_{i,m+1}, \]

\[ \Pr(T_i=m\mid T_i\ge m,L_{im},U_i) =\operatorname{logit}^{-1}(-1.2+0.55L_{im}+0.65U_i). \]

The initiation response is \((1.0,1.6,2.1)\) at horizons one through three. This violates level ignorability because \(U_i\) confounds initiation and \(Y_{i0}(0)\), but it satisfies the simulated conditional parallel-trends restriction.

Code 
rng = np.random.default_rng(91)
n, periods = 5_000, 6
truth = np.array([1.0, 1.6, 2.1])
history = rng.normal(size=(n, periods, 1))
u = rng.normal(size=n)
status = np.zeros((n, periods))
first = np.full(n, periods + 1, dtype=int)
untreated = np.zeros((n, periods + 1))
untreated[:, 0] = 2.0 * u + rng.normal(scale=0.4, size=n)

for m in range(periods):
    untreated[:, m + 1] = (
        untreated[:, m]
        + 0.35 * history[:, m, 0]
        + rng.normal(scale=0.5, size=n)
    )
    risk = first > m
    p = 1 / (1 + np.exp(-(-1.2 + 0.55 * history[:, m, 0] + 0.65 * u)))
    start = risk & (rng.uniform(size=n) < p)
    first[start] = m
    status[first <= m, m] = 1.0

y_forward = untreated.copy()
for i in range(n):
    if first[i] < periods:
        for h, effect in enumerate(truth, start=1):
            if first[i] + h <= periods:
                y_forward[i, first[i] + h] += effect

pt = cm.ParallelTrendsSNMM(
    max_horizon=3,
    treatment_mode="initiation",
    n_folds=3,
    seed=7,
)
pt.fit(y_forward, status, history)
pt_result = pt.summary()
pd.DataFrame({
    "horizon": pt_result["horizons"].astype(int),
    "truth": truth,
    "estimate": pt_result["coef"],
    "se": pt_result["se"],
})
horizon truth estimate se
0 1 1.0 0.983159 0.009227
1 2 1.6 1.589484 0.016232
2 3 2.1 2.103710 0.022914
diagnostic
moment rows 4.526500e+04
minimum fitted treatment mean 1.000000e-02
maximum fitted treatment mean 6.400915e-01
maximum absolute sample moment 4.106337e-15

Comparative simulation study

The three estimators target different primitive objects. A common scalar comparison is possible in a linear distributed-lag design. There are four treatment decisions, \(t=0,1,2,3\), and the target contrasts are

\[ d^0=(0,0,0,0),\qquad d^1=(0,1,1,1). \]

Because the final outcome is \(Y_4\), switching from \(d^0\) to \(d^1\) changes the treatments entering at horizons one, two, and three. The common target is therefore

\[ \tau=\mu_4(d^1)-\mu_4(d^0) =1+0.6+0.4=2. \]

For each unit, draw mutually independent standard-normal \(X_i,U_i,L_{i,-1},\eta_{it},\varepsilon_{it}\). Let \(\Lambda(v)=1/(1+e^{-v})\), set \(A_{i,-1}=A_{i,-2}=0\), and generate

\[ L_{it}=0.5X_i+\eta_{it}, \]

\[ A_{it}\mid X_i,L_{it},A_{i,t-1},U_i \sim\operatorname{Bernoulli}\left[ \Lambda\{-0.1+0.35X_i+0.35L_{it}+0.2A_{i,t-1}+\rho U_i\} \right], \]

\[ Y_{i0}=U_i+0.5X_i+0.4L_{i,-1}+\varepsilon_{i0}, \]

\[ Y_{i,t+1} =U_i+0.5X_i+0.15(t+1)+0.4L_{it} +A_{it}+0.6A_{i,t-1}+0.4A_{i,t-2}+\varepsilon_{i,t+1}. \]

The observed-selection design sets \(\rho=0\). The history \(H_{it}=(X_i,L_{it},L_{i,t-1},A_{i,t-1})\) then contains the treatment predictors. The fixed \(U_i\) shifts outcome levels but is independent of treatment.

Code 
flowchart LR
  X["baseline X"] --> L["state L_t"]
  X --> A["treatment A_t"]
  L --> A
  AP["past treatment A_(t-1)"] --> A
  X --> Y["outcome Y_(t+1)"]
  L --> Y
  AP --> Y
  A --> Y

flowchart LR
  X["baseline X"] --> L["state L_t"]
  X --> A["treatment A_t"]
  L --> A
  AP["past treatment A_(t-1)"] --> A
  X --> Y["outcome Y_(t+1)"]
  L --> Y
  AP --> Y
  A --> Y

The hidden-selection design sets \(\rho=0.6\). The unobserved \(U_i\) now causes both treatment and outcome levels. It does not enter the untreated trend because it cancels from \(Y_{i,t+1}(0)-Y_{it}(0)\). Sequential ignorability fails for the regression-blip and dynamic-balancing estimators, while the time-varying parallel-trends restriction remains plausible.

Code 
flowchart LR
  U["unobserved level U"] --> A["treatment A_t"]
  U --> Y["outcome Y_(t+1)"]
  X["baseline X"] --> L["state L_t"]
  X --> A
  L --> A
  AP["past treatment A_(t-1)"] --> A
  X --> Y
  L --> Y
  AP --> Y
  A --> Y

flowchart LR
  U["unobserved level U"] --> A["treatment A_t"]
  U --> Y["outcome Y_(t+1)"]
  X["baseline X"] --> L["state L_t"]
  X --> A
  L --> A
  AP["past treatment A_(t-1)"] --> A
  X --> Y
  L --> Y
  AP --> Y
  A --> Y

For each draw, RegressionBlip estimates \(\tau\) as \(\widehat\gamma_0+\widehat\gamma_1+\widehat\gamma_2\). ParallelTrendsSNMM estimates it as \(\widehat\psi_1+\widehat\psi_2+\widehat\psi_3\). DynamicCovariateBalance fits \(\widehat\mu_4(d^1)\) and \(\widehat\mu_4(d^0)\) separately and differences them.

Code 
import matplotlib.pyplot as plt

def simulate_comparison(seed, hidden_strength, n=1_200, periods=4):
    rng = np.random.default_rng(seed)
    x = rng.normal(size=n)
    u = rng.normal(size=n)
    treatment = np.zeros((n, periods))
    outcome = np.zeros((n, periods + 1))
    state = np.zeros((n, periods))
    lagged_state = np.zeros((n, periods))

    state_minus_one = 0.5 * x + rng.normal(size=n)
    outcome[:, 0] = (
        u + 0.5 * x + 0.4 * state_minus_one + rng.normal(scale=0.8, size=n)
    )
    previous_treatment = np.zeros(n)
    twice_lagged_treatment = np.zeros(n)
    previous_state = state_minus_one

    for time in range(periods):
        lagged_state[:, time] = previous_state
        state[:, time] = 0.5 * x + rng.normal(size=n)
        linear_probability = (
            -0.1
            + 0.35 * x
            + 0.35 * state[:, time]
            + 0.2 * previous_treatment
            + hidden_strength * u
        )
        propensity = 1.0 / (1.0 + np.exp(-linear_probability))
        treatment[:, time] = rng.binomial(1, propensity)
        outcome[:, time + 1] = (
            u
            + 0.5 * x
            + 0.15 * (time + 1)
            + 0.4 * state[:, time]
            + treatment[:, time]
            + 0.6 * previous_treatment
            + 0.4 * twice_lagged_treatment
            + rng.normal(scale=0.8, size=n)
        )
        twice_lagged_treatment = previous_treatment.copy()
        previous_treatment = treatment[:, time].copy()
        previous_state = state[:, time].copy()

    lagged_treatment = np.column_stack(
        [np.zeros(n), treatment[:, : periods - 1]]
    )
    pt_history = np.stack(
        [np.repeat(x[:, None], periods, axis=1), state, lagged_state, lagged_treatment],
        axis=2,
    )
    regression_treatment = np.zeros((n, periods + 1))
    regression_treatment[:, 1:] = treatment
    regression_history = np.zeros((n, periods + 1, 3))
    regression_history[:, 1:, 0] = x[:, None]
    regression_history[:, 1:, 1] = state
    regression_history[:, 1:, 2] = lagged_state
    dcb_history = np.stack(
        [
            np.repeat(x[:, None], periods, axis=1),
            state,
            lagged_state,
            lagged_treatment,
            outcome[:, :periods],
        ],
        axis=2,
    )

    regression = cm.RegressionBlip(max_lag=2)
    regression.fit(outcome, regression_treatment, regression_history)
    parallel_trends = cm.ParallelTrendsSNMM(
        max_horizon=3,
        treatment_mode="blip",
        n_folds=3,
        nuisance_penalty=1e-4,
        seed=seed,
    )
    parallel_trends.fit(outcome, treatment, pt_history)

    path_means = []
    dcb_success = []
    dcb_max_balance = []
    for path in ([0.0, 0.0, 0.0, 0.0], [0.0, 1.0, 1.0, 1.0]):
        dcb = cm.DynamicCovariateBalance(nuisance_penalty=1e-4)
        dcb.fit(outcome[:, -1], treatment, dcb_history, path)
        diagnostics = dcb.summary()
        path_means.append(dcb.potential_outcome)
        dcb_success.append(diagnostics["success"])
        dcb_max_balance.append(np.max(diagnostics["max_abs_balance"]))

    return {
        "Regression blip": regression.coef.sum(),
        "Parallel-trends SNMM": parallel_trends.coef.sum(),
        "Dynamic balance": path_means[1] - path_means[0],
        "dcb_success": all(dcb_success),
        "dcb_max_balance": max(dcb_max_balance),
    }


records = []
for design, hidden_strength in [
    ("Selection on observables", 0.0),
    ("Selection on unobservables", 0.6),
]:
    for replication in range(100):
        result = simulate_comparison(1_000 + replication, hidden_strength)
        for estimator in [
            "Regression blip",
            "Parallel-trends SNMM",
            "Dynamic balance",
        ]:
            records.append(
                {
                    "design": design,
                    "replication": replication,
                    "estimator": estimator,
                    "estimate": result[estimator],
                    "dcb_success": result["dcb_success"],
                    "dcb_max_balance": result["dcb_max_balance"],
                }
            )

comparison = pd.DataFrame(records)
truth = 2.0
monte_carlo = (
    comparison.groupby(["design", "estimator"])["estimate"]
    .agg(mean="mean", mc_sd="std")
    .reset_index()
)
monte_carlo["bias"] = monte_carlo["mean"] - truth
monte_carlo["rmse"] = (
    comparison.assign(squared_error=lambda frame: (frame["estimate"] - truth) ** 2)
    .groupby(["design", "estimator"])["squared_error"]
    .mean()
    .pow(0.5)
    .to_numpy()
)
monte_carlo.round(3)
design estimator mean mc_sd bias rmse
0 Selection on observables Dynamic balance 2.000 0.197 0.000 0.196
1 Selection on observables Parallel-trends SNMM 1.997 0.095 -0.003 0.095
2 Selection on observables Regression blip 1.996 0.076 -0.004 0.075
3 Selection on unobservables Dynamic balance 2.918 0.201 0.918 0.939
4 Selection on unobservables Parallel-trends SNMM 1.957 0.086 -0.043 0.095
5 Selection on unobservables Regression blip 2.636 0.075 0.636 0.641

Sampling distributions for the common final-period path contrast.

Monte Carlo bias and standard deviation separate robustness from efficiency.

Under selection on observables, all three sampling distributions are centered near two. The recursive regression is most efficient in this correctly specified low-dimensional outcome model; the parallel-trends estimator pays for orthogonalization and cross-fitting, while the path-specific dynamic-balancing contrast pays for estimating two increasingly selective path means.

When \(U_i\) enters treatment assignment, the regression-blip and dynamic-balancing estimates inherit the positive level confounding. The parallel-trends estimate remains close to the target because differencing removes \(U_i\). This is not generic robustness to hidden confounding: it relies specifically on \(U_i\) being time invariant and absent from untreated outcome trends.

strict_solver_success_rate maximum_retained_imbalance
design
Selection on observables 0.81 0.000130
Selection on unobservables 0.76 0.000053

BalancingWeights deliberately retains finite weights when its strict scaled-residual convergence gate is not met. DynamicCovariateBalance propagates that per-period status instead of relabeling the fit as converged. The simulation therefore reports both the strict solver-success rate and the maximum original-scale imbalance among retained fits; a production analysis should inspect both and tighten or relax the calibration design deliberately.

Interpretation and current boundaries

These classes do not estimate the same object under interchangeable assumptions.

Feature RegressionBlip ParallelTrendsSNMM DynamicCovariateBalance
Identifying restriction Sequential ignorability Time-varying conditional parallel trends Sequential ignorability plus potential projections
Primitive target Lag-specific additive blips Horizon-specific additive blips Mean final outcome under one path
Outcome timing \(A_t\) may affect \(Y_t\) \(A_m\) first affects \(Y_{m+1}\) History precedes \(A_t\); final outcome supplied separately
Nuisance strategy Recursive outcome regressions Cross-fitted treatment and trend regressions Recursive potential projections plus sequential calibration
Robustness Outcome-regression specification Doubly robust orthogonal moment Product of projection error and imbalance
Unobserved level confounding Generally not allowed Allowed when it does not confound untreated trends Generally not allowed
Inference Stagewise unit-clustered SE Unit-level influence-function sandwich Not yet implemented; use a unit bootstrap

The initial implementation is intentionally narrow:

  1. The two blip estimators use additive linear effects, with one coefficient per lag or horizon. Dynamic balance instead estimates one path mean at a time.
  2. ParallelTrendsSNMM uses ridge nuisance regressions and a conditional-mean treatment model. Flexible Python callbacks are not yet part of the Rust boundary.
  3. RegressionBlip reports unit-clustered uncertainty conditional on earlier recursive blip estimates. A unit block bootstrap is required for joint recursive inference.
  4. In ParallelTrendsSNMM, cross-fitting is over units, never unit-period rows.
  5. DynamicCovariateBalance implements the paper’s full-interaction recursion with ridge, exact balance, and the package’s quadratic calibration geometry. It does not yet implement the paper’s lasso tuning, coordinatewise approximate-balance program, or analytic longitudinal variance.
  6. None of the estimators makes parallel trends or sequential ignorability testable. The analyst still owns the history set, timing, treatment coding, positivity argument, path support, and effect specification.

References

  • Blackwell, Matthew, and Adam N. Glynn. 2018. “How to Make Causal Inferences with Time-Series Cross-Sectional Data under Selection on Observables.” American Political Science Review 112(4): 1067–1082. Paper.
  • Shahn, Zach, Oliver Dukes, Meghana Shamsunder, David Richardson, Eric Tchetgen Tchetgen, and James Robins. 2022. “Structural Nested Mean Models Under Parallel Trends Assumptions.” arXiv:2204.10291.
  • Viviano, Davide, and Jelena Bradic. 2026. “Dynamic Covariate Balancing: Estimating Treatment Effects over Time with Potential Local Projections.” First circulated 2021. arXiv:2103.01280.