Repository: statsmodels/statsmodels
Pull request: ENH: Added fft to ccovf and ccf
- [x] closes #5199
- [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). 


I have changed `ccovf` to use `scipy.signal.correlate` instead of `np.correlate` which can use FFT. The way to enable FFT is chosen to be similar to how `acovf` does it and not as `scipy.signal.correlate`. 

A change is done to how the adjusted `d` is calculated  to speed up the calculation, an example is given here:  
``` python
nobs = 100
d1 = np.arange(nobs, 0, -1)
xi = np.ones(nobs)
d2 = np.correlate(xi, xi, "full")[nobs - 1 :]
np.allclose(d1, d2)
```

At last, I have updated some docstrings with minor corrections, an if statement, and added a FIXME comment, but I can remove them if you want to have it in a separate PR.


<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>


Problem statement:
Speed up the cross-correlation function
Dears,

I noticed that the computation of the cross corrlation function is extremely slower than the autocorrelation function.

It seems to me that this is because of two reasons, but I am not an expert in the field: 
* acf function lets fix the nlags parameter and ccf does not, the latter compute every lag (commonly unneceary)
* acf 'does a trick' with FFT while ccf does not.

Is it possible to speed up the ccf?

Hints from issue discussion:
It's possible to add all the same methods as in adf to ccf.

I guess, without looking, that acf and ccf could share most of the code

Commit message: Fix linting

