Triple/Double-Debiased Lasso

Triple/Double-Debiased Lasso

\affil[1]{{University of California, Los Angeles}} \affil[2]{{University of Copenhagen}} \affil[3]{{Yale University}} \affil[4]{{Aarhus Center for Econometrics (ACE)}}

Introduction

Inference on a low-dimensional target parameter in the presence of many controls has become a central problem in econometrics because making a conditional exogeneity assumption plausible in empirical applications often requires using a rich set of covariates BCFH17. In such settings, in addition to the target parameter, we also have a high-dimensional nuisance parameter that has to be estimated via machine learning methods such as the Lasso estimator. Unfortunately, regularization underlying these machine learning methods generates non-regular first-stage bias that typically invalidates naive plug-in inference procedures for the target parameter. A major insight of the double/debiased machine learning literature is that this problem can often be overcome by replacing the original moment equation with the one based on a moment function whose first derivative with respect to the nuisance parameter vanishes at the true parameter values (first-order Neyman orthogonality condition), so that sufficiently accurate regularized estimators of the nuisance parameter matter only through second-order remainder terms and the target parameter can be estimated at the usual $1/\sqrt{n}$ rate CHS15. In this paper we show that, in the high-dimensional linear regression model, one can go one step further and construct an estimator that removes not only the first-order effect of nuisance parameter estimation but also the second-order effect. Because of the steps in the construction of our estimator, we refer to it as a triple Lasso estimator, or, equivalently, a double-debiased Lasso estimator.

Our starting point is the standard high-dimensional linear regression model in which the outcome is regressed on a scalar treatment variable of interest and a potentially high-dimensional vector of controls. In this model, the conventional double Lasso estimator BCH14 is based on a moment function that satisfies the first-order Neyman orthogonality condition and is often remarkably successful even when the underlying Lasso estimators are far from being $\sqrt{n}$-consistent. However, the moment function behind the double Lasso estimator does not satisfy the {\em second-order} Neyman orthogonality condition, and this means that the double Lasso estimator may perform poorly if second-order remainder terms are large enough to matter in empirically relevant samples, which happens when the first-stage Lasso estimators converge rather slowly. Our first main contribution is thus to construct a moment function that satisfies not only the first-order but also the second-order Neyman orthogonality condition. Our moment function is obtained from the double Lasso moment function by subtracting a carefully chosen quadratic correction term, weighted by the inverse Gram matrix of the controls, so that the second-order bias term is cancelled by construction. In turn, our second main contribution is to propose the triple Lasso estimator that solves a moment condition based on this moment function, with all underlying nuisance parameters being estimated by the Lasso-type methods.

We derive an asymptotic linear representation for our triple Lasso estimator and compare its remainder terms with those of the usual double Lasso estimator. We show that whenever the double Lasso estimator is already asymptotically normal, the remainder terms for the triple Lasso estimator are never larger and are often strictly smaller in order than the corresponding terms for the double Lasso estimator. As a result, the distribution of the triple Lasso estimator is often closer to the properly centered normal distribution than that of the double Lasso estimator. Moreover, in the special case of exact sparsity and consistent Lasso screening, the triple Lasso estimator can remain asymptotically normal along sequences of data-generating processes for which the double Lasso estimator fails to be asymptotically normal.

Our Monte Carlo simulations confirm that the implications of the asymptotic theory remain relevant in empirically meaningful finite-sample settings. Across designs that vary the sample size, the correlation structure of the controls, and the sparsity of the nuisance parameters, although the variance of the triple Lasso estimator is typically somewhat larger than that of the double Lasso estimator, the former estimator often has substantially smaller squared bias and the mean squared error, especially in the harder designs. More importantly, however, the triple Lasso estimator often leads to more precise inference. For example, in one of the designs, the empirical coverage of nominal $95\%$ confidence intervals rises from about $67\%$ under the double Lasso estimator to about $92\%$ under the triple Lasso estimator. This happens despite the fact that the triple Lasso confidence intervals are only slightly longer on average.

As our third main contribution, we extend the construction of the second-order Neyman orthogonal moment function to a general Z-estimation problem. In fact, we provide a general recursive formula that can be used to construct moment functions in the Z-estimation problem that satisfy the Neyman orthogonality condition to {\em any} desired order. We further illustrate this construction in an M-estimation problem with a single-index structure, delivering the results for the high-dimensional linear regression model as a special case.

Our paper is related to several strands of the literature. First, it builds directly on the foundational work on double/debiased Lasso and post-selection inference in high-dimensional linear regression models, including BCH14, JM14, GBRD14, and ZZ14, all of which established that valid inference on low-dimensional components is possible after regularization when the moment function is appropriately debiased. Closely related is the double machine learning literature of CHS15 and CCDDHNR18, from which we adopt the emphasis on orthogonal moment functions, while extending that perspective from first-order to higher-order orthogonality. An alternative approach to inference in high-dimensional linear regression models is developed by AKK20, who propose explicitly accounting for the bias induced by regularization rather than eliminating it through orthogonalization; an advantage of this approach is that it achieves validity under weaker structural conditions, although it requires the researcher to specify a bound on the magnitude of the nuisance coefficients. Second, our paper is connected to the broader literature on approximately sparse high-dimensional models, including BC11, belloni_sparse_2012, and G16, from which we borrow technical tools, in particular the convergence rates for Lasso estimators in linear regression models and for node-wise Lasso estimators of inverse Gram matrices. Third, our paper is related to the emerging literature on higher-order orthogonality, including MSZ18 and BJW24. In particular, MSZ18 construct higher-order Neyman orthogonal moment functions for the partially linear model under {\em non}-Gaussian first-stage errors, while BJW24 develop a general algorithm for constructing higher-order Neyman orthogonal moment functions in likelihood-based models. These results are not directly applicable in our setting, as we focus on regression-based estimation rather than likelihood models and do not restrict the distribution of noise. Also, our results do not contradict the negative result in MSZ18 on the existence of higher-order Neyman orthogonal moment functions in the partially linear model with Gaussian first-stage errors, because our moment function depends not only on a vector of observable random variables but also on an independent copy of this vector, and, moreover, our orthogonality requirement is formulated in terms of ordinary derivatives rather than the stronger functional-derivative notion studied there, reflecting the parametric nature of our regression model. Finally, our paper is also related to earlier work on higher-order influence functions and U-statistic-based estimators in semiparametric and nonparametric models, including RLTV08, RLMTV17, NR18, and LMNR23, which develop higher-order expansions of functionals and corresponding estimators that balance bias and variance in settings where first-order influence function methods are insufficient. This literature shows that higher-order corrections can be used to relax smoothness requirements and/or improve rates in semiparametric and nonparametric models.

