Merge 6ae72fc into e71a2fb

arch/tests/unitroot/test_unitroot.py

@@ -25,6 +25,7 @@
     mackinnoncrit,
     mackinnonp,
 )
+from arch.utility.exceptions import InfeasibleTestException
 
 DECIMAL_5 = 5
 DECIMAL_4 = 4
@@ -474,8 +475,40 @@ def test_adf_short_timeseries():
     # GH 262
     x = np.asarray([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0])
     adf = ADF(x)
-    assert_almost_equal(adf.stat, -2.236, decimal=3)
-    assert adf.lags == 1
+    with pytest.raises(InfeasibleTestException):
+        adf.stat
+
+
+def test_adf_buggy_timeseries1():
+    x = np.asarray([0])
+    adf = ADF(x)
+    # ValueError: maxlag should be < nobs
+    with pytest.raises(InfeasibleTestException, match="Fewer observations are"):
+        adf.stat
+
+
+def test_adf_buggy_timeseries2():
+    x = np.asarray([0, 0])
+    adf = ADF(x)
+    # IndexError: index 0 is out of bounds for axis 0 with size 0
+    with pytest.raises(InfeasibleTestException, match="Fewer observations are"):
+        adf.stat
+
+
+def test_adf_buggy_timeseries3():
+    x = np.asarray([1] * 1000)
+    adf = ADF(x)
+    # AssertionError: Number of manager items must equal union of block items
+    # # manager items: 1, # tot_items: 0
+    with pytest.raises(InfeasibleTestException, match="The maximum lag you are"):
+        adf.stat
+
+
+def test_kpss_buggy_timeseries1():
+    x = np.asarray([0])
+    adf = KPSS(x)
+    with pytest.raises(InfeasibleTestException):
+        adf.stat
 
 
 kpss_autolag_data = (
@@ -604,3 +637,36 @@ def test_invalid_trend():
 def test_nc_warning():
     with pytest.warns(FutureWarning, match='Trend "nc" is deprecated'):
         ADF(np.random.standard_normal(100), trend="nc")
+
+
+@pytest.mark.parametrize("nobs", np.arange(1, 11).tolist())
+@pytest.mark.parametrize(
+    "stat", [ADF, KPSS, DFGLS, PhillipsPerron, VarianceRatio, ZivotAndrews]
+)
+@pytest.mark.parametrize("trend", ["nc", "c", "ct", "ctt"])
+def test_wrong_exceptions(stat, nobs, trend):
+    y = np.random.standard_normal((nobs,))
+    try:
+        stat(y, trend=trend).stat
+    except InfeasibleTestException:
+        pass
+    except ValueError as ve:
+        if "trend" not in ve.args[0]:
+            raise
+
+
+@pytest.mark.parametrize("nobs", [2, 10, 100])
+@pytest.mark.parametrize(
+    "stat", [ADF, KPSS, DFGLS, PhillipsPerron, VarianceRatio, ZivotAndrews]
+)
+@pytest.mark.parametrize("trend", ["nc", "c", "ct", "ctt"])
+def test_wrong_exceptions_nearly_constant_series(stat, nobs, trend):
+    y = np.zeros((nobs,))
+    y[-1] = 1.0
+    try:
+        stat(y, trend=trend).stat
+    except InfeasibleTestException:
+        pass
+    except ValueError as ve:
+        if "trend" not in ve.args[0]:
+            raise

arch/unitroot/unitroot.py

@@ -3,7 +3,9 @@
 
 from numpy import (
     abs,
+    amax,
     amin,
+    any as npany,
     arange,
     argwhere,
     array,
@@ -19,6 +21,7 @@
     int32,
     int64,
     interp,
+    isnan,
     log,
     nan,
     ones,
@@ -30,7 +33,7 @@
     squeeze,
     sum,
 )
-from numpy.linalg import inv, matrix_rank, pinv, qr, solve
+from numpy.linalg import LinAlgError, inv, matrix_rank, pinv, qr, solve
 from pandas import DataFrame
 from scipy.stats import norm
 from statsmodels.iolib.summary import Summary
@@ -65,7 +68,12 @@
 from arch.unitroot.critical_values.zivot_andrews import za_critical_values
 from arch.utility import cov_nw
 from arch.utility.array import DocStringInheritor, ensure1d, ensure2d
-from arch.utility.exceptions import InvalidLengthWarning, invalid_length_doc
+from arch.utility.exceptions import (
+    InfeasibleTestException,
+    InvalidLengthWarning,
+    infeasible_test_error,
+    invalid_length_doc,
+)
 from arch.utility.timeseries import add_trend
 
 __all__ = [
@@ -108,6 +116,20 @@
 }
 
 
+def _is_reduced_rank(x: NDArray) -> bool:
+    """
+    Check if a matrix has reduced rank preferring quick checks
+    """
+    if x.shape[1] > x.shape[0]:
+        return True
+    elif npany(isnan(x)):
+        return True
+    elif sum(amax(x, axis=0) == amin(x, axis=0)) > 1:
+        return True
+    else:
+        return matrix_rank(x) < x.shape[1]
+
+
 def _select_best_ic(
     method: str, nobs: float, sigma2: NDArray, tstat: NDArray
 ) -> Tuple[float, int]:
@@ -357,6 +379,15 @@ def _df_select_lags(
         full_rhs = rhs
 
     start_lag = full_rhs.shape[1] - rhs.shape[1] + 1
+    if _is_reduced_rank(full_rhs):
+        raise InfeasibleTestException(
+            "The maximum lag you are considering results in an ADF regression with "
+            "a singular regressor matrix and so a specification test be run. This "
+            "may occur if your series have little variation and so is locally "
+            "constant, or may occur if you are attempting to test a very short "
+            "series. You can manually set maximum lag length to consider smaller "
+            "models."
+        )
     ic_best, best_lag = _autolag_ols(lhs, full_rhs, start_lag, max_lags, method)
     return ic_best, best_lag
 
@@ -443,6 +474,12 @@ def _repr_html_(self) -> str:
         """Display as HTML for IPython notebook."""
         return self.summary().as_html()
 
+    def _check_specification(self) -> None:
+        """
+        Check that the data are compatible with running a test.
+        """
+        raise NotImplementedError
+
     def _compute_statistic(self) -> None:
         """This is the core routine that computes the test statistic, computes
         the p-value and constructs the critical values.
@@ -460,6 +497,7 @@ def _compute_if_needed(self) -> None:
         needed
         """
         if self._stat is None:
+            self._check_specification()
             self._compute_statistic()
 
     @property
@@ -718,6 +756,12 @@ def _select_lag(self) -> None:
         self._ic_best = ic_best
         self._lags = best_lag
 
+    def _check_specification(self) -> None:
+        trend_order = len(self._trend) if self._trend != "nc" else 0
+        lag_len = 0 if self._lags is None else self._lags
+        if self._y.shape[0] < (3 + trend_order + lag_len):
+            raise InfeasibleTestException(infeasible_test_error)
+
     def _compute_statistic(self) -> None:
         if self._lags is None:
             self._select_lag()
@@ -845,6 +889,12 @@ def __init__(
         else:
             self._c = -13.5
 
+    def _check_specification(self) -> None:
+        trend_order = len(self._trend) if self._trend != "nc" else 0
+        lag_len = 0 if self._lags is None else self._lags
+        if self._y.shape[0] < (3 + trend_order + lag_len):
+            raise InfeasibleTestException(infeasible_test_error)
+
     def _compute_statistic(self) -> None:
         """Core routine to estimate DF-GLS test statistic"""
         # 1. GLS detrend
@@ -1028,6 +1078,12 @@ def __init__(
         self._test_name = "Phillips-Perron Test"
         self._lags = lags
 
+    def _check_specification(self) -> None:
+        trend_order = len(self._trend) if self._trend != "nc" else 0
+        lag_len = 0 if self._lags is None else self._lags
+        if self._y.shape[0] < (3 + trend_order + lag_len):
+            raise InfeasibleTestException(infeasible_test_error)
+
     def _compute_statistic(self) -> None:
         """Core routine to estimate PP test statistics"""
         # 1. Estimate Regression
@@ -1043,10 +1099,13 @@ def _compute_statistic(self) -> None:
         lhs = y[1:, None]
         if trend != "n":
             rhs = add_trend(rhs, trend)
-
+        if _is_reduced_rank(rhs):
+            raise InfeasibleTestException(infeasible_test_error)
         resols = OLS(lhs, rhs).fit()
         k = rhs.shape[1]
         n, u = resols.nobs, resols.resid
+        if lags > u.shape[0]:
+            raise InfeasibleTestException(infeasible_test_error)
         lam2 = cov_nw(u, lags, demean=False)
         lam = sqrt(lam2)
         # 2. Compute components
@@ -1186,6 +1245,12 @@ def __init__(
         self._alternative_hypothesis = "The process contains a unit root."
         self._resids: Optional[ArrayLike1D] = None
 
+    def _check_specification(self) -> None:
+        trend_order = len(self._trend) if self._trend != "nc" else 0
+        lag_len = 0 if self._lags is None else self._lags
+        if self._y.shape[0] <= (3 + trend_order + lag_len):
+            raise InfeasibleTestException(infeasible_test_error)
+
     def _compute_statistic(self) -> None:
         # 1. Estimate model with trend
         nobs, y, trend = self._nobs, self._y, self._trend
@@ -1199,6 +1264,8 @@ def _compute_statistic(self) -> None:
             else:
                 self._autolag()
         assert self._lags is not None
+        if u.shape[0] < self._lags + 1:
+            raise InfeasibleTestException(infeasible_test_error)
         lam = cov_nw(u, self._lags, demean=False)
         s = cumsum(u)
         self._stat = 1 / (nobs ** 2.0) * sum(s ** 2.0) / lam
@@ -1326,14 +1393,23 @@ def _quick_ols(endog: NDArray, exog: NDArray) -> NDArray:
         """
         Minimal implementation of LS estimator for internal use
         """
-        xpxi = inv(exog.T.dot(exog))
+        try:
+            xpxi = inv(exog.T.dot(exog))
+        except LinAlgError:
+            raise InfeasibleTestException(infeasible_test_error)
         xpy = exog.T.dot(endog)
         nobs, k_exog = exog.shape
         b = xpxi.dot(xpy)
         e = endog - exog.dot(b)
         sigma2 = e.T.dot(e) / (nobs - k_exog)
         return b / sqrt(diag(sigma2 * xpxi))
 
+    def _check_specification(self) -> None:
+        trend_order = len(self._trend) if self._trend != "nc" else 0
+        lag_len = 0 if self._lags is None else self._lags
+        if self._y.shape[0] < (3 + trend_order + lag_len):
+            raise InfeasibleTestException(infeasible_test_error)
+
     def _compute_statistic(self) -> None:
         """This is the core routine that computes the test statistic, computes
         the p-value and constructs the critical values.
@@ -1395,7 +1471,7 @@ def _compute_statistic(self) -> None:
             if bp == start_period + 1:
                 rank = matrix_rank(exog)
                 if rank < exog.shape[1]:
-                    raise ValueError(
+                    raise InfeasibleTestException(
                         "ZA: auxiliary exog matrix is not full rank.\n cols "
                         "{0}. rank = {1}".format(exog.shape[1], rank)
                     )
@@ -1549,6 +1625,12 @@ def debiased(self, value: bool) -> None:
         self._reset()
         self._debiased = bool(value)
 
+    def _check_specification(self) -> None:
+        trend_order = len(self._trend) if self._trend != "nc" else 0
+        lag_len = 0 if self._lags is None else self._lags
+        if self._y.shape[0] < (2 + trend_order + lag_len):
+            raise InfeasibleTestException(infeasible_test_error)
+
     def _compute_statistic(self) -> None:
         overlap, debiased, robust = self._overlap, self._debiased, self._robust
         y, nobs, q, trend = self._y, self._nobs, self._lags, self._trend

arch/utility/exceptions.py

@@ -80,3 +80,12 @@ class StudentizationError(RuntimeError):
 by func have not variability when resampled. The estimated covariance
 is:\n\n {cov}
 """
+
+
+class InfeasibleTestException(RuntimeError):
+    pass
+
+
+infeasible_test_error = """
+Fewer observations are available than are required to perform the test.
+"""

doc/source/changes/4.0.txt

@@ -4,6 +4,8 @@ Version 4
 
 Since 4.12
 ==========
+- Improved checks for unit root tests to help users understand when a test is infeasible
+  (:issue:`364`).
 - Added Phillips-Ouliaris (:func:`~arch.unitroot.cointegration.phillips_ouliaris`)
   cointegration tests (:issue:`360`).
 - Added three methods to estimate cointegrating vectors: