Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
59,585 characters · 23 sections · 57 citation commands
\spacingset{1}
{\it Keywords:} Doubly robust estimators, augmented balancing weights, time-varying treatment effects, surrogacy, mediation, recursive functionals
\spacingset{1.5}
Many of the counterfactual means estimated in empirical research are identified as recursions of regressions. For example, the mean outcome under a two-period treatment path is identified by two regressions iterated back to back: first, the outcome is regressed on the full history; then the predicted outcome, read at the counterfactual path, is regressed on the first period robins1986. The mean outcome using surrogate variables shares this structure: the surrogate index is the regression of the long-run outcome on the surrogates in one sample, which is then regressed on the treatment and covariates in another sample prentice1989surrogate. The natural indirect effect in mediation analysis is yet another example: first, the outcome is regressed on treatment, mediator and covariates; then the predicted outcome, evaluated at the opposite treatment, is regressed on the treatment and covariates robins1992identifiability, pearl2001direct. Across all of these examples, an initial regression generates a prediction, which is used as a generated outcome in a later regression. The classic identification results are stated as recursive regressions, i.e. with generated outcomes.
While the recursive regression formulation is in terms of the outcome mechanism, another formulation is in terms of the treatment mechanism. For example, the mean outcome under a two-period treatment path can also be expressed as the mean outcome after recursive reweighting by inverse propensity scores. These recursive weights may be viewed as recursively balancing the covariates of the treated and untreated subpopulations, using baseline covariates in the initial balancing and generated covariates in the subsequent balancing. Yet another formulation is doubly robust by augmenting the treatment mechanism formulation (in terms of generated regressors) with the outcome mechanism formulation (in terms of generated outcomes) bang2005doublyrobust. Phrased another way, generated regressors can be used to debias generated outcomes. In summary, empirical researchers studying time-varying treatment effects, surrogate analysis, and mediation analysis face a plethora of estimators.
In cross sectional causal inference, it is well known that certain regression, balancing weight, and doubly robust estimators are numerically equivalent robins_performance. Recent work has shown that, in cross sectional causal inference, debiasing a regularized regression estimator with a balancing weight effectively shrinks the regression towards least squares bruns2025augmented. In this paper, we ask: does the equivalence documented for cross sectional causal inference extend to longitudinal causal inference, where there is the additional richness of generated outcomes and generated regressors? Moreover, does the interpretation of debiasing as shrinkage towards least squares generalize?
Our primary contribution is a positive answer to the first question: estimation by recursively regressing generated outcomes is numerically equivalent to estimation by recursively balancing generated regressors. This equivalence is exact when, in each period, the recursive regressions and recursive balancing weights are linear in the same dictionary of basis functions and unregularized. The equivalence is approximate when, in each period, the recursive regressions and recursive balancing weights are linear in the same dictionary of basis functions but regularized. Several nonparametric estimators are linear in their basis functions, such as series, lasso, and kernel methods; see Appendix (ref).
Our secondary contribution is a partially negative answer to the second question: the shrinkage interpretation of debiasing vanishes as the number of time periods grows. While the doubly robust estimator with ridge-regularized balancing weights is still pulled towards OLS, the strength of that pull decays geometrically in the number of time periods. We conclude that the number of time periods is centrally important in understanding regularization. Only in the special case of cross sectional causal inference (with one time period) is it the case that debiasing and undersmoothing are two sides of the same coin.
Most directly, we build on previous work that characterizes numerical equivalences in cross-sectional causal inference rubin1980randomization,robins_performance, kline2011oaxaca, chattopadhyay2023implied. See bruns2025augmented for a comprehensive review and a general statement for the class of bounded linear functionals considered by chernozhukov2022automatic. Recently, rotnitzky2025onestep generalize this equivalence to the class of mixed bias linear functionals defined by rotnitzky2021mixedbias. We extend the equivalence in a different direction, to the class of recursive functionals of chernozhukov2022nested. Unlike earlier work, the class we study includes the canonical models of time-varying treatment effects, surrogate analysis, and mediation analysis, which are widely used in empirical research.
Whereas previous works on recursive balancing weights propose estimators and analyse their (probabilistic) convergence properties, in this work we analyse their (deterministic) algebraic equivalences. Recursive balancing criteria have been proposed for parametric bang2005doublyrobust, kernel kallus2021optimalbalancing, high dimensional linear viviano2021dynamic, and general machine learning chernozhukov2022nested function spaces. Our contribution is not to propose a new balancing criterion. Instead, we show when the balancing recursion with generated regressors is numerically identical to, or algebraically decomposable against, the regression recursion with generated outcomes.
More broadly, we contribute to a vast literature on doubly robust estimation, by shedding new light on how various estimators are interconnected. The literature establishes orthogonality, consistency, and asymptotic for recursive causal parameters molina2017multiple,luedtke2017sequential,rotnitzky2017multiply,chernozhukov2022nested. We study a complementary question: how do the regression, balancing weight, and doubly robust estimators relate algebraically? We find that the main equivalence from cross sectional causal inference generalizes to longitudinal causal inference, however the intuition that debiasing amounts to under-smoothing does not.
Section (ref) defines the class of recursive functionals, and details three leading examples: time-varying treatment, surrogate, and mediation analysis. Sections (ref), (ref), and (ref) study the cases with no regularization, ridge regularization, and general convex regularization, respectively. Section (ref) concludes.
We observe $n$ i.i.d.\ copies of a random vector $W$, the concatenation of all variables observed across time periods. The parameters we study are built from $T \geq 1$ stages of recursive conditional expectations, and we call $T$ the number of time periods. Following chernozhukov2022nested, each stage $t$ contains three time-varying objects: (i) the conditioning variables $Z_t$; (ii) an expectation $\mathbb{E}_t(\cdot)$, taken with respect to a (possibly) time-varying measure; and (iii) a linear formula $m_t$, defined below. A base expectation $\mathbb{E}_0(\cdot)$ averages over the population that defines the target. In the time-varying treatment example below $Z_t$ records the full observed history and expectation measures do not change. However, in the surrogacy example the treatment leaves the conditioning set between stages and the expectation measure changes from the observational group to the experimental group. Our setting nests both. When every stage shares the population measure, we drop the subscripts and write $\mathbb{E}(\cdot)$. The last stage, $t = T$, is the one with the outcome $Y$ as its response; the first stage, $t = 1$, defines the target.
That target, $\theta_0$, is an average of a counterfactual outcome: the mean outcome had treatment followed a counterfactual path. Units do not necessarily follow the full path, so $\theta_0$ is a functional of a distribution we never fully observe. Time-varying treatment effects, surrogate-identified treatment effects, and mediation effects all take this form, and they are our running examples. Setting $T = 1$ collapses the framework to the single-period problem of chernozhukov2022automatic,bruns2025augmented, which is our reference point throughout.
The recursive structure identifies $\theta_0$ by a backward recursion of conditional expectations. It is the g-computation formula of robins1986 extended to additional counterfactuals.
The recursion starts at the last stage, $t = T$, which regresses the outcome on the variables observed there:
where $Z_T$ collects those final-stage variables. Each outer stage regresses the previous stage's output on its own conditioning set. For $t = T-1, \ldots, 1$,
where $m_{t+1}(W; g)$ is a mean-square-continuous linear functional of $g$. We call $m_{t+1}(W; f_{t+1})$ a generated outcome: it is built from nuisance $f_{t+1}$ and evaluated at an argument that may be counterfactual. The generated outcome then serves as the response one stage further out, so identification iterates backward through the stages.
Taking the first-stage expectation identifies the target by the mean of the first-stage generated outcome,
Any parameter of the form (ref)--(ref) is a recursive functional in the sense of chernozhukov2022nested: a model in this class is specified by the number of time periods $T$ and the time-varying objects $\{Z_t, \mathbb{E}_t(\cdot), m_t\}$. The right-hand side of (ref) is a functional of the observed-data distributions, and it is the object every estimator in this paper targets. What gives the recursion a causal meaning is a separate matter, and the identifying assumptions differ across applications: time-varying treatment, surrogacy, and mediation each rest on their own conditions, stated with each running example in Section (ref). The algebra of this paper needs only the recursive form (ref)--(ref), whatever assumptions give it causal content.
Each outcome regression has a forward-pass Riesz representer $\alpha_t$, defined stage by stage through the Riesz representation theorem chernozhukov2022nested. The Riesz representer is a function, and its evaluations give a vector called the balancing weights. We take $\alpha_t$ to be the representer in the space of functions of $Z_t$ that are square-integrable under the stage-$t$ measure, so it is $Z_t$-measurable by construction, and we write $\alpha_t(Z_t)$ throughout.\footnote{If a representer on the full vector $W$ exists, $\alpha_t(Z_t)$ is its conditional expectation given $Z_t$ under the stage-$t$ measure, the minimal such representer; this convention is what makes the linear form $\alpha_t(Z_t) = \phi_t'\eta_t$ of Section (ref) well-posed, since any representer on $W$ is pinned down only up to additions with zero conditional mean given $Z_t$.} Set $\alpha_0 := 1$. The first representer $\alpha_1$ turns the functional $m_1$ into an inner product:
and for $t = 2, \ldots, T$, $\alpha_t$ represents the functional $L_t(g) := \mathbb{E}_{t-1}[\alpha_{t-1}(Z_{t-1})\, m_t(W; g)]$:
Existence of these representers is at the population level: each represented functional must be mean-square continuous, which in the running examples follows from sequential positivity and overlap between adjacent stage measures. The finite-sample identities of Sections (ref)--(ref) do not need it: they are statements about the coefficient vectors the stage programmes return, and they hold whether or not a population $\alpha_t$ exists. Once $\alpha_{t-1}$ is replaced by an estimate $\hat\alpha_{t-1}$, that estimate enters the stage-$t$ Riesz regression as a generated regressor, the forward-pass counterpart of the generated outcome.
The representers identify $\theta_0$ by a forward recursion. Start from $\theta_0 = \mathbb{E}_0[m_1(W; f_1)] = \mathbb{E}_1[\alpha_1(Z_1) f_1(Z_1)]$ and substitute $f_1(Z_1) = \mathbb{E}_1[m_2(W; f_2) \mid Z_1]$: \[ \theta_0 = \mathbb{E}_1[\alpha_1(Z_1)\, m_2(W; f_2)] = \mathbb{E}_2[\alpha_2(Z_2)\, f_2(Z_2)], \] where the second equality uses (ref). Iterating through all $T$ stages gives the inverse-probability-weighted (IPW) representation
a weighted average of the outcome under the last-stage measure. In this final expression, $\alpha_T(Z_T)$ may be viewed as the balancing weights for generated regressors. Specifically, the regressors being balanced at time $T$ are recursively generated by $\alpha_{t}(Z_t)$ for $t<T$.
Estimating $f_t$ can introduce regularisation bias. In the plug-in estimator $\hat{\mathbb{E}}_0[m_1(W; \hat{f}_1)]$, that bias accumulates from the last stage backward. First-order bias can be eliminated by a Neyman-orthogonal moment, built by applying the Riesz representation theorem recursively.
Write the stage-$t$ residual, i.e. the prediction error of $f_t$ for its generated outcome, as
We write $\hat\varepsilon_t$ for the same residual evaluated at estimated nuisances, with $\hat{f}_{t+1}$ and $\hat{f}_t$ in place of $f_{t+1}$ and $f_t$. Collecting the corrections gives the orthogonal moment chernozhukov2022nested:
Each correction $\mathbb{E}_t[\alpha_t(Z_t)\, \varepsilon_t(W)]$ is taken under its own stage measure and is zero at the true nuisances: $f_t$ is by definition the stage-$t$ conditional mean of its response, so $\mathbb{E}_t[\varepsilon_t(W) \mid Z_t] = 0$, and the weight $\alpha_t(Z_t)$ is a function of the same conditioning set. Orthogonality follows from the bias expansion below, in which every term is a product of two estimation errors. The estimator built from the orthogonal moment attains the semiparametric efficiency bound in many cases.
The recursive moment has a recursive mixed bias property chernozhukov2022nested, the recursive counterpart of the bias structure that rotnitzky2021mixedbias characterise for a single pair of nuisances. Take any alternative nuisances $(f_t^*, \alpha_t^*)$ with errors $\tilde{f}_t = f_t^* - f_t$ and $\tilde\alpha_t = \alpha_t^* - \alpha_t$. Evaluating the moment at the alternatives, the bias is
with the convention $m_{T+1}(W; \tilde{f}_{T+1}) := 0$. The inner error enters as the generated outcome $m_{t+1}(W; \tilde{f}_{t+1})$. Each summand is a product of a representer error and an outcome error, so it vanishes when $\tilde\alpha_t = 0$ (the stage-$t$ representer is correct) or when $\tilde{f}_{t+1} = \tilde{f}_t = 0$ (both adjacent outcome regressions are correct). This stage-wise robustness is the recursive analogue of double robustness: at $T = 1$, equation (ref) is the single term $-\mathbb{E}_1[\tilde\alpha_1 \tilde{f}_1]$, which vanishes if either $\alpha_1$ or $f_1$ is correct. rotnitzky2017multiply establish multiple robustness of this stage-wise kind for estimators of the time-varying treatment estimand of Example 2.
Two special cases of the moment (ref) are given below. At $T = 1$, with a single population measure, it is the standard doubly robust (AIPW) moment
our non-recursive baseline. With distinct base and stage measures it is instead the covariate-shift moment of chernozhukov2023covshifts. At $T = 2$ it expands to
the leading recursive case for the running examples.
We provide four examples that are nested by the recursive framework. The first ($T = 1$) is the non-recursive baseline. The other three ($T = 2$) show how recursion arises in time-varying treatment, surrogacy, and mediation.
\paragraph{Example 1: Single period ($T = 1$).} We observe $n$ i.i.d.\ copies of $W = (Z, Y)$, where $Z$ collects covariates and possibly a treatment indicator $D$. The single nuisance function is $f_1(z) = \mathbb{E}[Y \mid Z = z]$, the single Riesz representer $\alpha_1$ satisfies (ref), and $\theta_0 = \mathbb{E}[m_1(W; f_1)]$. The doubly robust moment reduces to the standard AIPW form (ref). This setting covers the average treatment effect robins1994estimating and other bounded linear functionals chernozhukov2022locallyrobust. bruns2025augmented show that when both the outcome model and the Riesz representer are estimated by ridge regression, the augmented estimator is algebraically equivalent to a single undersmoothed outcome regression; this equivalence is the direct precursor to our recursive results.
\paragraph{Example 2: $T = 2$ dynamic treatment effects.} We observe $n$ i.i.d.\ trajectories $W = (X_1, D_1, X_2, D_2, Y)$, with time-varying covariates $X_t$, binary treatments $D_t$, and final outcome $Y$. We set $Z_1 = (X_1, D_1)$ and $Z_2 = (X_1, D_1, X_2, D_2)$. The parameter of interest is the counterfactual mean $\theta_0(d_1, d_2) := \mathbb{E}[Y(d_1, d_2)]$. Three conditions at each stage $t$ identify $\theta_0$ via the g-computation formula robins1986. (i) If $(D_1,D_2)=(d_1,d_2)$ then $Y=Y(d_1,d_2)$. (ii) Potential outcomes are conditionally independent of the stage-$t$ treatment. (iii) Treatment propensities are bounded away from zero. Under these conditions the nuisance functions are:
where $m_2(W; g) = g(X_1, d_1, X_2, d_2)$ evaluates $g$ at the counterfactual treatment sequence $(d_1, d_2)$, and $\theta_0 = \mathbb{E}[f_1(X_1, d_1)]$, so $m_1(W; g) = g(X_1, d_1)$.
\paragraph{Example 3: Surrogacy.} We observe an experimental sample $(X, D, S)$ and an observational sample $(X, S, Y)$, with binary treatment $D$, surrogate $S$, and long-run outcome $Y$. The estimand is $\theta_0 = \mathbb{E}_0[Y(1)] - \mathbb{E}_0[Y(0)]$. Identification rests on the surrogate validity assumption of athey2025surrogate, which requires the conditional outcome mean to agree across the two populations, $\mathbb{E}_2[Y \mid X, S] = \mathbb{E}_1[Y \mid X, S]$, where the measures $\mathbb{E}_1$ and $\mathbb{E}_2$ in this case denote the experimental and observational samples, and $\mathbb{E}_0$ is the target population expectation. The nuisances are then
The stage functionals are evaluations: $m_2(W; g) = g(X, S)$ reads $g$ at the realised surrogate, with no counterfactual substitution, so the recursion content of the inner stage is the shift across samples and not a shift of the treatment. The outer functional is a contrast, $m_1(W; g) = g(X, 1) - g(X, 0)$. The two-sample structure means the inner and outer regressions are estimated on different datasets.
Two additional assumptions identify $\theta_0$. Common covariates, $p_1(X) = p_2(X)$, allows the surrogate distribution to be shifted across samples. Experimental unconfoundedness, $S(d) \perp D \mid X$ under measure $\mathbb{E}_1$, identifies the surrogate distribution under each arm.
\paragraph{Example 4: Mediation analysis.} We observe $n$ i.i.d.\ copies of $W = (X, D, M, Y)$, with covariates $X$, binary treatment $D$, mediator $M$, and outcome $Y$. The estimand is the cross-world mean $\theta_0 = \mathbb{E}[Y(d,\, M(1{-}d))]$, the mean outcome when the treatment is set to $d$ but the mediator is drawn as under treatment $1{-}d$ robins1992identifiability, pearl2001direct. Under a sequential ignorability assumption imai2010mediation, $\theta_0$ is identified by the mediation formula, a $T = 2$ recursive functional:
where $m_2(W; g) = g(X, d, M)$ evaluates $g$ at treatment $d$, and $\theta_0 = \mathbb{E}[f_1(X, 1{-}d)]$, so $m_1(W; g) = g(X, 1{-}d)$. The same formula arises in the alternative framework of interventional counterfactuals robins2010alternative. The two stages evaluate at different treatments. This sets mediation apart from Example 2, where one counterfactual sequence $(d_1, d_2)$ enters at both stages.
To make the algebra explicit, similar to bruns2025augmented, we work with estimators that are linear in a dictionary of basis functions. At each stage $t = 1, \ldots, T$, let $\phi_t : \mathcal{Z}_t \to \mathbb{R}^{k_t}$ be such a dictionary with each component square-integrable. Further, define $\phi_t^d$ as the image of the dictionary under the stage functional: $(\phi_t^d)_j := m_t(W;\, \phi_{t,j})$ for $j = 1, \ldots, k_t$. We take both nuisances to be linear in the dictionary, $f_t(Z_t) = \phi_t' \beta_t$ and $\alpha_t(Z_t) = \phi_t' \eta_t$, with coefficient vectors $\beta_t, \eta_t \in \mathbb{R}^{k_t}$. Since $m_t$ is linear in its functional argument, the definition gives the representation
so applying $m_t$ to a linear function swaps $\phi_t$ for $\phi_t^d$ and leaves the coefficients alone. The superscript $d$ records the leading case, in which $m_t$ substitutes a counterfactual treatment (or treatment path) and $\phi_t^d$ is simply $\phi_t$ with the treatment overwritten.
When the estimand is a contrast, as in the surrogacy estimand $\mathbb{E}_0[Y(1)] - \mathbb{E}_0[Y(0)]$, the outer functional is the combination $m_1(W; g) = g(X, 1) - g(X, 0)$, so $\phi_1^d = \phi_1^1 - \phi_1^0$, a difference of evaluations instead of a single one. The surrogacy inner stage does not substitute anything. Instead, $m_2(W; g) = g(X, S)$ reads $g$ at the realised surrogate, so $\phi_2^d = \phi_2$, and this stage shifts from the observational sample to the experimental sample.
Many nonparametric estimators have the property of linearity in $\phi_t$, e.g. series, lasso, and kernel methods, as well as certain neural networks and random forests. See bruns2025augmented for discussion.
We use the sample analogue of the time-varying $\mathbb{E}_t$ of Section (ref), denoted by $\hat{\mathbb{E}}_t$, which averages over the sample drawn from the stage-$t$ population. When there is only one measure, the single operator $\hat{\mathbb{E}}$ is used. Write $\hat{G}_t := \hat{\mathbb{E}}_t[\phi_t \phi_t']$ for the sample Gram matrix at stage $t$, and $\hat{M}_t := \hat{\mathbb{E}}_t[\phi_t\, (\phi_{t+1}^d)']$, $t = 1, \ldots, T-1$, for the cross-moment matrix linking adjacent stages, both under the outer stage-$t$ measure. Finally, at each stage $t$, let $P_t : \mathbb{R}^{k_t} \to \mathbb{R}$ be a convex penalty on the outcome coefficients, with parameter $\lambda_t \geq 0$, and $Q_t : \mathbb{R}^{k_t} \to \mathbb{R}$ a convex penalty on the Riesz coefficients, with parameter $\delta_t \geq 0$. In the notation of bruns2025augmented, $\hat{G}_1$ is their $\frac{1}{n}\Phi_p^\top\Phi_p$ and $\hat{\mathbb{E}}_0[\phi_1^d]$ their $\bar\Phi_q$.
The simplest estimator of $\theta_0$ is the plug-in estimator, which estimates each $f_t$ by an outcome regression and substitutes it into the identification formula. Each stage solves a regularised least-squares problem. The base case $t = T$ regresses $Y$ on the dictionary $\phi_T$:
where $P_T$ is convex and $\lambda_T \geq 0$ the regularisation parameter. At each stage $t = T{-}1, \ldots, 1$, the inner coefficient $\hat\beta_{t+1}$ produces the generated outcome $(\phi_{t+1}^d)'\hat\beta_{t+1}$, which is the target for stage $t$:
The first-order conditions of (ref)--(ref) define the backward pass:
The plug-in estimator is then \[ \hat\theta^{P} := \hat{\mathbb{E}}_0\!\left[m_1(W;\, \hat{f}_1)\right] = \hat{\mathbb{E}}_0\!\left[(\phi_1^d)'\right]\hat\beta_1. \] Here, the superscript $P$ records the outcome-regression penalty path $P = (P_1, \ldots, P_T)$ that produces $\hat\beta_1$.
The dual strategy estimates the representers $\alpha_t$ and reweights with them, instead of modelling the outcome functions. Each stage solves a regularised Riesz loss. The base case is stage $1$, which is a non-recursive problem expressed by chernozhukov2022automatic,singh2021kernel,chernozhukov2021riesz,chernozhukov2023covshifts as:
where $Q_1$ is the convex penalty and $\delta_1 \geq 0$ the regularisation parameter. At each stage $t = 2, \ldots, T$, the inner coefficient $\hat\eta_{t-1}$ produces the generated regressor $\hat\alpha_{t-1} = \phi_{t-1}'\hat\eta_{t-1}$, which enters the stage-$t$ Riesz loss for the functional $L_t(g) = \hat{\mathbb{E}}_{t-1}[\hat\alpha_{t-1}(Z_{t-1})\, m_t(W; g)]$ of chernozhukov2022nested:
The first-order conditions of (ref)--(ref) define the forward pass:
The balancing weight estimator is then \[ \hat\theta^{Q} := \hat{\mathbb{E}}_T\!\left[\hat\alpha_T \cdot Y\right] = \hat{\mathbb{E}}_T\!\left[Y\phi_T'\right]\hat\eta_T, \] the inverse-probability-weighted form of the Riesz regression estimator. Here, the superscript $Q$ denotes the Riesz penalty path $Q = (Q_1, \ldots, Q_T)$ that produces $\hat\eta_T$.
The recursive Riesz regression of the previous subsection has an equivalent reading as a recursive balancing procedure, in which $\hat\alpha_t = \phi_t'\hat\eta_t$ assigns unit $i$ at stage $t$ the weight $w_{i,t} = \phi_t(z_{t,i})'\hat\eta_t$. The reweighted stage-$t$ feature profile is then $\tfrac{1}{n}\sum_{i=1}^n w_{i,t}\,\phi_t(z_{t,i}) = \hat{G}_t\hat\eta_t$, and the forward-pass problems that produce $\hat\eta_t$ trade the variance of these weights against their imbalance relative to a target profile. This is the structure that bruns2025augmented, building on zubizarreta2015stable, athey2018residual, and hirshberg2021augmented, and many others, analyse for $T=1$. The change from one step to multiple steps is that the non-recursive target $\bar\phi^d = \hat{\mathbb{E}}[\phi^d]$ is replaced by a generated target $\hat\tau_t$.
\paragraph{Equivalence of recursive balancing and recursive Riesz regression.} The $T=1$ balancing-Riesz equivalence carries over to the recursive setting stage by stage. At each stage $t$, define the recursive target vector \[ \hat\tau_t :=
\] which replaces the simple target $\bar\phi^d$ of the non-recursive case.
Stage by stage, then, the forward-pass recursive Riesz regression can be viewed as a sequential balancing procedure. The constrained, penalised, and Riesz formulations of $\hat\eta_t$ all return the same weights.
\paragraph{Connection to dynamic covariate balancing.} viviano2021dynamic propose dynamic covariate balancing (DCB) for the time-varying treatment setting where at each period $t$, they solve a constrained programme minimising the $\ell_2$ norm of the weights subject to the dynamic balance condition
where $H_{i,t}$ is the observed history of unit $i$ at period $t$, $\gamma_{i,t}$ is the period-$t$ weight on unit $i$ being chosen, $\hat\gamma_{i,t-1}$ is the weight estimated at period $t-1$, the tolerance $\kappa_t$ bounds the imbalance, and $\hat\gamma_{i,0} = 1/n$ initialises the recursion. DCB also restricts each weight vector to the simplex (non-negative weights summing to one), caps the largest weight, and restricts the weights to the units on the target path, $\gamma_{i,t} = 0$ unless $D_{i,1:t} = d_{1:t}$. On the retained units the observed history equals the counterfactual one, which is what lets the same $H_{i,t}$ stand in for $\phi_t^d$ on both sides of (ref).
Under the identification $\hat\alpha_t(W_i) = n\,\gamma_{i,t}$ (the DCB weights sum to one while the representer averages to one), condition (ref) bounds the same imbalance as (ref) with $\|\cdot\|_* = \|\cdot\|_\infty$ and $\delta_t^{(C)} = \kappa_t$. The two programmes minimise the weight variance over different sets, form (ref) over the span of the dictionary on the full sample and DCB over free weight vectors that are zero off the target path, so for a generic dictionary the solutions differ even at exact balance. Without the simplex constraint and the weight cap, the two coincide when the stage-$t$ dictionary is saturated in the treatment path, meaning every entry either is $\mathbf{1}\{D_{1:t} = d_{1:t}\}$ times a function of the covariate history or vanishes on the path.
The doubly robust estimator takes the sample analogue of the orthogonal moment (ref), evaluated at the estimated nuisances:
where $\hat\varepsilon_t$ are the sample residuals of Section (ref): $\hat\varepsilon_t = (\phi_{t+1}^d)'\hat\beta_{t+1} - \phi_t'\hat\beta_t$ for $t < T$, by (ref), and $\hat\varepsilon_T = Y - \phi_T'\hat\beta_T$. In the linear framework, $\hat\theta^{DR}$ is estimated using the backward-pass coefficients $\hat\beta_t$ from (ref)--(ref) and the forward-pass coefficients $\hat\eta_t$ from (ref)--(ref).
When we estimate all $2T$ nuisance functions by OLS, with no penalty at any stage, the recursive plug-in estimator $\hat\theta^{P}$, the recursive balancing weight estimator $\hat\theta^{Q}$, and the doubly robust estimator $\hat\theta^{DR}$ coincide exactly. This is an algebraic identity, so it holds in any sample in which each Gram matrix $\hat{G}_t$ is invertible, and whether or not the linear models are correctly specified.
Throughout this section each Gram matrix $\hat{G}_t$ is invertible, that is, the stage-$t$ dictionary has full column rank in the sample, which requires $k_t \leq n_t$. OLS falls under our linear estimators where the regularisation penalties are set to zero, ($P_t = Q_t = 0$ for all $t$), such that the backward and forward passes (ref)--(ref) lead to FOCs:
and
The recursive plug-in and balancing weight estimators of Section (ref) specialise to:
The next proposition records the equivalence and the two facts behind it: every correction is individually zero, and every stage-wise bilinear form already equals the target.
Part (i) does not use OLS on the Riesz side. The proof uses only the outcome-side normal equations and the linearity of the weights, so the corrections vanish for any weights in the span of the regressors. This is a doubly robust property of OLS: adding any linear balancing weights to an OLS outcome model leaves the estimator unchanged. Equivalently, one can show that adding a linear outcome regression to OLS balancing weights does not change the estimator.
With ridge Riesz regressions and any outcome model, the doubly robust estimator can be represented as a single recursive outcome regression, but with augmented coefficients set by a backward recursion. Each stage shrinks towards OLS, but towards an OLS prediction of an already-augmented stage. The weight the estimator leaves on the pure OLS plug-in estimator decays geometrically in the number of time periods.
Each stage contains two penalties: $\lambda_t$ on the outcome regression and $\delta_t$ on the Riesz regression, each a positive definite diagonal matrix. Define \[ A_t := (\hat{G}_t + \delta_t)^{-1}\hat{G}_t, \qquad I - A_t = (\hat{G}_t + \delta_t)^{-1}\delta_t, \] where the second equation follows from $(\hat{G}_t + \delta_t) - \hat{G}_t = \delta_t$. The matrix $A_t$ is the ridge shrinkage operator of the time $t$ Riesz regression. For a stage-$t$ least-squares problem with a fixed target $b$, ridge predicts $\hat\beta^{R} = A_t \hat\beta^{OLS}$, so ridge is $A_t$ applied to the OLS coefficient, and $I - A_t$ measures the shrinkage between them.
For a scalar penalty $\delta_t = \delta I$ the two matrices commute, $A_t$ is symmetric, and $0 \prec A_t \prec I$. For a general diagonal $\delta_t$ the matrix $A_t$ need not be symmetric, but it is similar to the symmetric positive definite $(\hat{G}_t + \delta_t)^{-1/2}\hat{G}_t(\hat{G}_t + \delta_t)^{-1/2}$, so its eigenvalues are real and lie in $(0, 1)$ for any $\delta_t \succ 0$. Eigenvalues in $(0,1)$ are what make $A_t$ a contraction away from OLS.
The ridge Riesz first-order conditions generalise (ref) to:
with $\hat\alpha_t^R = \phi_t'\hat\eta_t^R$. Similarly, the ridge outcome regressions generalise (ref):
Here, the balancing weights estimator built from the ridge Riesz regression still equals the plug-in regression built from the ridge outcome pass, as long as both use the same penalty matrices, $\delta_t = \lambda_t$ at every stage. The estimators are not equal to the doubly robust estimator anymore, as was the case for OLS.
This equivalence is specific to ridge with equal penalties. At $T=1$, this equivalence reduces to kallus2020generalized,hirshberg2019minimaxlinear.
We now characterise the doubly robust estimator when the Riesz representers are estimated by ridge but the outcome regressions are estimated with a general linear estimator.
At time period $T$, the augmented coefficient is an affine blend of the penalised regression $\hat\beta_T^{Gen}$ and the final stage $T$ OLS regression $\hat\beta_T^{OLS}$, weighted by $I - A_T$ and $A_T$. This is equivalent to the single-regression result of bruns2025augmented. Each earlier stage has the same form but with a different target: it shrinks towards $\hat\beta_t^{Aug\text{-}OLS}$ instead of $\hat\beta_t^{OLS}$. Here, $\hat\beta_t^{Aug\text{-}OLS}$ are the coefficients of an OLS regression of the augmented generated outcomes from the stage before.
For the high-dimensional regime $k_t > n_t$ where $\hat{G}_t$ is not invertible, there the OLS targets $\hat\beta_t^{OLS}$ and $\hat\beta_t^{Aug\text{-}OLS}$ are read as minimum-norm solutions, with $\hat{G}_t^{-1}$ replaced by the Moore--Penrose pseudoinverse. Appendix (ref) shows equivalent results where the dictionaries are infinite-dimensional kernel features.
\paragraph{Comparison with bruns2025augmented.} For $T = 1$ the recursion collapses to a single equation:
recovering Proposition 3.2 of bruns2025augmented exactly. For $T \geq 2$, the last stage $T$ satisfies (ref) directly, while earlier stages have the augmented generated outcome as target, and therefore the final augmented plug-in estimator will have shrunk towards something else than recursive OLS.
\paragraph{Geometric decay of OLS weight.} The number of time periods determines how much weight on the pure OLS plug-in coefficients is retained. Take the scalar diagonal case where $\hat{G}_t = \sigma_t^2 I$, then $A_t = a_t I$ with $a_t = \sigma_t^2/(\sigma_t^2+\delta_t) \in (0,1)$. Unrolling the recursion gives:
and so on until stage 1. Every term but the last contains some regularisation. The one term without the general regression coefficients, which we call the OLS anchor, is $\bigl(\prod_{t=1}^T a_t\bigr)\hat\beta_1^{OLS}$. Note that this recursively shrunk OLS anchor is equal to a ridge plug-in estimator with the penalties from the Riesz regressions, which is $\hat\beta_1^{R}$ under penalties $\delta_t = \lambda_t$.
\paragraph{Shrinkage in the non-diagonal case.} For non-diagonal Gram matrices, we derive a contraction in the Euclidean norm. For a scalar penalty $\delta_t I > 0$ the matrix $A_t$ is symmetric, with spectral norm $\rho_t := \|A_t\|_2 = \max_j \mu_{t,j}/(\mu_{t,j}+\delta_t) < 1$, where $\mu_{t,j}$ are eigenvalues of $\hat{G}_t$.
For (non-scalar) diagonal penalty matrices, $A_t$ need not be symmetric and its Euclidean norm can exceed one even though its eigenvalues stay in $(0,1)$. However, the contraction then survives stage by stage in the $(\hat{G}_t+\delta_t)$-weighted norms, with $\rho_t$ replaced by the largest eigenvalue of $A_t$.
We now consider the doubly robust estimator where the outcome regression is also estimated with ridge penalties, and a shared penalty matrix at every stage, so $\delta_t = \lambda_t$.
The final stage ridge outcome coefficients can be written as $\hat\beta_T^R = A_T\hat\beta_T^{OLS}$, and substituting into (ref) gives \[ \hat\beta_T^{Aug} = \bigl[I - (I - A_T)^2\bigr]\hat\beta_T^{OLS}, \] so the augmented coefficients square the shrinkage towards OLS. When $\lambda_T$ and $\hat{G}_T$ commute, in particular for a scalar penalty $\lambda_T = \lambda I$, the squared gap is itself a ridge shrinkage. Like the undersmoothing characterisation of bruns2025augmented, $\hat\beta_T^{Aug}$ is a ridge estimator with the reduced effective penalty
which in the scalar case $\hat{G}_T = \sigma_T^2 I$ is $\Gamma_T = \lambda_T^2/(\sigma_T^2 + 2\lambda_T)$.
Each sequential stage then shrinks the outcome regression coefficients $\hat\beta_t^R$ not towards the OLS plug-in coefficients, but towards an OLS regression on the augmented generated outcome from the previous stage.
Finally, we consider $P_t$ and $Q_t$ general convex penalties. The stage coefficients have no closed form. We derive a telescoping identity relating the balancing weight estimator to OLS, and a penalty-gradient decomposition relating the doubly robust estimator to OLS.
Recall the Riesz target, the right-hand side of the forward pass (ref):
For any coefficient vectors $\hat\eta_1, \ldots, \hat\eta_T$, define the Riesz residual at stage $t$ as
When $\hat\eta_t$ solves the stage-$t$ Riesz regression (ref) with convex penalty $Q_t$ and parameter $\delta_t \geq 0$, convex optimality is the subgradient inclusion
so $c_t = \tfrac{\delta_t}{2}\hat{s}_t$ for the subgradient $\hat{s}_t \in \partial Q_t(\hat\eta_t)$ the solution selects, and $c_t = \tfrac{\delta_t}{2}\nabla Q_t(\hat\eta_t)$ when $Q_t$ is differentiable. Under OLS Riesz regressions ($\delta_t = 0$) there is exact balance, $\hat{G}_t\hat\eta_t = \hat\tau_t$, so $c_t = 0$ at every stage. Under ridge we obtain $c_t = \delta_t\hat\eta_t^R$.
The first identity writes the IPW estimator as OLS corrected by the Riesz residuals. At $T = 1$, bruns2025augmented show that any linear balancing weights estimate the same quantity as a regression on a shifted feature profile, $\hat\theta^{Q} = (\hat{\mathbb{E}}_0[\phi_1^d] - c_1)'\hat\beta_1^{OLS}$. The regularisation shifts the target $\hat{\mathbb{E}}_0[\phi_1^d]$ by $c_1$.
Here, we decompose the doubly robust estimator by combining Lemma (ref) with a similar telescoping identity for $\hat\theta^{DR}$ to $\hat\theta^{Q}$.
This result shows why the doubly robust estimator shrinks less to OLS as the number of time periods increases.
This paper takes the dissection of bruns2025augmented from a single penalised regression into the recursive setting of chernozhukov2022nested, in which time-varying treatment effects, surrogate identified treatment effects, and mediation effects use $T$ stages of regression.
The analysis characterises the resulting estimators under three forms of regularisation. Under OLS at every stage the plug-in, balancing weight, and doubly robust estimators coincide exactly. Under ridge penalties, the weight on the pure-OLS fit decays geometrically in the number of time periods. Under a general convex penalty the closed forms break down, but each stage still contributes one correction term, an inner product between the stage Riesz residual and the gap between the penalised and OLS outcome coefficients.
There are two takeaways. First, the main equivalence from cross sectional causal inference extends to longitudinal causal inference. Second, the interpretation of debiasing as undersmoothing weakens.
\spacingset{1.5}