The rest of the paper is organized as follows. Section (ref) introduces the triple Lasso moment function and the corresponding estimator in the high-dimensional linear regression model, proves second-order Neyman orthogonality of the proposed moment function, and derives the asymptotic linear representation and asymptotic normality result for the triple Lasso estimator. Section (ref) develops the general recursive formula for constructing higher-order Neyman orthogonal moment functions for the Z-estimation problem and applies it to M-estimators with a single-index structure, thereby showing that the logic behind the triple Lasso estimator is a part of a broader higher-order orthogonalization principle. Section (ref) reports Monte Carlo results showing that the triple Lasso estimator can substantially reduce bias and deliver more reliable inference in comparison with the double Lasso estimator under certain data-generating processes.

Notation

For any $p\in\mathbb N$, we denote $[p]:=\{1,\dots,p\}$ and $\boldsymbol{0}_p := (0,\dots,0)^\top \in\mathbb{R}^p$. Also, we let $\boldsymbol{0}_{p\times p}$ and $\mathbf{I}_p$ be the zero and identity matrices in $\mathbb{R}^{p\times p}$, respectively. In addition, for any finite set of positive integers $I$, we use $\mathrm{E}_I[f(\boldsymbol{X}_i)]$ to denote the average value of $f(\boldsymbol{X}_i)$ as $i$ varies over $I$, namely $\mathrm{E}_I[f(\boldsymbol{X}_i)] = |I|^{-1}\sum_{i\in I}f(\boldsymbol{X}_i)$. For numbers $a_n$ and $b_n$, $n\in\mathbb N$, we write $a_n\lesssim b_n$ if there exists a constant $C>0$ such that $|a_n|\leq C|b_n|$ for all $n\in \mathbb N$. For random variables $V_n$ and $R_n$, $n\in\mathbb N$, we write $V_n\lesssim_P R_n$ if $V_n = Y_n R_n$ for some $Y_n = O_P(1)$. Finally, for any matrix $\mathbf{A} = (A_{j,k})_{j,k=1}^p$, we denote $\|\mathbf{A}\|_{\mathrm{max}} := \mathrm{max}_{1\leq j,k\leq p}|A_{j,k}|$ and $\|\mathbf{A}\|_{\infty} := \mathrm{max}_{1\leq i\leq p}\sum_{j=1}^p |A_{j,k}|$.

Triple Lasso Estimator

In this section, we propose a moment function for the high-dimensional linear regression model that satisfies the second-order Neyman orthogonality condition and introduce the corresponding triple Lasso estimator. In addition, we derive the $\sqrt n$-consistency and asymptotic normality result for our triple Lasso estimator and compare the remainder terms of the corresponding asymptotic linear representation to those underlying the asymptotic linear representation for the double Lasso estimator.

Model, Second-Order Neyman Orthogonality, and Estimation

Consider the linear regression model

where $Y\in\mathbb{R}$ is the outcome, $D\in\mathbb{R}$ is a regressor of interest, $\boldsymbol{X} = (X_1,\dots,X_p)^\top\in\mathbb{R}^p$ is a vector of controls, $\beta_0\in\mathbb{R}$ is the parameter of interest, and $\boldsymbol{\theta}_0 = (\theta_{0,1},\dots,\theta_{0,p})^\top\in\mathbb{R}^p$ is a vector of nuisance parameters. Let $(\boldsymbol{X}_i,D_i,Y_i)$, $i=1,\dots,n$, be a random sample from the distribution of $(\boldsymbol{X},D,Y)$, where we denote $\boldsymbol{X}_i = (X_{i,1},\dots,X_{i, p})^\top$ for all $i\in[n]$. In this paper, we primarily focus on the high-dimensional setting in which the number of controls $p$ may be comparable to or larger than the sample size $n$.

A conventional way to estimate $\beta_0$ in model (ref) when $p$ is large is via the double Lasso estimator.\footnote{Various versions of the double/debiased Lasso estimator appeared in BCH14, JM14, GBRD14, and ZZ14.} To describe it, let the projection of $D$ on $\boldsymbol{X}$ be

where $\boldsymbol{\gamma}_0 = (\gamma_{0,1},\dots,\gamma_{0,p})^\top\in\mathbb{R}^p$ is a vector of parameters, and let the reduced-form regression of $Y$ on $\boldsymbol{X}$ be

where $\boldsymbol{\phi}_0 = (\phi_{0,1},\dots,\phi_{0,p})^\top := \boldsymbol{\theta}_0 + \beta_0 \boldsymbol{\gamma}_0$ and $e := \varepsilon + \beta_0\nu$. One version of the double Lasso estimator is constructed as follows. First, estimate $\boldsymbol{\phi}_0$ by running the Lasso regression of $Y$ on $\boldsymbol{X}$, yielding $\widehat\boldsymbol{\phi}$. Second, estimate $\boldsymbol{\gamma}_0$ by running the Lasso regression of $D$ on $\boldsymbol{X}$, yielding $\widehat\boldsymbol{\gamma}$. Third, estimate $\beta_0$ by running the OLS regression of $Y-\boldsymbol{X}^\top\widehat\boldsymbol{\phi}$ on $D-\boldsymbol{X}^\top\widehat\boldsymbol{\gamma}$, yielding the double Lasso estimator $\widehat\beta^{DL}$. Under certain conditions, the double Lasso estimator is $\sqrt n$-consistent and asymptotically normal, despite the fact that the underlying Lasso estimators $\widehat\boldsymbol{\phi}$ and $\widehat\boldsymbol{\gamma}$ are not $\sqrt n$-consistent.

To understand why this happens, introduce the moment function $\psi^{DL}\colon \mathbb{R}\times \mathbb{R}^p\times \mathbb{R}^p$ by \[ \psi^{DL}(\beta, \boldsymbol{\gamma}, \boldsymbol{\phi}) := \mathrm{E}[(Y - \boldsymbol{X}^\top\boldsymbol{\phi} - \beta(D - \boldsymbol{X}^\top\boldsymbol{\gamma}))(D - \boldsymbol{X}^\top\boldsymbol{\gamma})] \] and observe that the true value $\beta_0$ satisfies the moment equation

The double Lasso estimator $\widehat\beta^{DL}$ solves an empirical analogue of this equation obtained by replacing the population expectation with the empirical expectation and the true values $\boldsymbol{\gamma}_0$ and $\boldsymbol{\phi}_0$ with their estimators $\widehat\boldsymbol{\gamma}$ and $\widehat\boldsymbol{\phi}$. The key fact underlying the $\sqrt n$-consistency and asymptotic normality of the double Lasso estimator is that there is no first-order effect of replacing $\boldsymbol{\gamma}_0$ and $\boldsymbol{\phi}_0$ by $\widehat\boldsymbol{\gamma}$ and $\widehat\boldsymbol{\phi}$ on the moment function $\psi^{DL}$:

Indeed, provided that the estimators $\widehat\boldsymbol{\gamma}$ and $\widehat\boldsymbol{\phi}$ are sufficiently precise so that the second-order effects are asymptotically negligible, condition (ref) ensures that the estimator $\widehat\beta^{DL}$ is asymptotically equivalent to the infeasible estimator $\widetilde\beta^{DL}$ obtained by solving the empirical analogue of (ref) in which only the population expectation is replaced by the empirical expectation. The latter estimator is in turn $\sqrt n$-consistent and asymptotically normal by the classic $M$-estimation theory.

