Bootstrap Tests for Overidentification in Linear Regression Models

Abstract

We study the finite-sample properties of tests for overidentifying restrictions in linear regression models with a single endogenous regressor and weak instruments. Under the assumption of Gaussian disturbances, we derive expressions for a variety of test statistics as functions of eight mutually independent random variables and two nuisance parameters. The distributions of the statistics are shown to have an ill-defined limit as the parameter that determines the strength of the instruments tends to zero and as the correlation between the disturbances of the structural and reduced-form equations tends to plus or minus one. This makes it impossible to perform reliable inference near the point at which the limit is ill-defined. Several bootstrap procedures are proposed. They alleviate the problem and allow reliable inference when the instruments are not too weak. We also study their power properties.

1. Introduction

In recent years, there has been a great deal of work on the finite-sample properties of estimators and tests for linear regression models with endogenous regressors when the instruments are weak. Much of this work has focused on the case in which there is just one endogenous variable on the right-hand side, and numerous procedures for testing hypotheses about the coefficient of this variable have been studied. See, among many others, Staiger and Stock [1], Stock, Wright, and Yogo [2], Kleibergen [3], Moreira [4,5], Andrews, Moreira, and Stock [6], and Davidson and MacKinnon [7,8]. However, the closely related problem of testing overidentifying restrictions when the instruments are weak has been less studied. Staiger and Stock [1] provides some asymptotic results, and Guggenberger, Kleibergen, Mavroeidis, and Chen [9] derives certain asymptotic properties of tests for overidentifying restrictions in the context of hypotheses about the coefficients of a subset of right-hand side variables.

In the next section, we discuss the famous test of Sargan [10] and other asymptotic tests for overidentification in linear regression models estimated by instrumental variables (IV) or limited information maximum likelihood (LIML). We show that the test statistics are all functions of six quadratic forms defined in terms of the two endogenous variables of the model, the linear span of the instruments, and its orthogonal complement. In fact, they can be expressed as functions of a certain ratio of sums of squared residuals and are closely related to the test proposed by Anderson and Rubin [11]. In Section 3, we analyze the properties of these overidentification test statistics. We use a simplified model with only three parameters, which is nonetheless capable of generating statistics with exactly the same distributions as those generated by a more general model, and characterize the finite-sample distributions of the test statistics, not by providing closed-form solutions, but by providing unique recipes for simulating them in terms of just eight mutually independent random variables with standard distributions as functions of the three parameters.

In Section 4, we make use of these results in order to derive the limiting behavior of the statistics as the instrument strength tends to zero, as the correlation between the disturbances in the structural and reduced-form equations tends to unity, and as the sample size tends to infinity under the assumptions of weak-instrument asymptotics. Since the results of Section 4 imply that none of the tests we study is robust with weak instruments, Section 5 discusses a number of bootstrap procedures that can be used in conjunction with any of the overidentification tests. Some of these procedures are purely parametric, while others make use of resampling.

Then, in Section 6, we investigate by simulation the finite-sample behavior of the statistics we consider. Simulation evidence and theoretical analysis concur in strongly preferring a variant of a likelihood-ratio test to the more conventional forms of Sargan test. In Section 7, we look at the performance of bootstrap tests, finding that the best of them behave very well if the instruments are not too weak. However, as our theory suggests, they improve very little over tests based on asymptotic critical values in the neighborhood of the singularity that occurs where the instrument strength tends to zero and the correlation of the disturbances tends to one.

In Section 8, we analyze the power properties of the two main variants of bootstrap test. We obtain analytical results that generalize those of Section 3. Using those analytical results, we conduct extensive simulation experiments, mostly for cases that allow the bootstrap to yield reliable inference. We find that bootstrap tests based on IV estimation seem to have a slight power advantage over those based on LIML, at the cost of slightly greater size distortion under the null when the instruments are not too weak. Section 9 presents a brief discussion of how both the test statistics and the bootstrap procedures can be modified to take account of heteroskedasticity and clustered data. Finally, some concluding remarks are made in Section 10.

2. Tests for Overidentification

Although the tests for overidentification that we deal with are applicable to linear regression models with any number of endogenous right-hand side variables, we restrict attention in this paper to a model with just one such variable. We do so partly for expositional convenience and partly because this special case is of particular interest and has been the subject of much research in recent years. The model consists of just two equations,

$$\begin{array}{cc}\hfill {\mathit{y}}_1& =\beta {\mathit{y}}_2+\mathit{Z}\gamma +{\mathit{u}}_1,\text{and}\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill {\mathit{y}}_2& =\mathit{W}\mathit{\pi }+{\mathit{u}}_2.\hfill \end{array}$$

(2)

Here ${\mathit{y}}_1$ and ${\mathit{y}}_2$ are $n$-vectors of observations on endogenous variables, $\mathit{Z}$ is an $n\times k$ matrix of observations on exogenous variables, and $\mathit{W}$ is an $n\times l$ matrix of instruments such that $\mathcal{S}(\mathit{Z})\subset \mathcal{S}(\mathit{W})$, where the notation $\mathcal{S}(\mathit{A})$ means the linear span of the columns of the matrix $\mathit{A}$. The disturbances are assumed to be homoskedastic and serially uncorrelated. We assume that $l>k+1$, so that the model is overidentified.

The parameters of this model are the scalar β, the $k$-vector γ, the $l$-vector π, and the $2\times 2$ contemporaneous covariance matrix of the disturbances $u_{1i}$ and $u_{2i}$:

$$\Sigma \equiv \left[\begin{array}{cc}{\sigma }_1^2& \rho {\sigma }_1{\sigma }_2\\ \rho {\sigma }_1{\sigma }_2& {\sigma }_2^2\end{array}\right].$$

(3)

Equation (1) is the structural equation we are interested in, and Equation (2) is a reduced-form equation for the second endogenous variable ${\mathit{y}}_2$.

The model (1) and (2) implicitly involves one identifying restriction, which cannot be tested, and $q\equiv l-k-1$ overidentifying restrictions. These restrictions say, in effect, that if we append q regressors all belonging to $\mathcal{S}(\mathit{W})$ to Equation (1) in such a way that the equation becomes just identified, then the coefficients of these q additional regressors are zero.

The most common way to test the overidentifying restrictions is to use a Sargan test (Sargan, [10]), which can be computed in various ways. The easiest is probably to estimate Equation (1) by instrumental variables (IV), using the l columns of $\mathit{W}$ as instruments, so as to obtain the IV residuals

$${\widehat{\mathit{u}}}_1={\mathit{y}}_1-{\widehat{\beta }}_{\text{IV}}{\mathit{y}}_2-\mathit{Z}{\widehat{\gamma }}_{\text{IV}},$$

(4)

where ${\widehat{\beta }}_{\text{IV}}$ and ${\widehat{\gamma }}_{\text{IV}}$ denote the IV estimates of β and γ, respectively. The IV residuals ${\widehat{\mathit{u}}}_1$ are then regressed on $\mathit{W}$. The explained sum of squares from this regression divided by the IV estimate of ${\sigma }_1^2$ is the test statistic, and it is asymptotically distributed as ${\chi }^2(q)$.

The numerator of the Sargan statistic can be written as

$$({\mathit{y}}_1-\mathit{Z}{\widehat{\gamma }}_{\text{IV}}-{\widehat{\beta }}_{\text{IV}}{\mathit{y}}_2){}^{\top }{\mathit{P}}_{\mathit{W}}({\mathit{y}}_1-\mathit{Z}{\widehat{\gamma }}_{\text{IV}}-{\widehat{\beta }}_{\text{IV}}{\mathit{y}}_2),$$

(5)

where ${\widehat{\beta }}_{\text{IV}}$ and ${\widehat{\gamma }}_{\text{IV}}$ denote the IV estimates of β and γ, respectively, and ${\mathit{P}}_{\mathit{W}}\equiv \mathit{W}{(\mathit{W}{}^{\top }\mathit{W})}^{-1}\mathit{W}{}^{\top }$ projects orthogonally onto $\mathcal{S}(\mathit{W})$. We define ${\mathit{P}}_{\mathit{Z}}$ similarly, and let ${\mathit{M}}_{\mathit{W}}\equiv \mathit{I}-{\mathit{P}}_{\mathit{W}}$ and ${\mathit{M}}_{\mathit{Z}}\equiv \mathit{I}-{\mathit{P}}_{\mathit{Z}}$. Since $\mathit{Z}$ is orthogonal to the IV residuals,

$${\mathit{y}}_1-\mathit{Z}{\widehat{\gamma }}_{\text{IV}}-{\widehat{\beta }}_{\text{IV}}{\mathit{y}}_2={\mathit{M}}_{\mathit{Z}}({\mathit{y}}_1-\mathit{Z}{\widehat{\gamma }}_{\text{IV}}-{\widehat{\beta }}_{\text{IV}}{\mathit{y}}_2)={\mathit{M}}_{\mathit{Z}}({\mathit{y}}_1-{\widehat{\beta }}_{\text{IV}}{\mathit{y}}_2).$$

Then, since ${\mathit{P}}_{\mathit{W}}{\mathit{M}}_{\mathit{Z}}={\mathit{M}}_{\mathit{Z}}{\mathit{P}}_{\mathit{W}}={\mathit{P}}_{\mathit{W}}-{\mathit{P}}_{\mathit{Z}}={\mathit{M}}_{\mathit{Z}}-{\mathit{M}}_{\mathit{W}}$, the numerator of the Sargan statistic can also be written as

$$({\mathit{y}}_1-{\widehat{\beta }}_{\text{IV}}{\mathit{y}}_2){}^{\top }({\mathit{M}}_{\mathit{Z}}-{\mathit{M}}_{\mathit{W}})({\mathit{y}}_1-{\widehat{\beta }}_{\text{IV}}{\mathit{y}}_2).$$

(6)

Similarly, the denominator is just

$$\frac{1}{n}({\mathit{y}}_1-\mathit{Z}{\widehat{\gamma }}_{\text{IV}}-{\widehat{\beta }}_{\text{IV}}{\mathit{y}}_2){}^{\top }({\mathit{y}}_1-\mathit{Z}{\widehat{\gamma }}_{\text{IV}}-{\widehat{\beta }}_{\text{IV}}{\mathit{y}}_2)=\frac{1}{n}({\mathit{y}}_1-{\widehat{\beta }}_{\text{IV}}{\mathit{y}}_2){}^{\top }{\mathit{M}}_{\mathit{Z}}({\mathit{y}}_1-{\widehat{\beta }}_{\text{IV}}{\mathit{y}}_2).$$

(7)

In addition to being the numerator of the Sargan statistic, expression (6) is the numerator of the Anderson-Rubin, or AR, statistic for the hypothesis that $\beta ={\widehat{\beta }}_{\text{IV}}$; see Anderson and Rubin [11]. The denominator of this same AR statistic is

$$\frac{1}{n-l}({\mathit{y}}_1-{\widehat{\beta }}_{\text{IV}}{\mathit{y}}_2){}^{\top }{\mathit{M}}_{\mathit{W}}({\mathit{y}}_1-{\widehat{\beta }}_{\text{IV}}{\mathit{y}}_2),$$

(8)

which may be compared to the right-hand side of Equation (7). We see that the Sargan statistic estimates ${\sigma }_1^2$ under the null hypothesis, and the AR statistic estimates it under the alternative.

Of course, AR statistics are usually calculated for the hypothesis that β takes on a specific value, say ${\beta }_0$, rather than ${\widehat{\beta }}_{\text{IV}}$. Since by definition ${\widehat{\beta }}_{\text{IV}}$ minimizes the numerator (5), it follows that the numerator of the AR statistic is always no smaller than the numerator of the Sargan statistic. Even though the AR statistic is not generally thought of as a test of the overidentifying restrictions, it could be used as such a test, because it will always reject if the restrictions are sufficiently false, since that will cause the noncentrality parameter to be large.

It seems natural to modify the Sargan statistic by using (8) instead of (7) as the denominator, and this was done by Basmann [12]. The usual Sargan statistic can be written as

$$S=\frac{{\text{SSR}}_0-{\text{SSR}}_1}{{\text{SSR}}_0/n}=n\left(1-\zeta ({\widehat{\beta }}_{\text{IV}})\right)$$

(9)

