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.
87,804 characters · 15 sections · 95 citation commands
Identification-robust inference for the LATE with high-dimensional covariates
\interfootnotelinepenalty=10000
We propose an orthogonalized Anderson-Rubin (AR) test statistic with uniformly correct asymptotic size. The proposed method exhibits robustness against weak identification and high dimensionality in the LATE framework. Furthermore, we provide a practical guideline, including a step-by-step algorithm, for drawing inferences for the LATE with high-dimensional controls. This algorithm entails (1) inverting the proposed statistic to derive confidence intervals and (2) applying machine learning approaches to overcome the regularization bias and overfitting within the high-dimensional model. The objective is motivated by the persistent challenges of the weak-instrument problem in empirical research and the prevalence of rich data in the contemporary big data era.
In models where certain explanatory variables correlate with the error term, least squares estimators yield inconsistent coefficient estimates. To address this, instrumental variables are often employed, as they are uncorrelated with the error term but correlated with the endogenous explanatory variables. Nonetheless, if the correlation between the instruments and endogenous variables is weak, IV estimation becomes imprecise, resulting in unreliable tests and confidence intervals. This presents the weak-instrument problem, a notable concern in empirical practice.
Empirical researchers often seek to estimate the causal effects of endogenous regressors using instrumental variables regression. A prominent example is the influential study by angrist1991does, which uses quarter of birth as an instrument to estimate the returns to schooling. However, bound1995problems argue that Angrist and Krueger's results may be unreliable due to the weak correlation between one's quarter of birth and their education attainment. Moreover, the common practice of pretesting, with a rule-of-thumb F-statistic threshold of 10 proposed by staiger1997instrumental, is challenged by lee2022valid. In their paper, they introduce a novel critical value function and reveal that achieving a true 5 percent test with a critical value of 1.96 instead requires an F exceeding 104.7. Applying this criterion to their sample of 61 American Economic Review papers published between 2013 and 2019, they find that a quarter of the specifications initially presumed to be statistically significant are, in fact, insignificant.
imbens1994identification introduced a widely used framework for estimating the LATE. This parameter captures the treatment effect for the subgroup of compliers who take the treatment if and only if they are assigned to the treatment group. The use of instrumental variables to estimate LATE has received considerable attention in the literature. Within the LATE framework, weak identification emerges either when instruments correlate weakly with endogenous regressors or when the proportion of compliers is relatively small. We are particularly interested in exploring and addressing this weak identification issue within the LATE framework for several reasons. On the one hand, while compliers might be a minority, they often represent the population of critical interest to policymakers. Take the Vietnam-era draft lottery in angrist1990lifetime as an example. Even though compliers constituted a minority, approximately 0.10 to 0.16, their experiences provide valuable insights into the draft's direct consequences for individuals at the decision-making margin, thereby shedding light on the immediate effects of veteran status on civilian earnings. On the other hand, natural experiments often present challenges of weak identification, especially when researchers have no control over the size of the complier group, as highlighted by the instrument strength concerns in angrist1991does. Our objective is to yield reliable outcomes regardless of the proportion of individuals affected by the intervention, both theoretically and practically. Understanding this impact is crucial for expanding or scaling such interventions.
The weak-instrument literature has produced a range of econometric techniques for estimating and conducting inference on a structural parameter $\theta$ defined by moment conditions. In these models, one specifies a score $\psi(W;\theta)$ such that ${\mathrm{E}_{P}}[\psi(W;\theta_0)]=0$ at the true value $\theta_0$. We introduce an identification-robust test of the null hypothesis $H_0(\theta_0): {\mathrm{E}_{P}}[\psi(W;\theta_0)]=0$, which is equivalent to $H_0:\theta=\theta_0$. This problem has been addressed by, among others, anderson1949estimation, stock2000gmm, kleibergen2002pivotal, and andrews2016conditional. While these studies present methods tailored for inference about target parameters in the presence of weak identification, they overlook models equipped with high-dimensional covariates. Identification-robust testing procedures are particularly important in high-dimensional settings, as it is unclear how to pretest for weak identification in such cases. On the one hand, accounting for a large number of covariates can bolster the validity of the IV within the LATE framework. On the other hand, existing identification-robust methods face challenges in the presence of many covariates, especially due to size distortions when many covariates are included.
Our proposed method addresses these challenges by accommodating any arbitrary $N/p$ ratio, where $N$ is the sample size and $p$ is the number of covariates. This feature ensures that our method maintains correct size control regardless of the covariate dimensionality. Our technical innovation diverges from the typical approach of proposing a consistent LATE estimator. Instead, we introduce a stochastic process and its uniformly consistent estimator over the probability law under the null hypothesis. Building on this, we introduce a test statistic that exhibits uniformly correct asymptotic size. Our contribution generalizes prior research by developing an identification-robust statistic that employs machine learning techniques, enabling us to explore a broader set of controls than previously possible.
Our simulation results demonstrate that the proposed orthogonalized AR method performs reliably across identification strengths and covariate dimensionality, outperforming both the conventional AR test and the double/debiased machine learning (DML) method of CCDDHNR18. The conventional AR test exhibits severe size inflation as dimensionality increases and becomes infeasible when the number of covariates approaches or exceeds the sample size. The DML method, while valid for high-dimensional nuisance estimation, shows size deflation and over-coverage when instruments are weak or moderately weak, resulting in under-rejection. In contrast, the proposed method maintains empirical size close to the nominal level across sample sizes, identification strengths, and dimensionality. Under strong identification, its power is comparable to that of the DML method, while under weak identification it remains stable and correctly sized. These findings indicate that the proposed procedure combines the identification-robustness of the AR framework with regularized nuisance estimation, providing reliable inference across a wide range of data environments.
We further evaluate the proposed method in two empirical applications. The first, following hornung2015railroads, examines the effect of railroad access on city population growth in nineteenth-century Prussia. The second, following ambrus2020loss, investigates the long-term effect of cholera-related deaths on property values in nineteenth-century London. In both applications, the proposed method yields confidence intervals that are wider than those from the conventional AR test and the DML method, which is consistent with identification-robust inference when many covariates are included and instrument strength may be limited. Within our implementation, the resulting confidence intervals and inference outcomes are stable across regularization choices, including Lasso, Ridge, and Elastic Net, and across alternative control sets. Several effects that appear statistically significant under the AR test or the DML method are not significant under the proposed method. Taken together, these findings indicate that identification-robust inference with regularized nuisance estimation may alter statistical significance assessments in settings with many covariates and uncertain instrument strength.
This paper advances the well-established literature on weak identification by providing procedures for inference and the construction of confidence intervals for the LATE parameters in high-dimensional models. We develop a test statistic with uniformly correct asymptotic size. Furthermore, we provide a practical guideline, complete with a step-by-step algorithm, for drawing inferences and determining confidence intervals for the high-dimensional LATE using machine learning methods.
This paper contributes to the literature on weak identification and high-dimensional models by providing a test tailored for making inferences regarding the LATE in the presence of high-dimensional covariates.
Since staiger1997instrumental introduced the “local-to-zero” framework for weak instruments, a sequence of identification-robust tests has emerged.\footnote{See works by bound1995problems, kleibergen2002pivotal, andrews2006optimal, moreira2009tests, andrews2019identification, moreira2019optimal, mikusheva2022inference for various identification-robust inference methods developed over the past three decades. For detailed surveys on the weak identification literature, see stock2002survey, dufour2003identification, andrews2005inference, and andrews2019weak.} To test if the mean function equals zero at the true parameter value $\theta_0$, stock2000gmm pioneered the concept of weakly identified Generalized Method of Moments (GMM). They introduced the $S$ statistic as the quadratic form of the objective function, which is a generalized form of the AR test statistic as in anderson1949estimation and follows a $\chi^2$ asymptotic distribution under the null hypothesis. Later, kleibergen2005testing proposed the $K$ statistic, capitalizing on the asymptotic independence between the Jacobian estimator of the objective function and the sample average of the moment. moreira2003conditional sharpened power properties by constructing a conditional likelihood-ratio test that delivers exact size and locally most-powerful unbiased inference in linear IV models. andrews2016conditional further generalized the conditional approach by conditioning on the entire observed path of the sample moment process, obtaining procedures that remain valid regardless of identification strength and that achieve near-efficient power in both point-identified and set-identified designs. In the just-identified LATE setting with a single endogenous regressor and a single instrument, these statistics reduce algebraically to the AR test. We build on this equivalence and develop an orthogonalized AR procedure that retains identification-robust guarantees while accommodating high-dimensional nuisance learning through Neyman orthogonal scores and cross-fitting. This framework preserves the simplicity of AR, delivers size-correct inference under weak or strong first-stage relationships, and integrates naturally with modern control-selection methods introduced below.
Over the past decade, there has been a surge in the literature on machine learning-based econometric methods for high-dimensional models. belloni2015uniform advanced a Neyman orthogonal score for a Z-estimation framework in the presence of high-dimensional nuisance parameters. Subsequently, belloni2018uniformly constructed a confidence interval rooted in the Neyman orthogonality condition in the high-dimensional setting. In a series of contributions, Chernozhukov et al. (chernozhukov2013gaussian, chernozhukov2016empirical, chernozhukov2017central) established the Central Limit Theorem (CLT) for high-dimensional models using the Gaussian approximation approach. belloni2014high presented an overview of techniques for estimating and inferring in high-dimensional datasets. CCDDHNR18 introduced the DML methodology in the i.i.d. setting. They combined the Neyman orthogonality condition\footnote{We refer readers to pfanzagl1985contributions, bickel1993efficient, and newey1994asymptotic for the development of the Neyman orthogonal score.} and cross-fitting methods. Most recently, chernozhukov2022locally outlined a general construction for the doubly robust moment function, ensuring robustness against nonparametric or high-dimensional first steps. However, none of these papers on high-dimensional models consider weak identification issues.
This paper also relates to the literature on instrumental variables estimation of the LATE. imbens1994identification pioneered the introduction of a simple instrumental variables estimand for the average treatment effect among compliers. Motivated by angrist1991does, angrist1995two broadened the LATE framework to accommodate ordered treatments, such as years of schooling. A subsequent wave of research explores the incorporation of covariates into LATE estimation, including imbens1997bayesian, angrist2000interpretation, hirano2000assessing, yau2001inference, and abadie2003semiparametric, using either parametric or semiparametric estimation approaches. tan2006regression proposed a LATE estimator with robustness against the misspecification of either the propensity score model or the outcome regression model. hong2010semiparametric derived the semiparametric efficiency bounds for conditional and unconditional LATE. frolich2007nonparametric and ogburn2015doubly provided a fully nonparametric $\sqrt{N}$-consistent and efficient estimator for the LATE with confounding covariates. More recently, belloni2017program presented an efficient estimator alongside reliable confidence bands for the LATE with nonparametric/high-dimensional components, using the orthogonal moment condition and machine learning method. CCDDHNR18 incorporated their proposed DML method into the LATE framework, achieving an $\sqrt{N}$-consistent estimator for the LATE in the presence of high-dimensional covariates. angrist2022empirical underscored the importance of the LATE framework for causal inferences through empirical demonstrations.
Previous work on treatment effect estimation with many covariates included oprescu2019orthogonal, who proposed the orthogonal random forest to estimate heterogeneous treatment effects under unconfoundedness while allowing for nonparametric flexibility. qiu2021inference examined covariate-specific LATEs in high dimensional settings by imposing a parametric generalized linear model and applying variable selection and debiasing techniques to recover $\sqrt{N}$-consistency. sun2022high developed a doubly robust LATE estimator based on regularized, calibrated, and weighted $M$-estimators. boot2024inference studied the case of many and possibly weak instruments by introducing the saturated instrumental variable estimator, which debiased two-stage least squares in a fully interacted design. While their approach targeted a weighted, stratum-specific LATE, our method focused on the unconditional LATE under the assumption of strong monotonicity. Robustness was achieved through orthogonalization and the use of machine learning for estimating nuisance components, rather than relying on saturation. Although these methods enabled flexible, high-dimensional nuisance adjustment, their inferential guarantees depended on a well-identified first stage or correctly specified functional forms.
Empirical instrumental variables studies typically estimate the LATE conditional on a vector of pre-treatment covariates \(X\) and then average over the empirical distribution of \(X\) to obtain the unconditional LATE. Covariate adjustment relaxes the exogeneity condition to \(Z \perp\!\!\!\perp (Y(d), D(z)) \mid X\), which makes the independence and exclusion assumptions more plausible when the instrument is not fully randomized. It also absorbs heterogeneity, reducing residual variance and improving efficiency, and can enhance instrument sharpness. The findings of kennedy2020sharp demonstrated that a rich set of covariates predicting compliance can bolster causal identification even when instruments are weak by increasing instrument sharpness. These considerations explain the widespread practice of including many control variables in LATE analyses and motivate our high-dimensional framework, which delivers identification-robust inference uniformly in the number of controls \(p\). Accordingly, we target the unconditional LATE while allowing $p$ to grow with, or even exceed the sample size $N$. To the best of our knowledge, this paper is the first to provide identification-robust inference for the LATE under high-dimensional covariates, without imposing assumptions on instrument strength.
The rest of the paper is structured as follows. Section (ref) provides a practical guideline for implementing the proposed algorithm. Section (ref) presents the theoretical results. Section (ref) reports the Monte Carlo simulation findings. Section (ref) offers two empirical illustrations. Section (ref) concludes. The appendix contains all proofs of the theorems and lemmas, as well as an extension beyond the LATE framework to a general instrumental variables model.
In this section, we provide a brief overview of our proposed method without theories. This overview serves as a concise guideline in practice.
Consider the standard instrumental variable setup where the researcher has access to a dataset of N i.i.d. observations, represented as $\{W_i=(Y_i,D_i,Z_i,X_i')\}_{i=1}^N$. The outcome of interest for unit $i$ is denoted by $Y_i$. Let $D_i\in\{0,1\}$ be a binary indicator of the receipt of treatment for unit $i$. The instrumental variable $Z_i$ is also binary and can be interpreted, for example, as the offer of treatment. This instrument is randomly assigned conditional on the covariates. Let $X_i=(X_{i1},\cdots,X_{ip})'$ denote the $p$-dimensional vector of observed controls for unit $i$, where $X_{ij}$ represents the value of the $j$th covariate for unit $i$. Notably, the dimensionality $p$ can be substantially greater than the available sample size, $N$. For each unit \(i\), let \(D_i(z)\) denote the potential treatment received if the instrument takes the value \(Z_i=z\), and let \(Y_i(z,d)\) denote the potential outcome if \((Z_i,D_i)=(z,d)\). The causal parameter of interest is the LATE, $\theta_{\text{LATE}} = {\mathrm{E}_{P}}\!\left[Y_i(1,1)-Y_i(0,0)\mid D_i(1)>D_i(0)\right]$. Under the exclusion restriction, \(Y_i(z,d)=Y_i(d)\), the estimand simplifies to $\theta_{\text{LATE}} = {\mathrm{E}_{P}}\!\left[Y_i(1)-Y_i(0)\mid D_i(1)>D_i(0)\right].$ To streamline exposition, we adopt a binary representation for scalar $D$ and scalar $Z$. However, it is important to highlight that our framework can be extended to encompass broader contexts, including scenarios with ordered treatments like years of schooling, or when dealing with vector-valued $D$ and $Z$ as encountered in general instrumental variables models. Appendix (ref) develops the extension to general moment-restriction models.
Let $\{\mathcal{P}_N\}_N$ denote a sequence of probability laws associated with $\{W_i\}_i$. As the sample size $N$ grows, our analysis allows for an increasing dimensionality of $W_i$. Here, $P=P_N\in \mathcal{P}_N$ is defined with respect to a specific sample size $N$, and ${\mathrm{E}_{P}}$ stands for the expected value under the law $P$. For any set $B$, its complementary set is given by $B^c=\{1,\cdots,N\}\setminus B$, and $|B|$ represents the size or cardinality of $B$. We introduce the subsample expectation operator defined as $\mathbb{E}_B[\cdot]:=\frac{1}{|B|}\sum_{i\in B}[\cdot]$.
We model the random vector $W=(Y,D,Z,X')'$ as:
where $m_0$ is a function that maps the support of $(Z,X)$ to $(\varepsilon,1-\varepsilon)$, $g_0$ is a function that maps the support of $(Z,X)$ to $\mathbb{R}$, $p_0$ is a function that maps the support of $X$ to $(\varepsilon,1-\varepsilon)$ for some $\varepsilon\in (0,1/2)$, and $v, u, e$ are error terms. We do not impose any parametric assumptions\footnote{blandhol2022tsls demonstrated that 2SLS specifications can only have a LATE interpretation when controlling for “rich” covariates in a nonparametric manner.} on the form of the $m_0$, $g_0$, and $p_0$ functions here.
The population LATE proposed by tan2006regression is given by
The numerator includes the intent-to-treat (ITT) component along with augmentation terms that ensure robustness and orthogonality, while the denominator contains the compliance probability with analogous augmentation. The estimand in equation (ref) coincides with the causal LATE defined in Section (ref). Under the identification conditions for LATE, the augmentation terms have mean zero, implying that the numerator simplifies to the ITT effect and the denominator to the compliance probability. Let us define the compliance probability as $B_N={\mathrm{E}_{P}}[m_0(1,X)-m_0(0,X)]$. We retain the augmentation to construct a Neyman orthogonal score, which preserves the estimand but facilitates robust inference.
The standard normal distribution of the LATE estimator can be derived using the delta method, which linearizes the LATE estimator with respect to the estimators of the numerator and denominator in equation (ref). In line with the weak-instrument literature, we model weak identification by allowing the denominator to shrink so that the concentration parameter remains bounded as stated in Remark (ref). This corresponds to a small compliance probability. In such instances, the standard normal approximation fails because the LATE estimator is highly nonlinear with respect to the denominator estimator as the denominator approaches zero. To establish valid hypothesis tests and confidence sets for LATE without considering identification strength, we consider the function $\psi$ defined by
where $W=(Y,D,Z,X')'$, $\theta\in\Theta$ is our target parameter LATE, with $\Theta$ being a compact set on $\mathbb{R}$, and $\eta=(g,m,p)\in\mathcal{T}$\footnote{$\mathcal{T}$ is assumed to be a convex set because we want to ensure that $\psi(W;\theta_0,\eta_0+r(\eta-\eta_0))$ is well defined. Given that $\mathcal{T}$ is convex, $\eta_0+r(\eta-\eta_0)=(1-r)\eta_0+r\eta\in \mathcal{T}$ for all $r\in[0,1)$ and $\eta\in\mathcal{T}$.} are the nuisance parameters. In Section (ref), we demonstrate that the function $\psi$ adheres to the Neyman orthogonality condition. It is important to note that the score $\psi$ satisfies the moment condition ${\mathrm{E}_{P}}[\psi(W;\theta_0,\eta_0)]=0$, where $\theta_0$ and $\eta_0$ are the true values of $\theta$ and $\eta$, respectively. The orthogonal score is the same as that in Section 5.2 in CCDDHNR18, specialized to the LATE setting. With these properties, we can describe the function $\psi$ as an Anderson-Rubin-type (AR-type) Neyman orthogonal score function for the model (ref)-(ref).
We next introduce how to make inferences about the target parameter $\theta\in\Theta$. For each candidate value $\theta_0$, we test the moment condition $H_0(\theta_0): {\mathrm{E}_{P}}[\psi(W;\theta_0,\eta_0)]=0$ against $H_1(\theta_0): {\mathrm{E}_{P}}[\psi(W;\theta_0,\eta_0)]\neq 0$. Because the hypothesis is expressed directly in terms of the population moment, it imposes no restriction on the strength of identification. The test is valid whether the compliance probability is large, or small.
Initially, we estimate the first-stage nuisance parameters $\eta$, using some machine learning methods. With a fixed positive integer $K>1$, we randomly partition $\{1,\cdots,N\}$ into $K$ parts, denoted as $\{I_k\}_{k=1}^K$. For each $k\in\{1,\cdots,K\}$, the nuisance parameter estimate $\widehat{\eta}_k$ is computed using the subsample of those observations with index $i\in I_k^c$. Subsequently, we employ the cross-fitting/data-splitting method, as suggested by CCDDHNR18, to compute the covariance estimator of the process $\sqrt{N}\psi(W_i;\cdot,\eta_0)$, which is expressed as
for $\theta_1,\theta_2\in\Theta$. Observe that $\widehat{\Omega}(\theta_1,\theta_2)$ is computed using the sample of observations with index $i\in I_k$ and this computation is repeated $K$ times. For a candidate value $\theta_0$, define the cross-fitted sample moment $\widehat{q}_N(\theta_0)=\frac{1}{N}\sum_{k=1}^K\sum_{i\in I_k}\psi(W_i;\theta_0,\widehat\eta_k)$. Since $\widehat{\Omega}(\theta_0,\theta_0)$ is a scalar in the just-identified LATE setting, we form the statistic
and refer to it as the orthogonalized AR statistic. Under the null hypothesis, $AR(\theta_0)$ is asymptotically $\chi^2_1$ uniformly over the strength of identification. We reject $H_0(\theta_0)$ at level $\alpha$ whenever $AR(\theta_0)>\chi^2_{1,1-\alpha}$, where $\chi^2_{1,1-\alpha}$ denotes the $(1-\alpha)$-quantile of the $\chi^2_1$ distribution. Inverting the point-wise tests over $\theta\in\Theta$ yields a $(1-\alpha)$ confidence set that remains valid whether the compliance probability is large or small.
We specifically examine a logit model class where a binary outcome $D_i$, denoting individual $i$'s receipt of treatment, is determined by the treatment offer, $Z_i$, and a set of $p$-dimensional covariates, $X_i$. Moreover, we employ the logit model to estimate the propensity score and conduct linear regression analysis to estimate the outcome regression. The models can be expressed as:
where $\Lambda$ denotes the logistic CDF defined by $\Lambda(t)=\exp(t)/(1+\exp(t))$ for all $t\in\mathbb{R}$, and the true nuisance parameter vector $\eta_0=(\beta_{11}^0,\beta_{12}^0,\beta_{21}^0,\beta_{22}^0,\gamma^0)$. The log-likelihood functions for the logit model are $L_1(\beta_{11},\beta_{12})=\mathbb{E}_N[L_1(W_i;\beta_{11},\beta_{12})]$ and $L_2(\gamma)=\mathbb{E}_N[L_2(W_i;\gamma)]$, where $L_1(W_i;\beta_{11},\beta_{12})=\log(1+\exp(Z_i\beta_{11}+X_i'\beta_{12}))-D_i(Z_i\beta_{11}+X_i'\beta_{12})$ and $L_2(W_i;\gamma)=\log(1+\exp(X_i'\gamma))-Z_iX_i'\gamma$. We estimate the conditional mean functions $m_0(Z,X)$, $g_0(Z,X)$, and $p_0(X)$ introduced in (ref)--(ref) using OLS or a logistic link function, and denote the corresponding estimators by $\widehat{m}(Z,X), \widehat{g}(Z,X)$, and $\widehat{p}(X)$\footnote{Throughout, we use $\widehat{m}(Z,X), \widehat{g}(Z,X)$, and $\widehat{p}(X)$ to denote the estimated nuisance functions obtained via cross-fitting. Specifically, for each fold $k$, $\widehat{m}_k(Z, X)$, $\widehat{g}_k(Z, X)$, and $\widehat{p}_k(X)$ are estimated using the data excluding fold $k$, and applied to observations in fold $k$. For notational simplicity, we write $\widehat{m}$, $\widehat{g}$, and $\widehat{p}$ without the fold subscript except when additional clarity is required.}. The AR-type Neyman orthogonal score is then specified as
Note that within the score function, the logit model can be easily replaced by other models, such as the probit model or linear probability model. To illustrate the inference procedure, we outline a specific inference procedure in the subsequent algorithm. Although our algorithm primarily employs the lasso for illustration, other machine learning methods may be used in its place. We set penalty level $\lambda^k_1=\lambda^k_2=\lambda^k_3=1.1\sqrt{|I_k^c|}\Phi^{-1}(1-0.025/p)$ and construct penalty loading $\widehat\Psi_1^k, \widehat\Psi_2^k,$ and $\widehat\Psi_3^k$ based on Algorithm 6.1 from belloni2017program. Formally and theoretically justified choices of penalty loadings and penalty levels are elaborated in Algorithm (ref) and Lemma (ref) in Appendix (ref).
Section (ref) describes the algorithmic implementation of our orthogonalized AR test. We now establish its theoretical guarantees. First, we restate the score function and define the empirical process. Next, we present the regularity conditions required for the analysis. Finally, we prove two main results: a uniform functional central limit theorem and a uniform size result for the test statistic.
The test is built on the score $\psi(W;\theta,\eta)$ in equation (ref). At the true parameter value $(\theta_0,\eta_0)$, the score has zero mean and is Neyman orthogonal, which means the G\^ateaux derivative with respect to $\eta$ vanishes. Under these properties, estimation error in the high-dimensional nuisance functions contributes only an $o_p(N^{-1/2})$ remainder to the statistic. Appendix (ref) provides complete proofs.
For $\theta\in\Theta$, define $q_N(\theta)=N^{-1}\sum_{i=1}^N\psi(W_i;\theta,\eta_0)$ and let $S_N(\cdot)$ denote its expected value, given as $S_N(\cdot)={\mathrm{E}_{P}}[q_N(\cdot)]$. With this notation, we now define an empirical process $\mathbb{G}_N(\cdot)$ as
Later we show that under mild conditions, the process $\mathbb{G}_N(\cdot)$ weakly converges to a mean-zero Gaussian process $\mathbb{G}(\cdot)$ with a covariance function $\Omega(\theta_1,\theta_2)={\mathrm{E}_{P}}[\mathbb{G}(\theta_1)\mathbb{G}(\theta_2)]$. In practice $\eta_0$ is replaced by the cross-fitted estimator $\widehat{\eta}_k$ from Algorithm (ref). We propose an estimator for $\mathbb{G}_N(\theta)$ as
where the second term is a population expectation used solely for theoretical centring and is not evaluated in practice. An estimator of $q_N(\theta)$ is proposed as $\widehat{q}_N(\theta)=N^{-1}\sum_{k=1}^K\sum_{i\in I_k}\psi(W_i;\theta,\widehat\eta_k)$. An estimator of the covariance function is
Next we present the theoretical foundation for the orthogonalized AR test. The argument follows the empirical process framework of CCDDHNR18 but extends it to cover weak identification in the LATE setting. Whereas CCDDHNR18 establish the asymptotic normality of a $\sqrt{N}$-consistent estimator of $\theta$, we show the weak convergence result of our proposed empirical process under the null. We demonstrate that our test controls size uniformly over both weak- and strong-identification regimes.
To fix notation, for any finite-dimensional vector $\delta$, we write $\|\delta\|_1$ for the $\ell_1$ norm, $\|\delta\|_\infty$ for the sup norm, and $\|\delta\|_0$ for the number of nonzero components of $\delta$. For a function $f$, define the $L_q(P)$ norm by $\|f\|_{P,q} = ({\mathrm{E}_{P}}[f(W)^q])^{1/q}$. We define the sample expectation operator as $\mathbb{E}_N[\cdot]=\frac{1}{N}\sum_{i=1}^N[\cdot]$. The prediction norm of $\delta$ is given by $\|x_{ij}'\delta\|_{2,N}=\sqrt{\mathbb{E}_N[(x_{ij}'\delta)^2]}$. Let us define $\mathcal{T}_{N(i)}$ as the parameter space of the $i$-th parameter in $\eta=(\eta_1,\eta_2,\eta_3)$ with $i\in\{1,2,3\}$. The sequence $\{s_N\}_{N\geq 1}$ is a set of positive integers greater than 1. Let $c_0, c_1, C_1$ be some finite and positive constants. Let $a_N=p\vee N$. Let the sequence $\{M_N\}_{N\geq 1}$ be a set of positive constants such that $M_N\geq ({\mathrm{E}_{P}}[(Z_i\vee\|X_i\|_\infty)^{2\omega}])^{1/2\omega}$, where $\omega$ is a positive constant with $\omega>4$. Let $\{\Delta_N\}_{N\geq 1}$ be a sequence of positive constants that converges to zero at a speed at most polynomial in $N$. For any $T\subset[p+1]$, $\delta=(\delta_1,\cdots,\delta_{p+1})'\in\mathbb{R}^{p+1}$ with $\delta_{T,j}=\delta_j$ if $j\in T$ and $\delta_{T,j}=0$ if $j\notin T$. Define the minimum and maximum sparse eigenvalue by
Throughout, $B_N={\mathrm{E}_{P}}[m_0(1,X)-m_0(0,X)]$ is the compliance probability and $\kappa_N^2 = NB_N^2/{\mathrm{E}_{P}}[\text{Var}(D|Z,X)]$ is the concentration parameter analogue introduced in Section (ref). Recall that $\mathcal{P}_0$ denotes the collection of probability laws under the null.
Assumption (ref) establishes the identification assumption for the LATE framework. Assumption (ref)(i) requires that the instrument $Z$ is independent of potential outcome and potential treatment conditional on $X$, ensuring quasi-random assignment. Assumption (ref)(ii) stipulates that the instrument $Z$ affects the outcome solely through its impact on the treatment, with no direct effect on the outcome itself. Assumption (ref)(iii) rules out the presence of defiers by assuming the instrument cannot decrease the probability of treatment for any unit. Assumption (ref)(iv) is a standard overlap condition, indicating that for every value of the covariates $X$, there is a non-zero probability that a unit will either be treated or remain untreated.
Next assumption is to guarantee good performance of the target estimator under the linear and logistic link functions.
Assumption (ref)(i) specifies that the number of nonzero components in the high-dimensional nuisance parameter vector is controlled by the sparsity index $s_N$, which is further bounded in part (iv)(a). Assumption (ref)(ii) imposes the sparse eigenvalue condition, analogous to the RE condition in bickel2009simultaneous. Assumption (ref)(iii) imposes constraints on the error term $u$ from the reduced form equation (ref). Specifically, it establishes lower and upper bounds on the conditional second moment of $u$, as well as higher moment conditions on the covariates and outcomes. Assumption (ref)(iv) imposes regularity conditions on both the approximation errors and the empirical errors, in a manner analogous to Assumptions 6.1 and 6.2 of belloni2017program.
These conditions are sufficient for the high-level conditions invoked in Appendix (ref).
CCDDHNR18 establish $\sqrt{N}$-consistency and Wald-type inference for their debiased estimator under strong identification, which requires the smallest singular value of the Jacobian to be uniformly bounded away from zero. Assumption (ref)(v) relaxes this condition by allowing the concentration parameter to remain bounded, equivalently $\sqrt{N}B_N=O(1)$. By combining the Neyman orthogonal score with the AR statistic, we obtain weak convergence of the test under the null uniformly in instrument strength, achieving size-correct inference whether the first stage is strong or arbitrarily weak.
This shift from identification-dependent estimation to identification-robust inference yields a broader insight: by combining orthogonalization with AR test inversion, inference validity can be decoupled from first-stage strength, ensuring uniformly valid hypothesis testing even when consistent estimation is infeasible. In high-dimensional settings, conventional diagnostics for instrument relevance can be unreliable, as the inclusion or regularized selection of many covariates effectively partials out the identifying variation. As shown by belloni2014inference, such over-partialling can weaken the effective first stage and distort standard IV inference, highlighting the importance of orthogonalization-based procedures that safeguard inference when identification strength is uncertain.
This section describes the DGP, the simulation scenarios, and the finite-sample performance of the proposed orthogonalized AR test, in comparison with conventional alternatives.
For each replication we draw an i.i.d. sample $\{(Y_i,D_i,Z_i,X_i')\}_{i=1}^{N}$ with sample sizes $N\in\{50,100\}$ and covariate dimensions $p\in\{5,10,25,35,50,100\}$. The covariates are generated as $X_i\sim \mathcal{N}(0,\Sigma)$, where $\Sigma_{jk}=0.5^{|j-k|}$. The binary instrument is defined by
with $\gamma_0=-0.08$, $\gamma=(0.5,0.5^{2},\dots,0.5^p)'$, and $h(X_i)=X_{i1}^2-1$. The parameter $\rho_Z$ controls the degree of nonlinear misspecification in the instrument, and the intercept $\gamma_0$ is set so that the instrument is balanced, with equal probabilities of taking values zero and one. Instrument strength varies smoothly through a latent compliance model, with the complier share defined as $P_C=\kappa/\sqrt{N}$ with $\kappa\in\{1.5,3,4.5,6\}$. The shares of always-takers and never-takers are set to $P_{AT}=P_{NT}=(1-P_C)/2$. The implied concentration parameters range from $\kappa_N^2\in\{1.88,8.87,27.23,102.86\}$. We draw $\eta_i\sim\mathrm{Unif}(0,1)$. Potential treatments are given by $D_i(0)=\mathbbm 1\{\eta_i<P_{AT}\}$ and $D_i(1)=\mathbbm 1\{\eta_i<1-P_{NT}\}$, so the realized treatment is $D_i=(1-Z_i)D_i(0)+Z_iD_i(1)$. The outcome is generated as
with $\theta_0=1$, $\beta=(0.5,0.5^{2},\dots,0.5^p)'$, and $g(X_i)=0.5X_{i1}^2+\sin(X_{i2})$. The disturbance is heteroskedastic with $\sigma_i^2=(1+\rho_\sigma|X_{i1}|)^2$. Thus, $\rho_Y$ controls outcome nonlinearity, $\rho_\sigma$ governs heteroskedasticity, and $\rho_Z$ introduces nonlinear dependence of the instrument on covariates. Setting $(\rho_Y,\rho_\sigma,\rho_Z)=(0,0,0)$ recovers the baseline linear homoskedastic design. Positive values of $(\rho_Y,\rho_\sigma,\rho_Z)$ generate misspecification and heterogeneity to assess robustness.
Monte Carlo simulations are conducted with 3,000 iterations for each set. We compare three approaches in the simulation study: the conventional AR test, the DML method by CCDDHNR18, and the proposed orthogonalized AR (OAR) test. OAR and DML procedures use $K=5$ folds for cross-fitting. For the AR test, the nuisance functions, namely the propensity score, the first-stage, and the reduced-form, are estimated using logistic regression or OLS with all covariates included without model selection or regularization. This classical low-dimensional implementation is feasible only when the design matrix has full column rank and positive residual degrees of freedom. To satisfy these conditions, in the nominal $p\in\{50,100\}$ designs we set $p=48$ for $N=50$ and $p=98$ for $N=100$, ensuring that $X'X$ is nonsingular.
In high-dimensional designs, conventional regression either fails or produces severe overfitting, so the AR test cannot be applied. In particular, no AR results are reported for $N=50$ with $p=100$. By contrast, the DML method employs Lasso-regularized regression for nuisance estimation, which enables valid inference even when $p \gg N$. Confidence intervals are constructed using a Wald-type approach based on the asymptotic normality of the debiased estimator. Our proposed OAR method combines the identification-robust test inversion of the AR procedure with regularized nuisance estimation, extending AR inference to high-dimensional settings.
Table (ref) reports empirical size across identification strength, sample size, and dimensionality under the baseline linear and homoskedastic design $(\rho_Y,\rho_{\sigma},\rho_Z)=(0,0,0)$. The AR test shows pronounced size inflation that worsens as $p$ increases. Even at $N=50$ and $p=5$ the empirical size exceeds 5% and it rises sharply with $p$, which indicates high sensitivity to dimensionality. By contrast, the DML method shows size deflation under weak or moderately weak instruments at $\kappa_N^2 \in \{1.88, 8.87\}$ with rejection rates well below the 5% nominal level. It approaches nominal size only when instrument strength increases to $\kappa_N^2 \in \{27.23, 102.86\}$. The proposed OAR test stays near the nominal level across identification strengths and values of $p$, with only minor drift in the highest-dimensional designs. Overall, Table (ref) shows that OAR is stable in both identification strength and dimensionality, while AR over-rejects increasingly with dimensionality and the DML method under-rejects when instruments are weak.
Figures (ref) and (ref) plot power curves for OAR, AR, and DML under the baseline linear and homoskedastic design $(\rho_Y,\rho_{\sigma},\rho_Z)=(0,0,0)$. The AR test shows size inflation that worsens as $p$ increases. At $\theta=1$ the size lies well above the 5% nominal level across identification strengths. The DML method attains high power under strong identification, as seen in Figure (ref), but displays size deflation under weak or intermediate identification, as seen in Figure (ref), with size far below 5%. In contrast, OAR tracks the nominal size closely across $\kappa_N^2$ and $p$. Under strong identification its power is essentially on par with DML, while under weak identification it maintains size and improves on AR.
Table (ref) compares performance across instrument strengths $\kappa_N^2 \in \{1.88, 8.87, 27.23, 102.86\}$ for $N=50$, varying the dimensionality $p$ and considering three designs: (i) the baseline linear and homoskedastic specification $(\rho_Y,\rho_{\sigma},\rho_Z)=(0,0,0)$, (ii) outcome misspecification with heteroskedastic errors $(1,0.5,0)$, and (iii) instrument nonlinearity $(0,0,0.5)$. The conventional AR test exhibits pronounced size inflation that becomes more severe as $p$ increases, with the largest distortions under design (iii). The DML method performs well only under strong identification with $\kappa_N^2 \in \{27.23, 102.86\}$, whereas under weak or intermediate identification it displays systematic size deflation across all three designs, which indicates sensitivity to weak instruments despite tolerance to misspecification. In contrast, the proposed OAR test maintains near-nominal size across identification strengths, dimensionalities, and designs. Under strong identification, differences between OAR and DML are negligible. Overall, the simulations indicate that OAR is robust to misspecification, dimensionality, and instrument strength, that the AR test is highly sensitive to dimensionality and misspecification, and that the DML method, while resilient to misspecification, fails to deliver correct size under weak identification.
To demonstrate the methods outlined in the preceding sections, we revisit the instrumental variable analysis by hornung2015railroads concerning the impact of railroad access on city growth in 19th-century Prussia. In this study, straight-line corridors between major cities (nodes) are constructed, and whether a city is located on this line is used as an instrument for analysis. We compare the proposed orthogonalized AR test with the conventional AR test and the DML method to assess the effect of railroad access on city population growth. Our goal is to deepen our understanding of the conclusions presented in the literature. By conducting a new empirical analysis, we keep two econometric considerations in mind: 1. the inclusion of high-dimensional covariates to mitigate unobserved confoundedness, and 2. accounting for the weak identification issue in the data. For the first time, we report confidence intervals that are robust to weak identification and high dimensionality.
Consider the empirical model:
where $Y_{it}$ denotes the urban population growth rate in city $i$ at time period $t$, $D_{i}$ is a dummy variable indicating whether there is railroad access by 1848 in city $i$, and $Z_{i}$ denotes whether city $i$ was located within a straight-line corridor between junction stations (nodes) in 1848. The covariates $X_{i}$ include a lagged dependent variable, distance to the closest node of railroad lines, age composition, primary education of the urban population, county-level concentration of large landholdings, access to main roads, rivers, and ports, pre-railroad city growth from 1831-1837, and the size of the civilian and military population in 1849.
Within this study, the exclusion restriction condition would be violated if the location of the cities in the straight-line corridor were associated with urban population growth through a channel other than the railroad. The author asserts that the exclusion restriction is satisfied in Hornung (hornung2015railroads, pg. 714),
In the context of this study, it is noteworthy that the adoption of railroad technology by cities located on a straight line between two important cities was randomly assigned. This random assignment arises because the positioning of these cities along such lines was not intentionally controlled by any specific entity. In 19th-century Prussia, the government did not dictate railroad construction due to financial limitations. Instead, the decision fell to individual city councils negotiating with private railroad enterprises. Hence, each city had the autonomy to determine whether or not to proceed with railroad construction. Within this study, “compliers” refer to (1) cities situated on the straight line between two major cities AND eventually established a railroad station, and (2) cities NOT on such a line AND did NOT get a train station. The second part does not exist in this study as mentioned in Hornung (hornung2015railroads, pg. 731),
We implement the proposed method on the city-level railroad data from hornung2015railroads. As outlined in Table 5 of hornung2015railroads, the first-stage F-statistics vary between 26.46 and 38.29. This variation suggested instrument weakness based on the $tF$ critical value function proposed by lee2022valid. We conduct a re-analysis by incorporating the polynomial and interaction terms of the original covariates, and present the results in Table (ref).
Table (ref) summarizes the results. To emphasize the robustness of the proposed method within the high-dimensional framework, we report the confidence intervals for the LATE estimates and the lengths of the confidence intervals. Panel A reports results from the conventional AR test without a regularization step. Panels B–E report results from DML with Lasso and from the proposed OAR with Lasso, Ridge, and Elastic Net, computed using 5-fold cross-fitting. Panel F reports the low-dimensional results from Table 5 of hornung2015railroads. Panels A--E use an expanded, high-dimensional set of covariates. Different columns report the results for several dependent variables across diverse time periods. Specifically, Columns (1) and (2) capture outcomes spanning two main periods: 1831-37, and 1849-1871. Columns (3)-(9) depict findings across seven subperiods. To mitigate the uncertainty induced by sample splitting, we compute confidence intervals based on the median from 100 repetitions of resampled cross-fitting following CCDDHNR18.
Table (ref) presents a comparison of the three inference methods and illustrates the distinct properties of each approach. The DML estimator relies on a delta-method linearization of the LATE ratio, whose validity requires that the denominator be bounded away from zero. When this condition fails under weak identification, the ratio becomes highly nonlinear, and the Wald approximation produces confidence intervals that are systematically too narrow. Hence, the short DML intervals in Table (ref) do not imply higher precision but rather reflect the lack of robustness of Wald-type inference when the first stage is weak. In contrast, our proposed orthogonalized AR test directly tests the moment condition without relying on a linear approximation. It therefore achieves uniform size control irrespective of identification strength. The conventional AR test also has identification-robust properties in low dimensions, but when applied with a large number of controls and no regularization, it suffers from overfitting and instability in covariance estimation, leading to distorted inference. Our orthogonalized AR approach combines the identification-robustness of the AR framework with regularized nuisance estimation and cross-fitting, ensuring valid inference even when the covariate dimension is large relative to the sample size.
Empirically, the results align with these expectations. With the full high-dimensional control set, the OAR confidence intervals are wider than those obtained in DML, consistent with proper uncertainty quantification when the first stage may be weak, yet still informative. OAR continues to yield statistically significant positive effects in several economically meaningful subperiods (e.g., 1852–55, 1861–64, 1867–71, and for the aggregate 1849–71 period), where the lower bound of the 95% interval lies above zero. By contrast, a number of subperiod effects that appear significant under DML or under the unregularized AR approach become insignificant under OAR. This indicates that those earlier findings were sensitive to the linear approximation underlying Wald inference or to overfitting in the high-dimensional setting. Economically, these results suggest that the strong and precisely estimated effects in the original 2SLS specification in hornung2015railroads likely reflected a limited control set that overstated instrument strength and precision. When model uncertainty and weak identification are incorporated, the statistical evidence for a consistently positive effect across subperiods becomes weaker. Importantly, OAR produces the same qualitative conclusions across Lasso, Ridge, and Elastic Net regularization, illustrating that our identification-robust results do not hinge on any particular form of penalization.
Table (ref) reports OAR results using a reduced set of controls, with the number of covariates $p$ ranging from 66 to 90. The 95% confidence intervals are very similar to those in Table (ref). Trimming the control set does not alter the substantive conclusions, which suggests that our findings are not driven by the inclusion of specific covariates, provided the specification remains high dimensional.
In this subsection, we reexamine the instrumental variable estimation by ambrus2020loss concerning the long-term consequences of the 1854 cholera outbreak on housing prices. The authors investigate the impacts of a cholera epidemic in a neighbourhood of 19th-century London on property values in 1864, a decade following the outbreak. Page 479 of ambrus2020loss presents the background information,
In this context, $Y_i$ represents the log rental price of house $i$ in 1864. The variable $D_i$ is an indicator, set to 1 if house $i$ has at least one cholera death. ambrus2020loss use two instruments $Z_i$; we focus on their preferred instrument, which is an indicator for whether the property $i$ falls inside the Broad Street pump (BSP) catchment areas, which were the primary contaminated area during the outbreak. The controls $X_i$ comprise all house characteristic variables listed in Table 1 of ambrus2020loss, such as distance to the closest pump, distance to the fire station, distance to the urinal, sewer access, among a total of 26 variables.
In this study, the “compliers” refer to (1) houses situated within the cholera-affected contaminated areas AND witnessed at least one cholera-related death, and (2) houses outside these contaminated zone AND did not experience any cholera fatalities. However, the latter category does not actually exist in the study since the authors limit the sample to properties within a certain distance of the BSP boundary.
Table B2 in ambrus2020loss reports an IV estimate of \(-0.694\) with a 95% interval \([-1.333,-0.055]\), which excludes zero. The reported first-stage \(F\) statistics are close to 10, which signals possible weak identification under the \(tF\) critical-value function of lee2022valid. Table (ref) revisits this application with two control sets, \(p=52\) in Panel A and \(p=26\) in Panel B. In both panels the conventional AR intervals are negative and exclude zero. By contrast, the DML intervals straddle zero, and the OAR intervals also straddle zero and are wider than DML, which is the pattern expected from identification-robust inversion under potentially weak first-stage signal. The OAR conclusions are stable across regularization choices (Lasso, Ridge, and Elastic Net), and they are unchanged when moving from \(p=26\) to \(p=52\). Taken together, the reanalysis indicates that once inference is made identification-robust, the effect of cholera-related deaths on 1864 rents is not statistically distinguishable from zero in our specification, and this conclusion is not driven by the choice of controls or penalty.
Given these insights, we advocate for our proposed orthogonalized AR method in high-dimensional models, especially when weak identification may be a concern.
In this paper, we address the challenge of weak identification within the LATE framework, especially in the presence of high-dimensional covariates. Our primary contribution lies in the introduction of an identification-robust inference method for the high-dimensional LATE framework. This is paired with a user-friendly algorithm for inference and confidence interval construction for the LATE estimate. We validate the uniformly correct asymptotic size of our proposed method.
There are several potential directions for future research. First, while we rely on the typical microeconometric assumption of i.i.d. sampling, future work could address theoretical developments for data structures exhibiting complex dependence, such as separately and jointly exchangeable arrays. Second, whereas we focused on a single instrument scenario, it would be valuable to develop methods and theories for settings with many weak instruments, as discussed by mikusheva2022inference, mikusheva2024weak and matsushita2024jackknife. We leave these and other extensions for future research.