Formulas and Mathematical Detail¶

Notation¶

Interest is in recovering the parameter vector from the model

\[\begin{aligned} y_{i} & =\beta^{\prime}x_{i}+\epsilon_{i}\end{aligned}\]

where \(\beta\) is \(k\) by 1, \(x_{i}\) is a \(k\) by 1 vector of observable variables and \(\epsilon_{i}\) is a scalar error. \(x_{i}\) can be separated in two types of variables. The \(k_{1}\) by 1 set of variables \(x_{1i}\) are exogenous regressors in the sense that \(E\left[x_{1i}\epsilon_{i}\right]=0\). The \(k_{2}\) by 1 variables \(x_{2i}\) are endogenous. A set of \(p\) instruments is available that satisfy the requirements for validity where \(p\geq k_{2}\). The extended model can be written as

\[\begin{split}\begin{aligned} y_{i} & =\underset{\textrm{exogenous}}{\underbrace{\beta_{1}^{\prime}x_{1i}}}+\underset{\textrm{endogenous}}{\underbrace{\beta_{2}^{\prime}x_{2i}}}+\epsilon_{i}\\ x_{2i} & =\underset{\textrm{exogenous}}{\underbrace{\gamma_{1}^{\prime}z_{1i}}}+\underset{\textrm{instruments}}{\underbrace{\gamma_{2}^{\prime}z_{2i}}}+u_{i}\end{aligned}\end{split}\]

The model can be expressed compactly

\[\begin{split}\begin{aligned} Y & =X_{1}\beta_{1}+X_{2}\beta_{2}+\epsilon=X\beta+\epsilon\\ X_{2} & =Z_{1}\gamma_{1}+Z_{2}\gamma_{2}+u=Z\gamma+u\end{aligned}\end{split}\]

The vector of instruments \(z_{i}\) is \(p\) by 1. There are \(n\) observations for all variables. \(k_{c}=1\) if the model contains a constant (either explicit or implicit, i.e., including all categories of a dummy variable). The constant, if included, is in \(x_{1i}\). \(X\) is the \(n\) by \(k\) matrix of regressors where row \(i\) of \(X\) is \(x_{i}^{\prime}\). \(X\) can be partitioned into \(\left[X_{1}\;X_{2}\right]\). \(Z\) is the \(n\) by \(p\) matrix of instruments. The vector \(y\) is \(n\) by 1. Projection matrices for \(X\) is defined \(P_{X}=X\left(X^{\prime}X\right)^{-1}X^{\prime}\). The projection matrix for \(Z\) is similarly defined only using \(Z\). The annihilator matrix for \(X\) is \(M_{X}=I-P_{X}\).

Observation weights, \(w_{i}\), are scaled to have a mean of 1. When a model uses weights, estimation is performed on the weighted variables \(\sqrt{w_{i}}y_{i}\), \(\sqrt{w_{i}}x_{i}\) and \(\sqrt{w_{i}}z_{i}\), which take the place of \(y_{i}\), \(x_{i}\) and \(z_{i}\) in the formulas below. This includes the covariance estimators, the specification tests and the first-stage statistics. The parameters, covariance and test statistics of a weighted model are therefore those of the unweighted model that is estimated using the weighted variables.

Parameter Estimation¶

Two-stage Least Squares (2SLS)¶

The 2SLS estimator is

\[\begin{aligned} \hat{\beta}_{2SLS} & =\left(X^{\prime}P_{Z}X\right)^{-1}\left(X^{\prime}P_{Z}y\right)\end{aligned}\]

Limited Information Maximum Likelihood and k-class Estimators¶

The LIML or other k-class estimator is

\[\begin{aligned} \hat{\beta}_{\kappa} & =\left(X^{\prime}\left(I-\kappa M_{Z}\right)X\right)^{-1}\left(X^{\prime}\left(I-\kappa M_{Z}\right)y\right)\end{aligned}\]

where \(\kappa\) is the parameter of the class. When \(\kappa=1\) the 2SLS estimator is recovered. When \(\kappa=0\), the OLS estimator is recovered. The LIML estimator is recovered when \(\kappa\) set to

\[\hat{\kappa}=\min\mathrm{eig\left[\left(W^{\prime}M_{Z}W\right)^{-\frac{1}{2}}\left(W^{\prime}M_{X_{1}}W\right)\left(W^{\prime}M_{Z}W\right)^{-\frac{1}{2}}\right]}\]

where \(W=\left[y\:X_{2}\right]\) and \(\mathrm{eig}\) returns the eigenvalues.

Generalized Method of Moments (GMM)¶

The GMM estimator is defined as

\[\begin{aligned} \hat{\beta}_{GMM} & =\left(X^{\prime}ZWZ^{\prime}X\right)^{-1}\left(X^{\prime}ZWZ^{\prime}y\right)\end{aligned}\]

where \(W\) is a positive definite weighting matrix.

Continuously Updating Generalized Method of Moments (GMM-CUE)¶

The continuously updating GMM estimator solves the minimization problem

