diff --git a/linearmodels/iv/model.py b/linearmodels/iv/model.py
index 4eecedeab3..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
@@ -1005,6 +1006,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 +1115,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 +1124,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 +1150,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 +1167,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 +1444,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 +1465,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 +1476,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 +1485,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),
         }
 
@@ -1485,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)\\
@@ -1719,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,
         *,
@@ -1761,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.
@@ -1772,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
@@ -1795,14 +1899,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..2561ccb32c
--- /dev/null
+++ b/linearmodels/tests/iv/test_gmm_singular_covariance.py
@@ -0,0 +1,576 @@
+"""
+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.
+
+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.
+"""
+
+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)
+    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")
+
+
+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")
+
+
+# ----------------------------------------------------------------------------
+# 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)