De-Biased Machine Learning of Global and Local Parameters Using Regularized Riesz Representers

De-Biased Machine Learning of Global and Local Parameters Using Regularized Riesz Representers

Introduction

Many statistical objects of interest can be expressed as a linear functional of a regression function (or projection, more generally). Examples include global parameters: average treatment effects, policy effects from changing the distribution of or transporting regressors, and average directional derivatives, as well as their local versions defined by taking averages over regions of shrinking volume. This variety of important examples motivates the problem of learning linear functionals of regressions. Global parameters are typically regular (estimable at $1/\sqrt{n}$ rate), and local parameters are non-regular (estimable at slower than $1/\sqrt{n}$ rates). Global parameters can also be non-regular under weak identification (for example, in average treatment effects, when propensity scores accumulate mass near zero or one, along a given sequence of models).

Often the regression is high dimensional, depending on many variables such as covariates in a treatment effect model. Plugging a machine learner into a functional of interest can give a badly biased estimator. To avoid such bias, we use debiased/“double” machine learning (DML) based on Neyman orthogonal scores that have zero derivative with respect to each first step learner (e.g., neyman1959,belloni2014pivotal,belloni2015uniform,chernozhukov2016locally,chernozhukov2018double,foster2019orthogonal). Note that the word “double” emphasizes the connection to double robustness, a property which orthogonal scores have in this case. Such scores are constructed by adding a bias correction term: the average product of the regression residual with a learner of the functional's Riesz representer (RR). This construction builds upon and is directly inspired by newey1994asymptotic where such scores arise in the computation of the semi-parametric efficiency bound for regular functionals. We also remove overfitting bias (high entropy bias) by using cross-fitting, an efficient form of sample splitting, where we average over data observations different than those used by the nonparametric learners. See, e.g, schick1986asymptotically for early use and chernozhukov2018double for more recent use in the context of debiased machine learning.

Using closed-form solutions for Riesz representers in several examples, chernozhukov2016locally,chernozhukov2018double defined DML estimators in high dimensional settings and established their good properties. In comparison, the new approach proposed in this paper has the following advantages and some limitations:

This paper is a revised version of \citet*{chernozhukov2018robins} that gave an algorithm based on $\ell_1$ regularization for automatically estimating the Riesz representer. This version is distinguished from \citet*{chernozhukov2018robins}, chernozhukov2018double, chernozhukov2016locally, and \citet*{chernozhukov2018learning} in covering local objects that are estimated at a rate slower than $1/\sqrt{n}$. Providing debiased machine learning for such local objects is an important contribution of this paper.

Sections (ref) and (ref) present the main ideas for a general audience. In Section (ref), we define global, local, and perfectly localized linear functionals of the regression, and provide orthogonal representations for these functionals. In Section (ref), we present two empirical examples: local and global average treatment effects, and local and global average derivatives.

Sections (ref) and (ref) are theoretical. In Section (ref), we provide estimation theory, demonstrating concentration and approximate Gaussianity of the DML estimator with regression and Riesz representer estimated via regularized moment conditions. We provide rates of convergence for estimating the Riesz representer, giving both fast rates under approximate sparsity and slow rates under the dense model. In Section (ref), we demonstrate asymptotic consistency and Gaussianity of the DML estimator with regression estimated via general machine learning.

The supplement provides supporting material. In Section (ref), we give a detailed account of how our work relates to previous and contemporary work. In Section (ref), we review prelimaries of functional analysis. In Section (ref), we analyze the structure of the leading examples, providing bounds on operator norm, variance of the score, and kurtosis. Finally, we provide proofs for each section.

Overview of target functionals, orthogonal representation, estimation, and inference

Target functionals

We consider a random element $W$ with distribution $P$ taking values $w$ in its support $\mathcal{W}$. Denote the $L^q(P)$ norm of a measurable function $f: \mathcal{W} \to \Bbb{R}$ and also the $L^q(P)$ norm of random variable $f(W)$ by $\| f \|_{P,q} = \| f(W)\|_{P,q}$. For a differentiable map $x \mapsto g(x)$, from $\mathbb{R}^d$ to $\mathbb{R}^k$, we use $\partial_{x'} g$ to abbreviate the partial derivatives $(\partial/\partial x') g(x)$, and we use $\partial_{x'} g(x_0)$ to mean $\partial_{x'} g (x) \mid_{x = x_0}$, etc. We use $x'$ to denote the transpose of a column vector $x$.

Let $(Y,X)$ denote a random sub-vector of $W$ taking values in their support sets, $y \in \mathcal{Y} \subset \Bbb{R}$ and $x \in \mathcal{X} \subset \Bbb{R}^{d_x}$, where $d_x=\infty$ is allowed. Let $F$ denote the law of $X$. We define $$x \mapsto \gamma_0^\star(x) := {\mathrm{E}} [Y\mid X=x],$$ as the unknown regression function of $Y$ on $X$. We consider the convex parameter space $\Gamma_0$ for $\gamma_0^\star$ with elements $\gamma$. (Later, in the theoretical sections, we generalize and replace the regression function by a projection).

Our goal is to construct high-quality inference methods for real-valued linear functionals of $\gamma_0^\star$. To present examples below we need to endow $\gamma_0^\star$ with a causal interpretation, which requires us to assume that it is a structural function, invariant to the changes in the distribution of $X$ under policies described below. This property is not guaranteed for an arbitrary regression problem. For the reader who is unfamiliar with these concepts, we note that a simple sufficient condition for invariance is follows: given a stochastic process $x \mapsto Y(x)$, called potential outcomes or structural function, vector $X$ is generated to follow distribution $F$ independently of $x \mapsto Y(x)$ and $Y$ is generated as $Y = Y(X)$. In this case we have $\gamma_0^\star(x) = {\mathrm{E}} Y(x)$ for any $F$. This condition is conventionally called exogeneity in econometrics and random assignment in statistics. The measurability requirement here is that $(x, \omega)\mapsto Y(x,\omega)$ is a measurable map. We refer to imbens2015causal, hernan2019causal, and peters2017elements for the relevant formalizations that enable causal interpretation.

Here and below, a weighting function is a measurable function $x \mapsto \ell(x) $ such that $\int \ell dF = 1$ and $\int \ell^2 dF < \infty$. In this example, setting

