Program Evaluation and Causal Inference with High-Dimensional Data

Program Evaluation and Causal Inference with High-Dimensional Data

\enlargethispage*{\baselineskip}

\thispagestyle{empty}

Introduction

The goal of many empirical analyses is to understand the causal effect of a treatment, such as participation in a training program or a government policy, on economic and other outcomes. Such analyses are often complicated by the fact that treatments or policies are rarely randomly assigned. The lack of true random assignment has led to the adoption of a variety of quasi-experimental approaches to estimating treatment effects that are based on observational data. Such approaches include instrumental variable (IV) methods in cases where treatment is not randomly assigned but there is some other external variable, such as eligibility for receipt of a government program or service, that is either randomly assigned or the researcher is willing to take as exogenous conditional on the right set of control variables (or simply controls). Another common approach is to assume that the treatment variable itself may be taken as exogenous after conditioning on the right set of controls which leads to regression or matching based methods, among others, for estimating treatment effects.\footnote{There is a large literature about estimation of treatment effects. See, for example, the textbook treatments in AngristBook, wooldridge:text, and imbens:rubin:book.}

A practical problem empirical researchers face when trying to estimate treatment effects is deciding what conditioning variables to include. When the treatment variable or instrument is not randomly assigned, a researcher must choose what needs to be conditioned on to make the argument that the instrument or treatment is exogenous plausible. Typically, economic intuition will suggest a set of variables that might be important to control for but will not identify exactly which variables are important or the functional form with which variables should enter the model. While less crucial to identifying treatment effects, the problem of selecting controls also arises in situations where the key treatment or instrumental variables are randomly assigned. In these cases, a researcher interested in obtaining precisely estimated policy effects will also typically consider including additional controls to help absorb residual variation. As in the case where including controls is motivated by a desire to make identification of the treatment effect more plausible, one rarely knows exactly which variables will be most useful for accounting for residual variation. In either case, the lack of clear guidance about what variables to use presents the problem of selecting controls from a potentially large set including raw variables available in the data as well as interactions and other transformations of these variables.

In this paper, we consider estimation of the effect of an endogenous binary treatment, $D$, on an outcome, $Y$, in the presence of a binary instrumental variable, $Z$, in settings with very many potential controls, $f(X)$. Allowing many potential controls expressly covers both the case where there are simply many controls (where $f(X) = X$)) and the case where there are many technical controls $f(X)$ generated as transformations such as powers, b-splines, or interactions of raw controls,\footnote{See, e.g., koenker:jappliedeconometircs, newey:series, wasserman:npbook, chen:Chapter, and tsybakov:npbook.} $X$, along with combinations of the two cases. The notation $f(X)$ naturally accommodates these cases, and we call $f(X)$ the controls regardless of the case. We allow for fully heterogeneous treatment effects and thus focus on estimation of causal quantities that are appropriate in heterogeneous effects settings such as the local average treatment effect (LATE) or the local quantile treatment effect (LQTE). We focus our discussion on the endogenous case where identification is obtained through the use of an instrumental variable, but all results carry through to the exogenous case where the treatment is taken as exogenous unconditionally or after conditioning on sufficient controls by simply replacing the instrument with the treatment variable in the estimation and inference methods and in the formal results. In the latter case, LATE reduces to the average treatment effect (ATE) and LQTE to the quantile treatment effect (QTE).

The methodology for estimating treatment effects we consider allows for cases where the number of potential controls, $p := \dim f(X)$, is much larger than the sample size, $n$. Of course, informative inference about causal parameters cannot proceed allowing for $p \gg n$ without further restrictions. We impose sufficient structure through the assumption that reduced form relationships such as the conditional expectations ${\mathrm{E}}_P[D|X]$, ${\mathrm{E}}_P[Z|X]$, and ${\mathrm{E}}_P[Y|X]$ are approximately sparse. Intuitively, approximate sparsity imposes that these reduced form relationships can be represented up to a small approximation error as a linear combination, possibly inside of a known link function such as the logistic function, of a number $s \ll n$ of the variables in $f(X)$ whose identities are a priori unknown to the researcher. This assumption allows us to use methods for estimating models in high-dimensional sparse settings that are known to have good prediction properties to estimate the fundamental reduced form relationships. We may then use these estimated reduced form quantities as inputs to estimating the causal parameters of interest. Approaching the problem of estimating treatment effects within this framework allows us to accommodate the realistic scenario in which a researcher is unsure about exactly which confounding variables or transformations of these confounds are important and so must search among a broad set of controls.

Valid inference following model selection is non-trivial. Direct application of usual inference procedures following model selection does not provide valid inference about causal parameters even in low-dimensional settings, such as when there is only a single control, unless one assumes sufficient structure on the model that perfect model selection is possible. Such structure can be restrictive and seems unlikely to be satisfied in many economic applications. For example, a typical condition that allows perfect model selection is the “beta-min" condition, which requires that all but a small number of coefficients are exactly zero and that the non-zero coefficients are all large enough that they can be distinguished from zero with probability very near one in finite samples. Such a condition rules out the possibility that there may be some variables which have moderate, but non-zero, partial effects. Ignoring such variables may result in large omitted variables bias that has a substantive impact on estimation and inference regarding individual model parameters; see Leeb and P{\"o}tscher leeb:potscher:pms,leeb:potscher:review; Potscher2009; and Belloni, Chernozhukov, and Hansen BCH2011:InferenceGauss,BelloniChernozhukovHansen2011.

The first main contribution of this paper is providing inferential procedures for key parameters used in program evaluation that are theoretically valid within approximately sparse models allowing for imperfect model selection. Our procedures build upon BellChernHans:Gauss and BellChenChernHans:nonGauss, who were the first to demonstrate in a highly specialized context, that valid inference can proceed following model selection allowing for model selection mistakes under two conditions. We formulate and extend these two conditions to a rather general moment-condition framework (e.g., hansen_gmm and hansen-singleton_gmm) as follows. First, estimation should be based upon “orthogonal" moment conditions that are first-order insensitive to changes in the values of nuisance parameters that will be estimated using high-dimensional methods. Specifically, if the target parameter value $\alpha_0$ is identified via the moment condition

where $h_0$ is a function-valued nuisance parameter estimated via a model-selection or regularization method, one needs to use a moment function, $\psi$, such that the corresponding moment condition is orthogonal with respect to perturbations of $h$ around $h_0$. More formally, the moment condition should satisfy the Neyman orthogonality condition

where $\partial_h$ is a functional derivative operator with respect to $h$ restricted to directions of possible deviations of estimators of $h_0$ from $h_0$. Second, one needs to ensure that the model selection mistakes occurring in the estimation of nuisance parameters are uniformly “moderately" small with respect to the underlying model. Specifically, we will require that the nuisance parameter $h_0$ is estimated at the rate $o(n^{-1/4})$, which ensures small bias, and that the estimator takes values in a space whose entropy does not grow too fast, which ensures no overfitting. In this paper, we establish that building estimators based upon moment conditions with the orthogonality condition ((ref)) holding ensures that crude estimation of $h_0$ via post-selection or other regularization methods has an asymptotically negligible effect on the estimation of $\alpha_0$ in general frameworks. It then follows that we can form a regular, root-$n$ consistent estimator of $\alpha_0$, uniformly with respect to the underlying model.

In the endogenous treatment effects setting, we build moment conditions satisfying ((ref)) from the efficient influence functions for certain reduced form parameters, building upon hahn-pp. We illustrate how orthogonal moment conditions coupled with methods developed for forecasting in high-dimensional approximately sparse models can be used to estimate and obtain valid inferential statements about a wide variety of structural/treatment effects. We formally demonstrate the uniform validity of the resulting inference within a broad class of approximately sparse models including models where perfect model selection is theoretically impossible. An important feature of our main theoretical results is that they cover the use of variable selection for functional response data using $\ell_1$-penalized methods. Functional response data arises, for example, when one is interested in the LQTE at not just a single quantile but over a range of quantile indices. Considering this case then necessitates looking at the functional dependent variable $u \longmapsto 1(Y \leqslant u)$, where $u$ denotes various levels that $Y$ can cross. Treating such functional response data allows us to provide a unified inference procedure for interesting quantities such as the (local) distributional and quantile effects of the treatment, including simpler important parameters such as LQTE at a given quantile as a special case.

The second main contribution of this paper is providing a general set of results for uniformly valid estimation and inference methods in moment-condition problems, arising in structural analysis in econometrics and other data sciences. These results are useful not only for establishing the properties of treatment effects estimators developed here, but they are also useful for attacking a wide range of problems in structural econometrics. For example, CHS:PnP provide estimates of parameters characterizing a simple structural demand model based loosely on the analysis in BLP:Autos using the framework developed here; see also CHS:AnnRev. A key element to our establishing uniform validity of post-regularization inference is again the use of Neyman orthogonal moment conditions. In the general framework we consider, we may have (a continuum of) target parameters identified via (a continuum of) moment conditions that involve (a continuum of) nuisance functions that will be estimated via Lasso, Post-Lasso, or some other high-quality machine learning method. Our general theory expressly allows for a wide variety of traditional and machine learning methods, including those that do not rely on approximate sparsity, as long as the methods

