Robust Inference with High-Dimensional Instruments

Robust Inference with High-Dimensional Instruments

Introduction

Various recent surveys in leading economics journals suggest that weak instruments remain important concerns for empirical practice. For instance, I.Andrews-Stock-Sun(2019) survey 230 instrumental variable (IV) regressions from 17 papers published in the American Economic Review (AER). They find that many of the first-stage F-statistics (and non-homoskedastic generalizations) are in a range that raises such concerns, and virtually all of these papers report at least one first-stage F with a value smaller than 10. Similarly, in lee2021's (lee2021) survey of 123 AER articles involving IV regressions, 105 out of 847 specifications have first-stage Fs smaller than 10. Moreover, many IV applications involve a large number of instruments. For example, in their seminal paper, Angrist-Krueger(1991) study the effect of schooling on wages by interacting three base instruments (dummies for the quarter of birth) with state and year of birth, resulting in 180 instruments. Hansen-Hausman-Newey(2008) show that using the 180 instruments gives tighter confidence intervals than using the base instruments even after adjusting for the effect of many instruments. In addition, as pointed out by MS22, in empirical papers that employ the “judge design" (e.g., see maestas2013, sampat2019, and dobbie2018), the number of instruments (the number of judges) is typically proportional to the sample size, and the famous Fama-MacBeth two-pass regression in empirical asset pricing (e.g., see fama1973, shanken1992, and anatolyev2022) is equivalent to IV estimation with the number of instruments proportional to the number of assets. Similarly, belloni2012 consider an IV application involving more than one hundred instruments for the study of the effect of judicial eminent domain decisions on economic outcomes. carrasco2015 used many instruments in the estimation of the elasticity of intertemporal substitution in consumption. Furthermore, as pointed out by Goldsmith(2020), the shift-share or Bartik instrument (e.g., see bartik1991 and blanchard1992), which has been widely applied in many fields such as labor, public, development, macroeconomics, international trade, and finance, can be considered as a particular way of combining many instruments. For example, in the canonical setting of estimating the labor supply elasticity, the corresponding number of instruments is equal to the number of industries, which is also typically proportional to the sample size. Similar patterns of many IVs occur in wind-direction IVs, granular IVs, local average treatment effect estimation, and Mendelian randomization.\footnote{E.g., see deryugina2019mortality, bondy2020crime, gabaix2024, blandhol2022, boot2024, sloczynski2024should, davey2003, and davies2015.}

The problem of inference under many instruments has been well studied in the literature (e.g., see the literature review below). However, to our knowledge, there still does not exist a valid inference method that is robust to the high dimensionality of instruments and general error dependence, such as network and spatial dependence, simultaneously. To fill this gap, we propose in this paper a novel procedure for inference under high-dimensional IVs. Our test statistic takes a self-normalized form, and we establish the asymptotic validity of our test under general error dependence structure, by using random matrix theory (RMT). Additionally, we derive the power property of our test under many instruments. Monte Carlo simulations are conducted to assess the numerical performance of the test, confirming good size control and satisfactory testing power across a range of various error dependence structures.

