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..c85531cee9 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", @@ -603,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""" @@ -610,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:: @@ -630,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)] @@ -816,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 @@ -834,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""" @@ -965,6 +1072,79 @@ 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. + + 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 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`` + (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,46 +1169,47 @@ 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 ---------- 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. - include_constant : bool + 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. - 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 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. Raises ------ @@ -1037,8 +1218,9 @@ def wu_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 ----- @@ -1064,12 +1246,26 @@ def wu_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 ---------- @@ -1081,8 +1277,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.") @@ -1093,11 +1293,7 @@ def wu_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] @@ -1133,6 +1329,14 @@ def wu_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 @@ -1149,16 +1353,17 @@ def wu_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, ) - 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] @@ -1175,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/_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..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 @@ -488,7 +474,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 +484,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 +511,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 +529,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 +538,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 +549,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 +562,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 +572,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 +597,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 +608,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..3b8dfd5049 --- /dev/null +++ b/linearmodels/tests/panel/test_results_hausman.py @@ -0,0 +1,330 @@ +""" +Tests of RandomEffectsResults.hausman and HausmanTestResults: the return type and +the checks that the models are comparable. +""" + +import pickle +import warnings + +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, PooledOLS, RandomEffects +import linearmodels.panel.results +from linearmodels.panel.results import ( + HausmanTestResults, + PanelEffectsResults, + 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)) + + +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) 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)