By “not overfitting" we mean that the entropy of the function classes containing the realizations of the estimator of the nuisance function/parameter does not increase too rapidly with the sample size. This second condition can only be verified analytically, but can be avoided by the use of various data splitting methods. For example, we can set aside a vanishing fraction of the data to estimate the nuisance parameter, as in bickel:1982, or employ cross-fitting, as in Belloni et al. (2010, 2012) and CCDHM16. Either scheme ensures that there is no asymptotic efficiency loss from data-splitting. We refer the reader to CCDHM16 for a detailed discussion and analysis of cross-fitting in connection to inference on ATE and other causal parameters using machine learning methods for high-dimensional data.\footnote{Cross-fitting proceeds as follows: (1) split the sample into two equal parts, the auxiliary and main parts; (2) use the auxiliary part to estimate the nuisance parameter and the main part to estimate the target parameter, obtaining one estimator of the target parameter; (3) by reversing the roles of the main and auxiliary parts, obtain another estimator of the target parameter; and (4) average the two estimators of the target parameter to obtain the final estimator. The theorems established in Section 5 yield the properties of the final estimator; see CCDHM16.}

These results contain the results on treatment effects relevant for program evaluation, particularly the results for distributional and quantile effects, as a leading special case. These results are also immediately useful in other contexts such as nonseparable quantile models as in iqr:ema, CH06, C03, and IN09; semiparametric and partially identified models as in EZ2013; and many others. In our results, we first establish a functional central limit theorem for the continuum of target parameters and show that this functional central limit theorem holds uniformly in a wide range of data-generating processes $P$ with approximately sparse continua of nuisance functions. Second, we establish a functional central limit theorem for the multiplier bootstrap that resamples the first order approximations to the standardized estimators and demonstrate its uniform-in-$P$ validity. These uniformity results build upon and complement those given in Romano:Shaikh:AoS for the empirical bootstrap. Third, we establish a functional delta method for smooth functionals of the continuum of target parameters and a functional delta method for the multiplier bootstrap of these smooth functionals, both of which hold uniformly in $P$, using an appropriately strengthened notion of Hadamard differentiability. All of these results are new and are of independent interest outside of the treatment effects focus of this paper.

We illustrate the use of our methods by estimating the effect of 401(k) eligibility and 401(k) participation on measures of accumulated assets as in CH401k.\footnote{See also Poterba, Venti, and Wise pvw:94,pvw:95,pvw:nber96,pvw:01, abadie:401k, benjamin, and ORR:401k among others.} Similar to CH401k, we provide estimates of ATE and QTE of 401(k) eligibility and of LATE and LQTE of 401(k) participation. We differ from this previous work by using the high-dimensional methods developed in this paper to allow ourselves to consider a broader set of controls than has previously been considered. We find that 401(k) participation has a moderate impact on accumulated financial assets at low quantiles while appearing to have a much larger impact at high quantiles. Interpreting the quantile index as “preference for savings” as in CH401k, this pattern suggests that 401(k) participation has little causal impact on the accumulated financial assets of those with low desire to save but a much larger impact on those with stronger preferences for saving. It is interesting that these results are similar to those in CH401k despite allowing for a much richer set of controls.

Links to the literature

The Neyman orthogonality condition embodied in ((ref)) has a long history in statistics and econometrics. For example, this type of orthogonality was used by Neyman1979 in low-dimensional settings to deal with crudely estimated parametric nuisance parameters. See also newey90, andrews94, newey94, robins:dr, and linton96 for the use of this condition in semi-parametric problems.

To the best of our knowledge, BellChernHans:Gauss and BellChenChernHans:nonGauss were the first to use the orthogonality ((ref)) to expressly address the question of the uniform post-selection inference without imposing “beta-min" conditions, either in high-dimensional settings with $p \gg n$ or in low-dimensional settings with $p \ll n$. They applied it to the specific problem of the linear instrumental variables model with many instruments where the nuisance function $h_0$ is the optimal instrument estimated by Lasso or Post-Lasso methods and $\alpha_0$ is the coefficient of the endogenous regressor. BCH2011:InferenceGauss and BelloniChernozhukovHansen2011 also exploited this approach to develop a double-selection method that yields valid post-selection inference on the parameters of the linear part of a partially linear model and on average treatment effects when the treatment is binary and exogenous conditional on controls in both the $p \gg n$ and the $p \ll n$ setting.\footnote{Note that these results as well as results of this paper on the uniform post-selection inference in moment-condition problems are new for either $p\ll n$ or $p \gg n$ settings. The results also apply to arbitrary model selection devices, such as the Dantzig selector, Square-Root-Lasso, or Adaptive Lasso, that are able to select good sparse approximating models; and “moderate” model selection errors are explicitly allowed in the paper.} Subsequently, Farrell:JMP extended the results of BCH2011:InferenceGauss and BelloniChernozhukovHansen2011 to estimation of ATE when the treatment is multivalued and exogenous conditional on controls using group penalization for selection. Note that this previous work on treatment effects covers only the exogenous case and does not allow for functional responses which are necessary, for example, for working with distributional or quantile treatment effects.

Our work also contributes to the line of research on obtaining $\sqrt{n}$-consistent and asymptotically normal estimates for low-dimensional components within traditional semiparametric frameworks as in the important work by bickel:1982, robinson, newey90, vaart:1991, andrews94, newey94, ai:chen, AC2012, and CLK:EfficientSP. The major difference is that we allow for the use of modern high-dimensional methods, a.k.a. machine learning methods, for modeling and fitting the non-parametric (or high-dimensional) components of the model. In contrast to the former literature, we expressly allow for data-driven choice of the approximating model for the high-dimensional component, which addresses a crucial problem that arises in empirical work. Moreover, recent methods based on $\ell_1$-penalization, upon which we focus in this paper, allow for much more flexible modeling of the non-parametric/high-dimensional parts of the model.\footnote{See, for instance, BellChenChernHans:nonGauss and BCW-SqLASSO2 for a formalization of this claim in terms of rearranged Sobolev spaces where it is shown that traditional methods can fail to be consistent while $\ell_1$-penalized methods remain consistent and have good rates of convergence.} Our general theory in Section 5 also allows, in principle, for a wide variety of both traditional and machine learning methods.

The paper also generates a number of new results on sparse estimation with functional response data. These results are of independent interest in themselves, and they build upon the work of BC-SparseQR who provided rates of convergence for variable selection when one is interested in estimating the quantile regression process with exogenous variables. More generally, this theoretical work complements and extends the rapidly growing set of results for $\ell_1$-penalized estimation methods; see, for example, FF:1993; T1996; FanLi2001; Zou2006; CandesTao2007; vdGeer; HHS2008; BickelRitovTsybakov2009; MY2007; Bach2010; horowitz:lasso; BC-SparseQR; kato; BellChenChernHans:nonGauss; BC-PostLASSO; BCK-LAD; BCY-honest; CanerZhang:GMMEL; and the references therein.

Plan of the Paper

Section (ref) introduces the structural parameters for policy evaluation and relates these parameters to reduced form functions. Section (ref) describes a three step procedure to estimate and make inference on the structural parameters and functionals of these parameters, and Section (ref) provides asymptotic theory in the treatment effects setting. Section (ref) generalizes the setting and results to moment-condition problems with a continuum of structural parameters and a continuum of reduced form functions. Section (ref) derives general asymptotic theory for the Lasso and post-Lasso estimators for functional response data used in the estimation of the reduced form functions. Section (ref) presents the empirical application. We provide notation, proofs of key results, and details about implementation of the methods in the empirical example in Appendices (ref)--(ref). An on-line Supplementary Appendix provides all remaining proofs, additional technical material, and results from a small Monte Carlo simulation bcfh15sup.

The Treatment Effects Setting and Target Parameters

Observables and Reduced Form Parameters