\[\min_{\beta}n\bar{g}\left(\beta\right)^{\prime}W\left(\beta\right)^{-1}\bar{g}\left(\beta\right)\]

where \(\bar{g}\left(\beta\right)=n^{-1}\sum_{i=1}^{n}z_{i}\hat{\epsilon}_{i}\) and \(\hat{\epsilon}_{i}=y_{i}-x_{i}\beta\). Unlike standard GMM, the weight matrix, \(W\) directly depends on the model parameters \(\beta\) and so it is not possible to use a closed form estimator.

Jackknife Instrumental Variables (JIVE)¶

The 2SLS fitted value of a regressor, \(\hat{x}_{i}\), the \(i\)th row of \(\hat{X}=P_{Z}X\), is computed from a first stage that uses observation \(i\). This links \(\hat{x}_{i}\) to the error \(\epsilon_{i}\), which is the source of the bias of 2SLS towards OLS when there are many instruments. The jackknife IV estimator (JIVE1) of Angrist, Imbens and Krueger (1999) replaces the fitted values with leave-one-out predictions. Row \(i\) of \(\tilde{X}\) is the prediction of \(x_{i}\) from the first-stage regression of \(X\) on \(Z\) that is estimated without observation \(i\),

\[\tilde{x}_{i}^{\prime}=z_{i}^{\prime}\left(Z_{(-i)}^{\prime}Z_{(-i)}\right)^{-1}Z_{(-i)}^{\prime}X_{(-i)}\]

where the subscript \((-i)\) indicates that row \(i\) has been removed. The first stage does not have to be estimated \(n\) times, since

\[\tilde{x}_{i}=\frac{\hat{x}_{i}-h_{i}x_{i}}{1-h_{i}},\qquad h_{i}=z_{i}^{\prime}\left(Z^{\prime}Z\right)^{-1}z_{i}\]

where \(h_{i}\) is the leverage of observation \(i\), the \(i\)th diagonal element of \(P_{Z}\). The exogenous regressors are in \(Z\), so that they are their own predictions, \(\tilde{x}_{1i}=x_{1i}\), and only the endogenous regressors are replaced. The JIVE estimator is the IV estimator that uses \(\tilde{X}\) as the instruments for \(X\),

\[\hat{\beta}_{JIVE}=\left(\tilde{X}^{\prime}X\right)^{-1}\tilde{X}^{\prime}y\]

which replaces \(\hat{X}\) with \(\tilde{X}\) in \(\hat{\beta}_{2SLS}=\left(\hat{X}^{\prime}X\right)^{-1}\hat{X}^{\prime}y\). The estimator cannot be computed if \(h_{i}=1\) for an observation. This happens if the observation is the only one with a non-zero value of one of the instruments, for example a dummy variable that is 1 for a single observation.

Basic Statistics¶

Let \(\hat{\epsilon}=y-X\hat{\beta}\). The residual sum of squares (RSS) is \(\hat{\epsilon}^{\prime}\hat{\epsilon}\), the model sum of squares (MSS) is \(\hat{\beta}^{\prime}X^{\prime}X\hat{\beta}\) and the total sum of squares (TSS) is \(\left(y-k_{c}\bar{y}\right)^{\prime}\left(y-k_{c}\bar{y}\right)^{\prime}\)where \(\bar{y}\) is the sample average of \(y\). The model \(R^{2}\)is defined

\[\begin{aligned} R^{2} & =1-\frac{\hat{\epsilon}^{\prime}\hat{\epsilon}}{\left(y-k_{c}\bar{y}\right)^{\prime}\left(y-k_{c}\bar{y}\right)^{\prime}}=1-\frac{RSS}{TSS}\end{aligned}\]

and the adjusted \(R^{2}\) is defined

\[\begin{aligned} \bar{R}^{2} & =1-\left(1-R^{2}\right)\frac{N-k_{c}}{N-k}.\end{aligned}\]

The \(R^{2}\) is NaN, with no error or warning, if the total sum of squares is zero, up to rounding error, since it is not defined then, e.g., if the dependent variable does not vary in a model with a constant. The adjusted \(R^{2}\) is NaN as well, and it is NaN if \(N=k\).

The residual variance is \(s^{2}=n^{-1}\hat{\epsilon}^{\prime}\hat{\epsilon}\) unless the debiased flag is used, in which case a small sample adjusted version is estimated \(s^{2}=\left(n-k\right)^{-1}\hat{\epsilon}^{\prime}\hat{\epsilon}\). The model degree of freedom is \(k\) and the residual degree of freedom is \(n-k\).

The model F-statistic is defined

\[\begin{aligned} F & =\hat{\beta}_{-}^{\prime}\hat{V}_{-}^{-1}\dot{\hat{\beta}_{-}}\end{aligned}\]

where the notation \(\hat{\beta}_{-}\) indicates that the constant is excluded from \(\hat{\beta}\) and \(\hat{V}_{-}\) indicates that the row and column corresponding to a constant are excluded. [1] The test statistic is distributed as \(\chi_{k-k_{c}}^{2}.\) If the debiased flag is set, then the test statistic \(F\) is transformed as \(F/\left(k-k_{c}\right)\) and a \(F_{k-k_{c},n-k}\) distribution is used.