Because of (ref), the moment function $\psi^{DL}$ is said to satisfy the (first-order) {\em Neyman orthogonality} condition. To express this condition more compactly, let $\boldsymbol{\eta} := (\boldsymbol{\gamma}^\top,\boldsymbol{\phi}^\top)^\top$ and $\boldsymbol{\eta}_0 := (\boldsymbol{\gamma}_0^\top,\boldsymbol{\phi}_0^\top)^\top$. Then (ref) can be written equivalently as \[ \frac{\partial \psi^{DL}(\beta_0,\boldsymbol{\eta}_0) }{\partial \boldsymbol{\eta}}= \boldsymbol{0}_{2p}. \] Similarly, one can define the second-order Neyman orthogonality. Following MSZ18 and BJW24, we say that a moment function $\psi\colon\mathbb{R}\times\mathbb{R}^q \to\mathbb{R}$ satisfies the {\em second-order Neyman orthogonality} condition if \[ \frac{\partial \psi(\beta_0,\boldsymbol{\eta}_0)}{\partial\boldsymbol{\eta}} =\boldsymbol{0}_q \quad\text{and}\quad \frac{\partial^2 \psi(\beta_0,\boldsymbol{\eta}_0)}{\partial\boldsymbol{\eta}\partial\boldsymbol{\eta}^\top}=\boldsymbol{0}_{q\times q}, \] where $\boldsymbol{\eta}_0\in\mathbb{R}^q$ is a vector of nuisance parameters. Intuitively, using estimators based on moment functions satisfying the second-order Neyman orthogonality condition is beneficial because such moment functions eliminate not only the first-order effects of replacing the true values of the nuisance parameters by the corresponding estimators but also the second-order effects. Unfortunately, however, the moment function $\psi^{DL}$ underlying the double Lasso estimator does not have this property. Indeed, it is straightforward to verify that \[ \frac{\partial^2\psi^{DL}(\beta_0,\boldsymbol{\gamma}_0,\boldsymbol{\phi}_0)}{\partial\boldsymbol{\gamma}\partial\boldsymbol{\phi}} = \mathrm{E}[\boldsymbol{X}\boldsymbol{X}^\top] \neq \boldsymbol{0}_{p\times p}. \] It therefore follows that convergence of the double Lasso estimator to the properly centered normal distribution may fail or be slow if the Lasso estimators $\widehat\boldsymbol{\gamma}$ and $\widehat\boldsymbol{\phi}$ are not sufficiently precise and the second-order effects are non-negligible.

Motivated by this observation, we propose a novel moment function for estimating $\beta_0$ in model (ref) that satisfies not only the first-order Neyman orthogonality condition but also the second-order Neyman orthogonality condition. We then define a {\em triple Lasso} estimator that solves an empirical analogue of the moment condition implied by this function and derive its asymptotic properties.

In order to introduce our moment function, let $\boldsymbol{\Sigma}_0 := \mathrm{E}[\boldsymbol{X}\boldsymbol{X}^\top]$ denote the population Gram matrix of the controls $\boldsymbol{X}$, and let $\boldsymbol{\Theta}_0 := \boldsymbol{\Sigma}_0^{-1}$ denote its inverse. Define the function $\psi^{ADJ} \colon \mathbb{R}\times \mathbb{R}^p\times \mathbb{R}^p \times \mathbb{R}^{p\times p}$ by \[ \psi^{ADJ}(\beta, \boldsymbol{\gamma}, \boldsymbol{\phi}, \boldsymbol{\Theta}) := \mathrm{E}[(Y - \boldsymbol{X}^\top\boldsymbol{\phi} - \beta(D - \boldsymbol{X}^\top\boldsymbol{\gamma}))\boldsymbol{X}^\top \boldsymbol{\Theta} \tilde\boldsymbol{X}(\tilde D-\tilde \boldsymbol{X}^\top\boldsymbol{\gamma})], \] where $(\tilde\boldsymbol{X},\tilde D)$ is a copy of $(\boldsymbol{X},D)$ that is independent of $(\boldsymbol{X},D,Y)$. Our moment function $\psi^{TL}\colon \mathbb{R}\times\mathbb{R}^{2p+p^2}\to\mathbb{R}$ for estimating $\beta_0$ in model (ref) is then defined as

In the following lemma, we show that the true value $\beta_0$ solves the moment equation

where $\boldsymbol{\eta}_0 := (\boldsymbol{\gamma}^\top_0,\boldsymbol{\phi}^\top_0,\mathrm{vec}(\boldsymbol{\Theta}_0)^\top)^\top$ and that the moment function $\psi^{TL}$ satisfies the second-order Neyman orthogonality condition. The proof of this lemma can be found in the Appendix.

We next define the triple Lasso estimator as a solution to an empirical analogue of equation (ref), where $\boldsymbol{\gamma}_0$ and $\boldsymbol{\phi}_0$ are replaced by the corresponding Lasso estimators and $\boldsymbol{\Theta}_0$ is replaced by the node-wise Lasso estimator. For the reasons that will become clear shortly, we in fact use row-sparsified version of the node-wise Lasso estimator of $\boldsymbol{\Theta}_0$. In addition, since it is useful for the asymptotic normality result, we also combine the estimator with cross-fitting, which is standard in the literature on double machine learning; see CCDDHNR18.

