Double Machine Learning: Estimating a Causal Parameter Instead of Predicting Returns

Every modelling article on this blog so far poses the same shape of question: given features, predict a number, then validate that the prediction survives out of sample. Spread models predict spread. Fill models predict fill probability. The whole validation apparatus — purged walk-forward, deflated Sharpe, look-ahead taxonomy — exists to check whether a prediction is real.

This article poses a different shape of question, and it is the one piece of machinery the blog has never had: estimate a single scalar parameter that has a causal interpretation, and attach a standard error to it that survives the fact that a flexible ML model was used to get there.

That is not a rhetorical distinction. "Queue position predicts fill probability" is trivially true and operationally useless — of course it does, both are driven by depth and volatility. "Moving one position forward in the queue causes a change of $\theta$ in fill probability, holding market state fixed" is a number you can put into an order placement policy. The first is a regression fit. The second requires an estimator that does not exist in the standard ML toolbox, because regularization and overfitting in a flexible first stage bias exactly the coefficient you care about.

Double Machine Learning (Chernozhukov, Chetverikov, Demirer, Duflo, Hansen, Newey & Robins, 2018) is the estimator that fixes this. Every quantitative trader has heard it: correlation is not causation. DML is the part that comes after that sentence.

Why Naive Regression Fails

Set up the question in the blog's own territory. We want the causal effect of queue position on fill probability:

• $Y_i$: whether the resting limit order at observation $i$ filled within the horizon.
• $D_i$: queue position at placement (normalized by level size).
• $X_i$: the confounders — realized volatility, quoted spread, depth imbalance, level size, time of day, regime label. These are the market-state variables the blog already computes in spread modeling with machine learning.

The causal parameter is $\theta_0$ in

$Y_i = D_i \theta_0 + g_0(X_i) + \epsilon_i$

where $g_0(X_i)$ captures the (potentially complex, nonlinear) relationship between market state and fill outcome.

The treatment is not randomly assigned. You are near the front of the queue because the level was thin, or because you posted during a calm stretch, or because the book was imbalanced in your favour. The same conditions independently drive whether you fill. That is the confounding.

Approach 1: Ignore confounders. Regress $Y$ on $D$ alone. The estimate absorbs the effect of every confounder correlated with both. Textbook omitted variable bias: thin levels give you both a good queue slot and a high fill rate, so you overstate the value of the slot itself.

Approach 2: Linear regression with controls. Regress $Y$ on $D$ and $X$. This works only if $g_0(X)$ is truly linear. Book dynamics are not — the fill/volatility relationship has thresholds, the depth-imbalance effect flips sign by regime. Mis-specifying $g_0$ reintroduces bias.

Approach 3: ML prediction. Fit a gradient boosted model on $(D, X)$. You get good out-of-sample discrimination, and no causal interpretation whatsoever. The model captures every predictive pattern, causal or not; regularization shrinks the treatment's contribution in ways that bias $\theta$; and there is no standard error you can trust.

This is the core tension. ML is good at prediction, but naive application to a causal parameter produces biased, non-normal, unreliable estimates.

The Partially Linear Model

DML operates within a structural framework. The workhorse is the partially linear regression (PLR):

$Y = D\theta_0 + g_0(X) + U, \quad \mathbb{E}[U \mid X, D] = 0$

$D = m_0(X) + V, \quad \mathbb{E}[V \mid X] = 0$

• $\theta_0$ is the causal parameter of interest.
• $g_0(X)$ is a nuisance function — the part of the outcome explained by market state.
• $m_0(X)$ is another nuisance function — the conditional expectation of treatment given market state (the "propensity" in a continuous-treatment setting).
• $U$ and $V$ are structural residuals.

The key insight: $\theta_0$ is low-dimensional, but $g_0$ and $m_0$ can be arbitrarily complex. We want ML to handle the nuisance functions while still delivering valid inference on the scalar.

Why "Double"?

Two ML models, not one:

1. Outcome model: $\hat{\ell}(X) \approx \mathbb{E}[Y \mid X]$ — predict the outcome from market state alone.
2. Treatment model: $\hat{m}(X) \approx \mathbb{E}[D \mid X]$ — predict the treatment from market state alone.

Form residuals

$\tilde{Y}_i = Y_i - \hat{\ell}(X_i), \quad \tilde{D}_i = D_i - \hat{m}(X_i)$

and estimate $\theta_0$ by regressing $\tilde{Y}$ on $\tilde{D}$:

$\hat{\theta}_0 = \frac{\sum_i \tilde{D}_i \tilde{Y}_i}{\sum_i \tilde{D}_i^2}$

This is Frisch-Waugh-Lovell on steroids: partial out the confounders with ML instead of a linear projection, then read the treatment effect off the residual variation.