and so on. We can model small neighborhoods $N$ as shrinking in volume with the sample size. The local weighting and kernel weighting discussed below are applicable to all key examples. Moreover, they are combinable with other weighting functions so that, for example, we can target inference on local average treatment effects for the treated.

In this example, consider the case when $X=(D,Z)$ consists of continuous treatment variable $D$ and covariates $Z$. Further suppose $\ell(x)=\ell(d)$ and $t(x)=1$. Then the parameter of interest is $\theta_{0}^{\ast}={\mathrm{E}}[\ell(D)T(D)],$ where $T(d)={\mathrm{E}}[\partial_d\gamma_{0}^{\ast }(D,Z)|D=d].$ When $Y=Y(D)$ for a potential outcome process $Y(d)$ that is independent of treatment $D$ conditional on the covariates $Z$ and differentiable in $d,$ it was shown by altonji2005cross and florens2008identification that $T(d)={\mathrm{E}}[\partial_d Y(D)|D=d],$ which is an average treatment effect on the treated. Thus $\theta_{0}^{\ast}$ is a weighted average of the effect of treatment on the treated and would be equal to $T(d)$ for the perfectly localized $\ell(d)=1(D=d)/f_{D}(d),$ where $f_{D}(d)$ is the pdf of $D.$ Also for $\ell(d)\equiv1$, imbens2009identification showed that $\theta_{0}^{\ast}={\mathrm{E}}[T(D)]={\mathrm{E}}[\partial_d Y(D)]$, which is an average treatment effect. See also rothenhausler2019incremental.

In Example (ref), we consider the case where the variable of differentiation is also the variable of localization. As explained above, this case corresponds to effects of continuous treatments, and it turns out to require extra care in Section (ref). The other possible case is where the variable of differentiation is different than the variable of localization. Such a case turns out to be simpler and is handled by similar arguments as Examples (ref), (ref), and (ref) in Section (ref).

All of these statistical parameters play an important role in causal inference, counterfactual decompositions, and predictive analyses. Introduction of the weighting function $\ell(X)$ allows us to study subgroup effects and local effects, and these will be covered by our non-asymptotic results and asymptotic results. All of the above examples can be viewed as real-valued linear functionals of the regression function.

The linear operator $\gamma \mapsto \theta(\gamma)$ has the following generating function $m$ in these examples:

In these examples, we can recognize the dependency on the weighting function by writing $ m(w, \gamma; \ell). $ In Examples (ref), (ref), and (ref) we can decompose $m(w, \gamma;\ell) = m_0(w, \gamma)\ell(x).$

Estimation of some parameters of the form in Definition 2.1 is very straightforward, such as ${\mathrm{E}}[w(X)\gamma_0(X)]$ for a known function $w(x)$. These can be estimated as the sample mean of $w(X)Y$. Such simple estimation is not possible for the causal, counterfactual parameters in Examples (ref), (ref), (ref), and (ref). The approach of this paper provides estimators for these counterfactual parameters and can be used for many others.

Our local functionals are defined by using the weight function that localizes the functionals around value $d_0$ of a low-dimensional vector component $D$. Here $D$ is a $p_1$-dimensional component of vector $X$.

We consider the weighting function

where $K: \Bbb{R}^{p_1} \to \Bbb{R}$ is a kernel function of order $\mathsf{o}$ such that $\int K = 1$ and $$\int (\otimes^m u) K(u) du = 0, \quad \text{ for } m=1,..., \mathsf{o}-1,$$ with its support contained in the cube $[-1,1]^{p_1}$. The simplest example is the box kernel with $K(u) = \times_{j=1}^{p_1}1(-1 < u_j < 1)/2$, which is of order $\mathsf{o} = 2$.

To present the main results in the most clear way, we assume that $\ell_h$ is known, i.e. $\omega$ is known. Our main results also hold for one sided kernels. We leave to future work the application of this theory to settings with one sided limits, e.g. regression discontinuity design.

The difficulty in targeting localized functionals is that they are not pathwise differentiable. A key quantity in the analysis is the operator norm (the modulus of continuity) of $\gamma \mapsto \theta (\gamma)$ on $\Gamma$, defined as

We consider $L=\infty$ as non-regular cases, e.g. perfectly-localized functionals. We also consider cases where $L \to \infty$ as $n \to \infty$ as non-regular. Indeed, the latter case arises from approximating the functional with $L=\infty$ by functionals where $L \to \infty$, e.g. local functionals with $h \to 0$. The $L=\infty$ case also arises in triangular array asymptotics where $P$ changes with $n$. The asymptotic thought experiment, where $L \to \infty$, approximates non-asymptotic cases where $L$ is high. We emphasize that we derive both non-asymptotic results and their asymptotic corollaries (which lead to simplified statements conveying key qualitative features of non-asymptotic results).

Building an orthogonal representation of the target functional

Equation ((ref)) can be thought of as a direct formulation of the target parameter. Next we introduce a dual formulation and finally an orthogonal formulation. Towards this end, we define the Riesz representer $\alpha_0$.

A minimal linear representer exists if and only if $L< \infty$, as a consequence of the Riesz--Frechet theorem; see Lemma (ref) below. Therefore, when $L< \infty$, we define the following dual linear representation for the target parameter

To motivate the upcoming orthogonal representation, we note that either the direct or the dual identification strategies can be used for direct plug-in estimation, but this does not give good estimators, as explained in the following technical remark.

Therefore we proceed to construct another representation for $\theta^\star_0 $ that has the required Neyman orthogonality structure.

Unlike the direct or dual representations for the functional, this representation is Neyman orthogonal to perturbations $(\bar h, h)$ $\in \Gamma^2$ of $(\alpha_0^\star,\gamma_0^\star)$ such that

In fact, a stronger property holds

which implies ((ref)) as well as double robustness.

The Neyman orthogonality property states that the representation of the target parameter $\theta_0$ in terms of the nuisance parameters $(\alpha, \gamma)$ is invariant to the local perturbations of the values of the nuisance parameter. This property makes the orthogonal representation an excellent basis for constructing high quality point and interval estimators of $\theta^\star_0$ in modern high-dimensional settings when we will be plugging-in biased estimators in lieu of $\gamma_0^\star$ and $\alpha_0^\star$, where the bias occurs because of the regularization (see, e.g., chernozhukov2016locally,chernozhukov2018double).