The observed random variables consist of $( (Y_u)_{u \in \mathcal{U}}, X,Z,D)$. The outcome variable of interest $Y_u$ is indexed by $u \in \mathcal{U}$. We give examples of the index $u$ below. The variable $D \in \mathcal{D}=\{0,1\}$ is a binary indicator of the receipt of a treatment or participation in a program. It will typically be treated as endogenous; that is, we will typically view the treatment as assigned non-randomly with respect to the outcome. The instrumental variable $Z \in \mathcal{Z}=\{0,1\}$ is a binary indicator, such as an offer of participation, that is assumed to be randomly assigned conditional on the observable covariates $X$ with support $\mathcal{X}$.\footnote{Of course, by “randomly assigned" we mean independently of potential outcomes conditional on the covariates.} For example, we argue that 401(k) eligibility can be considered exogenous only after conditioning on income and other individual characteristics in the empirical application. The notions of exogeneity and endogeneity we employ are standard and thus omitted.\footnote{For completeness, we provide a review of these conditions as well as restate standard conditions that are sufficient for a causal interpretation of the target parameters in the Supplementary Appendix.}

The indexing of the outcome $Y_u$ by $u$ is useful to analyze functional data. For example, $Y_u$ could represent an outcome falling short of a threshold, namely $Y_u = 1(Y \leqslant u)$, in the context of distributional analysis; $Y_u $ could be a height indexed by age $u$ in growth charts analysis; or $Y_u$ could be a health outcome indexed by a dosage $u$ in dosage response studies. Our framework is tailored for such functional response data. The special case with no index is included by simply considering $\mathcal{U}$ to be a singleton set.

We make use of two key types of reduced form parameters for estimating the structural parameters of interest -- (local) treatment effects and related quantities. These reduced form parameters are defined as

where $z = 0$ or $z =1$ are the fixed values of $Z$.\footnote{The expectation that defines $\alpha_V(z)$ is well-defined under the standard support condition $0 < c < {\mathrm{P}}_P(Z=1 \mid X) < 1 - c < 1$ a.s. This condition is standard in treatment effects estimation; see, e.g., the supplementary appendix. We impose this condition in Assumption (ref).} The function $ g_{V}$ maps $\mathcal{ZX}$, the support of the vector $(Z,X)$, to the real line $\mathbb{R}$ and is defined as

We use $V$ to denote a target variable whose identity may change depending on the context such as $V=\mathbf{1}_d(D) Y_u$ or $V=\mathbf{1}_d(D)$ where $\mathbf{1}_d(D) := 1 (D=d)$ is the indicator function.

All the structural parameters we consider are smooth functionals of these reduced-form parameters. In our approach to estimating treatment effects, we estimate the key reduced form parameter $\alpha_V(z)$ using modern methods to deal with high-dimensional data coupled with orthogonal estimating equations. The orthogonality property allows us to deal with the “non-regular" nature of penalized and post-selection estimators which do not admit linearizations except under very restrictive conditions. The use of regularization by model selection or penalization is in turn motivated by the desire to accommodate high-dimensional data.

Target Structural Parameters -- Local Treatment Effects

The reduced form parameters defined in ((ref)) are key because the structural parameters of interest are functionals of these elementary objects. The local average structural function (LASF) defined as

underlies the formation of many commonly used treatment effects. Under standard assumptions, the LASF identifies average potential outcomes for the group of compliers, individuals whose treatment status may be influenced by variation in the instrument, in the treated and non-treated states; see, e.g. Abadie abadie:bstest,abadie:401k. The local average treatment effect (LATE) of imbens:angrist:94 corresponds to the difference of the two values of the LASF:

The term local designates that this parameter does not measure the effect on the entire population but rather measures the effect on the subpopulation of compliers.\footnote{The methods of the paper can be extended to analyze the marginal treatment effects of heckman:vytlacil, hv2005.}

When there is no endogeneity, formally when $D \equiv Z$, the LASF and LATE become the average structural function (ASF) and average treatment effect (ATE) on the entire population. Thus, our results cover this situation as a special case where the ASF and ATE simplify to

We also note that the impact of the instrument $Z$ itself may be of interest since $Z$ often encodes an offer of participation in a program. In this case, the parameters of interest are again simply the reduced form parameters $$ \alpha_{Y_u}(z) , \ \ \alpha_{Y_u}(1) - \alpha_{Y_u}(0).$$ Thus, the LASF and LATE are primary targets of interest in this paper, and the ASF and ATE are subsumed as special cases.

Local Distribution and Quantile Treatment Effects

Setting $Y_u = Y$ in ((ref)) and ((ref)) provides the conventional LASF and LATE. An important generalization arises by letting $Y_u = 1(Y \leqslant u)$ be the indicator of the outcome of interest falling below a threshold $u \in \mathbb{R}$. In this case, the family of effects

describe the local distribution treatment effects (LDTE). Similarly, we can look at the quantile left-inverse transform of the curve $u \longmapsto \theta_{Y_u}(d)$,

and examine the family of local quantile treatment effects (LQTE):

The LQTE identify the differences of quantiles between the distribution of potential outcomes in the treated and non-treated states for compliers.

Target Structural Parameters -- Local Treatment Effects on the Treated

We may also be interested in local treatment effects on the treated. The key object in defining these effects is the local average structural function on the treated (LASF-T) which is defined by its two values:

The LASF-T identifies average potential outcomes for the group of treated compliers in the treated and non-treated states under standard assumptions. The local average treatment effect on the treated (LATE-T) introduced in hong:nekipelov:2010 and frolich:melly is the difference of two values of the LASF-T:

The LATE-T may be of interest because it measures the average treatment effect for treated compliers, namely the subgroup of compliers that actually receive the treatment.

When the treatment is assigned randomly given controls so we can take $D=Z$, the LASF-T and LATE-T become the average structural function on the treated (ASF-T) and average treatment effect on the treated (ATE-T). In this special case, the ASF-T and ATE-T simplify to

and we can use our results to provide estimation and inference methods for these quantities.

Local Distribution and Quantile Treatment Effects on the Treated

Local distribution treatment effects on the treated (LDTE-T) and local quantile treatment effects on the treated (LQTE-T) can also be defined. As in Section 2.2.1, we let $Y_u = 1(Y \leqslant u)$ be the indicator of the outcome of interest falling below a threshold $u$. The family of treatment effects

then describes the LDTE-T. We can also use the quantile left-inverse transform of the curve $u \longmapsto \vartheta_{Y_u}(d)$, namely $ \vartheta^{\leftarrow}_Y(\tau, d) := \inf\{ u \in \mathbb{R} : \vartheta_{Y_u}(d) \geqslant \tau \},$ and define the LQTE-T:

Under conditional exogeneity LQTE and LQTE-T reduce to the quantile treatment effects (QTE) and quantile treatment effects on the treated (QTE-T) koenker:book.

Estimation of Reduced-Form and Structural Parameters in a Data-Rich Environment

The key objects used to define the structural parameters in Section 2 are the expectations

where $g_V(z,X)= {\mathrm{E}}_P[V|Z=z,X]$ and $V$ denotes a variable whose identity will change with the context. Specifically, we shall vary $V$ over the set $\mathcal{V}_u$:

It is clear that $g_V(z,X)$ will play an important role in estimating $\alpha_V(z)$. A related function that will also play an important role in forming a robust estimation strategy is the propensity score $m_Z: \mathcal{ZX} \longmapsto \mathbb{R}$ defined by

We will denote other potential values for the functions $g_{V}$ and $m_Z$ by the parameters $g$ and $m$, respectively. We can then estimate $\alpha_V(z)$ by estimating $g_V$ and $m_Z$ using high-dimensional modeling and estimation methods.\footnote{Note that there is an alternative approach based on decomposing $g_V$ as $g_{V}(z,x) = \sum_{d=0}^1 e_{V}(d,z,x) l_{D}(d,z,x)$ where the regression functions $e_{V}$ and $l_{D}$ map the support of $(D,Z,X)$, $\mathcal{DZX}$, to the real line and are defined by $e_{V}(d,z,x): = {\mathrm{E}}_P[V|D=d, Z=z,X=x]$ and $l_{D}(d,z,x): = {\mathrm{P}}_P[D=d|Z=z, X=x]$. We provide some discussion of this approach in the supplementary appendix.}

In the rest of this section, we describe the estimation of the reduced-form and structural parameters. The estimation method consists of 3 steps:

1) Estimate the predictive relationships $m_Z$ and $g_V$ using high-dimensional nonparametric methods with model selection.

2) Estimate the reduced form parameters $\alpha_V$ and $\gamma_V$ using orthogonal estimating equations to immunize the reduced form estimators to imperfect model selection in the first step.

3) Estimate the structural parameters and effects via the plug-in rule.

First Step: Modeling and Estimating $g_{V}$ and $m_{Z}$

In this section, we discuss estimation of the conditional expectation functions $g_{V}$ and $m_{Z}$. Since these functions are unknown and potentially complicated, we use a generalized linear combination of a large number of control terms

to approximate $g_{V}$ and $m_{Z}$. Specifically, we use

In these equations, $r_{V}(z,x)$ and $r_{Z}(x)$ are approximation errors, and the functions $ \Lambda_V (f(z,x) '\beta_{V})$ and $\Lambda_Z(f(x)'\beta_{Z})$ are generalized linear approximations to the target functions $g_{V}(z,x)$ and $m_{Z}(1,x)$. The functions $\Lambda_V$ and $\Lambda_Z$ are taken to be known link functions $\Lambda$. The most common example is the linear link $\Lambda(u) = u$. When the response variable is binary, we may also use the logistic link $\Lambda(u) = \Lambda_0 (u) = e^u/(1+ e^u)$ and its complement $1- \Lambda_0(u)$ or the probit link $\Lambda(u) = \Phi(u) = (2\pi)^{-1/2}\int_{-\infty}^u e^{-z^2/2}dz$ and its complement $1-\Phi(u)$. For clarity, we use links from the finite set $\mathcal{L}=\{ \mathrm{Id}, \Phi, 1-\Phi, \Lambda_0, 1-\Lambda_0\}$ where $\mathrm{Id}$ is the identity (linear) link.

As discussed in the Introduction, the dictionary of controls, denoted by $ f(X)$, can be “rich" in the sense that its dimension $p=p_n$ may be large relative to the sample size. Specifically, our results require only that $\log p = o(n^{1/3})$ along with other technical conditions. We also note that the functions $f$ forming the dictionary can depend on $n$, but we suppress this dependence.

