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.
48,980 characters · 16 sections · 15 citation commands
High-dimensional Penalized Linear IV Estimation & Inference using BRIDGE and Adaptive LASSO
\onehalfspacing
\affil{Northwestern University, USA}
\setcounter{dispsection}{-1} \setcounter{theorem}{0}
Using increasingly large datasets, economists and social scientists in general face new challenges that are connected with the high-dimensional nature of the problems they study. A particular instance of such a scenario is the class of applications that involve high-dimensional instruments. Valid instruments are often the holy grail of applied research; however, recently there is increased availability in rich datasets such as the GWAS Catalog Macarthur:2017, Sollis:2023, which includes information on how phenotypic variations affect biological traits and disease traits give new paths for research. This type of variables can provide a number of instrument alternatives for applications related to productivity, risk, welfare, actuarial science, etc Gui:2005, Ma:2007, Wang:2008. In order to be able to do estimation and inference in environments like this, the researcher has to address the statistical difficulties that occur.
I will focus on the usual linear model \[Y=X'\beta_0 +u\] where the usual exogeneity restriction, i.e. $E[u|X]= 0$ fails. Instead, I will assume that there exist observable variables $Z$ that can be used as instruments, such that $E[u|Z]= 0$, ${\rm Cov}[ZX]$ is full rank, and the conditional expectation $E[X|Z]$ is linear. I am allowing the number of both the endogenous covariates and the instruments to grow and I am making the high level assumption that the number of instruments is sufficient for the identification of the model.
Due to the dimensionality problem, if the number of covariates in each stage exceeds the sample size, the OLS estimator is not identified. In this case, one can obtain valid estimators by minimizing the penalized loss function, i.e. \(L_n(b)= \operatorname*{arg\,min}_b \frac{1}{2}\sum_{i=1}^{n}(Y_i-X_i'b)^2 +\lambda_n\sum_{j}^{p}p(b_j)\) where $p(b_j)$ is a penalty function. Due to the presence of endogeneity, instead of $X_i$, I will use an estimator of $E[X|Z]$. Under sparsity, the latter can be obtained by using any proper penalization function. Then, I will estimate the model choosing the following popular penalty functions: BRIDGE , with $p(b)=|b|^{\gamma}$, where $\gamma\in (0,1)$ and adaptive LASSO with $p(b)=w_b|b|$, where $w_b$ is a set of weights.
Related literature: The standard practice to estimate high-dimensional, sparse, linear models includes popular penalization functions, such as, for $\gamma=1$, LASSO Tibshirani:1996, adaptive LASSO Zou:2006, BRIDGE Knight:2000, the smoothly clipped absolute deviation (SCAD) Fan:2001, the minimax concave pentaly (MCP) Zhang:2010 etc. However, inference in these setups often remains a difficult problem even with the sole existence of exogenous covariates Van:2014. In environments that endogeneity is also an issue, the econometrician faces three distinct scenarios; the first is having a growing number of instruments, yet a small number of covariates of interest Belloni:2012,Caner:2015, Fan:2018, Hansen:2014. The second is allowing for a growing number of covariates but a low-dimensional set of them being endogenous, which permits for a low dimensional estimation in the first stage (eg. \citeasnoun{Fan:2014}). Finally, there is the case with high-dimensionality in both stages with a growing number of instruments and regressors. \citeasnoun{Lin:2015} allow for this setting and they prove model selection consistency and oracle efficiency using LASSO, SCAD, and MCP, under the restrictive assumption of normality. \citeasnoun{Caner:2014} propose a GMM estimator using adaptive elastic nets as a penalty function, allowing for a growing number of covariates but without exceeding the sample size.\citeasnoun{Gold:2020} follows a similar direction, providing general assumptions under which a penalty function would have the desirable properties but restricting the nature of the instruments to sub-Gaussian. Also, for LASSO, the required properties imply the presence of strong compatibility conditions to achieve the oracle property.
To address the need for a more flexible framework, a natural path would be to turn to a penalty function that ensures model selection consistency with minimal assumptions in the simple OLS case, and consider the extension to the endogeneity framework. BRIDGE and adaptive LASSO are natural candidates for this setup. There have been some attempts Bahador:2024 to apply BRIDGE in two-stage penalized least squares using a control function approach, yet there is no formal proof of model selection and they require normality for the error terms.
Contribution: The main contribution of this paper is showing that the estimator for $\beta_0$ using both penalty functions is model selection consistent and oracle efficient. For both methods, if the number of parameters grows faster than the sample size, the result is based only on sub-Gaussian tails of the error term and standard rate conditions. The argument for BRIDGE is an extension to the partial orthogonality condition in \citeasnoun{Huang:2008}. For adaptive LASSO, the proof is based on the requirement of a consistent estimator for the weights in the penalty function as stated for the exogenous case in \citeasnoun{Huang:2008b}. Moreover, I show that for a growing number of parameters but not higher than the sample size, the result for BRIDGE holds with just homoskedastic, mean 0 errors and the corresponding rate assumptions. Lastly, for BRIDGE, I provide a proposal of a consistent estimator for the standard errors as well as some computational evidence on the performance of the method. Corresponding exercises on adaptive LASSO are ongoing and available upon request.
Structure: The rest of the paper is organized as follows. Section 2 introduces the notation and the requirements for the 1st stage of the problem, regardless of the 2nd stage method. In Section 3, I present the results for BRIDGE when $p<n$. In Sections 4 and 5, I present the results for $p>n$ for BRIDGE and adaptive LASSO correspondingly. Finally, Section 6 contains some computational results comparing how BRIDGE and LASSO perform in this problem.
\setcounter{dispsection}{0} \setcounter{theorem}{0}
Consider the vector $(Y, X,Z )$, where $Y$ is the outcome random variable, $X$ is the random vector of observable covariates - possibly endogenous, and $Z$ is a fixed vector of instruments. The researcher observes a random sample of size $n$. Assume that the true conditional expectation $E[X|Z]$ is linear.
This paper is concerned with the following two-stage model:
where $X_i=(X_{i1}, \dots, X_{ip_{xn}})'\in \mathbb{R}^{p_{xn}}$ consists of the covariates in the $2^{nd}$ stage with $\beta_0\in \mathbb{R}^{p_{xn}}$ being the corresponding coefficients, and $Z_i\in \mathbb{R}^{p_{zn}}$ consists of the covariates in the $1^{st}$ stage with $\alpha_j\in \mathbb{R}^{p_{zn}}$ being the corresponding coefficients. Let $\alpha\in \mathbb{R}^{p_{zn}\times p_{xn}}$ the matrix of the stacked first stage coefficients. I define the respective conditional means as $d_{ij}:=E[X_{ij}|Z_i]=Z_i'\alpha_j $. Let the matrix of the second stage covariates be $X_n\in \mathbb{R}^{n\times p_{xn}}$. Accordingly, let the matrix of the first stage covariates be $Z_n\in \mathbb{R}^{n\times p_{zn}}$. Further, I define $D_n=E[X_n|Z_n]=Z\alpha\in \mathbb{R}^{n\times p_{xn}}$ , which rows are $d_i=(d_{i1}, \dots, d_{ip_x})'$.
Naturally, $u_i\ and\ \bm{v}_{i}:= (v_{i1}, v_{i2},\dots, v_{ip_x})'$ are the random noise terms in each stage, and satisfy $E[u_i|Z_i]=0,\ E[v_i|Z_i]=0$ but not necessarily $E[u_i|X_{ij}]=0$. Finally, note that the model above implies that $Y_i= (Z_i'\alpha)'\beta_0+v_i'\beta_0+u_i= d_i'\beta_0+\underbrace{v_i'\beta_0+u_i}$, so I define $\varepsilon_{i}=v_i'\beta_0+u$ the new error term that satisfies $E[\varepsilon_i|Z_{i}]=0$.
Define $X_{1n}\in\mathbb{R}^{n\times k_{xn}}$ the matrix of the relevant covariates in the 2nd stage and $X_{2n}\in\mathbb{R}^{n\times m_{xn}}$ the matrix of the rest. Let $\beta_{01}$ be the sub-vector of 2nd stage coefficients that are non $0$, and $\beta_{02}$ vector of the rest of elements of $\beta_0$. Similarly for the first stage; $Z_{1nj}\in \mathbb{R}^{n\times k_{znj}}$ is the matrix of relevant covariates corresponding to the 2nd stage covariate $j$ and $Z_{2nj}\in\mathbb{R}^{n\times m_{knj}}$ the matrix of the rest. Let $\alpha_{1j}$ be the sub-vector of the non-zero coefficients of $\alpha_j$ $\forall j$, and $\alpha_{2j}$ the sub-vector of the zero coefficients. The number of instruments and the number covariates are allowed to grow with the sample size but the model should be sparse in both stages; that is, $k_{xn}/n\to 0 $ and $\max k_{znj}/n\to 0$.
For simplicity of notation, let $Y_i$ be centered and the instruments be standardized, ie.
\[ \sum_{i=1}^{n}Y_i=0, \quad \sum_{i=1}^{n}Z_{ih}=0, \quad \frac{1}{n}\sum_{i=1}^{n}Z_{ih}^2=1, \] for all instruments $h=1,\dots,p_{zn}$.
Lastly, for a given vector $\delta$, $\left\lVert\delta\right\rVert$ is the Euclidean norm, and for a given matrix $\Delta$, $\left\lVert\Delta\right\rVert$ is the spectral norm.
Due to existing results, the 1st stage coefficients are easy to deal with. For their estimation, the researcher can use the method of choice that provides model selection consistency and oracle efficiency. However, the assumptions that are needed in order to use the said method will directly affect the distribution of the 2nd stage BRIDGE coefficients. Define $\Sigma_{zn}=n^{-1}Z_n'Z_n$. Let $\rho_{1n}^z$ be the smallest and $\rho_{2n}^{z}$ the largest eigenvalue of $\Sigma_{zn}$. Define $\Sigma_{1nj}=n^{-1}Z_{1nj}'Z_{1nj}$. Let $\tau_{1n}^z$ be the smallest and $\tau_{2n}^{z}$ the largest eigenvalue of $\Sigma_{1nj}$. A sufficient set of conditions for the 1st stage, in order to achieve the desirable properties in the 2nd stage, are the following:
For each $j$:
(A1) ensures the good behavior of the composite error term and allows for a complicated dependence structure between the regressors. The second assumption ensures the good behavior of the weighted Gram matrices $\Sigma_{zn},\ \Sigma_{1nj}$. Assumption (A3) ensures that the small coefficients are far enough from 0 and the large ones do not diverge. Finally, (A4) is the most high level assumption, ensuring knowing the limit distribution of $\hat{\alpha}$. Since this is a regular high-dimensional setting, this is not a difficult assumption to satisfy and it can be replaced with the set of assumptions of a specific choice of penalty function. For example, for BRIDGE, (A1-A3) and a set of rate assumption is sufficient for (A4) Huang:2008. For fewer instruments than the sample size, the assumption on the error term will be exactly the same, while for $p_{zn}>n$ the errors of the 1st stage should also be sub-Gaussian to satisfy the conditions of the aforementioned paper.
\setcounter{dispsection}{1} \setcounter{theorem}{0}
Having ensured a good estimator for the first stage, I only need to provide an estimator with good properties for the second stage. I first state the results for $p_{xn}<n$. In the case of $p_{xn}>n$, these results will be used after ensuring model selection consistency. That is, after picking the correct model, due to sparsity, I can apply the same results as in $p_{xn}<n$, and ensure consistency and oracle efficiency even in the most cumbersome case.
Let $\hat{\beta}_j\in\mathbb{R}^{p_{xn}}$ be the root of the following minimization problem: \[\hat{\beta}=\operatorname*{arg\,min}_{b}\underbrace{ \frac{1}{2}\sum_{i=1}^{n}(Y_i-\hat d_i'b)^2+\lambda_{xn}\sum_{j=1}^{p_{xn}}|b_{j}|^\gamma}_{L_n(b)}.\]
where $\lambda_{xn}$ is the tuning parameter. Define $\Sigma_{dn}=n^{-1}D_n'D_n$. Let $\rho_{1n}^d$ be the smallest and $\rho_{2n}^{d}$ the largest eigenvalue of $\Sigma_{dn}$. Let $D_{1n}\in\mathbb{R}^{n\times k_{xn}}$ be the matrix of the conditional expectations of the relevant covariates in the 2nd stage and $D_{2n}\in\mathbb{R}^{n\times m_{xn}}$ the matrix of the rest. Define $\Sigma_{1dn}=n^{-1}D_{1n}'D_{1n}$. Let $\tau_{1n}^d$ be the smallest and $\tau_{2n}^{d}$ the largest eigenvalue of $\Sigma_{1dn}$.
Also, define $\hat \Sigma_{dn}=n^{-1}\hat D_n'\hat D_n$. Let $\hat\rho_{1n}^d$ be the smallest and $\hat\rho_{2n}^{d}$ the largest eigenvalue of $\hat\Sigma_{dn}$. Let $\hat D_{1n}\in\mathbb{R}^{n\times k_{xn}}$ the matrix of the estimated conditional expectations of the relevant covariates in the 2nd stage and $\hat D_{2n}\in\mathbb{R}^{n\times m_{xn}}$ the matrix of the rest. Define $\hat \Sigma_{1dn}=n^{-1}\hat D_{1n}'\hat D_{1n}$. Let $\hat\tau_{1n}^d$ be the smallest and $\hat\tau_{2n}^{d}$ the largest eigenvalue of $\hat\Sigma_{1dn}$. Let $\hat d_{1i}:=Z_i'\hat\alpha $ a row of $\hat D_{1n}$ and the rest of the vectors accordingly.
Assumption (B.1) is of the same nature as (A.1). The 2nd stage error term is homoskedastic and the covariance with each one of the 1st stage error terms needs to be finite but is allowed to be non-zero. Combining the expressions of the two stages, the error term of interest is $\varepsilon_i=v_{i}'\beta_0+u_i$ which is i.i.d. with mean 0 and variance $\sigma_{\varepsilon}^2$. Assumption (B.2) ensures the invertibility of $\Sigma_{dn}$ but, contrary to the smallest eigenvalue of $\Sigma_{zn}$, $\rho_{1n}^{d}$ is not allowed to converge to 0. A direct consequence is that the convergence rates of the 2nd stage estimator will no longer depend on the eigenvalues of the corresponding Gram matrix in the way the 1st stage ones do. (B.3.b) is used for the convergence of the 2nd stage estimator. (B.3.b) ensures model selection consistency and the rest of the rate assumptions are used for the three upcoming results. Note that the rates are affected by the number of the instruments and the maximum number of relevant coefficients of the 1st stage separate equations. Assumption (B.4) is the same as (A.4) and common in high dimensional literature, as it ensures that the relevant coefficients are uniformly bounded away from 0 (and infinity). Assumption (B.5.a) is uniformly bounding the eigenvalues of $\Sigma_{1n}^d$ from $0$ and infinity, making it strictly positive definite. Assumption (B.5.b) is used in the asymptotic normality proof and is implied by (B.3.e) if the conditional expectations of the relevant covariates in the 2nd stage are uniformly bounded. (B.6) is used in the asymptotic normality proof.
Under the set of assumptions above, and the identification assumption that there are sufficiently many instruments for the endogenous covariates, I prove the following two results, the former for regarding the consistency of the 2nd stage BRIDGE estimator and the second regarding the model selection consistency and oracle efficiency that it achieves.
Note that (ref) utilizes all (A.1-A.4) as the oracle efficiency of the first stage coefficients is used in the proof to determine the rate of $\left\lVert\hat{\alpha}-\alpha\right\rVert$. Note that neither rate is affected by the actual number of the instruments but only on the maximum number of the relevant instruments among the endogenous regressors. It is not straightforward which of the rates is faster as, apart from $p_{xn}, k_{xn}$ and the sample size, they depend on different quantities. That is, $h_n$ also depends on the maximum number of relevant covariates in the separate 1st stage equations, while $h_n'$ depends on the tuning parameter $\lambda_{xn}$. Interestingly, the latter ensures consistency and is used as a middle step to prove the former rate which is essential for model selection consistency.
The first part of this result states that the method described gives exact zero estimated coefficient values for the sub-vector of $\hat{\beta}_n$ that corresponds to the irrelevant covariates. The second part, written in a similar fashion as the result in Huang:2008, states that the estimated non-zero coefficients are oracle efficient- that is, they asymptotically have the same distribution as they would if they were ex ante known and being the only ones used in the model. The result is written for the linear combination of the corresponding estimated sub-vector and it is straightforward to rewrite it for the marginal distributions. For each $j=1,\dots,k_{xn}$, I can pick $\delta_n= e_j$ where $e_j$ the $k_{xn}\times 1$ unit vector with $0$ elements everywhere apart from the $j$-th position. Also, define $s_{dnj}^2=\sigma_{\varepsilon}^2e_j'(\Sigma_{1n}^{d})e_j$. Then, for each $\hat{\beta}_{1nj}$, it holds that $n^{1/2}s_{dnj}^{-1}(\hat{\beta}_{1nj}-\beta_{01j})\xrightarrow{d}N(0,1)$.
\setcounter{dispsection}{2} \setcounter{theorem}{0}
Keeping (A1-A4) and adding some more structure to the error term, the following list of assumptions is sufficient for model selection consistency and oracle efficiency under the BRIDGE penalty function.
(C.1.a) is the same as in the previous case, while (b) relaxes the standard assumption of normality that is usual in this problem and is used to prove model selection consistency. (C.2) follows \citeasnoun{Huang:2008} and limits the correlation between relevant and irrelevant variables and ensures that the relevant ones are sufficiently important in the model. This is also a crucial element to prove selection consistency. For the same proof, I use the rate restrictions in (C.3); it is interesting to see that there is no direct restriction on the total number of the instruments, or the total number of covariates. (C.4) is the same as for BRIDGE when $p<n$, just avoiding the restriction on the eigenvalues of the full Gram matrix. (C.5) and (C.6) are used for the oracle efficiency proof, with the former used for the Lindeberg condition. The latter includes rate restrictions to ensure consistency and oracle efficiency. Since I have ensured model selection, I can now use the previous results for $p<n$, but using $k_{xn}$ as the number of covariates. Thus, these rates look very similar with the ones used in the previous section for the same purpose, by substituting $k_{xn}$ for $p_{xn}$. The following two theorems state the results formally:
\setcounter{dispsection}{3} \setcounter{theorem}{0}
In this section, the same notation will hold as well as the normalization of the variables. I now define the adaptive LASSO problem. I assume that an initial estimator $\tilde\beta_n$ is available and I define the weights and the relevant loss function:
Also, let $\Sigma_{n12}^d=(\Sigma_{n21}^{d})' =n^{-1}D_{1n}'D_{2n}$, as well as $H_n=D_{1n}(D_{1n'}D_{1n})^{-1}D_{1n}'$. Furthermore, for a vector $\delta$, let $sgn(\delta)=(sgn(\delta_1),sgn(\delta_2),\dots)$, and $\hat{\delta}=_s\delta$ if and only if $sgn(\hat{\delta})=sgn(\delta)$.
Continuing to assume (A1-A4) and adding some more structure to the error term, the following list of assumptions is sufficient for model selection consistency and oracle efficiency under the Adaptive LASSO penalty function.
where for (D.2) I use the following definition from the online appendix of \citeasnoun{Huang:2008b}:
(D.1) is crucial for the proof of model selection consistency, as it allows for the use of the relevant maximal inequalities. Assumption (D.2) states the strong prerequisite for the success of adaptive lasso; having a good initial estimator to calculate the weights. This can be weakened slightly by assuming a weaker consistency notion and adding a partial orthogonality condition in the same style as in BRIDGE. The exact condition is discussed formally in \citeasnoun{Huang:2008b}. (D.3) bounds the parameters that correspond to the relevant regressors away from 0 and away from drifting to infinity. The former is a stronger requirement as it does not allow for positive but arbitrarily small values. (D.4) includes all relevant rates conditions necessary. While the first three are used in multiple occasions, (d) and (e) are paired with the corresponding maximal inequalities that explore the $\psi_d$ norm of a given sub-gaussian term. (D.5) bounds the eigenvalues of the Gram matrix and the relevant Gram matrix. It is interesting to note that for BRIDGE, I could avoid assuming anything about the eigenvalues of the full matrix. Finally, as in the previous cases, the instruments are bounded in probability.
With this set of assumptions, I prove the following two results. It is noteworthy that the objective in the first theorem is slightly stronger than model selection consistency. Instead of only identifying the true zeros, the method also identifies the sign of the relevant covariates. The second theorem is the standard oracle result from the previous two cases.
\setcounter{dispsection}{3} \setcounter{theorem}{0}
In the two step problem, the researcher should first pick a tuning parameter $\gamma$ for each stage. In this exposition, in order to simplify the notation, I assumed the same $\gamma$ in both cases throughout the theoretical results. This comes without loss of generality analytically, but the implications of different choices will be exposed in a later subsection. After computing $\hat{\alpha}_n$ with BRIDGE, she should construct $\hat{D}_n$ and, plugging it in the 2nd stage, compute $\hat{\beta}_n$ with BRIDGE. Then, she can construct a consistent estimator of $s_{dn}$, which can be done in two parts: (1) estimate the Gram matrix of the estimated covariates with non-zero coefficients using $\hat{\alpha}_n$ which overlaps with $\hat{\Sigma}_{dn}$ as defined above, and (2) estimate $\sigma_{\varepsilon}^2$ using the residuals from the 2nd stage. The, $\hat{s}_{dnj}^{-1}$ gives an approximate standard error for $\hat{\beta}_j$ for each $j$.
Given that BRIDGE is often challenging to compute as a non-convex problem, I provide a short discussion on the algorithm I use for the simulations, proposed by Huang:2010:
In the simple case of the linear model, having only one-step estimators, we define the LS loss function as \(Q_n(\beta)=\frac{1}{2}\sum_{i=1}^{n}(Y_i-X_i\beta)^2,\) and the bridge penalized function as \(L_n(\beta)= Q_n(\beta)+\lambda \sum_{j=1}^{p}|\beta_j|^{\gamma},\ for\ 0<\gamma<1.\) I will follow the algorithm proposed by \citeasnoun{Huang:2010}, as an improvement to the algorithm of \citeasnoun{Huang:2008}. The former is more efficient as it does not require any approximation. To achieve this, the authors use the fact that the maximizer $\hat{\beta}_n$ of the following function: \[S_n(\beta,\theta)=Q_n(\beta)+\sum_{j=1}^{p}\theta_j^{1-\frac{1}{\gamma}}|\beta_j|+\tau_n\sum_{j=1}^{p}\theta_j\] is equal to the optimizer of $L_n(\beta)$, and vice versa, under the constraint that $\hat{\theta}_j\geq 0,\ for\ j=1,\dots, p$ and the tuning parameter $\lambda= \tau_n^{1-\gamma}\gamma^{\gamma}(1-\gamma)^{\gamma-1},\ where\ \tau_n$ is the penalty parameter of $S_n(\beta,\theta)$ Huang:2009.
Based on this result, they propose a simple iterative algorithm: First, I pick the initial value for $\beta^{(0)}$ to be the the corresponding LASSO estimate of the problem. Then, I compute the $\theta_j^{(s)}$ that appear in $S_n(\beta,\theta)$ as: \[\theta_j^{(s)}=\left(\frac{1-\gamma}{\tau_n\gamma}\right)^{\gamma}|\beta_j^{s-1}|^{\gamma},\ j=1,\dots,p.\]Using that value, I minimize \[Q_n(\beta)+\sum_{j=1}^{p}(\theta_j^{(s)})^{1-\frac{1}{\gamma}}|\beta_j|\] over $\beta$ to compute the new value of the vector $\beta^{(s)}$. Then, I repeat the last two steps until we achieve convergence, which is always attainable as $S_n(\beta,\theta)$ decreases in each step. To set the penalty parameters $\tau_n$, I use the functional form of $\lambda$ that the equivalence result requires. To compute $\lambda$ itself, I use cross-validation, splitting the sample in 5 folds.
Applying the general algorithm to my case, I compute the 1st stage coefficients, then estimate $\hat{D}_n$, I plug it in the 2nd stage and estimate $\hat{\beta}_n$ using the same algorithm.
Regarding the simulated DGP, starting from the first-stage parameters, i.e., the matrix $\alpha$, I draw its columns $\alpha_j,\text{ for}\ j=1,\dots, p_x$ as follows. For the column elements $\alpha_{jh}$, where $h=1,\dots,p_{zn}$, I pick $\alpha_{jh} \in (-5,5)$ only for $k_{xn}$ entries and $\alpha_{jh}=0$ otherwise. I add a small normal random noise to all the entries. The sparsity structure of the 1st stage coefficient matrix has to satisfy that the number of overall relevant instruments ($\sum_{j=1}^{p_{zn}}k_{znj}$) is at least as high as the number of the relevant covariates and, to avoid further identification issues, the same subset of instruments is only used for at most as many covariates as its cardinality. To keep the example clean, there is one relevant instrument per 2nd stage covariate with different degrees of relevance defined by the corresponding elements of $\alpha$. The 1st stage $\gamma$ for BRIDGE is set out to be $0.1$.
The instrumental variables are drawn from a multivariate normal, with a flexible correlation pattern among the relevant ones. That is, $Z_i\sim\mathcal{N}_{p_{zn}}(0,\Sigma_{\bm{z}}).$ We define the variance-covariance matrix $\Sigma_{z}$ to have a Toeplitz structure, i.e. $\Sigma_z|_{jk}=\rho^{|j-k|},\ j,k\in[k_{xn}],\ \rho=0.7$ for the principal submatrix referring to the relevant instruments and the identity matrix for the rest: \[\Sigma_{z} =
\]
Further, I pick a joint normal distribution for the i.i.d. error terms in the first and the second stage, $(\{v_{ij}\}_{j=1}^{p_x},u_i)$, and I choose the variance covaraince matrix to allow for non-trivial correlation terms between the two. That is, $(u_i, \{v_{ij}\}_{j=1}^{p_{xn}})\sim \mathcal{N}_{1+p_{xn}}(0,\Sigma_{uv})$, where \[\Sigma_{uv}=
\] with $\sigma_u=\sigma_{\bm{v}}=\sqrt{0.5}$ and for $\sigma_{uv}= (\sigma_{uv^1},\dots, \sigma_{uv^{p_x}})$, some elements are equal to 0.4 and some 0.15. These choices give us a positive-definite matrix. Lastly, The second stage vector of non-zero coefficients is also constructed from numbers in $(-5,5)$. Now, one can draw $X_{ij}$ as $X_{ij}=Z_i'\alpha_j+v_{ij}$ and the outcome variable as $Y_i=X_i'=\beta+u_i$.
The first set of simulations presented below is using very small sample sizes to examine the elementary dynamics between the number of covariates and $n$, as well as the role of $\gamma$ in the problem. The second set is examining a larger sample size $n=1000$ with increasing number of covariates from $100$ to $1000$, as a more realistic scenario on where the method might be useful.
The first table has data from 200 simulations with information on the average root mean squared error (RMSE), the median RMSE and the number of variables that the method picks to be non-zero. I set $p_{xn}=p_{zn}=30$, $k_{xn}=k_{zn}=6$. I present results for OLS as a baseline, even though it is expected to work poorly in this environment, for LASSO, as it is a popular solution once the researcher faces a high dimensional environment and BRIDGE for three values of $\gamma$, to examine its significance. I pick three sample sizes to observe whether the method performs well on the “difficult" case that the sample size is equal to the number of variables, on the case that the sample size is equal to the number of covariates plus the number of instruments, and on the case that the sample size is double this number and thus, an “easier" case.
The setup described in the previous subsection is an environment where LASSO is model selection consistent and oracle efficient Gold:2020. Even in this environment that does not take advantage of BRIDGE's validity under more flexible distributions, the two methods are fairly competitive. Looking at table (ref) below, one can see that OLS is performing very poorly even with slightly bigger sample size than the number of the covariates. Even though the mean MSE of LASSO and BRIDGE sharply drop once the sample size is slightly larger, OLS RMSE is marginally reduced. For a sample size of 30, the mean of the former two is significantly higher than the median, implying large outliers, but the two of them balance once I increase the sample size. LASSO picks, on average, a larger model for a small sample size, but once $n>30$, LASSO and BRIDGE have competitive values for RMSE and pick similar model sizes, close to the truth.
Let $\hat{S}$ to be the set of indices of the 2nd stage covariates chosen by each method to be non-zero. Table (ref) presents two values: the probability that $\hat{S}$ contains the set of indices of the non-zero variables of the true model and the probability that the two quantities overlap. As expected, OLS always chooses a very big model, failing to estimate any of the true zeros. For $n=30$, LASSO picks larger models than it should and BRIDGE smaller, even though BRIDGE has a higher probability of picking the exact correct one. As the sample size grows, LASSO and BRIDGE always pick a set that contains the true one and, with a persistently very high probability, they pick exactly the true one.
Further, it is interesting to examine further the role of $\gamma$ in the 2nd stage. In table (ref), I present results for $n=60$, $p_{xn}=30, p_{zn}=30$ and $k_{xn}=k_{znj}=6$ for 13 values of $\gamma$. For the first stage, the 1st stage $\gamma$ is set to be $0.1$ across all simulations. Note that as 2nd stage $\gamma$ grows, the penalty function approaches LASSO. It can be observed that even in the edge cases, close to $0$ and $1$, the median performance remains reasonably good, and the variables collected are the same with the true model. The only steep change is the increased mean RMSE for $\gamma$ around $0.6-0.75$, which may be incidental. There is no strong evidence that the choice of $\gamma$ significantly affects the performance of BRIDGE.
Regarding the set of simulations with a higher sample size ($n=1000$), there are data of up to 200 simulations, and I report the same summary statistics. I increase the number of the instruments to $100$ in every case, where the relevant ones remain 6. The number of covariates is either 100, 500, or 1000, creating approximately the same dynamics as the “smaller $n$" case.
By table (ref), one can observe that OLS significantly under-performs, even in the first case, with $p_{xn}=10\%$ of $n$.Even if picking the correct $0$s was not the researcher's objective, the mean RMSE is much higher than both LASSO and BRIDGE. BRIDGE seems to have many outliers in this case, increasing the mean RMSE and the average number of selected variables but it may be incidental, since it does very well on the next two cases. LASSO does worse in picking the correct model, but is competitive on the average error it makes.
A compatible story can be said about table (ref). OLS and LASSO always pick a model which contains the true one, but almost never exactly the true one. BRIDGE on the other hand picks exactly the correct model almost always. In the first case, potential outliers that involve a much larger model than the correct drop the $P(\hat S= True)$ to 73%, but for the other two cases it seems to be performing very well. It is important to iterate here that this is still an environment favorable to LASSO, since the errors are still following a normal distribution.
Forthcoming sets of simulations will include errors that follow a sub-Gaussian but not normal distribution as well as PDFs with fatter tails, as well as the comparative results with adaptive LASSO.
\FloatBarrier
\setcounter{dispsection}{3} \setcounter{theorem}{0}
This paper is an exposition of how BRIDGE and adaptive LASSO can be used in a very popular environment, the linear model under the presence of endogeneity, when the researcher faces high dimensionality in both stages of the process. Facing a larger class of problems compared to the usual analysis, i.e. replacing the assumption of normal with sub-gaussian errors, I prove that both methods are model selection consistent and oracle efficient even when the number of covariates exceeds the sample size. For BRIDGE, I also prove that if the former is lower than the latter, the same properties hold without sub-gaussian errors. BRIDGE requires a slightly weaker set of assumptions to have the desirable properties, while adaptive LASSO is expected to be much faster computationally, so the methods are competitive on different fronts and the one that is recommended depends on the researcher's resources.