To formally describe the resulting estimator, split at random the full sample into $K$ subsamples of approximately the same size for some small $K$ and let $I(1),\dots,I(K)$ be the sets of indices from $1$ to $n$ corresponding to the observations in these subsamples. Also, for each $k\in[K]$, let $I(-k) := [n]\setminus I(k)$. Then for each $k\in[K]$, let $\widehat\boldsymbol{\gamma}_k = (\widehat\gamma_{k,1},\dots,\widehat\gamma_{k,p})^\top$, $\widehat\phi_k = (\widehat\phi_{k,1},\dots,\widehat\phi_{k,p})^\top$, and $\widetilde\boldsymbol{\Theta}_k = (\widetilde\Theta_{k, j, l})_{j,l\in[p]}$ be the Lasso estimator of $\boldsymbol{\gamma}_0$, the Lasso estimator of $\boldsymbol{\phi}_0$, and the node-wise Lasso estimator of $\boldsymbol{\Theta}_0$, respectively, all defined on the subsample of observations with indices in $I(-k)$. Here, we use the row-wise version of the node-wise Lasso estimator. Specifically, to define the first row $\widetilde\boldsymbol{\Theta}_{k,1}$ of the matrix $\widetilde\boldsymbol{\Theta}_k$, we run the Lasso regression of $X_1$ on $\boldsymbol{X}_{-1} := (X_2,\dots,X_p)^\top$ and obtain the vector of slope coefficients, say $\widehat\boldsymbol{\psi}_1$. Then we set $\widehat\sigma_1^2 := \mathrm{E}_{I(-k)}[X_{i,1}(X_{i,1} - \boldsymbol{X}_{i,-1}^\top\widehat\boldsymbol{\psi}_1)]$ and $\widetilde\boldsymbol{\Theta}_{k,1} := (1, - \widehat\boldsymbol{\psi}_1^\top)/\widehat\sigma_1^2$, where $\boldsymbol{X}_{i,-1}:=(X_{i,2},\dots,X_{i,p})^\top$. All other rows of $\widetilde\boldsymbol{\Theta}_k$ are defined analogously. Next, let $\widehat T_{k,0} := \{j=1,\dots,p\colon \widehat\gamma_{k,j} \neq 0\}$ be the set of indices from $1$ to $p$ corresponding to the non-zero components of the vector $\widehat\boldsymbol{\gamma}_k$ and let $\widehat\boldsymbol{\Theta}_{k}$ be the $\mathbb{R}^{p\times p}$ matrix obtained from $\widetilde\boldsymbol{\Theta}_k$ by setting to zero all its rows whose indices are not in $\widehat T_k:=\widehat T_{k,0} \cup \widehat T_{k,1}$, where $\widehat T_{k,1}$ is an additional set of indices corresponding to the extra variables that the researcher might want to use as we discuss below. Throughout most of the paper, we assume that $\widehat T_{k,1} = \emptyset$. To define $\widehat\beta_k$, we now solve an empirical version of the equation (ref) on $I_k$ for $\beta$, namely we set $$ \widehat\beta_k := \frac{\mathrm{E}_{I(k)}[(Y_i - \boldsymbol{X}_i^\top \widehat\boldsymbol{\phi}_k)(D_i - \boldsymbol{X}_i^\top\widehat\boldsymbol{\gamma}_k)] - \mathrm{E}_{I(k)}[(Y_i - \boldsymbol{X}_i^\top\widehat\boldsymbol{\phi}_k)\boldsymbol{X}_i^\top]\widehat\boldsymbol{\Theta}_{k}\mathrm{E}_{I(k)}[\boldsymbol{X}_i(D_i - \boldsymbol{X}_i^\top\widehat\boldsymbol{\gamma}_k)]}{\mathrm{E}_{I(k)}[(D_i - \boldsymbol{X}_i^\top\widehat\boldsymbol{\gamma}_k)^2] - \mathrm{E}_{I(k)}[(D_i - \boldsymbol{X}_i^\top\widehat\boldsymbol{\gamma}_k)\boldsymbol{X}_i^\top]\widehat\boldsymbol{\Theta}_{k} \mathrm{E}_{I(k)}[\boldsymbol{X}_i(D_i - \boldsymbol{X}_i^\top\widehat\boldsymbol{\gamma}_k)]}. $$ Finally, we aggregate the subsample estimators $\widehat\beta_k$ to obtain the triple (or double-debiased) Lasso estimator: $$ \widehat\beta^{TL} := \frac{1}{K}\sum_{k=1}^K \widehat\beta_k. $$ In the next subsection, we will derive the asymptotic theory for the triple Lasso estimator $\widehat\beta^{TL}$ and compare it with that for the double Lasso estimator $\widehat\beta^{DL}$.

Asymptotic Theory

In this subsection, we establish the $\sqrt n$-consistency and asymptotic normality result for the triple Lasso estimator $\widehat\beta^{TL}$. To this end, we impose several regularity conditions. These conditions control the moments of the disturbances $(\varepsilon,\nu)$, the sample splitting scheme used in the construction of the estimator, the behavior of the Lasso estimators appearing in the algorithm, and the sparsity of the model.

Assumptions (ref)--(ref) impose basic moment and regularity conditions on the disturbances, regressors, and the sample splitting scheme. Assumption (ref) requires the disturbances $\varepsilon$ and $\nu$ to have bounded fourth moments and the variance of $\nu$ to be bounded away from zero. The latter serves as an identification condition ensuring that the leading term in the asymptotic expansion of the triple Lasso estimator has finite variance. Assumption (ref) imposes mild boundedness conditions on the conditional second moments of $\varepsilon$ and $\nu$ and requires the regressors to be uniformly bounded. These conditions simplify the derivations but we note that our results can be extended to allow for unbounded controls. We work with bounded controls to avoid technicalities and to keep the arguments transparent. Finally, Assumption (ref) formalizes the requirement that the subsamples used in cross-fitting contain a non-vanishing fraction of the observations in the full sample. Since the number of folds $K$ is fixed, this condition guarantees that the number of observations in each subsample is of order $n$.

Assumption (ref) imposes approximate sparsity and regularity conditions that are standard in the analysis of Lasso estimators. Parts (ref) and (ref) require that the parameters $\boldsymbol{\gamma}_0$ and $\boldsymbol{\theta}_0$ admit sparse approximations $\bar\boldsymbol{\gamma}_0$ and $\bar\boldsymbol{\theta}_0$ with sparsity indices $s_{\boldsymbol{\gamma}}$ and $s_{\boldsymbol{\theta}}$ and formalize the quality of these approximations by requiring the approximation errors to be small in both $\ell_2$ and $\ell_1$ norms as well as in the prediction norm induced by the controls $\boldsymbol{X}$. Comparable conditions appear, for example, in BCH14 and CHS15. Parts (ref) and (ref) impose the usual rate and sparsity properties of the Lasso estimators $\widehat\boldsymbol{\gamma}_k$ and $\widehat\boldsymbol{\phi}_k$ obtained in the auxiliary regressions. These conditions are satisfied under typically used assumptions on the design matrix, as explained for example in BC11, and ensure that the preliminary estimators are sufficiently accurate for the sparse approximations $\bar\boldsymbol{\gamma}_0$ and $\bar\boldsymbol{\phi}_0$. Part (ref) imposes a classical sparse eigenvalue condition on the population Gram matrix $\mathrm{E}[\boldsymbol{X}\boldsymbol{X}^\top]$, meaning that sparse subsets of the regressors are not excessively collinear. Finally, Part (ref) requires the sets $\widehat T_k$ used in the construction of the triple Lasso estimator to be sufficiently sparse. In light of Part (ref), this condition holds trivially if $\widehat T_{k,1} = \emptyset$.

Assumption (ref) imposes regularity conditions in terms of convergence rates on the node-wise Lasso estimators used to approximate the inverse Gram matrix $\boldsymbol{\Theta}_0 = \mathrm{E}[\boldsymbol{X}\boldsymbol{X}^\top]^{-1}$. These convergence rates are well-established in the literature on estimation of inverse Gram matrices; see, for example, G16. Part (ref) requires that the estimators $\widetilde\boldsymbol{\Theta}_k$ provide sufficiently accurate approximations to the inverse of the empirical Gram matrix computed on the subsample $I(-k)$. Part (ref) ensures that each row of the estimator $\widetilde\boldsymbol{\Theta}_k$ remains bounded in the $\ell_1$ norm. Assumption (ref) restricts the growth of the sparsity indices relative to the sample size. In particular, it requires that the sparsity levels $s_{\boldsymbol{\gamma}}$ and $s_{\boldsymbol{\phi}}$ do not grow too quickly compared to $n$ and $\log p$. Note also that Assumption (ref) is weaker than the corresponding conditions used to derive the $\sqrt n$-consistency and asymptotic normality result for the double Lasso estimator; see BCH14 for the case of the original double Lasso estimator and CCDDHNR18 for the case of the double Lasso estimator combined with cross-fitting.