Literature Review. The contributions in the present paper is related to the large literature on many (weak) instruments.\footnote{See, for example, Kunitomo1980, morimune1983, Bekker(1994), donald2001, chamberlain2004, Chao-Swanson(2005), stock2005, han2006, Andrews-Stock(2007), Hansen-Hausman-Newey(2008), Newey-Windmeijer(2009), anderson2010, kuersteiner2010, anatolyev2011, belloni2011, okui2011, belloni2012, carrasco2012, Chao(2012), Haus2012, K13, hansen2014, carrasco2015, Wang_Kaffo_2016, kolesar2018, EK18, solvsten2020, CNT23, boot2024, among others.} In the context of many instruments and heteroskedasticity, Chao(2012) and Haus2012 provide standard errors for Wald-type inferences that are based on JIVE and jackknifed versions of the limited information maximum likelihood (LIML) and Fuller(1977)'s (Fuller(1977)) estimators (HLIM and HFUL). These estimators are more robust to many instruments than the commonly used two-stage least squares (TSLS) estimator because they can correct the bias caused by the high dimension of IVs.\footnote{Specifically, the rate of growth of the concentration parameter, which measure the overall instrument strength, is denoted as $\mu_n^2$. JIVE, HLIM, and HFUL remain consistent with heteroskedastic errors even when instrument weakness is such that $\mu_n^2$ is slower than the number of instruments $K$, provided that $\mu_n^2/\sqrt{K} \rightarrow \infty$ as the number of observations $n \rightarrow \infty$ (Chao(2012), Haus2012). In contrast, TSLS is less robust to instrument weakness as it is shown to be consistent only under homoskedasticity if $\mu_n^2/K \rightarrow \infty$ Chao-Swanson(2005). In simulations derived from the data in Angrist-Krueger(1991), which is representative of empirical labor studies with many instrument concerns, Angrist-Frandsen2022 show that such bias-corrected estimators outperform the TSLS that is based on the instruments selected by the least absolute shrinkage and selection operator (LASSO) introduced in belloni2012 or the random forest-fitted first stage introduced in athey2019.} Furthermore, carrasco2012, carrasco2015, carrasco2016efficient, hansen2014, and carrasco2017 proposed regularization approaches for two-stage least squares, limited information maximum likelihood, and jackknife IV (Angrist(1999)) estimators. Under many weak moment asymptotics, Newey-Windmeijer(2009) provide new variance estimators for the jackknife GMM and the class of generalized empirical likelihood (GEL) estimators, which includes the continuous updating estimator (CUE) and EL estimator as special cases.\footnote{In the linear heteroskedastic IV model, consistency and asymptotic normality of CUE require $m^2/n \rightarrow 0$ and $m^3/n \rightarrow 0$, respectively, where $m$ and $n$ denote the number of moment conditions and the sample size (e.g., see p.689 of Newey-Windmeijer(2009)). Such conditions are needed to simultaneously control the estimation error for all the elements of the heteroskedasticity consistent weighting matrix. Somewhat stronger rate conditions are required for other GEL estimators.}

However, the Wald-type inference methods are invalid under weak identification, which occurs when the concentration parameter remains bounded as the sample size increases to infinity. In this case, all the estimators mentioned earlier become inconsistent, and there is no consistent test for the structural parameter of interest.\footnote{E.g., see Section 3 of MS22.} For weak-identification-robust inference based on the classical AR test, Andrews-Stock(2007) show its validity under many instruments, but their IV model is homoskedastic and requires the number of instruments to diverge slower than the cube root of the sample size ($K^3/n \rightarrow 0$). Newey-Windmeijer(2009) proposed a GMM-AR test under many (weak) moment conditions but imposed the same rate condition on $K$. anatolyev2011 constructed a modified AR test that allows the number of instruments to be proportional to the sample size but requires homoskedastic errors, and Kaffo-Wang(2017) proposed a bootstrap version of anatolyev2011's test. Furthermore, carrasco2016 first proposed a ridge-regularized AR test that allows for $K$ being larger than $n$ with homoskedastic errors. Bun-Farbmacher-Poldermans(2020) compared the centered and uncentered GMM-AR test and identified a missing degrees-of-freedom correction when $K/n \rightarrow 0$.

