Inference for High-Dimensional Sparse Econometric Models

Inference for High-Dimensional Sparse Econometric Models

Introduction

We consider linear, high dimensional sparse (HDS) regression models in econometrics. The HDS regression model allows for a large number of regressors, $p$, which is possibly much larger than the sample size, $n$, but imposes that the model is sparse. That is, we assume only $s \ll n$ of these regressors are important for capturing the main features of the regression function. This assumption makes it possible to estimate HDS models effectively by searching for approximately the right set of regressors. In this article, we review estimation methods for HDS models that make use of $\ell_1$-penalization and then provide a set of novel inference results. We also provide empirical examples that illustrate the potential wide applicability of HDS models and methods in econometrics.

The motivation for considering HDS models comes in part from the wide availability of data sets with many regressors. For example, the American Housing Survey records prices as well as a multitude of features of houses sold; and scanner data-sets record prices and numerous characteristics of products sold at a store or on the internet. HDS models are also partly motivated by the use of series methods in econometrics. Series methods use many constructed or series regressors -- regressors formed as transformation of elementary regressors -- to approximate regression functions. In these applications, it is important to have parsimonious yet accurate approximation of the regression function. One way to achieve this is to use the data to select a small of number of informative terms from among a very large set of control variables or approximating functions. In this article, we formally discuss doing this selection and estimating the regression function.

We organize the article as follows. In the next section, we introduce the concepts of sparse and approximately sparse regression models in the canonical context of modeling a conditional mean function and motivate the use of HDS models via an empirical and analytical examples. In Section (ref), we discuss some principal estimation methods and mention extensions of these methods to applications beyond conditional mean models. We discuss some key estimation results for HDS methods and mention various extensions of these results in Section (ref). We then develop HDS models and methods in instrumental variables models with many instruments in Section (ref) and a partially linear model with many series terms in Section (ref), with the main emphasis given to inference. Finally, we present two empirical examples which motivate the use of these methods in Section (ref).

Notation. We allow for the models to change with the sample size, i.e. we allow for array asymptotics. In particular we assume that $p=p_n$ grows to infinity as $n$ grows, and $s =s_n$ can also grow with $n$, although we require that $s \log p = o(n)$. Thus, all parameters are implicitly indexed by the sample size $n$, but we omit the index to simplify notation. We also use the following empirical process notation, ${\mathbb{E}_n}[f] = {\mathbb{E}_n}[f(z_i)] = \sum_{i=1}^n f(z_i)/n.$ The ${l}_2$-norm is denoted by $\|\cdot\|$, and the ${l}_0$-norm, $\|\cdot\|_0$, denotes the number of non-zero components of a vector. We use $\| \cdot \|_{\infty}$ to denote the maximal element of a vector. Given a vector $\delta \in {\Bbb{R}}^p$, and a set of indices $T \subset \{1,\ldots,p\}$, we denote by $\delta_T \in {\Bbb{R}}^p$ the vector in which $\delta_{Tj} = \delta_j$ if $j\in T$, $\delta_{Tj}=0$ if $j \notin T$. We use the notation $(a)_+ = \max\{a,0\}$, $a \vee b = \max\{ a, b\}$ and $a \wedge b = \min\{ a , b \}$. We also use the notation $a \lesssim b$ to denote $a \leqslant c b$ for some constant $c>0$ that does not depend on $n$; and $a\lesssim_P b$ to denote $a=O_P(b)$. For an event $E$, we say that $E$ wp $\to$ 1 when $E$ occurs with probability approaching one as $n$ grows.

Sparse and Approximately Sparse Regression Models

In this section we review the modeling foundations for HDS methods and provide motivating examples with emphasis on applications in econometrics. First, let us consider the following parametric linear regression model: $$y_i = x_i'\beta_0 + \epsilon_i, \ \ \epsilon_i \sim N(0, \sigma^2), \ \ \beta_0 \in \Bbb{R}^p, \ \ i=1,\ldots,n$$ $$ T = {\rm support}(\beta_0) \text{ has } s \text{ elements where } s< n, $$ where $p> n$ is allowed, $T$ is unknown, and regressors $X=[x_1,\ldots, x_n]'$ are fixed. We assume Gaussian errors to simplify the presentation of the main ideas throughout the article, but note that this assumption can be eliminated without substantially altering the results. It is clear that simply regressing $y$ on all $p$ available $x$ variables is problematic when $p$ is large relative to $n$ which motivates consideration of models that impose some regularization on the estimation problem.

The key assumption that allows effective use of this large set of covariates is sparsity of the model of interest. Sparsity refers to the condition that only $s \ll n$ elements of $\beta_0$ are non-zero but allows the identities of these elements to be unknown. Sparsity can be motivated on economic grounds in situations where a researcher believes that the economic outcome could be well-predicted by a small (relative to the sample size) number of factors but is unsure about the identity of the relevant factors. Note that we allow $s=s_n$ to grow with $n$, as mentioned in the notation section, although $s \log p = o(n)$ will be required for consistency. This simple sparse model substantially generalizes the classical parametric linear model by letting the identities, $T,$ of the relevant regressors be unknown. This generalization is useful in practice since it is problematic to assume that we know the identities of the relevant regressors in many examples.

The previous model is simple and allows us to convey the essential ideas of the sparsity-based approach. However, it is unrealistic in that it presumes exact sparsity or that, after accounting for $s$ main regressors, the error in approximating the regression function is zero. We shall make no formal use of the previous model, but instead use a much more general, approximately sparse or nonparametric model. In this model, all of the regressors potentially have a non-zero contribution to the regression function, but no more than $s$ unknown regressors are needed for approximating the regression function with a sufficient degree of accuracy.

We formally define the approximately sparse model as follows.

Condition ASM. We have data $\{(y_i,z_i), i=1,\ldots,n\}$ that for each $n$ obey the regression model:

where $y_i$ is the outcome variable, $z_i$ is a $k_z$-vector of elementary regressors, $f(z_i)$ is the regression function, and $\epsilon_i$ are i.i.d. disturbances. Let $x_i=P(z_i)$, where $P(z_i)$ is a vector of dimension $p=p_n$, that contains a dictionary of possibly technical transformations of $z_i$, including a constant. The values $x_1,\ldots, x_n$ are treated fixed, and normalized so that ${\mathbb{E}_n}[x_{ij}^2] = 1$ for $j=1,\ldots,p$. The regression function $f(z_i)$ admits the approximately sparse form, namely there exists $\beta_0$ such that

where $s=s_n=o(n/\log p)$ and $K$ is a constant independent of $n$. }

In the set-up we consider the fixed design case, which covers random sampling as a special case where $x_1,\ldots, x_n$ represent a realization of this sample on which we condition throughout. The vector $x_i=P(z_i)$ can include polynomial or spline transformations of the original regressors $z_i$ see, e.g., \citeasnoun{newey:series} and \citeasnoun{chen:Chapter} for various examples of series terms. The approximate sparsity can be motivated similarly to \citeasnoun{newey:series}, who assumes that the first $s=s_n$ series terms can approximate the nonparametric regression function well. Condition ASM is more general in that it does not impose that the most important $s=s_n$ terms in the approximating dictionary are the first $s$ terms; in fact, the identity of the most important terms is treated as unknown. We note that in the parametric case, we may naturally choose $x_i'\beta_0=f(z_i)$ so that $r_i=0$ for all $i=1,\ldots,n$. In the nonparametric case, we may think of $x_i'\beta_0$ as any sparse parametric model that yields a good approximation to the true regression function $f(z_i)$ in equation ((ref)) so that $r_i$ is “small” relative to the conjectured size of the estimation error. Given ((ref)), our target in estimation is the parametric function $x_i'\beta_0$, where we can call $$ T := {\rm support}(\beta_0) $$ the “true" model. Here we emphasize that the ultimate target in estimation is, of course, $f(z_i)$. The function $x_i'\beta_0$ is simply a convenient intermediate target introduced so that we can approach the estimation problem as if it were parametric. Indeed, the two targets, $f(z_i)$ and $x_i'\beta_0$, are equal up to the approximation error $r_i$. Thus, the problem of estimating the parametric target $x_i'\beta_0$ is equivalent to the problem of estimating the nonparametric target $f(z_i)$ modulo approximation errors.

One way to explicitly construct a good approximating model $\beta_0$ for ((ref)) is by taking $\beta_0$ as the solution to

