Repository: statsmodels/statsmodels
Pull request: BUG: Fix auto lag selection in acorr_ljungbox #8338
- [x] closes #8338
- [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). 

Corrected the implementation of optimal lag selection (Escanciano & Lobato, J. of Econometrics, 2009). `auto_lag = True` does not return a zero-size array ValueError when data is a non-autocorrelated series. Added test and corrected the previous test for this which did not flag `auto_lag = True` and therefore did not test this feature. Need additional test to compare to R results.

<details>

The `get_optimal_length` inner function of the `acorr_ljungbox` function calculated the optimal lag incorrectly. First, the optimal lag was initialized to `0`, which then caused an empty array to be initialized. The function barfed when the `max()` of this array was attempted, returning a "ValueError: zero-size array to reduction operation maximum which has no identity". Second, the least penalized LjungBox value was initialized at `0` and compared to subsequent values in a loop. However, the least penalized LB value can be less than zero, especially when the data series displays low autocorrelation. Thus, optimal lags were incorrectly calculated.

This fix eliminates `get_optimal_length` and instead performs an array-based calculation of optimal lag length. Returns a minimum lag of 1, as described in [1]. Further, it supports optimal lag calculation for low autocorrelation series (for e.g., a white noise process).

- No changes made to docstrings; some changes to inline comments
- New test for `auto_lag = True` when data is white noise process to replicate the problem and show resolution
- **IMPORTANT**: Need test to compare to results in R/Stata for additional confirmation

_[1] J. Carlos Escanciano, Ignacio N. Lobato “An automatic Portmanteau test for serial correlation”., J. of Econometrics Volume 151, 2009._

</details>

Problem statement:
auto_lag selection in acorr_ljungbox fails when there is no autocorrelation
#### Describe the bug

The `acorr_ljungbox` function has an `auto_lag` flag which determines the optimal lag for performing the Ljung-Box autocorrelations test. When there is no autocorrelation in the data, however, the function fails and returns a "ValueError: zero-size array to reduction operation maximum which has no identity". When there is low autocorrelation, it may be returning an incorrect optimal lag (for theoretical reasons - described under **Details**).

#### Code Sample, a copy-pastable example if possible

```python
import numpy as np
from statsmodels.stats.diagnostic import acorr_ljungbox

data = np.random.rand(1000) # white noise process with no autocorrelation
res = acorr_ljunbox(data, auto_lag=True)

ValueError: zero-size array to reduction operation maximum which has no identity
```
<details>

The `auto_lag` feature attempts to implement the lag selection criterion outlined in Escanciano & Lobato (J. of Econometrics, 2009). However, there are two errors in the implementation, which is performed by the `get_optimal_length` inner function of the` acorr_ljungbox` function. There is also one error in the associated pytest.

The first error is that the optimal_lag is initialized at `0`. When the rest of the routine does not update the optimal_lag (see second error below for why this may happen), the optimal lag length is set at `0`. This causes an empty array to be initialized, and when the function tries to find the `max()` of this array, the `ValueError` related to a zero-size array is returned. Even theoretically, [1] defines the putative optimal lag $\tilde{p}$ as `>=1`.

The second error is that `get_optimal_length` initializes the value of the least penalized test statistic ($L_p$ in eq. (4), pg. 143 in [1]) as `0`. This is incorrect --- there is no theoretical reason why the penalized value of the test statistic should be positive. This is especially true for series that exhibit low degrees of autocorrelation, when the unpenalized Box-Pierce (or Ljung-Box) test statistic can be a very low value. Thus, even if the least penalized test statistic occurs at a large lag, if it has a negative penalized score, `get_optimal_lag` would still return `0` as the optimal lag. 

Finally, the test associated with this function `test_ljungbox_auto_lag_selection()` (in `test_diagnostic.py`) calls `acorr_ljungbox()`, but omits to set `auto_lag = True`. So the test is invalid. 

I've initiated a pull request fixing the above issues and the test. This is my first bug report and pull request, so any guidance is greatly appreciated!

_[1] J. Carlos Escanciano, Ignacio N. Lobato “An automatic Portmanteau test for serial correlation”., J. of Econometrics Volume 151, 2009._

</details>


#### Expected Output

Return the result for the Ljung-Box test with the optimal lag as 1 if autocorrelation is absent, and the correct lag based on the Escanciano & Lobato procedure if autocorrelation is not absent.

#### Output of ``import statsmodels.api as sm; sm.show_versions()``

<details>

Output of ``import statsmodels.api as sm; sm.show_versions()`` 

INSTALLED VERSIONS
------------------
Python: 3.8.1.final.0

statsmodels
===========

