From 9ef80d937b3ae35c918ad6bbf450e569c271a43e Mon Sep 17 00:00:00 2001 From: Kevin Sheppard Date: Mon, 12 Oct 2026 12:52:37 +0100 Subject: [PATCH] BUG: Use weighted residuals in the covariance of weighted GLS and 3SLS The GLS fit of SUR and IV3SLS (method="gls", which is the default unless every equation has the same regressors) gave the covariance estimators the unweighted residuals y - x b together with the weighted regressors sqrt(w) x. The scores of the robust, kernel and clustered covariances were then sqrt(w) x e instead of w x e, so their standard errors were wrong whenever the weights were not all equal. The coefficients, the unadjusted covariance and the OLS fit (method="ols", which already used the weighted residuals) were not affected. With one equation, GLS is weighted least squares and the robust covariance is HC0. For a regression with 200 observations and chi-squared weights, the robust standard errors of SUR were 0.1299, 0.1823 and 0.1254, where statsmodels WLS(...).fit(cov_type="HC0") gives 0.1448, 0.1982 and 0.1324. For systems with several equations, 96 of the 140 weighted combinations of GLS or OLS, iterated GLS, full_cov, constraints, debiased and covariance type differed from a direct computation of the stacked sandwich estimator, by up to 28% of the largest covariance element, for SUR and for IV3SLS. After this change all of them agree to within 1e-7, and the unweighted fits are unchanged since the weights are all one. The covariance estimators now get the residuals scaled by the square root of the weights. The tests compare the weighted estimators with the unweighted estimators applied to the data scaled by the square root of the normalized weights, for every covariance type, and the robust standard errors of a weighted regression with the statsmodels WLS estimator. Seven of them fail without the change. --- linearmodels/system/model.py | 5 ++- linearmodels/tests/system/test_3sls.py | 32 ++++++++++++++++++ linearmodels/tests/system/test_sur.py | 47 ++++++++++++++++++++++++++ 3 files changed, 83 insertions(+), 1 deletion(-) diff --git a/linearmodels/system/model.py b/linearmodels/system/model.py index ad5f6a01f7..9fbe2ca992 100644 --- a/linearmodels/system/model.py +++ b/linearmodels/system/model.py @@ -1162,9 +1162,12 @@ def _gls_finalize( cov_estimator = COV_EST[cov_type] gls_eps = np.reshape(gls_eps, (k, gls_eps.shape[0] // k)).T eps = np.reshape(eps, (k, eps.shape[0] // k)).T + # The covariance combines the weighted regressors with the residuals, so the + # residuals must be weighted as well + weps = eps * np.sqrt(np.hstack(self._w)) cov_est = cov_estimator( self._wxhat, - eps, + weps, sigma, full_sigma, gls=True, diff --git a/linearmodels/tests/system/test_3sls.py b/linearmodels/tests/system/test_3sls.py index 7242e12081..1a856437d9 100644 --- a/linearmodels/tests/system/test_3sls.py +++ b/linearmodels/tests/system/test_3sls.py @@ -198,6 +198,38 @@ def test_weighted_against_transformed_data(common_exog): assert_allclose(np.diag(res.cov), np.diag(out.cov)) +@pytest.mark.parametrize("method", ["ols", "gls"]) +@pytest.mark.parametrize("cov_type", ["unadjusted", "robust", "kernel", "clustered"]) +def test_weighted_covariance_is_covariance_of_scaled_data(method, cov_type): + # Weighted estimation is estimation on the data scaled by the square root of the + # normalized weights, so the covariance must be the same for every estimator + data = generate_3sls_data( + n=250, + k=3, + p=3, + en=2, + instr=3, + const=True, + rho=0.8, + common_exog=False, + included_weights=True, + output_dict=True, + ) + scaled = {} + for key, eqn in data.items(): + w = np.sqrt(eqn["weights"] / eqn["weights"].mean()) + scaled[key] = {k: v * w for k, v in eqn.items() if k != "weights"} + kwargs = {} + if cov_type == "kernel": + kwargs["bandwidth"] = 3 + elif cov_type == "clustered": + kwargs["clusters"] = np.random.default_rng(0).integers(0, 25, 250) + res = IV3SLS(data).fit(method=method, cov_type=cov_type, **kwargs) + expected = IV3SLS(scaled).fit(method=method, cov_type=cov_type, **kwargs) + assert_allclose(res.params, expected.params) + assert_allclose(res.cov, expected.cov) + + def test_weighted_just_identified_against_iv2sls(): # With the same instruments in every equation and each equation exactly # identified, 3SLS reduces to equation-by-equation 2SLS, also with weights diff --git a/linearmodels/tests/system/test_sur.py b/linearmodels/tests/system/test_sur.py index 19552409fa..6c74885210 100644 --- a/linearmodels/tests/system/test_sur.py +++ b/linearmodels/tests/system/test_sur.py @@ -11,6 +11,7 @@ import pytest from scipy.sparse.linalg import inv as spinv import scipy.stats +import statsmodels.api as sm from linearmodels.iv.model import _OLS as OLS from linearmodels.shared.hypotheses import InvalidTestStatistic, WaldTestStatistic @@ -806,6 +807,52 @@ def test_iterated_gls_cov_uses_final_sigma(): assert_allclose(res.cov, expected, rtol=1e-6) +def test_weighted_gls_robust_covariance_against_wls(): + # With one equation, GLS is WLS and the robust covariance is HC0 + rng = np.random.default_rng(20261012) + nobs = 200 + x = np.column_stack([np.ones(nobs), rng.standard_normal((nobs, 2))]) + y = x @ np.array([1.0, 0.5, -0.3]) + rng.standard_normal(nobs) * ( + 1 + np.abs(x[:, 1]) + ) + w = rng.chisquare(5, nobs) / 5 + ref = sm.WLS(y, x, weights=w / w.mean()).fit(cov_type="HC0") + eqns = {"a": {"dependent": y[:, None], "exog": x, "weights": w[:, None]}} + for method in ("ols", "gls"): + res = SUR(eqns).fit(method=method, cov_type="robust") + assert_allclose(res.params, ref.params) + assert_allclose(res.std_errors, ref.bse) + + +def scale_by_sqrt_weights(eqns): + """The same equations without weights, with the data scaled by sqrt(weights)""" + scaled = {} + for key, eqn in eqns.items(): + w = np.sqrt(eqn["weights"] / eqn["weights"].mean()) + scaled[key] = {k: v * w for k, v in eqn.items() if k != "weights"} + return scaled + + +@pytest.mark.parametrize("method", ["ols", "gls"]) +@pytest.mark.parametrize("cov_type", ["unadjusted", "robust", "kernel", "clustered"]) +def test_weighted_covariance_is_covariance_of_scaled_data(method, cov_type): + # Weighted estimation is estimation on the data scaled by the square root of the + # normalized weights, so the covariance must be the same for every estimator + eqns = generate_data(k=3, p=3, const=True, rho=0.8, included_weights=True) + nobs = eqns["equ.0"]["dependent"].shape[0] + kwargs = {} + if cov_type == "kernel": + kwargs["bandwidth"] = 3 + elif cov_type == "clustered": + kwargs["clusters"] = np.random.default_rng(0).integers(0, 25, nobs) + res = SUR(eqns).fit(method=method, cov_type=cov_type, **kwargs) + scaled = SUR(scale_by_sqrt_weights(eqns)).fit( + method=method, cov_type=cov_type, **kwargs + ) + assert_allclose(res.params, scaled.params) + assert_allclose(res.cov, scaled.cov) + + def direct_gls(eqns, scale): y = [] x = []