Having very many controls $f(X)$ creates a challenge for estimation and inference. A useful condition that makes it possible to perform constructive estimation and inference in such cases is termed approximate sparsity or simply sparsity. Sparsity imposes that there exist approximations of the form given in ((ref))-((ref)) that require only a small number of non-zero coefficients to render the approximation errors small relative to estimation error. More formally, sparsity relies on two conditions. First, there must exist $\beta_{V}$ and $\beta_{Z}$ such that, for all $V \in \mathcal{V} := \{\mathcal{V}_u : u \in \mathcal{U}\},$

where $\|x \|_0$ is the number of non-zero components of vector $x$ and all other norms we use are defined in Appendix A. That is, there are at most $s=s_n \ll n$ components of $f(Z,X)$ and $f(X)$ with nonzero coefficient in the approximations to $g_V$ and $m_Z$. Second, the sparsity condition requires that the size of the resulting approximation errors is small compared to the conjectured size of the estimation error; namely, for all $V \in \mathcal{V}$,

Note that the size of the approximating model $s=s_n$ can grow with $n$ just as in standard series estimation, subject to the rate condition $$s^2 \log^2 (p\vee n) \log^2 n/n \to 0.$$ These conditions ensure that the functions $g_V$ and $m_Z$ are estimable at a $o(n^{-1/4})$ rate and are used to derive asymptotic normality results for the structural and reduced-form parameter estimators. They could be relaxed through the use of sample splitting methods as in BellChenChernHans:nonGauss.

The high-dimensional-sparse-model framework outlined above extends the standard framework in the program evaluation literature which assumes both that the identities of the relevant controls are known and that the number of such controls $s$ is small relative to the sample size.\footnote{For example, one would select a set of basis functions, $\{f_{j}(X)\}_{j=1}^{\infty}$, such as power series or splines and then use only the first $s \ll n$ terms in the basis under the assumption that $s^C/n \rightarrow 0$ for some number $C$ whose value depends on the specific context in a standard nonparametric approach using series.} Instead, we assume that there are many, $p$, potential controls of which at most $s$ controls suffice to achieve a desirable approximation to the unknown functions $g_V$ and $m_{Z}$; and we allow the identity and number of these controls to be unknown. Relying on this assumed sparsity, we use selection methods to choose approximately the right set of controls.

Current estimation methods that exploit approximate sparsity employ different types of regularization aimed at producing estimators that theoretically perform well in high-dimensional settings while remaining computationally tractable. Many widely used methods are based on $\ell_1$-penalization. The Lasso method is one such commonly used approach that adds a penalty for the weighted sum of the absolute values of the model parameters to the usual objective function of an M-estimator. A related approach is the Post-Lasso method which performs re-estimation of the model after selection of variables by Lasso. These methods are discussed at length in recent papers and review articles; see, for example, BCH2011:InferenceGauss.

In the following, we outline the general features of the Lasso and Post-Lasso methods focusing on estimation of $g_V$. Given the data $(\tilde Y_i, \tilde X_i )_{i=1}^n = (V_i, f(Z_i, X_i))_{i=1}^n$, the Lasso estimator $\widehat \beta_{V}$ solves

where $\widehat \Psi = {\rm diag}(\widehat l_1,\ldots,\widehat l_{\dim(\widetilde X)})$ is a diagonal matrix of data-dependent penalty loadings, $M(y, t) = (y-t)^2/2$ in the case of linear regression, and $M(y, t) = -\{1(y=1) \log \Lambda_V(t)+ 1(y=0) \log(1- \Lambda_V(t))\}$ in the case of binary regression. The penalty level, $\lambda$, and loadings, $\widehat l_j, \ j = 1,...,\dim(\widetilde X)$, are selected to guarantee good theoretical properties of the method. We provide further discussion of these methods for estimation of a continuum of functions in Section (ref), and we specify detailed implementation algorithms used in the empirical example in Appendix (ref). A key consideration in this paper is that the penalty level needs to be set to account for the fact that we will be simultaneously estimating potentially a continuum of Lasso regressions since our $V$ varies over the list $\mathcal{V}_u$ with $u$ varying over the index set $\mathcal{U}$.

The Post-Lasso method uses $\widehat \beta_V$ solely as a model selection device. Specifically, it makes use of the labels of the regressors with non-zero estimated coefficients, $ \widehat I_V := \mathrm{supp}(\widehat \beta_V). $ The Post-Lasso estimator is then a solution to