We can call ((ref)) the oracle problem,\footnote{By definition the oracle knows the risk function of any estimator, so it can compute the best sparse least square estimator. Under some mild condition the problem of minimizing prediction risk amongst all sparse least square estimators is equivalent to the problem written here; see, e.g., \citeasnoun{BC-LectureNotes}.} and so we can call $T = {\rm support}(\beta_0)$ the oracle model. Note that we necessarily have that $s=\|\beta_0\| \leqslant n$. The oracle problem ((ref)) balances the approximation error ${\mathbb{E}_n} [(f(z_i) - x_i'\beta)^2]$ over the design points with the variance term $\sigma^2 \|\beta\|_0/n$, where the latter is determined by the number of non-zero coefficients in $\beta$. Letting $ c^2_s:= {\mathbb{E}_n}[r^2_i] = {\mathbb{E}_n} [(f(z_i) - x_i'\beta_0)^2]$ denote the squared error from approximating values $f(z_i)$ by $x_i'\beta_0$, the quantity $ c^2_s + \sigma^2 s/n$ is the optimal value of ((ref)). In common nonparametric problems, such as the one described below, the optimal solution in ((ref)) would balance the approximation error with the variance term giving that $c_s \leqslant K \sigma \sqrt{s/n}.$ Thus, we would have $ \sqrt{c^2_s + \sigma^2 s/n} \lesssim \sigma \sqrt{s/n},$ implying that the quantity $\sigma \sqrt{s/n}$ is the ideal goal for the rate of convergence. If we knew the oracle model $T$, we would achieve this rate by using the oracle estimator, the least squares estimator based on this model. Of course, we do not generally know $T$ since we do not observe the $f(z_i)$'s and thus cannot attempt to solve the oracle problem ((ref)). Since $T$ is unknown, we will not generally be able to achieve the exact oracle rates of convergence, but we can hope to come close to this rate.

Before considering estimation methods, a natural question is whether exact or approximate HDS models make sense in econometric applications. In order to answer this question, it is helpful to consider the following two examples in which we abstract from estimation completely and only ask whether it is possible to accurately describe some structural econometric function $f(z)$ using a low-dimensional approximation of the form $P(z)'\beta_0$.

Example 1: Sparse Models for Earning Regressions. In this example we consider a model for the conditional expectation of log-wage $y_i$ given education $z_i$, measured in years of schooling. We can expand the conditional expectation of wage $y_i$ given education $z_i$:

using some dictionary of approximating functions $P(z_i) = (P_1(z_i),\ldots, P_p(z_i))'$, such as polynomial or spline transformations in $z_i$ and/or indicator variables for levels of $z_i$. In fact, since we can consider an overcomplete dictionary, the representation of the function using $P_1(z_i),\ldots, P_p(z_i)$ may not be unique, but this is not important for our purposes.

A conventional sparse approximation employed in econometrics is, for example,

where the $P_j$'s are low-order polynomials or splines, with typically one or two (linear or linear and quadratic) terms. Of course, there is no guarantee that the approximation error $\tilde r_i$ in this case is small or that these particular polynomials form the best possible $s$-dimensional approximation. Indeed, we might expect the function $E[y_i|z_i]$ to change rapidly near the schooling levels associated with advanced degrees, such as MBAs or MDs. Low-degree polynomials may not be able to capture this behavior very well, resulting in large approximation errors $\tilde r_i$.

A sensible question is then, “Can we find a better approximation that uses the same number of parameters?” More formally, can we construct a much better approximation of the sparse form

for some regressor indices $k_1,\ldots,k_s$ selected from $\{1,\ldots,p\}$? Since we can always include ((ref)) as a special case, we can in principle do no worse than the conventional approximation; and, in fact, we can construct ((ref)) that is much better, if there are some important higher-order terms in ((ref)) that are completely missed by the conventional approximation. Thus, the answer to the question depends strongly on the empirical context.

Consider for example the earnings of prime age white males in the 2000 U.S. Census see, e.g., \citeasnoun{ACF2006}. Treating this data as the population data, we can compute $f(z_i)=E[y_i|z_i]$ without error. Figure (ref) plots this function. We then construct two sparse approximations and also plot them in Figure (ref). The first is the conventional approximation of the form ((ref)) with $P_1, \ldots, P_s$ representing polynomials of degree zero to $s-1$ ($s=5$ in this example). The second is an approximation of the form ((ref)), with $P_{k_1}$, \ldots, $P_{k_s}$ consisting of a constant, a linear term, and three linear splines terms with knots located at 16, 17, and 19 years of schooling. We find the latter approximation automatically using the $\ell_1$-penalization or Lasso methods discussed below,\footnote{The set of functions considered consisted of 12 linear splines with various knots and monomials of degree zero to four. Note that there were only 12 different levels of schooling.} although in this special case we could construct such an approximation just by eye-balling Figure (ref) and noting that most of the function is described by a linear function with a few abrupt changes that can be captured by linear spline terms that induce large changes in slope near 17 and 19 years of schooling. Note that an exhaustive search for a low-dimensional approximation in principle requires looking at a very large set of models. Methods for HDS models, such as $\ell_1$-penalized least squares (Lasso), which we employed in this example, are designed to avoid this search. \qed

Example 2: Series approximations and Condition ASM. It is clear from the statement of Condition ASM that this expansion incorporates both substantial generalizations and improvements over the conventional series approximation of regression functions in \citeasnoun{newey:series}. In order to explain this consider the set $\{P_j(z), j\geqslant 1\}$ of orthonormal basis functions on $[0,1]^d$, e.g. orthopolynomials, with respect to the Lebesgue measure. Suppose $z_i$ have a uniform distribution on $[0,1]^d$ for simplicity.\footnote{The discussion in this example continues to apply when $z_i$ has a density that is bounded from above and away from zero on $[0,1]^d$.} Assuming ${\mathrm{E}}[f^2(z_i)] < \infty$, we can represent $f$ via a Fourier expansion, $ f(z) = \sum_{j=1}^\infty \delta_j P_j(z),$ where $\{\delta_j, j \geqslant 1\}$ are Fourier coefficients that satisfy $\sum_{j=1}^\infty \delta_j^2 < \infty$.

Let us consider the case that $f$ is a smooth function so that Fourier coefficients feature a polynomial decay $\delta_j \propto j^{-\nu}$, where $\nu$ is a measure of smoothness of $f$. Consider the conventional series expansion that uses the first $K$ terms for approximation, $f(z) = \sum_{j=1}^K \beta_{0j}P_j(z) + a_c(z)$, with $\beta_{0j} = \delta_j$. Here $a_c(z_i)$ is the approximation error which obeys $ \sqrt{{\mathbb{E}_n}[a^{2}_{c}(z_i)]} \lesssim_P \sqrt{{\mathrm{E}}[a^{2}_{c}(z_i)]} \lesssim K^{\frac{-2\nu+1}{2}}$. Balancing the order $K^{\frac{-2\nu+1}{2}}$ of approximation error with the order $\sqrt{K/n}$ of the estimation error gives the oracle-rate-optimal number of series terms $s = K \propto n^{1/2\nu}$, and the resulting oracle series estimator, which knows $s$, will estimate $f$ at the oracle rate of $n^{\frac{1-2\nu}{4\nu}}$. This also gives us the identity of the most important series terms $T_{} = \{1,...,s\}$, which are simply the first $s$ terms. We conclude that Condition ASM holds for the sparse approximation $f(z) = \sum_{j=1}^p \beta_{0j}P_j(z) + a_{}(z)$, with $\beta_{0j} = \delta_j$ for $j \leqslant s$ and $\beta_{0j} = 0$ for $s+ 1\leqslant j \leqslant p$, and $a_{}(z_i)=a_{c}(z_i)$, which coincides with the conventional series approximation above, so that $\sqrt{{\mathbb{E}_n}[a^2_{}(z_i)]} \lesssim_P \sqrt{s/n}$ and $\|\beta_{0}\|_0 \leqslant s$.

Next suppose that Fourier coefficients feature the following pattern $\delta_j = 0$ for $j \leqslant M$ and $\delta_j \propto (j-M)^{-\nu}$ for $j > M$. Clearly in this case the standard series approximation based on the first $K \leqslant M$ terms, $\sum_{j=1}^K \delta_jf_j(z)$, has no predictive power for $f(z)$, and the corresponding standard series estimator based on the first $K$ terms therefore fails completely.\footnote{This is not merely a finite sample phenomenon but is also accommodated in the asymptotics since we expressly allow for array asymptotics; i.e. the underlying true model could change with $n$. Recall that we omit the indexing by $n$ for ease of notation.} In contrast, Condition ASM is easily satisfied in this case, and the Lasso-based estimators will perform at a near-oracle level in this case. Indeed, we can use the first $p$ series terms to form the approximation $ f(z)= \sum_{j=1}^p \beta_{0j} P_j(z) + a_{}(z)$, where $ \beta_{0j} = 0$ for $j \leqslant M$ and $j > M + s$, $\beta_{0j} = \delta_j$ for $M+1 \leqslant j \leqslant M +s $ with $s \propto n^{1/2\nu}$, and $p$ such that $M + n^{1/2\nu} = o(p).$ Hence $\|\beta_{0}\|_0 = s$, and we have that $ \sqrt{{\mathbb{E}_n}[a^2_{}(z_i)]} \lesssim_P \sqrt{{\mathrm{E}}[a^2_{}(z_i)]} \lesssim \sqrt{s/n} \lesssim n^{\frac{1-2\nu}{4\nu}}$. \qed

