From c1a549ff87806dcb3a5e029b95cf29ea2f5dc9c9 Mon Sep 17 00:00:00 2001 From: Kevin Sheppard Date: Sun, 11 Oct 2026 22:47:16 +0100 Subject: [PATCH 1/3] BUG: Return an invalid J and C statistic if the moment covariance is singular Stacked on lm71-bug3b-iv-features (BUG3b), which has to be merged first. That branch adds the scale-free rank check of the covariance of the moment conditions to IVGMMResults.centered_j_stat. This commit moves the check to a function and also applies it to j_stat and c_stat. IVGMMResults.c_stat and IVGMMResults.j_stat gave meaningless numbers if the estimate of the covariance of the moment conditions was singular, which is the case if there are too few clusters. The weight matrix is the inverse of that estimate, so it is numerical noise, and so are the statistics. c_stat uses the weight type of the model since #744, and so it uses the clustered covariance. Reproduction: n = 120, a constant and an exogenous regressor, one endogenous regressor and two instruments, which are 4 moment conditions, and clusters = arange(n) % G. The model is IVGMM(..., weight_type="clustered", clusters=clusters).fit(cov_type="clustered", clusters=clusters). G c_stat() before c_stat() now 3 -3.73e15 invalid 4 -7.42e16 invalid 5 0.8760 0.8760 The model of the C statistic treats the endogenous regressor as exogenous, so it has 5 moment conditions, and its covariance is singular unless there are at least 5 clusters. 7.0 gave -0.049 for every G, since it ignored the clusters, which is not valid either. c_stat also raised a LinAlgError when NumPy found that the covariance is exactly singular, e.g., with 2 centered clusters. j_stat is affected in the same way, in 7.0 and in main. J of the same model, for the default of IVGMM (center=False) and for center=True, which is the default of IVGMMCUE: G center=False center=True now 2 -3.19e15 1.41 invalid, invalid 3 -6.87e15 1.73e14 invalid, invalid 4 0.764 1.52e15 valid, invalid 5 0.285 0.303 valid, valid A covariance that is estimated from G clusters is a sum of G outer products, so its rank is at most G. The centered cluster sums add to zero, so the rank is at most G - 1 if the moment conditions are centered. With center=True and G = 2 the statistic looks plausible, 1.41, but it is as meaningless as the others. The rank of a kernel or a robust covariance is bounded by the number of observations in the same way, and so these have singular cases too, e.g., 4 observations and 4 moment conditions with center=True gave J = 4.2e6 with a kernel and 6.96e9 with the robust estimator. The fix: * The check of centered_j_stat is a function, _is_singular. Its rank is that of the correlation matrix, so it does not depend on the units of the data. * IVGMM.fit and IVGMMCUE.fit pass the covariance that they inverted to the post-estimation, and j_stat is an InvalidTestStatistic if it is singular. It is not checked if no covariance is estimated, which is the case if iter_limit=1, since the weight matrix is then the inverse of Z'Z / n. * c_stat is an InvalidTestStatistic with a reason if the J statistic of the model that treats the tested variables as exogenous is invalid, or if fitting it raises a LinAlgError. The C statistic is the difference of two J statistics that use the covariance of this model. The covariance is tested and not its inverse. In 48,000 simulated clustered covariances, with the units of the instruments between 1e-3 and 1e3, 14,809 of which are singular by construction, the test of the covariance found all of them and had no false alarms. The same test of the inverse missed 780 of them (5%). The estimates are not changed. A GMM model with a singular covariance still uses the inverse of noise, so its parameters and standard errors are not meaningful, or it raises a LinAlgError. The statistics of every model that is not singular are unchanged: the reference tests of the J, centered J and C statistics are not modified. The new tests compute the expected statistics from the definitions, without the weight matrix classes, for the cases above and for models with two endogenous regressors and tested subsets, weights, continuously updating GMM, kernel and robust covariances with 4 to 20 observations, and instruments in units that differ by 1e16. They also check that the statistics are valid with exactly as many clusters as moment conditions, that j_stat is not invalid if no covariance is estimated, and the LinAlgError of c_stat. 54 of the 118 tests fail without the change. --- linearmodels/iv/model.py | 82 +++- linearmodels/iv/results.py | 57 ++- .../tests/iv/test_gmm_singular_covariance.py | 456 ++++++++++++++++++ 3 files changed, 569 insertions(+), 26 deletions(-) create mode 100644 linearmodels/tests/iv/test_gmm_singular_covariance.py diff --git a/linearmodels/iv/model.py b/linearmodels/iv/model.py index 4eecedeab3..a6510b7c9d 100644 --- a/linearmodels/iv/model.py +++ b/linearmodels/iv/model.py @@ -1005,6 +1005,36 @@ def from_formula( return mod +def _is_singular(cov: linearmodels.typing.data.Float64Array) -> bool: + """ + Test whether the covariance of the moment conditions is singular + + Parameters + ---------- + cov : ndarray + Estimate of the covariance of the moment conditions (ninstr by ninstr) + + Returns + ------- + bool + True if the covariance is not finite, has a variance that is not + positive or its correlation matrix does not have full rank. + + Notes + ----- + The rank is that of the correlation matrix of the moment conditions, so + that it does not depend on the units in which the instruments are + measured. The tolerance of the rank test is relative to the largest + eigenvalue, and so it treats the covariance as singular if the variances + of the moment conditions differ by a factor of about 1e16. + """ + variances = diag(cov) + if not (npall(isfinite(cov)) and npall(variances > 0)): + return True + std = sqrt(variances) + return bool(matrix_rank(cov / outer(std, std), hermitian=True) < cov.shape[0]) + + class _IVGMMBase(_IVModelBase): r""" Estimation of IV models using the generalized method of moments (GMM) @@ -1084,6 +1114,7 @@ def _gmm_post_estimation( params: linearmodels.typing.data.Float64Array, weight_mat: linearmodels.typing.data.Float64Array, iters: int, + moment_cov: linearmodels.typing.data.Float64Array | None = None, ) -> dict[str, Any]: """GMM-specific post-estimation results""" instr = self._instr_columns @@ -1092,7 +1123,7 @@ def _gmm_post_estimation( "weight_type": self._weight_type, "weight_config": self._weight.config, "iterations": iters, - "j_stat": self._j_statistic(params, weight_mat), + "j_stat": self._j_statistic(params, weight_mat, moment_cov), "centered_j_stat": self._centered_j_statistic(params), } @@ -1118,17 +1149,7 @@ def _centered_j_statistic( weight_matrix_estimator = WEIGHT_MATRICES[self._weight_type] centered = weight_matrix_estimator(**{**self._weight_kwargs, "center": True}) s_c = centered.weight_matrix(x, z, eps) - # The rank is that of the correlation matrix of the moment conditions, - # so that it does not depend on the units in which the instruments are - # measured. The tolerance of the rank test is relative to the largest - # eigenvalue, and so it treats the covariance as singular if the - # variances of the moment conditions differ by a factor of about 1e16 - variances = diag(s_c) - singular = not (npall(isfinite(s_c)) and npall(variances > 0)) - if not singular: - std = sqrt(variances) - singular = matrix_rank(s_c / outer(std, std), hermitian=True) < ninstr - if singular: + if _is_singular(s_c): # The centered cluster sums add to zero, so a clustered estimator # has rank at most one less than the number of clusters reason = ( @@ -1145,10 +1166,29 @@ def _j_statistic( self, params: linearmodels.typing.data.Float64Array, weight_mat: linearmodels.typing.data.Float64Array, - ) -> WaldTestStatistic: - """J stat and test""" + moment_cov: linearmodels.typing.data.Float64Array | None = None, + ) -> WaldTestStatistic | InvalidTestStatistic: + """ + J stat and test + + moment_cov is the estimate of the covariance of the moment conditions + that was inverted to compute weight_mat, if there is one. The statistic + is not defined if it is singular, since weight_mat is then numerical + noise. The test is applied to the covariance and not to weight_mat, + because the inverse of a singular covariance is not always found to be + singular. + """ y, x, z = self._wy, self._wx, self._wz nobs, nvar, ninstr = y.shape[0], x.shape[1], z.shape[1] + if moment_cov is not None and _is_singular(moment_cov): + reason = ( + "The covariance of the moment conditions that is inverted to " + "compute the weight matrix is singular, so the statistic is " + "not defined. This occurs, for example, when there are too " + "few clusters." + ) + name = "J-test of overidentifying restrictions" + return InvalidTestStatistic(reason, name=name) eps = y - x @ params g_bar = (z * eps).mean(0) stat = float(nobs * g_bar.T @ weight_mat @ g_bar.T) @@ -1403,9 +1443,11 @@ def fit( iters, norm = 1, 10 * tol + 1 vinv = eye(params.shape[0]) + moment_cov: linearmodels.typing.data.Float64Array | None = None while iters < iter_limit and norm > tol: eps = wy - wx @ params - wmat = inv(weight_matrix(wx, wz, eps)) + moment_cov = weight_matrix(wx, wz, eps) + wmat = inv(moment_cov) params = self.estimate_parameters(wx, wy, wz, wmat) delta = params - _params if iters == 1: @@ -1422,7 +1464,7 @@ def fit( ) results = self._post_estimation(params, cov_estimator, cov_type) - gmm_pe = self._gmm_post_estimation(params, wmat, iters) + gmm_pe = self._gmm_post_estimation(params, wmat, iters, moment_cov) results.update(gmm_pe) @@ -1433,6 +1475,7 @@ def _gmm_post_estimation( params: linearmodels.typing.data.Float64Array, weight_mat: linearmodels.typing.data.Float64Array, iters: int, + moment_cov: linearmodels.typing.data.Float64Array | None = None, ) -> dict[str, Any]: """GMM-specific post-estimation results""" instr = self._instr_columns @@ -1441,7 +1484,7 @@ def _gmm_post_estimation( "weight_type": self._weight_type, "weight_config": self._weight.config, "iterations": iters, - "j_stat": self._j_statistic(params, weight_mat), + "j_stat": self._j_statistic(params, weight_mat, moment_cov), "centered_j_stat": self._centered_j_statistic(params), } @@ -1795,14 +1838,15 @@ def fit( starting, wx, wy, wz, display, opt_options=opt_options ) eps = wy - wx @ params - wmat = inv(weight_matrix(wx, wz, eps)) + moment_cov = weight_matrix(wx, wz, eps) + wmat = inv(moment_cov) cov_config["debiased"] = debiased cov_estimator = IVGMMCovariance( wx, wy, wz, params, wmat, cov_type, **cov_config ) results = self._post_estimation(params, cov_estimator, cov_type) - gmm_pe = self._gmm_post_estimation(params, wmat, iters) + gmm_pe = self._gmm_post_estimation(params, wmat, iters, moment_cov) results.update(gmm_pe) return IVGMMResults(results, self) diff --git a/linearmodels/iv/results.py b/linearmodels/iv/results.py index ac2c848151..cd0fda3123 100644 --- a/linearmodels/iv/results.py +++ b/linearmodels/iv/results.py @@ -24,7 +24,7 @@ sqrt, squeeze, ) -from numpy.linalg import inv +from numpy.linalg import LinAlgError, inv from pandas import DataFrame, Series, concat, to_numeric from scipy import stats from statsmodels.iolib.summary import SimpleTable, fmt_2cols, fmt_params @@ -1667,8 +1667,12 @@ def j_stat(self) -> InvalidTestStatistic | WaldTestStatistic: Returns ------- - WaldTestStatistic - J statistic test of overidentifying restrictions + WaldTestStatistic or InvalidTestStatistic + 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 + few clusters. Notes ----- @@ -1685,6 +1689,17 @@ def j_stat(self) -> InvalidTestStatistic | WaldTestStatistic: :math:`\hat{\epsilon}_i` are the model residuals. :math:`W` is a consistent estimator of the variance of :math:`\sqrt{n}\bar{g}`. The degree of freedom is :math:`q = n_{instr} - n_{endog}`. + + The statistic is not defined if the estimate of the covariance of the + moment conditions that is inverted to compute the weight matrix is + singular. A clustered estimate that uses G clusters has rank at most + G, and at most G - 1 if the moment conditions are centered (``center`` + is True), because the 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 rank of a kernel or robust estimate is bounded in the + same way by the number of observations. Whether the estimate is + singular does not depend on the units in which the variables are + measured. """ return self._j_stat @@ -1763,7 +1778,9 @@ def centered_j_stat(self) -> InvalidTestStatistic | WaldTestStatistic: """ return self._centered_j_stat - def c_stat(self, variables: list[str] | str | None = None) -> WaldTestStatistic: + def c_stat( + self, variables: list[str] | str | None = None + ) -> InvalidTestStatistic | WaldTestStatistic: r""" C-test of endogeneity @@ -1775,8 +1792,11 @@ def c_stat(self, variables: list[str] | str | None = None) -> WaldTestStatistic: Returns ------- - WaldTestStatistic - Object containing test statistic, p-value, distribution and null + WaldTestStatistic or InvalidTestStatistic + Object containing test statistic, p-value, distribution and null. + An InvalidTestStatistic is returned if the covariance of the + moment conditions of the first model is singular, which can happen + if there are too few clusters. Notes ----- @@ -1821,6 +1841,16 @@ def c_stat(self, variables: list[str] | str | None = None) -> WaldTestStatistic: kernel, so that the statistic is valid under the same assumptions about the moment conditions. + The first model has a moment condition for each tested variable that + the model does not have. The estimate of the covariance of its moment + conditions has to be nonsingular, and the statistic is invalid if it is + not. A clustered estimate that uses G clusters has rank at most G, and + at most G - 1 if the moment conditions are centered (``center`` is + True), because the centered cluster sums add to zero. There have to be + at least as many clusters as moment conditions in the first model, and + more if they are centered. The rank of a kernel or robust estimate is + bounded in the same way by the number of observations. + See Hayashi (2000), pp. 218-221 and 232-234, or Baum, Schaffer and Stillman (2003), section 4.4, for details. """ @@ -1856,8 +1886,21 @@ def c_stat(self, variables: list[str] | str | None = None) -> WaldTestStatistic: weight_type=self.model._weight_type, **self.model._weight_kwargs, ) - res_e = mod.fit(cov_type=self.cov_type, **self.cov_config) + reason = ( + "The covariance of the moment conditions of the model that treats " + "the tested variables as exogenous is singular, so the statistic " + "is not defined. This occurs, for example, when there are too few " + "clusters, since that model has a moment condition for each tested " + "variable that this model does not." + ) + try: + res_e = mod.fit(cov_type=self.cov_type, **self.cov_config) + except LinAlgError: + # The covariance of the moment conditions is exactly singular + return InvalidTestStatistic(reason, name="C-statistic") assert isinstance(res_e, IVGMMResults) + if isinstance(res_e.j_stat, InvalidTestStatistic): + return InvalidTestStatistic(reason, name="C-statistic") j_e = res_e.j_stat.stat x = self.model._wx diff --git a/linearmodels/tests/iv/test_gmm_singular_covariance.py b/linearmodels/tests/iv/test_gmm_singular_covariance.py new file mode 100644 index 0000000000..0ed9da1b9e --- /dev/null +++ b/linearmodels/tests/iv/test_gmm_singular_covariance.py @@ -0,0 +1,456 @@ +""" +Tests of the J and C statistics of IVGMM and IVGMMCUE with a singular covariance + +The weight matrix of a GMM estimator is the inverse of an estimate of the +covariance of the moment conditions. That estimate is singular if there are too +few clusters. The covariance of K moment conditions that is estimated from G +clusters is a sum of G outer products of K-vectors, and so it has rank at most +G. The cluster sums of centered moment conditions add to zero, and so its rank +is at most G - 1 if the moment conditions are centered. The covariance of a +kernel or robust estimator has rank at most the number of observations, and one +less if the moment conditions are centered. + +The J statistic and the C statistic are functions of the inverse. When the +covariance is singular the inverse is numerical noise, and the statistics are +very large, negative or arbitrary numbers. They have to be invalid, as the +centered J statistic is. All of the expected values below are computed from the +definitions, without any of the weight matrix classes. + +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. +""" + +import numpy as np +from numpy.linalg import LinAlgError +from numpy.testing import assert_allclose +import pandas as pd +import pytest +from scipy import stats + +from linearmodels.iv import IVGMM, IVGMMCUE +import linearmodels.iv.model +from linearmodels.shared.hypotheses import InvalidTestStatistic, WaldTestStatistic + +# ---------------------------------------------------------------------------- +# Data +# ---------------------------------------------------------------------------- + + +def one_endog_data(nobs=120): + """2 exog, 1 endog and 2 instruments (4 moments), the model of the bug report""" + rng = np.random.default_rng(3) + z = pd.DataFrame(rng.standard_normal((nobs, 2)), columns=["z0", "z1"]) + 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 + + +def two_endog_data(nobs=120): + """1 exog, 2 endog and 3 instruments (4 moments), so that variables can be tested""" + rng = np.random.default_rng(5) + z = pd.DataFrame(rng.standard_normal((nobs, 3)), columns=["z0", "z1", "z2"]) + x = pd.DataFrame( + { + "x1": z.z0 + z.z1 + rng.standard_normal(nobs), + "x2": z.z1 - z.z2 + rng.standard_normal(nobs), + } + ) + exog = pd.DataFrame({"const": np.ones(nobs)}) + y = pd.Series(1 + 0.5 * x.x1 - 0.3 * x.x2 + rng.standard_normal(nobs), name="y") + return y, exog, x, z + + +def small_sample(nobs): + """A constant, 1 endogenous variable and 3 instruments: 4 moments""" + rng = np.random.default_rng(11) + z = pd.DataFrame(rng.standard_normal((nobs, 3)), columns=["z0", "z1", "z2"]) + x = pd.DataFrame({"x": z.sum(axis=1) + rng.standard_normal(nobs)}) + exog = pd.DataFrame({"const": np.ones(nobs)}) + y = pd.Series(1 + 0.5 * x.x + rng.standard_normal(nobs), name="y") + return y, exog, x, z + + +ONE_ENDOG = one_endog_data() +TWO_ENDOG = two_endog_data() + + +# ---------------------------------------------------------------------------- +# Definitions +# ---------------------------------------------------------------------------- + + +def gmm_params(x, y, z, w): + """GMM estimator that minimizes the quadratic form in the moments""" + xpz = x.T @ z + return np.linalg.solve(xpz @ w @ xpz.T, xpz @ w @ (z.T @ y)) + + +def j_value(x, y, z, w, params): + gbar = (z * (y - x @ params)[:, None]).mean(0) + return y.shape[0] * gbar @ w @ gbar + + +def clustered_cov(clusters, center): + """Sum of the outer products of the cluster sums of the moment conditions""" + + def cov(g): + if center: + g = g - g.mean(0) + sums = np.array([g[clusters == c].sum(0) for c in np.unique(clusters)]) + return sums.T @ sums / g.shape[0] + + return cov + + +def bartlett_cov(bandwidth, center): + """Heteroskedasticity and autocorrelation robust, with a Bartlett kernel""" + + def cov(g): + if center: + g = g - g.mean(0) + s = g.T @ g + for lag in range(1, bandwidth + 1): + op = g[lag:].T @ g[:-lag] + s += (1 - lag / (bandwidth + 1)) * (op + op.T) + return s / g.shape[0] + + return cov + + +def reference_j(x, y, z, cov): + """J statistic from two step GMM, using the covariance of the first step""" + 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 j_value(x, y, z, w, gmm_params(x, y, z, w)) + + +def reference_c(x, y, z_e, original, cov): + """ + C statistic, J_e - J_c + + z_e contains the moments of the model that treats the tested variables as + exogenous, and original are the columns of the moments of the model. + """ + first = gmm_params(x, y, z_e, np.linalg.inv(z_e.T @ z_e / y.shape[0])) + s_e = cov(z_e * (y - x @ first)[:, None]) + w_e = np.linalg.inv(s_e) + j_e = j_value(x, y, z_e, w_e, gmm_params(x, y, z_e, w_e)) + z = z_e[:, original] + w_c = np.linalg.inv(s_e[np.ix_(original, original)]) + return j_e - j_value(x, y, z, w_c, gmm_params(x, y, z, w_c)) + + +def arrays(data, tested=None, sqrt_weights=None): + """x, y, the model's moments and the moments with the tested variables added""" + y, exog, endog, instr = (np.asarray(v, dtype=float) for v in data) + y = y.reshape(-1) + names = list(data[2].columns) + tested = names if tested is None else tested + index = [names.index(name) for name in tested] + x = np.column_stack([exog, endog]) + z = np.column_stack([exog, instr]) + z_e = np.column_stack([exog, endog[:, index], instr]) + ninstr = instr.shape[1] + original = [*range(exog.shape[1]), *range(z_e.shape[1] - ninstr, z_e.shape[1])] + if sqrt_weights is not None: + x, z, z_e = (a * sqrt_weights[:, None] for a in (x, z, z_e)) + y = y * sqrt_weights + return x, y, z, z_e, original + + +def nonsingular(nmoments, nclusters, center): + """Rank condition of the covariance that is estimated from clusters""" + return nclusters - center >= nmoments + + +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. + """ + y, exog, endog, instr = data + clusters = np.arange(y.shape[0]) % nclusters + mod = model( + y, + exog, + endog, + instr, + weights=weights, + weight_type="clustered", + 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 + + +def check_invalid(stat, name): + assert type(stat) is InvalidTestStatistic + assert np.isnan(stat.stat) + assert np.isnan(stat.pval) + assert stat.critical_values is None + text = str(stat) + assert text.startswith("Invalid test statistic") + assert "singular" in text + assert name in text + + +def check_valid(stat, expected, df): + assert type(stat) is WaldTestStatistic + assert stat.df == df + assert np.isfinite(stat.stat) + assert_allclose(stat.stat, expected, rtol=1e-7, atol=1e-10) + assert_allclose(stat.pval, stats.chi2.sf(expected, df), rtol=1e-6, atol=1e-12) + + +J_NAME = "J-test of overidentifying restrictions" + +# ---------------------------------------------------------------------------- +# The model that the problem was reported with +# ---------------------------------------------------------------------------- + + +@pytest.mark.parametrize("nclusters", [3, 4]) +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) + check_invalid(res.c_stat(), "C-statistic") + + +def test_c_stat_with_enough_clusters_is_unchanged(): + res = fit(ONE_ENDOG, 5, center=False) + x, y, _, z_e, original = arrays(ONE_ENDOG) + cov = clustered_cov(np.arange(120) % 5, center=False) + check_valid(res.c_stat(), reference_c(x, y, z_e, original, cov), df=1) + # The value that the statistic has always had when it was valid + assert_allclose(res.c_stat().stat, 0.8759754606839069, rtol=1e-6) + + +# ---------------------------------------------------------------------------- +# The C statistic +# ---------------------------------------------------------------------------- + +C_STAT_MODELS = [ + pytest.param(ONE_ENDOG, None, id="one-endog"), + pytest.param(TWO_ENDOG, None, id="two-endog-all"), + pytest.param(TWO_ENDOG, ["x1"], id="two-endog-x1"), + pytest.param(TWO_ENDOG, ["x2"], id="two-endog-x2"), +] + + +@pytest.mark.parametrize("center", [False, True]) +@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: + return + c_stat = res.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): + check_invalid(c_stat, "C-statistic") + return + x, y, _, z_e, original = arrays(data, tested) + cov = clustered_cov(np.arange(y.shape[0]) % nclusters, center) + check_valid(c_stat, reference_c(x, y, z_e, original, cov), df=ntested) + assert c_stat.stat >= 0 + + +@pytest.mark.parametrize("nclusters", [5, 6]) +def test_c_stat_continuously_updating(nclusters): + # The covariance is centered unless center=False, so the model with 4 moments + # can be estimated with 5 clusters but the C statistic needs 6 + res = fit(ONE_ENDOG, nclusters, True, model=IVGMMCUE) + assert type(res.j_stat) is WaldTestStatistic + assert np.isfinite(res.j_stat.stat) + c_stat = res.c_stat() + if nclusters == 5: + check_invalid(c_stat, "C-statistic") + return + x, y, _, z_e, original = arrays(ONE_ENDOG) + cov = clustered_cov(np.arange(120) % nclusters, center=True) + check_valid(c_stat, reference_c(x, y, z_e, original, cov), df=1) + + +def test_c_stat_when_inverting_the_covariance_raises(monkeypatch): + # NumPy raises a LinAlgError if a matrix that it inverts is exactly singular, + # which happened in c_stat for the covariance of 2 centered clusters (it has + # rank 1) and was not caught. It depends on rounding, so raise it here, for + # the matrices of the model that treats the endogenous variable as exogenous, + # which has 5 moments and the model has 4 + res = fit(ONE_ENDOG, 6, center=False) + inv = linearmodels.iv.model.inv + + def raising_inv(a): + if a.shape[0] == 5: + raise LinAlgError("Singular matrix") + return inv(a) + + with monkeypatch.context() as patch: + patch.setattr(linearmodels.iv.model, "inv", raising_inv) + check_invalid(res.c_stat(), "C-statistic") + assert type(res.c_stat()) is WaldTestStatistic + + +# ---------------------------------------------------------------------------- +# The J statistic +# ---------------------------------------------------------------------------- + + +@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) + return + 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) + assert j_stat.stat >= 0 + + +@pytest.mark.parametrize("center", [False, True]) +@pytest.mark.parametrize("nclusters", [2, 3, 4, 5]) +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: + return + assert (type(res.j_stat) is WaldTestStatistic) == nonsingular(4, nclusters, center) + assert (type(res.centered_j_stat) is WaldTestStatistic) == (nclusters >= 5) + + +def test_j_stat_without_a_covariance_estimate_is_unchanged(): + # With one step the weight matrix is the inverse of the second moment of the + # instruments, which is not a covariance estimate that depends on the clusters + res = fit(ONE_ENDOG, 2, center=False, iter_limit=1) + x, y, z, _, _ = arrays(ONE_ENDOG) + w = np.linalg.inv(z.T @ z / 120) + expected = j_value(x, y, z, w, gmm_params(x, y, z, w)) + check_valid(res.j_stat, expected, df=1) + + +@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) + 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) + assert "Parameter Estimates" in str(res.summary) + assert np.all(np.isfinite(res.params)) + assert res.weight_type == "clustered" + + +# ---------------------------------------------------------------------------- +# Kernel and robust covariances, which are singular if there are too few observations +# ---------------------------------------------------------------------------- + + +@pytest.mark.parametrize("weight_type", ["kernel", "robust"]) +@pytest.mark.parametrize("nobs", [4, 5, 6, 8, 20]) +def test_small_samples_j_stat(weight_type, nobs): + # The centered moment conditions of nobs observations have rank at most + # 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) + return + 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)) + check_valid(j_stat, expected, df=2) + + +@pytest.mark.parametrize("center", [False, True]) +@pytest.mark.parametrize("weight_type", ["kernel", "robust"]) +@pytest.mark.parametrize("nobs", [5, 6, 8, 20]) +def test_small_samples_c_stat(weight_type, nobs, center): + # The C statistic needs 5 moments, so that nobs = 5 gives a singular + # covariance if the moments are centered + data = small_sample(nobs) + c_stat = IVGMM(*data, weight_type=weight_type, center=center).fit().c_stat() + if nobs == 5 and center: + check_invalid(c_stat, "C-statistic") + return + x, y, _, z_e, original = arrays(data) + bandwidth = nobs - 2 if weight_type == "kernel" else 0 + expected = reference_c(x, y, z_e, original, bartlett_cov(bandwidth, center)) + check_valid(c_stat, expected, df=1) + + +# ---------------------------------------------------------------------------- +# Weights and the units of the data +# ---------------------------------------------------------------------------- + + +@pytest.mark.parametrize("nclusters", [3, 4, 5, 6, 10]) +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: + return + 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) + if nclusters >= 5: + check_valid(res.c_stat(), reference_c(x, y, z_e, original, cov), df=1) + else: + check_invalid(res.c_stat(), "C-statistic") + + +@pytest.mark.parametrize("scale", [1e-8, 1e-4, 1e4, 1e8]) +@pytest.mark.parametrize("nclusters", [4, 12]) +def test_units_of_the_instruments_do_not_matter(nclusters, scale): + # The statistics do not depend on the units of the instruments. A covariance + # with variances that differ by 1e16 must not be mistaken for a singular one, + # which is why the rank is that of the correlation matrix + y, exog, endog, instr = ONE_ENDOG + scaled = (y, exog, endog, instr * [scale, 1.0]) + res = fit(ONE_ENDOG, nclusters, center=False) + res_scaled = fit(scaled, nclusters, center=False) + assert type(res_scaled.j_stat) is WaldTestStatistic + assert_allclose(res_scaled.j_stat.stat, res.j_stat.stat, rtol=1e-5) + if nclusters == 12: + assert type(res_scaled.c_stat()) is WaldTestStatistic + assert_allclose(res_scaled.c_stat().stat, res.c_stat().stat, rtol=1e-5) + else: + check_invalid(res_scaled.c_stat(), "C-statistic") From ef9ab91fdc1716891aa9515618aad0713c97fe9f Mon Sep 17 00:00:00 2001 From: Kevin Sheppard Date: Sun, 11 Oct 2026 22:47:17 +0100 Subject: [PATCH 2/3] BUG: Raise a ValueError that names the cause if IVGMMCUE has too few clusters IVGMMCUE with weight_type="clustered" raised numpy.linalg.LinAlgError: Singular matrix, from inside the objective function that scipy minimizes, if there were not more clusters than moment conditions. The message does not say what the problem is or what to do. It happens since #746, which makes the moment conditions centered by default, as the documentation says they should be. 7.0 fit a model with as many clusters as moment conditions, since it did not center them. Reproduction (reg1/repro_leads.py): n = 120, a constant and an exogenous regressor, one endogenous regressor and two instruments, so 4 moment conditions, and clusters = arange(n) % G. G center=True (the default) center=False 2 LinAlgError LinAlgError 3 LinAlgError LinAlgError 4 LinAlgError fits 5 fits fits The objective is a function of the inverse of the covariance of the moment conditions. The covariance that is estimated from G clusters is a sum of G outer products, and so it has rank at most G. The cluster sums of centered moment conditions add to zero, so its rank is at most G - 1. With 4 moment conditions it is singular for G = 4 and less if they are centered, and for G = 3 and less if they are not. Robust and kernel covariances have the same problem if there are too few observations, e.g., 4 observations and 4 moment conditions. fit now tests the rank of the covariance before it does anything else, with the function that the J and C statistics use, and raises a ValueError. For example, for G = 4: The clustered covariance of the 4 moment conditions is singular because there are too few clusters (4 clusters). Its rank is at most 3, one less than the number of clusters, since the cluster sums of centered moment conditions add to zero. 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. The covariance is estimated at the first step estimates of two-step GMM, which are 2SLS. Whether it is singular does not depend on the parameters, except for special values such as residuals that are all 0. The check is before the starting values are estimated, since that also inverts the covariance, and before the optimization starts, so there is no work to do before the error. The LinAlgError was raised after the optimizer started. A model with enough clusters is not changed, and its extra cost is one estimate of the covariance. The class and fit documentation say that the covariance has to be nonsingular. Not changed: * IVGMM, the two-step model, does not fail in the same way, and its estimation is not changed. With too few clusters it uses the inverse of a singular covariance, which is noise, and returns meaningless parameters, or it raises a LinAlgError if the covariance is exactly singular. For the model above and G = 3 the coefficient on the endogenous regressor is 1.70, and the true value is 0.5 and the estimate with G = 5 is 0.48. The J and C statistics of the previous commit are invalid. * A model with center=False and as many clusters as moment conditions still fits, as the report says. The covariance is not singular. The objective is then constant, and equal to the number of clusters, for every value of the parameters, so the estimates are the starting values. The new tests check the exact type of the error, since LinAlgError is a subclass of ValueError, the numbers of clusters and moment conditions in the message, the error with starting values, that the optimizer is not started, that models with enough clusters fit and have the objective of the definition as their J statistic, and the kernel and robust covariances with 4, 5 and 8 observations. 10 of the 20 new tests fail without the change. --- linearmodels/iv/model.py | 61 +++++++++ .../tests/iv/test_gmm_singular_covariance.py | 116 ++++++++++++++++++ 2 files changed, 177 insertions(+) diff --git a/linearmodels/iv/model.py b/linearmodels/iv/model.py index a6510b7c9d..38d72047ee 100644 --- a/linearmodels/iv/model.py +++ b/linearmodels/iv/model.py @@ -31,6 +31,7 @@ outer, sqrt, squeeze, + unique, ) from numpy.linalg import LinAlgError, eigvalsh, inv, matrix_rank, pinv, qr, solve from pandas import DataFrame, Series, concat @@ -1528,6 +1529,11 @@ class IVGMMCUE(_IVGMMBase): avoid starting value dependence. It is ``True`` unless ``center`` is passed in ``weight_config``. + 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``. + .. math:: \hat{\beta}_{cue} & =\min_{\beta}\bar{g}(\beta)'W(\beta)^{-1}g(\beta)\\ @@ -1762,6 +1768,51 @@ def _callback(params: ndarray[tuple[int], dtype[float64]]) -> None: return res.x[:, None], res.nit + 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. + """ + wy, wx, wz = self._wy, self._wx, self._wz + nobs, nmoment = wz.shape + 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) + def fit( self, *, @@ -1804,6 +1855,15 @@ def fit( IVGMMResults Results container + Raises + ------ + 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. + Notes ----- Starting values are computed by IVGMM. @@ -1815,6 +1875,7 @@ def fit( self._check_clusters(cov_config) wy, wx, wz = self._wy, self._wx, self._wz weight_matrix = self._weight.weight_matrix + self._check_moment_covariance() if starting is None: exog = None if self.exog.shape[1] == 0 else self.exog endog = None if self.endog.shape[1] == 0 else self.endog diff --git a/linearmodels/tests/iv/test_gmm_singular_covariance.py b/linearmodels/tests/iv/test_gmm_singular_covariance.py index 0ed9da1b9e..51853cf39f 100644 --- a/linearmodels/tests/iv/test_gmm_singular_covariance.py +++ b/linearmodels/tests/iv/test_gmm_singular_covariance.py @@ -16,6 +16,10 @@ centered J statistic is. All of the expected values below are computed from the definitions, without any of the weight matrix classes. +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 @@ -454,3 +458,115 @@ def test_units_of_the_instruments_do_not_matter(nclusters, scale): assert_allclose(res_scaled.c_stat().stat, res.c_stat().stat, rtol=1e-5) else: check_invalid(res_scaled.c_stat(), "C-statistic") + + +# ---------------------------------------------------------------------------- +# IVGMMCUE, which inverts the covariance in the objective that it minimizes +# ---------------------------------------------------------------------------- + + +def cue_model(data, nclusters, **config): + y, exog, endog, instr = data + clusters = np.arange(y.shape[0]) % nclusters + mod = IVGMMCUE( + y, + exog, + endog, + instr, + weight_type="clustered", + clusters=clusters, + **config, + ) + return mod, clusters + + +def cue_criterion(data, params, nclusters, center): + """The objective of continuously updating GMM at params, from the definition""" + x, y, z, _, _ = arrays(data) + g = z * (y - x @ params)[:, None] + gbar = g.mean(0) + cov = clustered_cov(np.arange(y.shape[0]) % nclusters, center)(g) + return y.shape[0] * gbar @ np.linalg.solve(cov, gbar) + + +# 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 +CUE_SINGULAR = [ + pytest.param(2, {}, id="default-2"), + pytest.param(3, {}, id="default-3"), + pytest.param(4, {}, id="default-4"), + pytest.param(4, {"center": True}, id="center-4"), + pytest.param(2, {"center": False}, id="uncentered-2"), + pytest.param(3, {"center": False}, id="uncentered-3"), +] + + +@pytest.mark.parametrize(("nclusters", "config"), CUE_SINGULAR) +def test_cue_too_few_clusters(nclusters, config): + # This was a LinAlgError (Singular matrix) from the objective function, which + # is a ValueError in recent versions of NumPy and says nothing about the cause + mod, clusters = cue_model(ONE_ENDOG, nclusters, **config) + with pytest.raises(ValueError, match="too few clusters") as exc_info: + mod.fit(cov_type="clustered", clusters=clusters) + assert exc_info.type is ValueError + message = str(exc_info.value) + assert f"{nclusters} clusters" in message + assert "4 moment conditions" in message + assert "more clusters" in message + assert "singular" in message + + +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"]) + with pytest.raises(ValueError, match="too few clusters") as exc_info: + mod.fit(starting=starting, cov_type="clustered", clusters=clusters) + assert exc_info.type is ValueError + + +def test_cue_too_few_clusters_is_found_before_the_optimization(monkeypatch): + def minimize(*args, **kwargs): + raise AssertionError("the objective is not defined, do not optimize it") + + monkeypatch.setattr(linearmodels.iv.model, "minimize", minimize) + mod, clusters = cue_model(ONE_ENDOG, 3) + with pytest.raises(ValueError, match="too few clusters"): + mod.fit(cov_type="clustered", clusters=clusters) + + +@pytest.mark.parametrize( + ("nclusters", "center"), + [(5, True), (6, True), (12, True), (4, False), (5, 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 + 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)) + params = res.params[["const", "w", "x"]].to_numpy() + expected = cue_criterion(ONE_ENDOG, params, nclusters, center) + assert_allclose(res.j_stat.stat, expected, rtol=1e-6) + # The default is to center the moment conditions + if center: + default, _ = cue_model(ONE_ENDOG, nclusters) + res_default = default.fit(cov_type="clustered", clusters=clusters) + assert_allclose(res_default.params, res.params) + + +@pytest.mark.parametrize("weight_type", ["kernel", "robust"]) +@pytest.mark.parametrize("nobs", [4, 5, 8]) +def test_cue_too_few_observations(weight_type, nobs): + # The rank of a covariance of centered moment conditions is at most nobs - 1 + data = small_sample(nobs) + mod = IVGMMCUE(*data, weight_type=weight_type) + if nobs > 4: + res = mod.fit() + assert np.isfinite(res.j_stat.stat) + return + 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) From f4f8a099a9c20909e3595d0a99ee42d601a37e65 Mon Sep 17 00:00:00 2001 From: Kevin Sheppard Date: Mon, 12 Oct 2026 09:33:50 +0100 Subject: [PATCH 3/3] TST: Allow the 3-cluster C-statistic test to see a LinAlgError The two Python 3.10 Azure jobs, which use NumPy 1.x, failed in test_c_stat_with_too_few_clusters_is_invalid[3] with "'NoneType' object has no attribute 'c_stat'". With 3 clusters the covariance of the model's own 4 moments is singular as well, and there NumPy raised LinAlgError when the model was estimated, which the fit helper of the test turns into None. The other tests in the file already return early in that case; this one did not. Checked by replacing inv() with one that raises for an ill-conditioned matrix, as NumPy did on those jobs: the test failed with the same AttributeError before this change and passes after it. It also passes with the real inv(), and the whole file passes (138 tests). --- linearmodels/tests/iv/test_gmm_singular_covariance.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/linearmodels/tests/iv/test_gmm_singular_covariance.py b/linearmodels/tests/iv/test_gmm_singular_covariance.py index 51853cf39f..2561ccb32c 100644 --- a/linearmodels/tests/iv/test_gmm_singular_covariance.py +++ b/linearmodels/tests/iv/test_gmm_singular_covariance.py @@ -229,6 +229,10 @@ def test_c_stat_with_too_few_clusters_is_invalid(nclusters): # 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 + return check_invalid(res.c_stat(), "C-statistic")