From 841342d438244b75021b1ac4ea62859abcda6922 Mon Sep 17 00:00:00 2001 From: Kevin Sheppard Date: Sun, 11 Oct 2026 19:35:25 +0100 Subject: [PATCH 1/3] ENH: Rename RandomEffectsResults.wu_hausman to hausman and return a named tuple RandomEffectsResults.wu_hausman is the fixed effects vs random effects Hausman test added in #740. It has not been released, so it is renamed without a deprecation. The old name is easily confused with IVResults.wu_hausman, which is a different test and keeps its name. hausman now returns HausmanTestResults, a typing.NamedTuple with two fields: - statistic: a WaldTestStatistic, or an InvalidTestStatistic - estimates: the DataFrame with the columns b0, b1, b0-b1 and Std. Err. Unpacking still works, stat, estimates = re_res.hausman(fe_res), and the fields can also be read by name. The str and the repr show the statistic (value, p-value, distribution and null hypothesis) followed by the table of estimates. Unlike the repr of the test statistics, they do not end with the id of the object, so the result displays cleanly as the last value in a notebook cell: Hausman specification test H0: No systematic difference in coefficients between models Statistic: 7.1909 P-value: 0.0274 Distributed: chi2(2) Coefficients (b0: consistent under H0 and H1, b1: efficient under H0) b0 b1 b0-b1 Std. Err. exper 150.088403 149.344223 0.744180 2.117138 expersq -7.261883 -7.332176 0.070293 0.193418 HausmanTestResults is added to panel.results.__all__ and to the panel reference. The parameters of hausman are documented as optional. The method is excluded from the attribute access helper in the panel tests under its new name. The existing Hausman tests are renamed. The 9 new tests in test_results_hausman.py check that the old name is gone and that IVResults.wu_hausman is untouched, the type, the fields, indexing and unpacking, the str and repr for valid and invalid statistics, and that the results survive pickling. They fail on main, where neither the method nor the class exists. --- doc/source/panel/reference.rst | 1 + linearmodels/panel/results.py | 93 +++++++--- linearmodels/tests/panel/_utility.py | 2 +- linearmodels/tests/panel/test_results.py | 44 ++--- .../tests/panel/test_results_hausman.py | 167 ++++++++++++++++++ 5 files changed, 262 insertions(+), 45 deletions(-) create mode 100644 linearmodels/tests/panel/test_results_hausman.py diff --git a/doc/source/panel/reference.rst b/doc/source/panel/reference.rst index f80e112d0d..7c8cfba504 100644 --- a/doc/source/panel/reference.rst +++ b/doc/source/panel/reference.rst @@ -41,6 +41,7 @@ Estimation Results PanelResults PanelEffectsResults RandomEffectsResults + HausmanTestResults PanelModelComparison compare diff --git a/linearmodels/panel/results.py b/linearmodels/panel/results.py index 15bd5bf6ef..0f96ccbd93 100644 --- a/linearmodels/panel/results.py +++ b/linearmodels/panel/results.py @@ -5,7 +5,7 @@ from collections.abc import Mapping import datetime as dt from functools import cached_property -from typing import Any, Union +from typing import Any, NamedTuple, Union import warnings from formulaic.utils.context import capture_context @@ -30,6 +30,7 @@ __all__ = [ "FamaMacBethResults", + "HausmanTestResults", "PanelEffectsResults", "PanelResults", "RandomEffectsResults", @@ -965,6 +966,48 @@ def _has_unadjusted_cov(res: PanelResults) -> bool: return type(estimator) is HomoskedasticCovariance +class HausmanTestResults(NamedTuple): + """ + Results of a Hausman specification test. + + Parameters + ---------- + statistic : WaldTestStatistic or InvalidTestStatistic + The test statistic with its p-value, distribution and null hypothesis. + An InvalidTestStatistic is used if the statistic is negative, which + means that the models do not meet the assumptions of the test. + estimates : DataFrame + The coefficients included in the test. The columns are ``b0`` (from + the model that is consistent under both hypotheses, ``other``), ``b1`` + (from the model that is efficient under the null), their difference + ``b0-b1`` and ``Std. Err.``, the standard error of the difference. + The standard error is NaN if the variance of the difference is + negative. + + See Also + -------- + RandomEffectsResults.hausman + Compute the Hausman specification test. + + Notes + ----- + The results are a named tuple, so the statistic and the estimates can be + unpacked as in ``statistic, estimates = results`` or read using the + ``statistic`` and ``estimates`` fields. Printing the results shows the + statistic followed by the table of estimates. + """ + + statistic: WaldTestStatistic | InvalidTestStatistic + estimates: DataFrame + + def __str__(self) -> str: + legend = "b0: consistent under H0 and H1, b1: efficient under H0" + return f"{self.statistic}\n\nCoefficients ({legend})\n{self.estimates}" + + def __repr__(self) -> str: + return self.__str__() + + class RandomEffectsResults(PanelResults): """ Results container for random effect panel data models @@ -989,15 +1032,15 @@ def theta(self) -> DataFrame: """Values used in generalized demeaning""" return self._theta - def wu_hausman( + def hausman( self, other: PanelResults, include_constant: bool = False, sigmamore: bool = False, sigmaless: bool = False, - ) -> tuple[InvalidTestStatistic | WaldTestStatistic, DataFrame]: + ) -> HausmanTestResults: r""" - Hausman specification test against results from a consistent model + Hausman specification test against results from a consistent model. Parameters ---------- @@ -1006,29 +1049,30 @@ def wu_hausman( hypothesis is true, typically a fixed effects (entity effects) model estimated using the same observations and regressors as this model. - include_constant : bool + include_constant : bool, optional Flag indicating whether to include the constant in the comparison. The default is False. - sigmamore : bool + sigmamore : bool, optional Flag indicating whether to base both covariance estimates on the error variance estimated by this model, which is efficient under - the null. Cannot be used with ``sigmaless``. - sigmaless : bool + the null. Cannot be used with ``sigmaless``. The default is False. + sigmaless : bool, optional Flag indicating whether to base both covariance estimates on the error variance estimated by ``other``. Cannot be used with - ``sigmamore``. + ``sigmamore``. The default is False. Returns ------- - WaldTestStatistic or InvalidTestStatistic - Object containing test statistic, p-value, distribution and null. - An InvalidTestStatistic is returned when the statistic is negative. - DataFrame - The coefficients included in the test. The columns are ``b0`` - (from ``other``), ``b1`` (from this model), their difference - ``b0-b1`` and ``Std. Err.``, the standard error of the difference. - The standard error is NaN if the variance of the difference is - negative. + HausmanTestResults + A named tuple with the fields ``statistic``, which contains the + test statistic, p-value, distribution and null as a + WaldTestStatistic, and ``estimates``, a DataFrame with the + coefficients included in the test. An InvalidTestStatistic is + returned as the statistic when the statistic is negative. The + columns of ``estimates`` are ``b0`` (from ``other``), + ``b1`` (from this model), their difference ``b0-b1`` and + ``Std. Err.``, the standard error of the difference. The standard + error is NaN if the variance of the difference is negative. Raises ------ @@ -1081,8 +1125,12 @@ def wu_hausman( >>> from linearmodels import PanelOLS, RandomEffects >>> fe_res = PanelOLS(y, x, entity_effects=True).fit() >>> re_res = RandomEffects(y, x).fit() - >>> test, estimates = re_res.wu_hausman(fe_res) - >>> print(test) + >>> results = re_res.hausman(fe_res) + >>> print(results) + + The statistic and the estimates can be unpacked. + + >>> test, estimates = results """ if not isinstance(other, PanelResults): raise TypeError("other must be the results of a panel model.") @@ -1156,9 +1204,10 @@ def wu_hausman( "assumptions of the Hausman test.", name=name, ) - return invalid, estimates + return HausmanTestResults(invalid, estimates) null = "No systematic difference in coefficients between models" - return WaldTestStatistic(stat, null, len(common), name=name), estimates + wald = WaldTestStatistic(stat, null, len(common), name=name) + return HausmanTestResults(wald, estimates) PanelModelResults = Union[PanelEffectsResults, PanelResults, RandomEffectsResults] diff --git a/linearmodels/tests/panel/_utility.py b/linearmodels/tests/panel/_utility.py index e71d7ed61f..d6df8f2219 100644 --- a/linearmodels/tests/panel/_utility.py +++ b/linearmodels/tests/panel/_utility.py @@ -250,7 +250,7 @@ def assert_frame_similar(result, expected): def access_attributes(result): d = dir(result) for key in d: - if not key.startswith("_") and key not in ("wald_test", "wu_hausman"): + if not key.startswith("_") and key not in ("wald_test", "hausman"): val = getattr(result, key) if callable(val): val() diff --git a/linearmodels/tests/panel/test_results.py b/linearmodels/tests/panel/test_results.py index cdcff2feb8..4ff225a3c2 100644 --- a/linearmodels/tests/panel/test_results.py +++ b/linearmodels/tests/panel/test_results.py @@ -488,7 +488,7 @@ def re_fe(hausman_data): @pytest.mark.parametrize("variant", ["", "sigmamore", "sigmaless"]) @pytest.mark.parametrize("constant", [False, True], ids=["", "include_constant"]) -def test_wu_hausman_stata(re_fe, constant, variant): +def test_hausman_stata(re_fe, constant, variant): re_res, fe_res = re_fe opts = { "include_constant": constant, @@ -498,14 +498,14 @@ def test_wu_hausman_stata(re_fe, constant, variant): if constant: # The covariance difference is not positive definite with the constant with pytest.warns(UserWarning, match="not positive definite") as record: - wald, estimates = re_res.wu_hausman(other=fe_res, **opts) + wald, estimates = re_res.hausman(other=fe_res, **opts) assert all(issubclass(warn.category, UserWarning) for warn in record) names = ["const", *HAUSMAN_NAMES] assert estimates["Std. Err."].isna().any() else: with warnings.catch_warnings(): warnings.simplefilter("error") - wald, estimates = re_res.wu_hausman(other=fe_res, **opts) + wald, estimates = re_res.hausman(other=fe_res, **opts) names = HAUSMAN_NAMES assert estimates["Std. Err."].notna().all() expected_stat, expected_pval = STATA_HAUSMAN[ @@ -525,10 +525,10 @@ def test_wu_hausman_stata(re_fe, constant, variant): assert_allclose(estimates["b0-b1"], estimates["b0"] - estimates["b1"]) -def test_wu_hausman_formula(re_fe): +def test_hausman_formula(re_fe): re_res, fe_res = re_fe expected = hausman_by_hand(re_res, fe_res) - wald, estimates = re_res.wu_hausman(fe_res) + wald, estimates = re_res.hausman(fe_res) assert_allclose(wald.stat, expected) assert_allclose(wald.pval, stats.chi2.sf(expected, 2)) assert wald.df == 2 @@ -543,7 +543,7 @@ def test_wu_hausman_formula(re_fe): @pytest.mark.parametrize( "include_constant", [False, True], ids=["", "include_constant"] ) -def test_wu_hausman_constants(hausman_data, re_constant, fe_constant, include_constant): +def test_hausman_constants(hausman_data, re_constant, fe_constant, include_constant): dependent, regressors = hausman_data re_exog = add_constant(regressors) if re_constant else regressors fe_exog = add_constant(regressors) if fe_constant else regressors @@ -552,7 +552,7 @@ def test_wu_hausman_constants(hausman_data, re_constant, fe_constant, include_co with warnings.catch_warnings(): # Warnings are only possible when the constant is part of the test warnings.simplefilter("ignore", UserWarning) - wald, estimates = re_res.wu_hausman(fe_res, include_constant=include_constant) + wald, estimates = re_res.hausman(fe_res, include_constant=include_constant) # The constant can only be compared when both models have one names = ( ["const"] if include_constant and re_constant and fe_constant else [] @@ -563,7 +563,7 @@ def test_wu_hausman_constants(hausman_data, re_constant, fe_constant, include_co assert_allclose(wald.stat, expected, atol=1e-12) -def test_wu_hausman_negative_statistic(): +def test_hausman_negative_statistic(): # In small samples Var(b0) - Var(b1) is often indefinite, and the # statistic can be negative, in which case the test is not valid rg = np.random.default_rng(3) @@ -576,7 +576,7 @@ def test_wu_hausman_negative_statistic(): fe_res = PanelOLS(y, x, entity_effects=True).fit() re_res = RandomEffects(y, add_constant(x)).fit() with pytest.warns(UserWarning, match="not positive definite"): - test, estimates = re_res.wu_hausman(fe_res) + test, estimates = re_res.hausman(fe_res) assert isinstance(test, InvalidTestStatistic) assert np.isnan(test.stat) assert np.isnan(test.pval) @@ -586,19 +586,19 @@ def test_wu_hausman_negative_statistic(): assert hausman_by_hand(re_res, fe_res, ["x0", "x1"]) < 0 -def test_wu_hausman_conflicting_options(re_fe): +def test_hausman_conflicting_options(re_fe): re_res, fe_res = re_fe with pytest.raises(ValueError, match="cannot both be True"): - re_res.wu_hausman(fe_res, sigmamore=True, sigmaless=True) + re_res.hausman(fe_res, sigmamore=True, sigmaless=True) @pytest.mark.parametrize( "other", [None, "fe", np.arange(3)], ids=["none", "str", "array"] ) -def test_wu_hausman_invalid_other(re_fe, other): +def test_hausman_invalid_other(re_fe, other): re_res, _ = re_fe with pytest.raises(TypeError, match="other must be the results of a panel model"): - re_res.wu_hausman(other) + re_res.hausman(other) @pytest.mark.parametrize("which", ["other", "self"]) @@ -611,7 +611,7 @@ def test_wu_hausman_invalid_other(re_fe, other): ("autocorrelated", {}), ], ) -def test_wu_hausman_requires_unadjusted_cov(hausman_data, which, cov_type, cov_config): +def test_hausman_requires_unadjusted_cov(hausman_data, which, cov_type, cov_config): dependent, regressors = hausman_data exog = add_constant(regressors) re_kwargs = fe_kwargs = {} @@ -622,36 +622,36 @@ def test_wu_hausman_requires_unadjusted_cov(hausman_data, which, cov_type, cov_c re_res = RandomEffects(dependent, exog).fit(**re_kwargs) fe_res = PanelOLS(dependent, exog, entity_effects=True).fit(**fe_kwargs) with pytest.raises(TypeError, match="unadjusted covariance estimator"): - re_res.wu_hausman(fe_res) + re_res.hausman(fe_res) -def test_wu_hausman_conventional_covariance(hausman_data): +def test_hausman_conventional_covariance(hausman_data): # "conventional" and "homoskedastic" are aliases of the unadjusted estimator dependent, regressors = hausman_data exog = add_constant(regressors) re_res = RandomEffects(dependent, exog).fit(cov_type="homoskedastic") fe_res = PanelOLS(dependent, exog, entity_effects=True).fit(cov_type="conventional") - wald, _ = re_res.wu_hausman(fe_res) + wald, _ = re_res.hausman(fe_res) assert_allclose(wald.stat, STATA_HAUSMAN[""][0], atol=1e-9) -def test_wu_hausman_different_observations(hausman_data): +def test_hausman_different_observations(hausman_data): dependent, regressors = hausman_data exog = add_constant(regressors) re_res = RandomEffects(dependent, exog).fit() fe_res = PanelOLS(dependent.iloc[:-8], exog.iloc[:-8], entity_effects=True).fit() with pytest.raises(ValueError, match="same observations"): - re_res.wu_hausman(fe_res) + re_res.hausman(fe_res) -def test_wu_hausman_no_common_coefficients(hausman_data, data): +def test_hausman_no_common_coefficients(hausman_data, data): dependent, _ = hausman_data data = data.set_index(["nr", "year"]) re_res = RandomEffects(dependent, add_constant(data[HAUSMAN_NAMES])).fit() fe_res = PanelOLS(dependent, data[["married", "union"]], entity_effects=True).fit() with pytest.raises(ValueError, match="coefficients in common"): - re_res.wu_hausman(fe_res) + re_res.hausman(fe_res) # Only the constant is shared fe_const = PanelOLS(dependent, add_constant(data[["married"]]), entity_effects=True) with pytest.raises(ValueError, match="coefficients in common"): - re_res.wu_hausman(fe_const.fit()) + re_res.hausman(fe_const.fit()) diff --git a/linearmodels/tests/panel/test_results_hausman.py b/linearmodels/tests/panel/test_results_hausman.py new file mode 100644 index 0000000000..c7a2d9d0f2 --- /dev/null +++ b/linearmodels/tests/panel/test_results_hausman.py @@ -0,0 +1,167 @@ +""" +Tests of RandomEffectsResults.hausman and HausmanTestResults: the return type and +the checks that the models are comparable. +""" + +import pickle + +import numpy as np +from numpy.testing import assert_allclose +import pandas as pd +from pandas.testing import assert_frame_equal +import pytest +from scipy import stats +from statsmodels.tools.tools import add_constant + +from linearmodels.datasets import wage_panel +from linearmodels.iv.results import IVResults +from linearmodels.panel import PanelOLS, RandomEffects +import linearmodels.panel.results +from linearmodels.panel.results import HausmanTestResults, RandomEffectsResults +from linearmodels.shared.hypotheses import InvalidTestStatistic, WaldTestStatistic + +HAUSMAN_NAMES = ["exper", "expersq"] + + +@pytest.fixture(scope="module") +def wage(): + return wage_panel.load().set_index(["nr", "year"]) + + +@pytest.fixture(scope="module") +def hausman_results(wage): + exog = add_constant(wage[HAUSMAN_NAMES]) + re_res = RandomEffects(wage["hours"], exog).fit() + fe_res = PanelOLS(wage["hours"], exog, entity_effects=True).fit() + return re_res, fe_res + + +def hausman_by_hand(re_res, fe_res, names=HAUSMAN_NAMES): + delta = (fe_res.params - re_res.params)[names].to_numpy() + diff = (fe_res.cov - re_res.cov).loc[names, names].to_numpy() + return float(delta @ np.linalg.pinv(diff) @ delta) + + +def negative_statistic_results(): + # In small samples Var(b0) - Var(b1) is often indefinite and the statistic + # can be negative + rg = np.random.default_rng(3) + index = pd.MultiIndex.from_product([np.arange(12), np.arange(3)]) + effects = np.repeat(rg.standard_normal(12), 3) + x = pd.DataFrame(rg.standard_normal((36, 2)), index=index, columns=["x0", "x1"]) + y = pd.Series( + x.sum(axis=1).to_numpy() + effects + rg.standard_normal(36), index=index + ) + fe_res = PanelOLS(y, x, entity_effects=True).fit() + re_res = RandomEffects(y, add_constant(x)).fit() + return re_res, fe_res + + +def test_hausman_replaces_wu_hausman(hausman_results): + re_res, _ = hausman_results + # The method was never released, so there is no alias + assert not hasattr(RandomEffectsResults, "wu_hausman") + assert not hasattr(re_res, "wu_hausman") + assert callable(re_res.hausman) + # IVResults.wu_hausman is a different test that keeps its name + assert callable(IVResults.wu_hausman) + + +def test_hausman_results_exported(): + assert "HausmanTestResults" in linearmodels.panel.results.__all__ + + +def test_hausman_named_tuple(hausman_results): + re_res, fe_res = hausman_results + result = re_res.hausman(fe_res) + assert isinstance(result, HausmanTestResults) + assert isinstance(result, tuple) + assert len(result) == 2 + assert result._fields == ("statistic", "estimates") + + # Field access and indexing + assert isinstance(result.statistic, WaldTestStatistic) + assert not isinstance(result.statistic, InvalidTestStatistic) + assert isinstance(result.estimates, pd.DataFrame) + assert result[0] is result.statistic + assert result[1] is result.estimates + assert_allclose(result.statistic.stat, hausman_by_hand(re_res, fe_res)) + assert list(result.estimates.columns) == ["b0", "b1", "b0-b1", "Std. Err."] + assert list(result.estimates.index) == HAUSMAN_NAMES + + # Unpacking + statistic, estimates = re_res.hausman(fe_res) + assert isinstance(statistic, WaldTestStatistic) + assert_frame_equal(estimates, result.estimates) + assert_allclose(statistic.stat, result.statistic.stat) + + +def test_hausman_str_repr(hausman_results): + re_res, fe_res = hausman_results + result = re_res.hausman(fe_res) + text = str(result) + assert repr(result) == text + # The statistic comes first, followed by the estimates + assert text.startswith(str(result.statistic)) + assert "Hausman specification test" in text + assert "H0: No systematic difference in coefficients between models" in text + assert f"Statistic: {result.statistic.stat:0.4f}" in text + assert f"P-value: {result.statistic.pval:0.4f}" in text + assert "Distributed: chi2(2)" in text + assert text.index("Statistic:") < text.index("Coefficients") + assert text.endswith(str(result.estimates)) + for label in [*result.estimates.columns, *result.estimates.index]: + assert label in text + # No object ids, which would make the repr change from run to run + assert "id:" not in text + assert "0x" not in text + + +def test_hausman_str_repr_invalid(): + re_res, fe_res = negative_statistic_results() + with pytest.warns(UserWarning, match="not positive definite"): + result = re_res.hausman(fe_res) + assert isinstance(result.statistic, InvalidTestStatistic) + text = str(result) + assert repr(result) == text + assert text.startswith(str(result.statistic)) + assert "Invalid test statistic" in text + assert "negative" in text + assert text.endswith(str(result.estimates)) + assert "x0" in text + assert "x1" in text + + +@pytest.mark.parametrize("protocol", [2, pickle.HIGHEST_PROTOCOL]) +def test_hausman_pickle(hausman_results, protocol): + re_res, fe_res = hausman_results + result = re_res.hausman(fe_res) + loaded = pickle.loads(pickle.dumps(result, protocol=protocol)) + assert type(loaded) is HausmanTestResults + assert isinstance(loaded.statistic, WaldTestStatistic) + assert loaded.statistic.stat == result.statistic.stat + assert loaded.statistic.pval == result.statistic.pval + assert loaded.statistic.df == result.statistic.df + assert loaded.statistic.dist_name == result.statistic.dist_name + assert loaded.statistic.null == result.statistic.null + assert_frame_equal(loaded.estimates, result.estimates) + assert str(loaded) == str(result) + + +def test_hausman_pickle_invalid(): + re_res, fe_res = negative_statistic_results() + with pytest.warns(UserWarning, match="not positive definite"): + result = re_res.hausman(fe_res) + loaded = pickle.loads(pickle.dumps(result)) + assert type(loaded) is HausmanTestResults + assert isinstance(loaded.statistic, InvalidTestStatistic) + assert np.isnan(loaded.statistic.stat) + assert np.isnan(loaded.statistic.pval) + assert_frame_equal(loaded.estimates, result.estimates) + assert str(loaded) == str(result) + + +def test_hausman_pval(hausman_results): + re_res, fe_res = hausman_results + result = re_res.hausman(fe_res) + assert_allclose(result.statistic.pval, stats.chi2.sf(result.statistic.stat, 2)) From 980dae373f351581bd54fd8146ed8af52f4dda2c Mon Sep 17 00:00:00 2001 From: Kevin Sheppard Date: Sun, 11 Oct 2026 19:35:50 +0100 Subject: [PATCH 2/3] BUG: Check that the models compared by the Hausman test use the same data hausman only checked that the two models have the same number of observations. Probing it on wage_panel with the probe scripts of the 7.1 release plan (probe_newfeat.py and probe_hausman2.py) found that it also accepted models that cannot be compared: - A weighted fixed effects model against an unweighted random effects model returned 118.8. The statistic is 112.1 when neither model is weighted and 109.6 when both are. The models then estimate different quantities, so the difference between their coefficients says nothing about the effects. A ValueError is now raised when the weights differ. The weights are compared after the rescaling to mean 1 that the models apply, so weights that only differ by scale are accepted. - A model of a different dependent variable (hours against lwage) was accepted. A ValueError is now raised when the dependent variables differ. The values are compared, not the names. - Models with the same number of observations but not the same observations, for example one that dropped the first entity and one that dropped the last, were accepted. A ValueError is now raised when the entities and periods of the observations differ. The probes print nan for the different dependent variable and for the random effects model against PooledOLS. That is InvalidTestStatistic.stat, which is always nan: both cases already returned an InvalidTestStatistic on main, with the reason that the statistic is negative. What was missing is the case where the statistic is not finite. A nan statistic is not less than 0, so a WaldTestStatistic with a nan statistic and a nan p-value was returned when a covariance or a coefficient is nan or infinite. hausman now returns an InvalidTestStatistic, with the reason that the coefficients or covariances are not finite or that the statistic is not finite, for these. A random effects model compared with a model that has time effects is still accepted, and this is now documented. RE has no time effects, so the null hypothesis that the test examines is then that the regressors are uncorrelated with the entity effects and with the time effects, and wage_panel gives 215.9 for FE with both effects against RE without time effects. Stata's hausman leaves this to the user. Rejecting it would also rule out the standard way to test the entity effects with year effects: include the year dummies in the random effects model and use the model with entity and time effects as other. The comparison then uses the coefficients that are in both models, and the tests cover this. The tests cover each ValueError, a renamed dependent variable, weights that differ by scale, nan and inf in the covariance and in the coefficients, the random effects model against PooledOLS, which is an InvalidTestStatistic on main and now, and the comparison with a model with time effects. 8 of the 12 tests that are added to test_results_hausman.py by this commit fail without the checks: the three tests of the weights, the different dependent variable, the different observations and the three tests of the nan and inf values. The other 4 describe behavior that does not change. --- linearmodels/panel/results.py | 94 +++++++--- .../tests/panel/test_results_hausman.py | 167 +++++++++++++++++- 2 files changed, 237 insertions(+), 24 deletions(-) diff --git a/linearmodels/panel/results.py b/linearmodels/panel/results.py index 0f96ccbd93..35b11d2efc 100644 --- a/linearmodels/panel/results.py +++ b/linearmodels/panel/results.py @@ -966,6 +966,36 @@ def _has_unadjusted_cov(res: PanelResults) -> bool: return type(estimator) is HomoskedasticCovariance +def _check_same_sample(res: PanelResults, other: PanelResults) -> None: + """Check that two sets of results use the same data and weights""" + if res.nobs != other.nobs: + raise ValueError( + "The models must be estimated using the same observations. " + f"This model uses {res.nobs} and other uses {other.nobs}." + ) + model, other_model = res.model, other.model + if not model.dependent.index.equals(other_model.dependent.index): + raise ValueError( + "The models must be estimated using the same observations. The " + "entities and time periods of the observations differ." + ) + if not np.allclose( + model.dependent.values2d, other_model.dependent.values2d, rtol=1e-10, atol=0.0 + ): + raise ValueError( + "The models must have the same dependent variable, but the " + "dependent variables of the models differ." + ) + # Weights are rescaled to have mean 1, so only the relative weights matter + if not np.allclose( + model.weights.values2d, other_model.weights.values2d, rtol=1e-10, atol=0.0 + ): + raise ValueError( + "The models must be estimated using the same weights, but the " + "weights of the models differ." + ) + + class HausmanTestResults(NamedTuple): """ Results of a Hausman specification test. @@ -974,8 +1004,9 @@ class HausmanTestResults(NamedTuple): ---------- statistic : WaldTestStatistic or InvalidTestStatistic The test statistic with its p-value, distribution and null hypothesis. - An InvalidTestStatistic is used if the statistic is negative, which - means that the models do not meet the assumptions of the test. + An InvalidTestStatistic is used if the statistic is negative or not + finite, which means that the models do not meet the assumptions of the + test. estimates : DataFrame The coefficients included in the test. The columns are ``b0`` (from the model that is consistent under both hypotheses, ``other``), ``b1`` @@ -1047,8 +1078,8 @@ def hausman( other : PanelResults Results from a model that is consistent whether or not the null hypothesis is true, typically a fixed effects (entity effects) - model estimated using the same observations and regressors as - this model. + model estimated using the same observations, dependent variable, + regressors and weights as this model. include_constant : bool, optional Flag indicating whether to include the constant in the comparison. The default is False. @@ -1068,8 +1099,8 @@ def hausman( test statistic, p-value, distribution and null as a WaldTestStatistic, and ``estimates``, a DataFrame with the coefficients included in the test. An InvalidTestStatistic is - returned as the statistic when the statistic is negative. The - columns of ``estimates`` are ``b0`` (from ``other``), + returned as the statistic when the statistic is negative or not + finite. The columns of ``estimates`` are ``b0`` (from ``other``), ``b1`` (from this model), their difference ``b0-b1`` and ``Std. Err.``, the standard error of the difference. The standard error is NaN if the variance of the difference is negative. @@ -1081,8 +1112,9 @@ def hausman( was not estimated using the unadjusted covariance estimator. ValueError If ``sigmamore`` and ``sigmaless`` are both True, if the models - use different numbers of observations, or if the models do not - have any coefficients in common. + were not estimated using the same observations, dependent variable + and weights, or if the models do not have any coefficients in + common. Notes ----- @@ -1108,12 +1140,26 @@ def hausman( The result is only meaningful when :math:`V_0 - V_1` is positive definite, which is guaranteed in large samples if the null is true. If it is not, a warning is issued and a generalized inverse is used. If - the resulting statistic is negative, the data do not satisfy the - assumptions of the test and an InvalidTestStatistic is returned. - - Both models must use the default unadjusted covariance estimator, and - the statistic requires that ``other`` is estimated using the same - observations as this model. + the resulting statistic is negative or not finite, the data do not + satisfy the assumptions of the test and an InvalidTestStatistic is + returned. This can happen, for example, if ``other`` is a pooled + model, since the unadjusted covariance of a pooled model ignores the + entity effects. + + Both models must use the default unadjusted covariance estimator. The + statistic also requires that the models are estimated using the same + observations, the same dependent variable and the same weights, which + are checked. The regressors are not checked: only coefficients with + the same name in both models are compared. + + Effects that are included in ``other`` but not in this model are not + tested separately. For example, if ``other`` includes time effects + and this model does not, the null hypothesis is that the regressors + are uncorrelated with the entity effects and with the time effects, + and the test can reject because of the time effects. To test only the + entity effects, include the time dummies as regressors in this model. + Only the coefficients found in both models are compared, so the + coefficients on the time dummies are not part of the test. References ---------- @@ -1141,11 +1187,7 @@ def hausman( "The Hausman test can only be used with models estimated using " "the unadjusted covariance estimator." ) - if self.nobs != other.nobs: - raise ValueError( - "The models must be estimated using the same observations. " - f"This model uses {self.nobs} and other uses {other.nobs}." - ) + _check_same_sample(self, other) in_other = set(other.params.index) common = [name for name in self.params.index if name in in_other] @@ -1181,6 +1223,14 @@ def hausman( diff = var_diff.to_numpy() delta = b_diff.to_numpy() + name = "Hausman specification test" + if not (np.all(np.isfinite(diff)) and np.all(np.isfinite(delta))): + invalid = InvalidTestStatistic( + "The coefficients or covariances of the models are not finite, " + "so the statistic cannot be computed.", + name=name, + ) + return HausmanTestResults(invalid, estimates) try: np.linalg.cholesky(diff) positive_definite = True @@ -1197,10 +1247,10 @@ def hausman( ) stat = float(delta @ np.linalg.pinv(diff, hermitian=True) @ delta) - name = "Hausman specification test" - if stat < 0: + if not np.isfinite(stat) or stat < 0: + problem = "not finite" if not np.isfinite(stat) else "negative" invalid = InvalidTestStatistic( - "The statistic is negative, so the models do not meet the " + f"The statistic is {problem}, so the models do not meet the " "assumptions of the Hausman test.", name=name, ) diff --git a/linearmodels/tests/panel/test_results_hausman.py b/linearmodels/tests/panel/test_results_hausman.py index c7a2d9d0f2..3b8dfd5049 100644 --- a/linearmodels/tests/panel/test_results_hausman.py +++ b/linearmodels/tests/panel/test_results_hausman.py @@ -4,6 +4,7 @@ """ import pickle +import warnings import numpy as np from numpy.testing import assert_allclose @@ -15,9 +16,13 @@ from linearmodels.datasets import wage_panel from linearmodels.iv.results import IVResults -from linearmodels.panel import PanelOLS, RandomEffects +from linearmodels.panel import PanelOLS, PooledOLS, RandomEffects import linearmodels.panel.results -from linearmodels.panel.results import HausmanTestResults, RandomEffectsResults +from linearmodels.panel.results import ( + HausmanTestResults, + PanelEffectsResults, + RandomEffectsResults, +) from linearmodels.shared.hypotheses import InvalidTestStatistic, WaldTestStatistic HAUSMAN_NAMES = ["exper", "expersq"] @@ -165,3 +170,161 @@ def test_hausman_pval(hausman_results): re_res, fe_res = hausman_results result = re_res.hausman(fe_res) assert_allclose(result.statistic.pval, stats.chi2.sf(result.statistic.stat, 2)) + + +def weighted_results(wage, re_weights, fe_weights): + n = wage.shape[0] + base = pd.Series(1.0 + (np.arange(n) % 5), index=wage.index) + exog = add_constant(wage[HAUSMAN_NAMES]) + re_res = RandomEffects(wage["hours"], exog, weights=re_weights(base)).fit() + fe_res = PanelOLS( + wage["hours"], exog, entity_effects=True, weights=fe_weights(base) + ).fit() + return re_res, fe_res + + +def unweighted(base): + return None + + +def weighted(base): + return base + + +@pytest.mark.parametrize( + ("re_weights", "fe_weights"), + [(unweighted, weighted), (weighted, unweighted)], + ids=["fe-weighted", "re-weighted"], +) +def test_hausman_different_weights(wage, re_weights, fe_weights): + re_res, fe_res = weighted_results(wage, re_weights, fe_weights) + with pytest.raises(ValueError, match="same weights"): + re_res.hausman(fe_res) + + +def test_hausman_different_weight_values(wage): + re_res, fe_res = weighted_results(wage, weighted, lambda base: base[::-1].array) + with pytest.raises(ValueError, match="same weights"): + re_res.hausman(fe_res) + + +def test_hausman_same_weights(wage): + re_res, fe_res = weighted_results(wage, weighted, weighted) + result = re_res.hausman(fe_res) + assert isinstance(result.statistic, WaldTestStatistic) + assert not isinstance(result.statistic, InvalidTestStatistic) + assert_allclose(result.statistic.stat, hausman_by_hand(re_res, fe_res)) + # The weights are only relevant up to scale + scaled = weighted_results(wage, weighted, lambda base: 3.0 * base) + scaled_result = scaled[0].hausman(scaled[1]) + assert_allclose(scaled_result.statistic.stat, result.statistic.stat, rtol=1e-8) + + +def test_hausman_different_dependent(wage): + exog = add_constant(wage[HAUSMAN_NAMES]) + fe_res = PanelOLS(wage["lwage"], exog, entity_effects=True).fit() + re_res = RandomEffects(wage["hours"], exog).fit() + with pytest.raises(ValueError, match="same dependent variable"): + re_res.hausman(fe_res) + # A transformed dependent variable is also a different dependent variable + re_res = RandomEffects(np.exp(wage["lwage"]), exog).fit() + with pytest.raises(ValueError, match="same dependent variable"): + re_res.hausman(fe_res) + + +def test_hausman_renamed_dependent(wage): + exog = add_constant(wage[HAUSMAN_NAMES]) + fe_res = PanelOLS(wage["hours"], exog, entity_effects=True).fit() + re_res = RandomEffects(wage["hours"].rename("other_name"), exog).fit() + result = re_res.hausman(fe_res) + assert isinstance(result.statistic, WaldTestStatistic) + assert_allclose(result.statistic.stat, hausman_by_hand(re_res, fe_res)) + + +def test_hausman_different_rows_same_count(wage): + # The first entity has 8 observations, as does the last entity + exog = add_constant(wage[HAUSMAN_NAMES]) + nobs = wage.shape[0] + first = RandomEffects(wage["hours"].iloc[8:], exog.iloc[8:]).fit() + last = PanelOLS( + wage["hours"].iloc[: nobs - 8], exog.iloc[: nobs - 8], entity_effects=True + ).fit() + assert first.nobs == last.nobs + with pytest.raises(ValueError, match="same observations"): + first.hausman(last) + + +@pytest.mark.parametrize("bad", [np.nan, np.inf]) +def test_hausman_non_finite_is_invalid(wage, bad): + exog = add_constant(wage[HAUSMAN_NAMES]) + re_res = RandomEffects(wage["hours"], exog).fit() + fe_res = PanelOLS(wage["hours"], exog, entity_effects=True).fit() + fe_res.cov = fe_res.cov * bad + with warnings.catch_warnings(): + warnings.simplefilter("error") + result = re_res.hausman(fe_res) + assert isinstance(result.statistic, InvalidTestStatistic) + assert np.isnan(result.statistic.stat) + assert np.isnan(result.statistic.pval) + assert "not finite" in str(result.statistic) + assert list(result.estimates.index) == HAUSMAN_NAMES + assert not np.isfinite(result.estimates["Std. Err."]).any() + assert "Invalid test statistic" in str(result) + + +def test_hausman_non_finite_coefficients_is_invalid(wage, monkeypatch): + exog = add_constant(wage[HAUSMAN_NAMES]) + re_res = RandomEffects(wage["hours"], exog).fit() + fe_res = PanelOLS(wage["hours"], exog, entity_effects=True).fit() + params = fe_res.params.copy() + params["exper"] = np.nan + monkeypatch.setattr(PanelEffectsResults, "params", property(lambda self: params)) + with warnings.catch_warnings(): + warnings.simplefilter("error") + result = re_res.hausman(fe_res) + assert isinstance(result.statistic, InvalidTestStatistic) + assert "not finite" in str(result.statistic) + + +def test_hausman_pooled_is_invalid(wage): + # The unadjusted covariance of the pooled model ignores the entity effects + # and it is not consistent under the alternative + exog = add_constant(wage[["expersq", "union", "married"]]) + re_res = RandomEffects(wage["lwage"], exog).fit() + pooled = PooledOLS(wage["lwage"], exog).fit() + with pytest.warns(UserWarning, match="not positive definite"): + result = re_res.hausman(pooled) + assert isinstance(result.statistic, InvalidTestStatistic) + assert np.isnan(result.statistic.stat) + assert np.isnan(result.statistic.pval) + assert "negative" in str(result.statistic) + assert "do not meet the assumptions of the Hausman test" in str(result.statistic) + assert list(result.estimates.index) == ["expersq", "union", "married"] + + +def test_hausman_two_way_fixed_effects(wage): + # A model with time effects is accepted as other. The test then also + # examines the time effects, so the dummies should be in this model. + years = pd.get_dummies( + wage.index.get_level_values("year"), prefix="y", drop_first=True, dtype=float + ) + years.index = wage.index + names = ["expersq", "union", "married"] + exog = add_constant(wage[names]) + fe_res = PanelOLS(wage["lwage"], exog, entity_effects=True, time_effects=True).fit() + + re_plain = RandomEffects(wage["lwage"], exog).fit() + plain = re_plain.hausman(fe_res) + assert isinstance(plain.statistic, WaldTestStatistic) + assert list(plain.estimates.index) == names + + re_dummies = RandomEffects(wage["lwage"], pd.concat([exog, years], axis=1)).fit() + with_dummies = re_dummies.hausman(fe_res) + assert isinstance(with_dummies.statistic, WaldTestStatistic) + # Only the coefficients that are in both models are compared + assert list(with_dummies.estimates.index) == names + assert with_dummies.statistic.df == len(names) + assert_allclose( + with_dummies.statistic.stat, hausman_by_hand(re_dummies, fe_res, names) + ) + assert not np.isclose(with_dummies.statistic.stat, plain.statistic.stat) From b7e2a72fd7d84eb9e5222e4fadc874890f4d961c Mon Sep 17 00:00:00 2001 From: Kevin Sheppard Date: Sun, 11 Oct 2026 19:36:19 +0100 Subject: [PATCH 3/3] BUG: Return InapplicableTestStatistic from pesaran_cd where the test is not valid PanelResults.pesaran_cd returned a statistic for every model. The Pesaran CD statistic is not valid for the models below, so pesaran_cd now returns an InapplicableTestStatistic whose reason says why: - BetweenOLS is fitted to the entity means, so it has no time-varying residuals, and the series implied by its coefficients is not the model's error. - FamaMacBeth: the errors of the model are the residuals of the period-by-period cross-sectional regressions. The series implied by the averaged coefficients contains the common, time-varying term (beta_t - beta_bar)'x_it. With coefficients that vary over time and independent errors, the statistic of that series rejects in 72.5% to 100% of the samples (mean statistic 4.6 to 51.8). The period-by-period residuals sum to zero in every period and have the problem of time effects below. Fama-MacBeth inference is also already robust to cross-sectional dependence. With constant coefficients the statistic has the right size (4.2% to 7.3%), so this exclusion rests on what the series is and not on its size. - Models with time effects. They demean the residuals across entities in every period, which makes the residuals negatively correlated across entities (Juodis and Reese, 2022, JBES 40(3)). With independent errors the statistic is centered at about -sqrt(T/2) with a standard deviation of 0.03 to 0.13, so it rejects in 98% to 100% of the samples when T is 10 or more and never rejects when T is 5. - Models with other effects that change over time for some entity and have a group of several entities that is missing from some periods. The panel example notebook specifies time effects as other_effects=time ids, which is the same model as time_effects=True. Industry-by-period effects have the same problem, and so do blocks of periods, also when they are combined with groups of entities (23% to 28% of rejections at T = 20 with blocks of 2 periods, and a statistic that is biased downwards for smaller T). Other effects that do not change over time for any entity, such as industry or region, do not change the correlations over time. Nor do effects that are present in every period, such as random groups, which reject as often as the other models in the simulation, or groups with a single entity. The other models were checked by simulation, with independent errors, y = 1 + 0.5 x + a_i + e_it, N = 20 and 100 entities, T = 5, 10 and 20 periods, 1,000 samples per cell (400 for the models with time effects) and the 5% test. Rejection rate in %, from scripts that are not part of this change: N=20 N=100 T=5 T=10 T=20 T=5 T=10 T=20 PooledOLS 7.2 5.0 5.1 5.6 5.6 4.2 PanelOLS, entity effects 7.5 4.8 5.2 5.6 5.7 4.2 RandomEffects 7.2 4.9 5.2 5.6 5.6 4.2 FirstDifferenceOLS 7.3 4.8 5.4 5.3 5.8 4.1 PanelOLS, other = groups of entities 7.3 5.0 5.2 5.7 5.6 4.3 PanelOLS, entity effects and random groups 7.3 5.0 5.2 5.4 5.6 4.2 BetweenOLS 7.0 5.5 4.9 5.8 5.8 5.8 FamaMacBeth, constant coefs 7.3 5.7 5.4 5.7 5.6 4.2 PanelOLS, time effects 0.0 98.2 100.0 0.0 99.2 100.0 PanelOLS, entity and time 0.0 98.0 100.0 0.0 99.0 100.0 PanelOLS, entity and industry-by-period effects 0.0 99.5 100.0 0.0 98.2 100.0 PanelOLS, entity effects and blocks of 2 periods 2.5 4.2 27.5 1.0 4.5 23.5 The applicable models have the same size because the correlation of two series over time removes anything that is constant over time, so the entity effects and the constant in the shocks do not matter. The size is 4% to 6% except for T = 5 and N = 20, where it is 7.2% to 7.5%. For pure noise the size is 5.9% there, and 6.6%, 7.1% to 7.5% and 8.6% to 11.7% for T = 4, 3 and 2: the variance of the statistic is T / (T - 1) for series that are demeaned, and the test needs more than a few periods. This is not changed. Not changed, because the size under independent errors is acceptable, but the docstring now says what is used: - FirstDifferenceOLS uses the shocks in levels (idiosyncratic) and not the residuals of the differenced equation. The size is 4.1% to 7.3% with independent errors and 6.4% to 8.1% with the differenced residuals. With random walk errors, for which first differences are efficient, it is 7.5% to 28.0% for the shocks in levels and 3.8% to 7.8% for the differenced residuals. (plm::pcdtest cannot be applied to an fd model: it fails with "argument 1 is not a vector". Run on a pooled regression of the differenced data it equals the statistic of the differenced residuals.) - The weights of a weighted model are not applied to the shocks. The size is 3.6% to 7.9% for weights that vary by observation or by entity, and up to 9.8% when the weights vary over time in the same way for all entities. Multiplying the shocks by the square root of the weights, which is what plm::pcdtest does for a pooled model, gives 3.4% to 4.8% in that case. - The statistic needs serially uncorrelated shocks. With random walk errors every model that uses the shocks in levels rejects in 25% to 28% of the samples at T = 20. The docstring explains each case, with the references (Pesaran 2021, Juodis and Reese 2022), and the memory use: the dense N by N matrices of correlations and overlapping observations take about 450 MB when N is 4,000. The tests check that the models above are applicable or not, the reasons, and that the statistic of an applicable model is the one computed by explicit loops over pairs of entities. test_pesaran_cd_all_estimators keeps the estimators where the test applies. 26 of the 39 tests in test_results_pesaran_cd.py fail on main, where every model that is not applicable returns a NormalTestStatistic, and the other 13 pass on both. --- linearmodels/panel/results.py | 135 +++++++++- linearmodels/tests/panel/test_results.py | 20 +- .../tests/panel/test_results_pesaran_cd.py | 251 ++++++++++++++++++ 3 files changed, 380 insertions(+), 26 deletions(-) create mode 100644 linearmodels/tests/panel/test_results_pesaran_cd.py diff --git a/linearmodels/panel/results.py b/linearmodels/panel/results.py index 35b11d2efc..c85531cee9 100644 --- a/linearmodels/panel/results.py +++ b/linearmodels/panel/results.py @@ -604,6 +604,19 @@ def wresids(self) -> Series: self._wresids.squeeze(), index=self._index, name="weighted residual" ) + def _pesaran_cd_inapplicable(self) -> str | None: + """Reason why the Pesaran CD test is not valid for the model, if any""" + from linearmodels.panel.model import BetweenOLS # noqa: PLC0415 + + if isinstance(self.model, BetweenOLS): + return ( + "The Pesaran CD test is not applicable to BetweenOLS. The model " + "is fitted to the entity means, so it has no time-varying " + "residuals, and the series implied by its coefficients is not " + "the model's error." + ) + return None + @cached_property def pesaran_cd(self) -> InapplicableTestStatistic | NormalTestStatistic: r""" @@ -611,17 +624,19 @@ def pesaran_cd(self) -> InapplicableTestStatistic | NormalTestStatistic: Returns ------- - NormalTestStatistic - Statistic value, distribution and p-value. + NormalTestStatistic or InapplicableTestStatistic + Statistic value, distribution and p-value. If the test is not + valid for the model, an InapplicableTestStatistic that explains + why in its reason. Notes ----- Tests the null hypothesis that the model's idiosyncratic shocks are - cross-sectionally independent. Let :math:`\hat{\rho}_{ij}` denote the - sample correlation of the estimated idiosyncratic shocks for entities - :math:`i` and :math:`j`, computed using their common time observations, - and let :math:`T_{ij}` denote the number of overlapping observations. - The statistic is + cross-sectionally independent [1]_. Let :math:`\hat{\rho}_{ij}` denote + the sample correlation of the estimated idiosyncratic shocks for + entities :math:`i` and :math:`j`, computed using their common time + observations, and let :math:`T_{ij}` denote the number of overlapping + observations. The statistic is .. math:: @@ -631,11 +646,55 @@ def pesaran_cd(self) -> InapplicableTestStatistic | NormalTestStatistic: overlapping observations and finite pairwise correlation estimates. Under the null, the statistic is asymptotically standard normal. + The shocks are the model's ``idiosyncratic`` values, so the weights of + a weighted model are not applied to them, and the shocks of + ``FirstDifferenceOLS`` are in levels rather than the residuals of the + differenced equation. Since the correlations are computed over time, + the test assumes that the shocks are not serially correlated, and it + over-rejects if they are persistent. + + The test is not valid for every model and an InapplicableTestStatistic + is returned for the following. + + * ``BetweenOLS``. The model is fitted to the entity means, so it has + no time-varying residuals, and the series implied by its + coefficients is not the model's error. + * ``FamaMacBeth``. The model's errors are the residuals of the + period-by-period cross-sectional regressions. The series implied by + the averaged coefficients contains the common, time-varying term + :math:`(\beta_t - \bar{\beta})^{\prime}x_{it}`, so it is + cross-sectionally dependent by construction. Fama-MacBeth inference + is also already robust to cross-sectional dependence. + * Models with time effects, and models with other effects that + change over time for some entity and have a group of several + entities that is missing from some periods, such as period, + industry-by-period or decade effects. These effects demean the + residuals across entities that are observed in the same periods, + which makes the residuals negatively correlated across entities + and biases the statistic downwards. For example, when the model + has time effects and the errors are independent, the statistic is + centered at about :math:`-\sqrt{T/2}` and its standard deviation + is far below 1, so it is not standard normal [2]_. Other effects + that do not change over time for an entity, or that are present + in every period, are not affected. + + The statistic is computed from dense N by N matrices of pairwise + correlations and numbers of overlapping observations, where N is the + number of entities. Computing them takes about 450 MB of memory when + N is 4,000, and the memory use is proportional to the square of N. + References ---------- - Pesaran, M. H. (2021). General diagnostic tests for cross-sectional - dependence in panels. Empirical Economics, 60(1), 13-50. + .. [1] Pesaran, M. H. (2021). General diagnostic tests for + cross-sectional dependence in panels. Empirical Economics, 60(1), + 13-50. + .. [2] Juodis, A. and Reese, S. (2022). The incidental parameters + problem in testing for remaining cross-section correlation. + Journal of Business & Economic Statistics, 40(3), 1191-1203. """ + reason = self._pesaran_cd_inapplicable() + if reason is not None: + return InapplicableTestStatistic(reason=reason, name="Pesaran CD test") eps = self.idiosyncratic.iloc[:, 0] wide = eps.unstack(level=0) wide = wide.loc[:, wide.notnull().any(axis=0)] @@ -817,6 +876,31 @@ def wald_test( ) +def _is_concentrated_in_time(groups: Any, entity_ids: Any, time_ids: Any) -> bool: + """ + Test whether groups of entities change over time and are missing in periods + + Removing the mean of a group that is the same in every period for an + entity does not change the correlations over time. Removing the mean of a + group that changes over time and is not observed in every period removes a + component that varies over time and is shared by the entities in the group. + A group with a single entity is not shared. + """ + frame = DataFrame( + { + "group": np.asarray(groups).ravel(), + "entity": np.asarray(entity_ids).ravel(), + "time": np.asarray(time_ids).ravel(), + } + ) + if frame.groupby("entity")["group"].nunique().max() <= 1: + return False + by_group = frame.groupby("group") + shared = by_group["entity"].nunique() > 1 + nperiod = frame["time"].nunique() + return bool((by_group["time"].nunique()[shared] < nperiod).any()) + + class PanelEffectsResults(PanelResults): """ Results container for panel data models that include effects @@ -835,6 +919,28 @@ def __init__(self, res: AttrDict) -> None: self._r2_ex_effects = res.r2_ex_effects self._effects = res.effects + def _pesaran_cd_inapplicable(self) -> str | None: + reason = ( + "The Pesaran CD test is not applicable to models with {effect}. " + "These effects demean the residuals across entities that are " + "observed in the same periods, which makes the residuals " + "negatively correlated across entities and biases the statistic " + "downwards, away from its null distribution (Juodis and Reese, " + "2022)." + ) + if self._time_effect: + return reason.format(effect="time effects") + if self._other_effect: + other = self.model._other_effect_cats.dataframe + entity_ids = self.model.dependent.entity_ids + time_ids = self.model.dependent.time_ids + for col in other: + if _is_concentrated_in_time(other[col], entity_ids, time_ids): + return reason.format( + effect="other effects that are concentrated in some periods" + ) + return super()._pesaran_cd_inapplicable() + @property def f_pooled(self) -> WaldTestStatistic: r""" @@ -1274,6 +1380,17 @@ def __init__(self, res: AttrDict): self._avg_r2 = res.avg_r2 self._avg_adj_r2 = res.avg_adj_r2 + def _pesaran_cd_inapplicable(self) -> str | None: + return ( + "The Pesaran CD test is not applicable to FamaMacBeth. The errors " + "of the model are the residuals of the period-by-period " + "cross-sectional regressions. The series implied by the averaged " + "coefficients contains the common, time-varying term " + "(beta_t - beta_bar)'x_it, so it is cross-sectionally dependent by " + "construction. Fama-MacBeth inference is already robust to " + "cross-sectional dependence." + ) + @property def all_params(self) -> DataFrame: """ diff --git a/linearmodels/tests/panel/test_results.py b/linearmodels/tests/panel/test_results.py index 4ff225a3c2..cfce9ce4a2 100644 --- a/linearmodels/tests/panel/test_results.py +++ b/linearmodels/tests/panel/test_results.py @@ -13,8 +13,6 @@ from linearmodels.iv.model import IV2SLS from linearmodels.panel.data import PanelData from linearmodels.panel.model import ( - BetweenOLS, - FamaMacBeth, FirstDifferenceOLS, PanelOLS, PooledOLS, @@ -227,30 +225,18 @@ def test_pesaran_cd_plm_wage_panel(data): @pytest.mark.parametrize( - "estimator", - [ - "pooled", - "entity", - "entity-time", - "random", - "between", - "first-difference", - "fama-macbeth", - ], + "estimator", ["pooled", "entity", "random", "first-difference"] ) def test_pesaran_cd_all_estimators(estimator): + # The estimators where the test is not applicable are in + # test_results_pesaran_cd.py frame = pesaran_panel(**PESARAN_PANELS["unbalanced-dependent"]) exog = add_constant(frame[["x"]]) models = { "pooled": lambda: PooledOLS(frame.y, exog), "entity": lambda: PanelOLS(frame.y, frame[["x"]], entity_effects=True), - "entity-time": lambda: PanelOLS( - frame.y, frame[["x"]], entity_effects=True, time_effects=True - ), "random": lambda: RandomEffects(frame.y, exog), - "between": lambda: BetweenOLS(frame.y, exog), "first-difference": lambda: FirstDifferenceOLS(frame.y, frame[["x"]]), - "fama-macbeth": lambda: FamaMacBeth(frame.y, exog), } res = models[estimator]().fit() cd = res.pesaran_cd diff --git a/linearmodels/tests/panel/test_results_pesaran_cd.py b/linearmodels/tests/panel/test_results_pesaran_cd.py new file mode 100644 index 0000000000..a4961f6abb --- /dev/null +++ b/linearmodels/tests/panel/test_results_pesaran_cd.py @@ -0,0 +1,251 @@ +""" +Tests of PanelResults.pesaran_cd: which models the test is applicable to. +""" + +from itertools import combinations + +import numpy as np +from numpy.testing import assert_allclose +import pandas as pd +import pytest +from scipy import stats +from statsmodels.tools.tools import add_constant + +from linearmodels.panel import ( + BetweenOLS, + FamaMacBeth, + FirstDifferenceOLS, + PanelOLS, + PooledOLS, + RandomEffects, +) +from linearmodels.shared.hypotheses import ( + InapplicableTestStatistic, + NormalTestStatistic, +) + + +def cd_panel(seed=0, nentity=30, nperiod=8, loading=0.0): + """Panel with independent errors, or a common factor if loading > 0""" + rg = np.random.default_rng(seed) + factor = rg.standard_normal(nperiod) + load = 1.0 + 0.5 * rg.standard_normal(nentity) + x = rg.standard_normal((nentity, nperiod)) + effects = rg.standard_normal((nentity, 1)) + shock = loading * load[:, None] * factor + rg.standard_normal((nentity, nperiod)) + index = pd.MultiIndex.from_product( + [np.arange(nentity), np.arange(nperiod)], names=["entity", "time"] + ) + y = 1.0 + 0.5 * x + effects + shock + return pd.DataFrame({"y": y.ravel(), "x": x.ravel()}, index=index) + + +def brute_force_cd(idiosyncratic): + """Pesaran CD using explicit loops over entity pairs""" + resid = idiosyncratic.iloc[:, 0] + by_entity = {key: grp.droplevel(0) for key, grp in resid.groupby(level=0)} + total = 0.0 + npairs = 0 + for first, second in combinations(by_entity, 2): + common = by_entity[first].index.intersection(by_entity[second].index) + if len(common) < 2: + continue + rho = np.corrcoef(by_entity[first][common], by_entity[second][common])[0, 1] + if np.isfinite(rho): + total += np.sqrt(len(common)) * rho + npairs += 1 + return total / np.sqrt(npairs) + + +def cd_models(frame): + """Functions that create the models and if the Pesaran CD test applies""" + y, x = frame.y, frame[["x"]] + xc = add_constant(x) + entity = pd.Series(frame.index.get_level_values(0), index=frame.index) + period = pd.Series(frame.index.get_level_values(1), index=frame.index) + nperiod = int(period.max()) + 1 + nentity = int(entity.max()) + 1 + rg = np.random.default_rng(1234) + random_cats = pd.Series(rg.integers(0, 4, size=frame.shape[0]), index=frame.index) + weights = pd.Series(1.0 + (np.arange(frame.shape[0]) % 4), index=frame.index) + group_period = (entity % 3) * nperiod + period + group_block = (entity % 3) * nperiod + period // 2 + # Groups with a single entity are not shared with other entities + entity_regime = entity * 2 + (period >= nperiod // 2).astype(int) + # Groups that change over time but are present in every period + states = np.empty((nentity, nperiod), dtype=int) + states[:, 0] = rg.integers(0, 3, size=nentity) + for i in range(1, nperiod): + stay = rg.random(nentity) < 0.6 + states[:, i] = np.where(stay, states[:, i - 1], rg.integers(0, 3, size=nentity)) + markov = pd.Series(states.ravel(), index=frame.index) + # A treatment indicator that is switched on at a different date for each entity + switch = pd.Series(rg.integers(1, nperiod, size=nentity)) + post = (period >= switch.loc[entity].to_numpy()).astype(int) + return { + "pooled": (lambda: PooledOLS(y, xc), True), + "entity": (lambda: PanelOLS(y, x, entity_effects=True), True), + "entity-weighted": ( + lambda: PanelOLS(y, x, entity_effects=True, weights=weights), + True, + ), + "random": (lambda: RandomEffects(y, xc), True), + "first-difference": (lambda: FirstDifferenceOLS(y, x), True), + # Time-invariant groups of entities, and groups unrelated to the periods + "other-group": (lambda: PanelOLS(y, xc, other_effects=entity % 3), True), + "entity-other-random": ( + lambda: PanelOLS(y, x, entity_effects=True, other_effects=random_cats), + True, + ), + "entity-other-markov": ( + lambda: PanelOLS(y, x, entity_effects=True, other_effects=markov), + True, + ), + "other-entity-regime": ( + lambda: PanelOLS(y, x, other_effects=entity_regime), + True, + ), + "between": (lambda: BetweenOLS(y, xc), False), + "fama-macbeth": (lambda: FamaMacBeth(y, xc), False), + "time": (lambda: PanelOLS(y, x, time_effects=True), False), + "entity-time": ( + lambda: PanelOLS(y, x, entity_effects=True, time_effects=True), + False, + ), + # Time effects specified as other effects + "other-period": (lambda: PanelOLS(y, xc, other_effects=period), False), + "entity-other-period": ( + lambda: PanelOLS(y, x, entity_effects=True, other_effects=period), + False, + ), + "entity-other-group-period": ( + lambda: PanelOLS(y, x, entity_effects=True, other_effects=group_period), + False, + ), + "entity-other-blocks": ( + lambda: PanelOLS(y, x, entity_effects=True, other_effects=period // 2), + False, + ), + # Groups of entities by blocks of periods, which are not in every period + "entity-other-group-block": ( + lambda: PanelOLS(y, x, entity_effects=True, other_effects=group_block), + False, + ), + "entity-other-post": ( + lambda: PanelOLS(y, x, entity_effects=True, other_effects=post), + False, + ), + "other-random-period": ( + lambda: PanelOLS( + y, + xc, + other_effects=pd.DataFrame({"random": random_cats, "period": period}), + ), + False, + ), + } + + +CD_MODELS = list(cd_models(cd_panel())) + + +@pytest.mark.parametrize("name", CD_MODELS) +def test_pesaran_cd_applicability(name): + create, applicable = cd_models(cd_panel())[name] + res = create().fit() + cd = res.pesaran_cd + if applicable: + assert isinstance(cd, NormalTestStatistic) + expected = brute_force_cd(res.idiosyncratic) + assert_allclose(cd.stat, expected, rtol=1e-10) + assert_allclose(cd.pval, 2 * stats.norm.sf(abs(expected)), rtol=1e-10) + return + assert isinstance(cd, InapplicableTestStatistic) + assert np.isnan(cd.stat) + assert np.isnan(cd.pval) + assert cd.critical_values is None + assert "Pesaran CD test is not applicable" in str(cd) + assert "Pesaran CD test" in str(cd) + + +@pytest.mark.parametrize( + ("name", "phrase"), + [ + ("between", "BetweenOLS"), + ("between", "entity means"), + ("fama-macbeth", "FamaMacBeth"), + ("fama-macbeth", "(beta_t - beta_bar)'x_it"), + ("fama-macbeth", "robust to cross-sectional dependence"), + ("time", "time effects"), + ("entity-time", "time effects"), + ("other-period", "concentrated in some periods"), + ("entity-other-period", "concentrated in some periods"), + ("entity-other-group-period", "concentrated in some periods"), + ("entity-other-blocks", "concentrated in some periods"), + ("entity-other-group-block", "concentrated in some periods"), + ("entity-other-post", "concentrated in some periods"), + ("other-random-period", "concentrated in some periods"), + ("time", "demean the residuals across entities"), + ], +) +def test_pesaran_cd_inapplicable_reason(name, phrase): + create, applicable = cd_models(cd_panel())[name] + assert not applicable + cd = create().fit().pesaran_cd + assert isinstance(cd, InapplicableTestStatistic) + assert phrase in str(cd) + + +def test_pesaran_cd_time_effects_mechanism(): + # Time effects remove the cross-sectional mean of the residuals in every + # period, which is why the correlation of independent errors is negative + frame = cd_panel(seed=7, nentity=40, nperiod=6) + y, x = frame.y, frame[["x"]] + two_way = PanelOLS(y, x, entity_effects=True, time_effects=True).fit() + period_means = two_way.idiosyncratic.iloc[:, 0].groupby(level=1).mean() + assert_allclose(period_means, 0.0, atol=1e-12) + entity = PanelOLS(y, x, entity_effects=True).fit() + period_means = entity.idiosyncratic.iloc[:, 0].groupby(level=1).mean() + assert np.abs(period_means).max() > 1e-3 + # The correlation of the demeaned residuals is about -1/(N-1) for every pair + wide = two_way.idiosyncratic.iloc[:, 0].unstack(level=0) + corr = wide.corr().to_numpy() + avg = corr[np.tril_indices(corr.shape[0], k=-1)].mean() + assert avg < -0.5 / (frame.index.get_level_values(0).nunique() - 1) + + +def test_pesaran_cd_still_detects_dependence(): + # The applicable models still reject in the presence of a common factor + frame = cd_panel(seed=11, nentity=60, nperiod=10, loading=1.0) + y, x = frame.y, frame[["x"]] + res = PanelOLS(y, x, entity_effects=True).fit() + assert isinstance(res.pesaran_cd, NormalTestStatistic) + assert res.pesaran_cd.pval < 0.01 + + +def test_pesaran_cd_time_invariant_other_effects_unbalanced(): + # A group that is only observed in some periods does not remove a component + # that varies over time when the group of every entity is the same in + # every period + frame = cd_panel(seed=5, nentity=30, nperiod=8) + entity = frame.index.get_level_values(0) + period = frame.index.get_level_values(1) + # Entities 20 to 29 are only observed in the first 4 periods + frame = frame.loc[~((entity >= 20) & (period >= 4))] + region = pd.Series(frame.index.get_level_values(0) // 10, index=frame.index) + assert region.groupby(level=1).nunique().min() < 3 + res = PanelOLS(frame.y, add_constant(frame[["x"]]), other_effects=region).fit() + cd = res.pesaran_cd + assert isinstance(cd, NormalTestStatistic) + assert_allclose(cd.stat, brute_force_cd(res.idiosyncratic), rtol=1e-10) + + +def test_pesaran_cd_first_difference_uses_levels(): + # The shocks are in levels and not the residuals of the differenced equation + frame = cd_panel(seed=3, nentity=40, nperiod=8) + res = FirstDifferenceOLS(frame.y, frame[["x"]]).fit() + cd = res.pesaran_cd + assert isinstance(cd, NormalTestStatistic) + assert_allclose(cd.stat, brute_force_cd(res.idiosyncratic), rtol=1e-10) + differenced = brute_force_cd(res.resids.to_frame()) + assert not np.isclose(cd.stat, differenced, rtol=1e-3)