Instrumental variables in Python: 2SLS, LIML and GMM

Ordinary least squares is consistent only when the regressors are uncorrelated with the error term. When a regressor is driven by something the model leaves out, such as ability in a wage equation, or when it is measured with error or determined jointly with the outcome, OLS is biased however large the sample. Instrumental variables (IV) recover the effect from the part of the regressor's variation that an outside variable, the instrument, creates.

What instrumental variables estimate

An instrument must satisfy two conditions. It must be relevant: it moves the endogenous regressor. It must be excluded: it affects the outcome only through that regressor, so it is uncorrelated with the error. Relevance can be checked in the data. Exclusion cannot; it has to be argued from how the instrument arises.

Two-stage least squares (2SLS) makes the idea concrete. The first stage regresses the endogenous variable on the instruments and the exogenous controls. The second stage regresses the outcome on the fitted values from the first stage and the same controls. Only instrument-driven variation is left in the fitted values, so the second-stage coefficient is consistent under the two conditions.

• LIML (limited-information maximum likelihood) is a k-class estimator that is less biased than 2SLS when instruments are weak or many. In an exactly identified model it equals 2SLS.
• GMM (generalized method of moments) weights the instrument moments by their estimated covariance. When the model has more instruments than endogenous regressors and the errors are heteroskedastic, two-step GMM is more efficient than 2SLS (Hansen 1982).

When the effect differs across people, IV with a binary instrument and a binary treatment estimates a local average treatment effect: the effect for those whose treatment the instrument changes (Imbens and Angrist 1994). That group need not be representative of the population.

Weak instruments are the main practical risk. If the instruments explain little of the endogenous regressor, 2SLS is biased towards OLS and its usual confidence intervals are unreliable (Staiger and Stock 1997). Andrews, Stock and Sun (2019) review the tests and confidence sets that remain valid in that case.

In OpenEconometrics

oe.ivregress fits all three estimators. The exogenous regressors go in x, the endogenous ones in endog and the excluded instruments in instruments. The exogenous regressors act as their own instruments and are not repeated. The model below corresponds to the Stata command ivregress 2sls lwage female (educ = distance tuition), vce(robust):

• method="2sls" (the default), "liml" or "gmm", plus "fuller" (adjusted LIML) and "kclass" with an explicit kappa.
• covariance="robust", "cluster" with cluster= (one or two columns) or "hac" with lags=. Inference uses z and chi-square statistics by default; small=True switches to t and F with N − K degrees of freedom.
• One fit reports the first-stage, weak-identification, overidentification and endogeneity statistics together, the results that Stata produces with estat firststage, estat overid and estat endogenous.

For weak-instrument-robust inference, oe.iv_weak_test runs Anderson–Rubin and conditional likelihood-ratio tests of a hypothesized coefficient, and oe.iv_ar_confidence_set inverts the Anderson–Rubin test. The IV method guide gives every formula and convention.

Supported options

Options of oe.ivregress

oe.iv_weak_test

Stata command: none listed.

oe.iv_ar_confidence_set

Stata command: none listed.

oe.stock_yogo

Stata command: none listed.

Generated from oe.capabilities() of openecon 0.3.20a1. Stata names identify the corresponding workflow, not numerical parity.

Example

The script simulates 5,000 workers whose unobserved ability raises both their schooling and their wages, so OLS overstates the return to schooling. Distance to the nearest college and local tuition shift schooling but, by construction, affect wages only through schooling. The true return to a year of schooling is 0.08.

import numpy as np
import pandas as pd
import openecon as oe


def simulate(n=5000, seed=7):
    """Log wages with unobserved ability: schooling is endogenous, the instruments are not."""
    rng = np.random.default_rng(seed)
    ability = rng.normal(size=n)
    distance = rng.exponential(scale=1.0, size=n)  # distance to the nearest college
    tuition = rng.normal(size=n)                    # local tuition at age 18
    female = rng.integers(0, 2, size=n)
    educ = 13 + 1.2 * ability - 0.8 * distance - 0.5 * tuition + rng.normal(size=n)
    lwage = 1.0 + 0.08 * educ - 0.10 * female + 0.4 * ability + rng.normal(scale=0.5, size=n)
    return pd.DataFrame({"lwage": lwage, "educ": educ, "female": female,
                         "distance": distance, "tuition": tuition})