Sparse Estimation Methods

$\ell_1$-penalized and post $\ell_1$-penalized estimation methods

In order to discuss estimation consider first, as a matter of motivation, the classical AIC/BIC type estimator Akaike1974,Schwarz1978 that solves the empirical (feasible) analog of the oracle problem: $$ \min_{\beta \in \Bbb{R}^p} {\mathbb{E}_n}[(y_i - x_i'\beta)^2] + \frac{\lambda}{n} \|\beta\|_{0}, $$ where $\lambda$ is a penalty level.\footnote{The penalty level $\lambda$ in the AIC/BIC type estimator needs to account for the noise since it observes $y_i$ instead of $f(z_i)$ unlike the oracle problem ((ref)).} This estimator has attractive theoretical properties. Unfortunately, it is computationally prohibitive since the solution to the problem may require solving $\sum_{k \leqslant n} \binom{p}{k}$ least squares problems.\footnote{Results on the computational intractability of this problem were established in \citeasnoun{Natarajan1995}, \citeasnoun{GeJiangYe2010} and \citeasnoun{ChenGeWangYe2011}.}

One way to overcome the computational difficulty is to consider a convex relaxation of the preceding problem, namely to employ a closest convex penalty -- the $\ell_1$ penalty -- in place of the $\ell_0$ penalty. This construction leads to the so called Lasso estimator $\widehat \beta$ T1996, defined as a solution for the following optimization problem:

where $\|\beta\|_{1} = \sum_{j=1}^p | \beta_j|$. The Lasso estimator is computationally attractive because it minimizes a convex function. A basic choice for penalty level suggested by \citeasnoun{BickelRitovTsybakov2009} is

where $c>1$ and $1-\gamma$ is a confidence level that needs to be set close to 1. The formal motivation for this penalty is that it leads to near-oracle rates of convergence of the estimator.

The penalty level specified above is not feasible since it depends on the unknown $\sigma$. \citeasnoun{BC-PostLASSO} propose to set

with $\widehat \sigma = \sigma + o_P(1)$ obtained via an iteration method defined in Appendix A, where $c>1$ and $1-\gamma$ is a confidence level.\footnote{Practical recommendations include the choice $c=1.1$ and $\gamma=.05$.} \citeasnoun{BC-PostLASSO} also propose the $X$-dependent penalty level:

where $$\Lambda(1-\gamma|X) = (1-\gamma)-\text{quantile of} \ \ n\|{\mathbb{E}_n}[x_ig_i]\|_\infty \mid X$$ where $X=[x_1,\ldots,x_n]'$ and $g_i$ are i.i.d. $N(0,1)$ , which can be easily approximated by simulation. We note that

so $\sqrt{2 n \log(2p/\gamma)}$ provides a simple upper bound on the penalty level. Note also that \citeasnoun{BellChenChernHans:nonGauss} formulate a feasible Lasso procedure for the case with heteroscedastic, non-Gaussian disturbances. We shall refer to the feasible Lasso method with the feasible penalty levels ((ref)) or ((ref)) as the Iterated Lasso. This estimator has statistical performance that is similar to that of the (infeasible) Lasso described above.

\citeasnoun{BCW-SqLASSO} propose a variant called the Square-root Lasso estimator $\widehat \beta$ defined as a solution to the following program:

with the penalty level

where $c>1$ and $$\widetilde \Lambda(1-\gamma|X)= (1-\gamma)-\text{quantile of } \ n\|{\mathbb{E}_n}[x_ig_i]\|_\infty/\sqrt{{\mathbb{E}_n}[g_i^2]} \mid X,$$ with $g_i \sim N(0,1)$ independent for $i=1,\ldots, n$. As with Lasso, there is also simple asymptotic option for setting the penalty level:

The main attractive feature of ((ref)) is that the penalty level $\lambda$ is independent of the value $\sigma$, and so it is pivotal with respect to that parameter. Nonetheless, this estimator has statistical performance that is similar to that of the (infeasible) Lasso described above. Moreover, the estimator is a solution to a highly tractable conic programming problem:

where the criterion function is linear in parameters $t$ and positive and negative components of $\beta$, while the constraint can be formulated with a second-order cone, informally known also as the “ice-cream cone".

There are several other estimators that make use of penalization by the $\ell_1$-norm. An important case includes the Dantzig selector estimator proposed and analyzed by \citeasnoun{CandesTao2007}. It also relies on $\ell_1$-regularization but exploits the notion that the residuals should be nearly uncorrelated with the covariates. The estimator is defined as a solution to:

where $\lambda = \sigma \Lambda(1-\gamma|X)$. In what follows we will focus our discussion on Lasso but virtually all theoretical results carry over to other $\ell_1$-regularized estimators including ((ref)) and ((ref)). We also refer to \citeasnoun{gautier:tsybakov} for a feasible Dantzig estimator that combines the square-root lasso method ((ref)) with the Dantzig method.

$\ell_1$-regularized estimators often have a substantial shrinkage bias. In order to remove some of this bias, we consider the post-model-selection estimator that applies ordinary least squares regression to the model $\widehat T$ selected by a $\ell_1$-regularized estimator $\widehat\beta$. Formally, set $$\widehat T = {\rm support}( \widehat \beta ) = \{ j \in \{1,\ldots,p\} \ : \ |\widehat\beta_j| > 0\},$$ and define the post model selection estimator $\widetilde \beta$ as

where $\widehat T^c = \{1,...,p\} \setminus \widehat T$. In words, the estimator is ordinary least squares applied to the data after removing the regressors that were not selected in $\widehat T$. When the $\ell_1$-regularized method used to select the model is Lasso (Square-root Lasso), the post-model-selection estimator is called Post-Lasso (Post-Square-root Lasso). If model selection works perfectly -- that is, $\widehat T = T$ -- then the post-model-selection estimator is simply the oracle estimator whose properties are well-known. However, perfect model selection is unlikely in many situations, so we are interested in the properties of the post-model-selection estimator when model selection is imperfect, i.e. when $\widehat T \neq T$, and are especially interested in cases where $T\nsubseteq \widehat T$. In Section (ref) we describe the formal properties of the Post-Lasso estimator.

Some Heuristics via Convex Geometry

Before proceeding to the formal results on estimation, it is useful to consider some heuristics for the $\ell_1$-penalized estimators and the choice of the penalty level. For this purpose we consider a parametric model, and a generic $\ell_1$-regularized estimator based on a differentiable criterion function $\widehat Q$:

where, e.g., $ \widehat Q(\beta) = {\mathbb{E}_n}[(y_i - x_i'\beta)^2]$ for Lasso and $ \widehat Q(\beta) = \sqrt{{\mathbb{E}_n}[(y_i - x_i'\beta)^2]}$ for Square-root Lasso. The key quantity in the analysis of ((ref)) is the score -- the gradient of $\widehat Q$ at the true value\footnote{In the case of a nonparametric model the score is similar to the gradient of $\widehat Q$ at $\beta_0$ but ignores the approximation errors $r_i$'s.}: $$ S = \nabla \widehat Q(\beta_0). $$ The score $S$ is the effective “noise" in the problem that should be dominated by the regularization. However we would like to make the regularization bias as small as possible. This reasoning suggests choosing the smallest penalty level $\lambda$ that is large enough to dominate the noise with high probability, say $1 - \gamma$, which yields

where $\Lambda$ is the maximal score scaled by $n$, and $c>1$ is a theoretical constant of \citeasnoun{BickelRitovTsybakov2009} that guarantees that the score is dominated. We note that the principle of setting $\lambda$ to dominate the score of the criterion function is a general principle that carries over to other convex problems with possibly non-differentiable criterion functions and that leads to the optimal -- near-oracle -- performance of $\ell_1$-penalized estimators. See, for instance, \citeasnoun{BC-SparseQR}.

It is useful to mention some simple heuristics for the principle ((ref)) which arise from considering the simplest case where none of the regressors are significant so that $\beta_0 =0$. We want our estimator to perform at a near-oracle level in all cases, including this case, but here the oracle estimator $ \beta^*$ sets $\beta^*= \beta_0 =0$. We also want $\widehat \beta = \beta_0 = 0$ in this case, at least with a high probability, say $1-\gamma$. From the subgradient optimality conditions for ((ref)), we must have $$ - S_j + \lambda/n > 0 \text{ and } S_j + \lambda/n > 0 \text{ for all } 1 \leqslant j \leqslant p$$ for this to be true. We can only guarantee this by setting the penalty level $\lambda/n$ such that $ \lambda > n \max_{1 \leqslant j\leqslant p} | S_j| = n \| S\|_{\infty}$ with probability at least $1-\gamma$. This is precisely the rule ((ref)) appearing above.

Finally, note that in the case of Lasso and Square-root Lasso we have the following expressions for the score: $$

$$ where $g_i$ are i.i.d. $N(0,1)$ variables. Note that the score for Square-root Lasso is pivotal, while the score for Lasso is not, as it depends on $\sigma$. Thus, the choice of the penalty level for Square-root Lasso need not depend on $\sigma$ to produce near-oracle performance for this estimator.

Beyond Mean Models

Most of the literature on high dimensional sparse models focuses on the mean regression model discussed above. Here we discuss methods that have been proposed to deal with quantile regression and generalized linear models in high-dimensional sparse settings. We assume i.i.d. sampling for $(y_i,x_i)$ in this subsection.

Quantile Regression

We consider a response variable $y_i$ and $p$-dimensional covariates $x_i$ such that the $u$-th conditional quantile function of $y_i$ given $x_i$ is given by

where $u \in (0,1)$ is quantile index of interest. Recall that the $u$-th conditional quantile $F^{-1}_{y_i|x_i}(u|x)$ is the inverse of the conditional distribution function $F_{y_i|x_i}(y|x)$ of $y_i$ given $x_i=x$. Suppose that the true model $\beta(u)$ has a sparse support: $$T_u = {\rm support}(\beta(u)) = \{ j \in \{1,\ldots,p\} \ : \ |\beta_j(u)|>0 \} $$ has only $s_u \leqslant s \leqslant n/\log(n \vee p)$ non-zero components.

The population coefficient $\beta(u)$ is known to be a minimizer of the criterion function

where $\rho_{u} (t) = (u - 1\{t\leqslant 0\})t$ is the asymmetric absolute deviation function; see \citeasnoun{Koenker:1978}. Given a random sample $(y_1,x_1),\ldots,(y_n,x_n)$, $\widehat\beta(u)$, the quantile regression estimator of $\beta(u)$, is defined as a minimizer of the empirical analog of ((ref)):

As before, in high-dimensional settings, ordinary quantile regression is generally not consistent, which motivates the use of penalization in order to remove all, or at least nearly all, regressors whose population coefficients are zero. The $\ell_1$-penalized quantile regression estimator $\widehat \beta (u)$ is a solution to the following optimization problem:

The criterion function in ((ref)) is the sum of the criterion function ((ref)) and a penalty function given by a scaled $\ell_1$-norm of the parameter vector.

In order to describe choice of the penalty level $\lambda$, we introduce the random variable

where $u_1,\ldots, u_n$ are i.i.d. uniform $(0,1)$ random variables, independently distributed from the regressors, $x_1,\ldots, x_n$. The random variable $\Lambda$ has a pivotal distribution conditional on $X= [x_1,\ldots, x_n]'$. Then, for $c>1$, \citeasnoun{BC-SparseQR} propose to set

and $1-\gamma$ is a confidence level that needs to be set close to 1.

The post-penalized estimator (post-$\ell_1$-QR) applies ordinary quantile regression to the model $\widehat T_u$ selected by the $\ell_1$-penalized quantile regression BC-SparseQR. Specifically, set $$\widehat T_u = {\rm support}( \widehat \beta(u) ) = \{ j \in \{1,\ldots,p\} \ : \ |\widehat\beta_j(u)| > 0\},$$ and define the post-penalized estimator $\widetilde \beta(u)$ as

which is just ordinary quantile regression removing the regressors that were not selected in the first step. \citeasnoun{BC-SparseQR} derive the basic properties of the estimators above; see also \citeasnoun{kato} for further important results in nonparametric setting, where group penalization is also studied.

Generalized Linear Models

From the discussion above, it is clear that $\ell_1$-regularized methods can be extended to other criterion functions $\widehat Q$ beyond least squares and quantile regression. $\ell_1$-regularized generalized linear models were considered in \citeasnoun{vdGeer}. Let $y \in {\Bbb{R}}$ denote the response variable and $x \in {\Bbb{R}}^p$ the covariates. The criterion function of interest is defined as $$\widehat Q(\beta) = \frac{1}{n}\sum_{i=1}^n h(y_i,x_i'\beta)$$ where $h$ is convex and $1$-Lipschitz with respect the second argument, $|h(y,t)-h(y,t')|\leqslant |t - t' |.$ We assume $h$ is differentiable in the second argument with derivative denoted $\nabla h$ to simplify exposition. Let the true model parameter be defined by $ \beta_0 \in \arg\min_{\beta\in {\Bbb{R}}^p} {\mathrm{E}}[ h(y_i,x_i'\beta) ]$, and consequently we have ${\mathrm{E}}[ x_i \nabla h(y_i,x_i'\beta_0) ] = 0$. The $\ell_1$-regularized estimator is given by the solution of $$ \min_{\beta\in {\Bbb{R}}^p} \widehat Q(\beta) + \frac{\lambda}{n} \|\beta\|_1. $$ Under high level conditions \citeasnoun{vdGeer} derived bounds on the excess forecasting loss, ${\mathrm{E}}[h(y_i,x_i'\widehat \beta)] - {\mathrm{E}}[h(y_i,x_i'\beta_0)]$, under sparsity-related assumptions, and also specialized the results to logistic regression, density estimation, and other problems.\footnote{Results in other norms of interest could also be derived, and the behavior of the post-$\ell_1$-regularized estimators would also be interesting to consider. This is an interesting venue for future work.} The choice of penalty parameter $\lambda$ derived in \citeasnoun{vdGeer} relies on using the contraction inequalities of \citeasnoun{LedouxTalagrandBook} in order to bound the score:

where $\xi_i$ are independent Rademacher random variables, $P(\xi_i=1)=P(\xi_i=-1)=1/2$. Then \citeasnoun{vdGeer} suggests further bounds on the right side of ((ref)). For efficiency reasons, we suggest simulating the $1-\gamma$ quantiles of the right side of ((ref)) conditional on regressors. In either way one can achieve the domination of “noise" $\lambda/n \geqslant c\|\nabla \widehat Q(\beta_0)\|_\infty$ with high probability. Note that since $h$ is 1-Lipschitz, this choice of the penalty level is pivotal.

Estimation Results for High Dimensional Sparse Models

Convergence Rates for Lasso and Post-Lasso

Having introduced Condition ASM and the target parameter defined via ((ref)), our task becomes to estimate $\beta_0$. We will focus on convergence results in the {\it prediction norm} for $\delta = \widehat \beta - \beta_0$, which measures the accuracy of predicting $x_i'\beta_0$ over the design points $x_1,\ldots,x_n$, $$ \|\delta\|_{2,n} := \sqrt{ {\mathbb{E}_n}[(x_i'\delta)^2] } = \sqrt{\delta '{\mathbb{E}_n}[ x_ix_i'] \delta }.$$