Both $\gamma_0^*$ and $\alpha_0^*$ are identified, $\gamma_0^*$ as ${\mathrm{E}}[Y|X]$ and $\alpha_0^*$ by virtue of the consistent estimator we give in Section (ref). Identification allows us to use the orthogonal representation to estimate target parameters.

The case of finite-dimensional linear regression

It is instructive to consider the case of linear finite-dimensional regression. Consider $x \mapsto b(x) = \{ b_j(x)\}_{j=1}^p$ as a $p$-dimensional dictionary of basis functions with $b_j \in L^2(F)$ for each $j=1,...,p$. The regression function is assumed to obey the linear functional form $\gamma^\star_0 = b'\beta_0$ for some $\beta_0$. Also define $$ G= {\mathrm{E}} b(X)b(X)', \quad M = {\mathrm{E}} m(W, b). $$

First, observe that for $\gamma = b'\beta$, $$\theta (\gamma) = {\mathrm{E}} m(W, b'\beta) = {\mathrm{E}} m(W, b)' \beta = M \beta.$$ For instance, in Examples (ref), (ref), (ref), and (ref):

Second, we make a guess that the linear representer $\alpha^\star_0$ to be of the form $ \alpha^\star_0(x)= b(x)' \rho_0, $ for $\rho_0$ defined below. We can define the parameters $\beta_0$ and $\rho_0$ as any minimal $\ell_1$-norm solutions to the system of equations:

In particular, if $G$ is full rank, the solutions are $\beta_0= G^{-1} {\mathrm{E}} b(X) Y$ and $\rho_0 = G^{-1} M$.