Recently, crudu2021 and MS22 proposed jackknifed versions of the AR test in a model with many instruments and general heteroskedasticity. DKM24 developed a ridge-regularized version of the jackknife AR test. tuvaandorj2024 established the validity of a permutation AR test under heteroskedasticity and diverging $K$, requiring $K^3/n \rightarrow 0$. boot-ligtenberg(2023) developed a dimension-robust AR test based on continuous updating but relied on an invariance assumption. Furthermore, LWZ2024dimension proposed a dimension-agnostic bootstrap AR test by deriving strong approximation results for both test statistic and its bootstrap version. Their bootstrap AR test is valid with heteroskedastic errors uniformly across a broad asymptotic regime for $K$, spanning from fixed to diverging faster than the sample size. In addition to the AR tests, Matsushita2020 propose a jackknife LM test. Lim2024 consider a linear combination of jackknife AR, jackknife LM, and orthogonalized jackknife LM tests and find that the resulting conditional linear combination (CLC) test has good power properties in a variety of scenarios. N23 proposed a dimension-robust version of Kleibergen(2002)'s K test, and his method relies on a sparse $\ell_1$-regularized estimation of $\rho(Z_i)$, the conditional correlation between the endogenous variable and the outcome error. Yap24 proposed a leave-three-out variance estimator, in the same spirit as that proposed by AS23, for inference under heterogeneous treatment effects. The above methods are robust to weak identification, many instruments, and heteroskedastic errors. However, the literature on inference procedures that are robust to more general error dependence structure remains sparse. ligtenberg2023 proposed a leave-cluster-out AR test that is robust to cluster dependence structure under many weak instruments. \nocite{Wang-Zhang2024}

The rest of this paper is organized as follows. In Section (ref), we formally define the observed intermediary, which is used to formulate our proposed statistic. Then, we establish their limiting distributions with power theories. Section (ref) reports the finite sample performance of our proposed statistics across multiple scenarios through Monte Carlo simulations. The Appendices are provided separately as the supplementary document, including proof of theorems and relevant lemmas.

Model and Test Statistic

We consider the linear IV regression with a scalar outcome $y_i$, a scalar endogenous variable $x_i$, and a $K \times 1$ vector of instruments $z_i$ such that

Following the literature on many instruments, we assume that there are no (low-dimensional) controls included in our model as they can be partialled out from $(Y_i,X_i,Z_i)$. The source of endogeneity is caused by a correlation between $\varepsilon_{i}$ and $v_{i}$, which represent the errors associated with the structural and the first-stage equation, respectively. The two sequences of regression errors, $\left\{ \varepsilon_{i}\right\} _{i=1}^{N}$ and $\left\{ v_{i}\right\} _{i=1}^{N}$, follow a mean-zero random process whose functional forms are unknown and potentially heteroskedastic. The unknown parameters are $\beta\in\mathbb{R}$ and $\pi\in\mathbb{R}^{K\times1}$, where $K$ is comparable or exceeds its sample size.

We are interested in testing the null hypothesis

even when the identification for the structural parameter of interest $\beta$ may be weak. We focus on the model with a scalar endogenous variable for two reasons. First, in many empirical applications of IV regressions, there is only one endogenous variable (as can be seen from the surveys by Andrews-Stock-Sun(2019) and lee2021). Second, the asymptotic results derived in the paper extend directly to the general case of full-vector inference with multiple endogenous variables. Additionally, for the subvector inference, one may use a projection approach Dufour-Taamouti(2005) after implementing our test on the whole vector of endogenous variables.\footnote{Alternative subvector inference methods for IV regressions (e.g., see GKMC(2012), Andrews(2017), and GKM(2019), GKM(2021)) provide a power improvement over the projection approach under fixed $K$. However, whether they can be applied to the current setting is unclear. Also, Wang-Doko(2018) and Wang(2020) show that bootstrap tests based on the standard subvector AR statistic may not be robust to weak identification even under fixed $K$ and conditional homoskedasticity.}

To proceed, let us denote $y^*_{i}=y_{i}-x_{i}\beta_{0}.$ Let a singular value decomposition of the normalized instrument matrix be expressed as

