Repository: statsmodels/statsmodels (python)

Issue:
`plot_pacf` default `lags` raises with a small sample
#### Describe the bug

When the parameter `lags` is not specified in `statsmodels.graphics.tsaplots.plot_pacf`, the value which is computed by default causes the function to raise a `ValueError` due to `nlags` being more than 50% of the sample size.

The limitation of 50% of the sample size should be applied when computing default `nlags`.

#### Code Sample, a copy-pastable example if possible


```python
import pandas as pd
from statsmodels.graphics.tsaplots import plot_pacf

df = df = pd.DataFrame(
    index=[pd.Timestamp.now() + pd.Timedelta(seconds=i) for i in range(10)],
    columns=["a"],
    data=list(range(10))
)
plot_pacf(df)
```
```traceback
File /tmp/venv/lib/python3.11/site-packages/statsmodels/graphics/tsaplots.py:353, in plot_pacf(x, ax, lags, alpha, method, use_vlines, title, zero, vlines_kwargs, **kwargs)
    351     acf_x = pacf(x, nlags=nlags, alpha=alpha, method=method)
    352 else:
--> 353     acf_x, confint = pacf(x, nlags=nlags, alpha=alpha, method=method)
    355 _plot_corr(
    356     ax,
    357     title,
   (...)
    364     **kwargs,
    365 )
    367 return fig

File /tmp/venv/lib/python3.11/site-packages/statsmodels/tsa/stattools.py:1009, in pacf(x, nlags, method, alpha)
   1007     nlags = min(int(10 * np.log10(nobs)), nobs // 2 - 1)
   1008 if nlags >= x.shape[0] // 2:
-> 1009     raise ValueError(
   1010         "Can only compute partial correlations for lags up to 50% of the "
   1011         f"sample size. The requested nlags {nlags} must be < "
   1012         f"{x.shape[0] // 2}."
   1013     )
   1015 if method in ("ols", "ols-inefficient", "ols-adjusted"):
   1016     efficient = "inefficient" not in method

ValueError: Can only compute partial correlations for lags up to 50% of the sample size. The requested nlags 9 must be < 5.
```
<details>

**Note**: As you can see, there are many issues on our GitHub tracker, so it is very possible that your issue has been posted before. Please check first before submitting so that we do not have to handle and close duplicates.

**Note**: Please be sure you are using the latest released version of `statsmodels`, or a recent build of `main`. If your problem has been fixed in an unreleased version, you might be able to use `main` until a new release occurs. 

**Note**: If you are using a released version, have you verified that the bug exists in the main branch of this repository? It helps the limited resources if we know problems exist in the current main branch so that they do not need to check whether the code sample produces a bug in the next release.

</details>


If the issue has not been resolved, please file it in the issue tracker.

#### Expected Output
The limitation of 50% of the sample size should be applied when computing default `nlags`.

See for example https://github.com/datamole-ai/edvart/pull/154 for a workaround.

#### Output of ``import statsmodels.api as sm; sm.show_versions()``

<details>

INSTALLED VERSIONS
------------------
Python: 3.11.3.final.0
OS: Linux 6.5.3-1-MANJARO #1 SMP PREEMPT_DYNAMIC Wed Sep 13 12:21:35 UTC 2023 x86_64
byteorder: little
LC_ALL: None
LANG: en_US.UTF-8

statsmodels
===========

Installed: 0.14.0 /tmp/venv/lib/python3.11/site-packages/statsmodels)

Required Dependencies
=====================

cython: Not installed
numpy: 1.24.4 (/tmp/venv/lib/python3.11/site-packages/numpy)
scipy: 1.10.1 (/tmp/venv/lib/python3.11/site-packages/scipy)
pandas: 1.5.3 (/tmp/venv/lib/python3.11/site-packages/pandas)
    dateutil: 2.8.2 (/tmp/venv/lib/python3.11/site-packages/dateutil)
patsy: 0.5.3 (/tmp/venv/lib/python3.11/site-packages/patsy)

Optional Dependencies
=====================

matplotlib: 3.7.3 (/tmp/venv/lib/python3.11/site-packages/matplotlib)
    backend: module://matplotlib_inline.backend_inline 
cvxopt: Not installed
joblib: 1.3.2 (/tmp/venv/lib/python3.11/site-packages/joblib)

Developer Tools
================

IPython: 8.12.2 (/tmp/venv/lib/python3.11/site-packages/IPython)
    jinja2: 3.1.2 (/tmp/venv/lib/python3.11/site-packages/jinja2)
sphinx: 7.1.2 (/tmp/venv/lib/python3.11/site-packages/sphinx)
    pygments: 2.16.1 (/tmp/venv/lib/python3.11/site-packages/pygments)