The prediction norm directly depends on the the Gram matrix ${\mathbb{E}_n} [x_ix_i']$. Whenever $p>n$, the empirical Gram matrix ${\mathbb{E}_n}[x_ix_i']$ does not have full rank and in principle is not well-behaved. However, we only need good behavior of certain moduli of continuity of the Gram matrix called sparse eigenvalues. We define the minimal $m$-sparse eigenvalue of a semi-definite matrix $M$ as

and the maximal $m$-sparse eigenvalue as

To assume that $\phi_{{\rm min}}(m)[{\mathbb{E}_n} [x_ix_i']] >0$ requires that all empirical Gram submatrices formed by any $m$ components of $x_i$ are positive definite. To simplify asymptotic statements for Lasso and Post-Lasso, we use the following condition:

Condition SE. There is $\ell_n \to \infty$ such that $$\kappa' \leqslant \phi_{{\rm min}}(\ell_n s)[{\mathbb{E}_n} [x_ix_i']] \leqslant \phi_{{\rm max}}(\ell_n s)[{\mathbb{E}_n} [x_ix_i']] \leqslant \kappa'',$$ where $0< \kappa' < \kappa'' < \infty$ are constants that do not depend on $n$.

The following theorem describes the rate of convergence for feasible Lasso in the Gaussian model under Conditions ASM and SE. We formally define the feasible Lasso estimator $\widehat\beta$ as either the Iterated Lasso with penalty level given by $X$-independent rule ((ref)) or $X$-dependent rule ((ref)) or Square-root Lasso with penalty level given by $X$-dependent rule ((ref)) or $X$-independent rule ((ref)), with the confidence level $1-\gamma$ such that

As mentioned before, $\ell_1$-regularized estimators have an inherent bias towards zero and Post-Lasso was proposed to remove this bias, at least in part. It turns out that we can bound the performance of Post-Lasso as a function of Lasso's rate of convergence and Lasso's model selection ability. For common designs, this bound implies that Post-Lasso performs at least as well as Lasso, and it can be strictly better in some cases. Post-Lasso also has a smaller shrinkage bias than Lasso by construction.

The following theorem applies to any Post-Lasso estimator $\widetilde \beta$ computed using the model $\widehat T = \text{support}(\widehat \beta)$ selected by a Feasible Lasso estimator $\widehat \beta$ defined before Theorem (ref).

Monte Carlo Example

In this section we compare the performance of various estimators relative to the ideal oracle linear regression estimator. The oracle estimator applies ordinary least square to the true model by regressing the outcome on only the control variables with non-zero coefficients. Of course, the oracle estimator is not available outside Monte Carlo experiments.

We considered the following regression model: $$ y = x'\beta_0 + \epsilon, \ \ \beta_0 =(1,1,1/2,1/3,1/4,1/5,0,\ldots,0)', $$ where $x = (1,z')'$ consists of an intercept and covariates $z \sim N(0,\Sigma)$, and the errors $\epsilon$ are independently and identically distributed $\epsilon \sim N(0,\sigma^2)$. The dimension $p$ of the covariates $x$ is $500$, and the dimension $s$ of the true model is $6$. The sample size $n$ is $100$. The regressors are correlated with $\Sigma_{ij} = \rho^{|i-j|}$ and $\rho = .5$. We consider the levels of noise to be $\sigma = 1$ and $\sigma=0.1$. For each repetition we draw new $x$'s and $\epsilon$'s.

We consider infeasible Lasso and Post-Lasso estimators, feasible Lasso and Post-Lasso estimators described in the previous section, all with X-dependent penalty levels, as well as (5-fold) cross-validated (CV) Lasso and Post-Lasso. We summarize results on estimation performance in Table (ref) which records for each estimator $\bar \beta$ the norm of the bias $\|{\mathrm{E}}[ \bar \beta - \beta_0]\|$ and also the empirical risk $\{{\mathrm{E}}[(x_i'( \bar \beta - \beta_0 ))^2]\}^{1/2}$ for recovering the regression function. In this design, infeasible Lasso, Square-root Lasso, and Iterated Lasso exhibit substantial bias toward zero. This bias is somewhat alleviated by choosing the penalty-level via cross-validation, though the remaining bias is still substantial. It is also apparent that, as intuition and theory would suggest, the post-penalized estimators remove a large portion of this shrinkage bias. We see that among the feasible estimators, the best performing methods are the Post-Square-root Lasso and Post-Iterated Lasso. Interestingly, cross-validation also produces a Post-Lasso estimator that performs nearly as well, although the procedure is much more expensive computationally. The Post-Lasso estimators perform better than Lasso estimators primarily due to a much lower shrinkage bias which is beneficial in the design considered.

{

}

Inference on Structural Effects with High-Dimensional Instruments

Methods and Theoretical Results

In this section, we consider the linear instrumental variable (IV) model with many instruments. Consider the Gaussian simultaneous equation model:

Here $y_{1i}$ is the response variable, $y_{2i}$ is the endogenous variable, $w_i$ is a $k_w$-vector of control variables, $z_i = (u_i',w_i')'$ is a vector of instrumental variables (IV), and $(\zeta_i, v_i)$ are disturbances that are independent of $z_i$. The function $f(z_i) = {\mathrm{E}}[y_{2i}|z_i]$, the optimal instrument, is an unknown, potentially complicated function of the elementary instruments $z_i$. The main parameter of interest is the coefficient on $y_{2i}$, whose true value is $\alpha_1$. We treat $\{z_i\}$ as fixed throughout.

Based on these elementary instruments, we create a high-dimensional vector of technical instruments, $x_i = P(z_i)$, with dimension $p$ possibly much larger than the sample size though restricted via conditions stated below. We then estimate the the optimal instrument $ f(z_i)$ by

where $\widehat \beta$ is a feasible Lasso or Post-Lasso estimator as formally defined in the previous section.

Sparse-methods take advantage of approximate sparsity and ensure that many elements of $\widehat\beta$ are zero when $p$ is large. In other words, sparse-methods will select a small subset of the available technical instruments. Let $A_i = (f(z_i), w_i')'$ be the ideal instrument vector, and let

be the estimated instrument vector. Denoting $d_i=(y_{2i},w_i')'$, we form the feasible IV estimator using the estimated instrument vector as

The main regularity condition is recorded as follows.

Condition ASIV. In the linear IV model ((ref))-((ref)) with technical instruments $x_i = P(z_i)$, the following assumptions hold: (i) the parameter values $\sigma_v$, $\sigma_\zeta$ and the eigenvalues of $Q_n={\mathbb{E}_n}[A_iA_i']$ are bounded away from zero and from above uniformly in $n$, (ii) condition ASM holds for ((ref)), namely for each $i =1,...,n$, there exists $\beta_0 \in \Bbb{R}^p$, such that $ f(z_i) = x_i'\beta_0+ r_i, \ \ \|\beta_0\| \leqslant s, \ \ \{ {\mathbb{E}_n}[r_i^2] \}^{1/2} \leqslant K \sigma_{v} \sqrt{s/n}, $ where constant $K$ does not depend on $n$, (iii) condition SE holds for ${\mathbb{E}_n}[x_ix_i']$, and (iv) $s^2\log^2 (p\vee n) = o(n)$.

The main inference result is as follows.

Weak Identification Robust Inference with Very Many Instruments

Consider the simultaneous equation model:

where $y_{1i}$ is the response variable, $y_{2i}$ is the endogenous variable, $w_i$ is a $k_w$-vector of control variables, $z_i = (u_i',w_i')'$ is a vector of instrumental variables (IV), and $\zeta_i$ is a disturbance that is independent of $z_i$. We treat $\{z_i\}$ as fixed throughout.

We would like to use a high-dimensional vector $x_i=P(z_i)$ of technical instruments for inference that is robust to weak identification. We propose a method for inference based on inverting pointwise tests performed using a sup-score statistic defined below. The procedure is similar in spirit to \citeasnoun{anderson:rubin} and \citeasnoun{ss:weakiv} but uses a very different statistics that is well-suited to cases with very many instruments.

In order to formulate the sup-score statistic, we first partial-out the effect of controls $w_i$ on the key variables. For an $n$-vector $\{u_i, i=1,...,n\}$, define $\tilde u_i = u_i - w_i'{\mathbb{E}_n}[w_i w_i']^{-1} {\mathbb{E}_n}[w_i u_i]$, i.e. the residuals left after regressing this vector on $\{w_i, i=1,...,n\}$. Hence $\tilde y_{1i}$, $\tilde y_{2i}$, and $\tilde x_{ij}$ are residuals obtained by partialling out controls. Also, let $\tilde x_i = (\tilde x_{i1},...,\tilde x_{ip})'$. In this formulation, we omit elements of $w_i$ from $\tilde x_{ij}$ since they are eliminated by partialling out. We then normalize without loss of generality

The sup-score statistic for testing the hypothesis $\alpha_1 = a $ takes the form: $$ \Lambda_a = \max_{1 \leqslant j \leqslant p} \frac{|n {\mathbb{E}_n} [(\tilde y_{1i} - \tilde y_{2i}a) \tilde x_{ij}]|}{\sqrt{{\mathbb{E}_n}[ (\tilde y_{1i} - \tilde y_{2i}a)^2 \tilde x^2_{ij} ]}}. $$ If the hypothesis $\alpha_1 = a$ is true, then the critical value for achieving level $\gamma$ is

where $W=[w_1,...,w_n]'$, $X=[x_1,...,x_n]'$, and $g_1,...,g_n$ are i.i.d. $N(0,1)$ variables independent of $W$ and $X$; $\tilde{g}_i$ denotes the residuals left after projecting $\{g_i\}$ on $\{w_i\}$ as defined above. We can approximate the critical value $\Lambda(1- \gamma|W,X)$ by simulation conditional on $X$ and $W$. It is also possible to use a simple asymptotic bound on this critical value of the form

for $c > 1$. The finite-sample $(1- \gamma)$ -- confidence region for $\alpha_1$ is then given by $$ \mathcal{C}:=\{ a \in \Bbb{R}: \Lambda_a \leqslant \Lambda(1- \gamma|W,X)\}, $$ while a large sample $(1- \gamma)$ -- confidence region is given by $ \mathcal{C}':= \{ a \in \Bbb{R}: \Lambda_a \leqslant \Lambda(1- \gamma)\}. $

The main regularity condition is recorded as follows.

Condition HDIV. Suppose the linear IV model ((ref)) holds. Consider the $p$-vector of instruments $x_i = P(z_i)$, $i=1,...,n$, such that $(\log p)/n \to 0$. Suppose further that the following assumptions hold uniformly in $n$: (i) the parameter value $\sigma_\zeta$ is bounded away from zero and from above, (ii) the dimension of $w_i$ is bounded and the eigenvalues of the Gram matrix ${\mathbb{E}_n}[w_i w_i']$ are bounded away from zero, (iii) $\|w_i\| \leqslant K$ and $|\tilde x_{ij}| \leqslant K$ for all $1\leqslant i \leqslant n$ and all $1 \leqslant j \leqslant p$, where $K$ is a constant, independent of $n$.

The main inference result is as follows.

Monte Carlo Example: Instrumental Variable Model

The theoretical results presented in the previous sections suggest that using Lasso to aid in fitting the first-stage regression should result in IV estimators with good estimation and inference properties. In this section, we provide simulation evidence on these properties of IV estimators using iterated Lasso to select instrumental variables for a second-stage estimator. We also considered Square-root Lasso for variable selection. The results were similar to those for iterated Lasso, so we report only the iterated Lasso results.

Our simulations are based on a simple instrumental variables model of the form $$

\ \ \ \left(

\right)\mid x_i \sim N\left(0,\left(

\right)\right) \ {\rm i.i.d.,} $$ where $\alpha=1$ is the parameter of interest, and $x_i = (x_{i1},...,x_{i100})' \sim N(0,\Sigma_X)$ is the instrument vector with $E[x_{ih}^2] = \sigma^2_x$ and $Corr(x_{ih},x_{ij}) = .5^{|j-h|}$. In all simulations, we set $\sigma^2_{\zeta} = 1$ and $\sigma^2_x = 1$. We also use $Corr(\zeta,v$) = .3.

We consider several different settings for the other parameters. We provide simulation results for sample sizes, $n$, of 100 and 500. In one simulation design, we set $\Pi = 0$ and $\sigma^2_v = 1$. In this case, the instruments have no information about the endogenous variable, so $\alpha$ is unidentified. We refer to this as the “No Signal” design. In the remaining cases, we use an “exponential” design for the first stage coefficients, $\Pi$, that sets the coefficient on $x_{ih} = .7^{h-1}$ for $h=1,...,100$ to provide an example of Lasso's performance in settings where the instruments are informative. This model is approximately sparse, since the majority of explanatory power is contained in the first few instruments, and obeys the regularity conditions put forward above. We consider values of $\sigma^2_v$ which are chosen to benchmark three different strengths of instruments. The three values of $\sigma^2_v$ are found as $\sigma^2_v = \frac{n \Pi'\Sigma_Z\Pi}{F^*\Pi'\Pi}$ for $F^*$ of 10, 40, or 160.

For each setting of the simulation parameter values, we report results from several estimation procedures. A simple possibility when presented with $p < n$ instrumental variables is to just estimate the model using 2SLS and all of the available instruments. It is well-known that this will result in poor-finite sample properties unless there are many more observations than instruments; see, for example, \citeasnoun{bekker}. Fuller's fuller estimator (FULL)\footnote{The Fuller estimator requires a user-specified parameter. We set this parameter equal to one which produces a higher-order unbiased estimator. See \citeasnoun{hhk:weakmse} for additional discussion.} is robust to many instruments as long as the presence of many instruments is accounted for when constructing standard errors and $p < n$; see \citeasnoun{bekker} and \citeasnoun{hhn:weakiv} for example. We report results for these estimators in rows labeled 2SLS(All) and FULL(All) respectively.\footnote{All models include an intercept. With $n = 100$, we randomly select 98 instruments to use for 2SLS(All) and FULL(All).} In addition, we report Fuller and IV estimates based on the set of instruments selected by Lasso with two different penalty selection methods. IV-Lasso and FULL-Lasso are respectively 2SLS and Fuller using instruments selected by Lasso with penalty obtained using the iterated method outlined in Appendix A. We use an initial estimate of the noise level obtained using the regression of $y_2$ on the instrument that has the highest simple correlation with $y_2$. IV-Lasso-CV and FULL-Lasso-CV are respectively 2SLS and Fuller using instruments selected by Lasso using 10-fold cross-validation to choose the penalty level. We also report inference results based on the Sup-Score test developed in Section 5.2.

In Table (ref), we report root-mean-squared-error (RMSE), median bias (Med. Bias), rejection frequencies for 5% level tests (rp(.05)), and the number of times the Lasso-based procedures select no instruments ($\|\widehat\Pi\|_0 = 0$). For computing rejection frequencies, we estimate conventional 2SLS standard errors for all 2SLS estimators, and the many instrument robust standard errors of \citeasnoun{hhn:weakiv} for the Fuller estimators. In cases where Lasso selects no instruments, the reported Lasso point estimation properties are based on the feasible procedure that enforces identification by lowering the penalty until one variable is selected. Rejection frequencies in cases where no instruments are selected are based on the feasible procedure that uses conventional IV inference using the selected instruments when this set is non-empty and otherwise uses the Sup-Score test.

The simulation results show that Lasso-based IV estimators is useful in situations with many instruments. As expected, 2SLS(All) does extremely poorly along all dimensions. FULL(All) also performs worse than the Lasso-based estimators in terms of estimator risk (RMSE) in all cases. The Lasso-based procedures do not dominate FULL(All) in terms of median bias, though all of the Lasso-based procedures have smaller median bias than FULL(All) when $n = 100$ and there is some signal in the instruments and are very similar with $n = 500$. In terms of size of 5% level tests, we see that the Sup-Score test uniformly controls size as indicated by the theory. IV-Lasso and FULL-Lasso using the iterated penalty selection method also do a very good job controlling size across all of the simulation settings with a worst-case rejection frequency of .064 (with simulation standard error of .01) and the majority of rejection frequencies below .05. Interestingly, when there is no signal in the instrument, the Lasso-based estimators using penalty selected by CV have substantial size-distortions when $n = 100$ which is due to the CV penalty being small enough that instruments are still selected despite there being no signal. The iterated penalty is such that, at least approximately, only instruments whose coefficients are outside of a $\sqrt{n}$ neighborhood of 0 are selected and thus overselection in cases with little signal is guarded against. Despite the problem with using CV when there is no signal, it is worth noting that the Lasso-based procedures with CV penalty produce tests with approximately correct size in all other parameter settings.

{

}

To further examine the properties of the inference procedures that appear to give small size distortions, we plot the power curves of 5% level tests using the Sup-Score test and IV-Lasso with the iterated and CV penalty choices with $n = 100$ in Figure (ref).\footnote{The power curves in the $n = 500$ case are qualitatively similar.} We see that both the Sup-Score test and IV-Lasso using the iterated procedure augmented with Sup-Score test when no instruments are selected appear to uniformly control size and have some power against alternatives when the model is identified. It is also clear that of these two procedures, the IV-Lasso has substantially more power than the Sup-Score test. The figures also show that IV-Lasso with iterated penalty has almost as much power as IV-Lasso using the CV penalty while avoiding the substantial size distortion and spurious power produced by using CV when there is no signal.

Overall, the simulation results are favorable to the Lasso-based IV methods. The Lasso-based estimators dominate the other estimators considered based on RMSE and have relatively small finite sample biases. The Lasso-based procedures also do a good job in producing tests with size close to the nominal level. There is some evidence that the Fuller-Lasso may do better than 2SLS-Lasso in terms of testing performance though these procedures are very similar in the designs considered. It also seems that tests based on IV-Lasso using the iterated penalty selection rule may perform better than tests based on IV-Lasso using cross-validation to choose the Lasso penalty levels, especially when there is little explanatory power in the instruments.

Inference on Treatment and Structural Effects Conditional on Observables

Methods and Theoretical Results

We consider the following partially linear model,

where $d_i$ is a policy/treament variable whose impact we would like to infer, and $z_i$ represents confounding factors on which we need to condition. This model is of interest in our international growth example discussed in the next section as well as in many empirical studies heckman:metricslabormarkets,imbens:review. The confounding factors affect the policy variable via $m(z_i)$. We assume that $m(z_i)$ and $g(z_i)$ each admit an approximately sparse form and use linear combinations of technical control terms $x_i = P(z_i)$ to approximate them.

There are at least three obvious strategies for inference:

Note that it is most natural not to penalize $\alpha_0$ since the goal is to quantify the impact of $d_i$. (The previous rate results derived in Theorems (ref) and (ref) for the regression function extend to the case where the coefficients on a fixed number of variables are not penalized.) In what follows, we shall refer to options (i), (ii), and (iii) respectively as Lasso, Post-Lasso, and Indirect Post-Lasso.

Regarding inference, “intuition" suggests that if $g$ can be estimated at faster than the $n^{1/4}$ rate then any of (i)-(iii) could be $\sqrt{n}$-consistent and asymptotically normal. It turns out that this “intuition" is often correct for options (ii) and (iii) but is wrong for option (i). Indeed, it is possible to show that under rather strong regularity conditions that

where $\sigma_{\zeta}^2 [\mathbb{E}_n v_i^2]^{-1}$ is the semi-parametric efficiency bound for estimating $\alpha_0$, for $\bar \alpha$ denoting the estimators (ii) or (iii) above. Unfortunately, the distributional result ((ref)) is not very robust to modest violations of regularity conditions and may provide a poor approximation to the finite-sample distributions of the estimators for $\alpha_0$. The reason is that Lasso applied to ((ref)) may miss important terms relating $d_i$ to $z_i$ through $m(z_i)$ and thus suffer from substantial omitted variables bias. On the other hand, Lasso applied only to ((ref)), even if successful in selecting adequate controls for the relationship between $d_i$ and $z_i$, may miss important terms in $g(z_i)$ and thus be highly inefficient. We illustrate this lack of robustness through a simulation experiment reported below.

Instead of using Lasso, Post-Lasso, or Indirect Post-Lasso, we advocate a “double-Post-Lasso” method. To define this estimator, we write the reduced form corresponding to ((ref))-((ref)):

Now we have two equations and hence can apply Lasso methods to each equation to select control terms. That is, we run Lasso regression of $y_{1i}$ on $x_i=P(z_i)$ and Lasso regression of $d_i$ on $x_i=P(z_i)$. Then we can run least squares of $y_{1i}$ on $d_i$ and the union of the controls selected in each equation to estimate and perform inference on $\alpha_0$. By using this procedure we increase the chances for successfully recovering terms that approximate the key control term $m(z_i)$, which results in improved robustness properties. Indeed, the resulting procedure is considerably more robust in computational experiments and requires much weaker regularity conditions than the obvious strategies outlined above.

Now we formally define the double-Post-Lasso estimator. Let $\widehat I_1 = {\rm support }(\widehat \beta_1)$ denote the control terms selected by a feasible Lasso estimator $\widehat \beta_1$ computed using data $(y_i,x_i) = (d_{i}, x_i), i =1,...,n$. Let $\widehat I_2 = {\rm support }(\widehat \beta_2)$ denote the control terms selected by a feasible Lasso estimator $\widehat \beta_2$ computed using data $(y_i,x_i) = (y_{1i}, x_i), i =1,...,n$. The double-Post-Lasso estimator $\check \alpha$ of $\alpha_0$ is defined as the least squares estimator obtained by regressing $y_{1i}$ on $d_i$ and the selected control terms $x_{ij}$ with $j \in \widehat I \supseteq \widehat I_1 \cup \widehat I_2$: $$ (\check \alpha, \check \beta) = \underset{ \alpha \in \Bbb{R}, \beta \in \Bbb{R}^p}{\rm argmin}\{ {\mathbb{E}_n}[(y_{1i} - d_i \alpha - x_i'\beta)^2] \ : \ \beta_j = 0, \forall j \not \in \widehat I \}. $$ The set $\widehat I$ can contain other variables with names $\widehat I_3$ that the analyst may think are important for ensuring robustness. Thus, $\widehat I = \widehat I_1 \cup \widehat I_2 \cup \widehat I_3$; let $\widehat s = |\widehat I|$ and $\widehat s_j = |\widehat I_j|$ for $j =1,2,3$.

Condition ASTE. (i) The data $(y_{1i}, d_i, z_i), i=1,...,n$, obeys model ((ref))-((ref)) for each $n$, and $x_i = P(z_i)$ is a dictionary of transformations of $z_i$. (ii) The parameter values $\sigma^2_v$ and $\sigma^2_{\zeta}$ are bounded from above by $\bar \sigma$ and away from zero, uniformly in $n$, and $|\alpha_0|$ is bounded uniformly in $n$. (iii) Regressor values $x_i, i=1,...,n$, obey the normalization condition ${\mathbb{E}_n}[x^2_{ij}]=1$ for all $j \in \{1,...,p\}$ and sparse eigenvalue condition SE. (iv) There exists $s \geqslant 1$ and $\beta_{m0}$ and $\beta_{g0}$ such that

where $K$ is an absolute constant, independent of $n$, but all other parameter values can depend $n$. (v) $s^2 \log^2 (p\vee n) = o(n)$ and $ \widehat s_3 \lesssim 1\vee \widehat s_1 \vee \widehat s_2$. }

Monte Carlo Example: Partially Linear Models

In this section, we compare the estimation strategies proposed above in the following model:

where the covariates $\tilde x \sim N(0,\Sigma)$, $\Sigma_{kj} = (0.5)^{|j-k|}$, and

with $\sigma_\zeta = \sigma_v =1$, and $\sigma_{\zeta v}=0$. The dimension $p$ of the covariates $x$ is $200$, and the sample size $n$ is $100$. We set $\alpha_0 =1$ and

We set $\lambda$ according to the $X$-dependent rule with $1-\gamma = .95$. For each repetition we draw new $x$'s, $\zeta$'s and $v$'s.

We summarize the inference performance of these methods in Table (ref) which illustrates mean bias, standard deviation, and rejection probabilities of 95% confidence intervals. As we had expected, Lasso and Post-Lasso exhibit a large mean bias which dominates the estimation error and results in poor performance of conventional inference methods. On the other hand, the Indirect Post-Lasso has a small bias relative to estimation error but is substantially more variable than double-Post-Lasso and produces a conservative test, a test with size much smaller than the nominal level. Notably, the double-Post-Lasso provides coverage that is close to the promised $5\%$ level and has the smallest mean bias and standard deviation.

{

}

Empirical Examples.

In this section, we illustrate the performance of sparse methods in two empirical examples. In the first, we revisit the classic \citeasnoun{AK1991}'s instrumental variables estimation of the returns to schooling. In this example, there are many instruments which can potentially be used in forming the IV estimator and there are concerns about the potential biases and inferential problems introduced from using many instruments. Our results show that sparse methods can be effectively used to alleviate these concerns. The second example concerns the use of $\ell_1$-penalized methods to select control variables for growth regressions in which there are many possible country level controls relative to the number of countries. Using Square-root Lasso to select control variables, we find that there is evidence in favor of the hypothesis of convergence.

Angrist and Krueger Example with 1530 instruments

We consider the \citeasnoun{AK1991} model $$

$$ where $y_{1i}$ is the log(wage) of individual $i$, $y_{2i}$ denotes education, $w_i$ denotes a vector of control variables, and $z_i$ denotes a vector of instrumental variables that affect education but do not directly affect the wage. The data were drawn from the 1980 U.S. Census and consist of 329,509 men born between 1930 and 1939. In this example, $w_i$ is a set of 510 variables: a constant, 9 year-of-birth dummies, 50 state-of-birth dummies, and 450 state-of-birth $\times$ year-of-birth interactions. As instruments, we use three quarter-of-birth dummies and interactions of these quarter-of-birth dummies with the set of state-of-birth and year-of-birth controls in $w_i$ giving a total of 1530 potential instruments. \citeasnoun{AK1991} discusses the endogeneity of schooling in the wage equation and provides an argument for the validity of $z_i$ as instruments based on compulsory schooling laws and the shape of the life-cycle earnings profile. We refer the interested reader to \citeasnoun{AK1991} for further details. The coefficient of interest is $\theta_1$, which summarizes the causal impact of education on earnings.

There are two basic options for estimating $\theta_1$ that have been used in the literature: one uses just the three basic quarter-of-birth dummies and the other uses 180 instruments corresponding to the three quarter-of-birth dummies and their interactions with the 9 main effects for year-of-birth and 50 main effects for state-of-birth. It is commonly-held that using the set of 180 instruments results in 2SLS estimates of $\theta_1$ that have a substantial bias, while using just the three quarter-of-birth dummies results in an estimator with smaller bias but a large variance; see, e.g., \citeasnoun{hhn:weakiv}. Another approach uses the 180 instruments and the Fuller estimator fuller (FULL) with an adjustment for the use of many instruments. Of course, using sparse methods for the first-stage estimation offers another option that could be used in place of any of the aforementioned approaches.

{

}

Table 5 presents estimates of the returns to schooling coefficient using 2SLS and FULL\footnote{We set the user-defined choice parameter in the Fuller estimator equal to one which results in a higher-order unbiased estimator.} and different sets of instruments. Given knowledge of the construction of the instruments, the first three rows of the table correspond to the natural groupings of the instruments into the three main quarter of birth effects, the three quarter-of-birth dummies and their interactions with the 9 main effects for year-of-birth and 50 main effects for state-of-birth, and the full set of 1530 potential instruments. The remaining two rows give results based on using Lasso to select instruments with penalty level given by the simple plug-in rule in Section 3 or by 10-fold cross-validation. Using the plug-in rule, Lasso selects only the dummy for being born in the fourth quarter; and with the cross-validated penalty level, Lasso selects 12 instruments which include the dummy for being born in the third quarter, the dummy for being born in the fourth quarter, and 10 interaction terms. The reported estimates are obtained using Post-Lasso.

The results in Table 5 are interesting and quite favorable to the idea of using Lasso to do variable selection for instrumental variables. It is first worth noting that with 180 or 1530 instruments, there are modest differences between the 2SLS and FULL point estimates that theory as well as evidence in \citeasnoun{hhn:weakiv} suggests is likely due to bias induced by overfitting the 2SLS first-stage which may be large relative to precision. In the remaining cases, the 2SLS and FULL estimates are all very close to each other suggesting that this bias is likely not much of a concern. This similarity between the two estimates is reassuring for the Lasso-based estimates as it suggests that Lasso is working as it should in avoiding overfitting of the first-stage and thus keeping bias of the second-stage estimator relatively small.

For comparing standard errors, it is useful to remember that one can regard Lasso as a way to select variables in a situation in which there is no a priori information about which of the set of variables is important; i.e. Lasso does not use the knowledge that the three quarter of birth dummies are the “main” instruments and so is selecting among 1530 a priori “equal” instruments. Given this, it is again reassuring that Lasso with the more conservative plug-in penalty selects the dummy for birth in the fourth quarter which is the variable that most cleanly satisfies \citeasnoun{AK1991}'s argument for the validity of the instrument set. With this instrument, we estimate the returns-to-schooling to be .0862 with an estimated standard error of .0254. The best comparison is FULL with 1530 instruments which also does not use any a priori information about the relevance of the instruments and estimates the returns-to-schooling as .1019 with a much larger standard error of .0422. One can be less conservative than the plug-in penalty by using cross-validation to choose the penalty level. In this case, 12 instruments are chosen producing a Fuller point estimate (standard error) of .0997 (.0139) or 2SLS point estimate (standard error) of .0982 (.0137). These standard errors are smaller than even the standard errors obtained using information about the likely ordering of the instruments given by using 3 or 180 instruments where FULL has standard errors of .0200 and .0143 respectively. That is, Lasso finds just 12 instruments that contain nearly all information in the first stage and, by keeping the number of instruments small, produces a 2SLS estimate that likely has relatively small bias. We believe that these empirical results are reliable. In particular, we note that the first stage $F$ statistic on the selected 12 instruments is approximately $20$; our computational experiments in the previous section employ designs with $F=10$ and $F=40$ to show that this method works well for both estimation and inference purposes.

As a final check, we report the 95% confidence interval obtained from the Sup-Score test of Section 5.2 based on the three natural groupings of 3, 180, and 1530 instruments. This test is robust to weak or non-identification and is simple to implement. For the three different sets of instruments, we obtain intervals that are much wider but roughly in line with the intervals discussed above. We note that our preferred method from the simulation section only makes use of the Sup-Score test when no instruments are selected, does a good job at controlling size in the simulation, and is more powerful than the Sup-Score test when the instruments contain signal about the endogenous variable. Using this procedure would lead us to use the much more precise IV-Lasso results.

Overall, these results demonstrate that Lasso instrument selection is feasible and produces sensible and what appear to be relatively high-quality estimates in this application. The results from the Lasso-based IV estimators are similar to those obtained from other leading approaches to estimation and inference with many-instruments and do not require ex ante information about which are the most relevant instruments. Thus, the Lasso-based IV procedures should provide a valuable complement to existing approaches to estimation and inference in the presence of many instruments.

Growth Example

In this section, we consider variable selection in an international economic growth example. We use the \citeasnoun{BarroLee1994} data consisting of a panel of 138 countries for the period of 1960 to 1985. We consider the national growth rates in GDP per capita as the dependent variable. In our analysis, we consider a model with $p=62$ covariates which allows for a total of $n=90$ complete observations. Our goal here is to provide estimates which shed light on the convergence hypothesis discussed below by selecting controls from among these covariates.\footnote{We can compare our results to those obtained in other standard models in the growth literature such as BarroSala1995,KoenkerMachado1999.}

One of the central issues in the empirical growth literature is the estimation of the effect of an initial (lagged) level of GDP per capita on the growth rates of GDP per capita. In particular, a key prediction from the classical Solow-Swan-Ramsey growth model is the hypothesis of convergence which states that poorer countries should typically grow faster than richer countries and therefore should tend to catch up with the richer countries over time. This hypothesis implies that the effect of a country's initial level of GDP on its growth rate should be negative. As pointed out in Barro and Sala-i-Martin BarroSala1995, this hypothesis is rejected using a simple bivariate regression of growth rates on the initial level of GDP. (In our case, regression yields a statistically insignificant coefficient of $.00132$.) In order to reconcile the data and the theory, the literature has focused on estimating the effect conditional on characteristics of countries. Covariates that describe such characteristics can include variables measuring education and science policies, strength of market institutions, trade openness, savings rates and others; see BarroSala1995. The theory then predicts that the effect of the initial level of GDP on the growth rate should be negative among otherwise similar countries.

Given that the number of covariates we can condition on is comparable to the sample size, covariate selection becomes an important issue in this analysis; see \citeasnoun{OneMillion}, \citeasnoun{TwoMillion}, \citeasnoun{Sala-i-MartinDoppelhoferMiller2004}. In particular, previous findings came under severe criticisms for relying upon ad hoc procedures for covariate selection; see, e.g., \citeasnoun{OneMillion}. Since the number of covariates is high, there is no simple way to resolve the model selection problem using only standard tools. Indeed the number of possible lower-dimensional model is very large, though see \citeasnoun{OneMillion}, \citeasnoun{TwoMillion} and \citeasnoun{Sala-i-MartinDoppelhoferMiller2004} for attempts to search over millions of these models. Here we use $\ell_1$-penalized methods to attempt to resolve this important issue.

We first present results for covariate selection using the different methods discussed in Section (ref): (a) a simple Post-Square-root-Lasso method which uses controls selected from applying the Square-root-Lasso to select controls in the regression of growth rates on log-GDP and other controls, and (b) the Post-double-selection method, which uses the controls selected by Square-root-Lasso in the regression of log-GDP on other controls and in the regression of growth rates on other controls. These were all based on Square-root Lasso to avoid the estimation of $\sigma$. We present the model selection results in Table (ref).

{

}

Square-root Lasso applied to the regression of growth rates on log-GDP and other controls selected only one control, the log of the black market premium which characterizes trade openness. The double selection method selected infant mortality rate, terms of trade shock, and several education variables (female gross enrollment for secondary education, percentage of “no schooling" in the female population, percentage of “higher school attained" in male population, and average schooling years in female population over the age of 25) to forecast log-GDP but no additional controls were selected to forecast growth. We refer the reader to \citeasnoun{BarroLee1994} and \citeasnoun{BarroSala1995} for a complete definition and discussion of each of these variables.

We then proceeded to construct confidence intervals for the coefficient on initial GDP based on each set of selected variables. We also report estimates of the effect of initial GDP in a model which uses the set of controls obtained from the double-selection procedure and additionally includes the log of the black market premium. We expressly allow for such amelioration strategy in our formal construction of the estimator. Table (ref) shows these results. We find that in all these models the linear regression coefficients on the initial level of GDP are negative. In addition, zero is excluded from the 90% confidence interval in each case. These findings support the hypothesis of (conditional) convergence derived from the classical Solow-Swan-Ramsey growth model. The findings also agree with and thus support the previous findings reported in \citeasnoun{BarroSala1995} which relied on ad-hoc reasoning for covariate selection.

{

}

Conclusion

There are many situations in economics where a researcher has access to data with a large number of covariates. In this article, we have presented results for performing analysis of such data by selecting relevant regressors and estimating their coefficients using $\ell_1$-penalization methods. We gave special attention to the instrumental variables model and the partially linear model, both of which are widely used to estimate structural economic effects. Through simulation and empirical examples, we have demonstrated that $\ell_1$ penalization methods may be usefully employed in these models and can complement tools commonly employed by applied researchers.

Of course, there are many avenues for additional research. The use of $\ell_1$-penalization is only one method of performing estimation with high-dimensional data. It will be interesting to consider and understand the behavior of other methods (e.g. \citeasnoun{HHS2008}, \citeasnoun{FanLi2001}, \citeasnoun{zhang:concave}, \citeasnoun{fan:liao}) for estimating structural economic objects. In addition, extending HDS models and methods to other types of economic models beyond those considered in this article will be interesting. An important problem in economics is the analysis of high-dimensional data in which there are many weak signals within the set of variables considered in which case the sparsity assumption may provide a poor approximation. The sup-score test presented in this article offers one approach to dealing with this problem, but further additional research dealing with this issue seems warranted. It would also be interesting to consider efficient use of high-dimensional data in cases in which scores are not independent across observations which is a much-considered case in economics. Overall, we believe the results in this article provide useful tools for applied economists but that there are still substantial and interesting topics in the use of high-dimensional economic data that warrant further investigation.