and the Basmann statistic as

$$S^{\prime }=\frac{{\text{SSR}}_0-{\text{SSR}}_1}{{\text{SSR}}_1/(n-l)}=(n-l)\left({\zeta }^{-1}({\widehat{\beta }}_{\text{IV}})-1\right),$$

(10)

where ${\text{SSR}}_0$ is the sum of squared residuals from regressing ${\mathit{y}}_1-{\widehat{\beta }}_{\text{IV}}{\mathit{y}}_2$ on $\mathit{Z}$, ${\text{SSR}}_1$ is the SSR from regressing ${\mathit{y}}_1-{\widehat{\beta }}_{\text{IV}}{\mathit{y}}_2$ on $\mathit{W}$, and $\zeta ({\widehat{\beta }}_{\text{IV}})\equiv {\text{SSR}}_1/{\text{SSR}}_0$. Observe that both test statistics are simply monotonic functions of $\zeta ({\widehat{\beta }}_{\text{IV}})$, the ratio of the two sums of squared residuals.

Another widely used test statistic for overidentification is the likelihood ratio, or LR, statistic associated with the LIML estimate ${\widehat{\beta }}_{\text{LIML}}$. This statistic is simply $nlog\kappa ({\widehat{\beta }}_{\text{LIML}})$, where

$$\kappa (\beta )=\frac{({\mathit{y}}_1-\beta {\mathit{y}}_2){}^{\top }{\mathit{M}}_{\mathit{Z}}({\mathit{y}}_1-\beta {\mathit{y}}_2)}{({\mathit{y}}_1-\beta {\mathit{y}}_2){}^{\top }{\mathit{M}}_{\mathit{W}}({\mathit{y}}_1-\beta {\mathit{y}}_2)}.$$

(11)

The LIML estimator ${\widehat{\beta }}_{\text{LIML}}$ minimizes $\kappa (\beta )$ with respect to β. Since $\widehat{\kappa }\equiv \kappa ({\widehat{\beta }}_{\text{LIML}})$ is just ${\zeta }^{-1}({\widehat{\beta }}_{\text{LIML}})$, we see that the LR statistic is $-nlog\zeta ({\widehat{\beta }}_{\text{LIML}})$. It is easy to show that, under both conventional and weak-instrument asymptotics, and under the null hypothesis, $\widehat{\kappa }-1=O_p(n^{-1})$ as the sample size n tends to infinity. Therefore, the LR statistic is asymptotically equivalent to the linearized likelihood ratio statistic

$${\text{LR}}^{\prime }\equiv (n-l)(\widehat{\kappa }-1)=(n-l)\widehat{\lambda }=(n-l)\left({\zeta }^{-1}({\widehat{\beta }}_{\text{LIML}})-1\right)$$

(12)

where $\widehat{\lambda }\equiv \widehat{\kappa }-1$. We define LR${}^{\prime }$ as $(n-l)\widehat{\lambda }$ rather than as $n\widehat{\lambda }$ by analogy with (10). In what follows, it will be convenient to analyze LR${}^{\prime }$ rather than LR.

We have seen that the Sargan statistic (9), the Basmann statistic (10), and the two likelihood ratio statistics LR and LR${}^{\prime }$ are all monotonic functions of the ratio of SSRs $\zeta (\widehat{\beta })$ for some estimator $\widehat{\beta }$. Both the particular function of $\zeta (\widehat{\beta })$ that is used and the choice of $\widehat{\beta }$ affect the finite-sample properties of an asymptotic test. For a bootstrap test, however, it is only the choice of $\widehat{\beta }$ that matters. This follows from the fact that it is only the rank of the actual test statistic in the ordered list of the actual and bootstrap statistics that determines a bootstrap P value; see Section 6 below and Davidson and MacKinnon [13]. Therefore, for any given bootstrap data-generating process (DGP) and any estimator $\widehat{\beta }$, bootstrap tests based on any monotonic transformation of $\zeta (\widehat{\beta })$ yield identical results.

3. Analysis Using a Simpler Model

It is clear from (6), (7), and (11) that all the statistics we have considered for testing the overidentifying restrictions depend on ${\mathit{y}}_1$ and ${\mathit{y}}_2$ only through their projections ${\mathit{M}}_{\mathit{Z}}{\mathit{y}}_1$ and ${\mathit{M}}_{\mathit{Z}}{\mathit{y}}_2$. We see also that $\zeta (\beta )$ is homogeneous of degree zero with respect to ${\mathit{M}}_{\mathit{Z}}{\mathit{y}}_1$ and ${\mathit{M}}_{\mathit{Z}}{\mathit{y}}_2$ separately, for any β. Thus the statistics depend on the scale of neither ${\mathit{y}}_1$ nor ${\mathit{y}}_2$. Moreover, the matrix $\mathit{Z}$ plays no essential role. In fact, it can be shown—see Davidson and MacKinnon [7], Section 3—that the distributions of the test statistics generated by the model (1) and (2) for sample size n are identical to those generated by the simpler model

$$\begin{array}{cc}\hfill {\mathit{y}}_1& =\beta {\mathit{y}}_2+{\mathit{u}}_1,\text{and}\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill {\mathit{y}}_2& =\mathit{W}\mathit{\pi }+{\mathit{u}}_2,\hfill \end{array}$$

(14)

where the sample size is now $n-k$, the matrix $\mathit{W}$ has $l-k$ columns, and ${\sigma }_1={\sigma }_2=1$.

In the remainder of the paper, we deal exclusively with the simpler model (13) and (14). By doing so, we avoid the additional algebra made necessary by the presence of $\mathit{Z}$ in Equation (1) without changing the results in any essential way. Of course, ${\mathit{y}}_1$, ${\mathit{y}}_2$, and $\mathit{W}$ in the simpler model are not the same as in the original model, although ρ is unchanged. For the original model, n and l in our results below would have to be replaced by $n-k$ and $l-k$, and ${\mathit{y}}_1$ and ${\mathit{y}}_2$ would have to be replaced by ${\mathit{M}}_{\mathit{Z}}{\mathit{y}}_1$ and ${\mathit{M}}_{\mathit{Z}}{\mathit{y}}_2$.

It is well known—see Mariano and Sawa [14]—that all the test statistics depend on the data generated by (13) and (14) only through the six quadratic forms

$$\begin{array}{cc}& P_{11}\equiv {\mathit{y}}_1{}^{\top }{\mathit{P}}_{\mathit{W}}{\mathit{y}}_1,P_{12}\equiv {\mathit{y}}_1{}^{\top }{\mathit{P}}_{\mathit{W}}{\mathit{y}}_2,P_{22}\equiv {\mathit{y}}_2{}^{\top }{\mathit{P}}_{\mathit{W}}{\mathit{y}}_2,\hfill \\ & M_{11}\equiv {\mathit{y}}_1{}^{\top }{\mathit{M}}_{\mathit{W}}{\mathit{y}}_1,M_{12}\equiv {\mathit{y}}_1{}^{\top }{\mathit{M}}_{\mathit{W}}{\mathit{y}}_2,\text{and}{M}_{22}\equiv {\mathit{y}}_2{}^{\top }{\mathit{M}}_{\mathit{W}}{\mathit{y}}_2.\hfill \end{array}$$

(15)

This is also true for the general model (1) and (2), except that the projection matrix ${\mathit{P}}_{\mathit{W}}$ must be replaced by ${\mathit{P}}_{\mathit{W}}-{\mathit{P}}_{\mathit{Z}}={\mathit{P}}_{\mathit{W}}{\mathit{M}}_{\mathit{Z}}$.

In this section and the next, we make the additional assumption that the disturbances ${\mathit{u}}_1$ and ${\mathit{u}}_2$ are normally distributed. Since the quadratic forms in (15) depend on the instruments only through the projections ${\mathit{P}}_{\mathit{W}}$ and ${\mathit{M}}_{\mathit{W}}$, it follows that their joint distribution depends on $\mathit{W}$ only through the number of instruments l and the norm of the vector $\mathit{W}\mathit{\pi }$. We can therefore further simplify Equation (14) as

$${\mathit{y}}_2=a\mathit{w}+{\mathit{u}}_2,$$

(16)

where the vector $\mathit{w}\in \mathcal{S}(\mathit{W})$ is normalized to have length unity, which implies that $a^2=\mathit{\pi }{}^{\top }\mathit{W}{}^{\top }\mathit{W}\mathit{\pi }$. Thus the joint distribution of the six quadratic forms depends only on the three parameters β, a, and ρ, and on the dimensions n and l; for the general model (1) and (2), the latter would be $n-k$ and $l-k$.

The above simplification was used in Davidson and MacKinnon [7] in the context of tests of hypotheses about β, and further details can be found there. The parameter a determines the strength of the instruments. In weak-instrument asymptotics, $a=O(1)$, while in conventional strong-instrument asymptotics, $a=O(n^{1/2})$. Thus, by treating a as a parameter of order unity, we are in the context of weak-instrument asymptotics; see Staiger and Stock [1]. The square of the parameter a is often referred to as the (scalar) concentration parameter; see Phillips [15] (p. 470), and Stock, Wright, and Yogo [2].

Davidson and MacKinnon [7] show that, under the assumption of normal disturbances, the six quadratic forms (15) can be expressed as functions of the three parameters β, a, and ρ and eight mutually independent random variables, the distributions of which do not depend on any of the parameters. Four of these random variables, which we denote by $x_1$, $x_2$, $z_P$, and $z_M$, are standard normal, and the other four, which we denote by $t_{11}^P$, $t_{22}^P$, $t_{11}^M$, and $t_{22}^M$, are respectively distributed as ${\chi }_{l-2}^2$, ${\chi }_{l-1}^2$, ${\chi }_{n-l}^2$, and ${\chi }_{n-l-1}^2$. In terms of these eight variables, we make the definitions

$$\begin{array}{cc}& Q_{11}\equiv x_1^2+z_P^2+t_{11}^P,Q_{12}\equiv x_1{x}_2+z_P\sqrt{t_{22}^P},Q_{22}\equiv x_2^2+t_{22}^P,\hfill \\ & N_{11}\equiv t_{11}^M,N_{12}\equiv z_M\sqrt{t_{11}^M},\text{and}{N}_{22}\equiv z_M^2+t_{22}^M.\hfill \end{array}$$

(17)

These quantities have simple interpretations: $Q_{ij}={\mathit{u}}_i{}^{\top }{\mathit{P}}_{\mathit{W}}{\mathit{u}}_j$, and $N_{ij}={\mathit{u}}_i{}^{\top }{\mathit{M}}_{\mathit{W}}{\mathit{u}}_j$, for $i=1,2$.

Theorem 1. 

The Basmann statistic $S^{\prime }$ defined in (10) can be expressed as

$$S^{\prime }=\frac{(n-l)(P_{11}{P}_{22}-P_{12}^2)}{M_{11}{P}_{22}-2P_{12}{M}_{12}+P_{12}^2{M}_{22}/P_{22}},$$

(18)

and the statistic LR${}^{\prime }\equiv (n-l)\widehat{\lambda }$ defined in (12) can be expressed as

$${\text{LR}}^{\prime }=\frac{(n-l)(P_{11}{M}_{22}-2P_{12}{M}_{12}+P_{22}{M}_{11}-{\Delta }^{1/2})}{2(M_{11}{M}_{22}-M_{12}^2)},$$

(19)

where

$$\begin{array}{cc}\hfill P_{11}& =Q_{11},P_{12}=a{x}_1+\rho Q_{11}+r{Q}_{12},M_{11}=N_{11},\hfill \\ \hfill P_{22}& =a^2+2a(\rho x_1+r{x}_2)+{\rho }^2{Q}_{11}+2r\rho Q_{12}+r^2{Q}_{22},\hfill \\ \hfill M_{22}& ={\rho }^2{N}_{11}+2r\rho N_{12}+r^2{N}_{22},\text{and}{M}_{12}=\rho N_{11}+r{N}_{12},\hfill \end{array}$$

(20)

with $r=\sqrt{1-{\rho }^2}$, and

$$\Delta ={(P_{11}{M}_{22}-2P_{12}{M}_{12}+P_{22}{M}_{11})}^2-4(M_{11}{M}_{22}-M_{12}^2)(P_{11}{P}_{22}-P_{12}^2).$$

