diff --git a/linearmodels/iv/model.py b/linearmodels/iv/model.py index 38d72047ee..a29dda3258 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,20 @@ 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 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:: * VCV: bootstrap @@ -1425,6 +1491,16 @@ 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. 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 -------- linearmodels.iv.gmm.IVGMMCovariance @@ -1448,6 +1524,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 @@ -1532,7 +1613,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:: @@ -1770,47 +1856,47 @@ 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, 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)): + 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 - 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." - ) + 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( @@ -1859,10 +1945,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/iv/results.py b/linearmodels/iv/results.py index cd0fda3123..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 @@ -1895,8 +1904,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..0f39e04492 100644 --- a/linearmodels/tests/iv/test_gmm_singular_covariance.py +++ b/linearmodels/tests/iv/test_gmm_singular_covariance.py @@ -18,13 +18,20 @@ 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. + +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 @@ -34,7 +41,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 +179,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 +192,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 +241,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 +274,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 +333,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 +350,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 +370,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 +401,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 +440,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: @@ -494,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"), @@ -521,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"]) @@ -541,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)) @@ -574,3 +651,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) 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)