A main contribution of this paper is establishing that the estimator $\widehat g_{V}(Z,X) = \Lambda_V(f(Z,X)'\bar \beta_V)$ of the regression function $g_{V}(Z,X)$, where $\bar \beta_V= \widehat \beta_V$ or $\bar \beta_V= \tilde \beta_V$, achieves the near oracle rate of convergence $\sqrt{(s \log p)/n} $ and maintains desirable theoretic properties while allowing for a continuum of response variables.

Estimation of $m_Z$ proceeds similarly. The Lasso estimator $\widehat \beta_Z$ and Post-Lasso estimator $\tilde \beta_Z$ are defined analogously to $\widehat \beta_V$ and $\tilde \beta_V$ using the data $(\tilde Y_i, \tilde X_i)_{i=1}^n$= $( Z_i, f(X_i))_{i=1}^n$. The estimator $\widehat m_Z(1, X) = \Lambda_Z(f(X)'\bar \beta_Z)$ of $m_{Z}(X)$, with $\bar \beta_Z= \widehat \beta_Z$ or $\bar \beta_Z= \tilde \beta_Z$, also achieves the near oracle rate of convergence $\sqrt{(s \log p)/n}$ and has other good theoretic properties. The estimator of $\widehat m_Z(0,X)$ is then formed as $1 - \widehat m_Z(1,X)$.

Second Step: Robust Estimation of the Reduced-Form Parameters $\alpha_V(z)$ and $\gamma_V$

Estimation of the key quantities $\alpha_V(z)$ will make heavy use of orthogonal moment functions as defined in ((ref)). These moment functions are closely tied to efficient influence functions, where efficiency is in the sense of locally minimax semi-parametric efficiency. The use of these functions will deliver robustness with respect to the non-regularity of the post-selection and penalized estimators needed to manage high-dimensional data. The use of these functions also automatically delivers semi-parametric efficiency for estimating and performing inference on the reduced-form parameters and their smooth transformations -- the structural parameters.

The efficient influence function and orthogonal moment function for $\alpha_V(z)$, $z \in \mathcal{Z} = \{0,1\}$, are given respectively by

This efficient influence function was derived by hahn-pp; it has recently been used by cattaneo2010efficient in the series context (with $p \ll n$) and rothe_firpo2013 in the kernel context. The efficient influence function and the moment function for $\gamma_V$ are trivially given by

We then define estimators of the reduced-form parameters $\alpha_V(z)$ and $\gamma_V(z)$ as solutions $\alpha = \widehat \alpha_V(z)$ and $\gamma= \widehat \gamma_V$ to the equations

where $\widehat g_V$ and $\widehat m_Z$ are constructed as in Section (ref). We apply this procedure to each variable name $V \in \mathcal{V}_u$ and obtain the estimator\footnote{By default notation, $(a_j)_{j \in \mathcal{J}}$ returns a column vector produced by stacking components together in some consistent order. }

The estimator and the parameter are vectors in $\mathbb{R}^{d_\rho}$ with dimension $d_{\rho} = 3 \times \dim \mathcal{V}_u = 15$.

In the next section, we formally establish a principal result which shows that

where $\mathcal{P}_n$ is a rich set of data generating processes $P$ which includes cases where perfect model selection is impossible theoretically. The notation “$Z_{n,P} \rightsquigarrow Z_P$ uniformly in $P \in \mathcal{P}_n$" is defined formally in Appendix (ref) and can be read as “$Z_{n,P}$ is approximately distributed as $Z_P$ uniformly in $P \in \mathcal{P}_n$." This usage corresponds to the usual notion of asymptotic distribution extended to handle uniformity in $P$.

We then stack all the reduced form estimators and parameters over $u \in \mathcal{U}$ as $$ \widehat \rho = (\widehat \rho_u)_{u \in \mathcal{U}} \ \text{ and } \ \ \rho = (\rho_u)_{u \in \mathcal{U}}, $$ giving rise to the empirical reduced-form process $\widehat \rho$ and the reduced-form function-valued parameter $\rho$. We establish that $\sqrt{n}(\widehat \rho - \rho)$ is asymptotically Gaussian: In $\ell^\infty(\mathcal{U})^{d_\rho}$,

where $\mathbb{G}_{P}$ denotes the $P$-Brownian bridge vdV-W. This result contains ((ref)) as a special case and again allows $\mathcal{P}_n$ to be a “rich" set of data generating processes $P$ that includes cases where perfect model selection is impossible theoretically. Importantly, this result verifies that the functional central limit theorem applies to the reduced-form estimators in the presence of possible model selection mistakes.

Since some of our objects of interest are complicated, inference can be facilitated by a multiplier bootstrap method as in GineZinn1984. We define $\widehat \rho^* = (\widehat \rho_u^*)_{u \in \mathcal{U}}$, a bootstrap draw of $\widehat \rho$, via

Here $(\xi_i)_{i=1}^n$ are i.i.d. copies of $\xi$ which are independently distributed from the data $(W_i)_{i=1}^n$ and whose distribution $P_\xi$ does not depend on $P$. We also impose that

Examples of $\xi$ include (a) $\xi = \mathcal{E}-1$, where $\mathcal{E}$ is a standard exponential random variable, (b) $\xi = \mathcal{N}$, where $\mathcal{N}$ is a standard normal random variable, and (c) $\xi = \mathcal{N}_1/\sqrt{2} + (\mathcal{N}_2^2-1)/2$, where $\mathcal{N}_1$ and $\mathcal{N}_2$ are mutually independent standard normal random variables.\footnote{We do not consider the nonparametric bootstrap, which corresponds to using multinomial multipliers $\xi$, to reduce the length of the paper; but we note that the conditions and analysis could be extended to cover this case.} The choices of (a), (b), and (c) correspond respectively to the Bayesian bootstrap (e.g., Hahn97 and chamberlain2003nonparametric), the Gaussian multiplier method (e.g, GineZinn1984 and vdV-W), and the wild bootstrap method (mammen1993:bootstrap).\footnote{ The motivation for method (c) is that it is able to match 3 moments since ${\mathrm{E}}[\xi^2] = {\mathrm{E}}[\xi^3]=1$. Methods (a) and (b) do not satisfy this property since ${\mathrm{E}}[\xi^2] = 1$ but ${\mathrm{E}}[\xi^3]\neq 1$ for these approaches. } $\widehat \psi^\rho_{u}$ in ((ref)) is an estimator of the influence function $\psi^{\rho}_u$ defined via the plug-in rule:

Note that this bootstrap is computationally efficient since it does not involve recomputing the influence functions $\widehat \psi_u^{\rho}$.\footnote{CH06 and HS2006 proposed related computationally efficient bootstrap schemes that resample the influence functions.} Each new draw of $(\xi_i)_{i=1}^n$ generates a new draw of $\widehat \rho^*$ holding the data and the estimates of the influence functions fixed. This method simply amounts to resampling the first-order approximations to the estimators. Here we build upon prior uses of this or similar methods in low-dimensional settings such as Hansen96 and KS12.

We establish that the bootstrap law of $\sqrt{n} (\widehat \rho^*- \widehat \rho)$ is uniformly asymptotically consistent: In the metric space $\ell^\infty(\mathcal{U})^{d_\rho}$, conditionally on the data,

where $\rightsquigarrow_B$ denotes weak convergence of the bootstrap law in probability, as defined in Appendix (ref).

Third Step: Robust Estimation of the Structural Parameters

All structural parameters we consider take the form of smooth transformations of the reduced-form parameters:

The structural parameters may themselves carry an index $q \in \mathcal{Q}$ that can be different from $u$; for example, the LQTE is indexed by a quantile index $q \in (0,1)$. This formulation includes as special cases all the structural functions of Section (ref). We estimate these quantities by the plug-in rule. We establish the asymptotic behavior of these estimators and the validity of the bootstrap as a corollary from the results outlined in Section 3.2 and the functional delta method (extended to handle uniformity in $P$).

For the application of the functional delta method, we require that the functional $\rho \longmapsto \phi(\rho)$ be Hadamard differentiable uniformly in $\rho \in \mathbb{D}_{\rho}$, where $\mathbb{D}_{\rho}$ is a set that contains the true values $\rho= \rho_P$ for all $P \in \mathcal{P}_n$, tangentially to a subset that contains the realizations of $Z_P$ for all $P \in \mathcal{P}_n$ with derivative map $h \longmapsto \phi'_\rho(h) = (\phi'_{\rho}(h)(q))_{q \in \mathcal{Q}}$.\footnote{We give the definition of uniform Hadamard differentiability in Definition (ref) of Appendix (ref).} We define the estimators of the structural parameters and their bootstrap versions via the plug-in rule as

We establish that these estimators are asymptotically Gaussian

and that the bootstrap consistently estimates their large sample distribution:

These results can be used to construct simultaneous confidence bands and test functional hypotheses on $\Delta$ using the methods described for example in CF15 and CFM.

Theory: Estimation and Inference on Local Treatment Effects Functionals

Consider fixed sequences of numbers $\delta_n \searrow 0$, $\epsilon_n \searrow 0$, $\Delta_n \searrow 0$, at a speed at most polynomial in $n$ (for example, $\delta_n \geqslant 1/n^c$ for some $c > 0$), $\ell_n := \log n$, and positive constants $c$, $C$, and $c' <1/2$. These sequences and constants will not vary with $P$. The probability $P$ can vary in the set $\mathcal{P}_n$ of probability measures, termed “data-generating processes", where $\mathcal{P}_n$ is typically a set that is weakly increasing in $n$, i.e. $\mathcal{P}_n \subseteq \mathcal{P}_{n+1}$. Other definitions and notation are collected in Appendix A.

Assumption (ref) is stated to deal with the measurability issues associated with functional response data. This assumption also implies that the set of functions $(\psi_u^{\rho})_{ u \in \mathcal{U}}$, where $$\psi^\rho_{u}:= (\{\psi^\alpha_{V,0}, \psi^\alpha_{V,1}, \psi^\gamma_{V}\})_{ V \in \mathcal{V}_u},$$ is $P$-Donsker uniformly in $\mathcal{P}$. That is, it implies

with $\mathbb{G}_{P}$ denoting the $P$-Brownian bridge vdV-W and with $Z_P$ having bounded, uniformly continuous paths uniformly in $P \in \mathcal{P}$:

We work with the sequence of constants defined prior to Assumption (ref).

Under the stated assumptions, the empirical reduced form process $ \widehat Z_{n,P} = \sqrt{n}(\widehat \rho - \rho)$ defined by ((ref)) obeys the following relations. We recall definitions of convergence uniformly in $P \in \mathcal{P}_n$ in Appendix (ref).

Another main result of this section shows that the bootstrap law of the process $$ \widehat Z^*_{n,P} := \sqrt{n} (\widehat \rho^* - \widehat \rho) := \frac{1}{\sqrt{n}} \sum_{i=1}^n \xi_i \widehat \psi^{\rho}_u(W_i), $$ where $\widehat \psi^\rho_u$ is defined in ((ref)), provides a valid approximation to the large sample law of $\sqrt{n}(\widehat \rho - \rho)$.

Next we consider inference on the structural functionals $\Delta$ defined in ((ref)). We derive the large sample distribution of the estimator $\widehat \Delta$ in ((ref)), and show that the multiplier bootstrap law of $\widehat \Delta^*$ in ((ref)) provides a consistent approximation to that distribution. We rely on the functional delta method in our derivations, which we modify to handle uniformity with respect to the underlying d.g.p. $P$. Our argument relies on the following assumption on the structural functionals.

The definition of uniform Hadamard differentiability is given in Definition (ref) of Appendix (ref). Assumption (ref) holds for all examples of structural parameters listed in Section 2. \\

The following corollary gives the large sample law of $\sqrt{n} (\widehat \Delta - \Delta)$, the properly normalized structural estimator. It also shows that the bootstrap law of $ \sqrt{n} (\widehat \Delta^* - \widehat \Delta), $ computed conditionally on the data, approaches the large sample law $\sqrt{n} (\widehat \Delta - \Delta)$. It follows from the previous theorems as well as from a more general result contained in Theorem (ref).

General Theory: Honest Inference in General Moment Condition Problems with Nuisance Functions Estimated by Machine Learning Methods

In this section, we consider a general moment condition framework, where possibly a continuum of target parameters is of interest and we use modern machine learning methods, with Lasso-type methods being a lead example, to estimate a continuum of high-dimensional nuisance functions. This setting covers a rich variety of modern moment-condition problems in econometrics including the treatment effects problem. We establish a functional central limit theorem for the estimators of the continuum of target parameters that holds uniformly in $P \in \mathcal{P}$, where $\mathcal{P}$ includes a wide range of data-generating processes with well-approximable continuums of nuisance functions. We also derive a functional central limit theorem for the multiplier bootstrap that resamples the first order approximations to the standardized estimators of the continuum of target parameters and establish its uniform validity. Moreover, we establish the uniform validity of the functional delta method and the functional delta method for the multiplier bootstrap for smooth functionals of the continuum of target parameters using an appropriate strengthening of Hadamard differentiability.

Setting

We are interested in function-valued target parameters indexed by $u \in \mathcal{U} \subset \mathbb{R}^{d_u}$. We denote the true value of the target parameter by $$\theta^0 = (\theta_u)_{u \in \mathcal {U}}, \text{ where } \theta_u \in \Theta_u \subset \Theta \subset \mathbb{R}^{d_\theta}, \text{ for each } u \in \mathcal{U}. $$ We assume that for each $u \in \mathcal{U},$ the true value $\theta_u$ is identified as the solution to the following moment condition:

where $W_u$ is a random vector that takes values in a Borel set $\mathcal{W}_u \subset \mathbb{R}^{d_w}$ and contains as a subcomponent the vector $Z_u$ taking values in a Borel set $\mathcal{Z}_u$, the moment function

is a Borel measurable map, and the function

is another Borel measurable map that denotes the possibly infinite-dimensional nuisance parameter. The sets $T_u(z)$ are assumed to be convex for each $u \in \mathcal{U}$ and $z \in \mathcal{Z}_u$. Finite-dimensional nuisance parameters that do not depend on $Z_u$ are treated as part of $h_u$ as well.

We assume that the continuum of nuisance functions $(h_{u})_{u \in \mathcal{U}}$ is well-approximable and can be well estimated by the modern generation of statistical and machine learning methods. In particular, our regularity conditions allow for approximately sparse nuisance functions, which can be modeled and estimated using methods such as Lasso and Post-Lasso. We let $\widehat h_u= (\widehat h_{um})_{m=1}^{d_t}$ denote the estimator of $h_u$, which we assume obeys the conditions in Assumption (ref). The estimator $\widehat \theta_u$ of $\theta_u$ is constructed as any approximate $\epsilon_n$-solution in $\Theta_u$ to a sample analog of the moment condition ((ref)), i.e.,

The Neyman Orthogonality or Immunization Condition

A key condition needed for regular estimation of $\theta_u$ is an orthogonality or immunization condition. The simplest to explain, yet strongest, form of this condition can be expressed as follows:

subject to additional technical conditions such as continuity ((ref)) and dominance ((ref)) stated below, where we use the symbol $\partial_t$ to abbreviate $\frac{\partial}{\partial t'}$. This condition holds in the previous setting of inference on treatment effects after interchanging the order of the derivative and expectation. The formulation here also covers certain non-smooth cases such as structural and instrumental quantile regression problems.

In the formal development, we use a more general form of the orthogonality condition.

The orthogonality condition ((ref)) reduces to ((ref)) when $\mathcal{H}_u$ can span all measurable functions $h: \mathcal{Z}_u \longmapsto T_u$ such that $\|h\|_{P,2} < \infty$ but is more general otherwise.

Regularity Conditions and Results

In what follows, we shall denote by $\delta$, $c_0$, $c$, and $C$ some positive constants. For a positive integer $d$, $[d]$ denotes the set $\{1,\ldots, d\}.$ We shall impose the following regularity conditions.

The conditions of Assumption (ref) are mild and standard in moment condition problems. Assumption (ref)(iv) encodes sufficient global and local identifiability to obtain a rate result. The suitably measurable condition, defined in Appendix (ref), is a mild condition satisfied in most practical cases.

Assumption (ref) imposes smoothness and integrability conditions on various quantities derived from $\psi_u$. It also imposes conditions on the complexity of the relevant function classes.

In what follows, let $\Delta_n \searrow 0$, $\delta_n \searrow 0$, and $\tau_n \searrow 0$ be sequences of constants approaching zero from above at a speed at most polynomial in $n$ (for example, $\delta_n \geqslant 1/n^c$ for some $c > 0$). \\

The following theorem is one of the main results of the paper:

We can estimate the law of $Z_P$ with the bootstrap law of

where $(\xi_i)_{i=1}^n$ are i.i.d. multipliers as defined in equation ((ref)), $ \widehat \psi_u(W_i)$ is the estimated score $$ \widehat \psi_u(W_i):= - \widehat J_u^{-1} \psi_u(W_{ui}, \widehat \theta_u, \widehat h_u(Z_{ui})), $$ and $\widehat J_u$ is a suitable estimator of $J_u$.\footnote{We do not discuss the estimation of $J_u$ since it is often a problem-specific matter. In Section 3, $J_u$ was equal to minus the identity matrix, so we did not need to estimate it.} The bootstrap law is computed by drawing $(\xi_i)_{i=1}^n$ conditional on the data.

The following theorem shows that the multiplier bootstrap provides a valid approximation to the large sample law of $\sqrt{n}(\widehat \theta_u- \theta_u)_{u \in \mathcal{U}}$.

We next derive the large sample distribution and validity of the multiplier bootstrap for the estimator $\widehat \Delta := \phi(\widehat \theta):= \phi( (\widehat \theta_u)_{u \in \mathcal{U}})$ of the functional $\Delta := \phi(\theta^0)= \phi( (\theta_u)_{u \in \mathcal{U}} )$ using the functional delta method. The functional $\theta^0 \longmapsto \phi(\theta^0)$ is defined as a uniformly Hadamard differentiable transform of $\theta^0 = (\theta_u)_{u \in \mathcal{U}}$. The following result gives the large sample law of $\sqrt{n} (\widehat \Delta - \Delta)$, the properly normalized estimator. It also shows that the bootstrap law of $ \sqrt{n} (\widehat \Delta^* - \widehat \Delta),$ computed conditionally on the data, is consistent for the large sample law of $\sqrt{n} (\widehat \Delta - \Delta)$. Here $\widehat \Delta^* := \phi(\widehat \theta^*) = \phi ( (\widehat \theta^*)_{u \in \mathcal{U}})$ is the bootstrap version of $\widehat \Delta$, and $\widehat \theta^*_u = \widehat \theta_u + n^{-1} \sum_{i=1}^n \xi_i \widehat \psi_u(W_i)$ is the multiplier bootstrap version of $\widehat \theta_u$ defined via equation ((ref)).

To derive Theorem (ref), we strengthen the usual notion of Hadamard differentiability to a uniform notion introduced in Definition (ref). Theorems (ref) and (ref) show that this uniform Hadamard differentiability is sufficient to guarantee the validity of the functional delta uniformly in $P$. These new uniform functional delta method theorems may be of independent interest.

Theory: Lasso and Post-Lasso for Functional Response Data

In this section, we provide results for Lasso and Post-Lasso estimators with function-valued outcomes and linear or logistic links. As these results are of interest beyond the context of estimation of nuisance functions for moment condition problems or treatment effects estimation, we present this section in a way that leaves it autonomous with respect to the rest of the paper.

The generic setting with function-valued outcomes

Consider a data generating process with a functional response variable $(Y_{u})_{u\in \mathcal{U}}$ and observable covariates $X$ satisfying for each $u\in \mathcal{U}$,

where $f:\mathcal{X}\to\mathbb{R}^p$ is a set of $p$ measurable transformations of the initial controls $X$, $\theta_u$ is a $p$-dimensional vector, $r_u$ is an approximation error, and ${\Lambda}$ is a fixed known link function. The notation in this section differs from the rest of the paper with $Y_u$ and $X$ denoting a generic response and a generic vector of covariates to facilitate the application of these results to other contexts. We only consider the linear link function, ${\Lambda}(t) = t$, and the logistic link function, ${\Lambda}(t)=\exp(t)/\{1+\exp(t)\}$, in detail.

Considering the logistic link is useful when the functional response is binary, though the linear link can be used in that case as well under some conditions. For example, it is useful for estimating a high-dimensional generalization of the distributional regression models considered in CFM where the response variable is the continuum $(Y_u = 1( Y \leqslant u))_{u \in \mathcal{U}}$. Even though we focus on these two cases we note that the principles discussed here apply to many other $M$-estimators with convex (or approximately convex) criterion functions. In the remainder of the section, we discuss and establish results for $\ell_1$-penalized and post-model selection estimators of $(\theta_u)_{u \in \mathcal{U}}$ that hold uniformly over $u\in\mathcal{U}$.

Throughout the section, we assume that $u\in \mathcal{U} \subset [0,1]^{d_u}$ and that we have $n$ i.i.d. observations from d.g.p.'s where ((ref)) holds, $\{( Y_{ui})_{u\in\mathcal{U}}, X_i)\}_{i=1}^n$, available for estimating $(\theta_u)_{u\in\mathcal{U}}$. For each $u\in\mathcal{U}$, penalty level $\lambda$, and diagonal matrix of penalty loadings $\widehat\Psi_u,$ we define the Lasso estimator as

where $M(y,t) = \frac{1}{2}(y-{\Lambda}(t))^2$ in the case of linear regression, and $M(y,t) = -\{1(y=1)\log {\Lambda}(t) + 1(y=0)\log(1-{\Lambda}(t))\}$ in the case of the logistic link function for binary response data. For each $u\in\mathcal{U}$, the Post-Lasso estimator based on a set of covariates $\widetilde T_u$ is then defined as

where the set $\widetilde T_u$ contains $\mathrm{supp}(\widehat\theta_u)$ and may also contain additional variables deemed as important.\footnote{The total number of additional variables $\widehat s_a$ should also obey the same growth conditions that $s$ obeys. For example, if the additional variables are chosen so that $\widehat s_a \lesssim \|\widehat\theta_u\|_0$ the growth condition is satisfied with probability going to one for the designs covered by Assumptions (ref) and (ref). See also BelloniChernozhukovHansen2011 for a discussion on choosing additional variables.} We will set $\widetilde T_u = \mathrm{supp}(\widehat\theta_u)$ unless otherwise noted.

The chief departure between the analysis when $\mathcal{U}$ is a singleton and the functional response case is that the penalty level needs to be set to control selection errors uniformly over $u\in\mathcal{U}$. To do so, we will set $\lambda$ so that with high probability

where $c>1$ is a fixed constant. When $\mathcal{U}$ is a singleton the strategy above is similar to BickelRitovTsybakov2009, BC-PostLASSO, and BCW-SqLASSO, who use an analog of ((ref)) to derive the properties of Lasso and Post-Lasso. When $\mathcal{U}$ is not a singleton, this strategy was first employed in the context of $\ell_1$-penalized quantile regression processes by BC-SparseQR.

To implement ((ref)), we propose setting the penalty level as

where ${d_u}$ is the dimension of $\mathcal{U}$, $1-\gamma$ with $\gamma = o(1)$ is a confidence level associated with the probability of event ((ref)), and $c>1$ is a slack constant.\footnote{When the set $\mathcal{U}$ is a singleton, one can use the penalty level in ((ref)) with ${d_u} = 0$. This choice corresponds to that used in BelloniChernozhukovHansen2011.} When implementing the estimators, we set $c=1.1.$ and $\gamma = .1/\log(n)$, which is theoretically motivated and practically tested in an extensive set of simulation experiments in BelloniChernozhukovHansen2011. In addition to the penalty parameter $\lambda$, we also need to construct a penalty loading matrix $\widehat\Psi_u = {\rm diag}(\{\widehat l_{u j}, j=1,\ldots,p\})$. This loading matrix can be formed according to the following iterative algorithm.

Properties of a Continuum of Lasso and Post-Lasso: Linear Link

We provide sufficient conditions for establishing good performance of the estimators discussed above when the linear link function is used. In the statement of the following assumption, $\delta_n\searrow 0$ and $\Delta_n\searrow 0$ are fixed sequences approaching zero from above at a speed at most polynomial in $n$ (for example, $\delta_n \geqslant 1/n^c$ for some $c > 0$), $\ell_n := \log n$, and $c, C, \kappa', \kappa''$ and $\nu \in (0,1]$ are positive finite constants.

Assumption (ref) is only a set of sufficient conditions. The finite sample results in the Supplementary Appendix allow for more general conditions (for example, ${d_u}$ can grow with the sample size). We verify that the more technical conditions in Assumption (ref)(iv)(b) hold in a variety of cases, see Lemma (ref) in Appendix (ref) in the Supplementary Appendix. Under Assumption (ref), we establish results on the performance of the estimators ((ref)) and ((ref)) for the linear link function case that hold uniformly over $u \in \mathcal{U}$ and $P \in \mathcal{P}_n$.

We note that the performance bounds are exactly of the type used in Assumption (ref) (see also Assumption (ref) in the Supplementary Appendix). Indeed, under the condition $s^2\log^2(p\vee n) \log^2 n \leqslant \delta_n n$, the rate of convergence established in Theorem (ref) yields $\sqrt{s\log(p\vee n)/n} \leqslant o( n^{-1/4})$.

Properties of Lasso and Post-Lasso Estimators: Logistic Link

We provide sufficient conditions to state results on the performance of the estimators discussed above for the logistic link function. Consider the fixed sequences $\delta_n\searrow 0$ and $\Delta_n\searrow 0$ approaching zero from above at a speed at most polynomial in $n$, $\ell_n := \log n$, and the positive finite constants $c$, $C$, $\kappa'$, $\kappa''$, and $\underline{c} \leqslant 1/2$.

The following result characterizes the performance of the estimators ((ref)) and ((ref)) for the logistic link function case under Assumption (ref).

Application: the Effect of 401(k) Participation on Asset Holdings

As a practical illustration of the methods developed in this paper, we consider estimation of the effect of 401(k) eligibility and participation on accumulated assets as in abadie:401k and CH401k. Our goal here is to illustrate the estimation results and inference statements and to make the following points that underscore our theoretical findings: 1) In a low-dimensional setting, where the number of controls is low and therefore there is no need for selection, our robust post-selection inference methods perform well. That is, the results of our methods agree with the results of standard methods that do not employ any selection. 2) In a high-dimensional setting, where there are (moderately) many controls, our post-selection inference methods perform well, producing well-behaved estimates and confidence intervals compared to the erratic estimates and confidence intervals produced by standard methods that do not employ selection as a means of regularization. 3) Finally, in a very high-dimensional setting, where the number of controls is comparable to the sample size, the standard methods break down completely, while our methods still produce well-behaved estimates and confidence intervals. These findings are in line with our theoretical results about uniform validity of our inference methods.