Installed: 0.13.2 (C:\Users\Vasudeva\AppData\Local\Programs\Python\Python38\lib\site-packages\statsmodels)

Required Dependencies
=====================

cython: 0.29.15 (C:\Users\Vasudeva\AppData\Local\Programs\Python\Python38\lib\site-packages\Cython)
numpy: 1.18.1 (C:\Users\Vasudeva\AppData\Local\Programs\Python\Python38\lib\site-packages\numpy)
scipy: 1.4.1 (C:\Users\Vasudeva\AppData\Local\Programs\Python\Python38\lib\site-packages\scipy)
pandas: 1.0.1 (C:\Users\Vasudeva\AppData\Local\Programs\Python\Python38\lib\site-packages\pandas)
    dateutil: 2.8.1 (C:\Users\Vasudeva\AppData\Local\Programs\Python\Python38\lib\site-packages\dateutil)
patsy: 0.5.2 (C:\Users\Vasudeva\AppData\Local\Programs\Python\Python38\lib\site-packages\patsy)

Optional Dependencies
=====================

matplotlib: 3.1.3 (C:\Users\Vasudeva\AppData\Local\Programs\Python\Python38\lib\site-packages\matplotlib)
    backend: TkAgg
cvxopt: Not installed
joblib: 1.1.0 (C:\Users\Vasudeva\AppData\Local\Programs\Python\Python38\lib\site-packages\joblib)

Developer Tools
================

IPython: 7.16.1 (C:\Users\Vasudeva\AppData\Local\Programs\Python\Python38\lib\site-packages\IPython)
    jinja2: 2.11.2 (C:\Users\Vasudeva\AppData\Local\Programs\Python\Python38\lib\site-packages\jinja2)
sphinx: 3.1.2 (C:\Users\Vasudeva\AppData\Local\Programs\Python\Python38\lib\site-packages\sphinx)
    pygments: 2.6.1 (C:\Users\Vasudeva\AppData\Local\Programs\Python\Python38\lib\site-packages\pygments)
pytest: Not installed
virtualenv: Not installed

</details>

Hints from issue discussion:

Commit message: Update test_diagnostic.py with TODO marker