Neyman Orthogonality: Why It Works

The naive partial-out approach (estimate $g_0$, subtract, regress) fails because ML estimation errors in $\hat{g}_0$ propagate straight into $\hat\theta_0$. The DML score is constructed to be Neyman orthogonal — insensitive to small perturbations in the nuisance functions.

The orthogonal score for PLR:

$\psi(W; \theta, \eta) = \big[Y - \ell(X) - \theta(D - m(X))\big] \cdot \big[D - m(X)\big]$

where $\eta = (\ell, m)$. The orthogonality condition is

$\left. \frac{\partial}{\partial \eta} \mathbb{E}\big[\psi(W; \theta_0, \eta)\big] \right|_{\eta = \eta_0} = 0$

Intuitively, the score uses only the variation in $D$ and $Y$ that is independent of $X$, and errors in one nuisance function are offset by the other. If $\hat{m}$ slightly over-predicts the treatment, $\tilde{D}$ is slightly too small, but the corresponding error in $\tilde{Y}$ from mis-estimating $\ell$ pushes in a compensating direction. The bias becomes second order — the product of two first-stage errors — rather than first order.

Formally, if both nuisance estimators converge at rate $n^{-1/4}$ (mild; most reasonable ML methods clear it), then

$\sqrt{N}(\hat{\theta}_0 - \theta_0) \xrightarrow{d} \mathcal{N}(0, \sigma^2)$

so $\hat\theta_0$ converges at the parametric rate and is asymptotically normal.

Cross-Fitting: Why It Is Mandatory Here

Orthogonality alone is not enough. If the nuisance models are fit on the same rows used to estimate $\theta_0$, first-stage overfitting contaminates the second stage — and the specific harm is worth stating precisely, because it is not the harm you are used to. Elsewhere, overfitting shows up as an inflated validation score: you notice it, you discount it, you move on. Here it shows up as a shifted point estimate for $\theta_0$, with a confidence interval that is still narrow and still centred on the wrong number. There is no score to be suspicious of. The estimator simply lies quietly.

Cross-fitting breaks the dependence: each observation's nuisance predictions come from a model trained without it, and $\theta_0$ is estimated from the pooled held-out residuals. The mechanics are ordinary K-fold machinery, covered in spread modeling with machine learning; what matters below is which folds you hand it.

The DML Algorithm Step by Step

Input: data $\{(Y_i, D_i, X_i)\}_{i=1}^N$, ML methods $\mathcal{M}_\ell$ and $\mathcal{M}_m$, folds $K$.

Step 1 — Partition: split $\{1, \ldots, N\}$ into $K$ disjoint folds.

Step 2 — Cross-fit nuisance models: for $k = 1, \ldots, K$, train $\hat{\ell}^{(-k)}$ and $\hat{m}^{(-k)}$ on the complement of fold $k$, then compute $\tilde{Y}_i = Y_i - \hat{\ell}^{(-k)}(X_i)$ and $\tilde{D}_i = D_i - \hat{m}^{(-k)}(X_i)$ for $i \in I_k$.

Step 3 — Estimate:

$\hat{\theta}_0 = \left(\sum_{i=1}^N \tilde{D}_i^2\right)^{-1} \sum_{i=1}^N \tilde{D}_i \tilde{Y}_i$

Step 4 — Inference:

$\hat{\sigma}^2 = \frac{1}{N} \sum_{i=1}^N \big(\tilde{Y}_i - \hat{\theta}_0 \tilde{D}_i\big)^2 \tilde{D}_i^2 \bigg/ \left(\frac{1}{N} \sum_{i=1}^N \tilde{D}_i^2\right)^2$

with intervals $\hat{\theta}_0 \pm z_{\alpha/2} \cdot \hat{\sigma} / \sqrt{N}$.

The confidence interval is valid for exactly one question

This is the caveat that decides whether a DML result is worth anything, and it is where most applied uses of the method quietly fall apart.

The asymptotic normality above is a statement about one pre-specified treatment, one pre-specified confounder set, one pre-specified score. Fix those in advance, run the estimator once, and the interval means what it says. Try three candidate treatments, or four confounder sets, or swap learners until the p-value looks better, and you are no longer doing inference — you are running a search, and the reported p-value is the p-value of a maximum, not of a draw.

The blog has already measured what that does. In the deflated Sharpe ratio study, searches over pure noise with zero true edge produce a naive false-discovery rate of 1.000 — the unadjusted test fires every single time — while the median naive p-value of the winner sits near 0.0007. Nothing about Neyman orthogonality protects you from this. Orthogonality fixes the bias from nuisance estimation; it says nothing about the bias from specification search. A DML p-value of 1e-05 obtained after trying six specifications deserves exactly the same Bonferroni/Holm/BHY treatment as any other winner pulled out of a grid, with $M$ set to the number of specifications you actually ran.

