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.
161,274 characters · 56 sections · 83 citation commands
lassopack: Model selection and prediction with regularized regression in Stata
\graphicspath{{fig/}}
\makeatletter \makeatother
\inserttype[st0001]{article} { Achim Ahrens\\The Economic and Social Research Institute\\Dublin, Ireland\\[email removed] \and Christian B. Hansen\\ University of Chicago\\[email removed] \and Mark E. Schaffer\\ Heriot-Watt University\\Edinburgh, United Kingdom \\ [email removed] }
\thispagestyle{empty}
Machine learning is attracting increasing attention across a wide range of scientific disciplines. Recent surveys explore how machine learning methods can be utilized in economics and applied econometrics Varian2014,Mullainathan2017,Athey2018,Kleinberg2018. At the same time, Stata offers to date only a limited set of machine learning tools. lassopack\xspace is an attempt to fill this gap by providing easy-to-use and flexible methods for regularized regression in Stata.\footnote{This article refers to version 1.2 of lassopack\xspace released on the 15th of January, 2019. For additional information and data files, see \url{https://statalasso.github.io/}.}
While regularized linear regression is only one of many methods in the toolbox of machine learning, it has some properties that make it attractive for empirical research. To begin with, it is a straightforward extension of linear regression. Just like ordinary least squares (OLS), regularized linear regression minimizes the sum of squared deviations between observed and model predicted values, but imposes a regularization penalty aimed at limiting model complexity. The most popular regularized regression method is the lasso---which this package is named after---introduced by Frank1993 and Tibshirani1996, which penalizes the absolute size of coefficient estimates.
The primary purpose of regularized regression, like supervised machine learning methods more generally, is prediction. Regularized regression typically does not produce estimates that can be interpreted as causal and statistical inference on these coefficients is complicated.\footnote{This is an active area of research, see for example Buhlmann2013,Meinshausen2009b,Weilenmann2017,Wasserman2009,Lockhart2014.} While regularized regression may select the true model as the sample size increases, this is generally only the case under strong assumptions. However, regularized regression can aid causal inference without relying on the strong assumptions required for perfect model selection. The post-double-selection methodology of Belloni2014b and the post-regularization approach of Chernozhukov2015 can be used to select appropriate control variables from a large set of putative confounding factors and, thereby, improve robustness of estimation of the parameters of interest. Likewise, the first stage of two-step least-squares is a prediction problem and lasso or ridge can be applied to obtain optimal instruments Belloni2012,Carrasco2012,Hansen2014. These methods are implemented in our sister package pdslasso\xspace pdslasso2018, which builds on the algorithms developed in lassopack\xspace.
The strength of regularized regression as a prediction technique stems from the bias-variance trade-off. The prediction error can be decomposed into the unknown error variance reflecting the overall noise level (which is irreducible), the squared estimation bias and the variance of the predictor. The variance of the estimated predictor is increasing in the model complexity, whereas the bias tends to decrease with model complexity. By reducing model complexity and inducing a shrinkage bias, regularized regression methods tend to outperform OLS in terms of out-of-sample prediction performance. In doing so, regularized regression addresses the common problem of overfitting: high in-sample fit (high $R^2$), but poor prediction performance on unseen data.
Another advantage is that the regularization methods of lassopack\xspace---with the exception of ridge regression---are able to produce sparse solutions and, thus, can serve as model selection techniques. Especially when faced with a large number of putative predictors, model selection is challenging. Iterative testing procedures, such as the general-to-specific approach, typically induce pre-testing biases and hypothesis tests often lead to many false positives. At the same time, high-dimensional problems where the number of predictors is large relative to the sample size are a common phenomenon, especially when the true model is treated as unknown. Regularized regression is well-suited for high-dimensional data. The $\ell_1$-penalization can set some coefficients to exactly zero, thereby excluding predictors from the model. The bet on sparsity principle allows for identification even when the number of predictors exceeds the sample size under the assumption that the true model is sparse or can be approximated by a sparse parameter vector.\footnote{Hastie2009 summarize the bet on sparsity principle as follows: `Use a procedure that does well in sparse problems, since no procedure does well in dense problems.'}
Regularized regression methods rely on tuning parameters that control the degree and type of penalization. lassopack\xspace offers three approaches to select these tuning parameters. The classical approach is to select tuning parameters using cross-validation in order to optimize out-of-sample prediction performance. Cross-validation methods are universally applicable and generally perform well for prediction tasks, but are computationally expensive. A second approach relies on information criteria such as the Akaike information criterion Zou2007,Zhang2010. Information criteria are easy to calculate and have attractive theoretical properties, but are less robust to violations of the independence and homoskedasticity assumptions Arlot2010. Rigorous penalization for the lasso and square-root lasso provides a third option. The approach is valid in the presence of heteroskedastic, non-Gaussian and cluster-dependent errors Belloni2012,Belloni2014c,Kozbur2015. The rigorous approach places a high priority on controlling overfitting, thus often producing parsimonious models. This strong focus on containing overfitting is of practical and theoretical benefit for selecting control variables or instruments in a structural model, but also implies that the approach may be outperformed by cross-validation techniques for pure prediction tasks. Which approach is most appropriate depends on the type of data at hand and the purpose of the analysis. To provide guidance for applied reseachers, we discuss the theoretical foundation of all three approaches, and present Monte Carlo results that assess their relative performance.
The article proceeds as follows. In Section (ref), we present the estimation methods implemented in lassopack\xspace. Section (ref)-(ref) discuss the aforementioned approaches for selecting the tuning parameters: information criteria in Section (ref), cross-validation in Section (ref) and rigorous penalization in Section (ref). The three commands, which correspond to the three penalization approaches, are presented in Section (ref), followed by demonstrations in Section (ref). Section (ref) presents Monte Carlo results. Further technical notes are in Section (ref).
\paragraph{Notation.} We briefly clarify the notation used in this article. Suppose $\bm{a}$ is a vector of dimension $m$ with typical element $a_j$ for $j=1,\ldots,m$. The $\ell_1$-norm is defined as $\left\lVert\bm{a}\right\rVert_1 = \sum_{j=1}^m |a_j|$, and the $\ell_2$-norm is $\left\lVert\bm{a}\right\rVert_2 = \sqrt{\sum_{j=1}^m |a_j|^2}$. The `$\ell_0$-norm' of $\bm{a}$ is denoted by $\left\lVerta\right\rVert_0$ and is equal to the number of non-zero elements in $\bm{a}$. $\mathbbm{1}\{.\}$ denotes the indicator function. We use the notation $b\vee c$ to denote the maximum value of $b$ and $c$, i.e., $\max(b,c)$.
This section introduces the regularized regression methods implemented in lassopack\xspace. We consider the high-dimensional linear model \[y_i=\bm{x}_i'\bm{\beta}+\varepsilon_i,\qquad i=1,\ldots, n,\] where the number of predictors, $p$, may be large and even exceed the sample size, $n$. The regularization methods introduced in this section can accommodate large-$p$ models under the assumption of sparsity: out of the $p$ predictors only a subset of $s\ll n$ are included in the true model where $s$ is the sparsity index \[s := \sum_{j=1}^p \mathbbm{1}\{\beta_j\neq 0\}=\left\lVert\bm\beta\right\rVert_0.\] We refer to this assumption as exact sparsity. It is more restrictive than required, but we use it here for illustrative purposes. We will later relax the assumption to allow for non-zero, but `small', $\beta_j$ coefficients. We also define the active set $\Omega=\{j \in \{1,\ldots,p \}: \beta_j\neq 0\}$, which is the set of non-zero coefficients. In general, $p$, $s$, $\Omega$ and $\bm{\beta}$ may depend on $n$ but we suppress the $n$-subscript for notational convenience.
We adopt the following convention throughout the article: unless otherwise noted, all variables have been mean-centered such that $\sum_i y_{i}=0$ and $\sum_ix_{ij}=0$, and all variables are measured in their natural units, i.e., they have not been pre-standardized to have unit variance. By assuming the data have already been mean-centered we simplify the notation and exposition. Leaving the data in natural units, on the other hand, allows us to discuss standardization in the context of penalization.
Penalized regression methods rely on tuning parameters that control the degree and type of penalization. The estimation methods implemented in lassopack\xspace, which we will introduce in the following sub-section, use two tuning parameters: $\lambda$ controls the general degree of penalization and $\alpha$ determines the relative contribution of $\ell_1$ vs. $\ell_2$ penalization. The three approaches offered by lassopack\xspace for selecting $\lambda$ and $\alpha$ are introduced in (ref).
The lasso takes a special position, as it provides the basis for the rigorous penalization approach (see Section (ref)) and has inspired other methods such as elastic net and square-root lasso, which are introduced later in this section. The lasso minimizes the mean squared error subject to a penalty on the absolute size of coefficient estimates:
The tuning parameter $\lambda$ controls the overall penalty level and $\psi_j$ are predictor-specific penalty loadings.
Tibshirani1996 motivates the lasso with two major advantages over OLS. First, due to the nature of the $\ell_1$-penalty, the lasso sets some of the coefficient estimates exactly to zero and, in doing so, removes some predictors from the model. Thus, the lasso serves as a model selection technique and facilitates model interpretation. Secondly, lasso can outperform least squares in terms of prediction accuracy due to the bias-variance trade-off.
The lasso coefficient path, which constitutes the trajectory of coefficient estimates as a function of $\lambda$, is piecewise linear with changes in slope where variables enter or leave the active set. The change points are referred to as knots. $\lambda=0$ yields the OLS solution and $\lambda\rightarrow\infty$ yields an empty model, where all coefficients are zero.
The lasso, unlike OLS, is not invariant to linear transformations, which is why scaling matters. If the predictors are not of equal variance, the most common approach is to pre-standardize the data such that $\frac{1}{n}\sum_i x^2_{ij}=1$ and set $\psi_j=1$ for $j=1,\ldots,p$. Alternatively, we can set the penalty loadings to $\hat\psi_j=(\frac{1}{n}\sum_ix_{ij}^2)^{-1/2}.$ The two methods yield identical results in theory.
Ridge regression Tikhonov1963,Hoerl1970 replaces the $\ell_1$-penalty of the lasso with a $\ell_2$-penalty, thus minimizing
The interpretation and choice of the penalty loadings $\psi_j$ is the same as above. As in the case of the lasso, we need to account for uneven variance, either through pre-estimation standardization or by appropriately choosing the penalty loadings $\psi_j$.
In contrast to estimators relying on $\ell_1$-penalization, the ridge does not perform variable selection. At the same time, it also does not rely on the assumption of sparsity. This makes the ridge attractive in the presence of dense signals, i.e., when the assumption of sparsity does not seem plausible. Dense high-dimensional problems are more challenging than sparse problems: for example, Dicker2016 shows that, if $p/n\rightarrow\infty$, it is not possible to outperform a trivial estimator that only includes the constant. If $p,n\rightarrow$ jointly, but $p/n$ converges to a finite constant, the ridge has desirable properties in dense models and tends to perform better than sparsity-based methods Hsu2014,Dicker2016,Dobriban2018.
Ridge regression is closely linked to principal component regression. Both methods are popular in the context of multicollinearity due to their low variance relative to OLS. Principal components regression applies OLS to a subset of components derived from principal component analysis; thereby discarding a specified number of components with low variance. The rationale for removing low-variance components is that the predictive power of each component tends to increase with the variance. The ridge can be interpreted as projecting the response against principal components while imposing a higher penalty on components exhibiting low variance. Hence, the ridge follows a similar principle; but, rather than discarding low-variance components, it applies a more severe shrinkage Hastie2009.
A comparison of lasso and ridge regression provides further insights into the nature of $\ell_1$ and $\ell_2$ penalization. For this purpose, it is helpful to write lasso and ridge in constrained form as
and to examine the shapes of the constraint sets. The above optimization problems use the tuning parameter $\tau$ instead of $\lambda$. Note that there exists a data-dependent relationship between $\lambda$ and $\tau$.
Figure (ref) illustrates the geometry underpinning lasso and ridge regression for the case of $p=2$ and $\psi_1=\psi_2=1$ (i.e., unity penalty loadings). The red elliptical lines represent residual sum of squares contours and the blue lines indicate the lasso and ridge constraints. The lasso constraint set, given by $|\beta_1| + |\beta_2| \leq \tau$, is diamond-shaped with vertices along the axes from which it immediately follows that the lasso solution may set coefficients exactly to 0. In contrast, the ridge constraint set, $\beta_1^2 + \beta_2^2 \leq \tau$, is circular and will thus (effectively) never produce a solution with any coefficient set to 0. Finally, $\hat\beta_0$ in the figure denotes the solution without penalization, which corresponds to OLS. The lasso solution at the corner of the diamond implies that, in this example, one of the coefficients is set to zero, whereas ridge and OLS produce non-zero estimates for both coefficients.
While there exists no closed form solution for the lasso, the ridge solution can be expressed as \[ \bm{\hat\beta}_{\textrm{ridge}} = (\bm{X}'\bm{X}+\lambda\bm{\Psi}'\bm{\Psi})^{-1}\bm{X}'\bm{y}. \] Here $\bm{X}$ is the $n \times p$ matrix of predictors with typical element $x_{ij}$, $\bm{y}$ is the response vector and $\bm{\Psi}=\textrm{diag}(\psi_1,\ldots,\psi_p)$ is the diagonal matrix of penalty loadings. The ridge regularizes the regressor matrix by adding positive constants to the diagonal of $\bm{X}'\bm{X}$. The ridge solution is thus well-defined generally as long as all the $\psi_j$ and $\lambda$ are sufficiently large even if $\bm{X}'\bm{X}$ is rank-deficient.
The elastic net of ZouHastie2005 combines some of the strengths of lasso and ridge regression. It applies a mix of $\ell_1$ (lasso-type) and $\ell_2$ (ridge-type) penalization:
The additional parameter $\alpha$ determines the relative to contribution of $\ell_1$ vs. $\ell_2$ penalization. In the presence of groups of correlated regressors, the lasso typically selects only one variable from each group, whereas the ridge tends to produce similar coefficient estimates for groups of correlated variables. On the other hand, the ridge does not yield sparse solutions impeding model interpretation. The elastic net is able to produce sparse solutions for some $\alpha$ greater than zero, and retains or drops correlated variables jointly.
The irrepresentable condition (IRC) is shown to be sufficient and (almost) necessary for the lasso to be model selection consistent Zhao2006,Meinshausen2006. However, the IRC imposes strict constraints on the degree of correlation between predictors in the true model and predictors outside of the model. Motivated by this non-trivial condition for the lasso to be variable selection consistent, Zou2006 proposed the adaptive lasso. The adaptive lasso uses penalty loadings of $\psi_j=1/|\hat\beta_{0,j}|^\theta$ where $\hat\beta_{0,j}$ is an initial estimator. The adaptive lasso is variable-selection consistent for fixed $p$ under weaker assumptions than the standard lasso. If $p<n$, OLS can be used as the initial estimator. Huang2008 prove variable selection consistency for large $p$ and suggest using univariate OLS if $p>n$. The idea of adaptive penalty loadings can also be applied to elastic net and ridge regression Zou2009.
The square-root lasso,
is a modification of the lasso that minimizes the root mean squared error, while also imposing an $\ell_1$-penalty. The main advantage of the square-root lasso over the standard lasso becomes apparent if theoretically grounded, data-driven penalization is used. Specifically, the score vector, and thus the optimal penalty level, is independent of the unknown error variance under homoskedasticity as shown by Belloni2011a, resulting in a simpler procedure for choosing $\lambda$ (see Section (ref)).
Penalized regression methods induce an attenuation bias that can be alleviated by post-estimation OLS, which applies OLS to the variables selected by the first-stage variable selection method, i.e.,
where $\tilde\beta_j$ is a sparse first-step estimator such as the lasso. Thus, post-estimation OLS treats the first-step estimator as a genuine model selection technique. For the case of the lasso, Belloni2013 have shown that the post-estimation OLS, also referred to as post-lasso, performs at least as well as the lasso under mild additional assumptions if theory-driven penalization is employed. Similar results hold for the square-root lasso Belloni2011a,Belloni2014c.
Since coefficient estimates and the set of selected variables depend on $\lambda$ and $\alpha$, a central question is how to choose these tuning parameters. Which method is most appropriate depends on the objectives and setting: in particular, the aim of the analysis (prediction or model identification), computational constraints, and if and how the i.i.d.\ assumption is violated. lassopack\xspace offers three approaches for selecting the penalty level of $\lambda$ and $\alpha$:
Information criteria are closely related to regularization methods. The classical Akaike's information criterion Akaike1974 is defined as $-2 \times \textnormal{log-likelihood} + 2 p$. Thus, the AIC can be interpreted as penalized likelihood which imposes a penalty on the number of predictors included in the model. This form of penalization, referred to as $\ell_0$-penalty, has, however, an important practical disadvantage. In order to find the model with the lowest AIC, we need to estimate all different model specifications. In practice, it is often not feasible to consider the full model space. For example, with only 20 predictors, there are more than 1 million different models.
The advantage of regularized regression is that it provides a data-driven method for reducing model selection to a one-dimensional problem (or two-dimensional problem in the case of the elastic net) where we need to select $\lambda$ (and $\alpha)$. Theoretical properties of information criteria are well-understood and they are easy to compute once coefficient estimates are obtained. Thus, it seems natural to utilize the strengths of information criteria as model selection procedures to select the penalization level.
Information criteria can be categorized based on two central properties: loss efficiency and model selection consistency. A model selection procedure is referred to as loss efficient if it yields the smallest averaged squared error attainable by all candidate models. Model selection consistency requires that the true model is selected with probability approaching 1 as $n\rightarrow\infty$. Accordingly, which information information criteria is appropriate in a given setting also depends on whether the aim of analysis is prediction or identification of the true model.
We first consider the most popular information criteria, AIC and Bayesian information criterion Schwarz1978:
where $\hat\sigma^2(\lambda ,\alpha)=n^{-1}\sum_{i=1}^n\hat\varepsilon_i^2$ and $\hat\varepsilon_i$ are the residuals. ${{df}}(\lambda,\alpha)$ is the effective degrees of freedom, which is a measure of model complexity. In the linear regression model, the degrees of freedom is simply the number of regressors. Zou2007 show that the number of coefficients estimated to be non-zero, $\hat{s}$, is an unbiased and consistent estimate of ${{df}}(\lambda)$ for the lasso ($\alpha=1$). More generally, the degrees of freedom of the elastic net can be calculated as the trace of the projection matrix, i.e., \[ \widehat{{df}}(\lambda,\alpha)=\textrm{tr}(\bm{X}_{\hat{\Omega}}(\bm{X}_{\hat{\Omega}}'\bm{X}_{\hat{\Omega}}+\lambda(1-\alpha)\bm{\Psi})^{-1}\bm{X}_{\hat{\Omega}}').\] where $\bm{X}_{\hat{\Omega}}$ is the $n \times \hat{s}$ matrix of selected regressors. The unbiased estimator of the degrees of freedom provides a justification for using the classical AIC and BIC to select tuning parameters Zou2007.
The BIC is known to be model selection consistent if the true model is among the candidate models, whereas AIC is inconsistent. Clearly, the assumption that the true model is among the candidates is strong; even the existence of the `true model' may be problematic, so that loss efficiency may become a desirable second-best. The AIC is, in contrast to BIC, loss efficient. Yang2005 shows that the differences between AIC-type information criteria and BIC are fundamental; a consistent model selection method, such as the BIC, cannot be loss efficient, and vice versa. Zhang2010 confirm this relation in the context of penalized regression.
Both AIC and BIC are not suitable in the large-$p$-small-$n$ context where they tend to select too many variables (see Monte Carlo simulations in Section (ref)). It is well known that the AIC is biased in small samples, which motivated the bias-corrected AIC Sugiura1978,Hurvich1989, \[ \mathrm{AIC}_c(\lambda,\alpha) = n\log\left(\hat\sigma^2(\lambda ,\alpha)\right) + 2{df}(\lambda,\alpha) \frac{n}{n-{df}(\lambda,\alpha)}. \] The bias can be severe if $df$ is large relative to $n$, and thus the AIC$_c$ should be favoured when $n$ is small or with high-dimensional data.
The BIC relies on the assumption that each model has the same prior probability. This assumptions seems reasonable when the researcher has no prior knowledge; yet, it contradicts the principle of parsimony and becomes problematic if $p$ is large. To see why, consider the case where $p=1000$ Chen2008: There are $1000$ models for which one parameter is non-zero ($s=1$), while there are $1000\times 999/2$ models for which $s=2$. Thus, the prior probability of $s=2$ is larger than the prior probability of $s=1$ by a factor of $999/2$. More generally, since the prior probability that $s=j$ is larger than the probability that $s=j-1$ (up to the point where $j=p/2$), the BIC is likely to over-select variables. To address this shortcoming, Chen2008 introduce the Extended BIC, defined as \[ \mathrm{EBIC}_\xi(\lambda,\alpha) = n\log\left(\hat\sigma^2(\lambda ,\alpha)\right) + {{df}}(\lambda,\alpha)\log(n)+2\xi{{df}}(\lambda,\alpha) \log(p), \] which imposes an additional penalty on the size of the model. The prior distribution is chosen such that the probability of a model with dimension $j$ is inversely proportional to the total number of models for which $s=j$. The additional parameter, $\xi\in[0,1]$, controls the size of the additional penalty.\footnote{We follow Chen2008 and use $\xi=1-\log(n)/(2\log(p))$ as the default choice. An upper and lower threshold is applied to ensure that $\xi$ lies in the [0,1] interval.} Chen2008 show in simulation studies that the EBIC$_\xi$ outperforms the traditional BIC, which exhibits a higher false discovery rate when $p$ is large relative to $n$.
The aim of cross-validation is to directly assess the performance of a model on unseen data. To this end, the data is repeatedly divided into a training and a validation data set. The models are fit to the training data and the validation data is used to assess the predictive performance. In the context of regularized regression, cross-validation can be used to select the tuning parameters that yield the best performance, e.g., the best out-of-sample mean squared prediction error. A wide range of methods for cross-validation are available. For an extensive review, we recommend Arlot2010. The most popular method is $K$-fold cross-validation, which we introduce in Section (ref). In Section (ref), we discuss methods for cross-validation in the time-series setting.
For $K$-fold cross-validation, proposed by Geisser1975, the data is split into $K$ groups, referred to as folds, of approximately equal size. Let $\mathcal{K}_k$ denote the set of observations in the $k$th fold, and let $n_k$ be the size of data partition $k$ for $k=1,...,K$. In the $k$th step, the $k$th fold is treated as the validation data set and the remaining $K-1$ folds constitute the training data set. The model is fit to the training data for a given value of $\lambda$ and $\alpha$. The resulting estimate, which is based on all the data except the observations in fold $k$, is $\bm{\hat\beta}_k(\lambda,\alpha)$. The procedure is repeated for each fold, as illustrated in Figure (ref), so that every data point is used for validation once. The mean squared prediction error for each fold is computed as \[\textrm{MSPE}_k(\lambda,\alpha)=\frac{1}{n_k}\sum_{i \in \mathcal{K}_k}\left(y_i - \bm{x}_i'\bm{\hat\beta}_k(\lambda,\alpha)\right)^2.\]
The $K$-fold cross-validation estimate of the MSPE, which serves as a measure of prediction performance, is \[\hat{\mathcal{L}}^{CV}(\lambda,\alpha)=\frac{1}{K}\sum_{k=1}^K\textrm{MSPE}_k(\lambda,\alpha).\] This suggests selecting $\lambda$ and $\alpha$ as the values that minimize $\hat{\mathcal{L}}^{CV}(\lambda,\alpha)$. An alternative common rule is to use the largest value of $\lambda$ that is within one standard deviation of the minimum, which leads to a more parsimonious model.
Cross-validation can be computationally expensive. It is necessary to compute $\hat{\mathcal{L}}^{CV}$ for each value of $\lambda$ on a grid if $\alpha$ is fixed (e.g.\ when using the lasso) or, in the case of the elastic net, for each combination of values of $\lambda$ and $\alpha$ on a two-dimensional grid. In addition, the model must be estimated $K$ times at each grid point, such that the computational cost is approximately proportional to $K$. \footnote{An exception is the special case of leave-one-out cross-validation, where $K=n$. The advantage of LOO cross-validation for linear models is that there is a closed-form expression for the MSPE, meaning that the model needs to be estimated only once instead of $n$ times.}
Standardization adds another layer of computational cost to $K$-fold cross validation. An important principle in cross-validation is that the training data set should not contain information from the validation dataset. This mimics the real-world situation where out-of-sample predictions are made not knowing what the true response is. The principle applies not only to individual observations, but also to data transformations such as mean-centering and standardization. Specifically, data transformations applied to the training data should not use information from the validation data or full dataset. Mean-centering and standardization using sample means and sample standard deviations for the full sample would violate this principle. Instead, when in each step the model is fit to the training data for a given $\lambda$ and $\alpha$, the training dataset must be re-centered and re-standardized, or, if standardization is built into the penalty loadings, the $\hat\psi_j$ must be recalculated based on the training dataset.
The choice of $K$ is not only a practical problem; it also has theoretical implications. The variance of $\hat{\mathcal{L}}^{CV}$ decreases with $K$, and is minimal (for linear regression) if $K=n$, which is referred to as leave-one-out or LOO CV. Similarly, the bias decreases with the size of the training data set. Given computational contraints, $K$ between 5 and 10 are often recommended, arguing that the performance of CV rarely improves for $K$ larger than 10 Hastie2009, Arlot2010.
If the aim of the researcher's analysis is model identification rather than prediction, the theory requires training data to be `small' and the evaluation sample to be close to $n$ Shao1993,Shao1997. The reason is that more data is required to evaluate which model is the `correct' one rather than to decrease bias and variance. This is referred to as cross-validation paradox Yang2006. However, since $K$-fold cross-validation sets the size of the training sample to approximately $n/K$, $K$-fold CV is necessarily ill-suited for selecting the true model.
Serially dependent data violate the principle that training and validation data are independent. That said, standard $K$-fold cross-validation may still be appropriate in certain circumstances. Bergmeir2018 show that $K$-fold cross-validation remains valid in the pure auto-regressive model if one is willing to assume that the errors are uncorrelated. A useful implication is that $K$-fold cross-validation can be used on overfit auto-regressive models that are not otherwise badly misspecified, since such models have uncorrelated errors.
Rolling $h$-step ahead CV is an intuitively appealing approach that directly incorporates the ordered nature of time series-data Hyndman2018.\footnote{Another approach is a variation of LOO cross-validation known as $h$-block cross-validation Burman1994, which omits $h$ observations between training and validation data.} The procedure builds on repeated $h$-step ahead forecasts. The procedure is implemented in lassopack\xspace and illustrated in Figure (ref)-(ref).
Figure (ref)(a) corresponds to the default case of 1-step ahead cross-validation. `$T$' denotes the observation included in the training sample and `$V$' refers to the validation sample. In the first step, observations 1 to 3 constitute the training data set and observation 4 is the validation point, whereas the remaining observations are unused as indicated by a dot (`.'). Figure (ref)(b) illustrates the case of 2-step ahead cross-validation. In both cases, the training window expands incrementally, whereas Table (ref) displays rolling CV with a fixed estimation window.
Since information-based approaches and cross-validation share the aim of model selection, one might expect that the two methods share some theoretical properties. Indeed, AIC and LOO-CV are asymptotically equivalent, as shown by Stone1977 for fixed $p$. Since information criteria only require the model to be estimated once, they are computationally much more attractive, which might suggest that information criteria are superior in practice. However, an advantage of CV is its flexibility and that it adapts better to situations where the assumptions underlying information criteria, e.g.\ homoskedasticity, are not satisfied Arlot2010. If the aim of the analysis is identifying the true model, BIC and EBIC provide a better choice than $K$-fold cross-validation, as there are strong but well-understood conditions under which BIC and EBIC are model selection consistent.
This section introduces the `rigorous' approach to penalization. Following hdm2016, we use the term `rigorous' to emphasize that the framework is grounded in theory. In particular, the penalization parameters are chosen to guarantee consistent prediction and parameter estimation. Rigorous penalization is of special interest, as it provides the basis for methods to facilitate causal inference in the presence of many instruments and/or many control variables; these methods are the IV-Lasso Belloni2012, the post-double-selection (PDS) estimator Belloni2014b and the post-regularization estimator (CHS) Chernozhukov2015; all of which are implemented in our sister package pdslasso\xspace pdslasso2018.
We discuss the conditions required to derive theoretical results for the lasso in Section (ref). Sections (ref)-(ref) present feasible algorithms for optimal penalization choices for the lasso and square-root lasso under i.i.d., heteroskedastic and cluster-dependent errors. Section (ref) presents a related test for joint significance testing.
There are three main conditions required to guarantee that the lasso is consistent in terms of prediction and parameter estimation.\footnote{For a more detailed treatment, we recommend Hastie2015 and Buhlmann2011.} The first condition relates to sparsity. Sparsity is an attractive assumption in settings where we have a large set of potentially relevant regressors, or consider various different model specifications, but assume that only one true model exists which includes a small number of regressors. We have introduced exact sparsity in Section 2, but the assumption is stronger than needed. For example, some true coefficients may be non-zero, but small in absolute size, in which case it might be preferable to omit them. For this reason, we use a weaker assumption:
\paragraph{Approximate sparsity.} Belloni2012 consider the approximate sparse model (ASM),
The elementary regressors $\bm{w}_i$ are linked to the dependent variable through the unknown and possibly non-linear function $f(\cdot)$. The aim of the lasso (and square-root lasso) estimation is to approximate $f(\bm{w}_i)$ using the target parameter vector $\bm{\beta}_0$ and the transformations $\bm{x}_i:=P(\bm{w}_i)$, where $P(\cdot)$ denotes a dictionary of transformations. The vector $\bm{x}_i$ may be large relative to the sample size, either because $\bm{w}_i$ itself is high-dimensional and $\bm{x}_i:=\bm{w}_i$, or because a large number of transformations such as dummies, polynomials, interactions are considered to approximate $f(\bm{w}_i)$.
The assumption of approximate sparsity requires the existence of a target vector $\bm{\beta}_0$ which ensures that $f(\bm{w}_i)$ can be approximated sufficiently well, while using only a small number of non-zero coefficients. Specifically, the target vector $\bm{\beta}_0$ and the sparsity index $s$ are assumed to meet the condition \[ \lVert\bm{\beta}_0\rVert_0:=s \ll n \qquad \mathrm{with} \quad\frac{s^2\log^2(p\vee n)}{n}\rightarrow 0,\] and the resulting approximation error $r_i=f(\bm{w}_i)-\bm{x}_i'\bm{\beta}_0$ is bounded such that
where $C$ is a positive constant. To emphasize the distinction between approximate and exact sparsity, consider the special case where $f(\bm{w}_i)$ is linear with $f(\bm{w}_i)=\bm{x}_i'\bm{\beta}^\star$, but where the true parameter vector $\bm{\beta}^\star$ violates exact sparsity so that $\lVert\bm{\beta}^\star\rVert_0 > n$. If $\bm{\beta}^\star$ has many elements that are negligible in absolute size, we might still be able to approximate $\bm{\beta}^\star$ using the sparse target vector $\bm{\beta}_0$ as long as $r_i=\bm{x}_i'(\bm{\beta}^\star-\bm{\beta}_0)$ is sufficiently small as specified in (ref).
\paragraph{Restricted sparse eigenvalue condition.} The second condition relates to the Gram matrix, $n^{-1}\bm{X}'\bm{X}$. In the high-dimensional setting where $p$ is larger than $n$, the Gram matrix is necessarily rank-deficient and the minimum (unrestricted) eigenvalue is zero, i.e., \[ \min_{\bm{\delta}\neq 0} \frac{\left\lVert\bm{X}\bm{\delta}\right\rVert_2}{\sqrt{n}\left\lVert\bm{\delta}\right\rVert_2}=0.\] Thus, to accommodate large $p$, the full rank condition of OLS needs to be replaced by a weaker condition. While the full rank condition cannot hold for the full Gram matrix if $p>n$, we can plausibly assume that sub-matrices of size $m$ are well-behaved. This is in fact the restricted sparse eigenvalue (RSEC) condition of Belloni2012. The RSEC formally states that the minimum sparse eigenvalues \[ \phi_{\min}(m) = \min_{1\leq\lVert\bm\delta\rVert_0 \leq m} \frac{\bm\delta'\bm{X}'\bm{X}\bm\delta}{\lVert\bm\delta\rVert^2_2} \qquad \mathrm{and} \qquad \phi_{\max}(m) = \max_{1\leq\lVert\bm\delta\rVert_0 \leq m} \frac{\bm\delta'\bm{X}'\bm{X}\bm\delta}{\lVert\bm\delta\rVert^2_2}\] are bounded away from zero and from above. The requirement $\phi_{\min}(m)>0$ implies that all sub-matrices of size $m$ have to be positive definite.\footnote{The RSEC is stronger than required for the lasso. For example, Bickel2009 introduce the restricted eigenvalue condition (REC). However, here we only present the RSEC which implies the REC and is sufficient for both lasso and post-lasso. Different variants of the REC and RSEC have been proposed in the literature; for an overview see Buhlmann2011.}
\paragraph{Regularization event.} The third central condition concerns the choice of the penalty level $\lambda$ and the predictor-specific penalty loadings $\psi_j$. The idea is to select the penalty parameters to control the random part of the problem in the sense that
with high probability. Here, $c>1$ is a constant slack parameter and $S_j$ is the $j$th element of the score vector $\bm{S} = \nabla \hat{Q}(\bm{\beta})$, the gradient of $\hat{Q}$ at the true value $\bm{\beta}$. The score vector summarizes the noise associated with the estimation problem.
Denote by $\Lambda = n \max_j | \psi_j^{-1} S_j |$ the maximal element of the score vector scaled by $n$ and $\psi_j$, and denote by $q_\Lambda(\cdot)$ the quantile function for $\Lambda$.\footnote{That is, the probability that $\Lambda$ is at most $a$ is $q_\Lambda(a)$.} In the rigorous lasso, we choose the penalty parameters $\lambda$ and $\psi_j$ and confidence level $\gamma$ so that
A simple example illustrates the intuition behind this approach. Consider the case where the true model has $\beta_j=0$ for $j=1,\ldots,p$, i.e., none of the regressors appear in the true model. It can be shown that for the lasso to select no variables, the penalty parameters $\lambda$ and $\psi_j$ need to satisfy $\lambda \geq 2 \max_j |\sum_i \psi_j^{-1} x_{ij} y_i|$.\footnote{See, for example, Hastie2015.} Because none of the regressors appear in the true model, $y_i=\varepsilon_i$. We can therefore rewrite the requirement for the lasso to correctly identify the model without regressors as $\lambda \geq 2 \max_j |\sum_i \psi_j^{-1} x_{ij} \varepsilon_i|$, which is the regularization event in (ref). We want this regularization event to occur with high probability of at least $(1-\gamma)$. If we choose values for $\lambda$ and $\psi_j$ such that $\lambda \geq q_\Lambda(1-\gamma)$, then by the definition of a quantile function we will choose the correct model---no regressors---with probability of at least $(1-\gamma)$. This is simply the rule in (ref).
The chief practical problem in using the rigorous lasso is that the quantile function $q_\Lambda(\cdot)$ is unknown. There are two approaches to addressing this problem proposed in the literature, both of which are implemented in \stcmd{rlasso}\xspace. The \stcmd{rlasso}\xspace default is the `asymptotic' or X-independent approach: theoretically grounded and feasible penalty level and loadings are used that guarantee that (ref) holds asymptotically, as $n\rightarrow\infty$ and $\gamma\rightarrow 0$. The X-independent penalty level choice can be interpreted as an asymptotic upper bound on the quantile function $q_\Lambda(.)$. In the `exact' or X-dependent approach, the quantile function $q_\Lambda(.)$ is directly estimated by simulating the distribution of $q_\Lambda(1-\gamma|\bm{X})$, the $(1-\gamma)$-quantile of $\Lambda$ conditional on the observed regressors $\bm{X}$. We first focus on the X-independent approach, and introduce the X-dependent approach in Section (ref).
Belloni2012 show---using moderate deviation theory of self-normalized sums from Jing2003---that the regularization event in (ref) holds asymptotically, i.e.,
if the penalty levels and loadings are set to
under homoskedasticity and heteroskedasticity, respectively. $c$ is the slack parameter from above and the significance level $\gamma$ is required to converge towards 0. \stcmd{rlasso}\xspace uses $c=1.1$ and $\gamma=0.1/\log(n)$ as defaults. \footnote{The parameters $c$ and $\gamma$ can be controlled using the options c({\it real}) and gamma({\it real}). Note that we need to choose $c$ greater than 1 for the regularization event to hold asymptotically, but not too high as the shrinkage bias is increasing in $c$.}\textsuperscript{,} \footnote{ An alternative X-independent choice is to set $\lambda = 2c\sigma \sqrt{2 n \log(2p/\gamma)}$. Since $\sqrt{n}\Phi^{-1}(1-\gamma/(2p)) \leq \sqrt{2n\log(2p/\gamma)}$, this will lead to a more parsimonious model, but also to a larger bias. To use the alternative X-independent, specify the lalt option. }
\paragraph{Homoskedasticity.} We first focus on the case of homoskedasticity. In the rigorous lasso approach, we standardize the score. But since $E(x_{ij}^2\varepsilon_i^2) = \sigma E(x_{ij}^2)$ under homoskedasticity, we can separate the problem into two parts: the regressor-specific penalty loadings $\psi_j = \sqrt{(1/n) \sum_i x_{ij}^2}$ standardize the regressors, and $\sigma$ moves into the overall penalty level. In the special case where the regressors have already been standardized such that $(1/n)\sum_i x_{ij}^2=1$, the penalty loadings are $\psi_j=1$. Hence, the purpose of the regressor-specific penalty loadings in the case of homoskedasticity is to accommodate regressors with unequal variance.
The only unobserved term is $\sigma$, which appears in the optimal penalty level $\lambda$. To estimate $\sigma$, we can use some initial set of residuals $\hat\varepsilon_{0,i}$ and calculate the initial estimate as $\hat\sigma_0=\sqrt{(1/n)\sum_i \hat\varepsilon_{0,i}^2}$. A possible choice for the initial residuals is $\hat\varepsilon_{0,i}=y_i$ as in Belloni2012 and Belloni2014b. rlasso uses the OLS residuals $\hat\varepsilon_{0,i}=y_i-\bm{x}_i[\mathcal{D}]'\hat{\bm{\beta}}_{OLS}$ where $\mathcal{D}$ is the set of 5 regressors exhibiting the highest absolute correlation with $y_i$.\footnote{This is also the default setting in hdm2016. The number of regressors used for calculating the initial residuals can be controlled using the corrnumber({\it integer}) option, where 5 is the default and 0 corresponds to $\hat\varepsilon_{0,i}=y_i$.} The procedure is summarized in Algorithm A:
The \stcmd{rlasso}\xspace default is to perform one further iteration after the initial estimate (i.e., $K=1$), which in our experience provides good performance. Both lasso and post-lasso can be used to update the residuals. \stcmd{rlasso}\xspace uses post-lasso to update the residuals.\footnote{The lassopsi option can be specified, if rigorous lasso residuals are preferred.}
\paragraph{Heteroskedasticity.} The X-independent choice for the overall penalty level under heteroskedasticity is $\lambda = 2c\sqrt{n}\Phi^{-1}(1-\gamma/(2p))$. The only difference with the homoskedastic case is the absence of $\varsigma$. The variance of $\epsilon$ is now captured via the penalty loadings, which are set to $\psi_j = \sqrt{\frac{1}{n}\sum_i x_{ij}^2\varepsilon_i^2}$. Hence, the penalty loadings here serve two purposes: accommodating both heteroskedasticity and regressors with uneven variance.
To help with the intuition, we consider the case where the predictors are already standardized. It is easy to see that, if the errors are homoskedastic with $\sigma=1$, the penalty loadings are (asymptotically) just $\psi_j=1$. If the data are heteroskedastic, however, the standardized penalty loading will not be $1$. In most practical settings, the usual pattern will be that $\hat\psi_j>1$ for some $j$. Intuitively, heteroskedasticity typically increases the likelihood that the term $\max_j |\sum_i x_{ij}\varepsilon_i|$ takes on extreme values, thus requiring a higher degree of penalization through the penalty loadings. \footnote{To get insights into the nature of heteroskedasticity, \stcmd{rlasso}\xspace also calculates and returns the standardized penalty loadings
which are stored in e(sPsi). }
The disturbances $\varepsilon_i$ are unobserved, so we obtain an initial set of penalty loadings $\hat\psi_j$ from an initial set of residuals $\hat{\varepsilon}_{0,i}$ similar to the i.i.d.\ case above. We summarize the algorithm for estimating the penalty level and loadings as follows:
\paragraph{Theoretical property.} Under the assumptions SEC, ASM and if penalty level $\lambda$ and the penalty loadings are estimated by Algorithm A or B, the lasso and post-lasso obey:\footnote{For the sake of brevity, we omit additional technical conditions relating to the moments of the error and the predictors. These conditions are required to make use of the moderate deviation theory of self-normalized sums Jing2003, which is employed to relax the assumption of Gaussian errors. See condition RF in Belloni2012 and condition SM in Belloni2014b.}
The first relation in (ref) provides an asymptotic bound for the prediction error, and the second relation in (ref) bounds the bias in estimating the target parameter $\bm\beta$. Belloni2012 refer to the above convergence rates as near-oracle rates. If the identity of the $s$ variables in the model were known, the prediction error would converge at the oracle rate $\sqrt{s/n}$. Thus, the logarithmic term $\log(p\vee n)$ can be interpreted as the cost of not knowing the true model.
The theory of the square-root lasso is similar to the theory of the lasso Belloni2011a,Belloni2014c. The $j$th element of the score vector is now defined as
To see why the square-root lasso is of special interest, we define the standardized errors $\nu_i$ as $\nu_i=\varepsilon_i/\sigma$. The $j$th element of the score vector becomes
and is thus independent of $\sigma$. For the same reason, the optimal penalty level for the square-root lasso in the i.i.d.\ case,
is independent of the noise level $\sigma$.
\paragraph{Homoskedasticity.} The ideal penalty loadings under homoskedasticity for the square-root lasso are given by formula (iv) in Table (ref), which provides an overview of penalty loading choices. The ideal penalty parameters are independent of the unobserved error, which is an appealing theoretical property and implies a practical advantage. Since both $\lambda$ and $\psi_j$ can be calculated from the data, the rigorous square-root lasso is a one-step estimator under homoskedasticity. Belloni2011a show that the square-root lasso performs similarly to the lasso with infeasible ideal penalty loadings.
\paragraph{Heteroskedasticity.} In the case of heteroskedasticity, the optimal square-root lasso penalty level remains (ref), but the penalty loadings, given by formula (v) in Table (ref), depend on the unobserved error and need to be estimated. Note that the updated penalty loadings using the residuals $\hat\varepsilon_i$ employ thresholding: the penalty loadings are enforced to be greater than or equal to the loadings in the homoskedastic case. The \stcmd{rlasso}\xspace default algorithm used to obtain the penalty loadings in the heteroskedastic case is analogous to Algorithm B.\footnote{The \stcmd{rlasso}\xspace default for the square-root lasso uses a first-step set of initial residuals. The suggestion of Belloni2014c to use initial penalty loadings for regressor $j$ of $\hat\psi_{0,j} = \max_i |x_{ij}|$ is available using the maxabsx option.} While the ideal penalty loadings are not independent of the error term if the errors are heteroskedastic, the square-root lasso may still have an advantage over the lasso, since the ideal penalty loadings are pivotal with respect to the error term up to scale, as pointed out above.
Kozbur2015 extend the rigorous framework to the case of clustered data, where a limited form of dependence---within-group correlation---as well as heteroskedasticity are accommodated. They prove consistency of the rigorous lasso using this approach in the large $n$, fixed $T$ and large $n$, large $T$ settings. The authors present the approach in the context of a fixed-effects panel data model, $y_{it}=\bm{x}_{it}'\bm{\beta}+\mu_i+\varepsilon_{it}$, and apply the rigorous lasso after the within transformation to remove the fixed effects $\mu_i$. The approach extends to any clustered-type setting and to balanced and unbalanced panels. For convenience we ignore the fixed effects and write the model as a balanced panel:
The intuition behind the Kozbur2015 approach is similar to that behind the clustered standard errors reported by various Stata estimation commands: observations within clusters are aggregated to create `super-observations' which are assumed independent, and these super-observations are treated similarly to cross-sectional observations in the non-clustered case. Specifically, define for the $i$th cluster and $j$th regressor the super-observation $u_{ij} := \sum_t x_{ijt}\varepsilon_{it}$. Then the penalty loadings for the clustered case are $$\psi_j = \sqrt{\frac{1}{nT}\sum_{i=1}^n u_{ij}^2},$$ which resembles the heteroskedastic case. The \stcmd{rlasso}\xspace default for the overall penalty level is the same as in the heteroskedastic case, $\lambda = 2c\sqrt{n}\Phi^{-1}(1-\gamma/(2p))$, except that the default value for $\gamma$ is $0.1/\log(n)$, i.e., we use the number of clusters $n$ rather than the number of observations $nT$. lassopack also implements the rigorous square-root lasso for panel data, which uses the overall penalty in (ref) and the penalty loadings in formula (vi), Table (ref).
There is an alternative, sharper choice for the overall penalty level, referred to as the X-dependent penalty. Recall that the asymptotic, X-independent choice in (ref) can be interpreted as an asymptotic upper bound on the quantile function of $\Lambda$, which is the scaled maximum value of the score vector. Instead of using the asymptotic choice, we can estimate by simulation the distribution of $\Lambda$ conditional on the observed $\bm{X}$, and use this simulated distribution to obtain the quantile $q_\Lambda(1-\gamma|\bm{X})$.
In the case of estimation by the lasso under homoskedasticity, we simulate the distribution of $\Lambda$ using the statistic $W$, defined as \[ W= 2 \hat \sigma \max_{1 \leq j \leq p} \left| \sum_{i=1}^n \psi_j^{-1} x_{ij} g_i \right|, \] where $\psi_j$ is the penalty loading for the homoskedastic case, $\hat \sigma$ is an estimate of the error variance using some previously-obtained residuals, and $g_i$ is an i.i.d.\ standard normal variate drawn independently of the data. The X-dependent penalty choice is sharper and adapts to correlation in the regressor matrix Belloni2011.
Under heteroskedasticity, the lasso X-dependent penalty is obtained by a multiplier bootstrap procedure. In this case the simulated statistic $W$ is defined as in formula (ii) in Table (ref). The cluster-robust X-dependent penalty is again obtained analogously to the heteroskedastic case by defining super-observations, and the statistic $W$ is defined as in formula (iii) in Table (ref). Note that the standard normal variate $g_i$ varies across but not within clusters.
The X-dependent penalties for the square-root lasso are similarly obtained from quantiles of the simulated distribution of the square-root lasso $\Lambda$, and are given by formulas (iv), (v) and (vi) in Table (ref) for the homoskedastic, heteroskedastic and clustered cases, respectively. \footnote{Note that since $g_i$ is standard normal, in practice the term $\frac{1}{\hat\sigma_g}$ that appears in the expressions for the square-root lasso $W$ will be approximately $1$.}
Inference using the lasso, especially in the high-dimensional setting, is a challenging and ongoing research area of research (see Footnote (ref)). A test that has been developed and is available in \stcmd{rlasso}\xspace corresponds to the test for joint significance of regressors using $F$ or $\chi^2$ tests that is common in regression analysis. Specifically, Belloni2012 suggest using the Kato2013 sup-score test to test for the joint significance of the regressors, i.e., \[ H_0: \beta_1 = \ldots = \beta_p = 0. \] As in the preceding sections, regressors are assumed to be mean-centered and in their original units.
If the null hypothesis is correct and the rest of the model is well-specified, including the assumption that the regressors are orthogonal to the disturbance $\varepsilon_i$, then $y_i = \varepsilon_i$ and hence $E(\bm{x}_{i}\varepsilon_i) = E(\bm{x}_i y_i) = 0$. The sup-score statistic is
where $\psi_j = \sqrt{\frac{1}{n}\sum_{i=1}^n (x_{ij} y_i)^2}$. Intuitively, the $\psi_j$ in (ref) plays the same role as the penalty loadings do in rigorous lasso estimation.
The $p$-value for the sup-score test is obtained by a multiplier bootstrap procedure simulating the distribution of $S_{SS}$ by the statistic $W$, defined as \[ W=\sqrt{n} \max_{1 \leq j \leq p} \left| \frac{1}{n}\sum_{i=1}^n \psi_j^{-1} x_{ij} y_i g_i \right|, \] where $g_i$ is an i.i.d. standard normal variate drawn independently of the data and $\psi_j$ is defined as in (ref).
The procedure above is valid under both homoskedasticity and heteroskedasticity, but requires independence. A cluster-robust version of the sup-score test that accommodates within-group dependence is
where $u_{ij} := \sum_{t=1}^T x_{ijt}y_{it}$ and $\psi_j = \frac{1}{nT}\sum_{i=1}^n u_{ij}^2$.
The $p$-value for the cluster version of $S_{SS}$ comes from simulating \[ W = \sqrt{nT} \max_{1 \leq j \leq p} \Big| \frac{1}{nT} \sum_{i=1}^n \psi_j^{-1} u_{ij} g_i \Big| \] where again $g_i$ is an i.i.d. standard normal variate drawn independently of the data. Note again that $g_i$ is invariant within clusters.
rlasso also reports conservative critical values for $S_{SS}$ using an asymptotically-justified upper bound: $CV = c\Phi^{-1}(1-\gamma_{SS}/(2p))$. The default value of the slack parameter is $c=1.1$ and the default test level is $\gamma_{SS}=0.05$. These parameters can be varied using the c({\it real}) and ssgamma({\it real}) options, respectively. The simulation procedure to obtain $p$-values using the simulated statistic $W$ can be computationally intensive, so users can request reporting of only the computationally inexpensive conservative critical values by setting the number of simulated values to zero with the ssnumsim(\it int) option.
The package lassopack\xspace consists of three commands: \stcmd{lasso2}\xspace, \stcmd{cvlasso}\xspace and \stcmd{rlasso}\xspace. Each command offers an alternative method for selecting the tuning parameters $\lambda$ and $\alpha$. We discuss each command in turn and present its syntax. We focus on the main syntax and the most important options. Some technical options are omitted for the sake of brevity. For a full description of syntax and options, we refer to the help files.
The primary purpose of \stcmd{lasso2}\xspace is to obtain the solution of (adaptive) lasso, elastic net and square-root lasso for a single value of $\lambda$ or a list of penalty levels, i.e., for $\lambda_1, \ldots, \lambda_r, \ldots, \lambda_R$, where $R$ is the number of penalty levels. The basic syntax of \stcmd{lasso2}\xspace is as follows:
The options \stcmd{\underbar{alp}ha({\it real})}, \stcmd{sqrt}, \stcmd{\underbar{ada}ptive} and ols can be used to select elastic net, square-root lasso, adaptive lasso and post-estimation OLS, respectively. The default estimator of \stcmd{lasso2}\xspace is the lasso, which corresponds to \stcmd{\underbar{alp}ha(1)}. The special case of \stcmd{alpha(0)} yields ridge regression.
The behaviour of \stcmd{lasso2}\xspace depends on whether the {\it numlist} in \stcmd{\underbar{l}ambda({\it numlist})} is of length greater than one or not. If {\it numlist} is a list of more than one value, the solution consists of a matrix of coefficient estimates which are stored in e(betas). Each row in e(betas) corresponds to a distinct value of $\lambda_r$ and each column to one of the predictors in indepvars. The `path' of coefficient estimates over $\lambda_r$ can be plotted using \stcmd{\underbar{plotp}ath({\it method})}, where {\it method} controls whether the coefficient estimates are plotted against lambda (`lambda'), the natural logarithm of lambda (`lnlambda') or the $\ell_1$-norm (`norm'). If the {\it numlist} in \stcmd{\underbar{l}ambda({\it numlist})} is a scalar value, the solution is a vector of coefficient estimates which is stored in e(b). The default behaviour of \stcmd{lasso2}\xspace is to use a list of 100 values.
In addition to obtaining the coefficient path, \stcmd{lasso2}\xspace calculates four information criteria (AIC, AIC$_c$, BIC and EBIC). These information criteria can be used for model selection by choosing the value of $\lambda_r$ that yields the lowest value for one of the four information criteria. The \stcmd{ic({\it string})} option controls which information criterion is shown in the output, where {\it string} can be replaced with `aic', `aicc', `bic or `ebic'. \stcmd{lic({\it string})} displays the estimation results corresponding to the model selected by an information criterion. It is important to note that \stcmd{lic({\it string})} will not store the results of the estimation. This has the advantage that the user can compare the results for different information criteria without the need to re-estimate the full model. To save the estimation results, the \stcmd{\underbar{postr}esults} option should be specified in combination with \stcmd{lic({\it string})}.
\hangpara \underbar{alp}ha({\it real})\ controls the elastic net parameter, $\alpha$, which controls the degree of $\ell_1$-norm (lasso-type) to $\ell_2$-norm (ridge-type) penalization. alpha(1) corresponds to the lasso (the default estimator), and alpha(0) to ridge regression. The {\it real} value must be in the interval [0,1]. \hangpara sqrt\ specifies the square-root lasso estimator. Since the square-root lasso does not employ any form of $\ell_2$-penalization, the option is incompatible with alpha({\it real}). \hangpara \underbar{ada}ptive\ specifies the adaptive lasso estimator. The penalty loading for predictor $j$ is set to $|\hat\beta_{j,0}|^{-\theta}$ where $\hat\beta_{j,0}$ is the OLS estimator or univariate OLS estimator if $p>n$. $\theta$ is the adaptive exponent, and can be controlled using the \texttt{adatheta({\it real})} option. \hangpara \texttt{\underbar{adal}oadings({\it matrix})}\ is a matrix of alternative initial estimates, $\hat\beta_{j,0}$, used for calculating adaptive loadings. For example, this could be the vector \texttt{e(b)} from an initial \stcmd{lasso2} estimation. The absolute value of $\hat\beta_{j,0}$ is raised to the power $-\theta$ (note the minus). \hangpara \texttt{\underbar{adat}heta({\it real})}\ is the exponent for calculating adaptive penalty loadings. The default is \texttt{\underbar{adat}heta(1)} \hangpara \texttt{ols}\ specifies that post-estimation OLS estimates are displayed and returned in \texttt{e(betas)} or \texttt{e(b)}.
\hangpara \underbar{l}ambda({{\it numlist}}) controls the penalty level(s), $\lambda_r$, used for estimation. {\it numlist}\ is a scalar value or list of values in descending order. Each $\lambda_r$ must be greater than zero. If not specified, the default list is used which, using Mata syntax, is defined by\\[.2cm] \centerline{exp(rangen(log(lmax),log(lminratio*lmax),lcount)),}\\[.2cm] where lcount, lminratio and lmax are defined below, exp() is the exponential function, \texttt{log()} is the natural logarithm, and \texttt{rangen({\it a},{\it b},{\it n})} creates a column vector going from \textit{a} to \textit{b} in \textit{n}-1 steps (see the \texttt{mf_range} help file). Thus, the default list ranges from \texttt{lmax} to \texttt{lminratio*lmax} and \texttt{lcount} is the number of values. The distance between each $\lambda_r$ in the sequence is the same on the logarithmic scale. \hangpara \texttt{\underbar{lc}ount({{\it integer}})} is the number of penalty values, $R$, for which the solution is obtained. The default is \texttt{lcount(100)}. \hangpara \texttt{lmax({{\it real}})} is the maximum value penalty level, $\lambda_1$. By default, $\lambda_1$ is chosen as the smallest penalty level for which the model is empty. Suppose the regressors are mean-centered and standardized, then $\lambda_1$ is defined as $\max_j \frac{2}{n\alpha}\sum_{i=1}^n |x_{ij}y_i|$ for the elastic net and $\max_j \frac{1}{n\alpha}\sum_{i=1}^n |x_{ij}y_i|$ for the square-root lasso Friedman2010. \hangpara \texttt{\underbar{lminr}atio({{\it real}})} is the ratio of the minimum penalty level, $\lambda_R$, to maximum penalty level, $\lambda_1$. {\it real}\ must be between 0 and 1. Default is \texttt{\underbar{lminr}atio(0.001)}.
\hangpara lic({\ststring}) specifies that, after the first \stcmd{lasso2}\xspace estimation using a list of penalty levels, the model that corresponds to the minimum information criterion will be estimated and displayer. `aic', `bic', `aicc', and `ebic' (the default) are allowed. However, the results are not stored in e(). \hangpara \underbar{postr}esults is used in combination with lic({\it string}). \underbar{postr}esults stores estimation results of the model selected by information criterion in e().\footnote{This option was called \texttt{postest} in earlier versions of \textbf{lassopack}\xspace.} \hangpara \texttt{ic({\ststring})} controls which information criterion is shown in the output of \stcmd{lasso2}\xspace when \texttt{lambda()} is a list. 'aic', 'bic', 'aicc', and 'ebic' (the default are allowed). \hangpara \texttt{\underbar{ebicx}i}({\it real}) controls the $\xi$ parameter of the EBIC. $\xi$ needs to lie in the [0,1] interval. $\xi=0$ is equivalent to the BIC. The default choice is $\xi=1-\log(n)/(2\log(p))$.
\hangpara \underbar{notp}en(\varlist) sets penalty loadings to zero for predictors in \varlist. Unpenalized predictors are always included in the model. \hangpara \underbar{par}tial(\varlist) specified that variables in \varlist\ are partialled out prior to estimation. \hangpara \underbar{pload}ings({{\it matrix}}) is a row-vector of penalty loadings, and overrides the default standardization loadings. The size of the vector should equal the number of predictors (excluding partialled-out variables and excluding the constant). \hangpara \underbar{unitl}oadings specifies that penalty loadings be set to a vector of ones; overrides the default standardization loadings. \hangpara prestd specifies that dependent variable and predictors are standardized prior to estimation rather than standardized “on the fly” using penalty loadings. See Section (ref) for more details. By default the coefficient estimates are un-standardized (i.e., returned in original units). \hangpara stdcoef returns coefficients in standard deviation units, i.e., do not un-standardize. Only supported with \texttt{prestd} option.
\hangpara fe within-transformation is applied prior to estimation. The option requires the data in memory to be xtset. \hangpara \underbar{nocons}tant suppress constant from estimation. Default behaviour is to partial the constant out (i.e., to center the regressors).
The replay syntax of \stcmd{lasso2}\xspace allows for plotting and changing display options, without the need to re-run the full model. It can also be used to estimate the model using the value of $\lambda$ selected by an information criterion. The syntax is given by:
\hangpara xb computes predicted values (the default). \hangpara \underbar{r}esiduals computes residuals. \hangpara ols specifies that post-estimation OLS will be used for prediction.
If the previous \stcmd{lasso2}\xspace estimation uses more than one penalty level (i.e. $R>1$), the following options are applicable:
\hangpara \underbar{l}ambda({\it real}) specifies that lambda value used for prediction. \hangpara lid({\it integer}) specifies the index of the lambda value used for prediction. \hangpara \underbar{appr}ox specifies that linear approximation is used instead of re-estimation. Faster, but only exact if coefficient path is piece-wise linear. \hangpara \underbar{noi}sily prompts display of estimation output if re-estimation required. \hangpara \underbar{postr}esults stores estimation results in e() if re-estimation is used.
\stcmd{cvlasso}\xspace implements $K$-fold and $h$-step ahead rolling cross-validation. The syntax of \stcmd{cvlasso}\xspace is:
The \underbar{alp}ha() option of \stcmd{cvlasso}\xspace option accepts a {\it numlist}, while \stcmd{lasso2}\xspace only accepts a scalar. If the {\it numlist}\ is a list longer than one, \stcmd{cvlasso}\xspace cross-validates over $\lambda_r$ with $r=1,\ldots,R$ and $\alpha_m$ with $m=1,\ldots,M$.
plotcv creates a plot of the estimated mean-squared prediction error as a function of $\lambda_r$, and plotopt({\ststring}) can be used to pass plotting options to Stata's line command.
Internally, \stcmd{cvlasso}\xspace calls \stcmd{lasso2}\xspace repeatedly. Intermediate \stcmd{lasso2}\xspace results can be stored using saveest({\ststring}).
\hangpara \underbar{nf}olds({\it integer}) is the number of folds used for $K$-fold cross-validation. The default is \underbar{nf}olds(10), or $K=10$. \hangpara \underbar{foldv}ar({\it varname}) can be used to specify what fold (data partition) each observation lies in. {\it varname}\ is an integer variable with values ranging from 1 to $K$. If not specified, the fold variable is randomly generated such that each fold is of approximately equal size. \hangpara \underbar{savef}oldvar({\it varname}) saves the fold variable variable in {\it varname}. \hangpara seed({\it integer}) sets the seed for the generation of a random fold variable.
\hangpara \underbar{roll}ing uses rolling $h$-step ahead cross-validation. The option requires the data to be \stcmd{tsset} or \stcmd{xtset}. \hangpara h({\it integer}) changes the forecasting horizon. The default is h(1). \hangpara \underbar{or}igin({\it integer}) controls the number of observations in the first training dataset. \hangpara \underbar{fixedw}indow ensures that the size of the training data set is constant.
\hangpara lopt specifies that, after cross-validation, \stcmd{lasso2}\xspace estimates the model with the value of $\lambda_r$ that minimizes the mean-squared prediction error. That is, the model is estimated with $\lambda=\hat\lambda_\texttt{lopt}$. \hangpara lse specifies that, after cross-validation, \stcmd{lasso2}\xspace estimates model with largest $\lambda_r$ that is within one standard deviation from $\hat\lambda_\texttt{lopt}$. That is, the model is estimated with $\lambda=\hat\lambda_\texttt{lse}$. \hangpara \underbar{postr}esults stores the \stcmd{lasso2}\xspace estimation results in e() (to be used in combination with lse or lopt).
Similar to \stcmd{lasso2}\xspace, \stcmd{cvlasso}\xspace also provides a replay syntax, which helps to avoid time-consuming re-estimations. The replay syntax of \stcmd{cvlasso}\xspace can be used for plotting and to estimate the model corresponding to $\hat\lambda_\texttt{lopt}$ or $\hat\lambda_\texttt{lse}$. The replay syntax of \stcmd{cvlasso}\xspace is given by:
\stcmd{rlasso}\xspace implements theory-driven penalization for lasso and square-root lasso. It allows for heteroskedastic, cluster-dependent and non-Gaussian errors. Unlike \stcmd{lasso2}\xspace and \stcmd{cvlasso}\xspace, \stcmd{rlasso}\xspace estimates the penalty level $\lambda$ using iterative algorithms. The syntax of \stcmd{rlasso}\xspace is given by:
\hangpara \underbar{pnotp}en(\varlist) specifies that variables in {\it varlist} are not penalized. \footnote{This option differs from that of \underbar{notp}en(\varlist) as used with cvlasso and lasso2; see the discussion in Section 9.} \hangpara \underbar{rob}ust specifies that the penalty loadings account for heteroskedasticity. \hangpara \underbar{cl}uster({\it varname}) specifies that the penalty loadings account for clustering on variable {\it varname}. \hangpara \texttt{center} center moments in heteroskedastic and cluster-robust loadings. \footnote{For example, the uncentered heteroskedastic loading for regressor $j$ is $\hat\psi_j = \sqrt{\frac{1}{n} \sum_i x_{ij}^2 \hat \varepsilon_i^2}$. In theory, $x_{ij}\varepsilon_i$ should be mean-zero. The centered penalty loading is $\hat\psi_j = \sqrt{\frac{1}{n} \sum_i (x_{ij} \hat \varepsilon_i-\hat\mu)^2}$ where $\hat\mu = \frac{1}{n} \sum_i x_{ij}\hat \varepsilon_i$.} \hangpara \texttt{lassopsi} use lasso or square-root lasso residuals to obtain penalty loadings. The default is post-estimation OLS.\footnote{The option was called \texttt{lassoups} in earlier versions.} \hangpara \texttt{\underbar{corrn}umber({\it integer})} number of high-correlation regressors used to obtain initial residuals. The default is \texttt{\underbar{corrn}umber(5)}, and \texttt{\underbar{corrn}umber(0)} implies that {\it depvar} is used in place of residuals. \hangpara \texttt{prestd} standardize data prior to estimation. The default is standardization during estimation via penalty loadings.
\hangpara \underbar{xdep}endent specifies that the X-dependent penalty level is used; see Section (ref). \hangpara numsim({\it integer}) is the number of simulations used for the X-dependent case. The default is 5,000. \hangpara \underbar{lalt}ernative specifies the alternative, less sharp penalty level, which is defined as $2c\sqrt{2n\log(2p/\gamma)}$ (for the square-root lasso, $2c$ is replaced with $c$). See Footnote (ref). \hangpara gamma({\it real}) is the `$\gamma$' in the rigorous penalty level (default ${\gamma}=1/{\log(n)}$; cluster-lasso default ${\gamma}=1/{\log(n_{clust})}$). See Equation (ref). \hangpara \texttt{c({\it real})} is the `$c$' in the rigorous penalty level. The default is \texttt{c(1.1)}. See Equation (ref).
\hangpara supscore reports the sup-score test of statistical significance. \hangpara testonly reports only the sup-score test without lasso estimation. \hangpara ssgamma({\it real}) is the test level for the conservative critical value for the sup-score test (default = 0.05, i.e., 5% significance level). \hangpara ssnumsim({\it integer}) controls the number of simulations for sup-score test multiplier bootstrap. The default is 500, while 0 implies no simulation.
\hangpara xb generate fitted values (default). \hangpara \underbar{r}esiduals generate residuals. \hangpara lasso use lasso coefficients for prediction (default is to use estimates posted in e(b) matrix). \hangpara ols use OLS coefficients based on lasso-selected variables for prediction (default is to use estimates posted in e(b) matrix).
In this section, we demonstrate the use of \stcmd{lasso2}\xspace, \stcmd{cvlasso}\xspace and \stcmd{rlasso}\xspace using one cross-section example (in Section (ref)) and one time-series example (in Section (ref)).
For demonstration purposes, we consider the Boston Housing Dataset available on the UCI Machine Learning Repository.\footnote{The dataset is available at \href{https://archive.ics.uci.edu/ml/machine-learning-databases/housing/housing.data}{https://archive.ics.uci.edu/ml/machine-learning-databases/housing/ housing.data}, or in CSV format via our website at \href{http://statalasso.github.io/dta/housing.csv}{http://statalasso.github.io/dta/housing.csv}.} The data set includes 506 observations and 14 predictors. \footnote{The following predictors are included: per capita crime rate (crim), proportion of residential land zoned for lots over 25,000 sq.ft. (zn), proportion of non-retail business acres per town (indus), Charles River dummy variable (chas), nitric oxides concentration (parts per 10 million) (nox), average number of rooms per dwelling (rm), proportion of owner-occupied units built prior to 1940 (\texttt{age}), weighted distances to five Boston employment centres (\texttt{dis}), index of accessibility to radial highways (\texttt{rad}), full-value property-tax rate per \$10,000 (\texttt{tax}), pupil-teacher ratio by town (\texttt{pratio}), $1000(Bk - 0.63)^2$ where Bk is the proportion of blacks by town (\texttt{b}), % lower status of the population (\texttt{lstat}), median value (\texttt{medv}).} The purpose of the analysis is to predict house prices using a set of census-level characteristics.
We first employ the lasso estimator:
The above lasso2 output shows the following columns:
Since lambda({\it real}) is not specified, \stcmd{lasso2}\xspace obtains the coefficient path for a default list of $\lambda_r$ values. The largest penalty level is 6858.99, in which case the model does only include the constant. Figure (ref) shows the coefficient path of the lasso for selected variables as a function of $\ln(\lambda)$.\footnote{Figure (ref) was created using the following command:\\ .\ lasso2 medv crim-lstat, plotpath(lnlambda) plotopt(legend(off)) plotlabel plotvar(rm chas rad lstat ptratio dis)}
lasso2 supports model selection using information criteria. To this end, we use the replay syntax in combination with the lic() option, which avoids that the full model needs to be estimated again. The lic() option can also be specified in the first \stcmd{lasso2}\xspace call. In the following example, the replay syntax works similar to a post-estimation command.
Two columns are shown in the output; one for the lasso estimator and one for post-estimation OLS, which applies OLS to the model selected by the lasso.
Next, we consider $K$-fold cross-validation.
The \stcmd{cvlasso}\xspace output displays four columns: the index of $\lambda_r$ (i.e., $r$), the value of $\lambda_r$, the estimated mean squared prediction error, and the standard deviation of the mean squared prediction error. The output indicates the value of $\lambda_r$ that corresponds to the lowest MSPE with an asterisk (\verb+*+). We refer to this value as $\hat\lambda_{\texttt{lopt}}$. In addition, the symbol \verb+^+ marks the largest value of $\lambda$ that is within one standard error of $\hat\lambda_{\texttt{lopt}}$, which we denote as $\hat\lambda_{\texttt{lse}}$.
The mean squared prediction is shown in Figure (ref), which was created using the plotcv option. The graph shows the mean squared prediction error estimated by cross-validation along with $\pm$ one standard error. The continuous and dashed vertical lines correspond to $\hat\lambda_{\texttt{lopt}}$ and $\hat\lambda_{\texttt{lse}}$, respectively.
To estimate the model corresponding to either $\hat\lambda_{\texttt{lopt}}$ or $\hat\lambda_{\texttt{lse}}$, we use the lopt or lse option, respectively. Similar to the lic() option of \stcmd{lasso2}\xspace, lopt and lse can either specified in the first \stcmd{cvlasso}\xspace call or after estimation using the replay syntax as in this example:
Lastly, we consider \stcmd{rlasso}\xspace. The program \stcmd{rlasso}\xspace runs an iterative algorithm to estimate the penalty level and loadings. In contrast to \stcmd{lasso2}\xspace and \stcmd{cvlasso}\xspace, it reports the selected model directly.
The supscore option prompts the sup-score test of joint significance. The $p$-value is obtained through multiplier bootstrap. The test statistic of 16.59 can also be compared to the asymptotic 5% critical value (here 3.18).
A standard problem in time-series econometrics is to select an appropriate lag length. In this sub-section, we show how lassopack\xspace can be employed for this purpose. We consider Stata's built-in data set lutkepohl2.dta, which includes quarterly (log-differenced) consumption (dln_consump), investment (dln_inv) and income (dln_inc) series for West Germany over the period 1960, Quarter 1 to 1982, Quarter 4. We demonstrate both lag selection via information criteria and by $h$-step ahead rolling cross-validation. We do not consider the rigorous penalization approach of rlasso due to the assumption of independence, which seems too restrictive in the time-series context.
After importing the data, we run the most general model with up to 12 lags of dln_consump, dln_inv and dln_inc using \stcmd{lasso2}\xspace with lic(aicc) option.
The output consists of two parts. The second part of the output is prompted since \stcmd{lic(aicc)} is specified. \stcmd{lic(aicc)} asks \stcmd{lasso2}\xspace to estimate the model selected by AIC$_c$, which in this case corresponds to $\lambda_{11}=0.207$.
In the next step, we consider $h$-step ahead rolling cross-validation.
The output indicates how the data set is partitioned into training data and the validation point. For example, the short-hand 13-80 (81) in the output above indicates that observations 13 to 80 constitute the training data set in the first step of cross-validation, and observation 81 is the validation point. The options \stcmd{fixedwindow}, h(integer) and origin(integer\texttt{)} can be used to control the partitioning of data into training and validation data. \texttt{h(}\textit{integer}\texttt{)} sets the parameter $h$. For example, \stcmd{h(2)} prompts 2-step ahead forecasts (the default is \stcmd{h(1)}). \stcmd{fixedwindow} ensures that the training data set is of same size in each step. If \stcmd{origin(50)} is specified, the first training partition includes observations 13 to 50, as shown in the next example:
The option plotcv creates the graph of the estimated mean squared prediction in Figure (ref). To estimate the model corresponding to $\hat\lambda_{\texttt{lse}}$, we can as in the previous examples use the replay syntax:
We point out that care should be taken when setting the parameters of $h$-step ahead rolling cross-validation. The default settings have no particular econometric justification.
We have introduced three alternative approaches for setting the penalization parameters in Sections (ref)-(ref). In this section, we present results of Monte Carlo simulations which assess the performance of these approaches in terms of in-sample fit, out-of-sample prediction, model selection and sparsity. To this end, we generate artificial data using the process
with $n=200$. We report results for $p=100$ and for the high-dimensional setting with $p=220$. The predictors $x_{ij}$ are drawn from a multivariate normal distribution with $\textnormal{corr}(x_{ij},x_{ir})=0.9^{|j-r|}$. We vary the noise level $\sigma$ between $0.5$ and $5$; specifically, we consider $\sigma=\{0.5,1,2,3,5\}$. We define the parameters as $\beta_j = \mathbbm{1}\{j\leq s\}$ for $j=1,\ldots,p$ with $s=20$, implying exact sparsity. This simple design allows us to gain insights into the model selection performance in terms of false positive (the number of variables falsely included) and false negative frequency (the number of variables falsely omitted) when relevant and irrelevant regressors are correlated. All simulations use at least 1,000 iterations. We report additional Monte Carlo results in Appendix (ref), where we employ a design in which coefficients alternate in sign.
Since the aim is to assess in-sample and out-of-sample performance, we generate $2n$ observations, and use the data $i=1,\ldots,n$ as the estimation sample and the observations $i=n+1,\ldots,2n$ for assessing out-of-sample prediction performance. This allows us to calculate the root mean squared error (RMSE) and root mean squared prediction error (RMSPE) as
where $\hat{y}_{i,n}$ are the predictions from fitting the model to the first $n$ observations.
Table (ref) and (ref) report results for the following estimation methods implemented in lassopack\xspace: the lasso with $\lambda$ selected by AIC, BIC, EBIC$_\xi$ and AIC$_c$ (as implemented in \stcmd{lasso2}\xspace); the rigorous lasso and rigorous square-root lasso (implemented in \stcmd{rlasso}\xspace), both using the X-independent and X-dependent penalty choice; and lasso with 5-fold cross-validation using the penalty level that minimizes the estimated mean squared prediction error (implemented in \stcmd{cvlasso}\xspace). In addition, we report post-estimation OLS results.
For comparison, we also show results of stepwise regression (for $p=100$ only) and the oracle estimator. Stepwise regression starts from the full model and iteratively removes regressors if the $p$-value is above a pre-defined threshold (10% in our case). Stepwise regression is known to suffer from overfitting and pre-testing bias. However, it still serves as a relevant reference point due to its connection with ad hoc model selection using hypothesis testing and the general-to-specific approach. The oracle estimator is OLS applied to the predictors included in the true model. Naturally, the oracle estimator is expected to show the best performance, but is not feasible in practice since the true model is not known.
We first summarize the main results for the case where $p=100$; see Table (ref). AIC and stepwise regression exhibit the worst selection performance, with around 18-20 falsely included predictors on average. While AIC and stepwise regression achieve the lowest RMSE (best in-sample fit), the out-of-sample prediction performance is among the worst---a symptom of over-fitting. It is interesting to note that the RMSE of AIC and stepwise regression are lower than the RMSE of the oracle estimator. The corrected AIC improves upon the standard AIC in terms of bias and prediction performance.
Compared to AIC$_c$, the BIC-type information criteria show similar out-of-sample prediction and better selection performance. While the EBIC performs only marginally better than BIC in terms of false positives and bias, we expect the relative performance of BIC and EBIC to shift in favour of EBIC as $p$ increases relative to $n$. 5-fold CV with the lasso behaves very similarly to the AIC$_c$ across all measures. The rigorous lasso, rigorous square-root lasso and EBIC exhibit overall the lowest false positive rates, whereas rigorous methods yield slightly higher RMSE and RMSPE than IC and CV-based methods. However, post-estimation OLS (shown in parentheses) applied to the rigorous methods improves upon first-step results, indicating that post-estimation OLS successfully addresses the shrinkage bias from rigorous penalization. The performance difference between X-dependent and X-independent penalty choices are minimal overall.
We also present simulation results for the high-dimensional setting in Table (ref). Specifically, we consider $p=220$ instead of $p=100$, while keeping the estimation sample size constant at $n=200$. With on average between 164 and 195 included predictors, it is not surprising that the AIC suffers from overfitting. The RMSPE of the AIC exceeds the RMSE by a factor of 5 or more. In comparison, AIC$_c$ and 5-fold cross-validation perform better as model selectors, with a false positive frequency between 6 and 8 predictors.
Despite the large number of predictors to choose from, EBIC and rigorous methods perform generally well in recovering the true structure. The false positive frequency is below 1 across all noise levels, and the false negative rate is zero if $\sigma$ is 1 or smaller. While the BIC performs similarly to the EBIC for $\sigma=0.5$ and $\sigma=1$, its performance resembles the poor performance of AIC for larger noise levels. The Monte Carlo results in Table (ref) highlight that EBIC and rigorous methods are well-suited for the high-dimensional setting where $p>n$, while AIC and BIC are not appropriate.
The computational costs of each method are reported in Table (ref). \stcmd{rlasso}\xspace with X-independent penalty is the fastest method considered. The run-time of lasso and square-root lasso with $p=100$ is 0.1s and 0.4s, respectively. The computational cost increased only slightly if $p$ is increased to $p=220$. \stcmd{rlasso}\xspace with X-dependent penalty simulates the distribution of the maximum value of the score vector. This process increases the computational cost of the rigorous lasso to 5.9s for $p=100$ (12.7s for $p=220$). With an average run-time of 3.1 seconds, \stcmd{lasso2}\xspace is slightly faster than \stcmd{rlasso}\xspace with X-dependent penalty if $p=100$, but slower in the high-dimensional set-up. Unsurprisingly, $K$-fold cross-validation is the slowest method as it requires the model to be estimated $K$ times for a range of tuning parameters.
lassopack implements the elastic net and square-root lasso using coordinate descent algorithms. The algorithm---then referred to as “shooting”---was first proposed by Fu1998 for the lasso, and by Van der Kooij (2007) for the elastic net. Belloni2011a and Belloni2014c employ the coordinate descent for the square-root lasso, and have kindly provided Matlab code.
Coordinate descent algorithms repeatedly cycle over predictors $j=1,...,p$ and update single coefficient estimates until convergence. Suppose the predictors are centered, standardized to have unit variance and the penalty loadings are $\psi_j=1$ for all $j$. In that case, the update for coefficient $j$ is obtained using univariate regression of the current partial residuals (i.e., excluding the contribution of predictor $j$) against predictor $j$. More precisely, the update for the elastic net is calculated as \[\tilde\beta_j \leftarrow \frac{\mathcal{S}\left(\sum_{i=1}^nx_{ij}(y_i-\tilde{y}^{(j)}_i),\lambda\alpha \right)}{1+\lambda(1-\alpha)}.\] where $\tilde\beta_j$ denotes the current coefficient estimate, $\tilde{y}^{(j)}_i=\sum_{\ell\neq j}x_{i\ell}\tilde\beta_\ell$ is the predicted value without the contribution of predictor $j$. Thus, since the predictors are standardized, $\sum_{i}x_{ij}(y_i-\tilde{y}^{(j)}_i)$ is the OLS estimate of regressing predictor $j$ against the partial residual $(y_i-\tilde{y}^{(j)}_i)$. The function $\mathcal{S}(a,b)$, referred to as soft-tresholding operator, \[ \mathcal{S}(a,b)=\left\{
\right. \] sets some of the coefficients equal to zero. The coordinate descent algorithm is spelled out for the square-root lasso in Belloni2014c.\footnote{Alexandre Belloni provides MATLAB code that implements the pathwise coordinate descent for the square-root lasso, which we have used for comparison.}
The algorithm requires an initial beta estimate for which the Ridge estimate is used. If the coefficient path is obtained for a list of $\lambda$ values, \stcmd{lasso2}\xspace starts from the largest $\lambda$ value and uses previous estimates as initial values (`warm starts'). See Friedman2007,Friedman2010, and references therein, for further information.
Since penalized regression methods are not invariant to scale, it is common practice to standardize the regressors $x_{ij}$ such that $\sum_ix_{ij}^2=1$ before computing the estimation results and then to un-standardize the coefficients after estimation. We refer to this approach as pre-estimation standardization. An alternative is to standardize on the fly by adapting the penalty loadings. The results are equivalent in theory. In the case of the lasso, setting $\psi_j=(\sum_ix_{ij}^2)^{1/2}$ yields the same results as dividing the data by $\sum_ix_{ij}^2$ before estimation. Standardization on-the-fly is the default in lassopack as it tends to be faster. Pre-estimation standardization can be employed using the prestd option. The prestd option can lead to improved numerical precision or more stable results in the case of difficult problems; the cost is (a typically small) computation time required to standardize the data. The \underbar{unitl}oadings option can be used if the researcher does not want to standardize data. In case the pre-estimation-standardization and standardization-\emph{on-the-fly} results differ, the user can compare the values of the penalized minimized objective function saved in \texttt{e(pmse)} (the penalized MSE, for the elastic net) or \texttt{e(prmse)} (the penalized root MSE, for the sqrt-lasso).
In many applications, theory suggests that specific predictors have an effect on the outcome variable. Hence, it might be desirable to always include these predictors in the model in order to improve finite sample performance. Typical examples are the intercept, a time trend or any other predictor for which the researcher has prior knowledge. lassopack offers two approaches for such situations:
It turns out that the two methods---zero-penalization and partialling out---are numerically equivalent. Formally, suppose we do not want to subject predictors $\ell$ with $\bar{p}>\ell\geq p$ to penalization. The zero-penalization and partialled-out lasso estimates are defined respectively as
where $\tilde{y}_i=y_i-\sum_{\ell=\bar{p}+1}^{p}x_{i\ell}\hat\delta_{y,\ell}$ and $\tilde{x}_{ij}=x_{ij}-\sum_{\ell=\bar{p}+1}^{p}x_{i\ell}\hat\delta_{j,\ell}$ are the residuals of regressing $y$ and the penalized regressors against the set of unpenalized regressors. The equivalence states that $\hat\beta_j=\tilde\beta_j$ for all $j=1,\ldots,\bar{p}$. The result is spelled out in Yamada2017 for the lasso and ridge, but holds for the elastic net more generally.
Either the \underbar{par}tial(\varlist) option or the \underbar{notp}en(\varlist) option can be used for variables that should not be penalized by the lasso. The options are equivalent in theory (see above), but numerical results can differ in practice because of the different calculation methods used. Partialling-out variables can lead to improved numerical precision or more stable results in the case of difficult problems compared to zero-penalization, but may be slower in terms of computation time.
The estimation of penalty loadings in the rigorous lasso introduces an additional complication that necessitates the \stcmd{rlasso}\xspace-specific option \underbar{pnotp}pen(\varlist). The theory for the \stcmd{rlasso}\xspace penalty loadings is based on the penalized regressors after partialling out the unpenalized variables. The \underbar{pnotp}en(\varlist) guarantees that the penalty loadings for the penalized regressors are the same as if the unpenalized regressors had instead first been partialled-out.
The fe fixed-effects option is equivalent to (but computationally faster and more accurate than) specifying unpenalized panel-specific dummies. The fixed-effects (`within') transformation also removes the constant as well as the fixed effects. The panel variable used by the fe option is the panel variable set by xtset. If installed, the within transformation uses the fast ftools package by Correia2016.
The prestd option, as well as the \underbar{notp}en(\varlist) and \underbar{pnotp}en(\varlist) options, can be used as simple checks for numerical stability by comparing results that should be equivalent in theory. The values of the penalized minimized objective function saved in e(pmse) for the elastic net and e(prmse) for the square-root lasso may also be used for comparison.
By default the constant, if present, is not penalized; this is equivalent to mean-centering prior to estimation. The \underbar{par}tial(\varlist) option also partials out the constant (if present). To partial out the constant only, we can specify \underbar{par}tial(_cons). Both \underbar{par}tial(\varlist) and fe mean-center the data; the \underbar{nocons}tant option is redundant in this case and may not be specified with these options. If the \underbar{nocons}tant option is specified an intercept is not included in the model, but the estimated penalty loadings are still estimated using mean-centered regressors (see the \texttt{center} option).
We thank Alexandre Belloni, who has provided MATLAB code for the square-root lasso, and Sergio Correia for supporting us with the use of ftools. We also thank Christopher F Baum, Jan Ditzen, Martin Spindler, as well as participants of the 2018 London Stata Conference and the 2018 Swiss Stata Users Group meeting for many helpful comments and suggestions. All remaining errors are our own.