From 727f80c6dab4d22af8546ac96fe5aa846ae116c5 Mon Sep 17 00:00:00 2001 From: Kevin Sheppard Date: Mon, 12 Oct 2026 11:30:48 +0100 Subject: [PATCH 1/5] BUG: Raise a ValueError if IVGMM has a singular moment covariance The weight matrix of two-step GMM is the inverse of an estimate of the covariance of the moment conditions. If there are fewer clusters than moment conditions (or too few observations) the estimate is singular, the inverse is numerical noise, and so are the estimates. #754 made j_stat and c_stat invalid and made IVGMMCUE raise a clear ValueError in this case, and left the estimates of IVGMM as they were: * In the model of the bug report with 3 clusters, the coefficient on x is 1.697 when it is 0.5, and its standard error is NaN (and 1.95e6 for the constant); it is a LinAlgError instead if NumPy finds that the matrix is exactly singular, which depends on rounding and so on the platform. * An exactly identified model does not need its weight matrix, but the noise in it is not harmless: with 3 clusters and 3 moment conditions the estimates are 9.7 away from the IV estimates, and the standard errors are not numbers. The first-step moment conditions of an exactly identified model add to zero, so that the covariance of 3 moment conditions is singular with 3 clusters even if they are not centered. IVGMM.fit now raises the same ValueError as IVGMMCUE.fit, with the message built by the same function: it says that the covariance is singular, that there are too few clusters (with their number and the rank that they allow) or too few observations, and to use more clusters than moment conditions. It is checked at every step before the covariance is inverted. It is not checked with iter_limit=1, which does not use a covariance. IVGMMResults.c_stat estimates a model with a moment condition for each tested variable, and returns an invalid statistic instead of an error if that model cannot be estimated, as it did for a LinAlgError. j_stat of an IVGMM is never invalid now, since the model is not estimated. The fits that are valid are not changed. Params, cov, std_errors, j_stat, centered_j_stat and c_stat of 600 random models (IVGMM and IVGMMCUE, the four weight types, observation weights, three covariances, units of the instruments from 1e-2 to 1e2, and cluster counts around the number of moment conditions) are the same to the last bit on origin/main and here for the 561 that are estimated on both. The 39 that raise here are clustered models with no more clusters than moment conditions: 19 raised on origin/main too, and the 20 that did not had an invalid J statistic, and NaN or enormous standard errors. Tests: the tests of #754 for the J and C statistics with too few clusters (and the small samples) accepted a result or a LinAlgError from fit, and returned if there was none. They test the ValueError and its message now, and the values of the statistics for the models that can be estimated. New tests check that the message is that of IVGMMCUE except for the sentence that says why the inverse is needed; that exactly identified models raise with 2 and 3 clusters and are the IV estimates with 4 or more; that every step is checked, and the initial weight does not avoid the check; that iter_limit=1 is not changed; and that the params, std errors, J, centered J and C of seven valid models are those of origin/main (and the params and J are the definitions). --- linearmodels/iv/model.py | 115 ++++-- linearmodels/iv/results.py | 6 +- .../tests/iv/test_gmm_singular_covariance.py | 381 +++++++++++++++--- 3 files changed, 412 insertions(+), 90 deletions(-) diff --git a/linearmodels/iv/model.py b/linearmodels/iv/model.py index 38d72047ee..4186866d3f 100644 --- a/linearmodels/iv/model.py +++ b/linearmodels/iv/model.py @@ -1196,6 +1196,58 @@ def _j_statistic( null = "Expected moment conditions are equal to 0" return WaldTestStatistic(stat, null, ninstr - nvar) + def _raise_if_singular( + self, + moment_cov: linearmodels.typing.data.Float64Array, + consequence: str, + ) -> None: + """ + Raise a ValueError if the covariance of the moment conditions is singular + + Parameters + ---------- + moment_cov : ndarray + Estimate of the covariance of the moment conditions (ninstr by + ninstr). + consequence : str + Sentence that says why the estimator needs the inverse of the + covariance. It ends with the period. + + Raises + ------ + ValueError + If the covariance is singular. The message says that there are too + few clusters or observations and how many are needed. + """ + if not _is_singular(moment_cov): + return + nobs, nmoment = self._wz.shape + if isinstance(self._weight, OneWayClusteredWeightMatrix): + nclusters = len(unique(asarray(self._weight_config["clusters"]))) + centered = self._weight_config["center"] + rank = nclusters - 1 if centered else nclusters + rank_reason = ( + "one less than the number of clusters, since the cluster sums of " + "centered moment conditions add to zero" + if centered + else "the number of clusters" + ) + msg = ( + f"The clustered covariance of the {nmoment} moment conditions is " + f"singular because there are too few clusters ({nclusters} " + f"clusters). Its rank is at most {rank}, {rank_reason}. " + f"{consequence} Use more clusters than moment conditions, or a " + "weight_type that does not use clusters." + ) + else: + msg = ( + f"The covariance of the {nmoment} moment conditions is singular, " + f"which occurs if there are too few observations ({nobs}) for the " + f"number of moment conditions. {consequence} Use more " + "observations or fewer moment conditions." + ) + raise ValueError(msg) + class IVGMM(_IVGMMBase): r""" @@ -1251,6 +1303,17 @@ class IVGMM(_IVGMMBase): paired by position. The same rule applies to ``clusters`` passed to ``fit`` or in ``weight_config``. + The weight matrix is the inverse of an estimate of the covariance of the + moment conditions, so that the estimate has to be nonsingular. A clustered + estimate from G clusters has rank at most G, and at most G - 1 if the + moment conditions are centered (``center`` is ``True``), since the + centered cluster sums add to zero. There have to be more clusters than + moment conditions. The rank of a kernel or robust estimate is bounded in + the same way by the number of observations. ``fit`` raises a + ``ValueError`` that says so if the estimate is singular, since the + estimates would be numerical noise. The estimate is not used, and so it + is not tested, if ``iter_limit`` is 1. + .. todo:: * VCV: bootstrap @@ -1425,6 +1488,15 @@ def fit( IVGMMResults Results container + Raises + ------ + ValueError + If the estimate of the covariance of the moment conditions is + singular, which occurs if there are too few clusters or + observations. The number of clusters has to be larger than the + number of moment conditions. The estimate is not used if + ``iter_limit`` is 1. + See also -------- linearmodels.iv.gmm.IVGMMCovariance @@ -1448,6 +1520,11 @@ def fit( while iters < iter_limit and norm > tol: eps = wy - wx @ params moment_cov = weight_matrix(wx, wz, eps) + self._raise_if_singular( + moment_cov, + "The weight matrix of the GMM estimator is the inverse of the " + "covariance, so it is not defined.", + ) wmat = inv(moment_cov) params = self.estimate_parameters(wx, wy, wz, wmat) delta = params - _params @@ -1779,39 +1856,13 @@ def _check_moment_covariance(self) -> None: estimates of two-step GMM. """ wy, wx, wz = self._wy, self._wx, self._wz - nobs, nmoment = wz.shape + nobs = wz.shape[0] params = IVGMM.estimate_parameters(wx, wy, wz, inv(wz.T @ wz / nobs)) - if not _is_singular(self._weight.weight_matrix(wx, wz, wy - wx @ params)): - return - if isinstance(self._weight, OneWayClusteredWeightMatrix): - nclusters = len(unique(asarray(self._weight_config["clusters"]))) - centered = self._weight_config["center"] - rank = nclusters - 1 if centered else nclusters - rank_reason = ( - "one less than the number of clusters, since the cluster sums of " - "centered moment conditions add to zero" - if centered - else "the number of clusters" - ) - msg = ( - f"The clustered covariance of the {nmoment} moment conditions is " - f"singular because there are too few clusters ({nclusters} " - f"clusters). Its rank is at most {rank}, {rank_reason}. The " - "objective function of the continuously updating GMM estimator " - "uses the inverse of the covariance, so it is not defined. Use " - "more clusters than moment conditions, or a weight_type that does " - "not use clusters." - ) - else: - msg = ( - f"The covariance of the {nmoment} moment conditions is singular, " - f"which occurs if there are too few observations ({nobs}) for the " - "number of moment conditions. The objective function of the " - "continuously updating GMM estimator uses the inverse of the " - "covariance, so it is not defined. Use more observations or fewer " - "moment conditions." - ) - raise ValueError(msg) + self._raise_if_singular( + self._weight.weight_matrix(wx, wz, wy - wx @ params), + "The objective function of the continuously updating GMM estimator " + "uses the inverse of the covariance, so it is not defined.", + ) def fit( self, diff --git a/linearmodels/iv/results.py b/linearmodels/iv/results.py index cd0fda3123..a11c936933 100644 --- a/linearmodels/iv/results.py +++ b/linearmodels/iv/results.py @@ -1895,8 +1895,10 @@ def c_stat( ) try: res_e = mod.fit(cov_type=self.cov_type, **self.cov_config) - except LinAlgError: - # The covariance of the moment conditions is exactly singular + except (LinAlgError, ValueError): + # The covariance of the moment conditions is singular: fit raises + # a ValueError, and NumPy raises a LinAlgError if it is exactly + # singular. LinAlgError is not a ValueError in NumPy 1.x return InvalidTestStatistic(reason, name="C-statistic") assert isinstance(res_e, IVGMMResults) if isinstance(res_e.j_stat, InvalidTestStatistic): diff --git a/linearmodels/tests/iv/test_gmm_singular_covariance.py b/linearmodels/tests/iv/test_gmm_singular_covariance.py index 2561ccb32c..0f186283c0 100644 --- a/linearmodels/tests/iv/test_gmm_singular_covariance.py +++ b/linearmodels/tests/iv/test_gmm_singular_covariance.py @@ -18,13 +18,15 @@ The objective function of IVGMMCUE is a function of the inverse of the covariance, and so it is not defined either. IVGMMCUE.fit has to raise a -ValueError that says why, before it optimizes the function. - -Whether NumPy raises a LinAlgError when it inverts a covariance that is -singular depends on rounding, and so on the platform. It does so if the matrix -is exactly singular, and otherwise returns the inverse of the noise. The -estimates are not changed by the statistics, and a model whose estimation fails -has no results to test. +ValueError that says why, before it optimizes the function. So does IVGMM.fit, +since the estimates of two-step GMM are numerical noise when the weight matrix +is the inverse of a singular covariance. Whether NumPy raises a LinAlgError when +it inverts a covariance that is singular depends on rounding, and so on the +platform: it does so if the matrix is exactly singular, and otherwise returns +the inverse of the noise. Both fits raise the same ValueError on every platform, +and a model that cannot be estimated has no statistics to test. The C statistic +estimates another model, which has more moment conditions, and is invalid if +that model cannot be estimated. """ import numpy as np @@ -34,7 +36,7 @@ import pytest from scipy import stats -from linearmodels.iv import IVGMM, IVGMMCUE +from linearmodels.iv import IV2SLS, IVGMM, IVGMMCUE import linearmodels.iv.model from linearmodels.shared.hypotheses import InvalidTestStatistic, WaldTestStatistic @@ -172,11 +174,7 @@ def nonsingular(nmoments, nclusters, center): def fit(data, nclusters, center, model=IVGMM, weights=None, **kwargs): - """ - Model of the data with the clusters arange(nobs) % nclusters - - Returns None if estimation fails because the covariance is exactly singular. - """ + """Fit the model of the data with the clusters arange(nobs) % nclusters""" y, exog, endog, instr = data clusters = np.arange(y.shape[0]) % nclusters mod = model( @@ -189,12 +187,22 @@ def fit(data, nclusters, center, model=IVGMM, weights=None, **kwargs): clusters=clusters, center=center, ) - try: - return mod.fit(cov_type="clustered", clusters=clusters, **kwargs) - except LinAlgError: - # Only a singular covariance can be the reason - assert not nonsingular(exog.shape[1] + instr.shape[1], nclusters, center) - return None + return mod.fit(cov_type="clustered", clusters=clusters, **kwargs) + + +def check_too_few_clusters(data, nclusters, center, **kwargs): + """The covariance of the moments of the model is singular and fit raises""" + nmoments = data[1].shape[1] + data[3].shape[1] + assert not nonsingular(nmoments, nclusters, center) + with pytest.raises(ValueError, match="too few clusters") as exc_info: + fit(data, nclusters, center, **kwargs) + # LinAlgError is a ValueError in recent versions of NumPy + assert exc_info.type is ValueError + message = str(exc_info.value) + assert f"({nclusters} clusters)" in message + assert f"{nmoments} moment conditions" in message + assert "singular" in message + assert "Use more clusters than moment conditions" in message def check_invalid(stat, name): @@ -228,12 +236,12 @@ def test_c_stat_with_too_few_clusters_is_invalid(nclusters): # With 3 and 4 clusters the statistic was -3.7e15 and -7.4e16. The model has 4 # moments and the model for the statistic has 5, and so the covariance of the # moments that it estimates is singular unless there are at least 5 clusters - res = fit(ONE_ENDOG, nclusters, center=False) - if res is None: - # With 3 clusters the covariance of the model's own 4 moments is singular - # too, and NumPy 1.x raised when the model was estimated + if nclusters == 3: + # The covariance of the model's own 4 moments is singular too, so that + # the model is not estimated + check_too_few_clusters(ONE_ENDOG, nclusters, center=False) return - check_invalid(res.c_stat(), "C-statistic") + check_invalid(fit(ONE_ENDOG, nclusters, center=False).c_stat(), "C-statistic") def test_c_stat_with_enough_clusters_is_unchanged(): @@ -261,10 +269,10 @@ def test_c_stat_with_enough_clusters_is_unchanged(): @pytest.mark.parametrize("nclusters", [2, 3, 4, 5, 6, 8]) @pytest.mark.parametrize(("data", "tested"), C_STAT_MODELS) def test_c_stat_few_clusters(data, tested, nclusters, center): - res = fit(data, nclusters, center) - if res is None: + if not nonsingular(data[1].shape[1] + data[3].shape[1], nclusters, center): + check_too_few_clusters(data, nclusters, center) return - c_stat = res.c_stat(tested) + c_stat = fit(data, nclusters, center).c_stat(tested) ntested = data[2].shape[1] if tested is None else len(tested) nmoments = data[1].shape[1] + data[3].shape[1] + ntested if not nonsingular(nmoments, nclusters, center): @@ -320,14 +328,12 @@ def raising_inv(a): @pytest.mark.parametrize("center", [False, True]) @pytest.mark.parametrize("nclusters", [2, 3, 4, 5, 6, 8]) def test_j_stat_few_clusters(nclusters, center): - res = fit(ONE_ENDOG, nclusters, center) - if res is None: - return - j_stat = res.j_stat if not nonsingular(4, nclusters, center): - # Between -3e15 and 2e15, or a plausible number in a few cases - check_invalid(j_stat, J_NAME) + # The statistic was between -3e15 and 2e15, or a plausible number in a + # few cases, and there was a LinAlgError in some + check_too_few_clusters(ONE_ENDOG, nclusters, center) return + j_stat = fit(ONE_ENDOG, nclusters, center).j_stat x, y, z, _, _ = arrays(ONE_ENDOG) cov = clustered_cov(np.arange(120) % nclusters, center) check_valid(j_stat, reference_j(x, y, z, cov), df=1) @@ -339,10 +345,11 @@ def test_j_stat_few_clusters(nclusters, center): def test_j_and_centered_j_stat_are_checked_separately(nclusters, center): # The statistic that is not centered is valid with as many clusters as moments, # whereas the centered statistic always centers the moments - res = fit(ONE_ENDOG, nclusters, center) - if res is None: + if not nonsingular(4, nclusters, center): + check_too_few_clusters(ONE_ENDOG, nclusters, center) return - assert (type(res.j_stat) is WaldTestStatistic) == nonsingular(4, nclusters, center) + res = fit(ONE_ENDOG, nclusters, center) + assert type(res.j_stat) is WaldTestStatistic assert (type(res.centered_j_stat) is WaldTestStatistic) == (nclusters >= 5) @@ -358,20 +365,20 @@ def test_j_stat_without_a_covariance_estimate_is_unchanged(): @pytest.mark.parametrize("iter_limit", [2, 3, 100]) def test_j_stat_iterated(iter_limit): - res = fit(ONE_ENDOG, 3, center=False, iter_limit=iter_limit) - if res is not None: - check_invalid(res.j_stat, J_NAME) + check_too_few_clusters(ONE_ENDOG, 3, center=False, iter_limit=iter_limit) res = fit(ONE_ENDOG, 20, center=False, iter_limit=iter_limit) assert type(res.j_stat) is WaldTestStatistic assert np.isfinite(res.j_stat.stat) -@pytest.mark.filterwarnings("ignore:invalid value encountered in sqrt:RuntimeWarning") def test_summary_and_attributes_with_an_invalid_statistic(): - res = fit(ONE_ENDOG, 3, center=False) - if res is None: - return - check_invalid(res.j_stat, J_NAME) + # The model can be estimated with 4 clusters, but not its C statistic or its + # centered J statistic + res = fit(ONE_ENDOG, 4, center=False) + check_invalid(res.c_stat(), "C-statistic") + check_invalid( + res.centered_j_stat, "Centered J-test of overidentifying restrictions" + ) assert "Parameter Estimates" in str(res.summary) assert np.all(np.isfinite(res.params)) assert res.weight_type == "clustered" @@ -389,14 +396,13 @@ def test_small_samples_j_stat(weight_type, nobs): # nobs - 1, and so the covariance of the 4 moments is singular if nobs = 4 data = small_sample(nobs) mod = IVGMM(*data, weight_type=weight_type, center=True) - try: - j_stat = mod.fit().j_stat - except LinAlgError: - assert nobs == 4 - return if nobs == 4: - check_invalid(j_stat, J_NAME) + with pytest.raises(ValueError, match="too few observations") as exc_info: + mod.fit() + assert exc_info.type is ValueError + assert "4 moment conditions" in str(exc_info.value) return + j_stat = mod.fit().j_stat x, y, z, _, _ = arrays(data) bandwidth = nobs - 2 if weight_type == "kernel" else 0 expected = reference_j(x, y, z, bartlett_cov(bandwidth, center=True)) @@ -429,16 +435,14 @@ def test_small_samples_c_stat(weight_type, nobs, center): def test_weights(nclusters): rng = np.random.default_rng(1234) weights = np.exp(0.5 * rng.standard_normal(120)) - res = fit(ONE_ENDOG, nclusters, center=False, weights=weights) - if res is None: + if nclusters == 3: + check_too_few_clusters(ONE_ENDOG, nclusters, False, weights=weights) return + res = fit(ONE_ENDOG, nclusters, center=False, weights=weights) sw = np.sqrt(weights / weights.mean()) x, y, z, z_e, original = arrays(ONE_ENDOG, sqrt_weights=sw) cov = clustered_cov(np.arange(120) % nclusters, center=False) - if nclusters >= 4: - check_valid(res.j_stat, reference_j(x, y, z, cov), df=1) - else: - check_invalid(res.j_stat, J_NAME) + check_valid(res.j_stat, reference_j(x, y, z, cov), df=1) if nclusters >= 5: check_valid(res.c_stat(), reference_c(x, y, z_e, original, cov), df=1) else: @@ -574,3 +578,268 @@ def test_cue_too_few_observations(weight_type, nobs): mod.fit() assert exc_info.type is ValueError assert "4 moment conditions" in str(exc_info.value) + + +# ---------------------------------------------------------------------------- +# IVGMM does not estimate a model with a singular covariance, and estimates the +# other models as it did +# ---------------------------------------------------------------------------- + +IVGMM_CONSEQUENCE = ( + "The weight matrix of the GMM estimator is the inverse of the covariance, so " + "it is not defined." +) +CUE_CONSEQUENCE = ( + "The objective function of the continuously updating GMM estimator uses the " + "inverse of the covariance, so it is not defined." +) + + +def exactly_identified_data(nobs=120): + """2 exog, 1 endog and 1 instrument: 3 moments and 3 parameters""" + rng = np.random.default_rng(3) + z = pd.DataFrame(rng.standard_normal((nobs, 1)), columns=["z0"]) + x = pd.DataFrame({"x": z.sum(axis=1) + rng.standard_normal(nobs)}) + exog = pd.DataFrame({"const": 1.0, "w": rng.standard_normal(nobs)}) + y = pd.Series(1 + 0.5 * x.x + rng.standard_normal(nobs), name="y") + return y, exog, x, z + + +EXACTLY_IDENTIFIED = exactly_identified_data() + + +@pytest.mark.parametrize("center", [False, True]) +@pytest.mark.parametrize("nclusters", [2, 3]) +def test_ivgmm_message_has_the_form_of_the_cue_message(nclusters, center): + # The same message, except for the sentence that says why the inverse is needed + with pytest.raises(ValueError, match="too few clusters") as ivgmm: + fit(ONE_ENDOG, nclusters, center) + with pytest.raises(ValueError, match="too few clusters") as cue: + fit(ONE_ENDOG, nclusters, center, model=IVGMMCUE) + message = str(ivgmm.value) + assert IVGMM_CONSEQUENCE in message + assert CUE_CONSEQUENCE in str(cue.value) + assert message.replace(IVGMM_CONSEQUENCE, CUE_CONSEQUENCE) == str(cue.value) + rank = nclusters - 1 if center else nclusters + assert f"Its rank is at most {rank}," in message + + +def test_ivgmm_message_with_too_few_observations_has_the_form_of_the_cue_message(): + data = small_sample(4) + with pytest.raises(ValueError, match="too few observations") as ivgmm: + IVGMM(*data, center=True).fit() + with pytest.raises(ValueError, match="too few observations") as cue: + IVGMMCUE(*data).fit() + message = str(ivgmm.value) + assert "(4) for the number of moment conditions" in message + assert message.replace(IVGMM_CONSEQUENCE, CUE_CONSEQUENCE) == str(cue.value) + + +@pytest.mark.parametrize("nclusters", [2, 3]) +def test_ivgmm_with_the_report_of_the_bug(nclusters): + # The estimate of the coefficient on x was 1.697 with 3 clusters when it is + # 0.5, the standard errors were NaN, and it was a LinAlgError on some + # platforms. There is nothing to estimate + check_too_few_clusters(ONE_ENDOG, nclusters, center=False) + + +@pytest.mark.parametrize("nclusters", [2, 3]) +@pytest.mark.parametrize("center", [False, True]) +def test_ivgmm_exactly_identified_with_too_few_clusters(nclusters, center): + # The weight matrix does not matter for the estimates of an exactly identified + # model, but the noise in it did: the estimates were 10 away from the IV + # estimates with 3 clusters, and the standard errors were not numbers. The + # moment conditions at the first step estimates of an exactly identified + # model add to zero, so that 3 clusters are not enough for its 3 moment + # conditions even if they are not centered + with pytest.raises(ValueError, match="too few clusters") as exc_info: + fit(EXACTLY_IDENTIFIED, nclusters, center) + assert exc_info.type is ValueError + assert "3 moment conditions" in str(exc_info.value) + assert f"({nclusters} clusters)" in str(exc_info.value) + + +@pytest.mark.parametrize("nclusters", [4, 5, 12]) +@pytest.mark.parametrize("center", [False, True]) +def test_ivgmm_exactly_identified_with_enough_clusters(nclusters, center): + res = fit(EXACTLY_IDENTIFIED, nclusters, center) + iv = IV2SLS(*EXACTLY_IDENTIFIED).fit() + assert_allclose(res.params, iv.params, rtol=1e-8) + assert np.all(np.isfinite(res.std_errors)) + + +# params, std_errors, J, centered J and C statistics of valid models, which are +# the same as before the models with a singular covariance were not estimated. +# nan is an invalid statistic: 4 clusters are not enough for the centered +# covariance of the 4 moments, or for the 5 moments of the C statistic +REFERENCE = { + "one_endog_4_uncentered": ( + ONE_ENDOG, + 4, + False, + {}, + None, + [1.1060294133733184, -0.09579152673627506, 0.5112972864437724], + [0.12262127936353251, 0.05938494098579946, 0.014495052398286558], + 0.7641884453729715, + np.nan, + np.nan, + ), + "one_endog_5_uncentered": ( + ONE_ENDOG, + 5, + False, + {}, + None, + [1.1686397759103972, 0.02314636167513907, 0.4759766716861976], + [0.02478733188229867, 0.10956469961382122, 0.043802173882083066], + 0.2854469335292345, + 0.2766732068378109, + 0.8759754606839069, + ), + "one_endog_5_centered": ( + ONE_ENDOG, + 5, + True, + {}, + None, + [1.170136490148785, 0.022642982069993402, 0.47318217390609696], + [0.02477994933455665, 0.10910357791629804, 0.04336862649295782], + 0.30272957956427843, + 0.31356562792598913, + np.nan, + ), + "one_endog_12_centered": ( + ONE_ENDOG, + 12, + True, + {}, + None, + [1.1416369745714467, -0.030890652538546615, 0.5142191592177071], + [0.08948919724364655, 0.0900933799632663, 0.03672622060720279], + 0.5092695078373961, + 0.511039165481266, + 0.3666652130379392, + ), + "one_endog_12_iterated": ( + ONE_ENDOG, + 12, + False, + {"iter_limit": 5}, + None, + [1.138693223505445, -0.03130105409944977, 0.514937395372159], + [0.08941806134708769, 0.09013853082731596, 0.036698020920592794], + 0.48939723718769423, + 0.5102753816307153, + 0.3276841075791217, + ), + "two_endog_6_uncentered": ( + TWO_ENDOG, + 6, + False, + {}, + ["x1"], + [0.9880803474466319, 0.4146283703294095, -0.3076604194012811], + [0.024254582559327085, 0.07723707544855464, 0.018531865035777875], + 0.7736276096934414, + 1.0570377720126551, + 0.12392028655991782, + ), + "two_endog_12_centered": ( + TWO_ENDOG, + 12, + True, + {}, + None, + [1.017400423332445, 0.3946521531745319, -0.32124858101494413], + [0.05108132199782605, 0.06125247442528457, 0.03134643143047266], + 0.4445148871191505, + 0.4406830165857778, + 5.354334781473466, + ), +} + + +def check_stat(stat, expected, name, df): + if np.isnan(expected): + check_invalid(stat, name) + else: + check_valid(stat, expected, df) + + +@pytest.mark.parametrize("name", REFERENCE) +def test_ivgmm_with_a_nonsingular_covariance_is_unchanged(name): + data, nclusters, center, kwargs, tested, params, std_errors, j, cj, c = REFERENCE[ + name + ] + res = fit(data, nclusters, center, **kwargs) + assert_allclose(res.params, params, rtol=1e-8) + assert_allclose(res.std_errors, std_errors, rtol=1e-8) + check_stat(res.j_stat, j, J_NAME, df=1) + check_stat( + res.centered_j_stat, cj, "Centered J-test of overidentifying restrictions", 1 + ) + ntested = data[2].shape[1] if tested is None else len(tested) + check_stat(res.c_stat(tested), c, "C-statistic", df=ntested) + + +@pytest.mark.parametrize("name", REFERENCE) +def test_ivgmm_with_a_nonsingular_covariance_is_the_definition(name): + data, nclusters, center, kwargs, _, params, _, j, _, _ = REFERENCE[name] + x, y, z, _, _ = arrays(data) + cov = clustered_cov(np.arange(y.shape[0]) % nclusters, center) + # Two steps unless iter_limit is given, which the definition does not do + if "iter_limit" not in kwargs: + assert_allclose(reference_params(x, y, z, cov), params, rtol=1e-8) + assert_allclose(reference_j(x, y, z, cov), j, rtol=1e-8) + + +def reference_params(x, y, z, cov): + """Two-step GMM estimates, from the definition""" + first = gmm_params(x, y, z, np.linalg.inv(z.T @ z / y.shape[0])) + w = np.linalg.inv(cov(z * (y - x @ first)[:, None])) + return gmm_params(x, y, z, w) + + +def test_ivgmm_one_step_does_not_use_the_covariance(): + # The weight matrix is (Z'Z/n)^{-1}, so there is nothing to be singular + res = fit(ONE_ENDOG, 2, center=False, iter_limit=1) + x, y, z, _, _ = arrays(ONE_ENDOG) + expected = gmm_params(x, y, z, np.linalg.inv(z.T @ z / y.shape[0])) + assert_allclose(res.params, expected, rtol=1e-8) + assert res.iterations == 1 + + +@pytest.mark.parametrize("iter_limit", [2, 5]) +def test_ivgmm_with_an_initial_weight_still_tests_the_covariance(iter_limit): + initial_weight = np.eye(4) + check_too_few_clusters( + ONE_ENDOG, 3, False, initial_weight=initial_weight, iter_limit=iter_limit + ) + + +def test_ivgmm_is_not_estimated_if_the_covariance_is_singular_in_a_later_step( + monkeypatch, +): + # The covariance is estimated at every step and not only at the first + calls = [] + original = linearmodels.iv.model._is_singular + + def second_call_is_singular(cov): + calls.append(1) + return len(calls) == 2 or original(cov) + + monkeypatch.setattr(linearmodels.iv.model, "_is_singular", second_call_is_singular) + with pytest.raises(ValueError, match="too few clusters"): + fit(ONE_ENDOG, 12, center=False, iter_limit=5, tol=0.0) + assert len(calls) == 2 + + +def test_c_stat_is_invalid_if_the_model_of_the_statistic_cannot_be_estimated(): + # It has 5 moment conditions and needs 5 clusters, whereas the model that was + # estimated has 4. The model of the statistic raises, as fit does, and the + # statistic does not + res = fit(ONE_ENDOG, 4, center=False) + c_stat = res.c_stat() + check_invalid(c_stat, "C-statistic") + assert "singular" in str(c_stat) From eea242037da137bd7ca496882cceba6ce8a95c1d Mon Sep 17 00:00:00 2001 From: Kevin Sheppard Date: Mon, 12 Oct 2026 11:32:44 +0100 Subject: [PATCH 2/5] BUG: Raise a ValueError if IVGMMCUE has as many clusters as moment conditions IVGMMCUE with weight_type="clustered" and center=False fitted a model with as many clusters as moment conditions without an error, but the objective function is constant. The cluster sums of the moment conditions are the rows of a square matrix Xi (G by K), the average of the moment conditions is Xi'1 / n and the covariance is Xi'Xi / n, so that the objective function is 1'Xi(Xi'Xi)^{-1}Xi'1, which is G whenever Xi is nonsingular, for every value of the parameters (divided by the small sample adjustment if the covariance is debiased). In the model of the bug report with 4 clusters the objective is 3.99999999997 at every parameter value that was tried, the optimizer stops after one iteration, and the "estimates" are the starting values, the estimates of IVGMM, with a J statistic equal to the number of clusters. #754 made IVGMMCUE raise if the covariance is singular, which needs more clusters than moment conditions if they are centered (the default) and at least as many if they are not. Now more clusters than moment conditions are needed with every value of center. The check is made at the same place, before the starting values are estimated and before the objective is optimized, and also if starting values are given. The message says that there are as many clusters as moment conditions, that the moment conditions are not centered, that the objective function is the same for every value of the parameters so that there is nothing to optimize, and to use more clusters than moment conditions or a weight_type that does not use clusters. With as many clusters as moment conditions a model with centered moment conditions, and an exactly identified model, are singular and raise the message of #754. Nothing else changes: of the 600 random models of the previous commit, the only one that is affected is the one with 4 clusters and 4 uncentered moment conditions, and its J statistic was 3.9999999999999654. Tests: the objective function of the model with 4 clusters is 4 at five random parameter vectors (and 4 / (the adjustment) if debiased), computed by IVGMMCUE.j and from the definition without the package's classes; fit raises before it optimizes, with and without starting values; the message; and a test of 2 to 8 clusters with both values of center, which raises for 4 or fewer clusters in IVGMMCUE, and for IVGMM exactly when the covariance is singular. The test of #754 that fitted a model with 4 uncentered clusters expects an error now. --- linearmodels/iv/model.py | 58 ++++++++++--- .../tests/iv/test_gmm_singular_covariance.py | 83 +++++++++++++++++-- 2 files changed, 123 insertions(+), 18 deletions(-) diff --git a/linearmodels/iv/model.py b/linearmodels/iv/model.py index 4186866d3f..cd1a2efbbc 100644 --- a/linearmodels/iv/model.py +++ b/linearmodels/iv/model.py @@ -1609,7 +1609,12 @@ class IVGMMCUE(_IVGMMBase): The covariance of the moment conditions has to be nonsingular. A clustered covariance that is estimated from G clusters has rank at most G, and at most G - 1 if the moment conditions are centered, so that there have to be - more clusters than moment conditions if ``center`` is ``True``. + more clusters than moment conditions if ``center`` is ``True``. If + ``center`` is ``False`` and there are as many clusters as moment + conditions, the covariance is nonsingular but the objective function is + the same for every value of the parameters, so that there is nothing to + optimize. There have to be more clusters than moment conditions whether or + not the moment conditions are centered. .. math:: @@ -1847,22 +1852,48 @@ def _callback(params: ndarray[tuple[int], dtype[float64]]) -> None: def _check_moment_covariance(self) -> None: """ - Raise a ValueError if the covariance of the moment conditions is singular - - The objective function uses the inverse of the covariance, so it is not - defined if the covariance is singular. The rank of the covariance does - not depend on the parameters, other than for special values such as - residuals that are all 0, and so it is tested at the first step - estimates of two-step GMM. + Raise a ValueError if the objective function is not defined or constant + + The objective function uses the inverse of the covariance of the moment + conditions, so it is not defined if the covariance is singular. The rank + of the covariance does not depend on the parameters, other than for + special values such as residuals that are all 0, and so it is tested at + the first step estimates of two-step GMM. + + If the covariance is clustered, the moment conditions are not centered + and there are as many clusters as moment conditions, the covariance is + not singular but the objective function is the same for all values of + the parameters. The cluster sums of the moment conditions are the rows + of a square matrix, Xi, the average of the moment conditions is + Xi'1 / n, and the covariance is Xi'Xi / n, so that the objective function + is 1'Xi(Xi'Xi)^{-1}Xi'1, which is the number of clusters whenever Xi is + nonsingular (it is divided by the small sample adjustment if debiased). """ wy, wx, wz = self._wy, self._wx, self._wz - nobs = wz.shape[0] + nobs, nmoment = wz.shape params = IVGMM.estimate_parameters(wx, wy, wz, inv(wz.T @ wz / nobs)) self._raise_if_singular( self._weight.weight_matrix(wx, wz, wy - wx @ params), "The objective function of the continuously updating GMM estimator " "uses the inverse of the covariance, so it is not defined.", ) + if not isinstance(self._weight, OneWayClusteredWeightMatrix): + return + nclusters = len(unique(asarray(self._weight_config["clusters"]))) + # Fewer clusters, or centered moment conditions, give a singular covariance + if nclusters > nmoment: + return + msg = ( + f"The clustered covariance of the {nmoment} moment conditions is " + f"estimated from as many clusters ({nclusters} clusters) as there " + "are moment conditions, and the moment conditions are not centered. " + "The objective function of the continuously updating GMM estimator " + "is then the same for every value of the parameters, so that there is " + "nothing to optimize and the estimates would be the starting values. " + "Use more clusters than moment conditions, or a weight_type that " + "does not use clusters." + ) + raise ValueError(msg) def fit( self, @@ -1910,10 +1941,11 @@ def fit( ------ ValueError If the covariance of the moment conditions is singular, which - occurs if there are too few clusters or observations. The - number of clusters has to be larger than the number of moment - conditions if the moment conditions are centered, which is the - default. + occurs if there are too few clusters or observations, or if the + moment conditions are not centered and there are as many clusters + as moment conditions, so that the objective function is constant. + The number of clusters has to be larger than the number of moment + conditions. Notes ----- diff --git a/linearmodels/tests/iv/test_gmm_singular_covariance.py b/linearmodels/tests/iv/test_gmm_singular_covariance.py index 0f186283c0..0f39e04492 100644 --- a/linearmodels/tests/iv/test_gmm_singular_covariance.py +++ b/linearmodels/tests/iv/test_gmm_singular_covariance.py @@ -27,6 +27,11 @@ and a model that cannot be estimated has no statistics to test. The C statistic estimates another model, which has more moment conditions, and is invalid if that model cannot be estimated. + +The objective function of IVGMMCUE is also constant, and not only undefined, if +the moment conditions are not centered and there are as many clusters as moment +conditions. IVGMMCUE.fit raises a ValueError then too, so that there have to be +more clusters than moment conditions with every value of center. """ import numpy as np @@ -498,8 +503,10 @@ def cue_criterion(data, params, nclusters, center): # The model has 4 moment conditions. Centered cluster sums add to zero, so that the -# covariance is singular unless there are more clusters than moment conditions, and -# at least as many are needed if the moment conditions are not centered +# covariance is singular unless there are more clusters than moment conditions. If +# the moment conditions are not centered the covariance is singular with fewer +# clusters than moment conditions, and the objective function is constant with +# as many (see test_cue_uncentered_with_as_many_clusters_as_moments) CUE_SINGULAR = [ pytest.param(2, {}, id="default-2"), pytest.param(3, {}, id="default-3"), @@ -525,6 +532,73 @@ def test_cue_too_few_clusters(nclusters, config): assert "singular" in message +@pytest.mark.parametrize("debiased", [False, True]) +def test_cue_uncentered_with_as_many_clusters_as_moments(debiased): + # The covariance is not singular, but the objective function is + # 1'Xi(Xi'Xi)^{-1}Xi'1 = 4 for all parameters, with Xi the matrix of the 4 + # cluster sums of the moment conditions, which is square. It is divided by + # the small sample adjustment of the clustered covariance if it is debiased. + # The model fitted without an error, and the estimates were the starting + # values + mod, clusters = cue_model(ONE_ENDOG, 4, center=False, debiased=debiased) + nobs, nvar = 120, 3 + scale = (nobs - 1) / (nobs - nvar) * 4 / 3 if debiased else 1.0 + rng = np.random.default_rng(0) + for params in rng.standard_normal((5, 3)) * 3 + [1.0, 0.0, 0.5]: + value = mod.j(params, mod._wx, mod._wy, mod._wz) + assert_allclose(value, 4.0 / scale, rtol=1e-9) + if not debiased: + assert_allclose(cue_criterion(ONE_ENDOG, params, 4, False), 4.0, rtol=1e-9) + with pytest.raises(ValueError, match="as many clusters") as exc_info: + mod.fit(cov_type="clustered", clusters=clusters) + assert exc_info.type is ValueError + message = str(exc_info.value) + assert "as many clusters (4 clusters) as there are moment conditions" in message + assert "4 moment conditions" in message + assert "not centered" in message + assert "the same for every value of the parameters" in message + assert "Use more clusters than moment conditions" in message + + +def test_cue_uncentered_with_as_many_clusters_as_moments_with_starting_values(): + mod, clusters = cue_model(ONE_ENDOG, 4, center=False) + starting = pd.Series([1.0, 0.0, 0.5], index=["const", "w", "x"]) + with pytest.raises(ValueError, match="as many clusters"): + mod.fit(starting=starting, cov_type="clustered", clusters=clusters) + + +def test_cue_uncentered_with_as_many_clusters_as_moments_is_found_first(monkeypatch): + def minimize(*args, **kwargs): + raise AssertionError("the objective is constant, do not optimize it") + + monkeypatch.setattr(linearmodels.iv.model, "minimize", minimize) + mod, clusters = cue_model(ONE_ENDOG, 4, center=False) + with pytest.raises(ValueError, match="as many clusters"): + mod.fit(cov_type="clustered", clusters=clusters) + + +@pytest.mark.parametrize("center", [False, True]) +@pytest.mark.parametrize("nclusters", [2, 3, 4, 5, 6, 8]) +def test_cue_and_ivgmm_clusters(nclusters, center): + # IVGMMCUE needs more clusters than the 4 moment conditions, with every + # value of center. IVGMM needs a nonsingular covariance, which is one cluster + # less if the moment conditions are not centered + if nclusters > 4: + mod, clusters = cue_model(ONE_ENDOG, nclusters, center=center) + res = mod.fit(cov_type="clustered", clusters=clusters) + assert np.all(np.isfinite(res.params)) + assert np.all(np.isfinite(res.std_errors)) + else: + mod, clusters = cue_model(ONE_ENDOG, nclusters, center=center) + with pytest.raises(ValueError, match="clusters") as exc_info: + mod.fit(cov_type="clustered", clusters=clusters) + assert exc_info.type is ValueError + if nonsingular(4, nclusters, center): + assert np.all(np.isfinite(fit(ONE_ENDOG, nclusters, center).params)) + else: + check_too_few_clusters(ONE_ENDOG, nclusters, center) + + def test_cue_too_few_clusters_with_starting_values(): mod, clusters = cue_model(ONE_ENDOG, 3) starting = pd.Series([1.0, 0.0, 0.5], index=["const", "w", "x"]) @@ -545,11 +619,10 @@ def minimize(*args, **kwargs): @pytest.mark.parametrize( ("nclusters", "center"), - [(5, True), (6, True), (12, True), (4, False), (5, False), (12, False)], + [(5, True), (6, True), (12, True), (5, False), (6, False), (12, False)], ) def test_cue_with_enough_clusters_is_not_changed(nclusters, center): - # The first three of these have more clusters than moment conditions. A model - # that is not centered can have as many + # These have more clusters than the 4 moment conditions mod, clusters = cue_model(ONE_ENDOG, nclusters, center=center) res = mod.fit(cov_type="clustered", clusters=clusters) assert np.all(np.isfinite(res.params)) From f88da073bab483c264a1d2f0336f89d6f2ea7522 Mon Sep 17 00:00:00 2001 From: Kevin Sheppard Date: Mon, 12 Oct 2026 11:36:27 +0100 Subject: [PATCH 3/5] DOC: Say that IVGMM and IVGMMCUE raise before j_stat can be invalid j_stat said that an InvalidTestStatistic is returned if the covariance of the moment conditions that is inverted is singular, which can happen with too few clusters. IVGMM.fit and IVGMMCUE.fit raise a ValueError in that case instead. The statistic can only be invalid now at the final estimates of IVGMMCUE, if the covariance is nonsingular at the first step estimates, where it is tested, and singular at the final estimates. --- linearmodels/iv/results.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/linearmodels/iv/results.py b/linearmodels/iv/results.py index a11c936933..2f78f5b13f 100644 --- a/linearmodels/iv/results.py +++ b/linearmodels/iv/results.py @@ -1671,7 +1671,9 @@ def j_stat(self) -> InvalidTestStatistic | WaldTestStatistic: J statistic test of overidentifying restrictions. An InvalidTestStatistic is returned if the estimate of the covariance of the moment conditions that was inverted to compute - the weight matrix is singular, which can happen if there are too + the weight matrix is singular. This can only happen at the final + estimates of IVGMMCUE, since IVGMM and IVGMMCUE raise a + ValueError if it is singular, which is the case if there are too few clusters. Notes @@ -1700,6 +1702,13 @@ def j_stat(self) -> InvalidTestStatistic | WaldTestStatistic: same way by the number of observations. Whether the estimate is singular does not depend on the units in which the variables are measured. + + ``IVGMM.fit`` raises a ``ValueError`` if the estimate is singular in + any step, and ``IVGMMCUE.fit`` raises one if it is singular at the + first step estimates, since the weight matrix and the objective + function are not defined. The statistic is invalid only if the + estimate of ``IVGMMCUE`` is singular at its final estimates and not at + the first step estimates. """ return self._j_stat From fe6b02c68f39a6ec454d7c145d1e10f8995cfc56 Mon Sep 17 00:00:00 2001 From: Kevin Sheppard Date: Mon, 12 Oct 2026 11:46:11 +0100 Subject: [PATCH 4/5] TST: Test the centered J with as many clusters as moment conditions using IVGMM test_centered_j_cue_with_too_few_clusters fitted an IVGMMCUE with center=False and 5 clusters for 5 moment conditions to check that the centered J statistic is invalid when the uncentered fit is valid. That model cannot be estimated any more: the objective function of the CUE is constant with as many clusters as moment conditions, and the CUE needs 6 clusters, for which the centered covariance is not singular either. The two-step estimator can still be estimated with 5 uncentered clusters, and it is the case where the J statistic is valid and the centered J statistic is not, so the test checks it with IVGMM. It checks that the CUE raises with 5 clusters, and that the CUE with 6 clusters has a valid centered J statistic. The test of the centered J statistic of the CUE with a singular covariance of #754 is unchanged. --- .../tests/iv/test_new_feature_hardening.py | 29 ++++++++++++------- 1 file changed, 18 insertions(+), 11 deletions(-) diff --git a/linearmodels/tests/iv/test_new_feature_hardening.py b/linearmodels/tests/iv/test_new_feature_hardening.py index d6c5bd9672..5da432c3a3 100644 --- a/linearmodels/tests/iv/test_new_feature_hardening.py +++ b/linearmodels/tests/iv/test_new_feature_hardening.py @@ -657,24 +657,31 @@ def test_centered_j_cue_without_centering(weight_type, weighted): def test_centered_j_cue_with_too_few_clusters(weighted): # The centered cluster sums add to zero, so that the covariance has at most # nclusters - 1 non-zero eigenvalues, and the model has 5 moment conditions. - # The uncentered covariance of the CUE has nclusters + # The uncentered covariance has nclusters. The objective function of the CUE + # is constant if there are as many clusters as moment conditions, so that the + # CUE needs 6 clusters, which are enough for the centered covariance too, and + # the case where only the uncentered covariance can be estimated is the + # two-step estimator y, exog, endog, z, w = simulate(23, ninstr=3, violation=0.3, weighted=weighted) outcomes = {} for nclusters in (5, 6): - mod = IVGMMCUE( - y, - exog, - endog, - z, - weights=w, - weight_type="clustered", - clusters=np.arange(y.shape[0]) % nclusters, - center=False, - ) + config = { + "weights": w, + "weight_type": "clustered", + "clusters": np.arange(y.shape[0]) % nclusters, + "center": False, + } + mod = IVGMM(y, exog, endog, z, **config) res = mod.fit() outcomes[nclusters] = res assert np.isfinite(res.j_stat.stat) assert res.j_stat.df == 2 + cue = IVGMMCUE(y, exog, endog, z, **config) + if nclusters == 5: + with pytest.raises(ValueError, match="as many clusters"): + cue.fit() + else: + assert type(cue.fit().centered_j_stat) is WaldTestStatistic invalid = outcomes[5].centered_j_stat assert isinstance(invalid, InvalidTestStatistic) assert np.isnan(invalid.stat) From 936799227f957a67d1a314549d89b7d12a280c25 Mon Sep 17 00:00:00 2001 From: Kevin Sheppard Date: Mon, 12 Oct 2026 11:46:12 +0100 Subject: [PATCH 5/5] DOC: Say how many clusters IVGMM needs The class and fit documentation of IVGMM said that there have to be more clusters than moment conditions. A clustered estimate from G clusters has rank at most G, and at most G - 1 if the moment conditions are centered, so that an overidentified model with center=False can be estimated with as many clusters as moment conditions, and one with center=True cannot. The first step moment conditions of an exactly identified model add to zero, so that it needs more clusters than moment conditions with every value of center. --- linearmodels/iv/model.py | 22 +++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/linearmodels/iv/model.py b/linearmodels/iv/model.py index cd1a2efbbc..a29dda3258 100644 --- a/linearmodels/iv/model.py +++ b/linearmodels/iv/model.py @@ -1307,12 +1307,15 @@ class IVGMM(_IVGMMBase): moment conditions, so that the estimate has to be nonsingular. A clustered estimate from G clusters has rank at most G, and at most G - 1 if the moment conditions are centered (``center`` is ``True``), since the - centered cluster sums add to zero. There have to be more clusters than - moment conditions. The rank of a kernel or robust estimate is bounded in - the same way by the number of observations. ``fit`` raises a - ``ValueError`` that says so if the estimate is singular, since the - estimates would be numerical noise. The estimate is not used, and so it - is not tested, if ``iter_limit`` is 1. + centered cluster sums add to zero. There have to be at least as many + clusters as moment conditions, and more if they are centered. The first + step moment conditions of an exactly identified model add to zero as well, + so that it needs more clusters than moment conditions with every value of + ``center``. The rank of a kernel or robust estimate is bounded in the same + way by the number of observations. ``fit`` raises a ``ValueError`` that + says so if the estimate is singular, since the estimates would be + numerical noise. The estimate is not used, and so it is not tested, if + ``iter_limit`` is 1. .. todo:: @@ -1493,9 +1496,10 @@ def fit( ValueError If the estimate of the covariance of the moment conditions is singular, which occurs if there are too few clusters or - observations. The number of clusters has to be larger than the - number of moment conditions. The estimate is not used if - ``iter_limit`` is 1. + observations. A clustered estimate needs at least as many clusters + as moment conditions, and more if the moment conditions are + centered or the model is exactly identified. The estimate is not + used if ``iter_limit`` is 1. See also --------