This has a direct practical consequence for the convenience features of the tooling. DoubleMLData accepts a list in d_cols and will happily return three treatment effects in one summary table:

dml_data_multi = dml.DoubleMLData(
    df, y_col='filled', d_cols=['queue_pos', 'toxicity', 'spread_at_post'],
    x_cols=confounder_cols,
)

Three rows, three p-values, and a multiple-testing problem the summary table does not mention. If you read all three, adjust all three. If only one was the pre-registered question, say so, and treat the other two explicitly as exploratory.

Cross-Fitting on Time Series: Purge, Embargo, Custom Folds

Standard DML assumes i.i.d. observations. Book data is not, and the failure mode is the one this blog has already documented in detail: adjacent rows share overlapping forward windows, so a plain TimeSeriesSplit still leaks the answer across the fold boundary. See spread modeling with machine learning for the purged, embargoed walk-forward implementation and the reason a gap of at least horizon rows is required, and the look-ahead bias taxonomy for the full catalogue of leaks and their measured magnitudes.

The genuinely DML-specific part is how you hand purged folds to the estimator, because set_sample_splitting has a contract that trips people up:

import numpy as np
import doubleml as dml

def purged_folds(n: int, n_splits: int, horizon: int):
    """Expanding-window folds with a purge/embargo gap of `horizon` rows.

    Same construction as the purged walk-forward CV in the spread-modeling
    article: the gap removes the overlap between a training row's forward
    window and a validation row's.
    """
    fold_size = n // (n_splits + 1)
    for k in range(1, n_splits + 1):
        train_end = fold_size * k
        val_start = train_end + horizon
        val_end = val_start + fold_size
        if val_end > n:
            break
        yield np.arange(0, train_end - horizon), np.arange(val_start, val_end)

folds = list(purged_folds(len(df), n_splits=5, horizon=HORIZON))

dml_plr = dml.DoubleMLPLR(dml_data, ml_l=ml_l, ml_m=ml_m)
dml_plr.set_sample_splitting([folds])
dml_plr.fit()

Two things to be aware of, neither of which is obvious from the library docs:

1. Purged walk-forward folds do not cover every row. The purge gaps and the initial training block are never anyone's test fold, so $\theta_0$ is estimated from strictly fewer than $N$ residuals. That is correct behaviour, not a bug, but it means the effective $N$ in the variance formula is the number of pooled test rows — check it rather than assuming.
2. n_rep repetitions are not free here. With random K-fold, repeating cross-fitting and averaging is a cheap variance reduction. With a deterministic time-ordered split there is only one splitting, so n_rep buys nothing and hides nothing; stability has to come from re-running on different data windows instead.

For panel structures (many symbols over the same period), DoubleML supports cluster-robust standard errors — cluster on symbol, not on time, and see multi-symbol validation for the blog's position on when a cross-instrument result is actually established.

The Measured Case: Queue Position and Fill Probability

This is the one causal question in the article where the project already has the data, and it should be run rather than posed. Queue position analysis already covers position estimation, FIFO mechanics, drain rates and time-to-fill on real book data; fill simulation covers fill probability modelling and the calibration loop against live fills. Both produce a predictive model of fills. DML turns the same inputs into a causal estimate.

The specification, pre-registered before looking at the estimate:

• Outcome $Y$: filled within HORIZON snapshots (binary).
• Treatment $D$: normalized queue position at placement.
• Confounders $X$: realized 1s volatility, quoted spread in bps, depth imbalance, level size at post, distance from mid in ticks, time-of-day encoding, regime label.

from sklearn.ensemble import GradientBoostingRegressor
import doubleml as dml

confounder_cols = [
    'rv_1s', 'spread_bps', 'depth_imbalance', 'level_size',
    'dist_from_mid_ticks', 'tod_sin', 'tod_cos', 'regime',
]

dml_data = dml.DoubleMLData(
    df, y_col='filled', d_cols='queue_pos_norm', x_cols=confounder_cols,
)

ml_l = GradientBoostingRegressor(n_estimators=300, max_depth=5)
ml_m = GradientBoostingRegressor(n_estimators=300, max_depth=5)

dml_plr = dml.DoubleMLPLR(dml_data, ml_l=ml_l, ml_m=ml_m, score='partialling out')
dml_plr.set_sample_splitting([folds])
dml_plr.fit()

print(dml_plr.summary)
dml_plr.sensitivity_analysis()
print(dml_plr.sensitivity_summary)