Patch to review:
diff --git a/statsmodels/tools/validation/validation.py b/statsmodels/tools/validation/validation.py
index 27727d34323..afbb5580e86 100644
--- a/statsmodels/tools/validation/validation.py
+++ b/statsmodels/tools/validation/validation.py
@@ -59,7 +59,7 @@ def array_like(
         of obj (if present) or uses NumPy to automatically detect the dtype
     ndim : {int, None}
         Required number of dimensions of obj. If None, no check is performed.
-        If the numebr of dimensions of obj is less than ndim, additional axes
+        If the number of dimensions of obj is less than ndim, additional axes
         are inserted on the right. See examples.
     maxdim : {int, None}
         Maximum allowed dimension.  Use ``maxdim`` instead of ``ndim`` when
diff --git a/statsmodels/tsa/stattools.py b/statsmodels/tsa/stattools.py
index b93639c1fe6..e6a40644193 100644
--- a/statsmodels/tsa/stattools.py
+++ b/statsmodels/tsa/stattools.py
@@ -15,6 +15,7 @@
 import pandas as pd
 from scipy import stats
 from scipy.interpolate import interp1d
+from scipy.signal import correlate
 
 from statsmodels.regression.linear_model import OLS, yule_walker
 from statsmodels.tools.sm_exceptions import (
@@ -397,7 +398,7 @@ def acovf(x, adjusted=False, demean=True, fft=True, missing="none", nlag=None):
         If True, then denominators is n-k, otherwise n.
     demean : bool, default True
         If True, then subtract the mean x from each element of x.
-    fft : bool, default None
+    fft : bool, default True
         If True, use FFT convolution.  This method should be preferred
         for long time series.
     missing : str, default "none"
@@ -586,7 +587,7 @@ def acf(
     qstat : bool, default False
         If True, returns the Ljung-Box q statistic for each autocorrelation
         coefficient.  See q_stat for more information.
-    fft : bool, default None
+    fft : bool, default True
         If True, computes the ACF via FFT.
     alpha : scalar, default None
         If a number is given, the confidence intervals for the given level are
@@ -731,8 +732,7 @@ def pacf_yw(x, nlags=None, method="adjusted"):
     nlags = int_like(nlags, "nlags", optional=True)
     nobs = x.shape[0]
     if nlags is None:
-        if nlags is None:
-            nlags = min(int(10 * np.log10(nobs)), nobs - 1)
+        nlags = min(int(10 * np.log10(nobs)), nobs - 1)
 
     method = string_like(method, "method", options=("adjusted", "mle"))
     pacf = [1.0]
@@ -1019,7 +1019,7 @@ def pacf(x, nlags=None, method="ywadjusted", alpha=None):
 
 
 @deprecate_kwarg("unbiased", "adjusted")
-def ccovf(x, y, adjusted=True, demean=True):
+def ccovf(x, y, adjusted=True, demean=True, fft=True):
     """
     Calculate the crosscovariance between two series.
 
@@ -1028,24 +1028,23 @@ def ccovf(x, y, adjusted=True, demean=True):
     x, y : array_like
        The time series data to use in the calculation.
     adjusted : bool, optional
-       If True, then denominators for autocovariance is n-k, otherwise n.
+       If True, then denominators for crosscovariance is n-k, otherwise n.
     demean : bool, optional
         Flag indicating whether to demean x and y.
+    fft : bool, default True
+        If True, use FFT convolution.  This method should be preferred
+        for long time series.
 
     Returns
     -------
     ndarray
         The estimated crosscovariance function.
-
-    Notes
-    -----
-    This uses np.correlate which does full convolution. For very long time
-    series it is recommended to use fft convolution instead.
     """
     x = array_like(x, "x")
     y = array_like(y, "y")
     adjusted = bool_like(adjusted, "adjusted")
     demean = bool_like(demean, "demean")
+    fft = bool_like(fft, "fft", optional=False)
 
     n = len(x)
     if demean:
@@ -1055,24 +1054,28 @@ def ccovf(x, y, adjusted=True, demean=True):
         xo = x
         yo = y
     if adjusted:
-        xi = np.ones(n)
-        d = np.correlate(xi, xi, "full")
+        d = np.arange(n, 0, -1)
     else:
         d = n
-    return (np.correlate(xo, yo, "full") / d)[n - 1 :]
+
+    method = "fft" if fft else "direct"
+    return correlate(xo, yo, "full", method=method)[n - 1:] / d
 
 
 @deprecate_kwarg("unbiased", "adjusted")
-def ccf(x, y, adjusted=True):
+def ccf(x, y, adjusted=True, fft=True):
     """
     The cross-correlation function.
 
     Parameters
     ----------
     x, y : array_like
-       The time series data to use in the calculation.
+        The time series data to use in the calculation.
     adjusted : bool
-       If True, then denominators for autocovariance is n-k, otherwise n.
+        If True, then denominators for cross-correlation is n-k, otherwise n.
+    fft : bool, default True
+        If True, use FFT convolution.  This method should be preferred
+        for long time series.
 
     Returns
     -------
@@ -1081,16 +1084,14 @@ def ccf(x, y, adjusted=True):
 
     Notes
     -----
-    This is based np.correlate which does full convolution. For very long time
-    series it is recommended to use fft convolution instead.
-
     If adjusted is true, the denominator for the autocovariance is adjusted.
     """
     x = array_like(x, "x")
     y = array_like(y, "y")
     adjusted = bool_like(adjusted, "adjusted")
+    fft = bool_like(fft, "fft", optional=False)
 
-    cvf = ccovf(x, y, adjusted=adjusted, demean=True)
+    cvf = ccovf(x, y, adjusted=adjusted, demean=True, fft=fft)
     return cvf / (np.std(x) * np.std(y))
 
 
diff --git a/statsmodels/tsa/tests/test_stattools.py b/statsmodels/tsa/tests/test_stattools.py
index 0e9ecb3ddfa..541024f3058 100644
--- a/statsmodels/tsa/tests/test_stattools.py
+++ b/statsmodels/tsa/tests/test_stattools.py
@@ -30,6 +30,7 @@
 from statsmodels.tsa.stattools import (
     acf,
     acovf,
+    ccovf,
     adfuller,
     arma_order_select_ic,
     breakvar_heteroskedasticity_test,
@@ -990,6 +991,30 @@ def test_acovf_fft_vs_convolution(demean, adjusted):
     assert_almost_equal(F1, F2, decimal=7)
 
 
+@pytest.mark.parametrize("demean", [True, False])
+@pytest.mark.parametrize("adjusted", [True, False])
+def test_ccovf_fft_vs_convolution(demean, adjusted):
+    np.random.seed(1)
+    x = np.random.normal(size=128)
+    y = np.random.normal(size=128)
+
+    F1 = ccovf(x, y, demean=demean, adjusted=adjusted, fft=False)
+    F2 = ccovf(x, y, demean=demean, adjusted=adjusted, fft=True)
+    assert_almost_equal(F1, F2, decimal=7)
+
+
+@pytest.mark.parametrize("demean", [True, False])
+@pytest.mark.parametrize("adjusted", [True, False])
+@pytest.mark.parametrize("fft", [True, False])
+def test_compare_acovf_vs_ccovf(demean, adjusted, fft):
+    np.random.seed(1)
+    x = np.random.normal(size=128)
+
+    F1 = acovf(x, demean=demean, adjusted=adjusted, fft=fft)
+    F2 = ccovf(x, x, demean=demean, adjusted=adjusted, fft=fft)
+    assert_almost_equal(F1, F2, decimal=7)
+
+
 @pytest.mark.smoke
 @pytest.mark.slow
 def test_arma_order_select_ic():
@@ -999,7 +1024,7 @@ def test_arma_order_select_ic():
     arparams = np.array([0.75, -0.25])
     maparams = np.array([0.65, 0.35])
     arparams = np.r_[1, -arparams]
-    maparam = np.r_[1, maparams]
+    maparam = np.r_[1, maparams]  # FIXME: Never used
     nobs = 250
     np.random.seed(2014)
     y = arma_generate_sample(arparams, maparams, nobs)


Code review comments:
diff_hunk: @@ -990,6 +991,30 @@ def test_acovf_fft_vs_convolution(demean, adjusted):
     assert_almost_equal(F1, F2, decimal=7)
 
 
+@pytest.mark.parametrize("demean", [True, False])
+@pytest.mark.parametrize("adjusted", [True, False])
+def test_ccovf_fft_vs_convolution(demean, adjusted):
+    np.random.seed(1)
+    x = np.random.normal(size=128)
+    y = np.random.normal(size=128)
+
+    F1 = ccovf(x, y, demean=demean, adjusted=adjusted, fft=False)
+    F2 = ccovf(x, y, demean=demean, adjusted=adjusted, fft=True)
+    assert_almost_equal(F1, F2, decimal=7)
+
+
+@pytest.mark.parametrize("demean", [True, False])
+@pytest.mark.parametrize("adjusted", [True, False])
+@pytest.mark.parametrize("fft", [True, False])
+def test_compare_acovf_vs_ccovf(demean, adjusted, fft):
path: statsmodels/tsa/tests/test_stattools.py
text: @user1:
pass `reset_randomstate` as an argument to the test function. Then you do not need to set `np.random.seed` internally.

https://github.com/statsmodels/statsmodels/blob/77bb1d276c7d11bc8657497b4307aa7575c3e65c/statsmodels/stats/tests/test_statstools.py#L86

This is a global fixture.