Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
77,386 characters · 14 sections · 47 citation commands
\titlerunning{HDSM in Econometrics} \authorrunning{Belloni and Chernozhukov}
\abstract{In this chapter we discuss conceptually high dimensional sparse econometric models as well as estimation of these models using $\ell_1$-penalization and post-$\ell_1$-penalization methods. Focusing on linear and nonparametric regression frameworks, we discuss various econometric examples, present basic theoretical results, and illustrate the concepts and methods with Monte Carlo simulations and an empirical application. In the application, we examine and confirm the empirical validity of the Solow-Swan model for international economic growth.}
We consider linear, high dimensional sparse (HDS) regression models in econometrics. The HDS regression model has a large number of regressors $p$, possibly much larger than the sample size $n$, but only a relatively small number $s < n$ of these regressors are important for capturing accurately the main features of the regression function. The latter assumption makes it possible to estimate these models effectively by searching for approximately the right set of the regressors, using $\ell_1$-based penalization methods. In this chapter we will review the basic theoretical properties of these procedures, established in the works of BickelRitovTsybakov2009,CandesTao2007,MY2007,Lounici2008,BC-PostLASSO,Koltchinskii2009,vdGeer,ZhaoYu2006,ZhangHuang2006, among others (see RigolletTsybakov2010,BC-PostLASSO for a detailed literature review). In this section, we review the modeling foundations as well as motivating examples for these procedures, with emphasis on applications in econometrics.
Let us first consider an exact or parametric HDS regression model, namely,
where $y_i$'s are observations of the response variable, $x_i$'s are observations of $p$-dimensional fixed regressors, and $\epsilon_i$'s are i.i.d. normal disturbances, where possibly $p \geqslant n$. The key assumption of the exact model is that the true parameter value $\beta_0$ is sparse, having only $s<n$ non-zero components with support denoted by
Next let us consider an approximate or nonparametric HDS model. To this end, let us introduce the regression model
where $y_i$ is the outcome, $z_i$ is a vector of elementary fixed regressors, $z \mapsto f(z)$ is the true, possibly non-linear, regression function, and $\varepsilon_i$'s are i.i.d. normal disturbances. We can convert this model into an approximate HDS model by writing
where $x_i=P(z_i)$ is a $p$-dimensional regressor formed from the elementary regressors by applying, for example, polynomial or spline transformations, $\beta$ is a conformable parameter vector, whose “true" value $\beta_0$ has only $s<n$ non-zero components with support denoted as in ((ref)), and $r_i:=r(z_i)= f(z_i) - x_i'\beta_0$ is the approximation error. We shall define the true value $\beta_0$ more precisely in the next section. For now, it is important to note only that we assume there exists a value $\beta_0$ having only $s$ non-zero components that sets the approximation error $r_i$ to be small.
Before considering estimation, 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 example, 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$. In particular, we are interested in improving upon the conventional low-dimensional approximations.
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. Since measured education takes on a finite number of years, we can expand the conditional expectation of wage $y_i$ given education $z_i$:
using some dictionary of approximating functions $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 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 $s=4$ or $5$ terms, but 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 exhibit oscillatory behavior near the schooling levels associated with advanced degrees, such as MBA or MD. Low-degree polynomials may not be able to capture this behavior very well, resulting in large approximation errors $\tilde r_i$'s.
Therefore, the question is: With the same number of parameters, can we find a much better approximation? In other words, can we find some higher-order terms in the expansion ((ref)) which will provide a higher-quality approximation? More specifically, can we construct an approximation
for some regressor indices $k_1,\ldots, k_s$ selected from $\{1,\ldots,p\}$, that is accurate and much better than ((ref)), in the sense of having a much smaller approximation error $r_i$?
Obviously the answer to the latter question depends on how complex the behavior of the true regression function ((ref)) is. If the behavior is not complex, then low-dimensional approximation should be accurate. Moreover, it is clear that the second approximation ((ref)) is weakly better than the first ((ref)), and can be much better if there are some important high-order terms in ((ref)) that are completely missed by the first approximation. Indeed, in the context of the earning function example, such important high-order terms could capture abrupt positive changes in earning associated with advanced degrees such as MBA or MD. 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., Angrist, Chernozhukov and Fernandez-Val ACF2006). Treating this data as the population data, we can then compute $f(z_i)=E[y_i|z_i]$ without error. Figure (ref) plots this function. (Of course, such a strategy is not generally available in the empirical work, since the population data are generally not available.) We then construct two sparse approximations and also plot them in Figure (ref): the first is the conventional one, of the form ((ref)), with $P_1, \ldots, P_s$ representing an $(s-1)$-degree polynomial, and 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 two linear splines terms with knots located at 16 and 19 years of schooling (in the case of $s=5$ a third knot is located at 17). In fact, we find the latter approximation automatically using $\ell_1$-penalization methods, 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 requires looking at a very large set of models. We avoided this exhaustive search by using $\ell_1$-penalized least squares (LASSO), which penalizes the size of the model through the sum of absolute values of regression coefficients. Table (ref) quantifies the performance of the different sparse approximations. (Of course, a simple strategy of eye-balling also works in this simple illustrative setting, but clearly does not apply to more general examples with several conditioning variables $z_i$, for example, when we want to condition on education, experience, and age.) \qed
The next two applications are natural examples with large sets of regressors among which we need to select some smaller sets to be used in further estimation and inference. These examples illustrate the potential wide applicability of HDS modeling in econometrics, since many classical and new data sets have naturally multi-dimensional regressors. For example, the American Housing Survey records prices and multi-dimensional features of houses sold, and scanner data-sets record prices and multi-dimensional information on products sold at a store or on the internet. \\
Example 2: Instrument Selection in Angrist and Krueger Data. The second example we consider is an instrumental variables model, as in Angrist and Krueger AK1991$$
$$ where, for person $i$, $y_{i1}$ denotes wage, $y_{i2}$ denotes education, $w_i$ denotes a vector of control variables, and $x_i$ denotes a vector of instrumental variables that affect education but do not directly affect the wage. The instruments $x_i$ come from the quarter-of-birth dummies, and from a very large list, total of $180$, formed by interacting quarter-of-birth dummies with control variables $w_i$. The interest focuses on measuring the coefficient $\theta_1$, which summarizes the causal impact of education on earnings, via instrumental variable estimators.
There are two basic options used in the literature: one uses just the quarter-of-birth dummies, that is, the leading 3 instruments, and another uses all 183 instruments. It is well known that using just 3 instruments results in estimates of the schooling coefficient $\theta_1$ that have a large variance and small bias, while using 183 instruments results in estimates that have a much smaller variance but (potentially) large bias, see, e.g., hhn:weakiv. It turns out that, under some conditions, by using $\ell_1$-based estimation of the first stage, we can construct estimators that also have a nearly efficient variance and at the same time small bias. Indeed, as shown in Table (ref), using the LASSO estimator induced by different penalty levels defined in Section (ref), it is possible to find just 37 instruments that contain nearly all information in the first stage equation. Limiting the number of the instruments from 183 to just 37 reduces the bias of the final instrumental variable estimator. For a further analysis of IV estimates based on LASSO-selected instruments, we refer the reader to BCCH-LASSOIV.
\qed
Example 3: Cross-country Growth Regression. One of the central issues in the empirical growth literature is estimating the effect of an initial (lagged) level of GDP (Gross Domestic Product) 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 and therefore should tend to catch up with the richer countries. Such a hypothesis implies that the effect of the initial level of GDP on the 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 this data set, linear regression yields an insignificant positive coefficient of $0.0013$.) In order to reconcile the data and the theory, the literature has focused on estimating the effect conditional on the pertinent 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 BarroSala1995. The theory then predicts that for countries with similar other characteristics the effect of the initial level of GDP on the growth rate should be negative (BarroSala1995). Thus, we are interested in a specification of the form:
where $y_i$ is the growth rate of GDP over a specified decade in country $i$, $G_i$ is the initial level of GDP at the beginning of the specified period, and the $X_{ij}$'s form a long list of country $i$'s characteristics at the beginning of the specified period. We are interested in testing the hypothesis of convergence, namely that $\alpha_1 <0$.
Given that in standard data-sets, such as Barro and Lee data BarroLee1994, the number of covariates $p$ we can condition on is large, at least relative to the sample size $n$, covariate selection becomes a crucial issue in this analysis (OneMillion, TwoMillion). In particular, previous findings came under severe criticism for relying on ad hoc procedures for covariate selection. In fact, in some cases, all of the previous findings have been questioned (OneMillion). Since the number of covariates is high, there is no simple way to resolve the model selection problem using only classical tools. Indeed the number of possible lower-dimensional models is very large, although OneMillion and TwoMillion attempt to search over several millions of these models. We suggest $\ell_1$-penalization and post-$\ell_1$-penalization methods to address this important issue. In Section (ref), using these methods we estimate the growth model ((ref)) and indeed find rather strong support for the hypothesis of convergence, thus confirming the basic implication of the Solow-Swan model. \qed
Notation. In what follows, all parameter values are indexed by the sample size $n$, but we omit the index whenever this does not cause confusion. In making asymptotic statements, we assume that $n \to \infty$ and $p=p_n \to \infty$, and we also allow for $s=s_n \to \infty$. We use the notation $(a)_+ = \max\{a,0\}$, $a \vee b = \max\{ a, b\}$ and $a \wedge b = \min\{ a , b \}$. The $\ell_2$-norm is denoted by $\|\cdot\|$ and the “$\ell_0$-norm" $\|\cdot\|_0$ denotes the number of non-zero components 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$ the vector in which $\delta_{Tj} = \delta_j$ if $j\in T$, $\delta_{Tj}=0$ if $j \notin T$. We also use standard notation in the empirical process literature, $$\mathbb{E}_n[f] = \mathbb{E}_n[f(w_i)] = \sum_{i=1}^n f(w_i)/n,$$ and we 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)$. Moreover, for two random variables $X, Y$ we say that $X=_dY$ if they have the same probability distribution. We also define the prediction norm associated with the empirical Gram matrix $\mathbb{E}_n[x_ix_i']$ as $$\|\delta\|_{2,n} = \sqrt{\mathbb{E}_n[(x_i'\delta)^2]}.$$
Throughout the rest of the chapter we consider the nonparametric model introduced in the previous section:
where $y_i$ is the outcome, $z_i$ is a vector of fixed regressors, and $\varepsilon_i$'s are i.i.d. disturbances. Define $x_i=P(z_i)$, where $P(z_i)$ is a $p$-vector of transformations of $z_i$, including a constant, and $f_i = f(z_i)$. For a conformable sparse vector $\beta_0$ to be defined below, we can rewrite ((ref)) in an approximately parametric form:
where $r_i := f_i - x_i'\beta_0$, $i=1,\ldots,n,$ are approximation errors. We note that in the parametric case, we may naturally choose $x_i'\beta_0=f_i$ so that $r_i=0$ for all $i=1,\ldots,n$. In the nonparametric case, we shall choose $x_i'\beta_0$ as a sparse parametric model that yields a good approximation to the true regression function $f_i$ in equation ((ref)).
Given ((ref)), our target in estimation will become the parametric function $x_i'\beta_0$. Here we emphasize that the ultimate target in estimation is, of course, $f_i$, while $x_i'\beta_0$ is a convenient intermediate target, introduced so that we can approach the estimation problem as if it were parametric. Indeed, the two targets are equal up to approximation errors $r_i$'s that will be set smaller than estimation errors. Thus, the problem of estimating the parametric target $x_i'\beta_0$ is equivalent to the problem of estimating the non-parametric target $f_i$ modulo approximation errors.
With that in mind, we choose our target or “true" $\beta_0$, with the corresponding cardinality of its support $$s= \| \beta_0\|_0,$$ as any solution to the following ideal risk minimization or oracle problem:
We call this problem the oracle problem for the reasons explained below, and we call $$ T = {\rm support}(\beta_0)$$ the oracle or the “true" model. Note that we necessarily have that $s \leqslant n$.
The oracle problem ((ref)) balances the approximation error $\mathbb{E}_n [(f_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_i - x_i'\beta_0)^2] $$ denote the average square error from approximating values $f_i$ by $x_i'\beta_0$, the quantity $ c^2_s + \sigma^2 s/n$ is the optimal value of ((ref)). Typically, the optimality in ((ref)) would balance the approximation error with the variance term so that for some absolute constant $K\geqslant 0$
so that $ \sqrt{c^2_s + \sigma^2 s/n} \lesssim \sigma \sqrt{s/n}.$ Thus, the quantity $\sigma \sqrt{s/n}$ becomes 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, but we in general do not know $T$, since we do not observe the $f_i$'s to attempt to solve the oracle problem ((ref)). Since $T$ is unknown, we will not be able to achieve the exact oracle rates of convergence, but we can hope to come close to this rate.
We consider the case of fixed design, namely we treat the covariate values $x_1,\ldots, x_n$ as fixed. This includes random sampling as a special case; indeed, in this case $x_1,\ldots, x_n$ represent a realization of this sample on which we condition throughout. Without loss of generality, we normalize the covariates so that
We summarize the setup as the following condition.
\\
Condition ASM. We have data $\{(y_i,z_i), i=1,\ldots,n\}$ that for each $n$ obey the regression model ((ref)), which admits the approximately sparse form ((ref)) induced by ((ref)) with the approximation error satisfying ((ref)). The regressors $x_i=P(z_i)$ are normalized as in ((ref)).
\\
Having introduced the model ((ref)) with the target parameter defined via ((ref)), our task becomes to estimate $\beta_0$. We will focus on deriving rate of convergence results in the {\it prediction norm}, 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}. $$ In what follows $ \delta$ will denote deviations of the estimators from the true parameter value. Thus, e.g., for $ \delta = \widehat \beta - \beta_0$, the quantity $\|\delta\|_{2,n}^2$ denotes the average of the square errors $x_i'\widehat\beta - x_i'\beta_0$ resulting from using the estimate $x_i'\widehat\beta$ instead of $x_i'\beta_0$. Note that once we bound $\widehat \beta - \beta_0$ in the prediction norm, we can also bound the empirical risk of predicting values $f_i$ by $x_i'\widehat \beta$ via the triangle inequality:
In order to discuss estimation consider first the classical ideal AIC/BIC type estimator (Akaike1974,Schwarz1978) that solves the empirical (feasible) analog of the oracle problem: $$ \min_{\beta \in \Bbb{R}^p} \widehat Q (\beta) + \frac{\lambda}{n} \|\beta\|_{0}, $$ where $ \widehat Q (\beta) = \mathbb{E}_n[(y_i - x_i'\beta)^2]$ and $ \| \beta\|_0 = \sum_{j=1}^p 1\{ | \beta_j| > 0 \}$ is the $\ell_0$-norm and $\lambda$ is the penalty level. This estimator has very attractive theoretical properties, but unfortunately it is computationally prohibitive, since the solution to the problem may require solving $\sum_{k \leqslant n} \binom{p}{k}$ least squares problems (generically, the complexity of this problem is NP-hard Natarajan1995,GeJiangYe2010).
One way to overcome the computational difficulty is to consider a convex relaxation of the preceding problem, namely to employ the closest convex penalty -- the $\ell_1$ penalty -- in place of the $\ell_0$ penalty. This construction leads to the so called LASSO estimator:\footnote{The abbreviation LASSO stands for Least Absolute Shrinkage and Selection Operator, c.f. T1996.}
where as before $ \widehat Q (\beta) = \mathbb{E}_n[(y_i - x_i'\beta)^2]$ and $\|\beta\|_{1} = \sum_{j=1}^p | \beta_j|$. The LASSO estimator minimizes a convex function. Therefore, from a computational complexity perspective, ((ref)) is a computationally efficient (i.e. solvable in polynomial time) alternative to AIC/BIC estimator.
In order to describe the choice of $\lambda$, we highlight that the following key quantity determining this choice: $$ S = 2\mathbb{E}_n[x_i\varepsilon_i], $$ which summarizes the noise in the problem. We would like to choose the smaller penalty level so that
where $1-\alpha$ needs to be close to one, and $c$ is a constant such that $c>1$. Following BC-PostLASSO and BickelRitovTsybakov2009, respectively, we consider two choices of $\lambda$ that achieve the above:
where $\alpha \in (0,1)$ and $c>1$ is constant, and $$ \Lambda(1-\alpha|X) := (1-\alpha)-\text{quantile of } n\|S/(2\sigma)\|_{\infty} , $$ $ \text{ conditional on } X=(x_1,\ldots,x_n)'$. Note that $$ \|S/(2\sigma)\|_{\infty} =_d \max_{1 \leqslant j \leqslant p} |\mathbb{E}_n [x_{ij} g_i]|, \text{ where $g_i$'s are i.i.d. } N(0,1), $$ conditional on $X$, so we can compute $\Lambda(1-\alpha|X)$ simply by simulating the latter quantity, given the fixed design matrix $X$. Regarding the choice of $\alpha$ and $c$, asymptotically we require $\alpha \to 0$ as $n \to \infty$ and $c>1$. Non-asymptotically, in our finite-sample experiments, $\alpha = .1$ and $c=1.1$ work quite well. The noise level $\sigma$ is unknown in practice, but we can estimate it consistently using the approach of Section 6. We recommend the $X$-dependent rule over the $X$-independent rule, since the former by construction adapts to the design matrix $X$ and is less conservative than the latter in view of the following relationship that follows from Lemma (ref):
Regularization by the $\ell_1$-norm employed in ((ref)) naturally helps the LASSO estimator to avoid overfitting the data, but it also shrinks the fitted coefficients towards zero, causing a potentially significant bias. In order to remove some of this bias, let us consider the Post-LASSO estimator that applies ordinary least squares regression to the model $\widehat T$ selected by LASSO. Formally, set $$\widehat T = {\rm support}( \widehat \beta ) = \{ j \in \{1,\ldots,p\} \ : \ |\widehat\beta_j| > 0\},$$ and define the Post-LASSO 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 by LASSO. If the model selection works perfectly -- that is, $\widehat T = T$ -- then the Post-LASSO estimator is simply the oracle estimator whose properties are well known. However, perfect model selection might be unlikely for many designs of interest, so we are especially interested in the properties of Post-LASSO in such cases, namely when $\widehat T \neq T$, especially when $T\nsubseteq \widehat T$.
In this section we discuss the intuition behind LASSO and Post-LASSO estimators defined above. We shall rely on a dual interpretation of the LASSO optimization problem to provide some geometrical intuition for the performance of LASSO. Indeed, it can be seen that the LASSO estimator also solves the following optimization program:
for some value of $\gamma \geqslant 0$ (that depends on the penalty level $\lambda$). Thus, the estimator minimizes the $\ell_1$-norm of coefficients subject to maintaining a certain goodness-of-fit; or, geometrically, the LASSO estimator searches for a minimal $\ell_1$-ball -- the diamond-- subject to the diamond having a non-empty intersection with a fixed lower contour set of the least squares criterion function -- the ellipse.
In Figure (ref) we show an illustration for the two-dimensional case with the true parameter value $(\beta_{01}, \beta_{02})$ equal $(1,0)$, so that $T=\text{support}(\beta_0)=\{1\}$ and $s=1$. In the figure we plot the diamonds and ellipses. In the top figure, the ellipse represents a lower contour set of the population criterion function $Q(\beta) = E[(y_i-x_i'\beta)^2]$ in the zero noise case or the infinite sample case. In the bottom figures the ellipse represents a contour set of the sample criterion function $\widehat Q(\beta)=\mathbb{E}_n[(y_i-x_i'\beta)^2]$ in the non-zero noise or the finite sample case. The set of optimal solutions $\widehat \beta$ for LASSO is then given by the intersection of the minimal diamonds with the ellipses. Finally, recall that Post-LASSO is computed as the ordinary least square solution using covariates selected by LASSO. Thus, Post-LASSO estimate $\widetilde \beta$ is given by the center of the ellipse intersected with the linear subspace selected by LASSO.
In the zero-noise case or in population (top figure), LASSO easily recovers the correct sparsity pattern of $\beta_0$. Note that due to the regularization, in spite of the absence of noise, the LASSO estimator has a large bias towards zero. However, in this case Post-LASSO $\widetilde \beta$ removes the bias and recovers $\beta_0$ perfectly.
In the non-zero noise case (middle and bottom figures), the contours of the criterion function and its center move away from the population counterpart. The empirical error in the middle figure moves the center of the ellipse to a non-sparse point. However, LASSO correctly sets $\widehat \beta_2 = 0$ and $\widehat\beta_1\neq 0$ recovering the sparsity pattern of $\beta_0$. Using the selected support, Post-LASSO $\widetilde \beta$ becomes the oracle estimator which drastically improves upon LASSO. In the case of the bottom figure, we have large empirical errors that push the center of the lower contour set further away from the population counterpart. These large empirical errors make the LASSO estimator non-sparse, incorrectly setting $\widehat \beta_2 \neq 0$. Therefore, Post-LASSO uses $\widehat T = \{1,2\}$ and does not use the exact support $T=\{1\}$. Thus, Post-LASSO is not the oracle estimator in this case.
All three figures also illustrate the shrinkage bias towards zero in the LASSO estimator that is introduced by the $\ell_1$-norm penalty. The Post-LASSO estimator is motivated as a solution to remove (or at least alleviate) this shrinkage bias. In cases where LASSO achieves a good sparsity pattern, Post-LASSO can drastically improve upon LASSO.
In both the parametric and non-parametric models described above, whenever $p>n$, the empirical Gram matrix $\mathbb{E}_n[x_ix_i']$ does not have full rank and hence it is not well-behaved. However, we only need good behavior of certain moduli of continuity of the Gram matrix called restricted sparse eigenvalues. We define the minimal restricted sparse eigenvalue
and the maximal restricted sparse eigenvalue as
where $m$ is the upper bound on the number of non-zero components outside the support $T$. To assume that $\kappa(m)>0$ requires that all empirical Gram submatrices formed by any $m$ components of $x_i$ in addition to the components in $T$ are positive definite. It will be convenient to define the following sparse condition number associated with the empirical Gram matrix:
In order to state simplified asymptotic statements, we shall also invoke the following condition. \\
Condition RSE. Sparse eigenvalues of the empirical Gram matrix are well behaved, in the sense that for $m=m_n=s\log n$
This condition holds with high probability for many designs of interest under mild conditions on $s$. For example, as shown in Lemma (ref), when the covariates are Gaussians, the conditions in ((ref)) are true with probability converging to one under the mild assumption that $s\log p = o(n)$. Condition RSE is likely to hold for other regressors with jointly light-tailed distributions, for instance log-concave distribution. As shown in Lemma (ref), the conditions in ((ref)) also hold for general bounded regressors under the assumption that $s (\log^4 n)\log (p\vee n) = o(n)$. Arbitrary bounded regressors often arise in non-parametric models, where regressors $x_i$ are formed as spline, trigonometric, or polynomial transformations $P(z_i)$ of some elementary bounded regressors $z_i$.
For proofs, see BC-PostLASSO; the first lemma builds upon results in ZhangHuang2006 and the second builds upon results in RudelsonVershynin2008.
In this section we discuss the rate of convergence of LASSO in the prediction norm; our exposition follows mainly BickelRitovTsybakov2009.
The key quantity in the analysis is the following quantity called “score": $$ S = S(\beta_0) = 2 \mathbb{E}_n [ x_i \varepsilon_i]. $$ The score is the effective “noise" in the problem. Indeed, defining $ \delta := \widehat \beta - \beta_0, $ note that by the H\"older's inequality
Intuition suggests that we need to majorize the “noise term" $\| S\|_{\infty}$ by the penalty level $\lambda/n$, so that the bound on $\| \delta\|^2_{2,n} $ will follow from a relation between the prediction norm $\| \cdot \|_{2,n}$ and the penalization norm $\|\cdot\|_{1}$ on a suitable set. Specifically, for any $c > 1$, it will follow that if $$\lambda \geqslant cn\|S\|_\infty$$ and $\|\delta\|_{2,n} \geqslant 2c_s$, the vector $\delta$ will also satisfy
where ${\bar c} = (c+1)/(c-1)$. That is, in this case the error in the regularization norm outside the true support does not exceed ${\bar c}$ times the error in the true support. (In the case $\|\delta\|_{2,n} \leqslant 2c_s$ the inequality ((ref)) may not hold, but the bound $\|\delta\|_{2,n} \leqslant 2c_s$ is already good enough.)
Consequently, the analysis of the rate of convergence of LASSO relies on the so-called restricted eigenvalue $\kappa_{\bar c}$, introduced in BickelRitovTsybakov2009, which controls the modulus of continuity between the prediction norm $\| \cdot \|_{2,n}$ and the penalization norm $\|\cdot\|_{1}$ over the set of vectors $\delta \in \Bbb{R}^p$ that satisfy ((ref)):
where $\kappa_{\bar c}$ can depend on $n$. The constant $\kappa_{\bar c}$ is a crucial technical quantity in our analysis and we need to bound it away from zero. In the leading cases that condition RSE holds this will in fact be the case as the sample size grows, namely
Indeed, we can bound $\kappa_{\bar c}$ from below by $$
$$ by Lemma \ref{Lemma:BoundKAPPA} stated and proved in the appendix. Thus, under the condition RSE, as $n$ grows, $\kappa_{\bar c}$ is bounded away from zero since $\kappa(s\log n)$ is bounded away from zero and $\phi(s\log n)$ is bounded from above as in (\ref{Eq:PRIMITIVE}). Several other primitive assumptions can be used to bound $\kappa_{\bar c}$. We refer the reader to \cite{BickelRitovTsybakov2009} for a further detailed discussion of lower bounds on $\kappa_{\bar c}$.
We next state a non-asymptotic performance bound for the LASSO estimator.
The proof of Theorem (ref) is given in the appendix. The theorem also leads to the following useful asymptotic bounds.
The non-asymptotic and asymptotic bounds for the empirical risk immediately follow from the triangle inequality:
Thus, the rate of convergence of $x_i'\widehat \beta$ to $f_i$ coincides with the rate of convergence of the oracle estimator $\sqrt{ c^2_s + \sigma^2 s/n}$ up to a logarithmic factor of $p$. Nonetheless, the performance of LASSO can be considered optimal in the sense that under general conditions the oracle rate is achievable only up to logarithmic factor of $p$ (see Donoho and Johnstone DonohoJohnstone1994 and Rigollet and Tsybakov RigolletTsybakov2010), apart from very exceptional, stringent cases, in which it is possible to perform perfect or near-perfect model selection.
The purpose of this section is, first, to provide bounds (sparsity bounds) on the dimension of the model selected by LASSO, and, second, to describe some special cases where the model selected by LASSO perfectly matches the “true" (oracle) model.
Although perfect model selection can be seen as unlikely in many designs, sparsity of the LASSO estimator has been documented in a variety of designs. Here we describe the sparsity results obtained in BC-PostLASSO. Let us define $$\widehat m := |\widehat T\setminus T| = \|\widehat \beta_{T^c}\|_0,$$ which is the number of unnecessary components or regressors selected by LASSO.
Under Conditions ASM and RSE, for $n$ sufficiently large we have $1/\kappa_{\bar c} \lesssim 1$, $c_s \lesssim \sigma\sqrt{s/n}$, and $\phi(s\log n) \lesssim 1$; and under the conditions of Corollary (ref), $\lambda \geq c \sigma\sqrt{n}$ with probability approaching one. Therefore, we have that $L\lesssim_P 1$ and $$ s\log n > s \phi(s\log n) \cdot 2L, \ \ \text{that is,} \ \ s\log n \in \mathcal{M} $$ with probability approaching one as $n$ grows. Therefore, under these conditions we have $$ \min_{m \in \mathcal{M}}\phi(m\wedge n)\lesssim_P 1.$$
Thus, using a penalty level that satisfies ((ref)) LASSO's sparsity is asymptotically of the same order as the oracle sparsity, namely
We note here that Theorem (ref) is particularly helpful in designs in which $ \min_{m \in \mathcal{M}}\phi(m)$ $\ll$ $\phi(n) $. This allows Theorem (ref) to sharpen the sparsity bound of the form $\widehat s \lesssim_P s \phi(n)$ considered in BickelRitovTsybakov2009 and MY2007. The bound above is comparable to the bounds in ZhangHuang2006 in terms of order of magnitude, but Theorem (ref) requires a smaller penalty level $\lambda$ which also does not depend on the unknown sparse eigenvalues as in ZhangHuang2006.
The purpose of this section is to describe very special cases where perfect model selection is possible. Most results in the literature for model selection have been developed for the parametric case only (MY2007,Lounici2008). Below we provide some results for the nonparametric models, which cover the parametric models as a special case.
Thus, we see from parts (1) and (2) that perfect model selection is possible under strong assumptions on the coefficients' separation away from zero. We also see from part (3) that the strong separation of coefficients can be considerably weakened in exchange for a strong assumption on the maximal pairwise correlation of regressors. These results generalize to the nonparametric case the results of Lounici2008 and MY2007 for the parametric case in which $c_s = 0$.
Finally, the following result on perfect model selection also requires strong assumptions on separation of coefficients and the empirical Gram matrix. Recall that for a scalar $v$, $ {\rm sign}(v)=v/|v|$ if $|v|>0$, and $0$ otherwise. If $v$ is a vector, we apply the definition componentwise. Also, given a vector $x \in {\Bbb{R}}^p$ and a set $T \subset \{1,...,p\}$, let us denote $x_{i}[T]:=\{x_{ij}, j \in T\}$.
The result follows immediately from the first order optimality conditions, see Wainright2006. ZhaoYu2006 and CandesPlan2009 provides further primitive sufficient conditions for perfect model selection for the parametric case in which $u_i = \varepsilon_i$. The conditions above might typically require a slightly larger choice of $\lambda$ than ((ref)) and larger separation from zero of the minimal non-zero coefficient $\min_{j\in T}|\beta_{0j}|$.
Next we study the rate of convergence of the Post-LASSO estimator. Recall that for $\widehat T = \text{ support } (\widehat \beta)$, the Post-LASSO estimator solves $$ \widetilde \beta \in \arg\min_{\beta \in \mathbb{R}^p} \ \widehat Q(\beta) : \beta_j = 0 \text{ for each } j \in \widehat T^c. $$ It is clear that if the model selection works perfectly (as it will under some rather stringent conditions discussed in Section (ref)), that is, $T = \widehat T$, then this estimator is simply the oracle least squares estimator whose properties are well known. However, if the model selection does not work perfectly, that is, $T \not = \widehat T$, the resulting performance of the estimator faces two different perils: First, in the case where LASSO selects a model $\widehat T$ that does not fully include the true model $T$, we have a specification error in the second step. Second, if LASSO includes additional regressors outside $T$, these regressors were not chosen at random and are likely to be spuriously correlated with the disturbances, so we have a data-snooping bias in the second step.
It turns out that despite of the possible poor selection of the model, and the aforementioned perils this causes, the Post-LASSO estimator still performs well theoretically, as shown in BC-PostLASSO. Here we provide a proof similar to BCCH-LASSOIV which is easier generalize to non-Gaussian cases.
This theorem provides a performance bound for Post-LASSO as a function of LASSO's sparsity characterized by $\widehat m$, 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, but it can be strictly better in some cases, and has a smaller shrinkage bias by construction.
It is also worth repeating here that finite-sample and asymptotic bounds in other norms of interest immediately follow by the triangle inequality and by definition of $\kappa(\widehat m)$:
The corollary above shows that Post-LASSO achieves the same near-oracle rate as LASSO. Notably, this occurs despite the fact that LASSO may in general fail to correctly select the oracle model $T$ as a subset, that is $T \not \subseteq \widehat T$. The intuition for this result is that any components of $T$ that LASSO misses cannot be very important. This corollary also shows that in some special cases Post-LASSO strictly improves upon LASSO's rate. Finally, note that Corollary (ref) follows by observing that under the stated conditions,
Our specification of penalty levels ((ref)) and ((ref)) require the practitioner to know the noise level $\sigma$ of the disturbances or at least estimate it. The purpose of this section is to propose the following method for estimating $\sigma$. First, we use a conservative estimate $\widehat \sigma^{0} = \sqrt{\text{Var}_n[y_i]}:= \sqrt{\mathbb{E}_n\left[(y_i - \bar y)^2\right]}$, where $\bar y = \mathbb{E}_n[y_{i}]$, in place of $\sigma^2$ to obtain the initial LASSO and Post-LASSO estimates, $\widehat \beta$ and $\widetilde \beta$. The estimate $\widehat \sigma^{0}$ is conservative since $\widehat \sigma^{0} = \sigma^{0} + o_P(1)$ where $\sigma^0 = \sqrt{\text{Var}[y_i]} \geqslant \sigma$, since $x_i$ contains a constant by assumption. Second, we define the refined estimate $\widehat \sigma$ as $$\widehat \sigma= \sqrt{\widehat Q(\widehat \beta)}$$ in the case of LASSO and $$\widehat \sigma = \sqrt{\frac{n}{n-\widehat s} \cdot \widehat Q(\widetilde \beta)} $$ in the case of Post-LASSO. In the latter case we employ the standard degree-of-freedom correction with $\widehat s = \|\widetilde \beta\|_0 = |\widehat T|$, and in the former case we need no additional corrections, since the LASSO estimate is already sufficiently regularized. Third, we use the refined estimate $\widehat \sigma^2$ to obtain the refined LASSO and Post-LASSO estimates $\widehat \beta$ and $\widetilde \beta$. We can stop here or further iterate on the last two steps.
Thus, the algorithm for estimating $\sigma$ using LASSO is as follows:
And the algorithm for estimating $\sigma$ using Post-LASSO is as follows:
We can also use $\lambda= 2c \widehat \sigma^k \sqrt{n} \Phi^{-1}( 1- \alpha/2p) $ in place of $X$-dependent penalty. We note that using LASSO to estimate $\sigma$ it follows that the sequence $\widehat \sigma^k$, $k\geqslant 2$, is monotone, while using Post-LASSO the estimates $\widehat \sigma^k$, $k\geqslant 1$, can only assume a finite number of different values.
The following theorem shows that these algorithms produce consistent estimates of the noise level, and that the LASSO and Post-LASSO estimators based on the resulting data-driven penalty continue to obey the asymptotic bounds we have derived previously.
In this section we compare the performance of LASSO, Post-LASSO, and the ideal oracle linear regression estimators. The oracle estimator applies ordinary least square to the true model. (Such an estimator is not available outside Monte Carlo experiments.)
We begin by considering the following regression model: $$ y = x'\beta_0 + \varepsilon, \ \ \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 $\varepsilon$ are independently and identically distributed $\varepsilon \sim N(0,\sigma^2)$. The dimension $p$ of the covariates $x$ is $500$, the dimension $s$ of the true model is $6$, and the sample size $n$ is $100$. We set $\lambda$ according to the $X$-dependent rule with $1-\alpha = 90\%$. The regressors are correlated with $\Sigma_{ij} = \rho^{|i-j|}$ and $\rho = 0.5$. We consider two levels of noise: Design 1 with $\sigma^2 = 1$ (higher level) and Design 2 with $\sigma^2=0.1$ (lower level). For each repetition we draw new vectors $x_i$'s and errors $\varepsilon_i$'s.
We summarize the model selection performance of LASSO in Figures (ref) and (ref). In the left panels of the figures, we plot the frequencies of the dimensions of the selected model; in the right panels we plot the frequencies of selecting the correct regressors. From the left panels we see that the frequency of selecting a much larger model than the true model is very small in both designs. In the design with a larger noise, as the right panel of Figure (ref) shows, LASSO frequently fails to select the entire true model, missing the regressors with small coefficients. However, it almost always includes the most important three regressors with the largest coefficients. Notably, despite this partial failure of the model selection Post-LASSO still performs well, as we report below. On the other hand, we see from the right panel of Figure (ref) that in the design with a lower noise level LASSO rarely misses any component of the true support. These results confirm the theoretical results that when the non-zero coefficients are well-separated from zero, the penalized estimator should select a model that includes the true model as a subset. Moreover, these results also confirm the theoretical result of Theorem (ref), namely, that the dimension of the selected model should be of the same stochastic order as the dimension of the true model. In summary, the model selection performance of the penalized estimator agrees very well with the theoretical results.
We summarize the results on the performance of estimators in Table (ref), which records for each estimator $\check \beta$ the mean $\ell_0$-norm $E[ \|\check \beta\|_0 ]$, the norm of the bias $\|E \check \beta - \beta_0\|$ and also the prediction error $E[\mathbb{E}_n[ |x_i'( \check \beta - \beta_0 )|^2]^{1/2}]$ for recovering the regression function. As expected, LASSO has a substantial bias. We see that Post-LASSO drastically improves upon the LASSO, particularly in terms of reducing the bias, which also results in a much lower overall prediction error. Notably, despite that under the higher noise level LASSO frequently fails to recover the true model, the Post-LASSO estimator still performs well. This is because the penalized estimator always manages to select the most important regressors. We also see that the prediction error of the Post-LASSO is within a factor $\sqrt{\log p}$ of the prediction error of the oracle estimator, as we would expect from our theoretical results. Under the lower noise level, Post-LASSO performs almost identically to the ideal oracle estimator. We would expect this since in this case LASSO selects the model especially well making Post-LASSO nearly the oracle.
The results above used the true value of $\sigma$ in the choice of $\lambda$. Next we illustrate how $\sigma$ can be estimated in practice. We follow the iterative procedure described in the previous section. In our experiments the tolerance was $10^{-8}$ times the current estimate for $\sigma$, which is typically achieved in less than 15 iterations.
We assess the performance of the iterative procedure under the design with the larger noise, $\sigma^2 = 1$ (similar results hold for $\sigma^2 = 0.1$). The histograms in Figure (ref) show that the model selection properties are very similar to the model selection when $\sigma$ is known. Figure (ref) displays the distribution of the estimator $\widehat \sigma$ of $\sigma$ based on (iterative) Post-LASSO, (iterative) LASSO, and the initial estimator $\widehat\sigma^0=\sqrt{{\rm Var}_n[y_i]}$. As we expected, estimator $\widehat \sigma$ based on LASSO produces estimates that are somewhat higher than the true value. In contrast, the estimator $\widehat \sigma$ based on Post-LASSO seems to perform very well in our experiments, giving estimates $\widehat \sigma$ that bunch closely near the true value $\sigma$.
In this section we apply LASSO and Post-LASSO to an international economic growth example. We use the Barro and Lee 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 a dependent variable $y$ for the periods 1965-75 and 1975-85.\footnote{The growth rate in GDP over a period from $t_1$ to $t_2$ is commonly defined as $\log (GDP_{t_2}/GDP_{t_1})$.} In our analysis, we will consider a model with $p=62$ covariates, which allows for a total of $n=90$ complete observations. Our goal here is to select a subset of these covariates and briefly compare the resulting models to the standard models used in the empirical growth literature (Barro and Sala-i-Martin BarroSala1995).
Let us now turn to our empirical results. We performed covariate selection using LASSO, where we used our data-driven choice of penalty level $\lambda$ in two ways. First we used an upper bound on $\sigma$ being $\widehat \sigma^0$ and decreased the penalty to estimate different models with $\lambda$, $\lambda/2$, $\lambda/3$, $\lambda/4$, and $\lambda/5$. Second, we applied the iterative procedure described in the previous section to define $\lambda^{it}$ (which is computed based on $\widehat\sigma^{it}$ obtained using the iterative Post-LASSO procedure).
The initial choice of the first approach led us to select no covariates, which is consistent with over-regularization since an upper bound for $\sigma$ was used. We then proceeded to slowly decrease the penalty level in order to allow for some covariates to be selected. We present the model selection results in Table (ref). With the first relaxation of the choice of $\lambda$, we select the black market exchange rate premium (characterizing trade openness) and a measure of political instability. With a second relaxation of the choice of $\lambda$ we select an additional set of variables reported in the table. The iterative approach led to a model with only the black market exchange premium. We refer the reader to BarroLee1994 and BarroSala1995 for a complete definition and discussion of each of these variables.
We then proceeded to apply ordinary linear regression to the selected models and we also report the standard confidence intervals for these estimates. Table (ref) shows these results. We find that in all models with additional selected covariates, the linear regression coefficients on the initial level of GDP is always negative and the standard confidence intervals do not include zero. We believe that these empirical findings firmly support the hypothesis of (conditional) convergence derived from the classical Solow-Swan-Ramsey growth model.\footnote{The inferential method used here is actually valid under certain conditions, despite the fact that the model has been selected; this is demonstrated in a work in progress.} Finally, our findings also agree with and thus support the previous findings reported in Barro and Sala-i-Martin BarroSala1995, which relied on ad-hoc reasoning for covariate selection.
{
}