The key problem in determining the effect of participation in 401(k) plans on accumulated assets is saver heterogeneity coupled with the fact that the decision to enroll in a 401(k) is non-random. It is generally recognized that some people have a higher preference for saving than others. It also seems likely that those individuals with high unobserved preference for saving would be most likely to choose to participate in tax-advantaged retirement savings plans and would tend to have otherwise high amounts of accumulated assets. The presence of unobserved savings preferences with these properties then implies that conventional estimates that do not account for saver heterogeneity and endogeneity of participation will be biased upward, tending to overstate the savings effects of 401(k) participation.

To overcome the endogeneity of 401(k) participation, abadie:401k and CH401k adopt the strategy detailed in Poterba, Venti, and Wise pvw:94,pvw:95,pvw:nber96,pvw:01 and benjamin, who used data from the 1991 Survey of Income and Program Participation and argue that eligibility for enrolling in a 401(k) plan in this data can be taken as exogenous after conditioning on a few observables of which the most important for their argument is income. The basic idea of their argument is that, at least around the time 401(k)'s initially became available, people were unlikely to be basing their employment decisions on whether an employer offered a 401(k) but would instead focus on income. Thus, eligibility for a 401(k) could be taken as exogenous conditional on income, and the causal effect of 401(k) eligibility could be directly estimated by appropriate comparison across eligible and ineligible individuals.\footnote{Poterba, Venti, and Wise pvw:94,pvw:95,pvw:nber96,pvw:01 and benjamin all focus on estimating the effect of 401(k) eligibility, the intention to treat parameter. Also note that there are arguments that eligibility should not be taken as exogenous given income; see, for example, engen and engen:gale.} {abadie:401k, CH401k, and ORR:401k} use this argument for the exogeneity of eligibility conditional on controls to argue that 401(k) eligibility provides a valid instrument for 401(k) participation and employ IV methods to estimate the effect of 401(k) participation on accumulated assets.