where $Z:=\left(z_{1},\ldots,z_{N}\right)^{'}\in\mathbb{R}^{N\times K},$ $\lambda_{\ell}$ denotes the $\ell$-th largest non-zero eigenvalue of the matrix $S_{N}:=\frac{1}{N}Z^{'}Z,$ $q_{\ell}$ and $w_{\ell}$ denotes the associated orthonormal eigenvector of the $\ell$-th largest eigenvalue of $\underline{S}_{N}:=\frac{1}{N}ZZ^{'}$ and $S_{N}$, respectively.

To conduct the hypothesis testing in (ref), the associated test statistic proposed in this paper can be obtained by equating the sample analogue of

where $Y^{*}=Y-X\beta_{0}$, $Y=\left(y_{1},\ldots,y_{N}\right)^{'}$, $X=\left(x_{1},\ldots,x_{N}\right)^{'}$, and $\gamma_{j}^{2}:=\mathbb{E}\varepsilon_{j}^{2}$. Similar (high-dimensional) moment condition was studied in Feng2024 to overcome the curse of dimensionality of coefficients in linear regression models with a presence of (unconditional) heteroskedasticity and autocorrelation of an unknown form. After substituting $\gamma_{j}^{2}$ with its plug-in estimator under the null, namely $\frac{Y^{*'}Y^{*}}{N}$, the statistic of interest is obtained by evaluating its associated sample counterpart:

where $\underline{Y}:=Y^{*}/\left(Y^{*'}Y^{*}\right)^{1/2}\in\left(0,1\right)$. Then, with a proper convergence rate (shown in Theorem (ref)), our (oracle) statistic of interest is given by

Unless stated otherwise, we denote $\underline{a}:a/\left(a^{'}a\right)^{1/2}$ a self-normalized version of any column vector $a.$

To build intuition, suppose for the moment that the errors are homoskedastic, i.e., $\mathbb{E}(\varepsilon \varepsilon' |Z) = \gamma^2I_N$, where $\varepsilon = (\varepsilon_1, ..., \varepsilon_N)'$. By the law of iterated expectation,

under the null hypothesis, where $\ell$ is randomly chosen from the set $\{1, ..., \min\{K,N\}\}$. We note that $q'_{\ell}q_{s}$ and $w'_{\ell}w_s$ are equal to one when $\ell = s$ for all $\ell, s = 1, ..., \min\{K,N\}$ but zero otherwise. On the other hand, under the alternative $\beta \neq \beta_0$,

where $\Delta = \beta - \beta_0$.

For general error dependence, i.e., $\mathbb{E}(\varepsilon \varepsilon' | Z) = \Psi$, it is expected that under the null hypothesis,

where $\Psi^2_{ii}$ denotes the $i$-th diagonal entry of $\Psi$. The proposed test becomes robust against general error dependence with high-dimensional IVs if there exists, for some $i \in \{1, ..., N\}$, the sequence $\{ (q_{\ell}' \varepsilon \cdot \Psi_{ii}^{-1})^2\}_{\ell=1}^{\min\{K, N\}}$ converges to one in probability as $K$ and $N$ grow, where the quantity $q_{\ell}'\varepsilon \cdot \Psi_{ii}^{-1}$ is the inner product of an eigen basis generated by the Gram matrix $\underline{S}_N$ and the normalized regression errors. We establish the asymptotic validity of our inference procedure under the following regularity conditions.

With rank deficiency generally observed in sample covariance matrices, this regime is widely regarded as a standard high-dimensional paradigm in the literature Ledoit2002, Bai2007, Onatski2013. This paper introduces a more adaptable framework at a general level of intensity beyond that of Khan2024, where the number of covariates may exceed its sample size.

Assumption (ref) is recognized as a typical design in RMT literature Karoui2008, dobriban2018high, hastie2022surprises, Li2024, farbmacher2024revisiting, Zhang2025,which we employ to characterize correlations within instruments. The instruments, therefore, possess a weak form of cross-sectional dependence defined in Chudik2011. Qualitative instruments can also be included, but require standardization or variable transformation to ensure zero mean, unit variance, and a finite fourth moment (see, Theorem (ref)). Assumption (ref) allows us to apply RMT, which helps to generalize the error dependence structure, in the current setting with high-dimensional instruments, as can be seen in Assumption (ref) below.

The permitted structure of errors are general and flexible, accommodating complicated network or spatial dependence structure into account. Assumption (ref) ensures the mean independence of $\varepsilon$ and $v$ against instruments $z$. Specifically, $\mathbb{E}\left(\left.f_{i}^{\delta}\right|\varepsilon_{i},v_{i}\right)=\mathbb{E}\left(f_{i}^{\delta}\right)$, for $0<\delta\leq4$.

Define the self-normalized vector $u:=\left(u_{1},\ldots,u_{N}\right)$ where $u_{i}=\varepsilon_{i}/\left(\varepsilon^{'}\varepsilon\right)^{1/2}.$ The following theorem provides the null limiting distribution of $Q_{N}$, the oracle statistic, under high-dimensional IVs with a presence of (unconditional) heteroskedasticity in error processes.

Replacing $\text{tr}\left(\Sigma^{2}\right)$ with the consistent estimator proposed by Li2012, we have the following asymptotic result for the feasible statistic.

Let $\Delta = \beta - \beta_0$. The alternative with local departure is defined as

where $h = \Delta \cdot N^{1/2}/ (2 tr(\Sigma^2))^{1/5}$ denotes a (deterministic) level of departure. By construction, the quantity $\left(\frac{\left(2\text{tr}\left(\Sigma^{2}\right)\right)^{2/5}}{N}\right)^{1/2}$ implies a rate of signal strength in the coefficient of interest that the test can detect.

Several comments are in order. First, the power of the feasible statistic is analogous to that of the oracle one, since one merely replaces the unknown quantity $\text{tr}\left(\Sigma^{2}\right)$ with its consistent estimator. Namely $\left(\widehat{Q}_{N}-Q_{N}\right)+Q_{N}-\varpi_{N}=Q_{N}-\varpi_{N}+o_{\mathbb{P}}\left(1\right).$ Second, Condition (ii) can be relaxed to a finite fourth moment with a cost of technical complication.

Monte Carlo Simulations

In this part, we present the power performance of our proposed statistics across a range of data generating processes. The statistical software adopted for numerical results is MATLAB 2023a with the default seed. We further note that all parameters and functional forms postulated in the error processes are arbitrarily chosen. To design regimes of high dimensionality, we fix the sample size to $400$, then vary the number of instruments so that the data intensities ($K/n$) are $1/4,1/2,1,2$ and $3$. The numerical size (power) of statistics is calculated as the probability of null rejections over 1,000 replications, given the data is generated under the null (alternative) hypothesis. Throughout the analysis, an intercept is set to $2$ and a level of significance is set to 0.05.

The errors in the structural equation, i.e., $\varepsilon$, are generated by one of the three following processes: (i) the network dependent process, (ii) the spatial dependent process, and (iii) the multiplicative heteroskedasticity of a mixture of two i.i.d. processes. Specifically, for the network dependence process, we follow the simulation design of the network formation model in Section 5 of Kojevnikov2021, and set $\gamma=0.5$ in (5.1) of their model. The processes (ii) and (iii) were also implemented in Feng2024. For the case with spatial dependence, the errors follow the spatially autoregressive error process of order one studied by kim2011spatial, that is,

where $\rho_s=0.8$ and the innovation $e$ is independently generated by a standardized $t$-distribution with 5 degrees of freedom. In addition, $W=[c_iw_{ij}]$ is a $N \times N$ spatial weighting matrix. To specify such a matrix, we set all diagonal entries to zero and quantify the off-diagonal entry $w_{ij}$ to one if and only if $w_{ij}=a>0.5$, where $a$ is drawn from a uniform distribution over an interval $[0,1]$, otherwise zero. Denote $c_i$ a scalar that standardizes each row of $W$ to sum to unity. For the multiplicative heteroskedasticity and autocorrelation of a mixture of two i.i.d. processes, we let

where $\zeta_i = 1 + \left( a i/N \right)$ and $\omega_i = \sqrt{1-0.7^2} \eta_{1,i} + 0.7 \eta_{2,i}$, $\eta_{1,i}$ is generated by i.i.d. standardized $t$-distribution with 5 degrees of freedom, and $\eta_{2,i}$ is generated by i.i.d. standardized chi-square distribution with 6 degrees of freedom. As a result, the underlying process $\omega_i$ is considered both heavy-tailed and highly-skewed. The scalar $s_1:= 289a^2/300 + (1+ a/2)^2$ and $s_2:= 1.2(a+2)$ normalizes the random variable $\varepsilon_i$ to mean zero and variance one. Under this design, the errors are equi-correlated with a correlation approximately equal to $0.52$ for $i, j=1, ..., N$. The impact of the skedastic component is set to $a=10$.

For brevity, the three processes are abbreviated as {[}NET-E{]}, {[}SPA-E{]}, and {[}MUL-E{]}, respectively. The errors in the first-stage equation are correlated with the other given by $v_{i}=\rho\varepsilon_{i}+\sqrt{1-\rho^{2}}\eta_{3,i}$, where $\eta_3$ is independently generated by a standardized $t$-distribution with $5$ degrees of freedom (abbreviated as $t_{5}$ hereafter). The two processes are strongly correlated, taking values of $\rho\in\left\{ -.9,.5,.9\right\} .$

Being consistent with our asymptotic theory, the instruments are generated according to a latent factor model studied in Section 2.1 of bhattacharya2011sparse. Namely,

where $\eta_{4,i}\in\mathbb{R}^{3}$ and $f_{i}\in\mathbb{R}^{K}$ are independently drawn from $t_{5}$ and the standard normal distribution, respectively. $\Lambda$ is a $K\times3$ loading matrix such that $\Lambda^{'}\Lambda=\text{diag}\left(6,5,3\right)$. The population covariance matrix $\Sigma$ is the Toeplitz matrix with the $\left(i,j\right)$-th elements equaling to $0.7^{\left|i-j\right|}$. The strength of IVs is set to be weak such that $\pi^{'}\pi=1$ throughout the simulations.

The finite sample performance of the feasible statistic is analogous to that of the oracle one. To save space, we thus only report the results for the feasible test. We highlight several findings below. First, the simulation results illustrate that the empirical size of our test is close to the nominal level across different settings of the number of instruments ($K/N$), the degree of endogeneity ($\rho$), and the error dependence structure. Second, we note that the power of our test improves as the departure $h$ increases and attains unity, ensuring complete distributional separation between hypotheses, when $h$ exceeds 5. Third, in line with our asymptotic theory, the correlation between the two error processes influences the power of the test. A negative contemporary dependence enhances the power overall, as well as its magnitude.

Conclusion

In this paper, we propose a weak-identification-robust test for linear instrumental variable (IV) regressions with high-dimensional instruments, whose number is allowed to exceed the sample size. In addition, our test is robust to general error dependence, such as network dependence and spatial dependence. The test statistic takes a self-normalized form and the asymptotic validity of the test is established by using RMT. Simulation studies are conducted to assess the numerical performance of the test, confirming good size control and satisfactory testing power across a range of various error dependence structures. For future research agenda, we note that although our new test is robust to high-dimensional instruments under general error dependence, it is not robust to the case with a fixed number of instruments. On the other hand, LWZ2024dimension recently proposed a dimension-agnostic bootstrap AR test under many weak instruments and (independent) heteroskedastic errors by using strong approximation. For dimension-robust inference under the current setting, it may be interesting to explore a bootstrap-based inference procedure and combine it with RMT.