We are now ready to establish the main result of this paper, which shows that the triple Lasso estimator $\widehat\beta^{TL}$ is $\sqrt n$-consistent and asymptotically normal under our assumptions. The proof of this theorem can be found in the Appendix.

The result in Theorem (ref) for the triple Lasso estimator can be compared with the corresponding result for the double Lasso estimator. In particular, it follows from the proof of Theorem (ref) that the cross-fitted double Lasso estimator $\widehat\beta^{DL}$ satisfies

which is essentially a well-known result; e.g., see CCDDHNR18. Both estimators are therefore $\sqrt n$-consistent and properly centered asymptotically normal if the corresponding remainder terms in asymptotic linear representations (ref) and (ref) are asymptotically vanishing.

In turn, to compare the remainder terms, suppose first that $R_4=o(1)$ so that the double Lasso estimator is $\sqrt n$-consistent and asymptotically normal. We claim that then the remainder terms $R_1$, $R_2$, and $R_3$ for the triple Lasso estimator are {\em always at most of the same order} as the remainder term $R_4$ for the double Lasso estimator and {\em often of smaller order}. Indeed, $$ \frac{R_1}{R_4} = \frac{1}{s_{\boldsymbol{\gamma}}^{1/4}\sqrt{\log p}} = o(1) $$ whenever either $p\to\infty$ or $s_{\boldsymbol{\gamma}}\to\infty$, and $$ \frac{R_2}{R_4} = \frac{s_{\boldsymbol{\gamma}}\sqrt{\log p}}{\sqrt n} \leq R_4 = o(1) $$ since $s_{\boldsymbol{\phi}}\ge s_{\boldsymbol{\gamma}}$ by the definition of $s_{\boldsymbol{\phi}}$. Also, in the {\em worst} case, the remainder term $R_3$ is of the same order as $R_4$: \[ |R_3| \le \sqrt{s_{\boldsymbol{\phi}}\log p} \left( \frac{1}{K}\sum_{k=1}^K \| \bar\boldsymbol{\gamma}_0-\widehat\boldsymbol{\gamma}_k \|_2 + \sqrt{\mathrm{E}[|\boldsymbol{X}^\top(\boldsymbol{\gamma}_0-\bar\boldsymbol{\gamma}_0)|^2]} + \| \boldsymbol{\gamma}_0-\bar\boldsymbol{\gamma}_0 \|_2 \right) \lesssim_P R_4, \] since $\widehat T_k$ contains the support of $\widehat\boldsymbol{\gamma}_k$, but {\em often} converges to zero faster than $R_4$. Indeed, the latter occurs whenever the approximation errors $\mathrm{E}[|\boldsymbol{X}^\top(\boldsymbol{\gamma}_0-\bar\boldsymbol{\gamma}_0)|^2]$ and $\|\boldsymbol{\gamma}_0-\bar\boldsymbol{\gamma}_0\|_2^2$ vanish faster than the upper bound $s_{\boldsymbol{\gamma}}/n$ imposed in Assumption (ref).(ref) and when most coefficients in the approximating vector $\bar\boldsymbol{\gamma}_0$ are not approaching zero too quickly. In the extreme case where $\boldsymbol{\gamma}_0$ is exactly sparse, so that $\boldsymbol{\gamma}_0-\bar\boldsymbol{\gamma}_0=\boldsymbol{0}_p$, and all nonzero coefficients in $\boldsymbol{\gamma}_0$ are sufficiently separated away from zero, so that the Lasso estimators $\widehat\boldsymbol{\gamma}_k$ perform consistent screening, the remainder term $R_3$ is actually {\em equal} to zero. This explains why the remainder terms in the asymptotic linear representation for the triple Lasso estimator often vanish faster than that of the double Lasso estimator, leading to more accurate inference. In Section (ref), we illustrate this phenomenon in Monte Carlo simulations, where we demonstrate that the triple Lasso estimator yields confidence intervals with a better coverage control. In fact, the improvement in coverage control is rather dramatic in some cases.

In addition, we note that the remainder term $R_3$ can be further reduced by enlarging the sets $\widehat T_k$. For example, in practice one can arrange the entries of \[ S_k := \widetilde\boldsymbol{\Theta}_k \mathrm{E}_{I(-k)}[\boldsymbol{X}_i(D_i-\boldsymbol{X}_i^\top\widehat\boldsymbol{\gamma}_k)] \] in decreasing order of absolute value and include the indices of, say, the $L$ largest entries in $\widehat T_{k,1}$, provided that they are not already contained in $\widehat T_{k,0}$. This strategy is motivated by Lemma (ref) and Assumption (ref), which imply that the vectors $S_k$ approximate the corresponding vectors $\boldsymbol{\gamma}_0-\widehat\boldsymbol{\gamma}_k$. Allowing $L$ to grow beyond what is permitted by Assumption (ref).(ref) would reduce $R_3$ at the expense of increasing $R_1$ and $R_2$, which reflects a familiar bias--variance trade-off. In this paper, however, we simply set $L=0$ and leave the question of how to choose $L$ so as to balance these effects for future work.

Finally, in the exactly sparse case with non-zero coefficients of $\boldsymbol{\gamma}_0$ being sufficiently separated away from zero, so that $R_3=0$, it is actually possible that the triple Lasso estimator remains asymptotically normal even when the double Lasso estimator fails to be asymptotically normal. This can occur when $s_{\boldsymbol{\gamma}}$ grows sufficiently slowly relative to $s_{\boldsymbol{\phi}}$, and $$ \frac{s_{\boldsymbol{\gamma}}s_{\boldsymbol{\phi}}(\log p)^2}{n} \to \infty\quad\text{but}\quad \frac{\sqrt{s_{\boldsymbol{\gamma}}}s_{\boldsymbol{\phi}}\log p}{n}\to 0\quad\text{and}\quad \frac{s_{\boldsymbol{\gamma}}^{3/2}\sqrt{s_{\boldsymbol{\phi}}}(\log p)^{3/2}}{n}\to 0 $$ since in this case we have $R_1\to0$ and $R_2\to0$ but $R_4\to\infty$ as $n\to\infty$.

Extensions

In this section, we provide a general recursive formula for constructing a moment function satisfying the Neyman orthogonality condition to any order in a Z-estimation problem. We then apply the general formula to derive the second-order Neyman orthogonal moment function in an M-estimation problem with a single-index structure.

General Formula

Consider a $Z$-estimation problem

where $\boldsymbol{Z}$ is an observable random vector, $\beta_0\in\mathbb{R}$ is the parameter of interest, $\boldsymbol{\theta}_0\in\mathbb{R}^p$ is a vector of nuisance parameters, and $(\bar f,\bar u)$ is a pair of known functions, with $\bar u$ being $\mathbb{R}^p$-valued. In this subsection, we describe a method to construct a $k$th order Neyman orthogonal moment function for estimating $\beta_0$, for any integer $k\geq 1$.

