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.
68,701 characters · 12 sections · 72 citation commands
Treatment Effects Inference with High-Dimensional Instruments and Control Variables
\doublespacing
Estimation of causal and treatment effects has provided a valuable method of statistical analysis and understanding of policy variable effects. This is especially true for program evaluation studies in economics, finance, and statistics, where these methods help to analyze how treatments or social programs affect the outcome variables of interest. In addition, with the increasing availability of richer data sets, the use of high-dimension methods in treatment effects models has become widespread in the literature (see, e.g., BelloniChernozhukovHansen14, chernozhukov2015post, chernozhukov2018debiased, and angrist2022machine).
Endogeneity of the variable of interest is a very pervasive issue in empirical applications, and one of the most important challenges faced by researchers. The instrumental variables (IV) analysis is commonly used in practice to compute treatment effects for endogenous regressors (see, e.g., ImbensRubin15 and references therein). This is notably relevant for estimating causal effects when instruments are exogenous conditional on observables. Recently, the use of high-dimensional linear IV models has proven to be crucial for correct and accurate estimation and inference in these studies (see, e.g., belloni2012sparse, hansen2014instrumental, zhu2018sparse, CattaneoJanssonMa19, gold2020inference, and QiuTaoZhou21). This scenario is prevalent when the number of variables is very large, or when interactions between variables should be considered. It has been common in the literature to separately allow for high-dimensional IV or high-dimensional control variables.\footnote{There is also an existing literature considering the case of many IV. In this case, it is common that the number of IV grows together with the sample size, at the same rate or slower, but it is not allowed to be larger than the sample size (see, e.g., hausman2012instrumental, ChaoSwansonWoutersen23, and mikusheva2020inference). In the high-dimension IV literature, the number of instruments is allowed to be larger than the sample size.} Ignoring the high dimension may result in biased inferences and misinterpretation of causal effects.
A typical solution to the high-dimension problem is to assume variable selection under the (approximate) sparsity, which effectively imposes that virtually all control coefficients are (very near) zero, and only a limited number of variables have substantial magnitude effects on the response. Many penalized-based regression methods have been proposed to obtain consistent estimators and confidence intervals for scalar parameters under the approximate sparsity assumption, including literature on selecting among high-dimension controls (see, e.g., zhang2014confidence, farrell2015robust, chernozhukov2015post, belloni2015uniform, and zhu2018sparse) and also in the high-dimension IV context (see, e.g., belloni2012sparse, and belloni2016post). Although the sparsity assumption is common for estimation and inference in many models in the presence of high-dimensional instruments and controls, it is unverifiable and may often be violated in empirical applications. Once sparsity is violated, the estimator may be severely biased, and the corresponding confidence interval for the low dimensional parameter is inconsistent. Thus, considering inference under a non-sparsity assumption is of practical importance.
Non-sparse structures of control variables have been studied in the literature. Recently, cha2023inference provide a general framework for inference in high-dimensional regression models without assuming sparsity, incorporating multiple estimation steps to account for complex dependence structures. A portion of it weakens the sparsity condition (see, e.g., he2000parameters, buhlmann2003boosting, belloni2014inference, hansen2014instrumental, kozbur2017testing, solvsten2020robust, kozbur2020analysis, mikusheva2020inference, ing2020model, and li2021LimitControl). A compelling strategy to address the high-dimensional variable problem in the absence of sparsity is to employ a two-step ridge regression. hansen2014instrumental use this framework and suggest a two-stage ridge jackknife instrumental variables estimator (RJIVE) under high-dimensional instruments. However, also considering high-dimensional control variables in this framework can add further challenges, since correlations between the two stages cannot be estimated/eliminated by conventional procedures.
In this paper, we contribute to the literature on IV models by proposing estimation and inference of treatment effects in the presence of both high-dimensional instruments and high-dimensional control variables with a non-sparse structure. We first show that, in the presence of high-dimensional control variables, unfortunately, the existing RJIVE estimator suffers from inconsistency, which is induced by the use of the jackknife strategy together with the ridge regression when controlling for the high-dimensional exogenous controls in the first stage. Hence, we introduce a novel two-step estimator that employs a ridge regularization to the instrumental variables coupled with a data-splitting strategy, as well as a ridge projection matrix to partially out exogenous regressors.\footnote{angrist1995split use a sample splitting strategy to solve the two-stage least squares bias. belloni2012sparse employ sample splitting in the context of sparse models.} In contrast to existing approaches, which often involve iterative model selection, boosting techniques, or multi-step procedures, our proposed method uses a two-step estimator, where the first stage uses ridge regression coupled with data-splitting, and the second stage uses a ridge-style projection matrix with simple least squares regression. This framework offers a simpler and computationally efficient alternative, and specifically addresses challenges in models with both high-dimensional IV and exogenous controls, introducing new features related to instrument-control correlations that our method is designed to handle.
The practical implementation of the estimator is simple. In particular, in the first step, one uses the first partition of the data to compute a ridge regression of the endogenous variable on the set of instruments. Then, given these estimates, one uses the second partition of the data to calculate the fitted values. In the second step, to compute the treatment effects estimator, one uses a projection matrix with a ridge penalty term to partial out the high-dimensional exogenous covariates together with a simple ordinary least squares estimator regressing fitted values from the first step on the outcome variable. Mild sufficient conditions are provided for this two-step estimator to have the desired asymptotic properties, namely, consistency and asymptotic normality. We highlight that this new framework is valid with or without sparsity assumptions. By avoiding variable selection procedures, the proposed methods are more robust than sparsity-based methods for causal inference with both high-dimensional IV and controls.
We develop statistical inference procedures for the proposed methods. It has been known that when the covariate dimension is large relative to the sample size, many classical methods become invalid for inference. liu2020estimation develop inference methods based on ridge regression that can be applied to both low- and high-dimensional models. We extend their results to our two-stage treatment effects estimator. In particular, we employ a ridge-type variance estimation technique, which offers improved efficiency compared to traditional methods, and is easy to compute. We formally establish the consistency of the variance estimator. Finally, since the parameter of interest is finite-dimensional and the weak limit of the estimator is standard, one is able to employ standard inference procedures.
We conduct Monte Carlo simulations to evaluate the finite sample performance of the proposed methods.\footnote{The practical implementation of the proposed method is available in the R package “HDRRTreat”.} We investigate models with and without the sparsity condition. Numerical simulation results document evidence that the new estimator outperforms existing sparsity-based as well as non-sparsity-based approaches across a variety of settings. The proposed estimator is approximately unbiased, and the empirical coverage of the confidence interval is close to the nominal one. Overall, simulation results are in line with the theory, and provide favorable numerical evidence for the suggested methods.
Finally, we provide an empirical application to illustrate the methods proposed in this paper. We revisit the classic example of angrist1991does and estimate the causal effect of schooling on earnings by addressing potential endogeneity through the use of high-dimension instrumental variables and high-dimension covariates. Results document evidence of imprecise estimates, with large confidence intervals.
The remainder of the paper is organized as follows. Section (ref) presents the model and the parameter of interest. In Section (ref), we discuss the two-step estimation procedure. The limiting statistical properties of the estimator are presented in Section (ref). Section (ref) presents inference procedures, and Section (ref) collects numerical Monte Carlo results. An empirical illustration of the methods is provided in Section (ref). Finally, Section (ref) concludes. All proofs are collected in the Appendix.
The main objective of this paper is to produce causal inference in a model that allows for endogeneity of the variable of interest, as well as high-dimensional instruments and high-dimensional exogenous covariates. The study of causal effects is an important tool for practitioners. We consider the following model
where $y_i$ is the response variable, $d_i $ is a scalar endogenous treatment variable, and $\alpha$ is the main treatment effects parameter of interest. The vectors $X_i$ and $Z_{1i}$ contain the exogenous controls and the instruments variables, respectively. We consider $X_i$ and $Z_{i}:=[Z_{1i},X_{i}]$ to have dimensions $p_x$ and $p_z$, respectively, with $p_{z} = p_{x} + p_{z_{1}}$. The function $\Upsilon_i$ is the unknown optimal instrument, and $f(\cdot)$ is the unknown function of exogenous controls. The endogeneity in the above system comes from the correlation of errors $\epsilon_i$ and $v_i$.
It is common in applications to specify and estimate simple linear models. We follow the literature and assume a linear specification for the controls $f(\cdot)$ and the optimal instruments $\Upsilon(\cdot)$. Although both of these functions might be unknown, and potentially complicated, the linear specification can be satisfied with a series expansion of the observable variables.\footnote{For example, one could have $X_{i} =\{p_{kK}(\zeta_{i})\}_{k=1}^{K}$ and $Z_{i} =\{p_{mM}(\upsilon_{i})\}_{m=1}^{M}$ for some set of basis functions $\{p_{jJ}(\cdot)\}$ such as orthogonal polynomials or splines formed from, respectively, a set covariates $\zeta_{i}$, and a set of `fundamental' instruments $\upsilon_{i}$.} Hence, we consider the following model
with $Z_i=[Z_{1i},X_{i}]$. We assume that $\textnormal{E}[v_{i}|Z_{i}]=0$ and that $\textnormal{E}[\epsilon_{i}|Z_{i}]=0$ for all $i$. Further, we assume that $\textnormal{Var}(Z_{i})\neq 0$ such that $Z$ is a valid instrument for $X$.
In this work, we allow for both the number of exogenous covariates and available instruments to be larger than the sample size. In particular, we assume that $p_x/n \rightarrow \tau_x$ and $p_z/n = (p_{z_{1}}+p_x)/n = (p_{z_{1}}/n+p_x/n)=\tau_{z_{1}}+\tau_x \rightarrow \tau_z$, where $\tau_{x}>0$ and $\tau_{z}>0$. The parameter of interest, $\alpha$, is presented here as a scalar for simplicity. The methods and results discussed can be readily generalized to accommodate a vector-valued $\alpha$, preserving the same theoretical properties and interpretations. Next, we introduce an estimator for the causal parameter of interest $\alpha$.
This section proposes a two-step estimator for the treatment effects parameter of interest. Estimating the parameters of the model in equations (ref)--(ref) in a high-dimensional non-sparse setting presents several challenges. First, it is computationally demanding to handle large numbers of instruments and exogenous control variables concurrently. Second, while sparsity generally facilitates the derivation of statistical properties for a penalized estimator, such as Lasso (see, e.g., chernozhukov2015post), it comes at the cost of complicated variance term derivations. Finally, and most crucially, simply using the ridge jackknife instrumental variables estimator (RJIVE) --- that is often used in the high-dimensional instrumental variable context (see, e.g., hansen2014instrumental) --- together with standard partialling out the exogenous regressors induces bias.
Before we describe the details of our proposed method, to fix ideas and build intuition, we discuss reasons why a simple extension of the existing two-stage RJIVE method would be inadequate for estimating a model with endogeneity and high-dimensional regressors and instruments.
Let us consider a simple potential extension of the RJIVE estimation procedure for the model in equations (ref)--(ref) above, which is inspired by a ridge regularization extension of the JIVE idea to the $n>p$ case. In the first-stage regression, the optimal instrument $\widehat{d}$ is estimated with a ridge regularized jackknife method as follows,
where $\lambda$ is the regularization parameter, and the subsamples $d_{-i}$ and $Z_{-i}=[Z_{-1i}, X_{-i}]$ consist of all but the $i$-th data point, and $I$ is a $p_{z} \times p_{z}$ identity matrix. This jackknife strategy introduced by angrist1999jackknife is able to eliminate the many instrument biases for an independent sample. In particular, it has been shown that, under the exogeneity of the instrument, and a small number of covariates, the JIVE strategy does not suffer from the many instrument biases, because of the following orthogonality
The last equality in (ref) is zero because both $\textnormal{E}\left[v_{-i} \epsilon_i \mid Z_{-i}\right]=0$ and $\textnormal{E}\left[Z_{-i}\gamma_z \epsilon_i \mid Z_{-i}\right]=0$. The former conditional expectation is zero due to $\epsilon_{i}$ being independent of $v_{j}$ if $j\neq i$, and the latter is zero due to the validity of the instrument.
Unfortunately, this jackknife strategy cannot accommodate scenarios with an additional large number of covariates, since equation (ref) will not be valid. In such high-dimensional regressor setting, a simple partialling out of exogenous variables is not feasible, and the resulting bias can no longer be disregarded. To see this, consider the residuals from equation (ref) after partialling out covariates using ridge regression:
where $A_{n}=X(X^{\top}X+\eta I)^{-1} X^{\top}$. The matrix $A_n$ denotes a feasible ridge projection matrix, regularized with the penalized term $\eta I$ (see, e.g., hansen2014instrumental for further discussion on the properties of this matrix).
Now, recall the optimal instrument in equation (ref), by taking the expectation of term $\widehat{\epsilon}_i \widehat{d}_i$, we obtain
It is important to note that the term $\textnormal{E}\left[v_{-i} \widehat{\epsilon}_i \mid Z_{-i}\right]$, in the above expression, is no longer equal to zero due to the fact that $\widehat{\epsilon}_i$ contains information from the entire sample because the matrix $A_n$ has information about every $i$, including the observation $-i$. Hence, $\widehat{\epsilon}_i$ is related to $v_{-i}$.
Consequently, the lack of orthogonality in (ref) induces inconsistency of this JIVE estimator. Intuitively, the inconsistency arises because the JIVE estimator relies on the assumption of independence between the error term and the constructed instrument, an assumption that is violated when the error term incorporates information from the entire sample. Thus, the JIVE strategy is not designed to accommodate the case with both high-dimensional instrumental variables and high-dimensional control variables.
To address the shortcomings of simultaneously accounting for a large number of instruments and control variables, we propose a novel two-step ridge regression approach for estimation under a non-sparsity assumption. To handle the high-dimension instrument and controls concurrently, we make use of two tools that are usually available in the literature. For estimation, we incorporate a data-splitting method in conjunction with a regularized (semi-)projection matrix for partialling out the exogenous regressors.
The two-step ridge regression estimation procedure is as follows. In the first step, we start by splitting the data set into two parts randomly. Denote the partitions of the data $\{1,\cdots,n\}$ into two subsets $\mathcal{S}_1 = \{1,\cdots,n_1\}$, $\mathcal{S}_2 = \{1,\cdots,n_2\}$ and $n=n_1+n_2$. Then, regarding the first-stage regression in equation (ref), we use the first partition of the data, and consider the following ridge regression estimation for the parameter $\gamma_{z}$:
where $\| \cdot \|_{2}$ is the Euclidean norm, and $\eta_{z}$ is a tuning parameter that we will discuss below. Thus, the estimator can be written as follows:
where $Z^{\left(1\right)}$ and $d^{\left(1\right)}$ denote the instruments and endogenous variables, respectively, with sub-samples indexed by $\mathcal{S}_1$. Now, using the second partition of the sample, together with the above estimates, we construct the fitted values as,
where $Z^{\left(2\right)}$ denotes the instruments for the sub-sample indexed by $\mathcal{S}_s$.
In the second step, we estimate the parameter $\alpha$ in equation (ref) by performing a simple ordinary least squares procedure after partialling out regressors as follows:
where
It is important to note that the notion of partialling out above uses the matrix $A_{n_{2}}$, which is a regularized (semi-)projection matrix with $\eta_{x}$ being a tuning parameter that we discuss further below.
We remark that, unlike hansen2014instrumental, our proposal incorporates a penalized term --- through the $A_{n_{2}}$ projection matrix --- in the second stage of estimation, rather than directly using the estimated values from the first stage. This approach enhances the robustness and stability of our estimator by regularizing the second stage, addressing high-dimensional issues in both stages. The two-step estimator in equation (ref) has two important features. First, the sample-splitting procedure allows for consistent estimation of the parameter of interest and is crucial for inference purposes. Technically, the procedure in equations (ref) and (ref) guarantees the independence of two leading terms in the linear representation of the estimator in the first step of estimation --- one term is the projection of the endogenous regressor on the set of instruments, and the other is the projection of the outcome variable on the set of exogenous controls. Hence, the sample splitting ensures that the limiting distribution of the estimator is well-behaved. Intuitively, this strategy allows for correct inference because the randomness inherent in the estimated treatment variable is independent of the dataset used for statistical inference. The data split strategy can also resolve the many instrument biases in the presence of controls due to the independent nature of the estimation and inference dataset.
The second feature is the partialling out of the high-dimensional control variables by the ridge projection matrix, $A_{n_{2}}$. As discussed previously, the matrix $A_{n_{2}}$ is a regularized (semi-) projection matrix, and the practical implementation of the estimator requires the selection of the penalty parameters $\eta_{x}$ and $\eta_{z}$. We follow the usual approach in the ridge regression literature (see, e.g., friedman2010regularization, and liu2020estimation) and set $\eta_{z}=c_{z} \max_{1 \leq j \leq p} \big| z_{(j)}^{\top} d \big| / (n_{1} p_{z})$, and $\eta_{x}=c_{x} \max_{1 \leq j \leq p} \big| x_{(j)}^{\top} Y \big| / (n_{2} p_{x})$, where $c_{z}$ and $c_{x}$ are preset constants. Here, $Y = (y_{1}, \ldots, y_{n})^{\top}$, $x_{(j)} = (x_{1j}, \ldots, x_{nj})^{\top}$ is the $j$-th column of $X \ (j = 1, \ldots, p_{x})$, and $z_{(j)} = (z_{1j}$, \ldots, $z_{nj})^{\top}$ is the $j$-th column of $Z \ (j = 1, \ldots, p_{z})$.
The tuning parameters $\eta_{z}$ and $\eta_{x}$ are used in the first and second stages, in equations (ref) and (ref), respectively. Parameters $\eta_{z}$ and $\eta_{x}$ capture the scale of the data, ensuring that the regularization terms are appropriately scaled relative to the dimensions of the control variables and instrumental variables, respectively. Specifically, they are designed to balance the regularization based on the maximum influence of individual predictors on the response variable $Y$ and the dependent variable $d$. In practice, we use $\eta = \min\{\eta_{x}, \eta_{z}\}$ for regularization to achieve a balanced approach.
In our proposed estimator, the tuning parameters control the extent of shrinkage applied to the coefficients rather than directly determining the sparsity of the model. In contrast to the $L_1$ norm requirements in chernozhukov2015post, where the tuning parameter is used to achieve varying levels of sparsity or to eliminate predictors in model selection, our method focuses on regularizing the magnitude of the coefficients, allowing them to be small but non-zero. It is worth noting that in this context, by employing ridge-regularization together with sample splitting, we can effectively control the complexity of the model while maintaining the inclusion of all predictors, thereby striking a balance between model interpretability and predictive performance.
A related contribution is chernozhukov2018debiased, who develop a debiased machine learning approach combining orthogonalization with sample splitting to deliver valid inference in high-dimensional settings. Our method differs in that it emphasizes shrinkage-based regularization rather than debiasing, aiming at stable estimation while retaining all predictors, and is thus, although related, more directly comparable to ridge-type procedures.
The next section shows that the two-step procedure discussed here has desirable statistical properties. In particular, we establish the consistency and asymptotic normality of the estimator. We will also discuss practical inference.
This section presents the statistical properties of the two-step estimator discussed in the previous section. We first discuss the assumptions, and then establish the consistency and asymptotic normality of the estimator.
Notation is as follows: We have a dataset $\{X_{i},Z_{i},y_{i},d_{i}\}_{i=1}^{n}$ of size $n$, where $d_{i}$ is the treatment, $X_{i}$ are control variables, $Z_{i}$ are high-dimensional instrument variables, and $y_{i}$ is the outcome variable. We randomly split the dataset into two parts, $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$, with size $n=n_1+n_2$. We denote the dataset $\mathcal{S}_{1}$ for the first partition $\{X_{i}^{\left(1\right)},Z_{i}^{\left(1\right)},y_{i}^{\left(1\right)},d_{i}^{\left(1\right)}\}_{i=1}^{n_1}$, and $\mathcal{S}_{2}$ for the second partition $\{X_{i}^{(2)},Z_{i}^{(2)},y_{i}^{(2)},d_{i}^{(2)}\}_{i=1}^{n_2}$. Finally, our goal is to estimate and conduct inference for the treatment effect parameter $\alpha$.
We start by stating and discussing the assumptions to establish the statistical properties of the estimator presented above. Let $\{X_{i},Z_{i},y_{i},d_{i}\}_{i=1}^{n}$ be independent and identically distributed (i.i.d.) observations. We also consider the following conditions.
Assumption (ref) describes the asymptotic representation property of optimal instruments. This condition decomposes the optimal instruments into two components: one residing in the column space of confounding variables, and another orthogonal to the confounding variables. This ensures stability and control over these terms as the sample size increases. Assumption (ref) is extensively employed and documented in the literature on many instruments (see, e.g., hausman2012instrumental and hansen2014instrumental).
Assumption (ref) ensures that the signal retained after regularization constitutes a non-negligible fraction of the original signal, in line with the literature on instrumental variables chernozhukov2015post,hansen2014instrumental. The boundedness of both the largest and smallest eigenvalues guarantees that the design matrices $X$ and $Z$ are well-conditioned, thereby ensuring numerical stability and invertibility of relevant sample moment matrices. The restrictions on $\|\gamma_x\|$ and $\|\gamma_z\|$ are bounds of the magnitude of the coefficients helping to prevent overfitting and to maintain stability of the regression estimates. The eigenvalue restrictions are imposed in a probabilistic sense, holding with probability approaching one as $n \to \infty$, and are satisfied under sub-Gaussian designs provided that $p_x/n \to \tau_x < 1$ and $p_z/n \to \tau_z < 1$ vershynin2010introduction.
In Assumption (ref), $B$ denotes the correlation structure between the instruments $Z$ and the covariates $X$. Bounding $\|B\|$ serves two purposes: it ensures that the projection of $Z$ onto the column space of $X$ remains well-behaved (avoiding arbitrarily large projections), avoiding unbounded projections, and excessively large coefficient estimates that may lead to numerical instability and obscure the model’s economic interpretation. The conditions on the covariance matrix $\Sigma_W$ guarantee that the component of $Z$ orthogonal to $X$ contains sufficient independent variation. Specifically, bounding the largest and smallest eigenvalues guarantees that the residualized instruments neither exhibit excessive variance nor approach degeneracy. This well-conditioned residual component makes the instruments informative beyond what is explained by $X$, preserving identification strength and promoting the numerical stability of subsequent estimators.
Notice that Assumption (ref) allows for both the number of instruments and covariates to grow with the sample size. In particular, these ratios could even diverge to infinity. This assumption is an extension of that in hansen2014instrumental.
Assumption (ref) imposes a homoskedasticity condition. This constraint on the variance structure simplifies the model's error term properties, facilitating more straightforward and easily implementable inference procedures. We leave extensions to the heteroskedastic case to future research.
In this section, we establish the asymptotic properties of the general two-step estimator described in equation (ref) above. In particular, we establish the consistency and the asymptotic normality. We will discuss details on inference procedures in the next section.
Theorem (ref) shows the consistency of the proposed estimator, a desired property for most estimators. The variable $Q$, as presented in the theorem statement, represents a specific transformation of the instrumental variables matrix $Z$. This transformation involves regularization by $\eta_{z}$ followed by projection using $I-A_{n}$. We assume that the projection $Q$ behaves in a predictable manner, and the impact of complex interactions diminishes. This assumption is natural since
is also a quadratic form. Note that this assumption pertains to the population $Z$, but in the proof below, only the subset $Z^{(1)}$ is required. Therefore, this assumption is sufficient.
Here, the tuning parameter $\eta_x$ controls the extent of shrinkage applied to the coefficients rather than directly determining the sparsity of the model. In contrast to the $L_1$ norm requirements in chernozhukov2015post, where the tuning parameter is used to achieve varying levels of sparsity or to eliminate predictors, our method focuses on regulating the magnitude of the coefficients, allowing them to be small but non-zero. It is worth noting that in this context, we are imposing a stronger assumption of normality on the regression coefficients. By leveraging this sufficient assumption, we can effectively control the complexity of the model while maintaining the inclusion of all predictors, thereby striking a balance between model interpretability and predictive performance.
The following result establishes the asymptotic normality of the two-step estimator, which is crucial for conducting statistical inference. The limiting distribution of the estimator is shown to be a standard normal, which provides the foundation for the inference procedures discussed in the subsequent section. In order to construct confidence intervals and perform hypothesis tests, an easy-to-use consistent estimator for the asymptotic variance term, denoted by $\sigma_{\alpha}$, is required. The derivation and properties of such a variance estimator will be addressed in the next section.
We close this section emphasizing the importance of the sample splitting for the properties of the estimator. When considering a high-dimensional structure, the dependence between the first and second stages can lead to biased and inconsistent estimates. The data splitting ensures that the estimates from the first stage are independent of the second stage errors, thereby preserving the validity of the asymptotic properties. Next, we develop inference procedures.
In this section, we turn our attention to inference procedures. Given the results on asymptotic normality of the estimator presented in Theorem (ref), inference for the treatment effects parameter is simple. The only remaining ingredient is a consistent estimator of the variance, $\sigma^2_{\alpha}$.
We construct an estimator for the variance by employing the $RidgeVar$ technique for inference (liu2020estimation). The proposed approach involves applying the ridge projection matrix to both the control variables (numerator) and the instrumental variables (denominator). In comparison to existing literature, our $RidgeVar$ method offers improved efficiency for conducting inference.
From the asymptotic normality in Theorem (ref), we wish to estimate the parameter $\sigma^2_{\alpha}$. Now, using the sample-splitting together with the homoskedasticity condition in Assumption (ref), we have that $\textnormal{E}_{X}\left[\epsilon^{(2)}\epsilon^{(2)\top}\right]=\sigma_{\epsilon}^{2}$ and $\textnormal{E}_{X}\left[v^{\left(1\right)}v^{\left(1\right)\top}\right]=\sigma_{v}^{2}$. Thus, the following sandwich-type variance estimator can be derived,
The approximation in the last line of the above display arises from the assumption that $(I-A_{n_{2}})^2 \approx (I-A_{n_{2}})$. This simplification is justified when $A_{n_{2}}$ is a projection matrix or when its eigenvalues are close to 0 or 1, leading to negligible differences between its square and itself.
Therefore, we propose to use the following estimator:
where the estimator for the variance of the error term is given by
with $P_{n_{2}}=S_{n_{2}}\left(S_{n_{2}}^{\top}S_{n_{2}}+\eta_{s}I\right)^{-1}S_{n_{2}}^{\top}$, and $S_{n_{2}}=\left[d^{(2)},X^{(2)}\right]$.
The practical estimation of the variance of the error term in equation (ref) is simple, and consequently, the estimation of the variance in equation (ref) is straightforward. This method relies on the methods developed in liu2020estimation. They suggest a consistent estimator of variance based on ridge regression and random matrix theory, which is valid under both low- and high-dimensional models. We extend their results to cases where both the regressors and instruments are of high dimension.
The following theorem formally establishes the consistency of the variance estimator $\widehat{\sigma}_{\alpha}^{2}$ in (ref).
Given such an estimator, it is possible to formulate a wide variety of tests and construct confidence intervals. General hypotheses on the vector $\alpha$ can be accommodated by standard Wald-type tests. These statistics and their associated limiting theory provide a natural foundation for testing the null hypothesis $H_{0}: \alpha = r$ when $r$ is known. Thus, the result in Theorem (ref) can be used to test the null hypothesis $H_{0}$.
Thus, the Wald statistic can be constructed as
where $\widehat{\sigma}_{\alpha}^{2}$ is given in equation (ref).
This section presents a Monte Carlo simulation study designed to assess the finite sample performance of the proposed estimation and inference methods. Overall, this simulation study aims to provide insights into the small sample behavior of the new procedures under different data generating process (DGP) specifications, taking into account the endogeneity of the variable of interest, and the presence of high dimension for both instrumental variables and exogenous regressors.
We first introduce the DGP of variables used in the simulations. The outcome variable $y_i$ is generated according to the following model:
where $d_{i}$ is the endogenous variable of interest, the vector $X_i$ contains the completely exogenous regressors, and $\epsilon_i$ is the error term. In all simulations, the true value of the parameter of interest, $\alpha$, is set to be 1. The endogeneity of the variable $d_i$ is generated according to the equation:
where $Z_i$ is a vector of instrumental variables, and $v_i$ is the error term. We will discuss the specific details about the choice of the parameter values for $\gamma_x$ and $\gamma_z$ below.
The error terms $\epsilon_i$ and $v_i$ are jointly normally distributed with a mean vector of zeros and a covariance matrix given by:
This covariance structure allows for correlation between the error terms, with a correlation coefficient $\sigma = 0.6$.
This design framework allows us to investigate the relative performance of various estimators under different sparsity patterns while controlling for the overall signal strength. The combination of these structures with correlation patterns (EC and AR(1)), described below, provides a comprehensive setting for examining the behavior of high-dimensional instrumental variable estimators under varying degrees of sparsity and correlation.
In the simulations that follow, we introduce two distinct DGPs, one designed for non-sparse settings, and another for sparse. We evaluate the proposed two-step ridge regression (TSRR) estimator for the parameter of interest discussed in equation (ref), as well as its corresponding variance estimator in (ref). For comparison, we also include results for high-dimensional methods (DML) presented in chernozhukov2015post and belloni2017program; and the regularized jackknife instrumental variable estimator (RJIVE) suggested in hansen2014instrumental.\footnote{The DML is implemented with the R package hdi, and the RJIVE is implemented by ourselves. All the codes are available upon request.} We present results for the following empirical statistics: bias (for point estimates), bias of the variance, mean squared error (MSE), the probability of the true parameter value being inside the confidence interval (P(cover)) for a 95% nominal coverage, the length of the confidence interval, and the computation time. Finally, in the simulations, the number of replications is set to $500$.
As discussed previously, $\eta = (\eta_x,\eta_z)$, the penalty parameter of TSRR, is chosen adaptively as the minimum of two data-driven values: one based on the correlation between controls and the outcome, i.e., $\eta_{x}=c_{x} \max_{1 \leq j \leq p} \big| x_{(j)}^{\top} Y \big| / (n_{2} p_{x})$; and another based on the correlation between controls and the endogenous variable, i.e., $\eta_{z}=c_{z} \max_{1 \leq j \leq p} \big| z_{(j)}^{\top} d \big| / (n_{1} p_{z})$, and , where $c_{z}$ and $c_{x}$ are preset constants. Both values are scaled by the sample size, dimension, and a tuning parameter $c_{\cdot} \in \{0.1,1\}$. This dual-equation approach ensures appropriate penalization for both the structural and first-stage relationships.
We evaluate the finite-sample performance of our proposed estimator through Monte Carlo simulations. Unless otherwise noted, the sample size is set to $n=500$, with $p_{x}=700$ control variables and $p_{z}=500$ instruments. All simulations are repeated 1000 times. We examine several data-generating processes (DGPs) that vary in both the strength and structure of the instruments as well as the correlation patterns among covariates. For each design, we report the average bias, root mean squared error (RMSE), and empirical coverage of the nominal 95% confidence intervals.
To examine the robustness of our estimator, we consider several distinct data-generating processes (DGPs) that differ in instrument strength and dependence structure, including an “all-weak” design, a “cutoff” structure, and a “sparse” case with an autoregressive AR(1) correlation pattern. The detailed specifications of these scenarios are given below.
The first DGP, which we label the “non-sparse” case, features an equi-correlation (EC) structure, where all variables exhibit a common, non-negligible level of pairwise correlation. Specifically, the correlation matrix, $\Sigma$, follows an equi-correlation structure given by $\Sigma_{i\neq j}^{EC}=\rho$ and $\Sigma_{i,i}^{EC}=1$. This structure ensures that all variables maintain a consistent correlation level $\rho$. We examine the case with equi-correlation structure $\rho = 0.04$.
The coefficient vector for the control variables is: $\gamma_x = (m,\ldots,m,c\xi_1,\ldots,c\xi_s,0,\ldots,0)^{\top}$, with 5 strong signal components of magnitude $m$, and 5% of the remaining coefficients are non-zero. These non-zero coefficients $\xi_j$ are independently drawn from $N(0,1)$ and scaled by constant $c$, which is calibrated to achieve a specified signal strength $\mu_x^2 \in \{300,600\}$. The scaling constant $c$ is chosen such that:
where $\mu_x^2$ represents the signal strength following chernozhukov2015post, and $\Sigma_x$ denotes the correlation matrix for control variables.
Moreover, in this setting, we have a moderately high-dimensional control variable space with $p_{x}=700$ covariates and a higher-dimensional instrument space of $p_{z}=500$ variables. The instrument coefficients maintain an “all-weak" structure, where $\gamma_z = (c\xi_1,\ldots,c\xi_s,0,\ldots,0)^{\top}$ with 50% of the coefficients being non-zero and scaled to achieve a signal strength of $\mu_z^2 = 600$, which measures the overall strength of the instruments. In another “cut-off" structure, the instrumental variables' coefficients follow a cutoff structure where exactly 70% of the instruments are non-zero, while the remaining 30% are zero coefficients. Specifically, we set the first 70% of the instruments to have equal coefficients (initially set to 1) and the remaining instruments to have zero coefficients. These coefficients are then rescaled by a constant to achieve the same concentration parameter ($\mu_z^2 = 600$) as the “all-weak" setting. It imposes a deterministic pattern where the relevance of instruments follows a strict ordering.
The results for the first DGP, the less sparse controls, are collected in Table (ref). The upper part of Table (ref), Panel A, presents results for the case where the instruments follow a cutoff structure with 70% active instruments, representing a scenario with strong identification. The middle panel, Panel B, examines a setting with a constant correlation structure among covariates and instruments that follow an all-weak pattern where 30% of instruments are active but with smaller magnitudes. The lower panel, Panel C, investigates the case where $c_{x}=c_{z}= 1$, where $c_{\cdot}$ is a constant tuning parameter that determines the relative magnitude of regularization $\eta$. With $c_{\cdot}= 1$, the penalty parameters are directly proportional to the maximum correlations in the data, scaled by sample size and dimension. This serves as a reference point for evaluating how different choices of tuning parameters affect the regularization strength.
In this non-sparse case, our proposed estimator demonstrates superior performance compared to both alternative methods in all cases. Our approach yields very small empirical biases for both cases. For Panel A the point estimate bias is $-0.0184$, while the bias of variance is $0.0252$. In Panel B the bias for point estimate and variance are $0.0077$ and $0.0038$, respectively. In Panel C bias of point estimate and variance for TSRR are the smallest among all estimators. In addition, the TSRR estimator shows minimal bias and maintains excellent empirical coverage properties, about 95.4%, 95.8%, and 96.6%, for Panels A, B, and C, respectively. In contrast, the RJIVE estimator produces biased results. In Panel A the point estimates of RJIVE are severely biased, while in Panels B and C the point estimates have relatively smaller bias, but large bias in the variance. This lack of precision reflects on the empirical coverage of the confidence intervals, which are only about 86% for the nominal 95%. Finally, the DML estimator shows deteriorating performance with higher bias, and its coverage probabilities (88%, 86%, and 85%) fall notably below the nominal level.
The second DGP, the “sparse" case, maintains the same dimensionality as the first setting, with $p_{x}=700$ control variables and $p_{z}=500$ instruments. However, the correlation structure now follows an autoregressive AR(1) pattern, where $\Sigma_{i,j}^{AR(1)}=\rho^{|i-j|}$. This implies a sparser dependence structure, where variables are only correlated with their immediate neighbors, and the correlation strength decays exponentially with distance. The control variable coefficients in this case follow below: $\gamma_x = (c\xi_1,\ldots,c\xi_s,0,\ldots,0)^{\top}$, with 30% of the remaining coefficients are non-zero but relatively weaker in magnitude. The instrument coefficient structure remains consistent with the first DGP, maintaining the “cutoff" specification with the same signal strength calibration. Specifically, in an AR(1) structure, each variable is only correlated with its immediate neighbors, and the strength of the correlation decreases exponentially as the distance between the variables increases.
Table (ref) collects the results for the second specification for the case with sparse (AR(1)) correlation structure ($\rho = 0.5$) and sparse control variables, where the first five coefficients have strong signals (magnitude of 2). There are three configuration in this table. Panel A examines the case where the control variables exhibit moderate sparsity, with 30% of the control coefficients being non-zero, alongside 5 strong signals of magnitude 2. Panel B presents a more sparse setting where only 5% of the control coefficients are non-zero, while maintaining the same 5 strong signals of magnitude 2. In both panels, the instruments follow a cutoff structure where the first 70% of instruments have non-zero coefficients. This design allows us to examine how the sparsity of control variables affects estimation while maintaining consistent instrument strength and correlation patterns. Panel C investigates the case where $c_{x}=c_{z} = 1$, where the penalization term in the objective function is neither amplified nor dampened, providing a natural benchmark for assessing the impact of regularization strength on the estimator's performance. In this setting, our proposed estimator demonstrates comparable performance to DML estimator, with both methods achieving empirical coverages close to the desired 95% nominal coverage probability. While DML shows lower MSE for all panels, both methods effectively estimate the true parameter value of 1, with TSRR presenting relatively small biases. Notably, both methods substantially outperform RJIVE, which exhibits considerably higher MSE for all three panels, and wider length for confidence intervals.
These results highlight two key strengths of our proposed estimator. First, it maintains competitive performance with state-of-the-art methods, such as DML in sparse settings. Second, and more importantly, it demonstrates superior robustness across different data generating processes, maintaining high accuracy and appropriate coverage even in non-sparse settings where alternative methods may struggle. This robustness to both sparse and non-sparse scenarios, combined with its computational efficiency, makes our method particularly attractive for empirical applications where the true data structure is unknown.
This section provides an empirical application to illustrate the proposed two-step ridge regression (TSRR) methods in this paper. Following angrist1991does and hansen2014instrumental, we revisit the classic example in the many-instrument literature, which also contains high-dimension covariates. This example focuses on estimating the causal effect of schooling on earnings by addressing the potential endogeneity of schooling through the use of quarters of birth and interactions as instrumental variables. While maintaining a similar methodological framework, our analysis differs in several key aspects from the original study. We use a subset of the angrist1991does dataset, focusing on individuals under 40 years of age, and employ a sample-splitting approach, where we use 5% of the total sample. This dataset exploits an unusual natural experiment to estimate the impact of compulsory schooling laws in the United States. The experiment stems from the fact that children born in different months of the year start school at different ages, while compulsory schooling laws generally require students to remain in school until their sixteenth or seventeenth birthday. In effect, the interaction of school entry requirements and compulsory schooling laws compel students born in certain months to attend school longer than students born in other months. Because one's birthday is unlikely to be correlated with personal attributes other than age at school entry, the quarter of birth generates exogenous variation in education that can be used to estimate the impact of compulsory schooling on education and earnings. We refer readers to angrist1991does for further details on the experiment and data.
Our instrumental variables framework follows the original linear specification for the structural equation:
with a standard linear first-stage regression
where $w_i$ is the weekly wage, $s_i$ represents years of education, $X_i$ is a vector of control variables, and $Z_i$ is the vector of instruments. Moreover, the error terms satisfy:
The control variables, $X_i$, include year-of-birth dummies, state-of-birth dummies, and their interactions, resulting in a dimension of $p_X = 510$. Our set of instruments, $Z_i$, comprises three quarter-of-birth dummies and their interactions with the control variables, yielding a total of $p_Z=1527$ instruments. We present results for a smaller configuration of the number of instruments as well. In addition, for comparison, we present results for the RJIVE estimator.\footnote{The RJIVE estimator for $\beta_1$ takes the form: $\widehat{\beta}_{RJIVE} = (D'P_{Z\lambda}D)^{-1}D'P_{Z\lambda}y$, where $P_{Z\lambda}$ is the regularized projection matrix with ridge penalty $\lambda$.}
The estimation results are collected in Table (ref), where we present the proposed two-step ridge regression (TSRR), and for completeness we present the RJIVE results from hansen2014instrumental. The main results for the TSRR model show that our estimates are positive across all specifications, with the point estimates about $0.17$. Most notably, when using the full set of instruments ($1527$ IVs), we obtain a statistically significant estimate of $0.18$ with a 95% confidence interval of $[0.01, 0.36]$. Economically, this result indicates a positive return to schooling, with an additional year of education associated with approximately 18% higher weekly wages. The estimates from the smaller instrument set (with 180 IVs) yield similar point estimates but with wider confidence intervals that include zero. Statistically, this result is very intuitive, since the proposed procedures are designed to accommodate a large number of instruments and exogenous covariates. In comparison, the RJIVE estimates show invariant results across different instrument specifications. With 180 instruments (including quarter-of-birth interactions with year-of-birth and state-of-birth main effects), they obtained an RJIVE estimate of 0.11 (SE: 0.02). In their most comprehensive specification with 1,527 instruments (including all possible interactions), the RJIVE estimate was 0.1067 (SE: 0.02).
Overall, our TSRR approach yields results that are somewhat larger in magnitude than those from the RJIVE method in the original paper, particularly when using the full set of instruments. Our estimate of 0.18 is approximately 69% higher than the RJIVE estimate of 0.1067, though there is some overlap in the confidence intervals. The statistical significance of our full-instrument model suggests that the TSRR method can effectively handle high-dimensional settings while producing economically meaningful estimates that are broadly consistent with the literature on returns to schooling.
This paper proposes novel ridge regularization-based methods for estimating treatment effects and conducting inference in the presence of both high-dimensional instrumental variables and high-dimensional control variables. An advantage of these methods is that they are valid with or without sparsity assumptions. To mitigate high-dimensional instrument biases, we suggest a data splitting strategy for constructing an optimal instrument. The data splitting ensures that the estimates from the first stage are independent of the second stage errors, thereby preserving the validity of the asymptotic properties of the estimator. We establish the statistical properties of the estimator, namely, consistency and asymptotic normality. We also propose an estimator for the variance and establish its consistency.
Moreover, we study the finite sample performance of the proposed methods using numerical simulations. Results document evidence that the proposed estimator has good finite sample properties, outperforming existing sparsity-based approaches across a variety of settings. The overall numerical results provide evidence of the resilience and adaptability of our method to changes in the sparsity structure, highlighting its robustness compared to existing approaches. Finally, we provide an empirical application to estimating the causal effect of schooling on earnings by addressing potential endogeneity through the use of high-dimension instrumental variables and high-dimension covariates.