Proof: 

See the Appendix.

Remark 1. 

By substituting the relations (20) and (17) into (18) and (19), realizations of $S^{\prime }$ and LR${}^{\prime }$ can be generated as functions of realizations of the eight independent random variables. Since these have known distributions, it follows that Theorem 1 provides a complete characterization of the distributions of the statistics $S^{\prime }$ and LR${}^{\prime }$ in finite samples for arbitrary values of the parameters β, a, and ρ. Although it seems to be impossible to obtain closed-form solutions for their distribution functions, it is of course possible to approximate them arbitrarily well by simulation.

4. Particular Cases

In this section, we show that no test of the overidentifying restrictions is robust to weak instruments. In fact, the distributions of $S^{\prime }$ and LR${}^{\prime }$ have a singularity at the point in the parameter space at which $a=0$ and $\rho =\pm 1$, or, equivalently, $a=r=0$. By letting $n\to \infty$ for arbitrary values of a and r, we may obtain asymptotic results for both weak and strong instruments. In this way, we show that the singularity is present both in finite samples and under weak-instrument asymptotics. For conventional asymptotics, it is necessary to let both a and n tend to infinity, and in that case the singularity becomes irrelevant.

Like Theorem 1, the theorems in this section provide finite-sample results. Asymptotic results follow from these by letting the sample size n tend to infinity.

Theorem 2. 

When the parameter a is zero, the linearized likelihood ratio statistic LR${}^{\prime }$, in terms of the random variables defined in (17), is

$$\frac{(n-l)(Q_{11}{N}_{22}-2Q_{12}{N}_{12}+Q_{22}{N}_{11}-{\Delta }^{1/2})}{2(N_{11}{N}_{22}-N_{12}^2)},$$

(21)

where the discriminant Δ has become

$${(Q_{11}{N}_{22}-Q_{22}{N}_{11})}^2+4(Q_{12}{N}_{11}-Q_{11}{N}_{12})(Q_{12}{N}_{22}-Q_{22}{N}_{12}),$$

(22)

independent of r. When $r=0$, LR${}^{\prime }$ is

$$(n-l)(Q_{11}-x_1^2)/N_{11},$$

(23)

independent of a.

Proof. 

See the Appendix.

Remark 2.1. 

As with Theorem 1, the expressions (21)–(23) constitute a complete characterization of the distribution of LR${}^{\prime }$ in the limits when $a=0$ and when $r=0$.

Remark 2.2. 

The limiting value (21) for $a=0$ is independent of r. Thus the distribution of LR${}^{\prime }$ in the limit of completely irrelevant instruments is independent of all the model parameters.

Remark 2.3. 

The limiting value (23) for $r=0$ is independent of a, and it tends to a ${\chi }_{l-1}^2$ variable as $n\to \infty$.

Remark 2.4. 

The singularity mentioned above is a consequence of the fact that the limit at $a=r=0$ is ill-defined, since LR${}^{\prime }$ converges to two different random variables as $r\to 0$ for $a=0$ and as $a\to 0$ for $r=0$. These random variables are quite different and have quite different distributions.

Remark 2.5. 

Simulation results show that the limit to which LR${}^{\prime }$ converges as it approaches the singularity with $a=\alpha r$ for some positive α depends on α. In other words, the limit depends on the direction from which the singularity is reached.

Remark 2.6. 

The limit as $r\to 0$ may appear to be of little interest. However, it may have a substantial impact on bootstrap inference. Davidson and MacKinnon [7] show that a and r are not consistently estimated with weak instruments, and so, even if the true value of r is far from zero, the estimate needed for the construction of a bootstrap DGP (see Section 5) may be close to zero.

Remark 2.7. 

In the limit of strong instruments, the distribution of LR${}^{\prime }$ is the same as its distribution when $r=0$, as shown by the following corollary.

Corollary 2.1. 

The limit of LR${}^{\prime }$ as $a\to \infty$, which is the limit when the instruments are strong, is

$$(n-l)(Q_{11}-x_1^2)/N_{11}.$$

(24)

Proof of Corollary 2.1. 

Isolate the coefficients of powers of a in (19), and perform a Taylor expansion for small $1/a$. The result (24) follows easily. ☐

As $n\to \infty$, $N_{11}/n\to 1$, and so the result (24) confirms that the asymptotic distribution with strong instruments is just ${\chi }_{l-1}^2$. In fact, Guggenberger et al. [9] demonstrate that, for our simple model, the asymptotic distribution is bounded by ${\chi }_{l-1}^2$ however weak or strong the instruments may be. Some of the simulation results in Section 6 for finite samples are in accord with this result.

Theorem 3. 

The Basmann statistic $S^{\prime }$ in the limit with $a=0$ becomes

$$(n-l)\frac{(Q_{11}{Q}_{22}-Q_{12}^2)({\rho }^2{Q}_{11}+2r\rho Q_{12}+r^2{Q}_{22})}{{\rho }^2{D}_0+2r\rho D_1+r^2{D}_2},$$

(25)

where

$$\begin{array}{cc}\hfill D_0& =Q_{12}^2{N}_{11}-2Q_{11}{Q}_{12}{N}_{12}+Q_{11}^2{N}_{22},\hfill \\ \hfill D_1& =Q_{12}{Q}_{22}{N}_{11}-N_{12}(Q_{11}{Q}_{22}+Q_{12}^2)+Q_{11}{Q}_{12}{N}_{22},\text{and}\hfill \\ \hfill D_2& =Q_{22}^2{N}_{11}-2Q_{12}{Q}_{22}{N}_{12}+Q_{12}^2{N}_{22}.\hfill \end{array}$$

The limit with $r=0$ is

$$\frac{(n-l)(Q_{11}-x_1^2)(a^2+2a{x}_1+Q_{11})}{N_{11}{(a+x_1)}^2}.$$

(26)

Proof. 

The proof is omitted. It involves tedious algebra similar to that in the proof of Theorem 2; see the Appendix. ☐

Remark 3.1. 

The expression (25) for $a=0$ depends on r, unlike the analogous expression for LR${}^{\prime }$. When $r\to 0$ with $a=0$, it is easy to see that $S^{\prime }$ tends to the limit

$$(n-l)\frac{Q_{11}(Q_{11}{Q}_{22}-Q_{12}^2)}{Q_{12}^2{N}_{11}-2Q_{11}{Q}_{12}{N}_{12}+Q_{11}^2{N}_{22}}.$$

(27)

Remark 3.2. 

Similarly, the expression (26) for $r=0$ depends on a. Its limit as $a\to 0$ is

$$(n-l)\frac{Q_{11}(Q_{11}-x_1^2)}{N_{11}{x}_1^2},$$

(28)

which is quite different from (27), where the order of the limits is inverted. As with LR${}^{\prime }$, the limit as the singularity is approached depends on the direction of approach.

Remark 3.3. 

As expected, the limit of $S^{\prime }$ as $a\to \infty$ is the same as that of LR${}^{\prime }$.

The fact that the test statistics $S^{\prime }$ and LR${}^{\prime }$ depend on the parameters a and ρ indicates that these statistics are not robust to weak instruments. Passing to the limit as $n\to \infty$ with weak-instrument asymptotics does not improve matters. Of the six quadratic forms on which everything depends, only the $M_{ij}$ depend on n. Their limiting behavior is such that $M_{11}/n\to 1$, $M_{22}/n\to 1$, and $M_{12}/n\to \rho$ as $n\to \infty$. In contrast, the $P_{ij}$ do not depend on n, and they do depend on a and ρ.

5. Bootstrap Tests

Every test statistic has a distribution which depends on the DGP that generated the sample from which it is computed. The “true” DGP that generated an observed realization of the statistic is in general unknown. However, according to the bootstrap principle, one can perform inference by replacing the unknown DGP by an estimate of it, which is called the bootstrap DGP. Because what we need for inference is the distribution of the statistic under DGPs that satisfy the null hypothesis, the bootstrap DGP must necessarily impose the null. This requirement by itself does not lead to a unique bootstrap DGP, and we will see in this section that, for an overidentification test, there are several plausible choices.

If the observed value of a test statistic τ is $\widehat{\tau }$, and the rejection region is in the upper tail, then the bootstrap P value is the probability, under the bootstrap distribution of the statistic, that τ is greater than $\widehat{\tau }$. To estimate this probability, one generates a large number, say B, of realizations of the statistic using the bootstrap DGP. Let the jth realization be denoted by ${\tau }_j^{*}$. Then the simulation-based estimate of the bootstrap P value is just the proportion of the ${\tau }_j^{*}$ greater than $\widehat{\tau }$:

$${\widehat{p}}^{*}(\widehat{\tau })=\frac{1}{B}\sum _{j=1}^B\mathbb{I}({\tau }_j^{*}>\widehat{\tau }),$$

where $\mathbb{I}(·)$ is the indicator function, which is equal to 1 when its argument is true and 0 otherwise. If this fraction is smaller than α, the level of the test, then we reject the null hypothesis. See Davidson and MacKinnon [13].

5.1. Parametric Bootstraps

The DGPs contained in the simple model defined by Equation (13) and Equation (16) are characterized by just three parameters, namely, β, a, and ρ. Since the value of β does not affect the distribution of the overidentification test statistics (see Lemma 1 in the Appendix), the bootstrap DGP for a parametric bootstrap (assuming normally distributed disturbances) is completely determined by the values of a and ρ that characterize it.

The test statistic $\widehat{\tau }$ itself may be any one of the overidentification statistics we have discussed. The model that is actually estimated in order to obtain $\widehat{\tau }$ is not the simple model, but rather the full model given by Equations (1) and (2). The parameters of this model include some whose values do not interest us for the purpose of defining a bootstrap DGP: β, since it has no effect on the distribution of the statistic, and γ, since the matrix $\mathit{Z}$ plays no role in the simple model, from which the bootstrap DGP is taken. There remain π, ρ, ${\sigma }_1$, and ${\sigma }_2$.

For Equation (16), the parameter a was defined as the square root of $\mathit{\pi }{}^{\top }\mathit{W}{}^{\top }\mathit{W}\mathit{\pi }$. However, this definition assumes that the vector $\mathit{w}$ has unit length and that all the variables are scaled so that the variance of the disturbances ${\mathit{u}}_2$ is 1. A more general definition of a that does not depend on these assumptions is

$$a=\sqrt{\mathit{\pi }{}^{\top }\mathit{W}{}^{\top }\mathit{W}\mathit{\pi }/{\sigma }_2^2}.$$

(29)

It follows from (29) that, in order to estimate a, it is necessary also to estimate ${\sigma }_2^2$.

Since the parameter ρ is the correlation of the disturbances, which are not observed, any estimate of ρ must be based on the residuals from the estimation of Equations (1) and (2). Let these residuals be denoted by ${\ddot{\mathit{u}}}_1$ and ${\ddot{\mathit{u}}}_2$. Then the obvious estimators of the parameters of the covariance matrix are

$${\ddot{\sigma }}_1^2=n^{-1}{\ddot{\mathit{u}}}_1{}^{\top }{\ddot{\mathit{u}}}_1,{\ddot{\sigma }}_2^2=n^{-1}{\ddot{\mathit{u}}}_2{}^{\top }{\ddot{\mathit{u}}}_2,\text{and}\ddot{\rho }=n^{-1}{\ddot{\mathit{u}}}_1{}^{\top }{\ddot{\mathit{u}}}_2/{({\ddot{\sigma }}_1^2{\ddot{\sigma }}_2^2)}^{1/2},$$

and the obvious estimator of a is given by

$${\ddot{a}}^2=n\ddot{\mathit{\pi }}{}^{\top }\mathit{W}{}^{\top }\mathit{W}{}^{\top }\ddot{\mathit{\pi }}/{\ddot{\mathit{u}}}_2{}^{\top }{\ddot{\mathit{u}}}_2,$$

(30)

where $\ddot{\mathit{\pi }}$ estimates π. For ${\ddot{\mathit{u}}}_1$, there are two obvious choices, the IV residuals and the LIML residuals from (1). For ${\ddot{\mathit{u}}}_2$, the obvious choice is the vector of OLS residuals from (2), possibly scaled by a factor of ${(n/(n-l))}^{1/2}$ to take account of the lost degrees of freedom in the OLS estimation. However, this obvious choice is not the only one, because, if we treat the model (1) and (2) as a system, the system estimator of π that comes with the IV estimator of β is the three-stage least squares (3SLS) estimator, and the one that comes with the LIML estimator of β is the full-information maximum likelihood (FIML) estimator. These system estimators give rise to estimators not only of π, but also of ${\mathit{u}}_2$, that differ from those given by OLS.