Parameter Covariance Estimation¶

Two-stage LS, LIML and k-class estimators¶

Four covariance estimators are available. The first is the standard homoskedastic covariance, defined as

\[\begin{aligned} \hat{\Sigma}=n^{-1}s^{2}\left(\frac{X^{\prime}\left(I-\kappa M_{z}\right)X}{n}\right)^{-1} & =n^{-1}s^{2}\hat{A}.\end{aligned}\]

Note that this estimator can be expressed as

\[\begin{aligned} \hat{\Sigma}=n^{-1}\hat{A}^{-1}\left\{ s^{2}\hat{A}\right\} \hat{A}^{-1} & =n^{-1}\hat{A}^{-1}\hat{B}\hat{A}^{-1}.\end{aligned}\]

All estimators take this form and only differ in how the asymptotic covariance of the scores, \(B\), is estimated. For the homoskedastic covariance estimator, \(\hat{B}=s^{2}\hat{A}.\) The score covariance in the heteroskedasticity robust covariance estimator is

\[\begin{aligned} \hat{B} & =n^{-1}\sum_{i=1}^{n}\hat{\epsilon}_{i}^{2}\hat{x}_{i}\hat{x}_{i}^{\prime}=n^{-1}\sum_{i=1}^{n}\hat{\xi}_{i}\hat{\xi}_{i}^{\prime}.\end{aligned}\]

where \(\hat{x_{i}}\) is row \(i\) of \(\hat{X}=P_{Z}X\) and \(\hat{\xi}_{i}=\hat{\epsilon}_{i}\hat{x}_{i}\).

The kernel covariance estimator is robust to both heteroskedasticity and autocorrelation and is defined as

\[\begin{split}\begin{aligned} \hat{B} & =\hat{\Gamma}_{0}+\sum_{i=1}^{n-1}K\left(\frac{i}{h}\right)\left(\hat{\Gamma}_{i}+\hat{\Gamma}_{i}^{\prime}\right)\\ \hat{\Gamma_{j}} & =n^{-1}\sum_{i=j+1}^{n}\hat{\xi}_{i-j}\hat{\xi}_{i}^{\prime}\end{aligned}\end{split}\]

where \(K\left(\frac{i}{h}\right)\) is a kernel weighting function where \(h\) is the kernel bandwidth.

The one-way clustered covariance estimator is defined as

\[\begin{aligned} n^{-1}\sum_{j=1}^{g}\left(\sum_{i\in\mathcal{G}_{j}}\hat{\xi}_{i}\right)\left(\sum_{i\in\mathcal{G}_{j}}\hat{\xi}_{i}\right)^{\prime}\end{aligned}\]

where \(\sum_{i\in\mathcal{G}_{j}}\hat{\xi}_{i}\) is the sum of the scores for all members in group \(\mathcal{G}_{j}\) and \(g\) is the number of groups.

If the debiased flag is used to perform a small-sample adjustment, all estimators except the clustered covariance are rescaled by \(\left(n-k\right)/n\). The clustered covariance is rescaled by \(\left(\left(n-k\right)\left(n-1\right)/n^{2}\right)\left(\left(g-1\right)/g\right)\). [2]

Standard Errors¶

Standard errors are defined as

\[s.e.\left(\hat{\beta}_{j}\right)=\sqrt{e_{j}^{\prime}\hat{\Sigma}e_{j}}\]

where \(e_{j}\) is a vector of 0s except in location \(j\) which is 1.

T-statistics¶

T-statistics test the null \(H_{0}:\beta_{j}=0\) against a 2-sided alternative and are defined as

\[z=\frac{\hat{\beta}_{j}}{s.e.\left(\hat{\beta}_{j}\right)}.\]

P-values¶

P-values are computes using 2-sided tests,

\[Pr\left(\left|z\right|>Z\right)=2-2\Phi\left(\left|z\right|\right)\]

If the covariance estimator was debiased, a Student’s t distribution with \(n-k\) degrees of freedom is used,

\[\begin{aligned} Pr\left(\left|z\right|>Z\right) & =2-2t_{n-k}\left(\left|z\right|\right)\end{aligned}\]

where \(t_{n-k}\left(\cdot\right)\) is the CDF of a Student’s T distribution.

Confidence Intervals¶

Confidence intervals are constructed as

\[CI_{i,1-\alpha}=\hat{\beta}_{i}\pm q_{\alpha/2}\times\hat{\sigma}_{\beta_{i}}\]

where \(q_{\alpha/2}\) is the \(\alpha/2\) quantile of a standard Normal distribution or a Student’s t. The Student’s t is used when a debiased covariance estimator is used.

Linear Hypothesis Tests¶

Linear hypothesis tests examine the validity of nulls of the form \(H_{0}:R\beta-r=0\) and are implemented using a Wald test statistic