To construct our moment function, we proceed recursively. Let $\boldsymbol{Z}^{(1)},\boldsymbol{Z}^{(2)},\dots$ be independent copies of $\boldsymbol{Z}$ and denote $\boldsymbol{Z}^k = (\boldsymbol{Z}^{(1)},\dots,\boldsymbol{Z}^{(k)})$ if $k\geq 1$ and $\boldsymbol{Z}^k = \boldsymbol{Z}$ if $k=0$. Fix any $k\geq 1$ and suppose that we have already constructed moment functions $F\colon\mathbb{R}\times\mathbb{R}^q\to\mathbb{R}$ and $U\colon\mathbb{R}\times\mathbb{R}^q\to\mathbb{R}^q$ of the form $F(\beta,\boldsymbol{\eta}) = \mathrm{E}[f(\boldsymbol{Z}^{k-1},\beta,\boldsymbol{\eta})]$ and $U(\beta,\boldsymbol{\eta}) = \mathrm{E}[u(\boldsymbol{Z}^{k-1},\beta,\boldsymbol{\eta})]$ such that the pair $(\beta_0,\boldsymbol{\eta}_0)$ solves the system of moment equations $$

$$ where $\boldsymbol{\eta}_0 = (\eta_{0,1},\dots,\eta_{0,q})^\top \in\mathbb{R}^q$ is a nuisance parameter, and $F$ satisfies the $(k-1)$th order Neyman orthogonality condition: $$ \nabla_{\boldsymbol{\eta}}^{m}F(\beta_0,\boldsymbol{\eta}_0) = \boldsymbol{0}_q^{\otimes m},\quadfor all m\in[k-1], $$ where $\nabla_{\boldsymbol{\eta}}^m F(\beta_0,\boldsymbol{\eta}_0)$ is a tensor in $(\mathbb{R}^q)^{\otimes m}:=\mathbb{R}^{q\times \dots \times q}$ defined by $$ (\nabla_{\boldsymbol{\eta}}^mF(\beta_0,\boldsymbol{\eta}_0))_{j_1,\dots, j_m} = \frac{\partial^mF(\beta_0,\boldsymbol{\eta}_0)}{\partial\eta_{j_1}\dots,\partial\eta_{j_m}},\quadfor all j_1,\dots,j_m \in [q], $$ and for any vector $\boldsymbol{a} = (a_1,\dots,a_q)^\top$ in $\mathbb{R}^q$, we use $\boldsymbol{a}^{\otimes m}$ to denote the tensor in $(R^q)^{\otimes m}$ defined by $$ (\boldsymbol{a}^{\otimes m})_{j_1,\dots, j_m} = a_{j_1}\dots a_{j_m},\quad for all j_1,\dots,j_m \in [q]. $$ Now, let

where the latter includes the mode-wise tensor-matrix product notation: for any matrix $\mathbf C = (C_{j_1,j_2})_{j_1,j_2=1}^q$ in $\mathbb{R}^{q\times q}$ and any tensor $\mathbf D = (D_{j_1,\dots,j_m})_{j_1,\dots,j_m=1}^q$ in $(\mathbb{R}^q)^{\otimes m}$, the product $\mathbf C^{\otimes m}\mathbf D$ is the tensor in $(\mathbb{R}^q)^{\otimes m}$ defined by $$ (\mathbf C^{\otimes m}\mathbf D)_{j_1,\dots,j_m} = \sum_{i_1,\dots,i_m=1}^q C_{j_1,i_1}\dots C_{j_m, i_m}D_{i_1,\dots,i_m},\quad\text{for all }j_1,\dots,j_m \in [q]. $$ Then define the moment functions $\tilde F\colon \mathbb{R}\times \mathbb{R}^{q + q^2 + q^k}\to\mathbb{R}$ and $\tilde U\colon \mathbb{R}\times \mathbb{R}^{q + q^2 + q^k}\to\mathbb{R}^{q + q^2 + q^k}$ by setting $$ \tilde F(\beta,\tilde\boldsymbol{\eta}) := F(\beta,\boldsymbol{\eta}) - \frac{1}{k!} \mathbf{B}[(U(\beta,\boldsymbol{\eta}))^{\otimes k}], $$ and $$ \tilde U(\beta,\tilde\boldsymbol{\eta}) := \left(

\right), $$ where $\tilde\boldsymbol{\eta} = (\boldsymbol{\eta}^\top,\mathrm{vec}(\mathbf{A})^\top,\mathrm{vec}(\mathbf{B})^\top)^\top$ and we used the contraction of tensors notation: for any tensors $\mathbf D = (D_{j_1,\dots,j_m})_{j_1,\dots,j_m=1}^q$ and $\mathbf E = (E_{j_1,\dots,j_m})_{j_1,\dots,j_m=1}^q$ in $(\mathbb{R}^q)^{\otimes m}$, the contraction $\mathbf D[\mathbf E]$ is a scalar in $\mathbb{R}$ defined by $$ \mathbf D[\mathbf E] = \sum_{j_1,\dots,j_m=1}^q D_{j_1,\dots,j_m}E_{j_1,\dots,j_m}. $$ Note in passing that $\tilde F$ and $\tilde U$ take the form $\tilde F(\beta,\tilde\boldsymbol{\eta}) = \mathrm{E}[\tilde f(\boldsymbol{Z}^{k},\beta,\tilde\boldsymbol{\eta})]$ and $\tilde U(\beta,\tilde \boldsymbol{\eta}) = \mathrm{E}[u(\boldsymbol{Z}^{k},\beta,\tilde\boldsymbol{\eta})]$, maintaining the recursion. The following lemma shows that $\tilde F$ is the desired moment function. The proof of this lemma can be found in the Appendix.

Application to M-Estimation with Single-Index Structure

We now apply the general formula obtained in the previous subsection to derive the second-order Neyman orthogonal moment function for estimating $\beta_0$ in the M-estimation problem with a single-index structure $$ (\beta_0,\boldsymbol{\theta}_0):=\operatornamewithlimits{argmin}\limits_{(\beta,\boldsymbol{\theta})\in\mathbb{R}\times\mathbb{R}^p}\mathrm{E}[m(D\beta+\boldsymbol{X}^\top\boldsymbol{\theta},\boldsymbol{Y})], $$ where $m\colon\mathbb{R}\times\mathcal Y\to\mathbb{R}$ is a known smooth loss function, $\boldsymbol{Y}\in\mathcal Y$ is one or more outcome variables, $D\in\mathbb{R}$ is a regressor of interest, $\boldsymbol{X}\in\mathbb{R}^p$ is a vector of controls, $\beta_0\in\mathbb{R}$ is the parameter of interest, and $\boldsymbol{\theta}_0\in\mathbb{R}^p$ is a vector of nuisance parameters.