The system estimators of π can be computed without actually performing a system estimation, by running the regression

$${\mathit{y}}_2=\mathit{W}\mathit{\pi }+\phi {\ddot{\mathit{u}}}_1\text{+}\text{residuals};$$

(31)

see Davidson and MacKinnon [7], where this matter is discussed in greater detail. If ${\ddot{\mathit{u}}}_1$ is the vector of IV residuals, then the corresponding estimator $\ddot{\mathit{\pi }}$ is the 3SLS estimator; if it is the vector of LIML residuals, then $\ddot{\mathit{\pi }}$ is the FIML estimator.

For the purpose of computation, it is worth noting that all these estimators can be expressed as functions of the six quadratic forms (15). A short calculation shows that the estimators of $a^2$ and ρ based on IV residuals, scaled OLS residuals, and the OLS estimator of π are

$${\widehat{a}}^2=(n-l)\frac{P_{22}}{M_{22}}\text{and}\widehat{\rho }=\frac{n^{-1}(M_{12}-\widehat{b}{M}_{22})}{{\widehat{\sigma }}_1{\widehat{\sigma }}_2},$$

(32)

where $\widehat{b}=P_{12}/P_{22}$ is the difference between the IV estimator of β and the true β of the DGP, and where

$${\widehat{\sigma }}_1^2=n^{-1}\left(Q_{11}+N_{11}-2\widehat{b}(P_{12}+M_{12})+{\widehat{b}}^2(P_{22}+M_{22})\right),\text{and}{\widehat{\sigma }}_2^2=M_{22}/(n-l).$$

From (17), we can express ${\widehat{a}}^2$ as

