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 = []