To this end, denote the first, second, and third derivatives of $m$ with respect to its first argument by $m_1'$, $m_{11}''$, and $m_{111}'''$, respectively, and let $\boldsymbol{\mu}_0\in\mathbb{R}^p$ be a solution to the moment equation $$ \mathrm{E}[m_{11}''(D\beta_0 + \boldsymbol{X}^\top\theta_0,\boldsymbol{Y})(D-\boldsymbol{X}^\top\boldsymbol{\mu})\boldsymbol{X}] = \boldsymbol{0}_p. $$ Then the first-order Neyman orthogonal moment function for estimating $\beta_0$ is $F\colon\mathbb{R}\times\mathbb{R}^{2p}\to\mathbb{R}$ given by $$ F(\beta,\boldsymbol{\eta}):=\mathrm{E}[m_1'(D\beta + \boldsymbol{X}^\top\boldsymbol{\theta},\boldsymbol{Y})(D-\boldsymbol{X}^\top\boldsymbol{\mu})],\quad \boldsymbol{\eta}:=(\boldsymbol{\theta}^\top,\boldsymbol{\mu}^\top)^\top; $$ e.g. see GBRD14 or chetverikov2025selecting. Thus, setting $$ U(\beta,\boldsymbol{\eta}):=\left(

\right), $$ we fit our problem into a setting of the previous subsection with $\boldsymbol{\eta}_0:= (\boldsymbol{\theta}_0^\top,\boldsymbol{\mu}_0^\top)^\top$.

Now, to apply the general formula from the previous subsection, let

Then $$ \frac{\partial U(\beta_0,\boldsymbol{\eta}_0)}{\partial\boldsymbol{\eta}^\top} =

, \qquad \nabla_{\boldsymbol{\eta}}^2 F(\beta_0,\boldsymbol{\eta}_0) =

. $$ Hence, with $\mathbf{A}_0$ and $\mathbf{B}_0$ defined by \eqref{eq: a0 and b0 definition}, which simplifies to \eqref{eq: a0 and b0 simplification} for $k=2$, we have \[ \mathbf{B}_0 =

. \] Therefore, the second-order Neyman orthogonal moment function is

where $\tilde\boldsymbol{\eta} := (\boldsymbol{\theta}^\top,\boldsymbol{\mu}^\top,\mathrm{vec}(\mathbf G)^\top,\mathrm{vec}(\mathbf H)^\top)^\top$ and $\tilde\boldsymbol{\eta}_0 := (\boldsymbol{\theta}_0^\top,\boldsymbol{\mu}_0^\top,\mathrm{vec}(\mathbf G_0)^\top,\mathrm{vec}(\mathbf H_0)^\top)^\top$. It is also straightforward to verify that $\beta_0$ solves the desired moment equation $\tilde F(\beta,\tilde\boldsymbol{\eta}_0) = 0$.

Monte Carlo Simulations

In this section, we study the finite-sample behavior of the triple Lasso estimator relative to the standard cross-fitted double Lasso estimator in the linear regression model from Section (ref). In our Monte Carlo exercises, we vary the sparsity of the nuisance vector $\boldsymbol{\gamma}_0$, the correlation structure of the controls $\boldsymbol{X}$, and the sample/problem size, while keeping the target parameter fixed at $\beta_0 = 1$.

Design

For each Monte Carlo replication, we generate i.i.d.\ observations $(\boldsymbol{X}_i,D_i,Y_i)$, $i=1,\dots,n$, from model (ref)--(ref) with $\beta_0 = 1$, $\nu\sim\mathcal{N}(0,1)$ and $\varepsilon\sim\mathcal{N}(0,1)$ independent of each other and of $\boldsymbol{X}$. The controls follow a Gaussian design, \[ \boldsymbol{X} \sim \mathcal{N}(\boldsymbol{0},\boldsymbol{\Sigma}_0), \qquad \boldsymbol{\Sigma}_0:=\boldsymbol{\Sigma}_0(\rho):=\sbr[1]{\rho^{|j-k|}}_{j,k=1}^p, \] so that the scalar $\rho \in \{0,0.2,0.4,0.6,0.8\}$ governs the correlation between controls. Note that with $\boldsymbol{X}\sim \mathcal{N}(\boldsymbol{0},\boldsymbol{\Sigma}_0)$ and $\boldsymbol{\Sigma}_0$ of the Toeplitz form given above, the inverse $\boldsymbol{\Omega}_0=\boldsymbol{\Sigma}_0^{-1}$ takes a tri-diagonal form. Joint normality therefore leads to any single control $X_j$ being conditionally normally distributed given the other controls $\boldsymbol{X}_{-j}$. Unpacking $\boldsymbol{\Omega}_0$, we get \[ \left.X_j \middle| \boldsymbol{X}_{-j}=\boldsymbol{x}_{-j}\right.\sim

\] In particular, the conditional means are sparse linear functions of the conditioning variables.

We consider the sample sizes $n \in \{500,1000,2000\}$ and limit attention to the high-dimensional regime $p=n/2$. The structure of the outcome nuisance vector $\boldsymbol{\theta}_0$ is held fixed across designs and is approximately sparse with geometrically decaying entries, \[ \theta_{0,j} = (0.5)^{\,j-1}, \quad j=1,\ldots,p. \] The treatment nuisance vector $\boldsymbol{\gamma}_0$ also has (at least) geometrically decaying entries \[ \gamma_{0,j} =

\] but its sparsity varies across three regimes:

For each of the $3\times5\times3=45$ design points $(n,\rho,s_{\boldsymbol{\gamma}})$, we run $2{,}000$ Monte Carlo replications.

The approximate sparsity of $\boldsymbol{\theta}_0$ is already challenging for the double lasso. Moving from exact to approximate sparsity of $\boldsymbol{\gamma}_0$ gradually weakens the environment in which even our second-order orthogonal procedure is expected to perform well since approximate sparsity design implies larger effective sparsity index $s_{\boldsymbol{\gamma}}$.

Estimation and Implementation

We compare two estimators. First, the double Lasso estimator is implemented as the cross-fitted residual-on-residual estimator using the DML1 aggregation rule of CCDDHNR18. In each fold, we estimate the treatment regression $D$ on $\boldsymbol{X}$ and the reduced-form regression $Y$ on $\boldsymbol{X}$ by lasso on the training sample, evaluate the fitted values on the holdout sample, and form the fold-specific orthogonal score based on the residuals.

Second, the triple Lasso estimator is implemented as the cross-fitted adjusted-score estimator. Relative to the double Lasso estimator, it augments the holdout score by the estimated second-order correction term built from estimates of selected rows of the inverse control correlation matrix. Operationally, after estimating $\widehat{\boldsymbol{\gamma}}_k$ and $\widehat{\boldsymbol{\phi}}_k$ on the training fold, we define the selected support \( \widehat T_k = \{j : \widehat\gamma_{k,j} \neq 0\} \), estimate only the rows of $\boldsymbol{\Theta}_0 = \boldsymbol{\Sigma}_0^{-1}$ indexed by $\widehat T_k$ using node-wise Lasso, and then evaluate the adjusted score on the holdout observations. We then use DML1 aggregation to produce the triple lasso estimate.