As a complement to the work cited above, we estimate various treatment effects of 401(k) participation on financial wealth using high-dimensional methods. A key component of the argument underlying the exogeneity of 401(k) eligibility is that eligibility may only be taken as exogenous after conditioning on income. Both abadie:401k and CH401k adopt this argument but control only for a small number of terms. One might wonder whether the small number of terms considered is sufficient to adequately control for income and other related confounds. At the same time, the power to learn anything about the effect of 401(k) participation decreases as one controls more flexibly for confounds. The methods developed in this paper offer one resolution to this tension by allowing us to consider a very broad set of controls and functional forms under the assumption that among the set of variables we consider there is a relatively low-dimensional set that adequately captures the effect of confounds. This approach is more general than that pursued in previous research which implicitly assumes that confounding effects can adequately be controlled for by a small number of variables chosen ex ante by the researcher.

We use the same data as CH401k. The data consist of 9,915 observations at the household level drawn from the 1991 SIPP. We use net financial assets as the outcome variable, $Y$, in our analysis. Our treatment variable, $D$, is an indicator for having positive 401(k) balances; and our instrument, $Z$, is an indicator for being eligible to enroll in a 401(k) plan. The vector of raw covariates, $X$, consists of age, income, family size, years of education, a married indicator, a two-earner status indicator, a defined benefit pension status indicator, an IRA participation indicator, and a home ownership indicator. Further details can be found in CH401k.

We present detailed results for three different sets of controls $f(X)$. The first specification uses indicators of marital status, two-earner status, defined benefit pension status, IRA participation status, and home ownership status, second order polynomials in family size and education, a third order polynomial in age, and a quadratic spline in income with six break points\footnote{Specifically, we allow for income, income-squared, and then interact these two variables with seven dummies for the categories formed by the cut points.} (Quadratic Spline specification). The second specification augments the Quadratic Spline specification by interacting all the non-income variables with each term in the income spline (Quadratic Spline Plus Interactions specification). The final specification forms a larger set of potential controls by starting with all of the variables from the Quadratic Spline specification and forming all two-way interactions between all of the non-income variables. The set of main effects and interactions of all non-income variables is then fully interacted with all of the income terms (Quadratic Spline Plus Many Interactions specification).\footnote{The specifications are motivated by the original specification used in abadie:401k, benjamin, and CH401k allowing for data-dependent accommodation of nonlinearity. We report results based on the exact specification used in previous papers in the Supplementary Appendix.} The dimensions of the set of controls are thus 35, 311, and 1756 for the Quadratic Spline, Quadratic Spline Plus Interactions, and Quadratic Spline Plus Many Interactions specification, respectively. For methods that do not use variable selection, we use 32, 272, and 1526 variables resulting from removing terms that are perfectly collinear. We refer to the specification without interactions as low-$p$, to the specification with only income interactions as high-$p$, and to the specification with all two-way interactions further interacted with income as very-high-$p$.

