Merge pull request #339 from bashtage/sur-breusch-pagan

ENH: Add Breusch-Pagan test

Author: bashtage

linearmodels/system/results.py

@@ -11,7 +11,7 @@
 
 import linearmodels
 from linearmodels.shared.base import _SummaryStr
-from linearmodels.shared.hypotheses import WaldTestStatistic
+from linearmodels.shared.hypotheses import InvalidTestStatistic, WaldTestStatistic
 from linearmodels.shared.io import _str, format_wide, param_table, pval_format
 from linearmodels.shared.utility import AttrDict
 from linearmodels.typing import ArrayLike, NDArray
@@ -516,6 +516,51 @@ def summary(self) -> Summary:
 
         return smry
 
+    def breusch_pagan(self) -> Union[WaldTestStatistic, InvalidTestStatistic]:
+        r"""
+        Breusch-Pagan LM test for no cross-correlation
+
+        Returns
+        -------
+        WaldTestStatistic
+            Test statistic for null all correlationsare zero.
+
+        Notes
+        -----
+        The null hypothesis is that the shock correlations are all 0, and
+        so there are no gains to using GLS estimation in the system estimator.
+        When the null is rejected, there should be efficiency gains to using
+        GLS as long the regressors are not common to all models.
+
+        The Breusch-Pagan test statistic is defined as
+
+        .. math::
+
+           LM = n \sum_{i=1}^k \sum_{j=i+1}^k \hat{\rho}_{ij}^2
+
+        where :math:`\hat{\rho}_{ij}` is the sample residual correlation
+        between series i and j. n is the sample size. It has an asymptotic
+        :math:`\chi^2_{k(k-1)/2}` distribution.
+        """
+        name = "Breusch-Pagan LM Test"
+        resids = self.resids
+        if resids.shape[1] == 1:
+            return InvalidTestStatistic(
+                "Cannot test correlation when the system contains a single "
+                "dependent variable.",
+                name=name,
+            )
+        r = np.corrcoef(resids.T)
+        k = r.shape[0]
+        distinct_corr = np.tril(r, -1)
+        stat = self.resids.shape[0] * (distinct_corr ** 2).sum()
+        return WaldTestStatistic(
+            stat,
+            "Residuals are uncorrelated",
+            k * (k - 1) // 2,
+            name=name,
+        )
+
 
 class SystemEquationResult(_CommonResults):
     """

linearmodels/tests/panel/test_model.py

@@ -2,7 +2,7 @@
 import pickle
 
 import numpy as np
-from numpy.testing import assert_allclose, assert_equal
+from numpy.testing import assert_allclose
 import pandas as pd
 import pytest
 

linearmodels/tests/system/test_sur.py

@@ -7,9 +7,10 @@
 from pandas import DataFrame, Series, concat
 from pandas.testing import assert_frame_equal, assert_series_equal
 import pytest
+import scipy.stats
 
 from linearmodels.iv.model import _OLS as OLS
-from linearmodels.shared.hypotheses import InvalidTestStatistic
+from linearmodels.shared.hypotheses import InvalidTestStatistic, WaldTestStatistic
 from linearmodels.shared.utility import AttrDict
 from linearmodels.system._utility import (
     blocked_column_product,
@@ -851,3 +852,29 @@ def test_tvalues_homogeneity(method, cov_type):
     if cov_type == "robust" and method == "gls":
         assert_allclose(res0.tstats, base_tstat)
         assert_allclose(res1.tstats, base_100_tstat)
+
+
+@pytest.mark.parametrize("k", [1, 3])
+def test_brequsch_pagan(k):
+    eqns = generate_data(k=k)
+    mod = SUR(eqns)
+    res = mod.fit()
+    stat = res.breusch_pagan()
+    if k == 1:
+        assert isinstance(stat, InvalidTestStatistic)
+        assert "Breusch-Pagan" in str(stat)
+        assert np.isnan(stat.stat)
+        return
+    rho = np.asarray(res.resids.corr())
+    nobs = res.resids.shape[0]
+    direct = 0.0
+    for i in range(3):
+        for j in range(i + 1, 3):
+            direct += rho[i, j] ** 2
+    direct *= nobs
+    assert isinstance(stat, WaldTestStatistic)
+    assert_allclose(stat.stat, direct)
+    assert stat.df == 3
+    assert_allclose(stat.pval, 1.0 - scipy.stats.chi2(3).cdf(direct))
+    assert "Residuals are uncorrelated" in stat.null
+    assert "Breusch-Pagan" in str(stat)