Both procedures use the same random $K=5$ folds in each replication. All Lasso penalties are determined by plug-in rules of the Bickel--Ritov--Tsybakov type bickel_simultaneous_2009. Specifically, let $n_K := (K-1)n/K$ denote the training-sample size and introduce the Belloni--Chen--Chernozhukov--Hansen baseline penalty level

where $c:=1.1$ and \( \alpha(n_{\text{obs}},p_{\text{pen}}):= 0.1/\log(\mathrm{max}\{p_{\text{pen}},n_{\text{obs}}\}), \) as suggested by belloni_sparse_2012, with $n_{\text{obs}}$ and $p_{\text{pen}}$ being placeholders for the numbers of observations and penalized parameters, respectively. The fold-invariant penalty levels are then set as

where \( \sigma^2_e := \beta_0^2 \sigma_\nu^2 + \sigma_\varepsilon^2 = 2 \) is the variance of the reduced-form error \( e_i = \beta_0 \nu_i + \varepsilon_i, \) and \[ \sigma^2_{X_j\mid \boldsymbol{X}_{-j}}(\rho):=\frac{1-\rho^2}{1+\boldsymbol{1}\{1<j<p\}\rho^2} \] is the conditional variance of control $X_j$ given all other controls $\boldsymbol{X}_{-j}$. Here, $\lambda_\gamma$ is used in the treatment Lasso, $\lambda_\phi$ in the outcome reduced-form Lasso, and the $\{\lambda_{\psi_j}(\rho)\}_{j=1}^p$ are used in the node-wise Lasso step of the triple Lasso estimator.\footnote{ The (infeasible) penalties in (ref) were chosen for computational simplicity. A more realistic comparison would use tuning parameter choices which are feasible in practice. One option is to replace the unknown standard deviations with (initally conservative) proxies to produce feasible analogs of the bickel_simultaneous_2009 penalties. One could refine such proxies in a possibly iterative manner borrowing ideas from belloni_sparse_2012. Alternatively, one could use the bootstrap after cross-validation method proposed in chetverikov2025selecting, which takes into account the correlation structure.}

All Lasso estimators are fit using glmnet in R with standardized controls.\footnote{The glmnet definition of Lasso is based on one half square loss, which cancels out a “2” in the original belloni_sparse_2012 baseline penalty, thus leading to our (ref).} The treatment and reduced-form Lasso regressions include intercepts; node-wise regressions are run without intercepts. The reported standard errors are heteroskedasticity-robust score-based standard errors computed from the corresponding cross-fitted influence function.

Performance Metrics

For each estimator $\widehat\beta$ we report the Monte Carlo

Error bars indicate $\pm 1.96$ Monte Carlo standard errors. To facilitate comparison, for the studentized statistics, we plot the kernel densities instead of histograms.\footnote{All kernel densities are created using the ggplot2 with geom_density. In expectation of an approximately normal distribution, we use a Gaussian kernel and the silverman_density_1986 rule-of-thumb bandwidth (both geom_density defaults).}

Results

Figures (ref)--(ref) show the performance of the triple and double Lasso estimators based on the stated performance metrics in turn.

A central pattern from the figures is that triple Lasso estimator substantially reduces bias precisely in the designs where one would expect second-order effects to matter most. This is most visible in Figure (ref), which plots the squared bias. For example, under approximate sparsity with $(n,p,\rho)=(1000,500,0)$, the squared bias falls from $0.00216$ for the double Lasso estimator to $0.00025$ for the triple Lasso estimator; with $(n, p,\rho)=(2000,1000,0)$, it falls from $0.00082$ to $0.00009$. Similar reductions appear throughout the intermediate-sparsity designs. Note that these findings do not contradict the asymptotic theory in Section (ref), which suggests that the gains of the triple Lasso estimator are largest in the exactly sparse designs with consistent screening. Indeed, our approximately sparse design corresponds to larger effective sparsity index relative to our exact sparsity design, which hurts the double Lasso estimator more than it hurts the triple Lasso estimator by construction.

The variance of the triple Lasso estimator is typically somewhat larger (Figure (ref)), which is the expected cost of estimating the additional score adjustment, but the increase is typically modest relative to the bias reduction in the difficult designs.

Figure (ref) shows that these bias gains translate into meaningful improvements in mean squared error whenever the double Lasso estimator is materially biased. The gains are especially pronounced at low and moderate $\rho$ values and under intermediate or approximate sparsity. For instance, with $(n,p,\rho,s_{\boldsymbol{\gamma}})=(500,250,0,\text{interm.})$, the MSE declines from $0.00630$ to $0.00275$, and with $(1000,500,0,\text{approx.})$ it declines from $0.00301$ to $0.00129$. As $n$ grows, the gap narrows, and in the easiest designs, where $\boldsymbol{\gamma}_0$ is exactly sparse and regressor correlation is high, the two methods become similar and the double Lasso estimator can occasionally have slightly smaller MSE because of its lower variance.

Figure (ref) shows the empirical coverage. The double Lasso estimator exhibits substantial undercoverage in the harder designs, again reflecting residual bias. For example, at $(500,250,0,\text{approx.})$, its empirical $95\%$ coverage is $0.630$; at $(1000,500,0,\text{approx.})$ it is $0.672$; and at $(2000,1000,0,\text{approx.})$ it is still only $0.750$. The triple Lasso estimator moves coverage much closer to the nominal level in those same cells, to $0.909$, $0.923$, and $0.930$, respectively. Even when the improvement in MSE becomes small, the coverage results indicate that the second-order adjustment remains useful for inference.

Figure (ref) shows the corresponding tradeoff in interval length. The triple Lasso confidence intervals are systematically longer than double lasso intervals, which is consistent with the modest increase in variance. The longer intervals should not be viewed as a drawback per se. Rather, they reflect the additional uncertainty captured by the second-order adjustment and are therefore the means by which coverage is brought closer to its nominal level in designs where the first-order double Lasso procedure understates uncertainty.

Finally, Figure (ref) examines the studentized statistics in the low-correlation designs $(\rho=0)$. The density of the double Lasso $t$-statistic is visibly shifted away from the standard normal benchmark, whereas the triple Lasso statistic appears more nearly centered and generally more aligned with the $\mathcal{N}(0,1)$ reference curve. The underlying Monte Carlo moments tell a similar story: for example, under approximate sparsity with $(n,p,\rho)=(1000,500,0)$, the mean of the double Lasso studentized statistic is $1.53$, while the mean for triple Lasso is $0.50$. This pattern is consistent with our second-order orthogonality result.

In Appendix (ref), we report the corresponding densities for $\rho\in\{0.2,0.4,0.6,0.8\}$. Across many of these designs, the triple Lasso studentized statistic also appears better centered and often closer to the $\mathcal{N}(0,1)$ benchmark than the double Lasso statistic, although the magnitude of the difference varies across cells.

Overall, the simulations support the main theoretical message of the paper. When nuisance estimation errors are sufficiently small, triple and double Lasso estimators behave similarly. When second-order terms are non-negligible, however, the triple Lasso estimator delivers a clear reduction in bias and a substantial improvement in inference.