The result to report is a four-column comparison, not a single coefficient: the naive OLS estimate of $Y$ on $D$ alone, the linear-control estimate of $Y$ on $D$ and $X$, the DML estimate with its standard error, and the robustness value from the sensitivity analysis — how strong an unobserved confounder would have to be to nullify the effect. The gap between the naive and DML columns is the quantity of real interest: it is how much of the apparent value of queue position was market state wearing a costume.

A null or negative $\theta$ here is a publishable result and fits this blog better than a clean positive one. If the causal effect of queue position collapses once volatility and level size are partialled out, that is a direct finding about order placement policy: the slot is not what is earning the fill, the conditions under which you got the slot are.

Causal Factor Analysis

The standard approach to factor investing is associational: sort by a characteristic, form long-short portfolios, observe that returns differ. DML enables a different test — estimate the direct causal effect of a characteristic on returns, netting out the confounding effects of the other characteristics. If the effect vanishes under DML, the factor is not independently causal; it is a proxy.

This is a second, independent form of factor skepticism, and it is worth being explicit about how it relates to the one the blog already publishes. The deflated Sharpe ratio attacks the factor zoo from the selection side: with enough trials, a factor can look significant purely because you looked many times. DML attacks it from the confounding side: a factor can look significant on a single honest test and still be a proxy for something else in the conditioning set. A factor has to survive both to be interesting, and the two failure modes are independent — passing one tells you nothing about the other.

What DML Cannot Do

1. It requires the confounders to be observed. If an unobserved variable drives both treatment and outcome, DML is biased, and no amount of ML sophistication fixes an identification problem. Sensitivity analysis bounds the risk; it does not remove it.

2. It estimates an average effect. If the effect of queue position varies sharply across regimes, the point estimate is the average over your sample's regime mix. For heterogeneity use the Interactive Regression Model (DoubleMLIRM) or a causal forest.

3. It assumes a structural model. The partially linear specification requires the treatment to enter the outcome equation in a particular way. If the true process is fundamentally different, DML is confidently wrong.

4. It does not discover causal structure. DML estimates the effect of a pre-specified treatment. It does not tell you which variables are causes.

5. It does not exempt you from multiple testing. Repeating the point above because it is the one most often skipped: orthogonality debiases nuisance estimation, not specification search.

Practical Notes

Sample size. DML needs the nuisance models to converge at $n^{-1/4}$, which in practice means enough rows for the ML models to approximate $g_0$ and $m_0$ at all. Rather than trust round-number thresholds, establish adequacy empirically the way multi-symbol validation does — check whether the estimate holds up across instruments and sub-periods, and treat instability as the signal it is.

Learner choice. The DML-specific fact is narrow but useful: given $n^{-1/4}$ convergence, the learner affects the efficiency of $\hat\theta_0$ (interval width), not its consistency. Which learners are worth reaching for on tabular market data, and why gradient boosting is the default, is already covered in spread modeling with machine learning. If $\hat\theta_0$ moves materially across learners, that is not a menu to choose from — it is evidence the nuisance functions are badly estimated, and per the section above, picking the friendliest one turns the exercise into a search.

Conclusion

DML gives the blog something it did not have: a way to state a market-microstructure claim as a causal parameter with a defensible standard error, instead of as a prediction with a good validation score.

The three load-bearing ideas are:

• Orthogonalize the score so first-stage errors cancel to second order.
• Cross-fit, with purged and embargoed folds on time-series data, so first-stage overfitting cannot shift $\theta_0$.
• Pre-specify, so the interval you report is the interval you actually earned.

The estimator is the easy part. The hard parts are unchanged: deciding which confounders matter, arguing that the identification assumptions hold, and resisting the urge to run the specification one more time.

References

DML's lineage is short and worth one line: it is Robinson's (1988) partially linear model with ML replacing the kernel estimators, achieving the semiparametric efficiency bound, with an orthogonality condition tracing to Neyman's C($\alpha$) test and a close cousin in the targeted-learning (TMLE) literature. Chernozhukov et al.'s contribution was showing this could be operationalized with arbitrary ML learners while retaining $\sqrt{N}$-consistent, asymptotically normal inference.

1. Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W., & Robins, J. (2018). Double/Debiased Machine Learning for Treatment and Structural Parameters. The Econometrics Journal, 21(1), C1-C68.
2. Robinson, P. M. (1988). Root-N-Consistent Semiparametric Regression. Econometrica, 56(4), 931-954.
3. Bach, P., Chernozhukov, V., Kurz, M. S., & Spindler, M. (2022). DoubleML — An Object-Oriented Implementation of Double Machine Learning in Python. Journal of Machine Learning Research, 23(53), 1-6.
4. Facure, M. (2022). Causal Inference for the Brave and True. Chapter 22: Debiased/Orthogonal Machine Learning.
5. Cahan, E., Bai, J., & Ng, S. (2024). Causal Factor Investing. Quantitative Finance.