pytest: 7.1.3 (/tmp/venv/lib/python3.11/site-packages/pytest)
virtualenv: Not installed

</details>


Pull request description:
ENH: Improve lag selection in pacf
Reduce default lag size when sample is small

closes #9015

- [X] closes #9015
- [X] tests added / passed. 
- [X] code/documentation is well formatted.  
- [X] properly formatted commit message. See 
      [NumPy's guide](https://docs.scipy.org/doc/numpy-1.15.1/dev/gitwash/development_workflow.html#writing-the-commit-message). 

<details>


**Notes**:

* It is essential that you add a test when making code changes. Tests are not 
  needed for doc changes.
* When adding a new function, test values should usually be verified in another package (e.g., R/SAS/Stata).
* When fixing a bug, you must add a test that would produce the bug in main and
  then show that it is fixed with the new code.
* New code additions must be well formatted. Changes should pass flake8. If on Linux or OSX, you can
  verify you changes are well formatted by running 
  ```
  git diff upstream/main -u -- "*.py" | flake8 --diff --isolated
  ```
  assuming `flake8` is installed. This command is also available on Windows 
  using the Windows System for Linux once `flake8` is installed in the 
  local Linux environment. While passing this test is not required, it is good practice and it help 
  improve code quality in `statsmodels`.
* Docstring additions must render correctly, including escapes and LaTeX.

</details>


Patch:
diff --git a/statsmodels/graphics/tsaplots.py b/statsmodels/graphics/tsaplots.py
index ddeec5293..c3358cd28 100644
--- a/statsmodels/graphics/tsaplots.py
+++ b/statsmodels/graphics/tsaplots.py
@@ -17,7 +17,7 @@ def _prepare_data_corr_plot(x, lags, zero):
     if lags is None:
         # GH 4663 - use a sensible default value
         nobs = x.shape[0]
-        lim = min(int(np.ceil(10 * np.log10(nobs))), nobs - 1)
+        lim = min(int(np.ceil(10 * np.log10(nobs))), nobs // 2)
         lags = np.arange(not zero, lim + 1)
     elif np.isscalar(lags):
         lags = np.arange(not zero, int(lags) + 1)  # +1 for zero lag
diff --git a/statsmodels/tsa/stattools.py b/statsmodels/tsa/stattools.py
index 6b30f2173..891c77634 100644
--- a/statsmodels/tsa/stattools.py
+++ b/statsmodels/tsa/stattools.py
@@ -8,7 +8,7 @@ from statsmodels.compat.pandas import deprecate_kwarg
 from statsmodels.compat.python import Literal, lzip
 from statsmodels.compat.scipy import _next_regular
 
-from typing import Tuple
+from typing import Union, List
 import warnings
 
 import numpy as np
@@ -40,6 +40,8 @@ from statsmodels.tsa._innovations import innovations_algo, innovations_filter
 from statsmodels.tsa.adfvalues import mackinnoncrit, mackinnonp
 from statsmodels.tsa.tsatools import add_trend, lagmat, lagmat2ds
 
+ArrayLike1D = Union[np.ndarray, pd.Series, List[float]]
+
 __all__ = [
     "acovf",
     "acf",
@@ -709,8 +711,11 @@ def acf(
     else:
         return acf, qstat, pvalue
 
-
-def pacf_yw(x, nlags=None, method="adjusted"):
+def pacf_yw(
+    x: ArrayLike1D,
+    nlags: int | None = None,
+    method: Literal["adjusted", "mle"] = "adjusted",
+) -> np.ndarray:
     """
     Partial autocorrelation estimated with non-recursive yule_walker.
 
@@ -747,8 +752,7 @@ def pacf_yw(x, nlags=None, method="adjusted"):
     nlags = int_like(nlags, "nlags", optional=True)
     nobs = x.shape[0]
     if nlags is None:
-        nlags = min(int(10 * np.log10(nobs)), nobs - 1)
-
+        nlags = max(min(int(10 * np.log10(nobs)), nobs - 1), 1)
     method = string_like(method, "method", options=("adjusted", "mle"))
     pacf = [1.0]
     with warnings.catch_warnings():
@@ -758,7 +762,9 @@ def pacf_yw(x, nlags=None, method="adjusted"):
     return np.array(pacf)
 
 
-def pacf_burg(x, nlags=None, demean=True):
+def pacf_burg(
+    x: ArrayLike1D, nlags: int | None = None, demean: bool = True
+) -> tuple[np.ndarray, np.ndarray]:
     """
     Calculate Burg"s partial autocorrelation estimator.
 
@@ -800,6 +806,7 @@ def pacf_burg(x, nlags=None, demean=True):
         x = x - x.mean()
     nobs = x.shape[0]
     p = nlags if nlags is not None else min(int(10 * np.log10(nobs)), nobs - 1)
+    p = max(p, 1)
     if p > nobs - 1:
         raise ValueError("nlags must be smaller than nobs - 1")
     d = np.zeros(p + 1)
@@ -818,14 +825,19 @@ def pacf_burg(x, nlags=None, demean=True):
         v[1:] = last_v[1:] - pacf[i] * last_u[:-1]
         d[i + 1] = (1 - pacf[i] ** 2) * d[i] - v[i] ** 2 - u[-1] ** 2
         pacf[i + 1] = 2 / d[i + 1] * v[i + 1 :].dot(u[i:-1])
-    sigma2 = (1 - pacf ** 2) * d / (2.0 * (nobs - np.arange(0, p + 1)))
+    sigma2 = (1 - pacf**2) * d / (2.0 * (nobs - np.arange(0, p + 1)))
     pacf[0] = 1  # Insert the 0 lag partial autocorrel
 
     return pacf, sigma2
 
 
 @deprecate_kwarg("unbiased", "adjusted")
-def pacf_ols(x, nlags=None, efficient=True, adjusted=False):
+def pacf_ols(
+    x: ArrayLike1D,
+    nlags: int | None = None,
+    efficient: bool = True,
+    adjusted: bool = False,
+) -> np.ndarray:
     """
     Calculate partial autocorrelations via OLS.
 
@@ -885,8 +897,9 @@ def pacf_ols(x, nlags=None, efficient=True, adjusted=False):
     adjusted = bool_like(adjusted, "adjusted")
     nobs = x.shape[0]
     if nlags is None:
-        nlags = min(int(10 * np.log10(nobs)), nobs - 1)
-
+        nlags = max(min(int(10 * np.log10(nobs)), nobs // 2), 1)
+    if nlags > nobs//2:
+        raise ValueError(f"nlags must be smaller than nobs // 2 ({nobs//2})")
     pacf = np.empty(nlags + 1)
     pacf[0] = 1.0
     if efficient:
@@ -903,14 +916,30 @@ def pacf_ols(x, nlags=None, efficient=True, adjusted=False):
             params = lstsq(xlags[:, :k], x0, rcond=None)[0]
             # Last coefficient corresponds to PACF value (see [1])
             pacf[k] = np.squeeze(params[-1])
-
     if adjusted:
         pacf *= nobs / (nobs - np.arange(nlags + 1))
-
     return pacf
 
 
-def pacf(x, nlags=None, method="ywadjusted", alpha=None):
+def pacf(
+    x: ArrayLike1D,
+    nlags: int | None = None,
+    method: Literal[
+        "yw",
+        "ywadjusted",
+        "ols",
+        "ols-inefficient",
+        "ols-adjusted",
+        "ywm",
+        "ywmle",
+        "ld",
+        "ldadjusted",
+        "ldb",
+        "ldbiased",
+        "burg",
+    ] = "ywadjusted",
+    alpha: float | None = None,
+) -> np.ndarray | tuple[np.ndarray, np.ndarray]:
     """
     Partial autocorrelation estimate.
 
@@ -996,7 +1025,7 @@ def pacf(x, nlags=None, method="ywadjusted", alpha=None):
         "ldb",
         "ldbiased",
         "ld_biased",
-        "burg"
+        "burg",
     )
     x = array_like(x, "x", maxdim=2)
     method = string_like(method, "method", options=methods)
@@ -1005,13 +1034,13 @@ def pacf(x, nlags=None, method="ywadjusted", alpha=None):
     nobs = x.shape[0]
     if nlags is None:
         nlags = min(int(10 * np.log10(nobs)), nobs // 2 - 1)
-    if nlags >= x.shape[0] // 2:
+    nlags = max(nlags, 1)
+    if nlags > x.shape[0] // 2:
         raise ValueError(
             "Can only compute partial correlations for lags up to 50% of the "
             f"sample size. The requested nlags {nlags} must be < "
             f"{x.shape[0] // 2}."
         )
-
     if method in ("ols", "ols-inefficient", "ols-adjusted"):
         efficient = "inefficient" not in method
         adjusted = "adjusted" in method
@@ -1031,7 +1060,6 @@ def pacf(x, nlags=None, method="ywadjusted", alpha=None):
         acv = acovf(x, adjusted=False, fft=False)
         ld_ = levinson_durbin(acv, nlags=nlags, isacov=True)
         ret = ld_[2]
-
     if alpha is not None:
         varacf = 1.0 / len(x)  # for all lags >=1
         interval = stats.norm.ppf(1.0 - alpha / 2.0) * np.sqrt(varacf)
@@ -1538,7 +1566,7 @@ def grangercausalitytests(x, maxlag, addconst=True, verbose=None):
     if x.shape[0] <= 3 * maxlag + int(addconst):
         raise ValueError(
             "Insufficient observations. Maximum allowable "
-            "lag is {0}".format(int((x.shape[0] - int(addconst)) / 3) - 1)
+            "lag is {}".format(int((x.shape[0] - int(addconst)) / 3) - 1)
         )
 
     resli = {}
@@ -1945,7 +1973,7 @@ def kpss(
     regression: Literal["c", "ct"] = "c",
     nlags: Literal["auto", "legacy"] | int = "auto",
     store: bool = False,
-) -> Tuple[float, float, int, dict[str, float]]:
+) -> tuple[float, float, int, dict[str, float]]:
     """
     Kwiatkowski-Phillips-Schmidt-Shin test for stationarity.
 
@@ -2032,7 +2060,7 @@ def kpss(
 
     # if m is not one, n != m * n
     if nobs != x.size:
-        raise ValueError("x of shape {0} not understood".format(x.shape))
+        raise ValueError(f"x of shape {x.shape} not understood")
 
     if hypo == "ct":
         # p. 162 Kwiatkowski et al. (1992): y_t = beta * t + r_t + e_t,
@@ -2106,8 +2134,8 @@ look-up table. The actual p-value is {direction} than the p-value returned.
         rstore.nobs = nobs
 
         stationary_type = "level" if hypo == "c" else "trend"
-        rstore.H0 = "The series is {0} stationary".format(stationary_type)
-        rstore.HA = "The series is not {0} stationary".format(stationary_type)
+        rstore.H0 = f"The series is {stationary_type} stationary"
+        rstore.HA = f"The series is not {stationary_type} stationary"
 
         return kpss_stat, p_value, crit_dict, rstore
     else:
@@ -2200,7 +2228,7 @@ def range_unit_root_test(x, store=False):
 
     # if m is not one, n != m * n
     if nobs != x.size:
-        raise ValueError("x of shape {0} not understood".format(x.shape))
+        raise ValueError(f"x of shape {x.shape} not understood")
 
     # Table from [1] has been replicated using 200,000 samples
     # Critical values for new n_obs values have been identified


Test patch:
diff --git a/statsmodels/graphics/tests/test_tsaplots.py b/statsmodels/graphics/tests/test_tsaplots.py
index 44c1b49d0..5ba24f02a 100644
--- a/statsmodels/graphics/tests/test_tsaplots.py
+++ b/statsmodels/graphics/tests/test_tsaplots.py
@@ -375,3 +375,14 @@ def test_predict_plot(use_pandas, model_and_args, alpha):
     res = model(y, **kwargs).fit()
     fig = plot_predict(res, start, end, alpha=alpha)
     assert isinstance(fig, plt.Figure)
+
+
+@pytest.mark.matplotlib
+def test_plot_pacf_small_sample():
+    idx = [pd.Timestamp.now() + pd.Timedelta(seconds=i) for i in range(10)]
+    df = pd.DataFrame(
+        index=idx,
+        columns=["a"],
+        data=list(range(10))
+    )
+    plot_pacf(df)
diff --git a/statsmodels/tsa/tests/test_stattools.py b/statsmodels/tsa/tests/test_stattools.py
index 720e2bedb..abbc07181 100644
--- a/statsmodels/tsa/tests/test_stattools.py
+++ b/statsmodels/tsa/tests/test_stattools.py
@@ -1533,7 +1533,7 @@ def test_acf_conservate_nanops(reset_randomstate):
 
 
 def test_pacf_nlags_error(reset_randomstate):
-    e = np.random.standard_normal(100)
+    e = np.random.standard_normal(99)
     with pytest.raises(ValueError, match="Can only compute partial"):
         pacf(e, 50)
 
@@ -1584,3 +1584,27 @@ def test_granger_causality_exception_maxlag(gc_data):
 def test_granger_causality_verbose(gc_data):
     with pytest.warns(FutureWarning, match="verbose"):
         grangercausalitytests(gc_data, 3, verbose=True)
+
+@pytest.mark.parametrize("size",[3,5,7,9])
+def test_pacf_small_sample(size,reset_randomstate):
+    y = np.random.standard_normal(size)
+    a = pacf(y)
+    assert isinstance(a, np.ndarray)
+    a, b = pacf_burg(y)
+    assert isinstance(a, np.ndarray)
+    assert isinstance(b, np.ndarray)
+    a = pacf_ols(y)
+    assert isinstance(a, np.ndarray)
+    a = pacf_yw(y)
+    assert isinstance(a, np.ndarray)
+
+
+def test_pacf_1_obs(reset_randomstate):
+    y = np.random.standard_normal(1)
+    with pytest.raises(ValueError):
+        pacf(y)
+    with pytest.raises(ValueError):
+        pacf_burg(y)
+    with pytest.raises(ValueError):
+        pacf_ols(y)
+    pacf_yw(y)