if __name__ == "__main__":
    df = simulate()
    ols = oe.ols(data=df, y="lwage", x=["educ", "female"], covariance="HC1")
    iv = oe.ivregress(data=df, y="lwage", x=["female"], endog=["educ"],
                      instruments=["distance", "tuition"], covariance="robust")
    print(f"OLS return to schooling: {ols.coefficients[1].estimate:.4f}")
    print(iv.summary())

With openecon 0.3.20a1 from PyPI the script prints:

OLS return to schooling: 0.2225
Instrumental-variables (2SLS) regression — lwage
Observations: 5000  |  Covariance: robust  |  Confidence: 95%

Term         Estimate  Std. error         z     P>|stat|   CI lower    CI upper
---------  ----------  ----------  --------  -----------  ---------  ----------
Intercept     1.03963    0.117923   8.81617  1.18439e-18   0.808503     1.27075
female     -0.0808504   0.0183723  -4.40068  1.07913e-05  -0.116859  -0.0448415
educ        0.0754119  0.00960766   7.84914  4.18889e-15  0.0565812   0.0942425

r_squared: 0.181967  |  adjusted_r_squared: 0.18164  |  rmse: 0.649729  |  df_model: 2
df_resid: 4997  |  n_instruments: 2  |  n_endogenous: 1
Wald chi2 test of the slopes: chi2(2) = 81.6302, p = 1.88e-18
Wooldridge robust score test of overidentifying restrictions: chi2(1) = 0.366141, p = 0.5451
Wooldridge robust score chi2 test of exogeneity: chi2(1) = 299.978, p = 3.33e-67
Robust regression-based F test of exogeneity: F(1, 4996) = 440.042, p = 1.047e-93
Cragg-Donald Wald F (minimum eigenvalue statistic): F(2, 4996) = 981.797
Kleibergen-Paap rk Wald F: F(2, 4996) = 950.273

Reading the output

• OLS is badly biased. It estimates 0.2225, almost three times the true 0.08, because ability sits in the error term and is correlated with schooling.
• 2SLS recovers the effect. The estimate is 0.0754 with a 95% confidence interval from 0.057 to 0.094, which contains 0.08. The interval is about twice as wide as the OLS one would be; IV pays for consistency with precision.
• The instruments are strong. The Cragg–Donald F of 981.8 and the Kleibergen–Paap rk Wald F of 950.3 measure how well the instruments predict schooling. The Kleibergen–Paap version is robust to heteroskedasticity. Neither has a p-value; compare them with weak-instrument critical values, not with 0.05.
• The overidentification test does not reject. With two instruments and one endogenous regressor there is one overidentifying restriction. Wooldridge's robust score test gives χ²(1) = 0.37 with p = 0.55. This shows only that the two instruments give similar answers. If both violated the exclusion restriction in the same way, the test would not notice.
• Schooling is endogenous. The robust endogeneity tests reject exogeneity decisively (χ²(1) = 300.0), so the OLS and IV estimates differ by more than sampling error.

Limitations

• Cluster-robust covariance for oe.ivregress allows at most 2 clustering dimensions.
• The Cragg–Donald and Kleibergen–Paap statistics are reported without p-values. oe.stock_yogo supplies Stock–Yogo reference values only for one endogenous regressor in unweighted 2SLS with iid errors.
• LIML fits report no endogeneity test, and Fuller and k-class fits report no overidentification or endogeneity test.
• A GMM weight matrix clustered on fewer clusters than instruments is singular and is rejected.
• The conditional likelihood-ratio test in oe.iv_weak_test covers one endogenous coefficient with iid errors; robust, cluster and HAC versions are rejected.
• No test can establish the exclusion restriction. It rests on the argument for the instrument.
• No numerical parity with Stata is claimed: the inventory's stata_parity_validated flag is False.

References

• Andrews, I., Stock, J. H. and Sun, L. (2019). Weak instruments in instrumental variables regression: Theory and practice. Annual Review of Economics, 11, 727–753. doi:10.1146/annurev-economics-080218-025643
• Hansen, L. P. (1982). Large sample properties of generalized method of moments estimators. Econometrica, 50(4), 1029–1054. doi:10.2307/1912775
• Imbens, G. W. and Angrist, J. D. (1994). Identification and estimation of local average treatment effects. Econometrica, 62(2), 467–475. doi:10.2307/2951620
• Staiger, D. and Stock, J. H. (1997). Instrumental variables regression with weak instruments. Econometrica, 65(3), 557–586. doi:10.2307/2171753