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.
83,567 characters · 15 sections · 89 citation commands
Jackknife Instrumental Variable Inference
\doublespacing
As emphasized by MS22, many contemporary empirical applications of instrumental variables (IV) estimation involve a large number of instruments, often with limited identifying content per instrument. Beyond the settings discussed by these authors, similar many-instrument environments arise in shift-share (Bartik) designs that exploit finely disaggregated shocks BorusyakHullJaravel22,GoldsmithPinkhamSorkinSwift20, in Mendelian randomization studies that use large sets of genetic variants as instruments Davies15,Sanderson22,PLB24, and in recent machine learning–assisted IV approaches that construct instruments from high-dimensional predictors BelloniChernozhukovHansen12,ChernozhukovChetverikovDemirer18. In such settings, instrument proliferation combined with weak first-stage relationships can undermine conventional IV inference, creating a demand for robust procedures that address not only estimation but also hypothesis testing. We address this need by developing a unified trinity of Wald, Lagrange multiplier, and distance-type tests for linear IV models with many potentially weak instruments.
The construction and asymptotic analysis of estimators and test statistics often rely on the definition of a given objective function. In this context, a typical analytical strategy consists of using some stochastic expansion argument and the application of a suitable central limit theorem (CLT) NM94. Classical likelihood theory Engle84 and inference based on the generalised method of moments NW87, NW09 exploit this strategy to derive the well-known {\it trinity} of tests, Wald ($W$), Lagrange multiplier ($LM$) and likelihood ratio ($LR$) in maximum likelihood, and the corresponding distance statistic ($D$) in GMM. In line with this idea, we develop an analogous trinity of tests for linear IV models with many weak instruments.
In particular, this paper proposes a unified approach to inference for linear regression models with multiple endogenous and exogenous variables, heteroskedastic disturbances and many (potentially weak) instrumental variables (IVs). A jackknife device is typically applied to remove the noise component due to the presence of heteroskedasticity while, ideally, retaining the whole signal (matrix) of the instruments HNWCS12,CHNSW14,BC15,CMS21,MS22,Yap23,MO22.\footnote{See PH77 and AIK99 for earlier applications of the jackknife to IV models and BD99 and DM06 for extensive simulation studies on the properties of jackknife IV estimators in finite samples.} In our setup, we allow the number of instruments, say, $k$ to grow proportionally with the sample size, say, $n$ and $k\le n$. This is a standard assumption in this literature. There are, however, some studies that allow the number of instruments to be larger than the sample size HK14,DKM24. Furthermore, while the number of instruments is allowed to diverge, we impose a restriction on the minimal amount of signal the instruments should collectively carry. As shown in MS22, this is a necessary condition to obtain consistent IV estimators and tests.\footnote{The sufficiency part was proven in CS05.} Our test statistics are derived from two types of objective functions. The first type is defined by quadratic forms, such as those underlying the jackknife IV estimators JIVE1 and JIVE2 AIK99, CMS21,MS22. The second type is defined by ratios of quadratic forms, such as the objective function for the limited information maximum likelihood (LIML) estimator AR49, Bekker94, HHN08 and the more general jackknife-based alternatives heteroskedasticity-robust LIML (HLIM) and symmetric JIVE (SJIVE) of HNWCS12 and BC15, respectively.\footnote{The LIML case is the object of a companion paper.} We show that the tests directly related to objective functions are asymptotically distributed as a weighted average of chi-squares, sometimes referred to as chi-bar-square distribution Vuong89, Hansen21. The weights of this distribution can be estimated and p values can be computed Farebrother84. However, by appropriately modifying the objective function, it is possible to obtain more conventional chi-square distributions in the limit.
Via extensive Monte Carlo experiments, we investigate the tests' finite-sample properties in terms of size and power. Although the simulations do not display a clear ranking of the tests, we can still identify specific scenarios where certain statistics outperform the others. For example, when instruments are weak, some $LM$ statistics suffer from significant loss of power against distant alternatives, whereas $W$ may have low power near the null. By contrast, $D$ tends to have good size and power properties across a range of designs. Furthermore, compared to the Anderson-Rubin ($AR$) tests of CMS21 and MS22, $D$ and $W$ exhibit higher power in our experiments. In the JIVE case, these statistics are relatively easy to compute (compared to SJIVE and HLIM) and, in some cases in our simulations, they outperform $AR$ both in terms of size and power. Moreover, our tests are computationally less burdensome in our implementations then their $AR$ counterparts.
We further illustrate the practical relevance of our methodology through an empirical application that studies the causal effect of alcohol consumption on body mass index using genetic variants as instrumental variables in the UK Biobank. This setting naturally involves a large number of instruments, a feature that raises concerns about instrument strength and the reliability of conventional IV inference. In this context, the proposed statistics can be implemented straightforwardly and allow us to conduct inference without relying on standard approximations that may perform poorly when instruments are numerous or weak. The empirical results therefore complement the simulation findings by showing that the proposed procedures can be effectively applied in a realistic setting with many potentially weak instruments.
The main contribution of this paper is fourfold: first, we introduce a trinity of tests for linear IV models robust to many weak instruments; second, we establish the limiting distribution of these tests under the null; third, we compare their finite-sample performance and highlighting scenarios of dominance; four, we illustrate the applicability of the proposed tests in a realistic setting with many potentially weak instruments through an empirical study that uses genetic variants to estimate the effect of alcohol consumption on body mass index. Alongside the main contributions, we provide two auxiliary results that could be of independent interest to the readership: first, we prove a general consistency result for constrained and unconstrained jackknife-based estimators; second, we introduce a cross-fit variance estimator for the $AR$ test of CMS21 along the lines of that proposed by MS22.\footnote{See Proposition (ref) in Appendix (ref) and equation (ref) in Section (ref), respectively.}
The results in this paper are largely new, particularly for tests based on ratios of quadratic forms and the jackknife approach of BC15. They also generalise earlier results, such as the $LM$ test of HHN08 and Bekker's LIML statistic Kleibergen02, which are developed under the assumption of homoskedasticity. There are, however, some interesting overlaps. In particular, the $LM$ test built on the JIVE2 objective function coincides with that of MO22. Moreover, the Wald test based on the JIVE2 objective function is related to the results in Yap23, while the Wald test based on the HLIM objective function was first introduced by HNWCS12. This unified treatment of $W$, $LM$, and $D$ tests (the latter being novel) under many weak instruments is, to our knowledge, new. Finally, the modification applied to the objective functions to obtain a standard chi-square limit is reminiscent of that used in Kleibergen02 to obtain the $K$ statistic BK03. Also in our case, the degrees of freedom of the chi-square distribution do not depend on the number of instruments.
Our analysis also relates to the weak instrument literature of the 1990s and early 2000s, which showed that standard Wald-type procedures can break down under weak identification and thereby spurred the development of robust alternatives NS90,BJB95,Dufour97,SS97,WZ98,Kleibergen02. The present paper complements this literature, and recent contributions such as KN23, by showing how analogous concerns can be addressed in linear IV models with many potentially weak instruments through a unified testing framework.
As noted, several studies have considered $W$ and $LM$ statistics, but, to the best of our knowledge, the $D$ statistic in this paper, both in the chi-bar-square version and in the chi-square version, is new. This fills a gap in the weak-instrument literature, which so far has mostly developed $AR$-, $W$- and $LM$-type tests but lacked an analogue of the $LR$ test. The $D$ statistic we propose provides such an analogue, with promising size and power properties when compared to its competitors.
There are clear similarities between the conditional $LR$ ($CLR$) statistic introduced by Moreira03 and its more recent versions designed to allow for instrument numerosity AMO22,LWZ22, in that both $D$ and $CLR$ are based on a discrepancy between objective functions evaluated at different values. However, there is also a difference, namely that $D$ does not rely on a conditional argument.
Moreover, to the best of our knowledge, no previous work in this setting explicitly considers general linear restrictions (including the important case of testing a subset of parameters). This may be an important aspect, as in a conventional linear regression model one may partial out variables, via the Frisch-Waugh-Lovell (FWL) theorem, and conduct inference on the remaining regressors. However, as pointed out by CSW23 and MS24, the FWL theorem prevents the jackknife device from removing the noise component MS25. Introducing an inferential method for linear restrictions avoids possible issues connected to the interplay between the jackknife approach and the FWL theorem.
The remainder of the paper is organised as follows. Section (ref) introduces the model and the test statistics for simple hypotheses and linear restrictions and shows the asymptotic behaviour of the tests. Section (ref) describes how to modify our statistics to have a chi square limiting distribution. Section (ref) studies the finite sample properties of the tests in a Monte Carlo experiment and in an empirical application using genetic instruments. Finally, Section (ref) offers some conclusions. Proofs and auxiliary results are deferred to the Appendix.
Throughout the paper we use the following notational conventions: unless differently stated $a$, $\bm a$ and $\bm A$ denote a scalar, a vector and a matrix, respectively; $a_i$ is the $i$-th element of vector $\bm a$ and $A_{ij}$ is the $(i,j)$-th element of matrix $\bm A$. $\bm a'$ and $\bm A'$ are the transposes of $\bm a$ and $\bm A$. For a square matrix $\bm A$, $\lambda(\bm A)$ denotes the vector of eigenvalues of matrix $\bm A$, while $\lambda_{\min}(\bm A)$ and $\lambda_{\max}(\bm A)$ are the smallest and largest eigenvalues of $\bm A$ respectively. Moreover, $\operatorname{rk}(\bm A)$ denotes the rank of matrix $\bm A$. Finally, $\bm I$ is an identity matrix and when we need to be more specific, we include a subindex. Thus, $\bm I_n$ denotes a $n\times n$ identity matrix; similarly, $\bm \iota$ and $\bm \iota_n$ denote vectors of ones. The symbol $*$ is used for elementwise or Hadamard multiplication and $^{(n)}$ denotes elementwise power. Hence,
Analogously, $\bm A^{(-1)}$ is used for the elementwise inverse of $\bm A$, where the generic $(i,j)$ entry of $\bm A^{(-1)}$ is equal to $A_{ij}^{-1}$.
Let us consider the model
where $\bm X$ is a $n\times g$ matrix containing potentially both endogenous and exogenous variables and $\bm Z$ is a $n\times k$ (nonstochastic) matrix of instruments and $\operatorname{E}[\bm X]=\bm Z\bm \varPi$, where the components of $\bm \varPi$ are allowed to vary with the sample size $n$. Such assumptions are made for convenience and may be generalised.\footnote{ We may, for example, consider $\bm Z$ to be stochastic and in this case $\operatorname{E}[\bm X]$ should be interpreted as a conditional expectation with respect to $\bm Z$. The linearity of $\operatorname{E}[\bm X]$ may also be relaxed as suggested in, e.g., Bekker94 and CHNSW14.} The rows of the disturbance couple $(\bm \varepsilon,\bm V)$, say $(\varepsilon_i,\bm V_i')$ $i=1,\dots,n$, are independent with zero mean and covariance matrices
In this section we are interested in test statistics for the null hypothesis
based on a certain objective function $Q_n(\bm \beta)$ that is assumed to produce consistent estimators in the many instruments sense and to be robust to the presence of heteroskedasticity. We consider two types of objective functions. The first type is based on the ratio of two quadratic forms as in HNWCS12 and BC15; specifically,
where, for the SJIVE estimator of BC15, we have
where $\bm D=\operatorname{diag}(\bm P)$ is a diagonal matrix that contains the diagonal elements of $\bm P=\bm Z(\bm Z'\bm Z)^{-1}\bm Z'$. We note that in this case $\operatorname{tr}(\bm B)=\operatorname{tr}(\bm P)=k$. For the HLIM estimator of HNWCS12 we specify
Such estimators mimic the ratio of quadratic forms structure typical of LIML. Another possibility would be to define an objective function that only considers the numerator of equation (ref):
The resulting estimator is the JIVE1 if $\bm C$ is chosen as in equation (ref) or the JIVE2 if $\bm C$ corresponds to the specification in equation (ref).\footnote{Both specifications of $\bm C$ remove the source of bias that would make an estimator inconsistent in the many instruments sense when heteroskedasticity is present. However, they treat the information contained in the IVs differently. In fact, for $\bm C$ defined as in equation (ref) $\operatorname{E}[\bm X'\bm C\bm X]=\bm \varPi'\bm Z'\bm Z\bm \varPi$, while the specification in equation (ref) produces $\operatorname{E}[\bm X'\bm C\bm X]=\bm \varPi'\bm Z'\bm Z\bm \varPi-\bm \varPi'\bm Z'\bm D\bm Z\bm \varPi$. See BC15 for a discussion.} \footnote{We label the JIVE case with $\bm C$ as in equation (ref) as JIVE1. However, it is important to notice that this is a symmetric version of the JIVE1 proposed in AIK99 and it mimics the structure of the SJIVE objective function in BC15. To the best of our knowledge, the corresponding estimator was first introduced in CMS21.} Whether we use the objective functions in equation (ref) or equation (ref), the estimator
is consistent for $\bm \beta$ under $H_0$.
To test the null hypothesis in equation (ref) we consider the following statistics
where $\bm H=\bm \varPi'\bm Z'\bm C\bm Z\bm \varPi$, $r_{\min}=\lambda_{\min}(\bm H)$ is the smallest eigenvalue of $\bm H$ and $\sigma^2$ is the top left entry of matrix $\bm \varSigma=\lim\frac{1}{\operatorname{tr}(\bm B)}\sum_{i=1}^nB_{ii}\bm \varSigma_i$ as $n\to\infty$.\footnote{We note that in the SJIVE case $\bm H=\bm \varPi'\bm Z'\bm Z\bm \varPi$ while in the HLIM case $\bm H=\bm \varPi'\bm Z'\bm Z\bm \varPi-\bm \varPi'\bm Z'\bm D\bm Z\bm \varPi$. Notice also that $\bm \varSigma$ depends on the choice of $\bm B$. To avoid notation clutter, we omit this dependency.} Furthermore,
When $Q_n(\bm \beta)$ is the objective function of a JIVE estimator (equation (ref)), $\sigma^2=1$ and $\lambda(\bm \beta)=0$. The limiting distribution of the tests in equations (ref) to (ref) is derived in Theorem (ref) below under a set of assumptions.
The assumptions we use match those in BC15 but different versions can be found in the literature. CMS21 provide a discussion of some of the assumptions; below we discuss the other assumptions. In what follows it is understood that $c_u$ is a generic positive constant that may take different values in different instances.
Let $r_{\max}=\lambda_{\max}(\bm H)$ denote the largest eigenvalue of $\bm H$.
Some comments are in order. Assumption (ref) and Assumption (ref) are standard conditions commonly found in the literature HNWCS12,BC15,CMS21: the former is a technical condition on the behaviour of the diagonal elements of the projection matrix $\bm P$ and the latter is a regularity condition on the boundedness of the second moments of the disturbances. Assumption (ref) formalises the fact that we are dealing with potentially weak instruments. Assumption (ref) does not allow for the possibility of including exogenous regressors in the model due to the fact that it requires $\bm \varSigma_{22}$ to be positive definite, which excludes that some components of the $\bm V_{i}$'s are zero. However, it is straightforward to modify this assumption for the case when some components of the $\bm V_{i}$'s are zero, so we do not provide it here. Assumption (ref) corresponds to Assumption 5 from HNWCS12. Notice that {
} Since $\sum_{i=1}^{n}\bm C_{i}\bm C_{i}^{\prime}=\bm C^{2}$ and $\bm C_{i}\bm C_{i}^{\prime}$ is positive semidefinite, $\bm C_{i}\bm C_{i}^{\prime}\leq \bm C^{2}$ for any $i$. This and $\bm Z^{\prime}\bm C^{2}\bm Z\leq c_{u}\bm Z^{\prime}\bm Z$ CMS21 imply that
so we obtain that $\frac{1}{k^{2}}\sum _{i=1}^{n}\left\Vert \bm \varPi^{\prime}\bm Z^{\prime}\bm C_{i}\right\Vert ^{4}\leq \frac{c_{u}}{k^{2}}\operatorname{tr}\left( \bm H^{2}\right) $. This goes to $0$ if $r_{\max }/k\rightarrow0$, so a primitive condition for Assumption (ref) is that the instruments are weak. Another primitive condition, which is also valid for strong instruments, is $\bm C_{i}\bm C_{i}^{\prime}\leq k^{-\alpha}\bm C^{2}$ for any $i$ and for some $\alpha>0$. This means that there is no observation $i$ for which $\bm C_{i}\bm C_{i}^{\prime}$ is excessively large. In this case $\frac{1}{k^{2}} \sum_{i=1}^{n}\left\Vert \bm \varPi^{\prime}\bm Z^{\prime}\bm C_{i}\right\Vert ^{4}\leq \frac{k^{-\alpha}}{k^{2}}\operatorname{tr}\left( \bm H^{2}\right) $, which goes to $0$ by Assumption (ref). Finally, Assumption (ref) corresponds to Assumption 6 from HNWCS12.
The next result shows that the test statistics specified above are distributed asymptotically as a weighted sum of chi squares.
The statistics $D$, $LM$, $W$ are not feasible as they are based on the unknown quantities $\bm H$ and $\sigma^{2}$. By replacing these by their feasible counterparts
respectively, we obtain the feasible statistics
where $\widehat{r}_{\min}=\lambda_{\min}(\widehat{\bm H})$ and $\widehat{\bm \psi} _{n}=\dfrac{1}{\sqrt{k}}\widehat{\bm H}(\widehat{\bm \beta}-\bm \beta_{0})$.
It is interesting to notice that test statistics based on quadratic forms such as Wald and $LM$ can be easily modified to obtain a more conventional chi square limiting distribution MO22. This seems not to be the case for $D$ statistics, unless we take into account the simple case with one endogenous variable ($g=1$). The following example illustrates that case.
Let us assume that the null hypothesis is defined as a set of linear equality restrictions
with alternative hypothesis $H_{1}:\bm A\bm \beta\neq \bm a$, where $\bm A$ is a $p\times g$ matrix with $g\geq p$ and $\operatorname{rk}(\bm A)=p$. The corresponding Lagrangian is \[ \mathcal{L}(\bm \beta,\bm \gamma)=Q_{n}(\bm \beta)+2\bm \gamma^{\prime}\left( \bm A\bm \beta -\bm a\right) , \] where $Q_{n}(\bm \beta)$ is defined as in equation (ref). The first order conditions (FOCs) are (see also footnote on page (ref))
where recall $\widehat{\bm C}(\bm \beta)=\bm C-\lambda(\bm \beta)\bm B$ and $\lambda(\bm \beta)=\frac{1} {k}Q_{n}(\bm \beta)$. Let us denote by $\widetilde{\bm \beta}$ and $\widetilde{\bm \gamma}$ the values of $\bm \beta$ and $\bm \gamma$ that solve the FOCs (i.e., first order conditions). The test statistics for the null in equation (ref) are
where
and
The next result shows that these test statistics are distributed asymptotically as a weighted sum of chi squares.
Under conditions presented in this section, it is possible to show that the proposed test statistics follow a standard chi square limit instead of a weighted sum of chi squares. In general, it is unclear which approximation provides better finite sample behaviour, yet the former is typically more commonly used by practitioners than the latter.\footnote{Section (ref) explores the finite sample properties of test statistics under either asymptotic distribution.} By using appropriate modifications, we can derive a set of test statistics that have a standard chi square limit. In some cases such modifications are relatively straightforward, as in the $LM$ case of MO22 and, in general, for tests defined by quadratic forms. When dealing with distance statistics, the adjustments are not so obvious.
Let us suppose we want to define a test statistic for the simple null hypothesis (ref). In order to find a chi square limit for the distance statistic we modify the objective function in the following way
where $\bm J=\bm C\bm X\bm \varPhi^{-1}\bm X^{\prime}\bm C$. Let $T^{\ast}$ denote a generic test statistic, either distance, Lagrange multiplier or Wald, that is obtained from equation (ref) and with $\widehat{\bm \beta}$ the estimator obtained from the optimisation of the objective function in equation (ref). Specifically, $T^{\ast}$ is one of the test statistics
The following theorem shows that under the null $T^{\ast}$ converges in distribution to a chi square.
If we are testing linear restrictions as defined in ((ref)), we modify the objective function to
where $\bm J_{a}=\bm C\bm X\bm \varGamma_{n}^{\prime}\left( \bm \varGamma_{n}\bm \varPhi\bm \varGamma_{n}^{\prime }\right) ^{+}\bm \varGamma_{n}\bm X^{\prime}\bm C$ with $\bm \varGamma_{n}=\bm A^{\prime} (\bm A\bm H^{-1}\bm A^{\prime})^{-1}\bm A\bm H^{-1}$ and $\left( \bm \varGamma_{n}\bm \varPhi\bm \varGamma_{n}^{\prime}\right) ^{+}$ is a generalised inverse of $\bm \varGamma_{n}\bm \varPhi\bm \varGamma_{n}^{\prime}$ that satisfies $\bm \varGamma_{n} \bm \varPhi\bm \varGamma_{n}^{\prime}\left( \bm \varGamma_{n}\bm \varPhi\bm \varGamma_{n}^{\prime}\right) ^{+}\bm \varGamma_{n}\bm \varPhi\bm \varGamma_{n}^{\prime}=\bm \varGamma_{n}\bm \varPhi\bm \varGamma_{n}^{\prime}$ and is bounded as $n\rightarrow\infty$.\footnote{Note that since $\bm \varGamma_{n}$ is idempotent, it is not of full rank, and therefore, $\bm \varGamma_{n}\bm \varPhi\bm \varGamma _{n}^{\prime}$ is not of full rank either.} The test statistics for the null hypothesis in equation (ref) are
The following result generalises Theorem (ref).
In this section, we first describe how to compute feasible versions of the proposed statistics by plugging in estimators for the unknown quantities. Moreover, we investigate the statistics' finite sample properties in terms of size and power in the context of two data generating processes (DGPs). DGP1 is similar to that in MO22. The model in DGP2 can be seen as a stylised version of a Cobb-Douglas production function with two endogenous variables corresponding to production factors, say labour and capital, as in NR08 CMS21. Finally, we discuss the application of our tests to a real data example that employs genetic variants as instruments to study the effect of alcohol consumption on body mass index. The results of our tests are compared against the AR tests of CMS21 and MS22 with both naive and cross-fitted variance estimators. The general expression of the AR test is
where $\bm \varepsilon=\bm y-\bm X\bm \beta$, while $\bm C$ in equation (ref) produces the AR test of CMS21 and $\bm C$ in equation (ref) produces the AR test of MS22. The naive variance estimator for $\omega(\bm \beta)$ is
while the cross-fitted alternative is
with
where $\mbox{diag}(\bm B)$ is a diagonal matrix containing the diagonal entries of $\bm B$ in the main diagonal. For the version of CMS21 $\bm B$ is as in equation (ref), while for the version of MS22, $\bm B=\bm I-\bm P$.
To make the tests operational, we first need to estimate the relevant quantities $\bm H$, $\sigma^2$, $\bm \varPhi$. Estimators for $\bm H$, $\sigma^2$, $\bm \varPhi$ are already defined in BC15 HNWCS12. Specifically, we have $\widehat\bm H(\bm \beta)=\bm X'\widehat\bm C(\bm \beta)\bm X$, $\widehat\sigma^2(\bm \beta)=\frac{1}{\operatorname{tr}(\bm B)}(\bm y-\bm X\bm \beta)'\bm B(\bm y-\bm X\bm \beta)$ and {
} with $\widetilde\bm X(\bm \beta)=\bm X-\frac{\bm \varepsilon\widehat\bm \sigma_{12}(\bm \beta)}{\widehat\sigma^2(\bm \beta)}$ and $\bm D_{\bm \varepsilon}$ a diagonal matrix with $\bm \varepsilon=\bm y-\bm X\bm \beta$ along the main diagonal. We further define the estimator for the variance covariance matrix $\bm \varSigma$ as
where $\widehat\bm \varOmega$ is an estimator of $\bm \varOmega=\operatorname{E}[(\bm y,\bm X)'\bm B(\bm y,\bm X)]$ BC15. Thus, $\widehat\bm \sigma_{12}(\bm \beta)$ and $\widehat\sigma^2(\bm \beta)$ correspond to the off-diagonal entry and the top left entry of $\widehat\bm \varSigma(\bm \beta)$, respectively. These quantities depend on the parameter vector $\bm \beta$, which needs to be replaced by a suitable plug-in estimator. Wald statistics use the unrestricted estimator $\widehat\bm \beta$, Lagrange multiplier statistics use either $\bm \beta_0$ or the restricted estimator $\widetilde\bm \beta$, depending on the type of null hypothesis being tested. Finally, distance statistics use the objective function evaluated at both the restricted and unrestricted estimators Engle84; however, the multiplicative coefficient $r_{\min}\sigma^2$ is evaluated at the unrestricted estimator, say, $\widehat r_{\min}\widehat \sigma^2(\widehat \bm \beta)$ with $\widehat r_{\min}=\lambda_{\min}\left(\widehat \bm H(\widehat\bm \beta)\right)$. The correction term in the distance statistic, when the test is chi-square distributed, is evaluated at $\widehat\bm \beta$.
When the asymptotic distribution is a weighted average of chi-squares, the calculation of the weights $\bm \varphi$ is based on $\widehat\bm H(\widehat\bm \beta)$ and $\widehat\bm \varPhi(\widehat\bm \beta)$ for distance and Wald statistics and on $\widehat\bm H(\bm \beta_0)$ and $\widehat\bm \varPhi(\bm \beta_0)$ or $\widehat\bm H(\widetilde\bm \beta)$ and $\widehat\bm \varPhi(\widetilde\bm \beta)$ for the Lagrange multiplier statistic. Similarly, for distance and Wald statistics, the minimum eigenvalue $r_{\min}$ is estimated, as seen above, by $\widehat r_{\min}=\lambda_{\min}\left(\widehat \bm H(\widehat\bm \beta)\right)$, while for the $LM$ statistic we find $\widetilde r_{\min}=\lambda_{\min}\left(\widehat \bm H(\bm \beta_0)\right)$ or $\widetilde r_{\min}=\lambda_{\min}\left(\widehat \bm H(\widetilde\bm \beta)\right)$, depending on the type of null hypothesis under study.
The plug-in estimators presented in this section are general in the sense that apply both to tests based on quadratic forms (JIVE1, JIVE2) and tests based on ratios of quadratic forms (SJIVE, HLIM). However, for the former case and in light of the discussion in Section (ref), some simplifications apply. First of all, $\widehat\bm H(\bm \beta)=\widehat\bm H=\bm X'\bm C\bm X$, $\widehat r_{\min}=\widetilde r_{\min}=\lambda_{\min}\left(\widehat\bm H\right)$ and $\widetilde\bm X(\bm \beta)=\bm X$. This implies that the variance covariance matrix $\widehat\bm \varPhi(\bm \beta)$ depends on $\bm \beta$ only through the matrix $\bm D_{\bm \varepsilon}$. Hence,
Tests that require a plug-in for $r_{\min}$ use $\widehat r_{\min}(=\widetilde r_{\min})$. Moreover, $\sigma^2=1$ and no other quantity related to $\bm \varSigma$ is required. It is also easy to see that the correction factors in the chi-square distributed distance statistics are zero. As explained in Section (ref), certain tests produce the same value. However, their corresponding p values are computed using different plug-ins for $\bm \varphi$, thus generating different rejection rates in the numerical exercises. Specifically, recall from Theorem (ref) that $\bm \varphi$ is the vector of eigenvalues of $\bm \varXi\bm \varPhi$. Thus, for the $LM$ statistic $\bm \varphi$ is computed by evaluating the estimators of $\bm \varXi$ and $\bm \varPhi$ at $\bm \beta_0$ or $\widetilde\bm \beta$. For distance and Wald statistics the estimators of $\bm \varXi$ and $\bm \varPhi$ are evaluated at $\widehat\bm \beta$. Furthermore, the chi-square tests for the general linear restrictions require the generalised inverse of $\bm \varGamma_n\bm \varPhi\bm \varGamma_n'$, for example, {
} Estimation of $\bm H$ and $\bm \varPhi$ follows the same guidelines as described above in this section. In Appendix (ref), we also consider the case in which the test statistics use a cross-fit estimator of the variance instead of the expression in equation (ref). The cross-fit estimator of $\bm \varPhi$ is adapted from MS24:
where $\bm v=\bm B\bm \varepsilon$. The definition of $\widetilde{\bm X}(\bm \beta)$ follows that of the standard case.\footnote{A proof of consistency for the estimator in equation (ref) is available upon request.}
Finally, to make our statistics operational, we need feasible counterparts of $\bm\vartheta_n$ and $\bm \xi_n$ (or $\bm \psi_n$ and $\bm \tau_n$). It is easy to see that also in this case we can use suitable plug-in estimators for $\bm H$ and $\bm \beta$.
This section investigates the finite sample properties of the proposed test statistics for four specifications of the objective function and in the context of two DGPs. DGP1 considers a model with one endogenous variable, while DGP2 focusses on a model with two endogenous variables. Since under the alternative we allow for the presence of estimable nuisance parameters, the test statistics used in the simulations are those described in Theorem (ref) and Theorem (ref).\footnote{To avoid notation clutter, in tables, figures and the discussion of the results we drop the subscript $a$.} We study both size and power properties. For size, each experiment uses 5000 replications, whereas power experiments use 1000 replications. For the sake of conciseness, the power comparisons in Figure (ref) to Figure (ref) involve tests based on SJIVE and JIVE1. As mentioned later in this section, the competing alternatives (based on HLIM and JIVE2) tend to be dominated by SJIVE and JIVE1 statistics, particularly in the case of Lagrange multiplier statistics. Since there are no noticeable differences between $W_1/W_1^*$ and $W_2/W_2^*$, Table (ref) and Table (ref) feature only $W_1/W_1^*$. There are some very small differences between $D_1^*$ and $D_2^*$ that do not change the interpretation of the results, so in the tables and figures of this section we drop $D_2^*$. It is worth noting that the results in Section (ref) on the behaviour of JIVE1/JIVE2-based tests find confirmation in our simulations. Furthermore, by direct inspection of the respective formulae, it is easy to see that, in general, the proposed statistics are computationally simpler to evaluate than AR statistics mainly due to the Hadamard square in the expression of the variance. The complete set of results can be found in Appendix (ref). In what follows, we provide a description of the DGPs and some comments on the simulation results.
This simulation setup considers the model
Here, $\bm x_1$ and $\bm y$ are $n$-dimensional vectors with $n=200$. Moreover, $\bm X_2=\bm Z_2$ is a $(n\times g_2)$-dimensional matrix with $g=g_1+g_2=1+g_2$ being the number of regressors in equation (ref). We choose $g_2=5$; the first column of $\bm X_2$ contains a vector of ones and the remaining columns are sampled independently from a standard normal. The $n\times k_1$ instrument matrix $\bm Z_1$ contains in the first to third column the variables $\bm z_1$, $\bm z_1^{(2)}$ and $\bm z_1^{(3)}$ where $\bm z_1\sim\mathcal N(\boldsymbol{0},\bm I_n)$. The remaining columns are sampled from a standard normal. The number of instruments $k=k_1+g_2$ with $k_1=\alpha n$. Clearly, the parameter $\alpha$ is the share of instruments with respect to the sample size and is chosen as $\alpha\in\{0.05, 0.1\}$. The disturbance couple $(\bm \varepsilon, \bm v)$ is defined by the relationship
with $\bm u_1$ and $\bm u_2$ being mutually independent and distributed as a standard normal. We set $\beta=1$, $\bm \beta_2=\bm \iota_5$, $\rho_1=0.3$ and $\rho_2=0.2$. We also define $\bm \pi=\pi\bm \iota_k$ for $\bm \iota_k$ a scalar $\pi$ chosen as
with $r\in\{32,64\}$ being a measure of instrument strength. The null hypothesis corresponds to the linear restrictions
with $\bm \beta=(\beta,\bm \beta_2')'$, $\bm A=(1,\boldsymbol{0}')$, a row vector, and $\bm a=\beta_{0}=1$.
The model in the second DGP consists of two endogenous variables and reads
where $\bm X_3$ is a $n\times 2$ matrix with the first column being a vector of ones and the second being sampled from a standard normal, $\bm \varepsilon=0.2\bm u+\bm e$ with $\bm e\sim\mathcal N(\boldsymbol{0},0.3^2\bm I_n)$ and $\bm u\sim\mathcal N(\boldsymbol{0},\bm I_n)$. Furthermore,
where $\delta_1=0.5$, $\delta_2=0.4$, $\bm v_j\sim\mathcal N(\boldsymbol{0}, \sigma_{v_j}^2\bm I_n)$ with $\sigma_{v_j}=0.3$. Furthermore, $\bm Z_j\sim\mathcal N(\boldsymbol{0},\bm I_n)$. Finally, we have $\bm \pi_j=\sqrt{\frac{r}{k_j(1-r)}}\bm \iota_{k_j}$ and $k_j=\alpha n$, $n=200$ and $\alpha\in\{0.05,0.1\}$. In this case, the null hypothesis
has $\bm A=(1,1,\boldsymbol{0}')$, $\bm \beta=(\beta_1,\beta_2,\bm \beta_3')'$ and $\bm a=1$. The true values for the parameters of the endogenous variables are $\beta_1=0.3$ and $\beta_2=0.7$, while $\bm \beta_3=\bm \iota_2$.
The performance of the tests is measured in terms of rejection rates compared to a 5% nominal rate. The asymptotically chi-square distributed tests follow a $\chi^2_1$, while the estimated version of $\bm \varphi$ has a positive number in the first position and $g-1$ zeros in the remaining entries. The p values for the chi-bar-squared tests are computed using the {\tt R} function {\tt farebrother} DL10. For power, we focus on weak instruments. Namely, $r=32$ for DGP1 and $r=0.1$ for DGP2. The power plots for experiments with stronger instruments can be found in Appendix (ref) along with results based on tests that use cross-fit variances. As previously mentioned in this section, since some test statistics show the same type of behaviour, they are excluded from the plots in Figures (ref) to (ref). Also in this case, the full set of simulation results can be found in Appendix (ref).
\paragraph{Size} The size results in Table (ref) and Table (ref) show similar scenarios. In general we notice that size distortions oscillate between about $+3\%$ to about $-4\%$ with respect to the $5\%$ nominal size. The former result is observed for the $D$ statistic while the latter is observed for the $AR$ statistic with naive variance estimator. Both cases pertain to DGP1.
More specifically, we notice that the $LM$ and $LM^*$ tests show the same rejection rates. They consistently maintain good size properties relative to other statistics under the SJIVE/HLIM specification, particularly in the case of DGP1; under DGP2, they have a slight tendency to overreject. Wald statistics perform similarly under both asymptotic approximations, yet they tend to overreject under the SJIVE/HLIM objective functions. The opposite occurs when JIVE1 or JIVE2 is used. As in the case of Lagrange multiplier tests both Wald statistics show the same rejection rates. In contrast, statistics $D$ and $D^*$ display noticeable differences: $D$ tends to overreject relative to the nominal 5% level, whereas $D^*$ remains conservative, but close to the 5% nominal level. The proposed tests are also compared to Anderson-Rubin statistics $AR_{naive}$ and $AR_{cf}$ MS22, CMS21. We observe that the Anderson-Rubin statistics consistently underreject. For DGP1, the tests with cross-fit variance ($AR_{cf}$) remain conservative but show some size improvements with respect to the naive counterpart.
\paragraph{Power} The power results are shown in Figure (ref) to Figure (ref). The figures include results for tests based on the SJIVE and JIVE1 objective functions only, as they are similar to and, in some cases, improve over tests based on HLIM and JIVE2. For similar reasons, we exclude tests that are asymptotically distributed as a chi-bar-square. In our experiments, these power differences are noticeable but generally small (see figures in Appendix (ref)), with the exception of the Lagrange multiplier tests. It is known that Lagrange multiplier tests may lose power away from the true value.\footnote{Andrews16 reports this phenomenon in the context of combination of tests.} In our simulations, we observe a significant power drop for $LM$ tests based on both the SJIVE and HLIM objective functions and for both asymptotic distributions (see panels (a) and (b) of the figures in Appendix (ref)). $LM$ tests based on JIVE1 and JIVE2 show some power loss but not as large as the aforementioned cases (see panels (c) and (d) of the figures in Appendix (ref)). Interestingly, the power loss observed in the case of JIVE1 and in particular in the case of SJIVE is visibly smaller than the corresponding cases for JIVE2 and HLIM, respectively. Intuitively, this may be attributed to the different jackknife scheme used for SJIVE/JIVE1 tests compared to HLIM/JIVE2 tests.
Of the three types of statistics in Figure (ref) to Figure (ref), we notice that the Wald tests tend to have the largest power in general as we move away from the truth. The Lagrange multiplier tests, despite having good size properties in some cases, tend to have the worst power performance. The distance statistics seem to be somewhere in the middle, showing in a number of cases power properties comparable to those of the Wald tests and good size properties in the case of $D^*$.
Interestingly, while Wald statistics generally demonstrate good size properties and high power away from the true parameter, their power may be low near the true value. This is particularly evident for Wald tests based on JIVE1 and JIVE2 objective functions (see panels (b) and (d) in Figure (ref) to Figure (ref). This phenomenon is much less visible for SJIVE and HLIM objective functions. In general, the $D^*$ statistic appears to strike a good balance between size control and power. When compared to Anderson-Rubin statistics, apart from the pathological case of $LM$ and $LM^*$, our proposed tests generally show clear improvements in terms of power.
It is useful to contrast our results with some of the findings of the weak instrument literature developed in the 1990s and early 2000s (and related recent contributions). A key result is that conventional IV estimation and inference can be unreliable under weak identification. In particular, under different definitions of weakness -- whether near non-identification as in Dufour97 or local-to-zero as in SS97 -- Wald statistics fail to provide valid inference. This has set the ground for the development of alternative statistics that are robust to the presence of weak instruments. Under the local-to-zero asymptotic framework of SS97, WZ98 show that certain LM and LR statistics can be boundedly pivotal, yielding conservative yet asymptotically valid inference. A complementary approach is to construct statistics with pivotal limiting distributions as in the case of Kleibergen's K statistic Kleibergen02 that also address the limitations of the AR test in the overidentified case. More recently, KN23 highlight additional failures of 2SLS-based t-tests and recommend robust alternatives such as the AR test in exactly identified models and LIML estimation and CLR tests in overidentified settings. Our contribution can be also interpreted as a complement to this literature in a many instruments context.
We apply our proposed test statistics to examine the causal effect of alcohol consumption on body mass index (BMI) using data from the UK Biobank (UKB). The UK Biobank is a large-scale biomedical database containing extensive genetic and health information for approximately 500,000 individuals aged 40--69 in the United Kingdom Sudlow2015, Bycroft2018.
Our starting sample consists of 337,260 participants who were included in the genetic principal component analysis and classified as genetically Caucasian, selected from the 501,932 individuals available in the UK Biobank. Mendelian randomisation exploits genetic variants as instrumental variables to identify causal effects. Because genetic variants are randomly allocated at conception, they provide a source of exogenous variation that mimics experimental randomisation and can help overcome confounding and reverse causality that typically bias observational estimates of alcohol's health effects DaveySmith2003,Burgess2015. In this application, genetic variants associated with alcohol consumption serve as instruments for alcohol intake.
We initially extract 98 single nucleotide polymorphisms (SNPs) previously identified as associated with alcohol consumption. To reduce linkage disequilibrium and ensure instrument independence, we perform a pruning procedure that removes SNPs with pairwise correlation exceeding 0.15, yielding a final set of 91 instruments.
Following Topiwala2022, alcohol consumption is constructed as total weekly intake in grams of ethanol. Participants report consumption either weekly or monthly across several beverage categories, including red wine, champagne or white wine, beer or cider, spirits, fortified wine, and other alcoholic beverages. Monthly quantities are converted to weekly consumption by dividing by 4.3. Beverage-specific unit conversions are then applied, using factors such as 1.7 units per glass for wine, 2.4 units per pint for beer or cider, and one unit per measure of spirits. Alcohol units are converted to grams using the standard conversion of one UK unit corresponding to eight grams of ethanol.
The dataset undergoes several cleaning steps. Observations with missing or unknown alcohol consumption information are removed, resulting in the exclusion of 74,717 participants. Outliers in alcohol consumption are then removed using Tukey's fences with parameter $k = 1.5$, where the first and third quartiles are $Q_1 = 54.4$ and $Q_3 = 200$. This procedure excludes an additional 12,991 observations. After dropping all remaining observations with missing values, the final estimation sample consists of 97,987 individuals.
The final dataset used for estimation includes the outcome variable BMI, the endogenous variable total weekly alcohol consumption measured in grams per week, the 91 genetic variants used as instruments, and the following control variables: sex, age at recruitment, age squared, indicators for the UK Biobank assessment center, and the first ten genetic principal components to account for population stratification.
The results in Table (ref) show uniform and extremely strong rejection of the null hypothesis across all methods and test statistics. All JIVI-based procedures yield p values at or below $0.0002$, irrespective of whether inference is conducted using the $\bar{\chi}^2$ or the standard $\chi^2$ reference distribution. This uniformity indicates that inference in the full sample is largely insensitive to both the choice of statistic and the choice of asymptotic approximation.
Among the jackknife-based estimators, SJIVE and HLIM deliver virtually identical outcomes across all statistics. Similarly, JIVE1 and JIVE2 produce indistinguishable results, with p values of $0.0001$ for all trinity statistics and $0.0000$ for both Anderson--Rubin implementations.
Comparisons across statistics within each method reveal negligible differences between $D$, $W_1$, $W_2$, and $LM$, as well as their starred counterparts. Switching from the $\bar{\chi}^2$ distribution to the standard $\chi^2$ distribution has no material effect on inference.
Finally, the Anderson--Rubin tests (columns 12--13) yield p values of $0.0000$ for JIVE1 and JIVE2, reinforcing the same qualitative conclusion of strong rejection. Overall, the table documents a high degree of robustness across estimators, statistics, and reference distributions in the full sample.\footnote{Under the null, we use the OLS estimator for computational convenience. Although OLS is consistent, it does not preserve the finite-sample equality among the test statistics mentioned in Section (ref) and established in Appendix (ref); this discrepancy is purely numerical and has no asymptotic consequence.}
To assess whether the uniformity of inference in the full sample is driven by large-sample size, we repeat the analysis using a random subsample of 10{,}000 observations.
Table (ref) shows that, in contrast to the full sample, nontrivial differences across test statistics and reference distributions emerge when inference is based on a smaller sample.
Among the jackknife-based estimators, SJIVE and HLIM continue to display very similar behavior. Both yield small p values for the $D$ statistic, while the remaining statistics, $W_1$, $W_2$, and $LM$ produce larger but still statistically significant p values in the 0.01--0.03 range. Under the standard $\chi^2$ distribution, all starred statistics remain significant, with $LM^*$ delivering the smallest p-values for both methods.
For JIVE1 and JIVE2, dispersion across statistics is more pronounced. The weighted-$D$ statistic yields p values around 0.04, whereas $W_2$ and $LM$ are close to 0.05, implying less pronounced rejection under conventional thresholds. Under the standard $\chi^2$ distribution, $D_1^*$, $W_1^*$, and $W_2^*$ cluster around 0.04, while $LM^*$ remains around 0.01--0.02.
The Anderson--Rubin tests (columns 12--13) yield p values close to $0.0001$ for both JIVE1 and JIVE2, which are substantially smaller than those produced by the trinity statistics. As in the full sample, this illustrates that Anderson--Rubin and trinity procedures can generate materially different degrees of rejection in finite samples.
Overall, the sub-sample results indicate that the near-complete uniformity observed in the full sample does not extend to smaller samples, and that inference can depend on the particular statistic and reference distribution employed.
In analogy with classical likelihood theory, we introduce a unified approach to inference for the parameters of linear regression models in the presence of endogeneity, heteroskedasticity and many, potentially weak, instruments. The test statistics are based on objective functions commonly encountered in the modern literature on many instruments HNWCS12,BC15. Our tests can be used for inference on simple hypotheses as well as on hypotheses expressed in terms of linear restrictions. We show that their asymptotic distribution is chi-bar-square but a more conventional chi square limit can be found after appropriate modifications, a result that, in our opinion, is more appealing to practitioners. In an extensive simulation experiment we find that, generally, our methods tend to have better finite sample properties than competing Anderson-Rubin statistics both in terms of size and power.
We also illustrate the practical relevance of our methodology through an empirical application that studies the causal effect of alcohol consumption on body mass index using genetic variants as instrumental variables in the UK Biobank. This setting naturally gives rise to a large number of instruments, some of which may be weak, making it well suited for the framework developed in this paper. The application demonstrates how the proposed tests can be implemented in a realistic empirical environment with many potentially weak instruments and heteroskedastic disturbances. The results highlight the feasibility of our procedures in large-scale genetic datasets and illustrate how the proposed statistics can provide reliable inference in settings where standard IV methods may perform poorly.
It is reasonable to believe that our methodology can be adapted to models with many instruments with fixed effects as in CSW23 as well as to models with cluster dependence FLM25. A further important extension is to develop confidence sets by inverting our tests. This is especially relevant under weak identification, where informative procedures may yield confidence sets with nonstandard geometry (e.g., unbounded or disconnected regions). An important direction is to characterise these shapes by adapting the geometric results of DT05. Some of these problems are part of our ongoing research agenda.