$$\frac{(n-l)(a^2+2a(\rho x_1+r{x}_2)+{\rho }^2{Q}_{11}+2r\rho Q_{12}+r^2{Q}_{22}}{{\rho }^2{N}_{11}+2r\rho N_{12}+r^2{N}_{22}}.$$

The weak-instrument asymptotic limit of this expression replaces the denominator divided by $n-l$ by 1. The expectation of the numerator without the factor of $n-l$ is $a^2+{\rho }^{2}l+r^{2}l=a^2+l$. Consequently, it may be preferable to reduce bias in the estimation of $a^2$ by setting ${\widehat{a}}^2=(n-l)max(0,P_{22}/M_{22}-l)$; see Davidson and MacKinnon [7].

The system estimator of the vector π, which we write as $\ddot{\mathit{\pi }}$, is computed by running regression (31). A tedious calculation shows that

$$\ddot{\mathit{\pi }}{}^{\top }\mathit{W}{}^{\top }\mathit{W}\ddot{\mathit{\pi }}=P_{22}+\frac{{\ddot{P}}_{11}{\ddot{M}}_{12}^2}{{\ddot{M}}_{11}^2}-2\frac{{\ddot{P}}_{12}{\ddot{M}}_{12}}{{\ddot{M}}_{11}},$$

(33)

with

$$\begin{array}{cc}\hfill {\ddot{P}}_{11}& ={\ddot{\mathit{u}}}_1{}^{\top }{\mathit{P}}_{\mathit{W}}{\ddot{\mathit{u}}}_1=P_{11}-2\ddot{b}{P}_{12}+{\ddot{b}}^2{P}_{22},\hfill \\ \hfill {\ddot{P}}_{12}& ={\ddot{\mathit{u}}}_1{}^{\top }{\mathit{P}}_{\mathit{W}}{\mathit{y}}_2=P_{12}-\ddot{b}{P}_{22},\hfill \\ \hfill {\ddot{M}}_{11}& ={\ddot{\mathit{u}}}_1{}^{\top }{\mathit{M}}_{\mathit{W}}{\ddot{\mathit{u}}}_1=M_{11}-2\ddot{b}{M}_{12}+{\ddot{b}}^2{M}_{22},\text{and}\hfill \\ \hfill {\ddot{M}}_{12}& ={\ddot{\mathit{u}}}_1{}^{\top }{\mathit{M}}_{\mathit{W}}{\mathit{y}}_2=M_{12}-\ddot{b}{M}_{22},\hfill \end{array}$$

where $\ddot{b}$ is $\widehat{b}={\widehat{\beta }}_{\text{IV}}-\beta$ and ${\ddot{\mathit{u}}}_1$ is the vector of IV residuals if the structural Equation (1) is estimated by IV, and $\ddot{b}$ is $\tilde{b}={\widehat{\beta }}_{\text{LIML}}-\beta$ and ${\ddot{\mathit{u}}}_1$ is the vector of LIML residuals if (1) is estimated by LIML. Then, if ${\sigma }_2^2$ is estimated using the residuals ${\mathit{y}}_2-\mathit{W}\ddot{\mathit{\pi }}$, we find that

$${\ddot{\sigma }}_2^2=\frac{1}{n-l}\left(M_{22}+\frac{{\ddot{M}}_{12}^2}{{\ddot{M}}_{11}^2}{\ddot{P}}_{11}\right).$$

(34)

Using the results (33) and (34) in (30) allows us to write

$${\ddot{a}}^2=\frac{(n-l)(P_{22}{\ddot{M}}_{11}^2+{\ddot{P}}_{11}{\ddot{M}}_{12}^2-2{\ddot{P}}_{12}{\ddot{M}}_{11}{\ddot{M}}_{12})}{{\ddot{M}}_{11}^2{M}_{22}+{\ddot{M}}_{12}^2{\ddot{P}}_{11}}.$$

For the parameter ρ, more calculation shows that

$$\ddot{\rho }={\ddot{M}}_{12}{\left[\frac{{\ddot{P}}_{11}+{\ddot{M}}_{11}}{{\ddot{M}}_{11}^2{M}_{22}+{\ddot{M}}_{12}^2{\ddot{P}}_{11}}\right]}^{1/2}.$$

We consider four different bootstrap DGPs. The simplest, which we call the IV-R bootstrap, uses the IV and OLS estimates of a and ρ given in (32). The IV-ER bootstrap is based on 3SLS estimation of the two-equation system, that is, on IV estimation of (1) and OLS estimation of regression (31) with ${\ddot{\mathit{u}}}_1$ the vector of IV residuals. Similarly, the LIML-ER bootstrap relies on FIML estimation of the system, that is, on LIML for (1) and on regression (31) with ${\ddot{\mathit{u}}}_1$ the LIML residuals. Finally, we also define the F(1)-ER bootstrap, which is the same as LIML-ER except that ${\widehat{\beta }}_{\text{LIML}}$ is replaced at every step by the modified LIML estimator of Fuller [16] with $\eta =1$ which is discussed in Section 6.

It is plain that, the closer the bootstrap DGP is to the true DGP, the better will be bootstrap inference; see Davidson and MacKinnon [17]. We may therefore expect that IV-ER should perform better than IV-R, and that LIML-ER should perform better than IV-ER. Between LIML-ER and F(1)-ER, there is no obvious reason a priori to expect that one of them would outperform the other. But, whatever the properties of these bootstraps may be when the true DGP is not in the neighborhood of the singularity at $a=0$, $\rho =1$, we cannot expect anything better than some improvement over inference based on asymptotic critical values, rather than truly reliable inference, in the neighborhood of the singularity.

5.2. Resampling

Any parametric bootstrap risks being unreliable if the strong assumptions used to define the null hypothesis are violated. Most practitioners would therefore prefer a more robust bootstrap method. The strongest assumption we have made so far is that the disturbances are normally distributed. It is easy to relax this assumption by using a bootstrap DGP based on resampling, in which the bivariate normal distribution is replaced by the joint empirical distribution of the residuals. The discussion of the previous subsection makes it clear that several resampling bootstraps can be defined, depending on the choice of residuals that are resampled.

The most obvious resampling bootstrap DGP in the context of IV estimation is

$$\begin{array}{cc}\hfill {\mathit{y}}_1^{*}& ={\widehat{\beta }}_{\text{IV}}{\mathit{y}}_2^{*}+\mathit{Z}{\widehat{\gamma }}_{\text{IV}}+{\widehat{\mathit{u}}}_1^{*}\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill {\mathit{y}}_2^{*}& =\mathit{W}\widehat{\mathit{\pi }}+{\widehat{\mathit{u}}}_2^{*},\hfill \end{array}$$

(36)

where ${\mathit{y}}_1^{*}$ and ${\mathit{y}}_2^{*}$ are $n$-vectors of bootstrap observations, ${\widehat{\mathit{u}}}_1^{*}$ and ${\widehat{\mathit{u}}}_2^{*}$ are $n$-vectors of bootstrap disturbances with typical elements ${\widehat{u}}_{1i}^{*}$ and ${\widehat{u}}_{2i}^{*}$, respectively, and $\widehat{\mathit{\pi }}$ is the OLS estimate from Equation (2). The bootstrap disturbances are drawn in pairs from the bivariate empirical distribution of the structural residuals ${\widehat{u}}_{1i}^{\text{IV}}$ and the rescaled reduced-form residuals ${\left(n/(n-l)\right)}^{1/2}{\widehat{u}}_{2i}^{\text{OLS}}$. That is,

$$\left[\begin{array}{c}{\widehat{u}}_{1i}^{*}\\ {\widehat{u}}_{2i}^{*}\end{array}\right]\sim \text{EDF}\left(\begin{array}{c}{\widehat{u}}_{1i}^{\text{IV}}\\ {\left(n/(n-l)\right)}^{1/2}{\widehat{u}}_{2i}^{\text{OLS}}\end{array}\right).$$

(37)

Here EDF stands for “empirical distribution function”.

The rescaling of the reduced form residuals ${\widehat{u}}_{2i}^{\text{OLS}}$ in (37) ensures that the distribution of the ${\widehat{u}}_{2i}^{*}$ has variance equal to the unbiased OLS variance estimator. The fact that the bootstrap disturbances are drawn in pairs ensures that the bootstrap DGP mimics the correlation between the two sets of residuals. Of course, whether that correlation consistently estimates the correlation between the actual disturbances depends on whether or not $\widehat{\rho }$ can be estimated consistently. Under weak-instrument asymptotics, $\widehat{\rho }$ cannot be estimated consistently; see Davidson and MacKinnon [7].

Since all of the overidentification test statistics are invariant to the values of β and γ, we may replace the bootstrap DGP for ${\mathit{y}}_1^{*}$ given by (35) by

$${\mathit{y}}_1^{*}={\widehat{\mathit{u}}}_1^{*}.$$

(38)

The bootstrap statistics generated by (38) and (36) are identical to those generated by (35) and (36). We will refer to the bootstrap DGP given by (36)–(38), as the IV-R resampling bootstrap. It is a semiparametric bootstrap, because it uses parameter estimates of the reduced-form equation, but it does not assume a specific functional form for the joint distribution of the disturbances. The empirical distribution of the residuals has a covariance matrix which is exactly that used to estimate a and ρ by the IV-R parametric bootstrap; hence our nomenclature.

The IV-ER resampling bootstrap draws pairs from the joint EDF of the IV residuals ${\widehat{\mathit{u}}}_1^{\text{IV}}$ from Equation (1) and the residuals ${\mathit{y}}_2-\mathit{W}{\widehat{\mathit{\pi }}}_{\text{IV}}$ computed by running regression (31) with ${\widehat{\mathit{u}}}_1^{\text{IV}}$ replacing ${\ddot{\mathit{u}}}_1$. It also uses the resulting estimator ${\widehat{\mathit{\pi }}}_{\text{IV}}$ in (36) instead of the OLS estimator $\widehat{\mathit{\pi }}$. Note that the residuals ${\mathit{y}}_2-\mathit{W}{\widehat{\mathit{\pi }}}_{\text{IV}}$ are not the residuals from regression (31), but rather those residuals plus $\ddot{\phi }{\widehat{\mathit{u}}}_1^{\text{IV}}$.

The LIML-ER resampling bootstrap is very similar to the IV-ER one, except that it uses ${\widehat{\mathit{u}}}_1^{\text{LIML}}$ both directly and in regression (31). Formally, the resampling procedure draws pairs from the bivariate empirical distribution of

$$\left[\begin{array}{c}{\widehat{u}}_{i1}^{\text{LIML}}\\ {\widehat{u}}_{i2}^{\text{LIML}}\end{array}\right]=\left[\begin{array}{c}{y}_{i1}-{\widehat{\beta }}_{\text{LIML}}{y}_{2i}-{\mathit{Z}}_i{\widehat{\gamma }}_{\text{LIML}}\\ y_{i2}-{\mathit{W}}_i{\widehat{\mathit{\pi }}}_{\text{LIML}}\end{array}\right].$$

(39)

Similarly, for the F(1)-ER resampling bootstrap, the structural Equation (1) is estimated by Fuller’s method, and the residuals from this are used both for resampling and in regression (31); see the next section.

A word of caution is advisable here. Although the values of overidentification test statistics are invariant to β, thereby allowing us to use (38) instead of (35) in the bootstrap DGP, the residuals from which we resample in (37) and (39) do depend on the estimate of β, as does the estimate of π if it is based on any variant of Equation (31). Nevertheless, the test statistics depend on the estimate of β only through the residuals and the estimate of π.

6. Performance of Asymptotic Tests: Simulation Results

The discussion in Section 4 was limited to the statistics $S^{\prime }$ and LR${}^{\prime }$. For bootstrap tests, it is enough to consider just these two, since all other statistics mentioned in Section 2 are monotonic transforms of them. Of course, the different versions of the Sargan test and the LR test have different properties when used with (strong-instrument) asymptotic critical values. In this section, we therefore present some Monte Carlo results on the finite-sample performance of five test statistics, including S, $S^{\prime }$, LR, and LR${}^{\prime }$.

The fifth test statistic we examine is based on the estimator proposed by Fuller [16]. Like the IV and LIML estimators, Fuller’s estimator is a K-class estimator for model (1) and (2). It takes the form

$$\left[\begin{array}{c}\widehat{\beta }\\ \widehat{\gamma }\end{array}\right]=(\mathit{X}{}^{\top }(\mathit{I}-K{\mathit{M}}_{\mathit{W}}){\mathit{X})}^{-1}\mathit{X}{}^{\top }(\mathit{I}-K{\mathit{M}}_{\mathit{W}}){\mathit{y}}_1.$$

(40)

Setting $K=\widehat{\kappa }$, the minimized value of the variance ratio (11), in Equation (40) gives the LIML estimator, while setting $K=1$ gives the IV estimator. Fuller’s estimator sets $K=\widehat{\kappa }+\eta /(n-l)$ for some fixed number $\eta >0$ independent of the sample size n. We set $\eta =1$. With this choice, Fuller’s estimator ${\widehat{\beta }}_{\mathrm{F}}$ has all moments (except when the sample size is very small) and is approximately unbiased. The corresponding test statistic is simply $-nlog\zeta ({\widehat{\beta }}_{\mathrm{F}})$, which has the same form as the LR statistic. We will refer to this as the LRF test.

The DGPs used for our simulations all belong to the simplified model (13) and (16). The disturbances are generated according to the relations

$${\mathit{u}}_1={\mathit{v}}_1\text{and}{\mathit{u}}_2=\rho {\mathit{v}}_1+r{\mathit{v}}_2,$$

where ${\mathit{v}}_1$ and ${\mathit{v}}_2$ are $n$-vectors with independent standard normal elements, and $r\equiv {(1-{\rho }^2)}^{1/2}$. Of course, it is quite unnecessary to generate simulated samples of n observations, as it is enough to generate the six quadratic forms (15) as functions of eight mutually independent random variables, using the relations (17) and (20). The sample size n affects only the degrees of freedom of the two ${\chi }^2$ random variables $t_{11}^M$ and $t_{22}^M$ that appear in (17). Although any DGP given by (13) and (16) involves no explicit overidentifying restrictions, the test statistics are computed for the model (1) and (2), for which there are $q\equiv l-k-1$ of them, with $k=0$ in the experiments.

The first group of experiments is intended to provide guidance on the appropriate sample size to use in the remaining experiments. Our objective is to mimic the common situation in which the sample size is reasonably large and the instruments are quite weak. Since our simulation DGPs embody weak-instrument asymptotics, we should not expect any of the test statistics to have the correct size asymptotically. However, for any given a and ρ, the rejection frequency converges as $n\to \infty$ to that given by the asymptotic distribution of the statistic used; these asymptotic distributions were discussed at the end of Section 4. In the experiments, we use sample sizes of 20, 28, 40, 56, and so on, up to 1810. Each of these numbers is larger than its predecessor by approximately $\sqrt{2}$. Each experiment used ${10}^6$ replications.

The results of four sets of experiments are presented in Figure 1, in which we plot rejection frequencies in the experiments for a nominal level of 5%. In the top two panels, $a=2$, so that the instruments are very weak. In the bottom two panels, $a=8$, so that they are reasonably strong. Recall that the concentration parameter is $a^2$. In the two panels on the left, $\rho =0.5$, so that there is moderate correlation between the structural and reduced form disturbances. In the two panels on the right, $\rho =0.9$, so that there is strong correlation. Note that the vertical axis differs across most of the panels.

Figure 1. Rejection frequencies for asymptotic tests as functions of n for $q=8$. (a) $a=2$, $\rho =0.5$; (b) $a=2$, $\rho =0.9$; (c) $a=8$, $\rho =0.5$; (d) $a=8$, $\rho =0.9$. Note that the scale of the vertical axis differs across panels.

Figure 1. Rejection frequencies for asymptotic tests as functions of n for $q=8$. (a) $a=2$, $\rho =0.5$; (b) $a=2$, $\rho =0.9$; (c) $a=8$, $\rho =0.5$; (d) $a=8$, $\rho =0.9$. Note that the scale of the vertical axis differs across panels.

It is evident that the various asymptotic tests can perform very differently and that their performance varies greatly with the sample size. The Basmann test $S^{\prime }$ always rejects more often than LR${}^{\prime }$. This must be the case, because the two test statistics are exactly the same function of ζ—compare Equation (10) and Equation (12)—but the LIML estimator minimizes ${\zeta }^{-1}$ and the IV estimator does not. Thus ${\zeta }^{-1}({\widehat{\beta }}_{\text{LIML}})$ must be closer to unity than ${\zeta }^{-1}({\widehat{\beta }}_{\text{IV}})$.

The Sargan (S) and Basmann ($S^{\prime }$) tests perform almost the same for large samples but very differently for small ones, with the latter much more prone to overreject than the former. This could have been predicted from Equation (9) and Equation (10). If it were not for the different factors of n and $n-l$, $S^{\prime }$ would always be greater than S, because ${\zeta }^{-1}-1>1-\zeta$. For $\zeta \approx 1$, the difference between ${\zeta }^{-1}-1$ and $1-\zeta$ is $O({\zeta }^{-2})$. Therefore, under the null hypothesis, the difference between the two test statistics is $O(n^{-1})$. This explains their rapid convergence as n increases.

For $a=2$, the LR test and its linearized version LR${}^{\prime }$ perform somewhat differently in small samples but almost identically once $n\ge 200$, with LR rejecting more often than LR${}^{\prime }$. The difference between the two LIML-based statistics is much smaller than the difference between the two IV-based statistics for two reasons. First, because the LIML estimator minimizes ${\zeta }^{-1}$, the two LR test statistics are being evaluated at values of ζ that are closer to unity than the two IV-based statistics. Second, $({\zeta }^{-1}-1)+log(\zeta )$ is almost exactly half the magnitude of $({\zeta }^{-1}-1)-(1-\zeta )$. Based on the fact that $({\zeta }^{-1}-1)>-log(\zeta )$, we might expect LR${}^{\prime }$ to reject more often than LR. However, the difference between $n-l$ and n evidently more than offsets this, at least for $l=9$, the case in the figure.

In this case, the Fuller variant of the LR test performs somewhat differently from both LR and LR${}^{\prime }$ for all sample sizes. In contrast, for $a=8$, LR and LRF are so similar that we did not graph LR to avoid making the figure unreadable. LR, LR${}^{\prime }$, and LRF perform almost identically, and very well indeed, for large sample sizes, even though they overreject severely for small sample sizes.

As expected, all of the rejection frequencies seem to be converging to constants as $n\to \infty$. Moreover, in every case, it appears that the (interpolated) results for $n=400$ are very similar to the results for larger values up to $n=1810$. Accordingly, we used $n=400$ in all the remaining experiments.

In the second group of experiments, the number of overidentifying restrictions q is varied. The four panels in Figure 2 correspond to those of Figure 1. In most cases, performance deteriorates as q increases. Sometimes, rejection frequencies seem to be converging, but by no means always. In the remaining experiments, we somewhat arbitrarily set $q=8$. Choosing a smaller number would generally have resulted in smaller size distortions.

In the third group of experiments, the results of which are shown in Figure 3, we set $n=400$ and $q=8$, and we vary ρ between 0.0 and 0.99 at intervals of 0.01 for four values of a. The vertical axis is different in each of the four panels, because the tests all perform much better as a increases. For clarity, rejection frequencies for LR are not shown in the figure, because they always lie between those for LR${}^{\prime }$ and LRF. They are very close to those for LR${}^{\prime }$ when a is small, and they are very close to those for LRF when a is large.

For the smaller values of a, all of the tests can either overreject or underreject, with rejection frequencies increasing in ρ. The Sargan and Basmann tests overreject very severely when a is small and ρ is large. The LR${}^{\prime }$, LR, and LRF tests underreject severely when a is small and ρ is not large, but they overreject slightly when a is large. Based on Figure 1 and on the analysis of the previous section, we expect that this slight overrejection vanishes for larger samples.

Figure 2. Rejection frequencies for asymptotic tests as functions of q for $n=400$. (a) $a=2$, $\rho =0.5$; (b) $a=2$, $\rho =0.9$; (c) $a=8$, $\rho =0.5$; (d) $a=8$, $\rho =0.9$. Note that the scale of the vertical axis differs across panels.

Figure 2. Rejection frequencies for asymptotic tests as functions of q for $n=400$. (a) $a=2$, $\rho =0.5$; (b) $a=2$, $\rho =0.9$; (c) $a=8$, $\rho =0.5$; (d) $a=8$, $\rho =0.9$. Note that the scale of the vertical axis differs across panels.

Figure 3. Rejection frequencies for asymptotic tests as functions of ρ for $q=8$ and $n=400$. (a) Very weak instruments: $a=2$; (b) Weak instruments: $a=4$; (c) Moderately strong instruments: $a=8$; (d) Very strong instruments: $a=16$. Note that the scale of the vertical axis differs across panels.

Figure 3. Rejection frequencies for asymptotic tests as functions of ρ for $q=8$ and $n=400$. (a) Very weak instruments: $a=2$; (b) Weak instruments: $a=4$; (c) Moderately strong instruments: $a=8$; (d) Very strong instruments: $a=16$. Note that the scale of the vertical axis differs across panels.

Although the performance of all the overidentification tests is quite poor when a is small, it is worth noting that the Sargan tests are not as unreliable as t tests of the hypothesis that β has a specific value, and the LR tests for overidentification are not as unreliable as LR tests for that hypothesis; see Davidson and MacKinnon [7,8].

6.1. Near the Singularity

From Figure 1, Figure 2 and Figure 3, we see that the rejection probabilities of all the tests vary considerably with the parameters a and ρ as they vary in the neighborhood of the singularity at $a=0$, $\rho =1$. Further insight into this phenomenon is provided by Figure 4 and Figure 5. These are contour plots of rejection frequencies near the singularity for tests at the 0.05 level with a and ρ on the horizontal and vertical axes, respectively. Figure 4 is for the Basmann statistic $S^{\prime }$, and Figure 5 is for the LR${}^{\prime }$ statistic. Both figures are for the case dealt with in Figure 3, for which $n=400$ and $q=8$. The rejection frequencies are, once again, estimated using ${10}^6$ replications.

Figure 4. Contours of rejection frequencies for $S^{\prime }$ tests with $q=8$ and $n=400$.

Figure 4. Contours of rejection frequencies for $S^{\prime }$ tests with $q=8$ and $n=400$.

Figure 5. Contours of rejection frequencies for likelihood ratio (LR)${}^{\prime }$ tests with $q=8$ and $n=400$.

Figure 5. Contours of rejection frequencies for likelihood ratio (LR)${}^{\prime }$ tests with $q=8$ and $n=400$.

It is clear from these figures that rejection frequencies tend to be greatest as the singularity is approached by first setting $r=0$ and then letting a tend to zero. In this limit, $S^{\prime }$ is given by expression (28), and LR${}^{\prime }$ is given by expression (23). For extremely small values of a, $S^{\prime }$ actually underrejects. But, as a rises to values that are still very very small, rejection frequencies soar, sometimes to over 0.80. In contrast, LR${}^{\prime }$ underrejects severely for small values of a, values which do not have to be nearly as small as in the case of $S^{\prime }$. In much of the figure, however, the rejection frequencies for LR${}^{\prime }$ are just a little greater than 0.05.

The 95% quantile of the distribution of expression (28) has the huge value of 16,285, as estimated from 9,999,999 independent realizations. In contrast, recall that the 95% quantile of the ${\chi }_q^2$ distribution for $q=8$ is 15.5073. Since the distribution of $S^{\prime }$ for arbitrary a and ρ is stochastically bounded by that of (28), $S^{\prime }$ is boundedly pivotal. However, basing inference on the distribution of (28) is certain to be extremely conservative.

7. Performance of Bootstrap Tests

In principle, any of the bootstrap DGPs discussed in the previous section can be combined with any of the test statistics discussed in Section 2. However, there is no point considering both S and $S^{\prime }$, or both LR and LR${}^{\prime }$, because in each case one test statistic is simply a monotonic transformation of the other. If both the statistics in each pair are bootstrapped using the same bootstrap DGP, they must therefore yield identical results.

All of our experiments involve 100,000 replications for each set of parameter values, and the bootstrap tests mostly use $B=399$. This is a smaller number than should generally be used in practice, but it is perfectly satisfactory for simulation experiments, because experimental randomness in the bootstrap p values tends to average out across replications. Although the disturbances of the true DGPs are taken to be normally distributed, the bootstrap DGPs we investigate in the main experiments are resampling ones, because we believe they are the ones that will be used in practice.

Figure 6, Figure 7 and Figure 8 present the results of a large number of Monte Carlo experiments. Figure 6 concerns Sargan tests, Figure 7 concerns LR tests, and Figure 8 concerns Fuller LR tests. Each of the figures shows rejection frequencies as a function of ρ for 34 values of ρ, namely, 0.00, 0.03, 0.06, …, 0.99. The four panels correspond to $a=2$, 4, 6, and 8. Note that the scale of the vertical axis often differs across panels within each figure and across figures for panels corresponding to the same value of a. It is important to keep this in mind when interpreting the results.

As we have already seen, for small and moderate values of a, Sargan tests tend to overreject severely when ρ is large and to underreject modestly when it is small. It is evident from Figure 6 that, for $a=2$, using either the IV-R or IV-ER bootstrap improves matters only slightly. However, both these methods do provide a more and more noticeable improvement as a increases. For $a=8$, the improvement is very substantial. If we were increasing n as well as a, it would be natural to see this as evidence of an asymptotic refinement.

Figure 6. Rejection frequencies for Sargan tests as functions of ρ for $q=8$ and $n=400$. (a) $a=2$; (b) $a=4$; (c) $a=8$; (d) $a=16$. Note that the scale of the vertical axis differs across panels.

Figure 6. Rejection frequencies for Sargan tests as functions of ρ for $q=8$ and $n=400$. (a) $a=2$; (b) $a=4$; (c) $a=8$; (d) $a=16$. Note that the scale of the vertical axis differs across panels.

There seems to be no advantage to using IV-ER rather than IV-R. In fact, the latter always works a bit better when ρ is very large. This result is surprising in the light of the findings of Davidson and MacKinnon [7,8] for bootstrapping t tests on β. However, the bootstrap methods considered in those papers imposed the null hypothesis that $\beta ={\beta }_0$, while the ones considered here do not. Apparently, this makes a difference.

Using the LIML-ER and F(1)-ER bootstraps with the Sargan statistic yields entirely different results. The former underrejects very severely for all values of ρ when a is small, but the extent of the underrejection drops rapidly as a increases. The latter always underrejects less severely than LIML-ER (it actually overrejects for large values of ρ when $a=2$), and it performs surprisingly well for $a\ge 6$. Of course, it may seem a bit strange to bootstrap a test statistic based on IV estimation using a bootstrap DGP based on LIML or its Fuller variant.

Figure 7. Rejection frequencies for LR tests as functions of ρ for $q=8$ and $n=400$. (a) $a=2$; (b) $a=4$; (c) $a=8$; (d) $a=16$. Note that the scale of the vertical axis differs across panels.

Figure 7. Rejection frequencies for LR tests as functions of ρ for $q=8$ and $n=400$. (a) $a=2$; (b) $a=4$; (c) $a=8$; (d) $a=16$. Note that the scale of the vertical axis differs across panels.

In Figure 7, we see that, in contrast to the Sargan test, the LR test generally underrejects, often very severely when both ρ and a are small. Its performance improves rapidly as a increases, however, and it actually overrejects slightly when ρ and a are both large. All of the bootstrap methods improve matters, and the extent of the improvement increases with a. For $a=8$, all the bootstrap methods work essentially perfectly. For small values of a, the IV-R bootstrap actually seems to be the best in many cases, although it does lead to modest overrejection when ρ is large.

Figure 8. Rejection frequencies for Fuller LR tests as functions of ρ for $q=8$ and $n=400$. (a) $a=2$; (b) $a=4$; (c) $a=8$; (d) $a=16$. Note that the scale of the vertical axis differs across panels.

Figure 8. Rejection frequencies for Fuller LR tests as functions of ρ for $q=8$ and $n=400$. (a) $a=2$; (b) $a=4$; (c) $a=8$; (d) $a=16$. Note that the scale of the vertical axis differs across panels.

In Figure 8, we see that the Fuller LR test never underrejects as much as the LR test, and it actually overrejects quite severely when ρ is large and $a=2$. However, that is the only case in which it overrejects much. This is the only test for which its own bootstrap DGP, namely, F(1)-ER, is arguably the best one to use. Except when the asymptotic test already works perfectly, using that bootstrap method almost always improves the performance of the test. The bottom two panels of Figure 8 look very similar to the corresponding panels of Figure 7, except that the bootstrapped Fuller test tends to underreject just a bit. It is evident that, as a increases, the LR test and its Fuller variant become almost indistinguishable.

Figure 6, Figure 7 and Figure 8 provide no clear ranking of tests and bootstrap methods. There seems to be a preference for the LR and Fuller LR tests, and for the LIML-ER and F(1)-ER bootstrap DGPs. In no case does any combination of those tests and those bootstrap DGPs overreject anything like as severely as the Sargan test bootstrapped using IV-R or IV-ER. Provided the instruments are not very weak, any of these combinations should yield reasonably accurate, but perhaps somewhat conservative, inferences in most cases.

The rather mixed performance of the bootstrap tests can be understood by using the concept of “bootstrap discrepancy,” which is a function of the nominal level of the test, say α. The bootstrap discrepancy is simply the actual rejection rate for a bootstrap test at level α minus α itself. Davidson and MacKinnon [18] shows that the bootstrap discrepancy at level α is a conditional expectation of the random variable

$$q(\alpha )\equiv R\left(Q(\alpha ,{\mu }^{*}),\mu \right)-\alpha ,$$

(41)

where $R(\alpha ,\mu )$ is the probability, under the DGP μ, that the test statistic is in the rejection region for nominal level α, and $Q(\alpha ,\mu )$ is the inverse function that satisfies the equation

$$R\left(Q(\alpha ,\mu ),\mu \right)=\alpha =Q\left(R(\alpha ,\mu ),\mu \right).$$

Thus $Q(\alpha ,\mu )$ is the true level-α critical value of the asymptotic test under μ. The random element in (41) is ${\mu }^{*}$, the bootstrap DGP. If ${\mu }^{*}=\mu$, then we see clearly that $q(\alpha )=0$, and the bootstrap discrepancy vanishes. For more detail, see Davidson and MacKinnon [18].

Figure 9. Contours of rejection frequencies for instrumental variables (IV)-R bootstrap Sargan tests.

Figure 9. Contours of rejection frequencies for instrumental variables (IV)-R bootstrap Sargan tests.

Suppose now that the true DGP ${\mu }_0$ is near the singularity. The bootstrap DGP can reasonably be expected also to be near the singularity, but most realizations are likely to be farther away from the singularity than ${\mu }_0$ itself. If ${\mu }_0$ were actually at the singularity, then any bootstrap DGP would necessarily be farther away. If the statistic used is S, then we see from Figure 4 that rejection frequencies fall as the DGP moves away from the singularity in most, but not all, directions. Thus, for most such bootstrap DGPs, $Q(\alpha ,{\mu }^{*})$ is smaller than $Q(\alpha ,{\mu }_0)$ for any α, and so the probability mass $R\left(Q(\alpha ,{\mu }^{*}),{\mu }_0\right)$ in the distribution generated by ${\mu }_0$ is greater than α. This means that $q(\alpha )$ is positive, and so the bootstrap test overrejects. However, if the statistic used is LR, the reverse is the case, as we see from Figure 5, and the bootstrap test underrejects. This is just what we see in Figure 6, Figure 7 and Figure 8.

Figure 9 and Figure 10 are contour plots similar to Figure 4 and Figure 5, but they are for bootstrap rather than asymptotic tests. The IV-R parametric bootstrap is used for the Sargan test in Figure 9, and the LIML-ER parametric bootstrap is used for the LR test in Figure 10. In both cases, there are 100,000 replications, and $B=199$. Figure 9 looks remarkably like Figure 4, with low rejection frequencies for extremely small values of a, then a ridge where rejection frequencies are very high for slightly larger values of a. The ridge is not quite as high as the one in Figure 4, and the rejection frequencies diminish more rapidly as a increases.

Figure 10. Contours of rejection frequencies for limited information maximum likelihood (LIML)-ER bootstrap LR tests.

Figure 10. Contours of rejection frequencies for limited information maximum likelihood (LIML)-ER bootstrap LR tests.

Similarly, Figure 10 looks like Figure 5, but the severe underrejection in the far left of the figure occurs over an even smaller region, and there is an area of modest overrejection nearby. Both these size distortions can be explained by Figure 5. When a is extremely small, the estimate used by the bootstrap DGP tends on average to be larger, so the bootstrap critical values tend, on average, to be overestimates. This leads to underrejection. However, there is a region where a is not quite so small in which the bootstrap DGP uses estimates of a that are sometimes too small and sometimes too large. The former causes overrejection, the latter underrejection. Because of the curvature of the rejection probability function, the net effect is modest overrejection; see Davidson and MacKinnon [17]. This is actually the case for most of the parameter values shown in the figure, but the rejection frequencies are generally not much greater than 0.05.

8. Power Considerations

Overidentification tests are performed in order to check whether some of the assumptions for the two-equation model (1) and (2) to be correctly specified are valid. Those assumptions are not valid if the DGP for Equation (1) is actually

$${\mathit{y}}_1=\mathit{Z}{\gamma }_1+{\mathit{W}}_1\mathit{\delta }+\beta {\mathit{y}}_2+{\mathit{u}}_1,$$

(42)

where the columns of the matrix ${\mathit{W}}_1$ are in the span of the columns of the matrix $\mathit{W}$ and are linearly independent of those of $\mathit{Z}$. As in Section 3, we can eliminate $\mathit{Z}$ from the model, replacing all other variables and the disturbances by their projections onto the orthogonal complement of the span of the columns of $\mathit{Z}$. The simpler model of Equation (13) and Equation (16) becomes

$$\begin{array}{cc}\hfill {\mathit{y}}_1& =\beta {\mathit{y}}_2+\delta {\mathit{w}}_p+{\mathit{v}}_1\hfill \end{array}$$

(43)

$$\begin{array}{cc}\hfill {\mathit{y}}_2& =a{\mathit{w}}_1+{\mathit{u}}_2=a{\mathit{w}}_1+\rho {\mathit{v}}_1+r{\mathit{v}}_2.\hfill \end{array}$$

(44)

The vector $\mathit{W}\mathit{\pi }$ is now written as $a{\mathit{w}}_1$ instead of $a\mathit{w}$, and the vector ${\mathit{W}}_1\mathit{\delta }$ is written as $\delta {\mathit{w}}_p$. As before, we make the normalizations that $\parallel {\mathit{w}}_1{\parallel }^2=1$ and $a^2=\mathit{\pi }{}^{\top }\mathit{W}{}^{\top }\mathit{W}\mathit{\pi }$. In addition, we normalize so that $\parallel {\mathit{w}}_p{\parallel }^2=1$ and ${\delta }^2=\mathit{\delta }{}^{\top }{\mathit{W}}_1{}^{\top }{\mathit{W}}_1\mathit{\delta }$.

The disturbance vector ${\mathit{u}}_1$ of (13) is written as ${\mathit{v}}_1$. In the rightmost expression of Equations (44), the vector ${\mathit{u}}_2$ has been replaced by $\rho {\mathit{v}}_1+r{\mathit{v}}_2$, where $r\equiv {(1-{\rho }^2)}^{1/2}$. Since we are assuming, as in Section 3, that the disturbance vectors ${\mathit{u}}_1$ and ${\mathit{u}}_2$ are normally distributed, the vectors ${\mathit{v}}_1$ and ${\mathit{v}}_2$ are independent $\mathbf{N}(\mathbf{0},\mathbf{I})$. Further, since ${\mathit{w}}_1$ and ${\mathit{w}}_p$ are not in general orthogonal, we write

$${\mathit{w}}_p=\theta {\mathit{w}}_1+t{\mathit{w}}_2,$$

where ${\mathit{w}}_2{}^{\top }{\mathit{w}}_1=0$, and $t={(1-{\theta }^2)}^{1/2}$.

The Basmann statistic $S^{\prime }$ is still given by Equation (18), which is simply an algebraic consequence of the definition (10). Since the DGP for ${\mathit{y}}_2$ is unchanged, the quantities denoted in (10) by $P_{22}$ and $M_{22}$ are the same under the alternative as under the null. Since the DGP for ${\mathit{M}}_{\mathit{W}}{\mathit{y}}_1$ is also the same under the null and the alternative, so are $M_{11}$ and $M_{12}$. Thus only $P_{11}$ and $P_{12}$ differ from the expressions for them in Equations (20). It is easy to check that neither the numerator nor the denominator of $S^{\prime }$ in (18) depends on β under the alternative, and so in our computations we set $\beta =0$ without loss of generality.

In order to analyze the asymptotic power of the Sargan test in Basmann form, we seek to express its limiting asymptotic distribution as a chi-squared variable that is non-central under the alternative. As usual, in order for the non-centrality parameter (NCP) to have a finite limit, we invoke a Pitman drift. With our normalization of ${\mathit{w}}_p$, this just means that δ is constant as the sample size n tends to infinity. Again, we cannot expect to find a limiting chi-squared distribution with weak-instrument asymptotics, and so our asymptotic construction supposes that $a\to \infty$ as $n\to \infty$.

Under the null and the alternative, the denominator of (18), divided by $(n-l)P_{22}$, is simply an estimate of the variance of ${\mathit{v}}_1$. For the purposes of the asymptotic analysis of the simpler model, it can therefore be replaced by 1. The quantity of which the limiting distribution is expected to be chi-squared is therefore $P_{11}-P_{12}^2/P_{22}$. Recall that this is just the numerator of both the S and $S^{\prime }$ statistics.

With $\beta =0$, we compute as follows:

$$\begin{array}{cc}\hfill P_{11}& ={\mathit{y}}_1{}^{\top }{\mathit{P}}_{\mathit{W}}{\mathit{y}}_1={\delta }^2+2\delta \theta x_1+2\delta t{z}_1+{\mathit{v}}_1{}^{\top }{\mathit{P}}_{\mathit{W}}{\mathit{v}}_1,\hfill \\ \hfill P_{12}& ={\mathit{y}}_1{}^{\top }{\mathit{P}}_{\mathit{W}}{\mathit{y}}_2=a(x_1+\delta \theta )+O_p(1),\text{and}\hfill \\ \hfill P_{22}& ={\mathit{y}}_2{}^{\top }{\mathit{P}}_{\mathit{W}}{\mathit{y}}_2=a^2+2a(\rho x_1+r{x}_2)+O_p(1),\hfill \end{array}$$

where the symbol $O_p(1)$ means of order unity as $a\to \infty$. As before, we let $x_i={\mathit{w}}_1{}^{\top }{\mathit{v}}_i$ for $i=1,2$, and we also let $z_i={\mathit{w}}_2{}^{\top }{\mathit{v}}_i$. Thus the limit as $a\to \infty$ of $P_{11}-P_{12}^2/P_{22}$ is

$${\delta }^2+2\delta \theta x_1+2\delta t{z}_1+{\mathit{v}}_1{}^{\top }{\mathit{P}}_{\mathit{W}}{\mathit{v}}_1-{(x_1+\delta \theta )}^2={\mathit{v}}_1{}^{\top }{\mathit{P}}_{\mathit{W}}{\mathit{v}}_1-x_1^2+{\delta }^2{t}^2+2\delta t{z}_1.$$

(45)

In Equation (17), we introduced the quantity $Q_{11}$, equal to ${\mathit{v}}_1{}^{\top }{\mathit{P}}_{\mathit{W}}{\mathit{v}}_1$ and distributed as ${\chi }_l^2$. It was expressed as the sum of three mutually independent random variables, $x_1^2$, $z_P^2$, and $t_{11}^P$. Now we separate out both the terms $x_1^2$ and $z_1^2$ to obtain

$$Q_{11}=x_1^2+z_1^2+z_P^2+t_{11}^{P_0},$$

(46)

where all four random variables above are independent, with $x_1$, $z_1$, and $z_P$ standard normal, and $t_{11}^{P_0}$ distributed as ${\chi }_{l-3}^2$. The random variable $t_{11}^{P_0}$ is not to be confused with $t_{11}^P$ in Equations (17), which is distributed as ${\chi }_{l-2}^2$. It is legitimate to write $Q_{11}$ in this way because it can be constructed as the sum of the squares of the l independent N(0,1) variables ${\mathit{w}}_j{}^{\top }{\mathit{v}}_1$, where the ${\mathit{w}}_j$ form an arbitrary orthonormal basis of the span of the columns of $\mathit{W}$.

Using (46), the right-hand side of Equation (45) can be written as

$$z_P^2+t_{11}^{P_0}+{(z_1+\delta t)}^2.$$

This is the sum of three independent random variables. The first is ${\chi }_1^2$, the second is ${\chi }_{l-3}^2$, and the last is noncentral ${\chi }_1^2({\delta }^2{t}^2)$. It follows that, when $a^2$ and the sample size both tend to infinity, which implies that the instruments are not weak, the numerator of the test statistic is distributed as ${\chi }_{l-1}^2({\delta }^2{t}^2)$. Note that, if $\theta =1$, so that ${\mathit{w}}_p={\mathit{w}}_1$, the NCP ${\delta }^2{t}^2$ vanishes.

For the general model (1) and (2), with DGP given by Equation (42), it can be shown that the NCP is

$$\frac{1}{{\sigma }_1^2}\mathit{\delta }{}^{\top }{\mathit{W}}_1{}^{\top }{\mathit{M}}_{\mathit{Z}}{\mathit{W}}_1\mathit{\delta }-\frac{1}{{\sigma }_1^2}\frac{{(\mathit{\pi }{}^{\top }\mathit{W}{}^{\top }{\mathit{M}}_{\mathit{Z}}{\mathit{W}}_1\mathit{\delta })}^2}{\mathit{\pi }{}^{\top }\mathit{W}{}^{\top }{\mathit{M}}_{\mathit{Z}}\mathit{W}\mathit{\pi }}.$$

(47)

For the simpler model given by Equations (43) and (44), the first term here collapses to ${\delta }^2$ and the second term, which arises because β has to be estimated, collapses to $-{\theta }^2{\delta }^2$. Therefore, expression (47) as a whole corresponds to ${\delta }^2{t}^2$ for the simpler model.

8.1. Finite-Sample Concerns

The asymptotic result that $S^{\prime }$ follows the ${\chi }_{l-1}^2({\delta }^2{t}^2)$ distribution strongly suggests that S, LR, and LR${}^{\prime }$ must do so as well, because all these statistics are asymptotically equivalent. In fact, a more tedious calculation than that in Equation (45) and Equation (46) shows that the limiting distribution of ${\text{LR}}^{\prime }$ as both n and a tend to infinity is the same as for $S^{\prime }$, namely ${\chi }_{l-1}^2({\delta }^2{t}^2)$. Because these results are only asymptotic, however, it is necessary to resort to simulation to investigate behavior under the alternative in finite samples.

Under the null, we were able in Section 3 to express all six quantities, the $P_{ij}$ and the $M_{ij}$, for $i,j=1,2$, in terms of eight independent random variables. Under the alternative, we require ten of these variables. For the $M_{ij}$, there is no need to change the expressions for them in (20), where we use the three variables $t_{11}^M$, $t_{22}^M$, and $z_M$, distributed respectively as ${\chi }_{n-l}^2$, ${\chi }_{n-l-1}^2$, and N(0,1). These represent the projections of ${\mathit{v}}_1$ and ${\mathit{v}}_2$ onto the orthogonal complement of the span of the instruments. For the $P_{ij}$, however, we decompose as follows:

$$\begin{array}{cc}\hfill Q_{11}& =x_1^2+z_1^2+z_P^2+t_{11}^{P_0},\hfill \\ \hfill Q_{12}& =x_1{x}_2+z_1{z}_2+z_P\sqrt{t_{22}^{P_0}},\text{and}\hfill \\ \hfill Q_{22}& =x_2^2+z_2^2+t_{22}^{P_0}.\hfill \end{array}$$

Here $x_i$, $z_i$, $i=1,2$, and $z_P$ are standard normal, $t_{11}^P$ is ${\chi }_{l-3}^2$, and $t_{22}^P$ is ${\chi }_{l-2}^2$, all seven variables being mutually independent. We can simulate both $S^{\prime }$ and LR${}^{\prime }$ very cheaply, by drawing ten random variables, independently of either the sample size n or the degree of overidentification $l-1$, because all the statistics are deterministic functions of the $P_{ij}$ and the $M_{ij}$, and, of course, n and l. The relations in (20) hold except those for $P_{11}$ and $P_{12}$. These are replaced by

$$P_{11}=Q_{11}+{\delta }^2+2\delta \theta x_1+2\delta t{z}_1,$$

and

$$P_{12}=a{x}_1+\rho Q_{11}+r{Q}_{12}+\delta (a\theta +\rho \theta x_1+\rho t{z}_1+r\theta x_2+rt{z}_2).$$

These equations differ from the corresponding ones in (20) only by terms proportional to a positive power of δ.

8.2. Simulation Evidence

Since we have seen that the ${\text{LR}}^{\prime }$ test often has much better finite-sample properties than the $S^{\prime }$ test, even when both are bootstrapped, it is important to see whether the superior performance of ${\text{LR}}^{\prime }$ comes at the expense of power. In this section, we employ simulation methods to do so.

Given the considerable size distortion of the asymptotic tests for most of that part of the parameter space considered in Section 7, we limit attention to parametric bootstrap tests. In this, we follow Horowitz and Savin [19], which argues that the best way to proceed, as long as the rejection probability of a test is far removed from its nominal level, is to consider a bootstrap test. But that proposition is based on the assumption that the bootstrap discrepancy is small enough to be ignored, which is not the case for the overidentification tests we have considered in the neighborhood of the singularity. Because of that, and because it is unreasonable to expect that there is much in the way of usable power near the singularity, it is primarily of interest to investigate power for situations in which the instruments are not too weak.

As before, all the simulation results are presented graphically. These results are based on 200,000 replications with 399 bootstrap repetitions. The same random variables are used for every set of parameter values. These experiments would have been extremely computationally demanding without the theoretical results of Section 6 and the first part of this section, which allow us to calculate everything very cheaply after we have generated and stored $200,000\times 10$ plus $200,000\times 399\times 8$ random variables. The first set of random variables is used to calculate the actual test statistics and the estimates of a and ρ, and the second set is used to calculate the bootstrap statistics.

Figure 11. Power of bootstrap tests as functions of δ with $n=400$, $q=8$, $\rho =0.5$, $\theta =0$, and $t=1$. (a) $a=2$; (b) $a=4$; (c) $a=8$; (d) $a=6$.

Figure 11. Power of bootstrap tests as functions of δ with $n=400$, $q=8$, $\rho =0.5$, $\theta =0$, and $t=1$. (a) $a=2$; (b) $a=4$; (c) $a=8$; (d) $a=6$.

We report results only for $S^{\prime }$ bootstrapped using the IV parameter estimates and for LR${}^{\prime }$ bootstrapped using the LIML estimates. Recall from Section 2 that the former results apply to S as well as $S^{\prime }$, and the latter apply to LR as well as LR${}^{\prime }$, because the test statistics in each pair are monotonically related.

Figure 11 shows power functions for $q=8$, $\rho =0.5$, and four values of a. When $a=2$, LR${}^{\prime }$ rejects much less frequently than $S^{\prime }$, both under the null and under the alternative. Both power functions level out as δ becomes large, and it appears that neither test rejects with probability one as $\delta \to \infty$. As a increases, the two power functions converge, and both tests do seem to reject with probability one for large δ.

The top two panels of Figure 12 are comparable to the top two panels of Figure 11, but with $q=2$. When $a=2$, $S^{\prime }$ now rejects less often that it did before, but LR${}^{\prime }$ rejects more often. When $a=4$, LR${}^{\prime }$ rejects very much more often than it did before, and the two power functions are quite close. We also obtained results for $a=6$, $a=8$, and $a=16$, which are not shown. For $a=6$, the power functions for $S^{\prime }$ and LR${}^{\prime }$ are extremely similar, and for $a\ge 8$ they are visually indistinguishable.

Figure 12. Power of bootstrap tests as functions of δ with $n=400$, $q=2$, $\theta =0$, and $t=1$. (a) $a=2$, $\rho =0.5$; (b) $a=4$, $\rho =0.5$; (c) $a=4$, $\rho =0.1$; (d) $a=4$, $\rho =0.9$.

Figure 12. Power of bootstrap tests as functions of δ with $n=400$, $q=2$, $\theta =0$, and $t=1$. (a) $a=2$, $\rho =0.5$; (b) $a=4$, $\rho =0.5$; (c) $a=4$, $\rho =0.1$; (d) $a=4$, $\rho =0.9$.

The bottom two panels of Figure 12 are comparable to the top right panel, except that $\rho =0.1$ or $\rho =0.9$ instead of $\rho =0.5$. It is evident that the shapes of the power functions depend on ρ, but for most values of δ the dependence is moderate. This justifies our use of $\rho =0.5$ in most of the experiments. Using other values of ρ would not change the main results.

When one power function is always above another, as is the case in all the panels of Figure 11 and Figure 12, it is difficult to conclude that one test is genuinely more powerful than the other. Perhaps greater power is just an artifact of greater rejection frequencies whether or not the null hypothesis is true.

One way to compare such tests is to graph rejection frequencies under the alternative against rejection frequencies under the null. Each point on such a “size-power curve” corresponds to some nominal level for the bootstrap test, with levels running from 0 to 1. The abscissa is the rejection frequency when the DGP satisfies the null, the ordinate the rejection frequency when the DGP belongs to the alternative. For a level of 0, the test never rejects, since bootstrap P values cannot be negative. If the level is 1, the test always rejects. As the nominal level increases from 0 to 1, we expect power (on the vertical axis) to increase more rapidly than the rejection frequency under the null (on the horizontal axis). See Davidson and MacKinnon [20].

The top two panels of Figure 13 show size-power curves for $q=2$, $a=4$, and four values of δ. Perhaps surprisingly, the curves for LR${}^{\prime }$ in the left-hand panel look remarkably similar to the ones for $S^{\prime }$ in the right-hand panel. The apparently greater power of $S^{\prime }$, which is evident in the top right panel of Figure 12, seems to be almost entirely accounted for by its greater tendency to reject under the null.

The bottom two panels of Figure 13 show size-power curves for $q=2$, $\delta =4$, and four values of a. It is clear that power increases with a, but at a decreasing rate. As $a\to \infty$, the curves converge to the one given by asymptotic theory, where the distribution under the null is central ${\chi }_{l-1}^2$ and the one under the alternative is noncentral ${\chi }_{l-1}^2({\delta }^2{t}^2)$. This curve is graphed in the figure and labelled $a=\infty$.

The asymptotic result that the test statistics follow the ${\chi }_{l-1}^2({\delta }^2{t}^2)$ distribution suggests that only the product $\delta t=\delta {(1-{\theta }^2)}^{1/2}$ influences power, and that, in particular, there should be no power beyond the level of the test when $\theta =1$. In finite samples, things turn out to be more complicated, as can be seen from Figure 14, which plots power against θ for $\delta =4$. The top two panels show results for $a=2$ and $a=4$. The $S^{\prime }$ test has substantial power when $\theta =1$ and $a=2$, which presumably reflects its tendency to overreject severely under the null when the instruments are weak. Those panels also show, once again, that $S^{\prime }$ can reject far more often than LR${}^{\prime }$ when the instruments are weak. This is much less evident in the bottom two panels, which show results for larger values of a, namely, 6 and 8.

One surprising feature of Figure 14 is that, in all cases, power initially increases as θ increases from 0, even though $\delta {(1-{\theta }^2)}^{1/2}$ declines. This is true even for quite large values of a, such as $a=16$, although, of course, it is not true for extremely large values.

Figure 13. Size-power curves for $q=2$, $n=400$, $\rho =0.5$, and $\theta =0$. (a) LR${}^{\prime }$: $a=4$ and several values of δ; (b) $S^{\prime }$: $a=4$ and several values of δ; (c) LR${}^{\prime }$: $\delta =2$ and several values of a; (d) $S^{\prime }$: $\delta =2$ and several values of a.

Figure 13. Size-power curves for $q=2$, $n=400$, $\rho =0.5$, and $\theta =0$. (a) LR${}^{\prime }$: $a=4$ and several values of δ; (b) $S^{\prime }$: $a=4$ and several values of δ; (c) LR${}^{\prime }$: $\delta =2$ and several values of a; (d) $S^{\prime }$: $\delta =2$ and several values of a.

Figure 14. Power as a function of θ for $n=400$, $\rho =0.5$, and $\delta =4$. (a) $q=2$; (b) $q=8$; (c) $q=2$; (d) $q=8$

Figure 14. Power as a function of θ for $n=400$, $\rho =0.5$, and $\delta =4$. (a) $q=2$; (b) $q=8$; (c) $q=2$; (d) $q=8$

9. Relaxing the IID Assumption

The resampling bootstraps that we looked at in Section 7 do not implicitly make the assumption that the disturbances are normal. They do, however, assume that the disturbances are pairwise IID. If instead the disturbances are heteroskedastic, then the covariance matrix of their bivariate distribution may be different for each observation. In that case, all the test statistics we have studied have distributions that depend on the pattern of heteroskedasticity, and so they are no longer approximately pivotal for the model (1) and (2) under either weak-instrument or strong-instrument asymptotics.

Andrews, Moreira, and Stock [21] proposes heteroskedasticity-robust versions of test statistics for tests about the value of β that are robust to weak instruments. Note that, although Andrews, Moreira, and Stock [6] is based on the 2004 paper and has almost the same title, it does not contain this material. However, this work cannot be applied here, because, as we have seen, the overidentification tests are not robust to weak instruments.

The role of the denominators of the statistics S, $S^{\prime }$, and LR${}^{\prime }$ is simply to provide non-robust estimates of the scale of the numerators. In order to make those statistics robust to heteroskedasticity, we have to provide robust measures instead. The numerators of all three statistics can be written as

$${\widehat{\mathit{u}}}_1{}^{\top }{\mathit{P}}_{\mathit{W}}{\widehat{\mathit{u}}}_1={\widehat{\mathit{u}}}_1{}^{\top }\mathit{W}{(\mathit{W}{}^{\top }\mathit{W})}^{-1}\mathit{W}{}^{\top }{\widehat{\mathit{u}}}_1,$$

(48)

where the vector ${\widehat{\mathit{u}}}_1$ denotes either ${\mathit{y}}_1-\mathit{Z}{\widehat{\gamma }}_{\text{IV}}-{\widehat{\beta }}_{\text{IV}}{\mathit{y}}_2$, in the case of S and $S^{\prime }$, or ${\mathit{y}}_1-\mathit{Z}{\widehat{\gamma }}_{\text{LIML}}-{\widehat{\beta }}_{\text{LIML}}{\mathit{y}}_2$, in the case of LR${}^{\prime }$. Expression (48) is a quadratic form in the $l$-vector $\mathit{W}{}^{\top }{\widehat{\mathit{u}}}_1$. The usual estimate of the covariance matrix of that vector is $\mathit{W}{}^{\top }\widehat{\Omega }\mathit{W}$, where $\widehat{\Omega }=\text{diag }{\widehat{u}}_{1i}^2$. Thus the heteroskedasticity-robust variant of all three test statistics is the quadratic form

$${\widehat{\mathit{u}}}_1{}^{\top }\mathit{W}{(\mathit{W}{}^{\top }\widehat{\Omega }\mathit{W})}^{-1}\mathit{W}{}^{\top }{\widehat{\mathit{u}}}_1.$$

(49)

There would be no point in using a heteroskedasticity-robust statistic along with a bootstrap DGP that imposed homoskedasticity. The natural way to avoid doing so is to use the wild bootstrap. In Davidson and MacKinnon [8], the wild bootstrap is shown to have good properties when used with tests about the value of β. The disturbances of the wild bootstrap DGP are given by

$$\left[\begin{array}{c}{u}_{1i}^{*}\\ u_{2i}^{*}\end{array}\right]=\left[\begin{array}{c}{\widehat{u}}_{1i}{\nu }_i^{*}\\ {\widehat{u}}_{2i}{\nu }_i^{*}\end{array}\right],$$

(50)

where ${\nu }_i^{*}$ is an auxiliary random variable with expectation 0 and variance 1. The easiest choice for the distribution of the ${\nu }_i^{*}$ is the Rademacher distribution, which sets ${\nu }_i^{*}$ to $+1$ or $-1$, each with probability one half. This is also probably the best choice in most cases; see Davidson and Flachaire [22].

The IID assumption can, of course, be relaxed in other ways. In particular, it would be easy to modify the test statistic (49) to allow for clustered data by replacing the middle matrix with one that resembles the middle matrix for the usual cluster robust covariance matrix. We could then use a variant of the wild cluster bootstrap of Cameron, Gelbach, and Miller [23] that allows for simultaneity. The Rademacher random variable associated with each cluster, the analog of ${\nu }_i^{*}$ in Equation (50), would then multiply the residuals for all observations within that cluster for both equations.

10. Concluding Remarks

We have shown that the well-known Sargan test for overidentification in a linear simultaneous-equations model estimated by instrumental variables often overrejects severely when the instruments are weak. In the same circumstances, the likelihood ratio test often underrejects severely. In Section 3 and Section 4, we provided a finite-sample analysis that explains these facts. The distributions of the test statistics we consider have a singularity when the concentration parameter vanishes and the absolute value of the correlation between the disturbances of the structural and reduced-form equations tends to one. Thus it can be risky to use asymptotic tests in this situation. In addition, we proposed a new test based on Fuller’s modified LIML estimator, which often outperforms the ordinary LR test.

We have also proposed four bootstrap methods which can be applied to all three of these tests. Although bootstrapping does not help much when the instruments are extremely weak, especially when the disturbances of the two equations are highly correlated, it does help substantially when the instruments are only moderately weak. In particular, using a bootstrap DGP based on Fuller’s estimator generally leads to much more accurate inferences than simply using asymptotic theory in this case.

There is a cost in terms of power to using a bootstrap test based on any version of the likelihood ratio statistic relative to a test based on the conventional Sargan or Basmann statistics. This cost generally seems to be very modest, except when the instruments are very weak.