\[W=\left(R\hat{\beta}-r\right)^{\prime}\left[R^{\prime}\hat{\Sigma}R\right]^{-1}\left(R\hat{\beta}-r\right)\sim\chi_{q}^{2}\]

where \(q\) is the \(rank\left(R\right)\) which is usually the number of rows in \(R\) . If the debiased flag is used, then \(W/q\) is reported and critical and p-values are taken from a \(F_{q,n-k}\) distribution.

GMM Covariance estimators¶

GMM covariance depends on the weighting matrix used in estimation and the assumed covariance of the scores. In most applications these are the same and so the inefficient form,

\[\hat{\Sigma}=n^{-1}\left(\left(\frac{X'Z}{n}\right)W\left(\frac{Z'X}{n}\right)\right)^{-1}\left(\left(\frac{X'Z}{n}\right)WSW\left(\frac{Z'X}{n}\right)\right)\left(\left(\frac{X'Z}{n}\right)W\left(\frac{Z'X}{n}\right)\right)^{-1}\]

will collapse to the simpler form

\[\hat{\Sigma}=n^{-1}\left(\left(\frac{X'Z}{n}\right)W\left(\frac{Z'X}{n}\right)\right)^{-1}\]

when \(W=S^{-1}\). When an unadjusted (homoskedastic) covariance is used,

\[\hat{S}=\tilde{s}^{2}n^{-1}\sum_{j=1}^{n}z_{j}^{\prime}z_{j}\]

where \(\tilde{s}^{2}=n^{-1}\sum_{i=1}^{n}\left(\epsilon_{i}-\bar{\epsilon}\right)^{2}\) subtracts the mean which may be non-zero if the model is overidentified. Like previous covariance estimators, if the debiased flag is used, \(n^{-1}\) is replaced by \(\left(n-k\right)^{-1}\). When a robust (heteroskedastic) covariance is used, the estimator of \(S\) is modified to

\[\hat{S}=n^{-1}\sum_{i=1}^{n}\hat{\epsilon}_{i}^{2}z_{i}^{\prime}z_{i}.\]

If the debiased flag is used, \(n^{-1}\) is replaced by \(\left(n-k\right)^{-1}\).

Kernel covariance estimators of \(S\) take the form

\[\begin{split}\begin{aligned} \hat{S} & =\hat{\Gamma}_{0}+\sum_{i=1}^{n-1}k\left(i/h\right)\left(\hat{\Gamma}_{i}+\hat{\Gamma}_{i}^{\prime}\right)\\ \hat{\Gamma_{j}} & =n^{-1}\sum_{i=j+1}^{n}\hat{\epsilon}_{i-j}\hat{\epsilon}_{i}z_{i-j}^{\prime}z_{i}\end{aligned}\end{split}\]

and \(k\left(\cdot\right)\) is a kernel weighting function with bandwidth \(h\). If the debiased flag is used, \(n^{-1}\) is replaced by \(\left(n-k\right)^{-1}\).

The one-way clustered covariance estimator is defined as

\[\hat{S}=n^{-1}\sum_{j=1}^{g}\left(\sum_{i\in\mathcal{G}_{j}}\hat{\epsilon}_{i}z_{i}\right)^{\prime}\left(\sum_{i\in\mathcal{G}_{j}}\hat{\epsilon}_{i}z_{i}\right)\]

where \(\sum_{i\in\mathcal{G}_{j}}\hat{\epsilon}_{i}z_{i}\) is the sum of the moment conditional for all members in group \(\mathcal{G}_{j}\) and \(g\) is the number of groups. If the debiased flag is used, the \(n^{-1}\) term is replaced by

\[\left(n-k\right)^{-1}\frac{n-1}{n}\frac{g}{g-1}.\]

GMM Weight Estimators¶

The GMM optimal weight estimators are identical to the the estimators of \(S\) with two notable exceptions. First, they are never debiased and so always use \(n^{-1}\). Second, if the center flag is true, the demeaned moment conditions defined as \(\tilde{g}_{i}=z_{i}\hat{\epsilon}_{i}-\overline{z\epsilon}\) are used in-place of \(g_{i}\) in the robust, kernel and clustered estimators. The unadjusted estimator is always centered, and so this option has no effect.

JIVE Covariance¶

JIVE is a just identified IV estimator that solves the moment conditions \(n^{-1}\sum_{i=1}^{n}\tilde{x}_{i}\hat{\epsilon}_{i}=0\), where \(\hat{\epsilon}_{i}=y_{i}-x_{i}^{\prime}\hat{\beta}_{JIVE}\) is computed using \(X\) and not \(\tilde{X}\). Its covariance is

\[\hat{\Sigma}=n^{-1}\left(\frac{\tilde{X}^{\prime}X}{n}\right)^{-1}\hat{S}\left(\frac{X^{\prime}\tilde{X}}{n}\right)^{-1}\]