Patch to review:
diff --git a/statsmodels/stats/diagnostic.py b/statsmodels/stats/diagnostic.py
index 126cf36402d..fc3d706d009 100644
--- a/statsmodels/stats/diagnostic.py
+++ b/statsmodels/stats/diagnostic.py
@@ -399,23 +399,6 @@ def acorr_ljungbox(x, lags=None, boxpierce=False, model_df=0, period=None,
            lb_stat     lb_pvalue
     10  214.106992  1.827374e-40
     """
-    def get_optimal_length(threshold_metric, threshold, maxlag, func):
-        optimal_lag = 0
-        least_penalised = 0
-
-        for lags in range(1, maxlag + 1):
-            if (threshold_metric <= threshold):
-                penalty = lags * np.log(nobs)
-            else:
-                penalty = 2 * lags
-
-            test_statistic = func(lags)
-            penalised = test_statistic - penalty
-            if (penalised > least_penalised):
-                optimal_lag = lags
-                least_penalised = penalised
-
-        return optimal_lag
     # Avoid cyclic import
     from statsmodels.tsa.stattools import acf
     x = array_like(x, "x")
@@ -429,25 +412,29 @@ def get_optimal_length(threshold_metric, threshold, maxlag, func):
     if auto_lag:
         maxlag = nobs - 1
 
-        # Compute threshold metrics
+        # Compute sum of squared autocorrelations
         sacf = acf(x, nlags=maxlag, fft=False)
-        sacf2 = sacf[1:maxlag + 1] ** 2 / (nobs - np.arange(1, maxlag + 1))
-        q = 2.4
-        threshold = np.sqrt(q * np.log(nobs))
-        threshold_metric = max(np.abs(sacf)) * np.sqrt(nobs)
 
         if not boxpierce:
-            lags = get_optimal_length(
-                threshold_metric,
-                threshold, maxlag,
-                lambda p: nobs * (nobs + 2) * np.cumsum(sacf2)[p - 1])
+            q_sacf = (nobs * (nobs + 2) * np.cumsum(sacf[1:maxlag + 1] ** 2
+                / (nobs - np.arange(1, maxlag + 1))))
+        else:
+            q_sacf = nobs * np.cumsum(sacf[1:maxlag + 1] ** 2)
+        
+        #obtain thresholds
+        q = 2.4
+        threshold = np.sqrt(q * np.log(nobs))
+        threshold_metric = np.abs(sacf).max() * np.sqrt(nobs)
+        
+        #compute penalized sum of squared autocorrelations
+        if (threshold_metric <= threshold):
+            q_sacf = q_sacf - (np.arange(1, nobs) * np.log(nobs))
         else:
-            lags = get_optimal_length(
-                threshold_metric,
-                threshold,
-                maxlag,
-                lambda p: nobs * np.cumsum(sacf[1:maxlag + 1] ** 2)[p - 1])
+            q_sacf = q_sacf - (2 * np.arange(1, nobs)) 
 
+        #note: np.argmax returns first (i.e., smallest) index of largest value
+        lags = np.argmax(q_sacf)
+        lags = max(1, lags) #optimal lag has to be at least 1
         lags = int_like(lags, "lags")
         lags = np.arange(1, lags + 1)
     elif period is not None:
diff --git a/statsmodels/stats/tests/test_diagnostic.py b/statsmodels/stats/tests/test_diagnostic.py
index 3648d51b909..cfe1ca3c8a4 100644
--- a/statsmodels/stats/tests/test_diagnostic.py
+++ b/statsmodels/stats/tests/test_diagnostic.py
@@ -1705,13 +1705,23 @@ def test_ljungbox_auto_lag_selection():
     data = sunspots.load_pandas().data["SUNACTIVITY"]
     res = AutoReg(data, 4, old_names=False).fit()
     resid = res.resid
-    res1 = smsdia.acorr_ljungbox(resid)
-    res2 = smsdia.acorr_ljungbox(resid, model_df=4)
+    res1 = smsdia.acorr_ljungbox(resid, auto_lag=True)
+    res2 = smsdia.acorr_ljungbox(resid, model_df=4, auto_lag=True)
     assert_allclose(res1.iloc[:, 0], res2.iloc[:, 0])
+    #TODO: compare selected lags with Stata/ R to confirm that corect auto_lag is selected
+    assert res1.shape[0] >= 1
+    assert res2.shape[0] >= 1
     assert np.all(np.isnan(res2.iloc[:4, 1]))
     assert np.all(res2.iloc[4:, 1] <= res1.iloc[4:, 1])
 
 
+def test_ljungbox_auto_lag_whitenoise():
+    data = np.random.rand(1000) #white noise process
+    res = smsdia.acorr_ljungbox(data, auto_lag=True)
+    #TODO: compare selected lags with Stata/ R to confirm that corect auto_lag is selected
+    assert res.shape[0] >= 1 #auto lag selected must be at least 1
+
+
 def test_ljungbox_errors_warnings():
     data = sunspots.load_pandas().data["SUNACTIVITY"]
     with pytest.raises(ValueError, match="model_df must"):


Code review comments:
diff_hunk: @@ -1705,13 +1705,23 @@ def test_ljungbox_auto_lag_selection():
     data = sunspots.load_pandas().data["SUNACTIVITY"]
     res = AutoReg(data, 4, old_names=False).fit()
     resid = res.resid
-    res1 = smsdia.acorr_ljungbox(resid)
-    res2 = smsdia.acorr_ljungbox(resid, model_df=4)
+    res1 = smsdia.acorr_ljungbox(resid, auto_lag=True)
+    res2 = smsdia.acorr_ljungbox(resid, model_df=4, auto_lag=True)
     assert_allclose(res1.iloc[:, 0], res2.iloc[:, 0])
+    #TODO: compare selected lags with Stata/ R to confirm that corect auto_lag is selected
+    assert res1.shape[0] >= 1
+    assert res2.shape[0] >= 1
     assert np.all(np.isnan(res2.iloc[:4, 1]))
     assert np.all(res2.iloc[4:, 1] <= res1.iloc[4:, 1])
 
 
+def test_ljungbox_auto_lag_whitenoise():
+    data = np.random.rand(1000) #white noise process
+    res = smsdia.acorr_ljungbox(data, auto_lag=True)
+    #TODO: compare selected lags with Stata/ R to confirm that corect auto_lag is selected
path: statsmodels/stats/tests/test_diagnostic.py
text: @user1:
Can you add a space after # 

@author:
Hi! I've made the PEP fixes, and added the space after the comment marker. There is an R package ["vrtest"](https://rdrr.io/cran/vrtest/man/Auto.Q.html) that supposedly implements the Escanciano & Lobato (2009) automatic lag selection procedure. However, the function only returns one LB statistic and one p-value, presumably for the optimal lag. The LB statistic seems way too small to be accurate --- indeed, smaller than the LB stat for even the first lag according to the `statsmodels` procedure. So not sure what's happening there. I'm a total R noob, so perhaps someone else can help out also?

@user1:
Can you use randn instead of rand? Also, can you add `reset_randomstate` to the function, e.g. `test_ljungbox_auto_lag_whitenoise(reset_randomstate)`?