We report a variety of results for each specification. Under the maintained assumption that 401(k) eligibility may be taken as exogenous after controlling for the variables defined in the preceding paragraph, we can use the methods of this paper to estimate intention to treat effects of 401(k) eligibility by setting 401(k) eligibility as $D = Z$. We report the estimated average intention to treat and average intention to treat on the treated as the ATE and ATE-T, and we report estimates of quantile intention to treat and quantile intention to treat on the treated effects as QTE and QTE-T. We also directly apply the results of this paper to estimate effects of 401(k) participation, reporting estimates of the LATE, LATE-T, LQTE, and LQTE-T for each specification.\footnote{We note that because of one-sided compliance the local effects for the treated actually coincide with population effects for the treated; see frolich:melly.} For comparison, we also report estimates of the eligibility effect from the linear model without selection and with selection using the approach of BelloniChernozhukovHansen2011 and estimates of the participation effect from linear instrumental variables estimation without selection and with selection as in CHS:PnP.

Estimation of all these treatment effects depends on first-stage estimates of reduced form functions as detailed in Section (ref). We estimate reduced form functions where the outcome is continuous using ordinary least squares when no model selection is used or Post-Lasso when selection is used. We estimate reduced form functions where the outcome is binary by logistic regression when no model selection is used or Post-$\ell_1$-penalized logistic regression when selection is used. We only report selection-based estimates in the very-high-$p$ setting.\footnote{The estimated propensity score shows up in the denominator of the efficient moment conditions. As is conventional, we use trimming to keep the denominator bounded away from zero with trimming set to $10^{-12}$. Trimming occurs in the Quadratic Spline Plus Interactions (12 observations trimmed) and Quadratic Spline Plus Many Interactions specifications (9915 observations trimmed) when selection is not done. Trimming never occurs in the selection-based estimates in this example. We choose not to report unregularized estimates in the very-high-$p$ specification since all observations are trimmed and, in fact, have estimated propensity scores of either 0 or 1.} We refer to Appendix (ref) for detailed discussion of implementing our approach in this example.

Estimates of the ATE, ATE-T, LATE and LATE-T as well as the coefficient on 401(k) eligibility from the linear model and coefficient on 401(k) participation in the linear IV model are given in Table 1. In this table, we provide point estimates for each of the three sets of controls with and without variable selection. We report conventional heteroscedasticity consistent standard error estimates for the linear model and linear IV coefficient. For the ATE, ATE-T, LATE, and LATE-T, we report both analytic and multiplier bootstrap standard errors. The bootstrap standard errors are based on 500 bootstrap replications with mammen1993:bootstrap weights as multipliers.

Looking first at the two sets of standard error estimates for the average treatment effect estimates, we see that the bootstrap and analytic standard errors are quite similar and that one would not draw substantively different conclusions from using one versus the other. We also see that estimates of the effect of 401(k) eligibility using the linear model and estimates of the effect of 401(k) participation using the linear IV model are broadly consistent with each other across all specifications and regardless of whether or not variable selection is done. We also have that the estimates of the ATE, ATE-T, LATE, and LATE-T are very similar regardless of whether selection is used in the low-p Quadratic Spline specification. The ATE and ATE-T both indicate a positive and significant average effect of 401(k) eligibility; and the LATE and LATE-T suggest positive and significant effects of 401(k) participation for compliers. The similarity in the low-p case is reassuring as it illustrates that there is little impact of variable selection relative to simply including everything in a low-dimensional setting.\footnote{In the low-dimensional setting, using all available controls is semi-parametrically efficient and allows uniformly valid inference. Thus, the similarity between the results in this case is an important feature of our method which results from our reliance on low-bias moment functions and sensible variable selection devices to produce semi-parametrically efficient estimators and uniformly valid inference statements following model selection.}

We observe somewhat different results in the Quadratic Spline Plus Interactions specification. For both the ATE and the LATE in the Quadratic Spline Plus Interactions case, we see a substantially larger point estimate without selection than with selection, with the selection results being similar to those obtained in the low-p case. Along with the larger point estimate, we also see that the estimated standard errors in the no selection case for the ATE and LATE are roughly three times larger than the standard errors in the selection case. For the ATE-T and LATE-T in the Quadratic Spline Plus Interactions case, point estimates following selection are notably smaller than without selection but estimated standard errors after selection are somewhat larger. We note that one might suspect estimated standard errors for all of the estimators without selection to be substantially downward biased in this case due to the use of many control variables without regularization as in CJN:PLMStandardError. Finally, we see a large difference in the Orthogonal Polynomials Plus Many Interactions Specifications as estimates cannot even be computed reliably without selection due to severe overfitting: The estimated propensity score is either 0 or 1 for every observation.

We provide estimates of the QTE and QTE-T in Figure 1 and estimates of the LQTE and LQTE-T in Figure 2. The left column of Figure 1 gives results for the QTE, and the right column displays the results for the QTE-T. Similarly, the left and right columns of Figure 2 provide the LQTE and LQTE-T respectively. We give the results for the Quadratic Spline, Quadratic Spline Plus Interactions, and Quadratic Spline Plus Many Interactions specification in the top row, middle row, and bottom row respectively. In each graphic, we use solid lines for point estimates and report uniform 95% confidence intervals with dashed lines.

Looking across the figures, we see a similar pattern to that seen for the estimates of the average effects in that the selection-based estimates are stable across all specifications and are very similar to the estimates obtained without selection from the baseline low-$p$ Quadratic Spline specification. In the more flexible Quadratic Spline plus Interactions specification, the estimates that do not make use of selection behave somewhat erratically. This erratic behavior is especially apparent in the estimated LQTE of 401(k) participation where we observe that small changes in the quantile index may result in large swings in the point estimate of the LQTE and estimated standard errors are quite large. Again, this erratic behavior is likely due to overfitting due to the large set of variables considered. As with the average effects, estimated quantile effects without selection in the Quadratic Spline Plus Many Interactions specification are not reported as the estimated propensity score is always 0 or 1.

If we focus on the LQTE and LQTE-T estimated from variable selection methods, we find that 401(k) participation has a small impact on accumulated net total financial assets at low quantiles while appearing to have a larger impact at high quantiles. Looking at the uniform confidence intervals, we can see that this pattern is statistically significant at the 5% level and that we would reject the hypothesis that 401(k) participation has no effect and reject the hypothesis of a constant treatment effect more generally.

It is also worth discussing the results of the variable selection briefly as well. Due to the number of models and variable selection steps taken, especially in computing quantile effects, it is not practical to give a complete accounting of the selected variables here. Rather, we note that for the linear model, linear IV, ATE, and LATE results, we select between two and 22 variables depending on the specification of controls and left-hand-side variable. The median number of variables selected for the QTE and LQTE results, where the median is taken across index values $u$, across the different specifications of controls and left-hand-side variables varies between one and 11. There is considerable variability in the number of variables selected across $u$ though, ranging from a minimum of no variables selected to a maximum of 237 selected variables.\footnote{Having more than 100 variables selected occurs in the very high dimensional setting when the outcome in the penalized regression is $\mathbf{1}_0(D) Y_u$ for the six lowest values of $u$ among the subset of households eligible for 401(k)'s and for the six highest values of $u$ among the subset of households that are not eligible for 401(k)'s.} The selected variables themselves mostly correspond to capturing the effect of income. For example, the union of the variables selected in forming each of the reduced form quantities used for estimating the LATE in the Quadratic Spline Plus Many Interactions specification consists of 36 variables, only four of which do not include income.\footnote{Let $i_1$ be the indicator for income in the first income category, and define $i_2-i_7$ similarly. Let $db$ be the defined benefit dummy, $ira$ be the IRA dummy, $hown$ be the home ownership dummy, $mar$ be the married dummy, $te$ be the two-earner household dummy, $ed$ be years of schooling, and $fsize$ be family size. The exact identities of the variables selected for modeling any reduced form quantity used in estimating the LATE in the very-high-dimensional case are $i_1$, $i_2$, $i_3$, $income*i_3$, $income^2*i_6$, $db$, $ira*hown$, $age*ira$, $ed*ira$, $i_1*fsize$, $i_1*fsize^2*db$, $i_2*fsize$, $i_2*fsize^2*db$, $i_3*age^3$, $i_3*fsize$, $i_3*mar$, $i_3*fsize*te$, $i_3*fsize*mar$, $i_4*fsize$, $i_4*te$, $i_4*fsize*te$, $i_4*mar*te$, $i_4*ed*fsize$, $i_5*ed^2*te$, $i_5*fsize*te$, $income*ira$, $income*hown$, $income*mar*hown$, $income*te*hown$, $income*fsize*ira$, $income*fsize^2*ira$, $income*i_1*fsize$, $income^2*ira*hown$, $income^2*ed^2*te$, $income^2*i_3*fsize$, and $income^2*i_6*hown$.} This pattern of largely selecting terms that are direct income effects or interactions of income with other variables holds up across the specifications considered.

It is interesting that our results are similar to those in CH401k despite allowing for a much richer set of controls. The fact that we allow for a rich set of controls but produce similar results to those previously available lends further credibility to the claim that previous work controlled adequately for the available observables.\footnote{Of course, the estimates are still not valid causal estimates if one does not believe that 401(k) eligibility can be taken as exogenous after controlling for income and the other included variables.} Finally, it is worth noting that this similarity is not mechanical or otherwise built in to the procedure. For example, applications in BellChenChernHans:nonGauss and BelloniChernozhukovHansen2011 use high-dimensional variable selection methods and produce sets of variables that differ substantially from intuitive baselines.