We now verify the representation property for our guess: $$ {\mathrm{E}} \gamma(X) \alpha^\star_0(X) = {\mathrm{E}} \beta' b(X)b(X)'\rho_0 = \beta' G \rho_0 = \beta'M = \theta(\gamma), $$ for all $\beta$'s and hence all $\gamma$'s. The operator norm of $\theta (\gamma) = M'\beta$ is given by $$ L= \sup_{\beta \in \Bbb{R}^p\setminus \{0\}}\frac{ | M'\beta| }{\sqrt{\beta'G\beta}} = \sup_{\beta \in \Bbb{R}^p\setminus \{0\}} \frac{|\beta'G\rho_0|}{\sqrt{\beta'G\beta}} = \sqrt{\rho_0' G \rho_0} < \infty. $$

We conclude that direct, dual, and orthogonal representations are given by $$ \theta(\gamma) = M'\beta; \quad \theta(\alpha) = \rho' {\mathrm{E}} b(X)Y ; \quad \theta(\gamma, \alpha) = M'\beta + \rho' {\mathrm{E}} b(X)Y - \rho' G\beta, $$ where $\beta$ is $\gamma$'s parameter and $\rho$ is $\alpha$'s parameter. These representations appear to be both novel and useful.

The case of infinite-dimensional regression

In the infinite-dimensional case, we can employ the Riesz--Frechet representation theorem and Hahn--Banach extension theorem to establish existence of the linear Riesz representer.

The first part of the lemma shows (implicit) existence of a linear representer when $L < \infty$. Our estimation results will rely only on the existence of minimal representers. In some cases, however, we may utilize the closed-form solutions for linear representers (see, e.g., Section (ref) for the key examples), to improve the basis functions for estimating the minimal representers. There is also an efficiency reason to work with minimal representers rather than any linear representer, as highlighted in Section (ref) analyzing semi-parametric efficiency.

Informal preview of estimation and inference results

Our estimation and inference will exploit empirical analogs of both the orthogonal representation of the parameter ((ref)) and the equation defining the RR property ((ref)).

To approximate the regression function and the RR, we consider the $p$-vector of dictionary functions $b$, where the dimension $p$ of the dictionary can be large, potentially much larger than $n$. We approximate $\alpha^\star_0$ by a linear form $b'\rho_0$, and we approximate $\gamma^\star_0$ by a linear form $b'\beta_0$, and estimate the parameters using the algorithms below.

We remark that the RR estimator in step 2 is of Dantzig selector type, but is not exactly the Dantzig selector, requiring some new analysis. We use the GDS rather than series or spline estimation to accommodate high dimensional specifications for the regression and Riesz representer.

The dictionary $b(x)$ is very important for the GDS estimator. This dictionary should be chosen so that linear combinations of $b(x)$ can approximate in mean square any element of $\Gamma$. For example if $\Gamma$ is the set of linear combinations of an infinite sequence of regressors, as for a high dimensional regression, then $b(x)$ could be chosen as the first $p$ elements of that sequence. Also $p$ can be chosen flexibly, because $p$ will be allowed to grow faster than the sample size, as specified in the asymptotic theory to follow. In practice multiple choices of $p$ could be tried.

Next, we state the key concentration and approximate Gaussianity results informally. Key quantities in the analysis are the “true” score and its moments: $$ \psi_{0}^\star(W) := \theta_0^\star-m(W, \gamma_0^\star)-\alpha^\star _0(X)(Y- \gamma_0^\star(X)),\ \ \sigma^2 := {\mathrm{E}} \psi^2_{0}(W), \ \ \kappa^3 := {\mathrm{E}} |\psi^3_{0}(W)|. $$

We establish conditions under which

These include a bound on the $\ell_1$ norm of coefficients and that either the regression function or the RR is approximately sparse with the effective dimension $s$ less than $\sqrt{n}$. This allows either nuisance parameter to be completely dense.

Given that ((ref)) holds, we establish that the resulting de-biased (or “double”) machine learning (DML) estimator $\hat \theta$ is adaptive, namely it is approximated up to the error $o(\sigma/\sqrt{n})$ by the oracle estimator $$ \bar \theta := \theta^\star_0 - n^{-1} \sum_{i =1}^n \psi_0(W_i), $$ where the oracle knows the scores $\psi_0$. Hence the approximate deviation of $\hat \theta$ from $\theta_0^\star$ is determined by $\|\psi_{0}\|_{P,2}/\sqrt{n}$, which is the standard deviation of the oracle estimator.

Consequently, $\hat \theta$ concentrates in a $\sigma/\sqrt{n}$ neighborhood of the target with deviations controlled by the normal laws, $$ \sup_{t \in \Bbb{R} }\Big |{\mathrm{P}} ( \sqrt{n}(\hat \theta - \theta^\star_0)/ \sigma \leq t) - {\mathrm{P}}( N(0,1) \leq t) \Big | \leq A (\kappa/\sigma)^3/\sqrt{n} + \mathrm{error}_n \to 0, $$ where the $\mathrm{error}_n$ bound is non-asymptotic and tends to zero as $n \to \infty$. Of course, $\sigma/ \sqrt{n} \to 0$ is required for concentration. The non-asymptotic bound automatically implies the uniform validity of results over large classes of probability laws $P$ for $W$.

There are two cases to consider:

As we show in Section (ref), in the case of local functionals, the latter condition can be more succinctly stated as $$ (\kappa/\sigma) \lesssim \sigma \asymp L , \quad L/\sqrt{n} \to 0. $$

Finally, we establish that we can transfer learning and inference guarantees for local functionals to those for the (perfectly) localized functionals if the localization bias is sufficiently small, namely $$ \sqrt{n} (\theta(\gamma_0^\star; \ell_h) - \theta(\gamma_0^\star; \ell_0)/\sigma \to 0. $$

We think it is remarkable that a single inference theory covers both regular and non-regular cases, and provides uniform validity over large classes of $P$.

Applications

Global and local effects of 401(k) eligibility on net financial assets

First, we use our method to answer a question in household finance: what is the average treatment effect of 401(k) eligibility on net financial assets (over a horizon of about two years)? We follow the identification strategy of poterba1994 and poterba1995, who assume selection on observables. The authors assume that when 401(k) was introduced, workers ignored whether a given job offered 401(k) and instead made employment decisions based on income and other observable job characteristics; after conditioning on income and job characteristics, 401(k) eligibility was exogenous at the time. This empirical question corresponds to Example (ref).

We use data from the 1991 US Survey of Income and Program Participation (chernozhukov2018double_data), using sample selection and variable construction as in abadie2003semiparametric and chernozhukov2004effects. The outcome $Y$ is net financial assets defined as the sum of IRA balances, 401(k) balances, checking accounts, US saving bonds, other interest-earning accounts, stocks, mutual funds, and other interest-earning assets minus non-mortgage debt. The treatment $D$ is an indicator of eligibility to enroll in a 401(k) plan. The raw covariates $X$ are age, income, years of education, family size, marital status, two-earner status, benefit pension status, IRA participation, and home-ownership. We impose common support of the propensity score for the treated and untreated groups based on these covariates, yielding $n=9869$ observations. We consider the fully-interacted specification $b(D,X)$ of chernozhukov2018double with $p=277$ including polynomials of continuous covariates, interactions among all covariates, and interactions between covariates and treatment status.

Tables (ref) and (ref) summarize results for the entire population and for each quintile of the income distribution. We use $K=5$ folds in cross-fitting. To estimate the RR, we use the generalized Dantzig selector (GDS) procedure introduced in the present work. To estimate the regression, we use GDS, Lasso, random forest, or neural network. GDS is implemented using the tuning procedure described in Section (ref). Lasso is implemented using the tuning procedure described in \citet*{chernozhukov2018learning}. Random forest and neural network are implemented with the same settings as chernozhukov2018double, i.e. with 1000 trees or a single hidden layer of eight neurons, respectively. We find ATE of $7608$ $(1395)$ using GDS for both the RR and the regression. This ATE estimate is stable across different choices of regression estimator. We find that localized ATE is not statistically significant for the second quintile, and it is statistically significant, positive, and strongly heterogeneous for the other quintiles. Interpreting the relatively high effect of 401(k) eligibility for the first quintile is a question for future research.

For comparison, chernozhukov2018double report ATE of $7170$ $(1398)$ by DML, which estimates the RR by estimating the propensity score and plugging it into the RR functional form. Though these two estimators are asymptotically equivalent under correct specification, our estimator avoids the estimated propensity score in the denominator which could cause numerical instability. The ATE results are broadly consistent with poterba1994 and poterba1995, who use a simpler specification motivated by economic reasoning. The localized ATE estimates by income quintile group appear to be new empirical results and are of interest in their own right. In Section (ref) we report analogous estimates without debiasing. Without debiasing, the GDS and Lasso estimates of ATE are attenuated due to regularization. The bias is smaller for the random forest and neural network estimates of ATE.

Global and local price elasticity of gasoline demand

Second, we use our method to estimate the average price elasticity of household gasoline demand: the percentage change in demand due to a unit percentage change in price. This parameter is critical for assessing the welfare consequences of tax changes, and it has been studied in hausman1995nonparametric,schmalensee1999household,yatchew2001household,blundell2012measuring. Formally, the parameter of interest is the average derivative of log demand with respect to log price holding income and demographic characteristics fixed. The exact version of this empirical question corresponds to Example (ref). The approximate version of this empirical question corresponds to Example (ref).

We use data from the 1994-1996 Canadian National Private Vehicle Use Survey (semenova2020debiased_data), using sample selection and variable construction as in yatchew2001household and belloni2019conditional. The outcome $Y$ is log gasoline consumption. The variable $D$ with respect to which we differentiate is log price per liter. The raw covariates $X$ are log age, log income, and log distance as well as geographical, time, and household composition dummies. In total we have $n=5001$ observations. We consider the specification $b(D,X)$ previously considered by semenova2020debiased augmented with additional interactions. The semenova2020debiased specification includes polynomials of continuous covariates, and interactions of log price (and its square) with time and household composition dummies. We further include interactions of log price (and its square) with log age, log age squared, log income, and long income squared to allow for heterogeneity. Altogether, $p=99$.

Table (ref) summarizes results for the entire population and for each quintile of the income distribution. We use $K=5$ folds in cross-fitting. To estimate the RR, we use the generalized Dantzig selector (GDS) procedure introduced in the present work. To estimate the regression, we use GDS, Lasso, random forest, or neural network. Again, GDS is implemented using the tuning procedure described in Section (ref), Lasso is implemented using the tuning procedure described in \citet*{chernozhukov2018learning}, and random forest and neural network are implemented with the same settings as chernozhukov2018double. We find average price elasticity of $-0.28$ $(0.06)$ using GDS for both the RR and the regression. Lasso gives similar results. Note that random forest is not differentiable, and the derivative of a neural network may be difficult to extract from a black-box implementation. When using these estimators, we implement a partial difference approximation of the derivative, detailed in Section (ref). We conjecture that this approximation explains why the results using random forest appear attenuated and why the results using neural network appear positive or statistically insignificant. Using GDS, we find that localized average price elasticity is statistically significant and negative in each income quintile, with substantial heterogeneity.

For comparison, OLS regression of log consumption on log price, log age, log income, and log distance as well as geographical, time, and household composition dummies yields an estimate of $0.14$ $(0.06)$. The linear specification leads to a positive elasticity estimate, contradicting economic intuition (since it says there would be more gasoline consumption when prices are higher). Our localized average price elasticity results using GDS are broadly consistent with semenova2020debiased, who more explicitly consider the relationship between average price elasticity and income. In Section (ref) we report analogous estimates without debiasing. Without debiasing, the GDS and Lasso estimates of quintile elasticities are attenuated due to regularization. The bias is smaller for the random forest and neural network estimates of quintile elasticities.

Estimation and inference for high dimensional approximately linear models

Best linear approximations for the regression function and the Riesz representer

To approximate the regression function, we consider the $p$-vector of dictionary functions $$x \mapsto b(x) = (b_j(x))_{j=1}^p, \quad b_j \in L^2(F).$$ The dimension $p$ of the dictionary can be large, potentially much larger than $n$. Let $\Gamma_b$ be the linear subspace of $L^2(F)$ generated by $b$. We assume that as $n \to \infty$ we have that $p \to \infty$ and $ \Gamma_b \to \bar \Gamma:= \text{closure}(\Gamma), $ where $\bar \Gamma$ is a linear subspace of $L^2(F)$ with the basis functions $\{\tilde b_j\}_{j=1}^{\infty}$. Here convergence means that any convergent sequence in $\Gamma_b$ has its limit in $\bar \Gamma$ and for each $\gamma \in \bar \Gamma$ we have a sequence in $\Gamma_b$ converging to it, with respect to the $L^2(F)$ norm. Note that this setup allows the dictionary $b= b_n$ to change with $n$, as for example with $b$-splines.

Here we define $\gamma^\star_0$ as a projection of $Y$ onto $\bar \Gamma$, i.e. $\gamma^\star_0$ is the projection of $Y$ on the infinite set of variables $\{\tilde b_j(X)\}_{j=1}^{\infty}$. This setup is slightly more general than in the introduction, where $\gamma^\star_0$ was the conditional expectation function. Of course, if the latter is an element of $\bar \Gamma$, it automatically coincides with $\gamma^\star_0$.

We approximate $\gamma_0^\star$ by the finite-dimensional best linear predictor (BLP) $\gamma_0$ via $$ \gamma^\star_0 = \gamma_0 + r_\gamma := b'\beta_0 + r_\gamma: {\mathrm{E}} [ b(X) r_\gamma(X) ] =0, $$ where $r_\gamma$ is the approximation error, and $ \gamma_0 := b'\beta_0$ is the best linear predictor of $Y$ and best linear approximation to $\gamma^\star_0$. We define $\beta_0$ as a minimal $\ell_1$-norm solution to the system of equations $$ \min \|\beta\|_1 : {\mathrm{E}} [ b(X) ( \gamma^\star_0(X) - b(X)'\beta)] = 0, $$ when $G = {\mathrm{E}} b(X)b(X)'$ is not full rank.

Similarly, we approximate the Riesz representer $\alpha^\star_0$, which exists by Lemma (ref) whenever $L<\infty$, via the best linear approximation $\alpha_0$: $$ \alpha^\star_0 = \alpha_0 + r_\alpha = b'\rho_0 + r_\alpha: {\mathrm{E}} [r_\alpha (X) b(X) ] = 0. $$ We define $\rho_0$ as a minimal $\ell_1$-norm solution to the system of equations $$ \min \| \rho \|_1 : {\mathrm{E}} [(\alpha^\star_0(X) - b(X)'\rho)b(X)] = 0. $$ Using that ${\mathrm{E}} \alpha_0^\star(X) b(X) = {\mathrm{E}} m(W, b)$, we note that

Hence $\alpha_0$ is the Riesz representer for ${\mathrm{E}} m(W,\gamma)$ for each $\gamma \in \Gamma_b$. Here we can interpret $\Gamma_b$ as the collection of test functions on which the representation property ((ref)) holds.

If the approximation errors are such that

our inference will target the ultimate parameter. In the non-regular setup, the second order error condition $\int r_{\alpha}r_{\gamma}dF\leq \sigma/\sqrt{n}$ is weaker than what is usually required for pathwise differentiable functionals (since $\sigma\rightarrow\infty$ is the non-regular case); there is a lower bar for oracle rates in non-regular problems. This phenomenon was also noted by foster2019orthogonal and kennedy2020optimal. Otherwise our inference will target an interpretable penultimate parameter. We shall formally refer to the latter case as the misspecified case.

The second claim of the lemma is immediate from the definition of $(\beta_0, \rho_0)$ and the first follows from elementary calculations. The orthogonality property above says that the score function is invariant to small perturbations of the nuisance parameters $\rho$ and $\beta$ around their “true values” $\rho_0$ and $\beta_0$. This invariance property plays a crucial role in removing the impact of biased estimation of nuisance parameters $\rho_0$ and $\beta_0$ on the estimation of the main parameters $\theta_0$.

Estimators

Estimation will be carried out using the following Dantzig selector-type estimators (candes2007dantzig). In a follow-up work, \citet*{chernozhukov2018learning} consider Lasso-type estimators.

Here we record the possibility of convex restrictions on the parameter space by placing $t$ in a convex parameter space $T$. If parameter restrictions are correct, then this can potentially improve theoretical guarantees by weakening the requirements on $G$ and other primitives.

In this setting, our estimator specializes to the original Dantzig selector. In practice, we use $T_\beta = \Bbb{R}^p$, although when we are interested in average derivative functionals, it is theoretically helpful to impose the convex restrictions of the sort $T = \{ t \in \Bbb{R}^p: \sup_{x \in \mathcal{X}} | \partial_d b(x)'t | \leq B\}$, where $B$ is some a priori known upper bound on the derivative. Ideally, $D_{\beta}$ is chosen such that $\mathrm{diag}( Var( D^{-1}_{\beta} (\hat G \beta_0 - \hat M) ) = I$. Our practical algorithm given in Section (ref) estimates $D_\beta$ from the data.

In this setting, our estimator is a generalization of the original Dantzig selector. In practice, we are using $T_\rho =\Bbb{R}^p$, even though it is possible to exploit some structured restrictions on the problem motivated by the nature of the universal Riesz representers. Ideally, $D_{\rho}$ is chosen such that $\mathrm{diag}( Var( D^{-1}_{\rho} (\hat G \rho_0 - \hat M) ) = I$. Our practical algorithm given in Section (ref) estimates $D_\rho$ from the data.

We now define the DML estimator with Riesz Representers, which makes use of cross-fitting.

Properties of DML: Main result

We provide a single non-asymptotic result that allows us to cover both global and local functionals, implying uniformly valid rates of concentration and normal approximations over large sets of $P$.

Consider the oracle estimator based upon the true score functions: $$ \bar \theta := \theta_0 - n^{-1} \sum_{i =1}^n \psi_0(W_i), \quad \psi_0(W) :=\psi(W, \theta_0; \beta_0, \rho_0). $$ We seek to establish minimal conditions under which the DML estimator approximates the oracle estimator, and is approximately normal with distribution $$ N(0, \sigma^2/n), \quad \sigma := \| \psi_0\|_{P,2}. $$ For regular functionals $\sigma$ is bounded, giving $1/\sqrt{n}$ concentration around $\theta_0$, and for non-regular functionals $\sigma \propto L \to \infty$ requring $L/\sqrt{n} \to 0$ to get concentration. Our normal approximation is accurate if kurtosis of $\psi_0$ does not grow too fast: $$ (\kappa/\sigma)^3 /\sqrt{n} \text{ is small}, \ \ \kappa:= \| \psi_0\|_{P,3}. $$ In the regular case $(\kappa/\sigma)^3$ is bounded, but for the non-regular cases it can scale as fast as $L$, again requiring $L/\sqrt{n} \to 0$.

Fix all of these sequences and the constants. Define the guarantee set: $$ \mathsf{S} = \left\{

\right \}, $$ We will take $u=\hat{\beta}_k-\beta_0$ and $v=\hat{\rho}_k-\rho_0$. As such, $r_1$ measures the non-asymptotic mean square rate for the BLP; $r_2$ measures the non-asymptotic mean square rate for the RR; and $r_3$ measures how the estimation errors interact. Note the presence of $\sigma$ acting on $r_2$ and $r_3$, which accommodates non-regular functionals. We will instantiate $(r_1,r_2,r_3)$ as fast and slow rates by analyzing the GDS estimator, in Theorem (ref) below.

Next, define $\mu$ to be the smallest modulus of continuity such that on $(u,v) \in \mathsf{S} $ $$ \sqrt{Var} ( (-m(W, b) + {\rho_0}' b(X) b(X)) 'u) \leq \mu \sigma \| b'u\|_{P,2}, \sqrt{Var} ( (Y - b(X)'\beta_0) b(X) 'v) \leq \mu \| b'v\|_{P,2}, \\ $$$$ \sqrt{Var} (u 'b(X) b(X)' v) \leq \mu ( \| b'u\|_{P,2}+\| b'v\|_{P,2}). $$ In typical applications, the modulus of continuity $\mu$ is bounded. Indeed, if elements of the dictionary are bounded with probability one, $ \| b(X) \|_\infty \leq C$, then we can select $\mu = C B$ for many functionals of interest, so the assumption is plausible. If $b(X) = X$ are sub-Gaussian, then this assumption is also easily satisfied; however, this case is not of central interest to us. See chernozhukov2021simple for a more general discussion.

Consider $P$ that satisfies the following conditions.

$R(\delta)$ is a requirement on how the sequences $(r_1,r_2,r_3)$ evolve relative to $(\sigma,\mu,m)$. We will formally verify $R(\delta)$ for the approximately sparse setting, in Corollary (ref) below. $R(\delta)$ is the key condition for our main result, Theorem (ref).

The conclusions of this result are distinguished from those of chernozhukov2018double and \citet*{chernozhukov2018learning} in applying to local, nonparametric objects, in providing finite sample bounds, and in being uniform over the parameter space. The conclusions are similar to this previous work in relying on a rate condition that is the product of rates of estimation for two distinct functions, here the regression and the Riesz representer.

The constants can be chosen to yield an asymptotic result.

Hence the DML estimator of the linear functionals of the BLP function $\gamma_0$ enjoys good properties under the stated regularity conditions. This result does not distinguish between inference on global functionals from inference on local functionals, as long as the latter are not perfectly localized. We state a separate result for perfectly localized functionals below.

The approximation bias for the ultimate target can be plausibly small due to the fact that many rich function classes admit regularized linear approximations with respect to conventional dictionaries $b$. For instance, tsybakov2012introduction and belloni2014pivotal show small approximation bias using Fourier bases as dictionaries, and using Sobolev and rearranged Sobolev balls, respectively, as the function classes.

Semi-parametric efficiency

Below we use concepts from semi-parametric efficiency, as presented in bickel1993efficient and vaart; we do not recall them here for brevity.

The DML estimator $\hat \theta$ will be asymptotically efficient for estimating $\theta^\star_0$, defined as a functional of $\gamma^\star_{0}$, the projection of $Y$ on $\bar{\Gamma}$. The distribution of a data observation is unrestricted in this case, so that there will only be one influence function for each functional of interest, and the estimator is asymptotically linear with that influence function. The standard semiparametric efficiency results then imply that our estimator will have the smallest asymptotic concentration among estimators that are locally regular; see bickel1993efficient and vaart.

Our formal result stated below only implies efficiency for the regular case, where the operator norm of the function $L$ is bounded, holding $P$ fixed. We expect that a similar result continues to hold with $L \to \infty$, by developing an appropriate formalization that handles $P$ changing with $n$ and rules out super-efficiency phenomena. However, this formalization requires a separate major development, which we leave to future research. In what follows, the notation $\gamma^\star_{0,P}$ emphasizes the dependence of the projection $\gamma^\star_0$ on $P$.

Properties of GDS estimators

Our goal is to verify that the guarantee $R(\delta)$ holds. In particular we have to analyze $(r_1,r_2,r_3)$ by bounding the population prediction norm $ v \mapsto \sqrt{v' G v}. $ This is a more nuanced problem than bounding the empirical prediction norm $v \mapsto \sqrt{v' \hat G v}$, which has been accomplished in a variety of prior analyses done on Dantzig-type and Lasso-type estimators.

We begin with the following condition, which only controls the max of error rates and controls the $\ell_1$ norm of true parameters:

The bounds on $\ell_1$ norm of coefficients are naturally motivated, for example, by working in Sobolev or rearranged Sobolev spaces (see, tsybakov2012introduction and belloni2014pivotal, respectively). Rearranged Sobolev spaces allow the first $p$ regression coefficients in the series expansion to be arbitrarily rearranged, allowing a much greater degree of oscillatory behaviors than in the original Sobolev spaces. The complexity of these function classes are also different. Sobolev spaces are Donsker sets under sufficient smoothness, whereas rearranged Sobolev spaces have the covering entropy bounded below by $\log p$ and are not Donsker if $p \to \infty$.

At the core of this approach is the restricted set $$ S (t_0, \nu) := \{ \delta: \| G \delta \|_\infty \leq \nu , \| t_0+ \delta\|_1 \leq \|t_0\|_1, t + \delta \in T\}, $$ where $\nu$ is the noise level. As demonstrated in the proof, the GDS estimator belongs to this set with high probability $1- \epsilon$ for the noise level $\nu= 4B \bar \lambda,$ where $\lambda$ is the penalty level of GDS ($\nu$ scales like $\sqrt{ \log (p \vee n) }/\sqrt{n}$ in our problems).

The effective dimension is defined in terms of the population (rather than sample) covariance matrix $G$, which makes it easy to verify regularity conditions. Note that if $G=I$ and $\|t_0\|_0 =s$, then $s(t_0) \leq s.$ More generally, $s(t_0)$ measures the effective difficulty of estimating $t_0$ in the prediction norm, created by design $G$ and the structure of $t_0$. The condition imposes no conditions on the restricted or sparse eigenvalues of $G$. For example, take $G = 11'$, a rank 1 matrix, and suppose $\|t_0\|_0 = 1$. Then $s(t_0) \leq 1$ holds in this case, giving useful and intuitive performance bounds, while the standard restricted eigenvalues and cone invertibility factors are all zero in this case, yielding no bounds on the performance in the population prediction norm. This type of example illustrates the possibility of accommodation of overcomplete (multiple or amalgamated) dictionaries in $b$, whose use in conjunction with $\ell_1-$ penalization has been advocated by donoho2005stable. Of course, the bounds on effective dimension follow from the bounds on cone-invertibility factors and restricted eigenvalues.

Given a vector $\delta \in \Bbb{R}^p$, let $\delta_A$ denote a vector with the $j$-th component set to $\delta_j$ if $j \in A$ and $0$ if $j \not \in A$.

The cone invertibility factor is a generalization of the restricted eigenvalue condition of bickel2009simultaneous, proposed by ye2010rate. The concept of the effective dimension does not split $t_0$ into a sparse component and a small dense component, as is done in the now standard analysis of $\ell_1$-regularized estimators of approximately sparse $t_0$. The effective dimension is simply stated in terms of $t_0$ alone.

Since approximate sparsity is a simple condition that implies a bound on effective dimension, we pause and interpret approximate sparsity in the context of a motivating example from causal inference. In particular, we revisit ATE (Example (ref)). For simplicity, consider the global parameter and assume that the function ${\mathrm{E}}[Y|D,Z]$ is an element of $\bar{\Gamma}$, so that $\gamma_0^{\star}(D,Z)={\mathrm{E}}[Y|D,Z]$ and $\alpha_0^{\star}(D,Z)=D/\pi_0^{\star}(Z)-(1-D)/(1-\pi_0^{\star}(Z))$ where $\pi_0^{\star}(Z)={\mathrm{E}}[D|Z]$ is the propensity score. Consider the dictionary $b(d,z)=(dq(z)',(1-d)q(z)')'$ where $\{q_j(z)\}_{j=1}^{p/2}$ are the initial $p/2$ elements of a sequence of basis functions that approximates the functions ${\mathrm{E}}[Y|1,Z]$, ${\mathrm{E}}[Y|0,Z]$, $1/\pi_0^{\star}(Z)$, and $1/(1-\pi_0^{\star}(Z))$.

Suppose the minimal $\ell_1$-norm mean square projections of ${\mathrm{E}}[Y|1,Z]$ and ${\mathrm{E}}[Y|0,Z]$ onto $\{q_j(z)\}_{j=1}^{p/2}$ are approximately sparse after rescaling appropriately by $D^{-1}_{\beta}$. (Note that if ${\mathrm{E}}[Y|1,Z]$ and ${\mathrm{E}}[Y|0,Z]$ are already approximately sparse then so are their projections.) It follows that the minimal $\ell_1$-norm mean square projection of $\gamma_0^{\star}$ is approximately sparse and $s_{\beta}:=s(D^{-1}_\beta \beta_0;\nu)$ is small.

Suppose instead that the minimal $\ell_1$-norm mean square projections of $1/\pi_0^{\star}(Z)$ and $1/(1-\pi_0^{\star}(Z))$ onto $\{q_j(z)\}_{j=1}^{p/2}$ are approximately sparse after rescaling appropriately by $D^{-1}_{\rho}$. (Note that if $1/\pi_0^{\star}(Z)$ and $1/(1-\pi_0^{\star}(Z))$ are already approximately sparse then so are their projections.) It follows that the minimal $\ell_1$-norm mean square projection of $\alpha_0^{\star}$ is approximately sparse and $s_{\rho}:=s(D^{-1}_\rho \rho_0;\nu)$ is small.

The bound is a minimum of what is called the “fast rate bound” and the “slow rate” bound. This result tightens the result in chatterjee2013assumptionless who established a “slow rate” bound (in the context of Lasso) that applies under no assumptions on $G$. If the effective dimension is not too big, as in the examples above, the “fast rate” $s(t_0;\nu) \nu^2$ provides a tighter bound under weak assumptions on $G$. It is important to emphasize that the result is stated in terms of the population prediction norm rather than the empirical norm.

We now apply this result to GDS estimators of the Riesz representer and the BLP. We impose the following conditions. Let $\Bbb{G}_A$ denote the empirical process over $f \in \mathcal{F}: \mathcal{W} \to \Bbb{R}^p$ and $ i \in A$, namely $$\mathbb{G}_A f := \mathbb{G}_A f (W) := | I | ^{-1/2} \sum_{ i \in A} (f(W_i) - P f), \quad Pf := P f(W) := \int f(w) d P(w).$$

The following is a sufficient condition that will deliver the guarantee $R(\delta)$ for $\delta \to 0$. Let $\tilde \ell$ denote a positive constant (that increases to $\infty$ as $n \to \infty$ in the asymptotic results).

SC(a) records a restriction on the $\ell_1$ norm of $\beta_0$ and $\rho_0$. For instance, in Examples (ref), (ref), and (ref), $D_\rho \asymp \sigma I \asymp L I$, which requires the $\ell_1$-norm of $\rho_0$ to increase at most at the speed $L\asymp \sigma$.

SC(b) is a weak assumption: the bound $\bar \lambda $ and the penalty level $ \lambda$ can be chosen proportionally to $\sqrt{ \log (p \vee n) }/\sqrt{n}$, that is $$\tilde \ell \asymp \sqrt{\log (p \vee n)}$$ using self-normalized moderate deviation bounds (jing2003self,belloni2014pivotal) or high-dimensional central limit theorems (chernozhukov2017central), under mild moment conditions, without requiring sub-Gaussianity. For instance, belloni2014pivotal employ these tools to show that, for the bounded design case $\| b\|_\infty \leq C$, $\lambda$ can be chosen as in the Gaussian error case, provided that errors follow $t(2+\delta)$ distribution (having above $2$ bounded moments), and get the error bounds similar to the Gaussian case. Here we state a general condition as our working assumption, instead of focusing on more specific condition that get us Gaussian-type conclusions.

In other words, we have instantiated $(r_1,r_2,r_3)$ for approximately sparse models in the guarantee set $\mathsf{S}$. We have the following corollary, which verifies $R(\delta)$ for approximately sparse models and hence provides sufficient conditions for Theorem (ref).

Estimation and inference using general regression learners

In this section we generalize the previous analysis to allow for any regression learner $\hat{\gamma}$ of ${\mathrm{E}}[Y|X]$ to be used in the construction of the estimator. As we have done in preceding sections we continue to include local functionals in our analysis, so that the results apply to nonregular objects as well as regular ones that can be estimated $\sqrt{n}$-consistently.

Compared to the global case, the local case may have smaller regularization and model selection biases relative to the variance. Nonetheless, bias correction is important for inference in theory and in practice. Theoretically, the local case begins to resemble the global case as the number of dimensions being integrated increases. Empirically, we provide local estimates without bias correction in Section (ref). The differences can be substantial.

The only conditions we will impose on the regression learner are certain $L^{2}$ convergence properties that we will specify in this section. These properties will allow for a wide variety of learners, including GDS, Lasso, neural nets, boosting, and others. Thus we provide estimators of local functions that can be constructed using many regression learners.

We continue to consider estimators that use cross-fitting and have the form $$ \hat{\theta}=\frac{1}{n}\sum_{k=1}^{K}\sum_{i\in I_{k}}\left\{ m(W_{i} ,\hat{\gamma}_{k})+\hat{\alpha}_{k}(X_{i})[Y_{i}-\hat{\gamma}_{k} (X_{i})]\right\} . $$ where $\hat{\gamma}_{k}$ denotes the regression learner computed from observations not in $I_{k}$ and $\hat{\alpha}_{k}(x)=b(x)^{\prime}\hat{\rho }_{k}$ is the GDS learner of the Riesz representer described in previous sections.

To allow for as many regression learners as possible under as weak conditions as possible we focus on asymptotic analysis in this section. The fundamental property we will require of $\hat{\gamma}$ is that it have some mean square convergence rate as an estimator of the true conditional mean $\gamma _{0}^{\star}.$ Specifically we require that there is $r_{1}^{\star}$ converging to zero such that for each $k,$

For purposes of formulating regularity conditions it is also useful to work with $\alpha_{0}^{\star}$ rather than $\alpha_{0}$. We will also require that

In the regular case these conditions generally follow from the mean square consistency of $\hat{\gamma}_{k}$ under boundedness of $\alpha_{0}^{\star}.$ In nonregular cases they may impose additional conditions. For example, under the conditions of Lemma (ref) it will be sufficient for these conditions to hold that $$ h^{-p_{1}/2}\left\Vert \hat{\gamma}_{k}-\gamma_{0}^{\star}\right\Vert _{P,2}\to_{{\mathrm{P}}}0. $$ This condition will hold as long as $h$ grows slowly enough relative to the mean-square convergence rate of each $\hat{\gamma}_{k}.$

Recall from Theorem (ref) that $r_{2}$ is the convergence rate of $\sigma ^{-1}\left\Vert \hat{\alpha}_{k}-\alpha_{0}\right\Vert _{P,2}$. Let $r_{2}^{\star}=r_{2}+\sigma^{-1}\left\Vert \alpha_{0}-\alpha_{0}^{\star }\right\Vert _{P,2}$ and $$ \psi_{0}^{\star}(W)=\theta_{0}-m(W,\gamma_{0}^{\star})-\alpha_{0}^{\star }(X)[Y-\gamma_{0}^{\star}(X)]. $$

This result shows that asymptotic linearity of the estimator $\hat{\theta}$ will result if $r_{2}^{\star}\rightarrow0$ fast enough relative to $r_{1}^{\star}$. Asymptotic linearity implies asymptotic Gaussian inference by standard central limit theorem arguments. As in regular doubly robust estimation problems it allows for a tradeoff between the speed of convergence $r_{2}^{\star}$ of the Riesz representer and $r_{1}^{\star}$ of the regression. It only requires a mean square convergence rate for the regression learner $\hat{\gamma}_{k}$ and so allows for a wide variety of first step machine learning estimators.

We could also formulate a non-asymptotic analog to this asymptotic result. This would depend on the availability of non-asymptotic results for the learner $\hat{\gamma}_{k}$. To the best of our knowledge such results are not available for many learners, such as neural nets and random forests. To allow the results of this section to include as many first steps as possible we focus here on the asymptotic result and reserve the non-asymptotic result to future work.

Acknowledgements

The National Science Foundation provided partial financial support via grants 1559172 and 1757140. Rahul Singh thanks the Jerry Hausman Dissertation Fellowship.