where \(\hat{S}\) is an estimator of the covariance of the moment conditions \(\tilde{x}_{i}\hat{\epsilon}_{i}\). The estimators are those of the GMM covariance, with \(\tilde{x}_{i}\) in place of \(z_{i}\). The unadjusted estimator is \(\hat{S}=s^{2}n^{-1}\tilde{X}^{\prime}\tilde{X}\) with \(s^{2}=n^{-1}\hat{\epsilon}^{\prime}\hat{\epsilon}\), which is not centered. The robust estimator is \(\hat{S}=n^{-1}\sum_{i=1}^{n}\hat{\epsilon}_{i}^{2}\tilde{x}_{i}\tilde{x}_{i}^{\prime}\). The kernel and one-way clustered estimators are defined as for the GMM covariance. The bandwidth of the kernel estimator is selected using the data if it is not provided. Two-way clustering is not available. If the debiased flag is used, \(\hat{S}\) is rescaled as in the GMM covariance estimators.

Inference with many instruments

These covariance estimators treat \(\tilde{X}\) as if it were a fixed set of instruments. Chao, Swanson, Hausman, Newey and Woutersen (2012) derive the distribution of JIVE when the number of instruments, \(K\), grows with the sample size. JIVE is consistent and asymptotically normal if \(\sqrt{K}/r_{n}\to0\), where \(r_{n}\) is the concentration parameter, which measures the strength of the instruments. The asymptotic variance is the sum of the term that is estimated above and a second term of order \(K/r_{n}\). The second term depends on the variances of the errors of the first stage and of the model and on their covariance. It can be ignored if \(K/r_{n}\to0\), but not if \(K\) is large relative to the strength of the instruments. The covariance estimators above omit it, and so they can understate the uncertainty in the latter case. In a simulation with \(n=200\), \(K=100\) instruments that each have a first-stage coefficient of 0.05, and a correlation of 0.95 between the errors of the first stage and of the model, a 95% confidence interval that used the robust covariance contained the true value in about 85% of 2,000 replications. The variance estimator that includes the second term is not implemented.

Davidson and MacKinnon (2006) compare two variants of JIVE with 2SLS and LIML in Monte Carlo experiments. They find that JIVE is always more dispersed than 2SLS, often very much so, that it is almost always inferior to LIML, and that it performs particularly badly when the instruments are weak.

Post-estimation¶

A model that uses weights is tested using the weighted variables, as described in the notation. Sargan’s, Basmann’s, Anderson-Rubin, Basmann’s F, Durbin’s and the Wu-Hausman tests assume that the errors are homoskedastic, and so they do not use the covariance estimator that was used to fit the model. Wooldridge’s score tests of exogeneity and of the overidentifying restrictions are robust to heteroskedasticity, and they also do not depend on the covariance estimator. Wooldridge’s regression test of exogeneity uses the covariance estimator that was used to fit the model, and so it has the robustness of that estimator. The J-statistic and the C-statistic of the GMM estimators use the estimator of the covariance of the moment conditions that is selected by the weight type, and not the covariance estimator that is selected in the call to fit.

2SLS and LIML Post-estimation Results¶

Sargan

Sargan’s test of over-identifying restrictions examines whether the variance of the IV residuals is similar to that of the OLS residuals. The test statistic is computed

\[s=n(1-\hat{\epsilon}^{\prime}M_{Z}\hat{\epsilon}/\hat{\epsilon}^{\prime}\hat{\epsilon})\sim\chi_{v}^{2}\]

where \(\hat{\epsilon}\) are the IV residuals and \(M_{Z}\) is the annihilator matrix using all exogenous variables.\(\nu\) is the number of overidentifying restrictions, which is the number of instruments minus the number of endogenous variables, \(p-k_{2}\).

Basmann

Basmann’s test is a small-sample corrected version of Sargan’s test of over-identifying restrictions. It has the same distribution. Let \(s\) be Sargan’s test statistic, then Basmann’s test statistic is

\[s(n-p)/(n-s)\sim\chi_{v}^{2}\]

Wooldridge score

Wooldridge’s score test of exogeneity examines the magnitude of the correlation between the OLS residuals and a an appropriately constructed set of residuals of the instruments. Define \(e=M_{X}Y\) and \(\nu=M_{X}M_{Z}X_{2}\). Then the test statistic is computed from the regression

\[1=\gamma_{1}\hat{\epsilon}_{1}\hat{v}_{1,i}+\ldots+\gamma_{p}\hat{\epsilon}_{1}\hat{v}_{p,i}+\eta_{i}\]

as \(nR^{2}\sim\chi_{k_{2}}^{2}\).

Wooldridge regression

Wooldridge’s regression test of exogeneity is closely related to the score test and is generally more powerful. It also inherits robustness to heteroskedasticity and/or autocorrelation the comes from the choice of covariance estimator used in the model. Define \(R=M_{Z}X_{2}\). The test is implemented in a regression of

\[Y=X\beta+R\gamma+\eta\]

\[\hat{\gamma}^{\prime}\hat{\Sigma}_{\gamma}^{-1}\gamma^{\prime}\sim\chi_{k_{2}}^{2}\]

