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.
Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.
Finite-Sample Guarantees for High-Dimensional DML
abstractAbstract \quad Debiased machine learning (DML) offers an attractive way to estimate treatment effects in observational settings, where identification of causal parameters requires a conditional independence or unconfoundedness assumption, since it allows to control flexibly for a potentially very large number of covariates. This paper gives novel finite-sample guarantees for joint inference on high-dimensional DML, bounding how far the finite-sample distribution of the estimator is from its asymptotic Gaussian approximation. These guarantees are useful to applied researchers, as they are informative about how far off the coverage of joint confidence bands can be from the nominal level. There are many settings where high-dimensional causal parameters may be of interest, such as the ATE of many treatment profiles, or the ATE of a treatment on many outcomes. We also cover infinite-dimensional parameters, such as impacts on the entire marginal distribution of potential outcomes. The finite-sample guarantees in this paper complement the existing results on consistency and asymptotic normality of DML estimators, which are either asymptotic or treat only the one-dimensional case.
Introduction
A recent strand of literature in econometrics has considered estimation of treatment effects and causal or structural parameters using machine learning (ML) methods chernozhukov2018automatic, belloni2017program, athey2019machine, farrell2021deep. In many observational settings, the treatment or policy whose impact we wish to quantify was not randomly assigned, so a simple comparison of treatment and control groups is confounded by factors that are correlated both with the outcome and with the treatment. We may still be able to identify the causal effect of the treatment, however, if we are willing to make an unconfoundedness assumption, i.e., that the treatment is exogenous conditional on an appropriate set of controls. There are several ways in which ML can be useful in such quasi-experimental research designs. On the one hand, it allows to control flexibly for a large number of covariates, as are typically available in modern datasets. On the other hand, it opens a wide range of possibilities in terms of which types of data can be used as controls, e.g., textual or image data.\footnote{This is increasingly relevant in applied economics. For example, see dube2020monopsony, who use data from job ads in Amazon MTurk to study whether online employers have monopsonistic labor market power, controlling for job characteristics in the form of features learned by a random forest trained on the ad descriptions. Similarly, an IGC project by olken2017evidence, uses Google Street View images to infer J-PAL survey responses in Indonesia; one could imagine using those same images as a proxy for socio-demographic controls.}
The first section of this paper considers inference on a set of parameters that can be expressed as averages of a functional,
\[\theta_{0j} = \mathrm{E}_{P}\left[ m_j(W, \gamma_{0j}) \right], \quad j = 1, \, \ldots, \, p, \]
where $\gamma_{0j}$ is an infinite-dimensional nuisance parameter (for example, a conditional expectation or regression function), and $p$ is potentially large. This setting encompasses, for example, joint inference on the effects of many different treatments or on the impact of a treatment on many different outcomes.
Modern ML methods perform very well in predictive settings, typically by trading off variance and bias through some form of explicit or implicit regularization (and so they offer an attractive methodology to estimate $\gamma_{0j}$ when it is some form of “best predictor”). At the same time, that trade-off means that their convergence rates are slower than the parametric $\sqrt{n}$-rate, causing a first-order bias in the estimation of the target parameters $\theta_{0j}$ that does not, in general, vanish asymptotically. chernozhukov2018double and subsequent work show how to construct estimators of $\theta_{0j}$ that are correctly centered asymptotically by “debiasing” the moment conditions above, making them first-order robust\footnote{More precisely, Neyman-orthogonal, as defined below.} to the ML estimation error. These are known as double or debiased machine learning (DML) estimators.
An applied research may be interested in conducting joint inference about a high-dimensional parameter $\{\theta_{0j}\}_{j=1, \, \ldots, \, p}$. That can be done by constructing simultaneous confidence bands that cover the true parameter with a pre-specified probability (approximately) $1-\alpha$. Based on DML estimators $\{\hat{\theta}_{j}\}_{j = 1, \, \ldots, p}$ of $\{\theta_{0j}\}_{j = 1, \, \ldots, p}$, and an estimator $\hat{\sigma}_{j}^2$ of $\sigma_j^2 = \mathrm{Var}(\sqrt{n}(\hat{\theta}_{j} - \theta_{0j}))$ for each $j = 1,\ldots,p$, we will consider a joint confidence band of the form $\times_{j=1}^p [\hat{\theta}_j \mp n^{-1/2}\hat{\sigma}_jc_{\alpha}]$, where the critical value $c_{\alpha}$ is chosen so that the sup-$t$-statistic satisfies: \[\mathrm{P}_P\left( \max_{1 \leq j \leq p} \sqrt{n}\frac{|\hat{\theta}_{j} - \theta_{0j}|}{\hat{\sigma}_j} \leq c_\alpha \right) \approx 1 - \alpha.\] One way to choose $c_\alpha$ is to rely on the joint asymptotic normality of the $t$-statistics, which is a well-established fact in the literature (chernozhukov2018double, chernozhukov2018automatic, chernozhukov2021automatic), even in the high-dimensional case and for continua of parameters belloni2017program, belloni2018uniformly.
The goal of this paper is to provide finite-sample guarantees for the normal approximation above, in the form of a bound on the Kolmogorov distance between the finite-sample distribution of the sup-$t$-statistic and $\max_{1 \leq j \leq p} Z_j$ for a suitable jointly Gaussian distribution $(Z_1, \ldots, Z_p) \sim \mathcal{N}(0, \Sigma)$, i.e.: \[\sup_{t \in \mathbb{R}} \left| \mathrm{P}_P\left( \max_{1 \leq j \leq p} \sqrt{n}\frac{(\hat{\theta}_j - \theta_{0j})}{\sigma_j} \leq t \right) - \mathrm{P} \left( \max_{1 \leq j \leq p} Z_j \leq t \right)\right|.\] The dependence of the bound on $n$ and $p$ will be made explicit. These guarantees are useful to applied researchers, as they are informative about how far off the coverage of joint confidence bands can be from the nominal level for a given sample size and complexity of the problem. The closest existing result to ours is chernozhukov2021simple, who provide finite-sample guarantees in the single-parameter case, $p = 1$, when using pure sample splitting.\footnote{A precise definition of what we mean by sample splitting will be provided in (ref).} We extend their results in two ways, allowing for a potentially large $p$ and for no sample splitting. To that end, we leverage results from the literature on high-dimensional estimation, including a maximal inequality for empirical processes chernozhukov2014gaussian and new normal approximation results for high-dimensional vectors chernozhukov2021nearly.
In the second section of this paper, we consider a more general moment problem,
\[\mathrm{E}_{P}\left[ \psi_u(W, \theta_{0u}, \eta_{0u}) \right] = 0, \quad u \in \mathcal{U},\] where $\mathcal{U}$ is potentially an uncountable set. In the continuum of target parameters case, the $\{\theta_{0u}\}_{u \in \mathcal{U}}$ could represent, for example, the marginal distribution of an outcome under a treatment, which allows to derive many other interesting statistics (e.g., quantile treatment effects, Gini coefficients, Oaxaca-Blinder decompositions of distributional shifts, etc.). Again, we wish to construct a simultaneous confidence band for $\{\theta_{0u}\}_{u \in \mathcal{U}}$ using DML, based on a normal approximation to the sup-$t$-statistic. We provide finite-sample guarantees for that normal approximation also in this setting. To the best of our knowledge, ours is the first paper to provide non-asymptotic guarantees for the DML estimator of a continuum of parameters, which extend and complement the asymptotic results of belloni2017program, belloni2018uniformly.
\paragraph{Notation} Throughout the paper, we use the following notation. For a random variable $W \in \mathcal{W}$, distributed according given probability measure $P$ on $\mathcal{W}$, we denote by $\mathrm{E}_{P}\left[ \cdot \right]$ the expectation with respect to $P$, i.e., $\mathrm{E}_{P}\left[ f(W) \right] = \int f(w) \mathrm{d}P(w)$ for a suitably measurable and integrable $f$. We denote by $\mathbb{E}_{n}\left[ \cdot \right]$ the average of a sample $\{W_i\}_{i=1}^n$ of size $n$, i.e., $\mathbb{E}_{n}\left[ f(W) \right] = n^{-1}\sum_{i=1}^n f(W_i)$. We use $\mathbb{G}_n [\cdot] = \mathbb{G}_{n,p} [\cdot]$ for an empirical process $\sqrt{n}(\mathbb{E}_n[\cdot] - \mathrm{E}_P[\cdot])$ over a class $\mathcal{F}$ of suitable measurable and integrable functions $f : \mathcal{W} \rightarrow \mathbb{R}$, i.e., \[\mathbb{G}_n[f] = \mathbb{G}_n[f(W)] = \frac{1}{\sqrt{n}}\sum_{i=1}^n \left(f(W_i) - \mathrm{E}_{P}\left[ f(W) \right] \right).\]
We denote by $L^q = L^q(P)$ the space of functions with finite $q$-th absolute moments with respect to $P$, $\{f: \mathcal{W} \rightarrow \mathbb{R} \, : \, \mathrm{E}_{P}\left[ |f(W)|^q \right] < \infty\}$. For $f \in L^q(P)$, we denote by $\Vert f \Vert_{P,q} = (\mathrm{E}_{P}\left[ |f(W)|^q \right])^{1/q}$, the $L^q$ norm. For a bounded function $f : \mathcal{W} \rightarrow \mathbb{R}$, we denote by $\Vert \cdot \Vert_{\infty}$ the sup norm, i.e., $\Vert f \Vert_{\infty} = \sup_{w \in \mathcal{W}} |f(w)|$.
For a function $f : \mathbb{R} \rightarrow \mathbb{R}$, we denote by $\partial_r f(r) = f'(r)$ the derivative with respect to $r$. For a functional $F : \mathcal{F} \rightarrow \mathbb{R}$ over a class of functions $\mathcal{F}$, we define the Gateaux derivative of $F$ at $f$ in the direction $u \in \mathcal{F}$ as $\partial_r F(f + r u) |_{r = 0}$. We say that $F$ is Gateaux differential at $F$ if that derivative exists for all $u \in \mathcal{F}$.
For a function class $\mathcal{F}$ endowed with a norm $\Vert \cdot \Vert_{\mathcal{F}}$, and for any $\varepsilon > 0$, we define the covering number $N(\varepsilon, \mathcal{F}, \Vert \cdot \Vert_{\mathcal{F}})$ as the smallest number of closed balls with radius $\varepsilon$ that could cover $\mathcal{F}$. Denote by $F$ a measurable envelope for $\mathcal{F}$, i.e., a function such that $F \geq \sup_{f\in \mathcal{F}}|f|$. The uniform entropy number is, for any $\varepsilon > 0$, $\log \sup_Q N(\varepsilon, \mathcal{F}, \Vert \cdot \Vert_{Q,2})$, where the supremum is taken over any finitely discrete probability measure $Q$ such that $\Vert F \Vert_{Q,2} > 0$.
Averages of Many Linear Functionals
In this section we extend the results of chernozhukov2021simple to the high-dimensional case. Suppose we have access to an i.i.d. sample $\{(Y_i, W_i)\}_{i=1}^n$ from a probability distribution $P$, where $Y \in \mathcal{Y}$ denotes outcomes of interest and $W \in \mathcal{W}$ are other observed data. Our goal is to construct a simultaneous confidence band for the set of (scalar) parameters $\{\theta_{0j}\}_{j=1, \, \ldots, \, p}$, satisfying the following moment condition:
equation[equation omitted — 131 chars of source]
where the moment functional $m_j(\cdot, \cdot)$ may depend on the unknown true value of an infinite-dimensional nuisance parameter $\gamma_{0j} \in \Gamma$, e.g., a conditional expectation. The set $\Gamma \subset L^2$ is assumed to be a linear function space, and can be used to encode restrictions on $\gamma_{0j}$, such as smoothness chernozhukov2021simple. We assume that an ML estimator $\hat{\gamma}_j$ of $\gamma_{0j}$ can be obtained from the same data.
Below we give some examples where this setting may apply.
example[ATE of many treatments]
Researchers may be interested in estimating the average treatment effect (ATE) of many different treatments, combinations of treatments or dosages. In observational settings, where treatments are not randomly assigned, ATEs may still be identified under an unconfoundedness assumption rosenbaum1983central.
Let $D \in \mathcal{D}$ denote the treatment variable, where $\mathcal{D} = \{d_0, \ldots, d_p\}$ is the set of different treatment profiles, with $d_0$ denoting no treatment. Let $Y(d)$ for $d \in \mathcal{D}$ denote potential outcomes, so that the observed outcome is $Y = Y(D)$. Suppose that the researcher has access to a set of control variables $X$ such that $D$ is independent of $Y(d)$ given $X$ for all $d \in \mathcal{D}$. In that case, $\mathrm{E}_{P}\left[ Y(d) \right] = \mathrm{E}_{P}\left[ \gamma_0(d_j, X) \right]$, where $\gamma_0(d, x) = \mathrm{E}_{P}\left[ Y \mid D = d, X = x \right]$. (Notice that, with this formulation, $\gamma_0$ does not depend on $j$.) The ATE of treatment profile $d_j$ with respect to no treatment is: \[\theta_{0j} = \mathrm{E}_{P}\left[ m_j(W, \gamma_0) \right], \quad m_j(w, \gamma_0) = \gamma_0(d_j, x) - \gamma_0(d_0, x), \quad j = 1, \, \ldots, p.\]
In modern datasets, the number of available controls $X$ may be large, so that ML methods may be especially well suited to estimate the nuisance parameter $\gamma_{0j}$.
example[ATE on many outcomes]
Our framework also provides guarantees for uniform inference on the ATE of a binary treatment $D \in \{0,1\}$ on a large set of outcomes $Y = (Y_1, \ldots, Y_p)$ using the sup-$t$-statistic. In this case, under the appropriate unconfoundedness assumption, we have: \[\theta_{0j} = \mathrm{E}_{P}\left[ m_j(W, \gamma_{0j}) \right], \quad m_j(w, \gamma_{0j}) = \gamma_{0j}(1, x) - \gamma_{0j}(0, x), \quad j = 1, \, \ldots, p,\] where now $\gamma_{0j}(d,x) = \mathrm{E}_{P}\left[ Y_j \mid D = d, X = x \right]$.
example[Policy optimization]
Consider again a binary treatment $D \in \{0,1\}$ and an outcome of interest $Y$. Policymakers may wish to select the best treatment assignment policy based on a set of characteristics $X$, where a treatment assignment policy is a mapping $\pi : X \rightarrow \{0,1\}$. Suppose we only have access to observational data, collected under an unknown treatment policy. We want to evaluate and compare the effects of an alternative set of candidate policies $\{\pi_1, \ldots, \pi_p\}$, where the average effect of policy $\pi_j$ is: \[\theta_{0j} = \mathrm{E}_{P}\left[ m_j(W,\gamma_0) \right], \quad m_j(w, \gamma_{0j}) = \gamma_{0}(0, x) + \pi_j(X)(\gamma_{0}(1, x) - \gamma_{0}(0, x) ), \quad j = 1, \, \ldots, p.\] Here, $\gamma_0(d,x) = \mathrm{E}_{P}\left[ Y \mid D = d, X = x \right]$ again does not depend on $j$. Other recent literature has also considered doubly-robust approaches to policy optimization based on observational data (e.g., athey2021policy).
Given an estimator $\hat{\gamma}_j$ of $\gamma_{0j}$, it may seem natural to estimate $\theta_{0j}$ by the empirical analog of (ref), \[\Check{\theta}_j = \mathbb{E}_{n}\left[ m_j(W, \hat{\gamma}_j) \right], \quad j = 1, \, \ldots, \, p.\]
However, when $\hat{\gamma}_j$ is obtained using modern ML methods, $\hat{\gamma}_j$ typically converges to $\gamma_{0j}$ more slowly than the parametric $\sqrt{n}$-rate. In this formulation, chernozhukov2018double show that the bias from using $\hat{\gamma}_j$ instead of $\gamma_{0j}$ is of first-order magnitude, so that $\Check{\theta}_j$ is not asymptotically centered around the true value $\theta_{0j}$. Inference based on $\Check{\theta}_j$ will thus be incorrect if one fails to account for that.
Following chernozhukov2018double and subsequent work, we proceed by adjusting the moment condition (ref) to make it immune, to a first order, against the estimation error in $\hat{\gamma}_j$. This method is known as double or debiased machine learning (DML). Below we collect some existing results (Propositions (ref) and (ref)) that show how to “debias” the moment condition (ref) in the particular case of linear, mean-square continuous functionals of $\gamma$.
assumption[Linearity and mean-square continuity] For all $j = 1, \, \ldots, \, p$, the moment functional $\gamma \mapsto m_j(w, \gamma)$ is linear and mean-square continuous, i.e., there exists $\Bar{Q} < \infty$ such that \[\mathrm{E}_{P}\left[ m_j(W,\gamma)^2 \right] \leq \Bar{Q}^2 \mathrm{E}_{P}\left[ \gamma(W)^2 \right] \quad \text{for all } \gamma \in \Gamma, j = 1, \, \ldots, \, p.\]
proposition[Riesz representation] Suppose (ref) holds. Then, there exists a unique $\alpha_{0j} \in \mathrm{cl}(\mathrm{span} (\Gamma))$ such that \[\mathrm{E}_{P}\left[ m_j(W,\gamma) \right] = \mathrm{E}_{P}\left[ \alpha_{0j}(W)\gamma(W) \right] \quad \text{for all } \gamma \in \Gamma, j = 1, \, \ldots, \, p.\]
proofThis is a consequence of the Riesz representation theorem. For a proof under more general conditions in the context of classic semiparametric theory, see, e.g., newey1994asymptotic, ichimura2022influence. See also, e.g., chernozhukov2018automatic for a more detailed discussion of the role of the Riesz representer in DML.
remarkIt is easy to see that the functionals in examples (ref) to (ref) are linear. One can also show mean-square continuity under appropriate regularity conditions (e.g., an overlap condition, $0 < \underline{p} \leq \mathrm{P}_P\left( D = d \mid X \right) \leq \overline{p} < 1$ a.s. for all possible treatments or treatment profiles $d$), see e.g., chernozhukov2018automatic. As a consequence, the Riesz representer exists in examples (ref) to (ref).
In general, the true Riesz representer $\alpha_{0j}$ will be unknown. We assume that an estimator $\hat{\alpha}_j$ can be obtained from the same data. In some cases, the explicit form of $\alpha_{0j}$ can be derived. For instance, it can be shown that \[\alpha^*_{0j}(D,X) = \frac{\mathbf{1}\{D = d_j\}}{\mathrm{P}_P\left( D = d_j \mid X \right)} - \frac{\mathbf{1}\{D = d_0\}}{\mathrm{P}_P\left( D = d_0 \mid X \right)}\] is a Riesz representer for the functional (ref).\footnote{With a restricted semiparametric model $\Gamma$, it is possible that $\alpha^*_{0j} \notin \mathrm{cl}(\mathrm{span} (\Gamma))$, in which case the unique or minimal Riesz representer of (ref) would be $\alpha_{0j} = \mathrm{Proj}(\alpha^*_{0j} \mid \mathrm{cl}(\mathrm{span} (\Gamma)))$.} An estimator $\hat{\alpha}_j$ can be then constructed by plugging in a non-parametric estimate of the propensity scores $\mathrm{P}_P\left( D = d \mid X \right)$. A more recent strand of literature considers automatic estimation of $\alpha_{0j}$, where knowledge of the explicit form of $\alpha_{0j}$ is not required chernozhukov2018automatic, chernozhukov2020adversarial, chernozhukov2021automatic.
We consider a point estimator for $\{\theta_{0j}\}_{j=1, \, \ldots, \, p}$ based on the following augmented moment condition:
equation[equation omitted — 170 chars of source]
For the ATE examples, this is the augmented inverse propensity-weighted estimator (AIPW) of robins1994estimation. The addition of the term $\alpha_{0j}(W)(Y - \gamma_{0j}(W))$ can be seen as “debiasing” the moment condition (ref), since it makes it robust to estimation errors in $\hat{\gamma}_{0j}$ and $\hat{\alpha}_{0j}$ in a sense made explicit by the proposition below.
proposition[Neyman orthogonality and double robustness] Let $Z = (Y, W)$, and $\psi_j(Z, \theta, \gamma, \alpha)$ denote the augmented score: \[\psi_j(Z, \theta, \gamma, \alpha) = m_j(W, \gamma) + \alpha(W)(Y - \gamma(W)) - \theta.\] We have: \begin{enumerate}[(i)]
• (Neyman orthogonality) The Gateaux derivative maps of $\mathrm{E}_{P}\left[ \psi_j(Z, \theta, \gamma, \alpha) \right]$ with respect to $\gamma$ and $\alpha$ are 0 at $(\theta_{0j}, \gamma_{0j}, \alpha_{0j})$:
\begin{align*}
\partial_r \mathrm{E}_{P}\left[ \psi_j(Z, \theta_{0j}, \gamma_{0j} + r(\gamma - \gamma_{0j}), \alpha_{0j}) \right] |_{r = 0} &= 0 \quad for all \gamma \in \Gamma,\\ \partial_r \mathrm{E}_{P}\left[ \psi_j(Z, \theta_{0j}, \gamma_{0j}, \alpha_{0j} + r(\alpha - \alpha_{0j})) \right] |_{r = 0} & = 0 \quad for all \alpha \in \Gamma.
\end{align*}
• (Double robustness) Moreover, \[\mathrm{E}_{P}\left[ \psi_j(Z, \theta_{0j}, \gamma, \alpha) \right] = - \mathrm{E}_{P}\left[ (\alpha(W) - \alpha_{0j}(W))(\gamma(W) - \gamma_{0j}(W)) \right]\] so that the augmented score is mean zero for all $\alpha \in \Gamma$ whenever $\gamma = \gamma_{0j}$, or for all $\gamma \in \Gamma$ whenever $\alpha = \alpha_{0j}$.
\end{enumerate}
proofSee, e.g., chernozhukov2018automatic.
remark(ref) (ii) hints at a trade-off between the quality of the estimates for $\gamma_{0j}$ and $\alpha_{0j}$. In situations where $\gamma_{0j}$ can be estimated very well, it may be possible to achieve $\sqrt{n}$-convergence and asymptotic normality even when the rate of convergence of $\hat{\alpha}_j$ is slow, and vice versa.
Consider a simultaneous confidence band for $\{\theta_{0j}\}_{j=1, \, \ldots, \, p}$ constructed as $\times_{j=1}^p [\hat{\theta}_j \mp n^{-1/2}\hat{\sigma}_jc_{\alpha}]$, where the point estimates $\{\hat{\theta}_{j}\}_{j=1, \, \ldots, \, p}$ are based on an empirical analog of (ref), \[\hat{\theta}_{0j} = \mathbb{E}_{n}\left[ m_j(W, \hat{\gamma}_{j}) + \hat{\alpha}_{j}(W)(Y - \hat{\gamma}_{j}(W)) \right], \quad j = 1, \, \ldots, \, p,\]
and $\hat{\sigma}_{j}^2$ is an estimator of $\sigma_j^2 = \mathrm{Var}(\sqrt{n}(\hat{\theta}_{j} - \theta_{0j}))$. The critical value $c_{\alpha}$ will be chosen such that \[\mathrm{P}_P\left( \max_{1 \leq j \leq p} \sqrt{n}\frac{|\hat{\theta}_{j} - \theta_{0j}|}{\hat{\sigma}_j} \leq c_\alpha \right) \approx 1 - \alpha,\] using a normal approximation for the sup-$t$-statistic. The goal of this section is to provide finite-sample guarantees for this normal approximation in the form of a bound on the Kolmogorov distance between the finite-sample distribution of the sup-$t$-statistic and $\max_{1 \leq j \leq p} Z_j$ for a suitable jointly Gaussian distribution $(Z_1, \ldots, Z_p) \sim \mathcal{N}(0, \Sigma)$: \[\sup_{t \in \mathbb{R}} \left| \mathrm{P}_P\left( \max_{1 \leq j \leq p} \sqrt{n}\frac{(\hat{\theta}_j - \theta_{0j})}{\sigma_j} \leq t \right) - \mathrm{P} \left( \max_{1 \leq j \leq p} Z_j \leq t \right)\right|.\]
Towards that goal, we list a set of sufficient regularity conditions in the following assumptions.
assumption[Moment conditions]
Suppose the following moment conditions hold:
\begin{enumerate}[(i)]
• (Bounded heteroskedasticity of the outcome) $\mathrm{E}_{P}\left[ (Y - \gamma_{0j}(W))^2 \mid W \right] \leq \Bar{\sigma}$ for all $j = 1, \, \ldots, \, p$.
• (Variance bounded away from 0) Let $\overline{\psi}_{0j}(Z) = m_j(W, \gamma_{0j}) + \alpha_{0j}(Y - \gamma_{0j}(W)) - \theta_{0j}$ denote the oracle score (that is, the score evaluated at the true value of the parameters). We assume $\sigma_j^2 = \mathrm{E}_{P}\left[ \overline{\psi}_{0j}(Z)^2 \right] \geq \sigma_{\min}^2 > 0$ for all $j = 1, \, \ldots, \, p$. Moreover, let $\Sigma$ denote the correlation matrix of the $\psi_{0j}(Z)$, with $(j,k)$-th entry given by \[\Sigma_{jk} = \mathrm{E}_{P}\left[ \frac{\overline{\psi}_{0j}(Z) \overline{\psi}_{0k}(Z)}{\sigma_j \sigma_k} \right] \quad j,k = 1, \, \ldots, \, p.\] We assume that the smallest eigenvalue of $\Sigma$ is bounded below by some $\lambda_{\min} \geq 0$.
• (Higher-order moments) For some $q \geq 4$, there exists $b_n < \infty$ such that \[\Vert \max_{1 \leq j \leq p} |\overline{\psi}_{0j}(Z)/\sigma_j|\Vert_{P,q} \leq b_n,\] and, for all $j = 1, \, \ldots, \, p$, \[\mathrm{E}_{P}\left[ (\overline{\psi}_{0j}(Z)/\sigma_j)^4 \right] \leq b_n^2.\]
\end{enumerate}
remark[On the eigenvalue condition] The assumption that the minimum eigenvalue of $\Sigma$ is bounded below allows us to obtain nearly-optimal rates with respect to the sample size in the normal approximation we use chernozhukov2021nearly. In practice, it implies that the identifying moments do not become perfectly correlated as the number of parameters grows. This restriction precludes certain applications, e.g., using a grid of $p$ points to approximate the CDF of an outcome, with the grid becoming dense asymptotically. The case where there can possibly be a continuum of parameters will be covered in the next section.
assumption[Nuisance parameters] Suppose the following:
\begin{enumerate}[(i)]
• (RMSE convergence rates) We have $\Vert \hat{\gamma}_{j} - \gamma_{0j} \Vert_{P,2} \leq \mathcal{R}_n(\hat{\gamma})$ and $\Vert \hat{\alpha}_{j} - \alpha_{0j} \Vert_{P,2} \leq \mathcal{R}_n(\hat{\gamma})$ for all $j = 1, \, \ldots, \, p$.
• (Boundedness of Riesz Representer) We have $\Vert \alpha_{0j} \Vert_{\infty} \leq \Bar{\alpha}$, $\Vert \hat{\alpha}_{j} \Vert_{\infty} \leq \Bar{\alpha}$ for all $j = 1, \, \ldots, \, p$.
• (Envelope and entropy conditions) The class of functions
\begin{align*}
\mathcal{F} = \{ (Y,W) \mapsto & (m_j(W, \gamma) + \alpha(W)(Y - \gamma(W)) - m_j(W, \gamma_{0j}) - \alpha_{0j}(W)(Y - \gamma_{0j}(W)) \, : \\ & \Vert \gamma - \gamma_{0j} \Vert_{P,2} \leq \mathcal{R}_n(\hat{\gamma}), \Vert \alpha - \alpha_{0j} \Vert_{P,2} \leq \mathcal{R}_n(\hat{\alpha}), j = 1, \, \ldots, \, p \}.
\end{align*} is suitably measurable, with a measurable envelope $F \geq \sup_{f \in \mathcal{F}} |f|$ that satisfies $\Vert F \Vert_{P,2+\delta} \leq M_n$ for some $\delta\geq 0$ and some sequence of constants $M_n$. There exist sequences $v_n \geq 1$, $a_n \geq n \vee M_n$, such that the uniform entropy numbers of $\mathcal{F}$ obey \begin{align*}
\log \sup_Q N(\varepsilon \Vert F \Vert_{Q,2}, \mathcal{F}, \Vert \cdot \Vert_{Q,2}) \leq v_n \log(a_n/\varepsilon), \qquad for all 0 < \varepsilon \leq 1.
\end{align*}
\end{enumerate}
remark[On entropy conditions]
The goal of entropy conditions is to control the complexity of the class of functions used to estimate nuisance parameters. On the one hand, this class needs to be rich enough for it to be possible to obtain good MSE convergence rates. On the other hand, if the class is too complex, it may lead to an overfitting bias when using the same data to estimate the nuisance parameters $\gamma_{0j}$, $\alpha_{0j}$ and the target parameter $\theta_{0j}$. An alternative to restricting the entropy of the class of functions considered is using some form of sample splitting, as discussed below in (ref).
The following is one of the main theoretical results of this paper. It provides a bound on the Kolmogorov distance between the finite-sample distribution of the sup-$t$-statistic of the DML estimators and $\max_{1 \leq j \leq p} Z_j$ for the Gaussian limit distribution $(Z_1, \ldots, Z_p) \sim \mathcal{N}(0, \Sigma)$, where $\Sigma$ is as defined in Assumption (ref).
theoremSuppose Assumptions (ref), (ref) and (ref) hold. Then,
\[\sup_{t \in \mathbb{R}} \left| \mathrm{P}_P\left( \max_{1 \leq j \leq p} \sqrt{n}\frac{(\hat{\theta}_j - \theta_{0j})}{\sigma_j} \leq t \right) - \mathrm{P} \left( \max_{1 \leq j \leq p} Z_j \leq t \right)\right| \leq \varrho(n, p),\] where \begin{align*}
\varrho(n, p) & = C(q) \left\lbrace \frac{b_n (\log p)^{3/2} \log n}{\sqrt{n} \lambda_{\min}} + \frac{b_n^2 (\log p)^{2} \log n}{n^{1-2/q} \lambda_{\min}} + \left[\frac{b_n^q (\log d)^{3q/2-4} \log n \log (pn)}{n^{q/2 - 1} (\lambda_{\min})^{q/2}} \right]^\frac{1}{q-2}\right\rbrace \tag{A}\\ & + \frac{6\sqrt{\log p}}{\sigma_{\min}} \left\lbrace \Delta_{1n} + \Delta_{2n} \right\rbrace \tag{B} \\ & + \frac{c}{\log n}
\tag{C} \\
\Delta_{1n} & = K\left(2+\delta,\frac{c}{3}\right)\Bigg( \left[\left((2+\sqrt{2})\Bar{\alpha} + \sqrt{2} \Bar{Q} \right)\mathcal{R}_n(\hat{\gamma}) + \Bar{\sigma}\mathcal{R}_n(\hat{\alpha}) \right]\sqrt{3v_n\log(3a_n)} \\ & \quad + 3v_n n^{\frac{1}{2+\delta}-\frac{1}{2}} 5M_n \log(3a_n) \Bigg), \\
\Delta_{2n} & = \sqrt{n} \mathcal{R}_n(\hat{\gamma}) \mathcal{R}_n(\hat{\alpha}),
\end{align*} for any $c > 0$, some constant $C(q) > 0$ depending on $q$, and some constant $K\left(2+\delta,\frac{c}{3}\right) > 0$ that may depend on $\delta$ and $c$.
proofThe full proof is in Appendix (ref). Here we discuss the heuristics, which may help understand each of the terms (A), (B) and (C).
The first step of the proof is a decomposition of $\hat{\theta}_j$ into an oracle estimator, \[\Bar{\theta}_j = \mathbb{E}_{n}\left[ m(W, \gamma_{0j}) + \alpha_{0j}(W)(Y - \gamma_{0j}(W)) \right]\] (i.e., the sample average of the augmented moment condition (ref) if $\gamma_{0j}$ and $\alpha_{0j}$ were known), and a deviation from that oracle estimator, $\mathbb{E}_n[m(W, \hat{\gamma}_j) + \hat{\alpha}_{j}(W)(Y - \hat{\gamma}_{j}(W)) - m(W, \gamma_{0j}) - \alpha_{0j}(W)(Y - \gamma_{0j}(W))]$.
On the one hand, the oracle estimator satisfies a high-dimensional version of a Berry-Esseen type of inequality, which allows us to quantify the Kolmogorov distance between its finite-sample distribution and the distribution of the corresponding multivariate normal distribution. In particular, we use the nearly-optimal rates in chernozhukov2021nearly, which yield term (A). We can obtain refinements on this term by assuming stronger conditions on the higher-order moments of $\overline{\psi}_{0j}(Z)$, as discussed in (ref).
On the other hand, the deviation from the oracle estimator can be bound using empirical process techniques. In the Appendix, we show that, with probability no more than $c/\log n$ (C), the deviation is upper bounded by $\Delta_{1n} + \Delta_{2n}$ (B). Improvements on (B) can be obtained by using some form of sample splitting, as discussed in (ref).
remark[Stronger tail conditions]
We could obtain a simpler bound for term (A) by assuming stronger tail conditions than the ones in (ref) (iii), as made clear in (ref), which collects results from chernozhukov2021nearly. In particular,
\begin{enumerate}[(i)]
• Assuming that $\overline{\psi}_{0j}(Z)$ is sub-Gaussian with Orlicz norm upper-bounded by $b_n$, (A) could be replaced by:
\[C\left\lbrace \frac{b_n (\log p)^{3/2} \log n}{\sqrt{n} \lambda_{\min}} + \frac{b_n^2 (\log p)^{2}}{\sqrt{n\lambda_{\min}}} \right\rbrace\] for an absolute constant $C > 0$.
• Assuming that $\overline{\psi}_{0j}(Z)$ is almost-surely bounded by $b_n$, (A) could be replaced by:
\[C\frac{b_n (\log p)^{3/2} \log n}{\sqrt{n} \lambda_{\min}}\] for an absolute constant $C > 0$.
\end{enumerate}
These stronger assumptions may be satisfied in certain economic applications, for example if outcomes are binary or naturally bounded (e.g., hours worked in a labor supply example).
remark[Removing entropy conditions by sample splitting]
The role of the entropy conditions in (ref) is to prevent overfitting bias, due to the same data being used in the nuisance parameter estimators $\hat{\gamma}_j$, $\hat{\alpha}_j$ and the estimator of the target parameter $\hat{\theta}_j$. Another way to overcome this overfitting problem is, as discussed in chernozhukov2018double or newey2018cross, to use sample splitting.
\begin{enumerate}
• With pure sample splitting, observations 1 to $n$ are divided randomly into $L$ data folds of roughly equal size, $I_\ell$, $\ell = 1, \, \ldots, L$. For a given $\ell$, estimators $\hat{\gamma}_{j\ell}$, $\hat{\alpha}_{j\ell}$ of $\gamma_{0j}$ and $\alpha_{0j}$ are constructed using the data not in $\ell$, $I_{\ell}^c$. An estimator for $\theta_{0j}$ is then constructed as: \[\hat{\theta}_{j} = \frac{1}{n} \sum_{\ell = 1}^L \sum_{i \in I_\ell} [m_j(W_i, \hat{\gamma}_{j\ell}) + \hat{\alpha}_{j\ell}(W_i)(Y_i - \hat{\gamma}_{j\ell})].\]
In that case, within the $\ell$-th fold and after conditioning on $I_{\ell}^c$, the summands are i.i.d., and so we can set $v_n = 1$, $a_n = e$ in (ref) (iii).
• An alternative is to use a “dirty” version of sample splitting, in which data in the $\ell$-th fold is used to select amongst a finite set of estimators $\{(\hat{\gamma}^{(1)}_{j\ell}, \hat{\alpha}^{(1)}_{j\ell}), \ldots, (\hat{\gamma}^{(r)}_{j\ell}, \hat{\alpha}^{(r)}_{j\ell})\}$ trained on data not in $\ell$. In that case, because a covering number for a finite class of functions is at most its cardinality, we can set $v_n = 1$, $a_n = e \vee r$. \qedhere
\end{enumerate}
Finally, the asymptotic validity of the simultaneous confidence band follows as a corollary of (ref) under two additional assumptions. This is not a new result, as it could be seen as a particular case of belloni2018uniformly, but we present it here for completeness.
assumption[Consistent variance estimation] Suppose that we have an estimator $\hat{\sigma}_j$ of $\sigma_j$ for all $j = 1, \, \ldots, \, p$ such that $\max_{1\leq j \leq p} (\hat{\sigma}_j/\sigma_j) \overset{p}{\rightarrow} 1$.
assumption[Growth conditions] Suppose the following growth conditions as $n, p \rightarrow \infty$:
\begin{enumerate}[(i)]
• (Nuisance parameters converge fast enough) $\sqrt{\log (p) n} \mathcal{R}_n(\hat{\gamma})\mathcal{R}_n(\hat{\alpha}) \rightarrow 0$.
• (Complexity characteristics do not grow too fast) \[\sqrt{\log (p) v_n\log(a_n)} [\mathcal{R}_n(\hat{\gamma})\vee\mathcal{R}_n(\hat{\alpha})] \rightarrow 0 \quad \text{and} \quad \sqrt{\log (p)} v_n n^{\frac{1}{2+\delta} - \frac{1}{2}}M_n \log(a_n) \rightarrow 0.\]
• (Moment bounds do not grow too fast)
\[\frac{b_n (\log p)^{3/2} \log n}{\sqrt{n} \lambda_{\min}} + \frac{b_n^2 (\log p)^{2} \log n}{n^{1-2/q} \lambda_{\min}} + \left[\frac{b_n^q (\log d)^{3q/2-4} \log n \log (pn)}{n^{q/2 - 1} (\lambda_{\min})^{q/2}} \right]^\frac{1}{q-2} \rightarrow 0.\]
\end{enumerate}
remarkAs we pointed out in (ref), there is a tradeoff between the RMSE convergence rate of $\hat{\gamma}$ and of $\hat{\alpha}$, which is made explicit in (i). In the case of a single parameter, $p = 1$, a sufficient condition for (i) is that both $\hat{\gamma}$ and $\hat{\alpha}$ converge faster than $n^{-1/4}$, a rate that is typically attainable by non-parametric estimators chernozhukov2018double. Note that the dependence of $p$ in the bound is only logarithmic, allowing for very high dimensional cases (potentially $p \gg n$).
corollary[Validity of the simultaneous confidence band]
Under Assumptions (ref) to (ref), we have \[\mathrm{P}_P\left( \hat{\theta}_j - n^{-1/2}\hat{\sigma}_j c_{\alpha} \leq \theta_{0j} \leq \hat{\theta}_j + n^{-1/2}\hat{\sigma}_j c_{\alpha}, \, \forall 1 \leq j \leq p \right) \rightarrow 1 - \alpha.\]
Continua of Parameters
In this section, we consider the same setting as belloni2017program. Again, suppose we have access to an i.i.d. sample $\{W_i\}_{i=1}^n$ from a probability distribution $P$ on $\mathcal{W}$. Now, we are interested in constructing a simultaneous confidence band for the set of (scalar) parameters $\{\theta_{0u}\}_{u \in \mathcal{U}}$, satisfying the following moment condition:
equation[equation omitted — 133 chars of source]
for a possibly uncountable set $\mathcal{U}$, where $\eta_{0u} \in \mathcal{T}_u$ is the unknown true value of an infinite-dimensional nuisance parameter. Here we also assume that an ML estimator $\hat{\eta}_u$ of $\eta_{0u}$ can be obtained from the same data.
Below we give a leading example where this framework may be appropriate.
example[Distributional treatment effects] Consider a setting where we want to evaluate the effect of a binary treatment $D \in \{0,1\}$ on an outcome $Y$. We work with observational data and, as in Examples $\ref{ex:1}$ and $\ref{ex:2}$, we suppose that we have access to a rich enough set of controls $X$ such that an unconfoundedness assumption holds.
In some cases, features of the marginal distributions of potential outcomes beyond the mean may be of interest. For example, policymakers may care about how the treatment impacts inequality. In other cases, economic theory will make predictions about how different regions of the outcome distribution should be affected by the treatment, so looking at features other than the mean can be used to probe or validate the theory.
Under the unconfoundedness assumption, the marginal distribution of $Y(d)$ can be identified as \[\theta_{0u} = F_{Y(d)}(u) = \mathrm{E}_{P}\left[ \gamma_{0u}(d, X) \right],\] where $\gamma_{0u}(d, x) = F_{Y}(u \mid D = d, X = x) = \mathrm{P}_P\left( Y \leq u \mid D, X \right)$ is the conditional distribution of $Y$ given $D$ and $X$ at point $u$. A non-parametric estimator of $\gamma_{0u}$ may be constructed using different techniques, for example, distribution regression chernozhukov2013inference.
Having access to estimates of $\{\theta_{0u}\}_{u \in \mathcal{U}}$ for a suitable range $\mathcal{U}$ allows to construct many interesting statistics, such as quantile treatment effects, Gini indices, and Oaxaca-Blinder type of decompositions.
Again, our goal in this section is to give finite-sample guarantees for a simultaneous confidence band for $\{\theta_{0u}\}_{u \in \mathcal{U}}$ using a set of DML point estimates $\{\hat{\theta}_{u}\}_{u \in \mathcal{U}}$ based on an empirical analog of (ref). We assume that each $\theta_{0u} \in \Theta_u$ for some $\Theta_u \subset \mathbb{R}$, and that we can find, for each $u \in \mathcal{U}$, an approximate solution to the empirical analog of (ref), i.e., a $\hat{\theta}_u$ such that
equation[equation omitted — 113 chars of source]
for some sequence of $\epsilon_n > 0$ such that $n^{-1/2}\epsilon_n \rightarrow 0$ as $n \rightarrow \infty$. Again, we want to choose a critical value $c_{\alpha}$ such that \[\mathrm{P}_P\left( \sup_{u \in \mathcal{U}} \sqrt{n}\frac{|\hat{\theta}_{u} - \theta_{0u}|}{\sigma_u} \leq c_\alpha \right) \approx 1 - \alpha,\] using a normal approximation for the sup-$t$-statistic. As in the previous section, our objective is to provide finite-sample guarantees for this normal approximation in the form of a bound on the Kolmogorov distance between the finite-sample distribution of the sup-$t$-statistic and the supremum of a suitable Gaussian process. We begin by giving some sufficient regularity conditions towards this result.
assumption[Moment problem]
For all $u \in \mathcal{U}$, the following conditions hold.
\begin{enumerate}[(i)]
• The true parameter satisfies $\theta_{0u} \in \mathrm{int} \, \Theta_u$.
• The map $(\theta, \eta) \mapsto \mathrm{E}_{P}\left[ \psi_u(W, \theta, \eta) \right]$ is twice continuously Gateaux-differentiable on $\Theta_{u} \times \mathcal{T}_u$.
• (Neyman orthogonality) For $\Bar{r} \in [0,1)$, let $D_{\Bar{r}u}[\eta - \eta_{0u}]$ denote the Gateaux derivative map of $\mathrm{E}_{P}\left[ \psi_u(W, \theta, \eta) \right]$ with respect to $\eta$ at $(\theta_{0u}, \eta_{0u})$ in the direction $\eta - \eta_{0u}$, \[D_{\Bar{r}u}[\eta - \eta_{0u}] = \partial_r \mathrm{E}_{P}\left[ \psi_u(W, \theta_{0u}, \eta_{0u} + r(\eta - \eta_{0u})) \right] |_{r = \Bar{r}}.\] Then, $D_{0u}[\eta - \eta_{0u}] = 0$ for all $\eta \in \mathcal{T}_u$.
• (Bounded derivatives with respect to $\theta$) Let $J_{0u} = \partial_{\theta} \mathrm{E}_{P}\left[ \psi_u(W, \theta, \eta_{0u}) \right] |_{\theta = \theta_{0u}}$. Then $c_0 \leq |J_{0u}| \leq C_0$. Moreover, $|\partial_{\theta} \mathrm{E}_{P}\left[ \psi_u(W, \theta, \eta_{0u}) \right]| > c_1$ for all $\theta \in \Theta_u$.
• (Lipschitz-continuity at the true parameters) For all $\theta \in \Theta_u$ and $\eta \in \mathcal{T}_u$, \[\mathrm{E}_{P}\left[ (\psi_u(W, \theta, \eta) - \psi_u(W, \theta_{0u}, \eta_{0u}))^2 \right] \leq C_0 (|\theta-\theta_{0u}|\vee\Vert\eta - \eta_{0u}\Vert_{P,2})^\omega.\]
• (Bounded derivatives with respect to $\eta$) For all $r \in [0,1)$, $\theta \in \Theta_u$ and $\eta \in \mathcal{T}_u$, \[|\partial_r \mathrm{E}_{P}\left[ \psi_u(W, \theta, \eta_{0u} + r(\eta - \eta_{0u})) \right]| \leq B_{1n}\Vert\eta - \eta_{0u}\Vert_{P,2}. \]
• (Bounded second derivatives) For all $r \in [0,1)$, $\theta \in \Theta_u$ and $\eta \in \mathcal{T}_u$,
\[|\partial^2_r \mathrm{E}_{P}\left[ \psi_u(W, \theta_{0u} + r(\theta - \theta_{0u}), \eta_{0u} + r(\eta - \eta_{0u})) \right]| \leq B_{2n}(|\theta-\theta_{0u}|\vee\Vert\eta - \eta_{0u}\Vert_{P,2})^2. \]
\end{enumerate}
remark[On the Neyman orthogonality condition] As opposed to the previous section, here we take Neyman orthonality (iii) as a primitive condition, and hence we assume that $\eta_{0j}$ contains all nuisance parameters needed to make the moment functional Neyman-orthogonal. In the case of a linear, mean-square continuous functional of a regression, the same construction as in (ref) is valid, and so $\eta_{0u} = (\gamma_{0u}, \alpha_{0u})$ for $\alpha_{0u}$ the Riesz representer. This is true, for instance, in (ref). More generally, chernozhukov2018automatic discuss how to orthogonalize non-linear functionals, and chernozhukov2018double, belloni2018uniformly cover many other important models, such as conditional moment restrictions.
assumption[Nuisance parameters] Suppose the following:
\begin{enumerate}[(i)]
• (RMSE convergence rates) For all $u \in \mathcal{U}$ we have $\Vert\hat{\eta}_u - \eta_{0u}\Vert_{P,2} \leq \mathcal{R}_n(\hat{\eta})$.
• (Envelope and entropy conditions) The class of functions \[\mathcal{F} = \{W \mapsto \psi_u(W, \theta, \eta) \, : \, \theta \in \Theta_u, \eta \in \mathcal{T}_u, u \in \mathcal{U}\}\] is suitably measurable, with a measurable envelope $F \geq \sup_{f \in \mathcal{F}} |f|$ that satisfies $\Vert F \Vert_{P,2+\delta} \leq M_n$ for some $\delta\geq 0$ and some sequence of constants $M_n$. There exist sequences $v_n \geq 1$, $a_n \geq n \vee M_n$, such that the uniform entropy numbers of $\mathcal{F}$ obey \begin{align*}
\log \sup_Q N(\varepsilon \Vert F \Vert_{Q,2}, \mathcal{F}, \Vert \cdot \Vert_{Q,2}) \leq v_n \log(a_n/\varepsilon), \qquad for all 0 < \varepsilon \leq 1.
\end{align*} Finally, for all $f \in \mathcal{F}$, we have $c_0 \leq \Vert f \Vert_{P,2} \leq C_0$.
\end{enumerate}
assumption[Entropy and moments of the score at the truth] Suppose $\sigma_u^2 = J_{0u}^{-2} \mathrm{E}_P[\psi_u(W, \allowbreak \theta_{0u}, \eta_{0u})^2] \geq C_0^{-2}c_0^2$, and let $\overline{\psi}_{0u}(W) = -(\sigma_{u}J_{0u})^{-1} \psi_u(W, \allowbreak \theta_{0u}, \eta_{0u})$ denote the re-scaled score evaluated at the true values $(\theta_{0u}, \eta_{0u})$. Then:
\begin{enumerate}[(i)]
• (Envelope and entropy conditions) The class of functions \[\mathcal{F}_0 = \{W \mapsto \overline{\psi}_{0u}(W) \, : \,u \in \mathcal{U}\}\] is suitably measurable, with a measurable envelope $F_0 \geq \sup_{f \in \mathcal{F}_0} |f|$ that satisfies $\Vert F_0 \Vert_{P,q} \leq b_n$ for some $q\geq 4$ and some sequence of constants $b_n$. There exist sequences $V_n \geq 1$, $A_n \geq n$, such that the uniform entropy numbers of $\mathcal{F}_0$ obey \begin{align*}
\log \sup_Q N(\varepsilon \Vert F_0 \Vert_{Q,2}, \mathcal{F}_0, \Vert \cdot \Vert_{Q,2}) \leq V_n \log(A_n/\varepsilon), \qquad for all 0 < \varepsilon \leq 1.
\end{align*}
• (Moments) For all $f \in \mathcal{F}_0$ and $k = 3,4$, $\mathrm{E}_{P}\left[ |f(W)|^k \right] \leq C_0 b_n^{k-2}$.
\end{enumerate}
The following theorem is the second main result of this paper. Again, it provides a bound on the Kolmogorov distance between the finite-sample distribution of the sup-$t$-statistic of the DML estimators and its limiting distribution. In the statement of the theorem, $G_P$ denotes a tight mean-zero Gaussian process indexed by the class of functions in (ref), with covariance function $\mathrm{E}_{}\left[ G_P [\overline{\psi}_{0u}] G_P [\overline{\psi}_{0u'}] \right] = \mathrm{E}_{P}\left[ \overline{\psi}_{0u}(W) \overline{\psi}_{0u'}(W) \right]$ for all $u,u' \in \mathcal{U}$.
theoremSuppose Assumptions (ref), (ref) and (ref) hold. Then,
\begin{multline*}
\sup_{t \in \mathbb{R}} \left| \mathrm{P}_P\left( \sup_{u \in \mathcal{U}} \sqrt{n}\frac{(\hat{\theta}_u - \theta_{0u})}{\sigma_u} \leq t \right) - \mathrm{P} \left( \sup_{u \in \mathcal{U}} G_P [\overline{\psi}_{0u}] \leq t \right)\right| \leq \\ \kappa r_{1n} \left(\chi \sqrt{V_n \log(A_n b_n)} + \sqrt{1 \vee \log(1 / r_{1n})}\right) + r_{2n},
\end{multline*}
where $\kappa, \chi > 0$ are universal constants,
\[r_{1n} = c_0^{-1}\epsilon_n + \Delta_{1n} + \Delta_{2n} + \Delta_{3n}\] for:
\begin{align*}
\Delta_{1n} & = C_0^{-1}K\left(2+\delta,\frac{c}{2}\right)\Bigg(
\sqrt{C_0} [\mathcal{R}_n^\vee (\hat{\eta})]^{\omega/2}\sqrt{2v_n\log(2a_n)} + 2v_n n^{\frac{1}{2+\delta}-\frac{1}{2}} 2M_n \log(2a_n) \Bigg). \\
\Delta_{2n} & = C_0^{-1} \tfrac{1}{2} \sqrt{n} B_{2n} [\mathcal{R}_n^\vee (\hat{\eta})]^2. \\
\Delta_{3n} & = \frac{b_nL_n}{\gamma^{1/2} n^{1/2 - 1/q}} + \frac{(b_n)^{1/2}L_n^{3/4}}{\gamma^{1/2} n^{1/4}} + \frac{(b_nL_n^2)^{1/3}}{\gamma^{1/3} n^{1/6}} \\
\mathcal{R}_n^\vee (\hat{\eta}) & = \Big\lbrace c_1^{-1} n^{-1/2}\epsilon_n + c_1^{-1} n^{-1/2}K\left(2 + \delta, \frac{c}{2}\right) \left(C_0 \sqrt{v_n \log(a_n)} + v_n n^{\frac{1}{2+\delta} - \frac{1}{2}}M_n\log(a_n) \right) \\
& + c_1^{-1} B_{1n}\mathcal{R}_n(\hat{\eta})\Big\rbrace \vee \mathcal{R}_n(\hat{\eta}), \\
L_n & = d(q) V_n(\log n \vee \log(A_n b_n)),
\end{align*}
and \[r_{2n} = D(q)\left(\gamma + \log n / n\right) + c/\log n,\] where $d(q),D(q)$ are constants depending only on $q$, for any $c>0$ and $\gamma \in (0,1)$.
proofThe full proof is in Appendix (ref). As above, we discuss the heuristics here. First, $\mathcal{R}_n^\vee (\hat{\eta})$ is the maximum of two objects: the rate $\mathcal{R}_n(\hat{\eta})$ for $\hat{\eta}$ and a preliminary rate for $\hat{\theta}$ (an upper bound for how far $\hat{\theta}$ can be from $\theta_0$ based only on the smoothness conditions). Notice that this step becomes unnecessary whenever the moment function $\psi_{u}$ is linear in $\theta$, which will be the case in many applications, including (ref).
Second, the Kolmogorov distance-based statement of the theorem is related by Lemma (ref) to another kind of approximation, of the form: \[\mathrm{P}_P\left( \left| \sup_{u \in \mathcal{U}} \sqrt{n}\frac{(\hat{\theta}_u - \theta_{0u})}{\sigma_u} - G_P [\overline{\psi}_{0u}] \right| > r_{1n} \right) \leq r_{2n}.\]
As an intermediate step, we first approximate $\sqrt{n}(\hat{\theta}_u - \theta_{0u})/\sigma_u$ by the empirical process on the re-scaled oracle score, $\mathbb{G}_n[\Bar{\psi}_{0u}]$. With high probability, the distance between these two objects is bounded by $\Delta_{1n} + \Delta_{2n}$. The first term quantifies the size of the deviation between $\psi_u(W, \hat{\theta}, \hat{\eta})$ and $\psi_u(W, \theta_{0u}, \eta_{0u})$ when $\psi_u(W, \hat{\theta}, \hat{\eta})$ is in the class of functions $\mathcal{F}$. The second term bounds the error that we incur by linearizing the score. In particular, if the score is linear in both $\theta$ and $\eta$, this term can be ignored.
Finally, the term $\Delta_{3n}$ arises when approximating the supremum of the empirical process on the oracle score, $\mathbb{G}_n[\Bar{\psi}_{0u}]$, by the supremum of the corresponding Gaussian process, $G_P[\Bar{\psi}_{0u}]$, and it is a consequence of (ref).
As in the previous section, we give two additional conditions for the asymptotic validity of the simultaneous confidence band. This is also not a new result, but was shown in belloni2017program, belloni2018uniformly. We present it below for completeness.
assumption[Consistent variance estimation] Suppose that we have an estimator $\hat{\sigma}_u$ of $\sigma_u$ for all $u \in \mathcal{U}$ such that $\sup_{u \in \mathcal{U}} (\hat{\sigma}_u/\sigma_u) \overset{p}{\rightarrow} 1$.
assumption[Growth conditions] Suppose the following growth conditions as $n \rightarrow \infty$:
\begin{enumerate}[(i)]
• (Nuisance parameters converge fast enough) $\sqrt{n} [\mathcal{R}_n(\hat{\eta})]^2 \rightarrow 0$.
• (Complexity characteristics do not grow too fast) \[\sqrt{v_n\log(a_n)} [\mathcal{R}_n^\vee (\hat{\eta})]^{\omega/2}, \quad v_n n^{\frac{1}{2+\delta} - \frac{1}{2}}M_n \log(a_n) \rightarrow 0 \quad \text{and} \quad r_{1n} \sqrt{V_n \log(A_n b_n)} \rightarrow 0.\]
• (Moment bounds do not grow too fast)
\[\frac{b_nL_n}{\gamma^{1/2} n^{1/2 - 1/q}} + \frac{(b_n)^{1/2}L_n^{3/4}}{\gamma^{1/2} n^{1/4}} + \frac{(b_nL_n^2)^{1/3}}{\gamma^{1/3} n^{1/6}} \rightarrow 0.\]
\end{enumerate}
corollary[Validity of the simultaneous confidence band]
Under Assumptions (ref) to (ref), we have \[\mathrm{P}_P\left( \hat{\theta}_u - n^{-1/2}\hat{\sigma}_u c_{\alpha} \leq \theta_{0u} \leq \hat{\theta}_u + n^{-1/2}\hat{\sigma}_u c_{\alpha}, \, \forall u \in \mathcal{U} \right) \rightarrow 1 - \alpha.\]
Conclusions
In many applications, researchers are interested in the causal impact of a treatment or policy that was not randomly assigned. Inference in such non-experimental settings is still possible by controlling for a set of covariates, conditional on which the treatment becomes plausibly exogenous. In modern settings, DML offers an alternative way to leverage a large number of potential controls with regularization and model selection, which does not bias the estimates of the target parameters thanks to a Neyman orthogonality condition. Often, we want to make joint inference on a high-dimensional set of parameters, such as the ATE of many treatments or combinations thereof, the ATE of a treatment on many outcomes, or effects on the entire marginal distribution of the potential outcomes. In this paper we have complemented existing asymptotic results for high-dimensional DML belloni2017program, belloni2018uniformly with finite-sample guarantees. These finite sample guarantees can be useful to applied researchers, as they are informative of how far off the coverage of simultaneous confidence bands can be from the nominal levels.
There is one natural extension of the paper that would be an interesting avenue for future research. Our guarantees are based on the true standard error, $\sigma_j$ or $\sigma_u$, respectively, which in general will not be known and will have to be estimated. For our asymptotic corollaries, we have simply assumed that estimators $\hat{\sigma}_j$ or $\hat{\sigma}_u$ exist. We leave it to future work to provide finite-sample guarantees for variance estimation in the high dimensional case, although we note that such guarantees are available for one-dimensional DML with sample splitting chernozhukov2021simple.