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.