where \(\hat{\Sigma}_{\gamma}\) is the block of the covariance matrix corresponding to the \(\gamma\) parameters. \(\hat{\Sigma}\) is estimated using the same covariance estimator as the model fit.

Wooldridge’s Test of Overidentifying restrictions

Wooldridge’s test is a score test examining whether the component of the instrument that is uncorrelated with both the included exogenous and the fitted exogenous is uncorrelated with the IV residuals. Define \(\tilde{Z}=M_{\left[X_{1}\:\hat{X}_{2}\right]}Z_{2,1:q}\) where \(\hat{X_{2}}\) are the fitted values from the first stage regression of the endogenous on all exogenous variables and \(Z_{2,1:q}\) contains any \(q\) columns of \(Z_{2}\), \(q=p-k_{2}\) . The test statistic is computed using a regression of 1s on the test functions \(\hat{\epsilon}_{i}\tilde{z}_{i,j}\) for \(j=1,\ldots,q\) which should have expected value 0.

\[1=\gamma_{1}\hat{\epsilon}_{i}\tilde{z}_{i,1}+\ldots+\gamma_{q}\hat{\epsilon}_{i}\tilde{z}_{i,q}\]

The test statistic is \(nR^{2}\sim\chi_{q}^{2}\) from the regression.

Anderson-Rubin

The Anderson-Rubin test of overidentification examines the magnitude of the LIML \(\hat{\kappa}\)which should be close to unity when the model is not overidentified.

\[n\ln(\hat{\kappa})\sim\chi_{q}^{2}\]

where \(q=p-k_{2}\).

Basmann’s F

Basmann’s F test of overidentification also examines the magnitude of the LIML \(\hat{\kappa}\). The test statistic is

\[\hat{\kappa}(n-p)/q\sim F_{q,n-n_{instr}}\]

where \(q=p-k_{2}\).

Durbin and Wu-Hausman

Durbin’s and the Wu-Hausman tests of exogeneity test of exogeneity is depends on the variance of the residuals when come endogenous variables are treated as exogenous against the case where they are treated as endogenous. Durbin’s test statistic is

