SUR regression and three-stage least squares in Python

Some questions involve several regressions on the same units. A household's spending shares on food and housing respond to the same income shocks; a firm decides how many workers to hire and how much capital to buy at the same time; price and quantity in a market are set jointly by supply and demand. Fitting each equation separately ignores what the equations share. A system estimator uses the correlation between their errors to sharpen the estimates, allows tests and constraints that span equations, and, in the simultaneous case, deals with regressors that are outcomes of another equation.

What SUR and 3SLS estimate

Seemingly unrelated regression (SUR) fits M linear equations, each with its own regressors, whose errors are correlated across equations for the same unit but independent across units. Zellner (1962) proposed feasible generalized least squares: fit each equation by OLS, estimate the M × M error covariance Σ from the residuals, then estimate all equations together:

β̂ = [X′(Σ̂⁻¹ ⊗ I)X]⁻¹ X′(Σ̂⁻¹ ⊗ I)y

Here X is block diagonal with one block of regressors per equation and y stacks the outcomes. Every regressor must be exogenous in every equation, so OLS on each equation is already consistent; SUR is about efficiency. The gain is larger when the errors are strongly correlated and when the equations have different regressors. If the errors are uncorrelated, or if every equation has the same regressors, SUR reproduces OLS equation by equation (Zellner 1962). The Breusch–Pagan test checks whether the error correlations are zero (Breusch and Pagan 1980). Iterating the two steps to convergence gives the Gaussian maximum-likelihood estimator.

Three-stage least squares (3SLS) extends this to simultaneous equations, where some right-hand-side variables are endogenous (Zellner and Theil 1962). It runs two-stage least squares on each equation, estimates Σ from those residuals, and applies the GLS step with the instruments' projection in place of the identity. Like 2SLS, it needs excluded instruments for each equation. The trade-off is that 3SLS is more efficient than 2SLS when the errors are correlated and the equations are overidentified, but a misspecified equation then contaminates the estimates of the others. When every equation is exactly identified, 3SLS equals 2SLS.

In OpenEconometrics

oe.sureg and oe.reg3 take a list of equations. Each is a dictionary with the outcome "y", the regressors "x", an optional label "name" (the outcome by default) and an optional "constant": False. Coefficients are named label:term, such as food:kids, and oe.test and oe.lincom accept those names for cross-equation hypotheses. The estimation sample is casewise over every column of the system.

The SUR model in the example corresponds to the Stata command sureg (food lninc kids) (housing lninc urban), corr. Options of oe.sureg:

• iterate=True iterates to the maximum-likelihood estimate (Stata's isure).
• dfk=True uses the small-sample divisor for Σ̂; small=True reports t and F statistics instead of z and χ².
• constraints=[{"terms": {"food:lninc": 1, "housing:lninc": 1}, "value": 0}] imposes linear restrictions within or across equations.
• weights= with weight_type="aweight" or "fweight".

The Breusch–Pagan test of independent equations is in fit.tests["breusch_pagan"]. When every equation has the same regressors, oe.mvreg gives the OLS coefficients with their joint covariance.

oe.reg3 takes the same equation list. Every equation's outcome is endogenous; endogenous= adds right-hand-side variables that are no equation's outcome, exogenous= adds instruments that appear in no equation, and instruments= gives the complete list of exogenous variables instead. The system in the example corresponds to reg3 (demand: q p income pop) (supply: q p cost rain), endog(p). method="3sls" is the default; "2sls" and "ols" estimate equation by equation, "sure" treats every regressor as exogenous, and ireg3=True iterates the covariance. Each equation must meet the rank condition; otherwise the fit stops with underidentified. The systems family guide gives every formula and convention.

Supported options

Options of oe.sureg

Options of oe.reg3

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

Example

The script has two parts. The first simulates 600 households whose food share falls with log income (true slope −0.06) and rises with the number of children (0.015), and whose housing share rises with log income (0.04) and urban residence (0.08). The two errors have a correlation of −0.6. The second part simulates 400 regional markets. Demand is q = 10 − 1.0 · p + 0.8 · income + 0.5 · pop + u_d and supply is q = 2 + 0.5 · p − 0.6 · cost + 0.4 · rain + u_s, with error correlation 0.5. Price clears the market, so it is correlated with both errors; income and population shift only demand, cost and rain only supply.

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


def budget_shares(rng, n=600):
    """Food and housing budget shares with correlated errors and different regressors."""
    lninc = rng.normal(10, 0.5, size=n)
    kids = rng.integers(0, 4, size=n)
    urban = rng.integers(0, 2, size=n)
    e = rng.multivariate_normal([0, 0], [[0.04**2, -0.6 * 0.04 * 0.03],
                                         [-0.6 * 0.04 * 0.03, 0.03**2]], size=n)
    food = 0.90 - 0.06 * lninc + 0.015 * kids + e[:, 0]
    housing = -0.20 + 0.04 * lninc + 0.08 * urban + e[:, 1]
    return pd.DataFrame({"food": food, "housing": housing, "lninc": lninc,
                         "kids": kids, "urban": urban})


def market(rng, n=400):
    """Price and quantity set jointly by demand and supply in n regional markets."""
    income, pop, cost, rain = rng.normal(size=(4, n))
    u = rng.multivariate_normal([0, 0], [[1, 0.5], [0.5, 1]], size=n)
    demand_shift = 10 + 0.8 * income + 0.5 * pop + u[:, 0]  # q = shift - 1.0 p
    supply_shift = 2 - 0.6 * cost + 0.4 * rain + u[:, 1]    # q = shift + 0.5 p
    p = (demand_shift - supply_shift) / 1.5
    q = supply_shift + 0.5 * p
    return pd.DataFrame({"q": q, "p": p, "income": income, "pop": pop,
                         "cost": cost, "rain": rain})


if __name__ == "__main__":
    rng = np.random.default_rng(2026)
    shares = budget_shares(rng)
    sur = oe.sureg(data=shares, equations=[
        {"y": "food", "x": ["lninc", "kids"]},
        {"y": "housing", "x": ["lninc", "urban"]},
    ])
    print(sur.summary())
    ols_food = oe.ols(data=shares, y="food", x=["lninc", "kids"])
    ols_housing = oe.ols(data=shares, y="housing", x=["lninc", "urban"])
    print(f"OLS se: kids {ols_food.coefficients[2].std_error:.5f}, "
          f"urban {ols_housing.coefficients[2].std_error:.5f}")
    print()

    df = market(rng)
    eqs = [{"y": "q", "x": ["p", "income", "pop"], "name": "demand"},
           {"y": "q", "x": ["p", "cost", "rain"], "name": "supply"}]
    fits = {"OLS": oe.reg3(data=df, equations=eqs, method="ols"),
            "2SLS": oe.reg3(data=df, equations=eqs, endogenous=["p"], method="2sls"),
            "3SLS": oe.reg3(data=df, equations=eqs, endogenous=["p"])}
    print("Price coefficient   demand (true -1.0)    supply (true 0.5)")
    for label, fit in fits.items():
        c = {k.term: k for k in fit.coefficients}
        d, s = c["demand:p"], c["supply:p"]
        print(f"{label:<6}{d.estimate:>18.4f} ({d.std_error:.4f}){s.estimate:>14.4f} "
              f"({s.std_error:.4f})")

With openecon 0.3.20a1 from PyPI the script prints:

Seemingly unrelated regression — food
Observations: 600  |  Covariance: nonrobust  |  Confidence: 95%

Term                 Estimate  Std. error         z      P>|stat|    CI lower    CI upper
-----------------  ----------  ----------  --------  ------------  ----------  ----------
[food]
food:Intercept       0.926035   0.0313318   29.5558  5.53118e-192    0.864626    0.987444
food:lninc         -0.0623661  0.00311562  -20.0172    3.8978e-89  -0.0684726  -0.0562596
food:kids           0.0152176   0.0011689   13.0188   9.57235e-39   0.0129266   0.0175086
[housing]
housing:Intercept   -0.224107    0.023576  -9.50573   1.98653e-21   -0.270315   -0.177899
housing:lninc       0.0422899  0.00234423     18.04    9.4464e-73   0.0376953   0.0468845
housing:urban       0.0773994  0.00195821   39.5257             0   0.0735614   0.0812374

n_equations: 2  |  food:rmse: 0.0400975  |  food:r_squared: 0.452828  |  housing:rmse: 0.0301455
housing:r_squared: 0.678225  |  log_likelihood: 2468.23  |  df_model: 4  |  df_resid: 597
Wald chi2 test of all slopes of the system: chi2(4) = 2156.14, p = 0
Wald chi2 test of the slopes of equation food: chi2(2) = 573.895, p = 2.401e-125
Wald chi2 test of the slopes of equation housing: chi2(2) = 1846.73, p = 0
Breusch-Pagan test of independent equations: chi2(1) = 223.075, p = 1.93e-50
OLS se: kids 0.00148, urban 0.00248

Price coefficient   demand (true -1.0)    supply (true 0.5)
OLS              -0.5624 (0.0573)        0.0897 (0.0532)
2SLS             -0.9771 (0.0955)        0.4521 (0.0793)
3SLS             -0.9765 (0.0950)        0.4527 (0.0789)

Reading the output

• SUR recovers the budget-share coefficients. Every 95% interval contains its true value: −0.0624 for log income and 0.0152 for children in the food equation, 0.0423 for log income and 0.0774 for urban residence in the housing equation.
• The correlated errors buy precision. The standard error of the children coefficient falls from 0.00148 with OLS to 0.00117 with SUR, and that of urban residence from 0.00248 to 0.00196, both about 20% smaller. These regressors appear in only one equation, so the other equation's residuals carry information about their errors. With identical regressors in both equations there would be no gain.
• The Breusch–Pagan test confirms the correlation. χ²(1) = 223.1 rejects independent equations decisively, as expected with a true correlation of −0.6.
• OLS is badly biased in the market system. Because price responds to both demand and supply shocks, OLS gives −0.5624 for the demand slope and 0.0897 for the supply slope, against true values of −1.0 and 0.5.
• 2SLS and 3SLS recover the slopes; the system step adds almost nothing here. 2SLS gives −0.9771 and 0.4521, 3SLS −0.9765 and 0.4527. The 3SLS standard errors are about 0.5% smaller, which is roughly what the N − k divisor of method="2sls" accounts for on its own (√(400/396) ≈ 1.005). With one overidentifying instrument per equation and a moderate error correlation, there is little extra information for 3SLS to use.

Limitations

• oe.sureg and oe.reg3 report the conventional covariance only; robust, cluster, bootstrap and jackknife variances are not available.
• oe.reg3 has no equivalent of Stata's allexog, corr(), dfk2 or first options, and oe.sureg none of dfk2.
• small=True changes only the reference distributions, using the first equation's degrees of freedom for t statistics; the covariance divisor changes only with dfk=True.
• Iterated SUR is refused with singular_sigma when an outcome is also a regressor in another equation, because the likelihood is then unbounded; use the two-step estimator or oe.reg3.
• oe.reg3 reports structural residuals, computed with the observed endogenous regressors, so an equation's R² can be negative.
• No numerical parity with Stata is claimed: the inventory's stata_parity_validated flag is False.

References

• Breusch, T. S. and Pagan, A. R. (1980). The Lagrange multiplier test and its applications to model specification in econometrics. Review of Economic Studies, 47(1), 239–253. doi:10.2307/2297111
• Zellner, A. (1962). An efficient method of estimating seemingly unrelated regressions and tests for aggregation bias. Journal of the American Statistical Association, 57(298), 348–368. doi:10.1080/01621459.1962.10480664
• Zellner, A. and Theil, H. (1962). Three-stage least squares: Simultaneous estimation of simultaneous equations. Econometrica, 30(1), 54–78. doi:10.2307/1911287