\[\begin{split}\begin{aligned} \delta= & \hat{\epsilon}'_{e}P_{[z,w]}\hat{\epsilon}_{e}-\hat{\epsilon}'_{c}P_{z}\hat{\epsilon}_{c}\\ D= & \delta/(\hat{\epsilon}'_{e}\hat{\epsilon}_{e})/n\sim\chi_{q}^{2}\end{aligned}\end{split}\]

and the Wu-Hausman test statistic is

\[\begin{aligned} WH= & \frac{\delta/q}{(\hat{\epsilon}'_{e}\hat{\epsilon}_{e}-\delta)/v}\sim F_{q,\nu}\end{aligned}\]

where \(\hat{\epsilon}_{e}\) treats the selected set of endogenous variables as exogenous (efficient estimate) and \(\hat{\epsilon}_{c}\) is a consistent estimator if these variables are endogenous.\(P_{\left[Z\,W\right]}\) is the projection matrix containing all exogenous variables including the instrument as well as the variables being tested for endogeneity \(\left(W\right)\).\(q\) is the number of variables being tested for exogeneity and \(\nu=n-k1-k2-q\).

GMM Post-estimation Results¶

J-stat

The J-statistic tests whether the model is over-identified, and is defined

\[n\bar{g}'W^{-1}\bar{g}\sim\chi_{q}^{2}\]

where \(\bar{g}=n^{-1}\sum\hat{\epsilon}_{i}z_{i}\) and \(W\) is a consistent estimator of the variance of \(\sqrt{n}\bar{g}\). The degree of freedom is \(q=p-k_{2}\).

Centered J-stat

The estimators of the covariance of the moment conditions that are used to compute the weight matrix subtract the mean of the moment conditions only if the center flag is true, and the unadjusted estimator always does. Let \(g_{i}=z_{i}\hat{\epsilon}_{i}\), and let \(W=n^{-1}\sum_{i=1}^{n}g_{i}g_{i}^{\prime}\) and \(W_{c}=n^{-1}\sum_{i=1}^{n}\left(g_{i}-\bar{g}\right)\left(g_{i}-\bar{g}\right)^{\prime}\) be the robust estimators that do not center and that center the moment conditions, so that \(W=W_{c}+\bar{g}\bar{g}^{\prime}\). If the model is correctly specified, the moment conditions have mean zero at the true parameters and the two estimators are asymptotically equivalent. If it is not, \(\bar{g}\) does not converge to zero, and \(W\) also contains \(\bar{g}\bar{g}^{\prime}\). When both are evaluated using the same parameters, the Sherman-Morrison formula shows that the J-statistic and the centered J-statistic, \(J_{c}=n\bar{g}^{\prime}W_{c}^{-1}\bar{g}\), are related by

\[J=n\bar{g}^{\prime}W^{-1}\bar{g}=\frac{J_{c}}{1+J_{c}/n},\qquad J_{c}=\frac{J}{1-J/n}\]

The J-statistic is therefore smaller than \(n\) however large the violation of the moment conditions, while the centered J-statistic is not bounded. The relationship is exact when the weight matrix is the robust estimator that does not center, evaluated at the final parameters, as in iterated GMM that has converged. In two-step GMM the weight matrix is computed from the first-step parameters, and the relationship is approximate. The centered J-statistic has the same asymptotic distribution as the J-statistic when the moment conditions hold, and it is available as the attribute centered_j_stat of the results (Hall, 2000; Hansen and Lee, 2021). It uses the estimator of the covariance that is selected by the weight type of the model, applied to the centered moment conditions and evaluated at the final parameter estimates. The cluster sums of centered moment conditions add to zero, and so a clustered estimate that uses \(g\) clusters has rank at most \(g-1\). The statistic is not defined if the estimate is singular, which requires more clusters than moment conditions.

C-stat

The C-statistic tests the exogeneity of \(m\) of the endogenous variables by treating them as exogenous, so that they are included in the moment conditions, \(z_{i}^{e}=\left(z_{i}^{\prime},x_{t,i}^{\prime}\right)^{\prime}\), where \(x_{t,i}\) are the tested variables. A likelihood-ratio-like test statistic is computed as

\[C=J_{e}-J_{c}\sim\chi_{m}^{2}\]

where \(J_{e}\) is the J-statistic of the model that treats the tested variables as exogenous and \(J_{c}\) is the J-statistic of the model that leaves them as endogenous. Both J-statistics use the same estimate of the covariance of the moment conditions. Let \(S_{e}\) be the estimate of the covariance of \(z_{i}^{e}\hat{\epsilon}_{i}\) computed from the first-step residuals of the model that treats the variables as exogenous. It is computed using the estimator that is selected by the weight type of the model, for example robust, kernel or clustered. The weight matrix of that model is \(S_{e}^{-1}\), and \(J_{e}\) is its J-statistic. The moment conditions of the original model are a subset of the moment conditions of this model. Let \(S_{e,c}\) be the sub-matrix of \(S_{e}\) that contains the rows and columns of the original moment conditions. The original model is estimated using the weight matrix \(S_{e,c}^{-1}\), and \(J_{c}\) is its J-statistic. Note that \(S_{e,c}^{-1}\) is the inverse of the sub-matrix of \(S_{e}\), which is not the sub-matrix of \(S_{e}^{-1}\).

Using a common estimate of the covariance ensures that the test statistic is not negative, since \(g_{e}^{\prime}S_{e}^{-1}g_{e}\geq g_{c}^{\prime}S_{e,c}^{-1}g_{c}\) for any parameters, where \(g_{e}\) and \(g_{c}\) are the average moment conditions of the two models. If the original model is just identified then \(J_{c}=0\) and the C-statistic is \(J_{e}\). The statistic is not defined if \(S_{e}\) is singular, which requires more clusters than the moment conditions of the model that treats the variables as exogenous.

First-stage Estimation Analysis¶

Partial R2 and Partial F-statistic

The \(R^{2}\) is reported after orthogonalizing the instruments from included exogenous variables, so that the model estimated is

\[x_{2i}=\gamma_{0}+\tilde{z}_{2i}\gamma+\eta_{i}\]

where \(\tilde{Z}_{2}=M_{X_{1}}\tilde{Z}\). The partial \(F\)-statistic tests \(H_{0}:\gamma=0\) in this regression. It is the Wald statistic \(W=\hat{\gamma}^{\prime}\hat{V}_{\gamma}^{-1}\hat{\gamma}\), where \(\hat{V}_{\gamma}\) is the covariance of \(\hat{\gamma}\) estimated using the covariance estimator of the model. If the covariance is unadjusted, the statistic is reported as \(W/p_{2}\) with an \(F_{p_{2},n-p}\) distribution. This is the classical \(F\)-statistic only if the debiased flag is used. Otherwise, the covariance uses \(n^{-1}\) in place of \(\left(n-p\right)^{-1}\) and the statistic is \(n/\left(n-p\right)\) times larger. In all other cases, \(W\) is reported with an asymptotic \(\chi_{p_{2}}^{2}\) distribution. It is not divided by \(p_{2}\), which is how other software reports the \(F\)-statistic of a robust test.

Shea’s R2

Shea’s \(R^{2}\) is defined as the ratio of the parameter variances under OLS and 2SLS estimation standardized by the unexplained variance under each,

\[\frac{\frac{\hat{\sigma}_{OLS,\beta_{i}}^{2}}{1-R_{OLS}^{2}}}{\frac{\hat{\sigma}_{IV,\beta_{i}}^{2}}{1-R_{IV}^{2}}}=\frac{\hat{\sigma}_{OLS,\beta_{i}}^{2}}{\hat{\sigma}_{IV,\beta_{i}}^{2}}\times\frac{1-R_{IV}^{2}}{1-R_{OLS}^{2}}\]

If the estimator under 2SLS was as good as under OLS, both ratios would be 1 and Shea’s \(R^{2}=1\). On the other hand, the worse the \(IV\) fit in terms of either \(R^{2}\) or the parameter variances, the lower the value reported by Shea’s \(R^{2}\).

Kernel Weights and Bandwidth Selection¶

Kernel weights

In all formulas, \(m\) is the kernel bandwidth parameter.

Bartlett

\[w_{i}=1-\frac{i}{m+1},\,i<m\]
Parzen

\[\begin{split}\begin{aligned} z_{i} & =\frac{i}{m+1},\,i<m\\ w_{i} & =1-6z_{i}^{2}+6z_{i}^{3},z\leq0.5\\ w_{i} & =2(1-z_{i})^{3},z>0.5\end{aligned}\end{split}\]
Quadratic-Spectral

\[\begin{split}\begin{aligned} z_{i} & =\frac{6\pi i}{5m}\\ w_{0} & =1\\ w_{i} & =3(\sin(z_{i})/z_{i}-\cos(z_{i}))/z_{i}^{2},\:i\geq1\end{aligned}\end{split}\]

Optimal BW selection

If the bandwidth of a kernel estimator is not provided, it is selected using the data with the plug-in method of Newey and West (1994). The method uses a scalar series, \(g_{i}\), that is a sum of the scores or moment conditions. In the 2SLS and LIML covariances, \(g_{i}\) is the sum of \(\hat{x}_{ij}\hat{\epsilon}_{i}\) over the regressors other than the constant, where \(\hat{x}_{ij}\) are the first-stage fitted values. In the GMM weight matrix and covariance and in the JIVE covariance, the moment conditions are divided by their standard deviations before they are summed. The moment conditions are \(z_{ij}\hat{\epsilon}_{i}\), and \(\tilde{x}_{ij}\hat{\epsilon}_{i}\) in JIVE. Let

\[\hat{\sigma}_{j}=T^{-1}\sum_{i=j+1}^{T}g_{i}g_{i-j},\qquad\hat{s}^{(0)}=\hat{\sigma}_{0}+2\sum_{j=1}^{m^{*}}\hat{\sigma}_{j},\qquad\hat{s}^{(q)}=2\sum_{j=1}^{m^{*}}j^{q}\hat{\sigma}_{j}\]

where \(T\) is the number of observations. The number of autocovariances, \(m^{*}\), and the constants \(q\) and \(c_{q}\) depend on the kernel. For the Bartlett kernel, \(q=1\), \(c_{q}=1.1447\) and \(m^{*}=\lceil4\left(T/100\right)^{2/9}\rceil\). For the quadratic spectral kernel, \(q=2\), \(c_{q}=1.3221\) and \(m^{*}=\lceil4\left(T/100\right)^{2/25}\rceil\). For the Parzen kernel, \(q=2\), \(c_{q}=2.6614\) and \(m^{*}=\lceil4\left(T/100\right)^{4/25}\rceil\). The bandwidth is

\[\hat{m}=\left\lceil c_{q}\left[\left(\frac{\hat{s}^{(q)}}{\hat{s}^{(0)}}\right)^{2}\right]^{\frac{1}{2q+1}}T^{\frac{1}{2q+1}}\right\rceil\]

or \(T-1\), if that is smaller.

The method needs \(m^{*}<T-1\). Otherwise \(\hat{s}^{(0)}\) is the sum of all of the autocovariances, \(T^{-1}\left(\sum_{i}g_{i}\right)^{2}\), which is zero, up to rounding error, for a series with a mean of zero, such as the scores of a least squares regression, and the method divides by it. The condition holds with at least 4 observations for the Bartlett kernel (\(m^{*}=2\)), 5 for the Parzen kernel (\(m^{*}=3\)) and 6 for the quadratic spectral kernel (\(m^{*}=4\)). An error is raised for fewer observations, and if the bandwidth is not finite, which happens if \(\hat{s}^{(0)}\) is zero or the series has values that are not finite.

Constant Detection¶

Whether a model contains a constant or equivalent is tested using three tests. These are executed in order and so once a constant is detected, the others are not executed. The simplest method to ensure that a constant is correctly detected is to include a columns of 1s.

A column with only 1.0s
A column with a maximum minus minimum equal to 0 and that is not all 0s.
Whether the column of 1s, \(1_{N}\), is in the column space of \(X\). If it is, then the model contains a constant equivalent (e.g., dummies for all categories). \(X\) can have any rank. The column of 1s is regressed on the columns of \(X\), which are scaled to have a largest absolute value of 1, using a singular value decomposition. It is in the column space if the norm of the residual, relative to the norm of \(1_{N}\), is at most the square root of the machine precision (about \(1.5\times10^{-8}\)). The test is not used if the rank of \(X\) is \(N\), since the column space is then everything.

References

Sources used in writing the code include [AngImbKru], [Baltagi], [BaumEtAl] and [CamEtAl], [CamTri05], [CamTri09], [ChaoEtAl], [DavMac06], [Greene], [Hall00], [HansenLee], [NewWes94], [Stata], [Wool10] and [Wool12].