Repository: alkaline-ml/pmdarima (language: python)

Issue:
SARIMA prediction returns differenced time series when simple_differencing = True
Hi, thank you for the useful library in Python.
I'm working on a long time series with n ~ 2000 with period ~ 300.
I use auto_arima() for SARIMA with simple_differencing = True. Otherwise computation takes forever.

**Describe the bug**
With simple_differencing, predict() function returns a differenced time series.
Then I have to use diff_inv() for point forecast outside pmdarima.
Furthermore, it is not very straightforward to undifference confidence interval.
I'm wondering if pmdarima could return the undifferenced forecast and confidence interval.

**To Reproduce**  

```python
import pmdarima as pm
wineind = pm.datasets.load_wineind()

def do_fit(simple_differencing):
    return pm.auto_arima(wineind, start_p=1, start_q=1, max_p=1,
                         max_q=2, m=12, start_P=0,
                         seasonal=True,
                         d=1, D=1, stepwise=True,
                         error_action='ignore',
                         sarimax_kwargs={
                             'simple_differencing': simple_differencing
                         },
                         maxiter=2,
                         trace=True)

# simple_differencing = False
seasonal_fit = do_fit(False)
prediction = seasonal_fit.predict(n_periods=60, return_conf_int=True)
pred_mid = prediction[0]
pred_conf_int = prediction[1]

# simple_differencing = True
simple = do_fit(True)
prediction_simple = simple.predict(n_periods=60, return_conf_int=True)
pred_simple_mid = prediction_simple[0]
pred_simple_conf_int = prediction_simple[1]
```
The plot is
![image](https://user-images.githubusercontent.com/16064801/113437443-8ae49580-9421-11eb-874a-6ed357a668fc.png)

Differenced time series is returned because it inherits SARIMAX in statsmodels:
[https://www.statsmodels.org/stable/generated/statsmodels.tsa.statespace.sarimax.SARIMAX.html](https://www.statsmodels.org/stable/generated/statsmodels.tsa.statespace.sarimax.SARIMAX.html)

> If simple_differencing = True is used, then the endog and exog data are differenced prior to putting the model in state-space form. This has the same effect as if the user differenced the data prior to constructing the model, which has implications for using the results:

> - _Forecasts and predictions will be about the differenced data, not about the original data. (while if simple_differencing = False is used, then forecasts and predictions will be about the original data)._

**Versions**
<!--
Please run the following snippet and paste the output below.
import pmdarima; pmdarima.show_versions()

For pmdarima versions <1.5.2 use this:
import platform; print(platform.platform())
import sys; print("Python", sys.version)
import pmdarima; print("pmdarima", pmdarima.__version__)
import numpy; print("NumPy", numpy.__version__)
import scipy; print("SciPy", scipy.__version__)
import sklearn; print("Scikit-Learn", sklearn.__version__)
import statsmodels; print("Statsmodels", statsmodels.__version__)
-->
python: 3.7.6
machine: Windows-10-10.0.18362-SP0
statsmodels: 0.11.0
pmdarima: 1.8.0

**Expected behavior**
predict() function returns undifferenced time series. This is more useful than a differenced series. Also, this behavior is more consistent with auto.arima() in R.  The R function returns undifferenced series for both approximation = true/false.

**Solution**
My change is [https://github.com/drxyzw/pmdarima/commit/4853b6267892f6595250200dddf09893c3a42ee5](https://github.com/drxyzw/pmdarima/commit/4853b6267892f6595250200dddf09893c3a42ee5)
- Adding a routine inside _seasonal_prediction_with_confidence() in arima.py.
- For forecast mid, use diff_inv() to undifference the forecast.
- Confidence interval is computed with coefficients of SARIMA and sigma2. I derive an equivalent MA model from SARIMA. With the coefficients of the equivalent MA and sigma2, we can compute variance and confidence interval.

The result is
![image](https://user-images.githubusercontent.com/16064801/113646693-47f61c80-96c4-11eb-9b4f-fea932d89844.png)
Indeed, results of simple_differencing = true/false are different because the selected SARIMA orders are (0,1,2)x(0,1,1)[12] / (1,1,2)x(0,1,1)[12] respectively.

If it is fine for you, I will make a pull request from this change after cleaning tests.


Discussion / hints:

Fix patch:
diff --git a/pmdarima/arima/arima.py b/pmdarima/arima/arima.py
index e7649e21b..26eedb360 100644
--- a/pmdarima/arima/arima.py
+++ b/pmdarima/arima/arima.py
@@ -12,6 +12,7 @@
 
 from scipy.stats import gaussian_kde, norm
 import numpy as np
+import numpy.polynomial.polynomial as np_polynomial
 import warnings
 
 from . import _validation as val
@@ -23,15 +24,92 @@
 from ..compat import pmdarima as pm_compat
 from ..utils import if_has_delegate, is_iterable, check_endog, check_exog
 from ..utils.visualization import _get_plt
+from ..utils.array import diff_inv, diff
 
 # Get the version
 import pmdarima
 
 __all__ = [
-    'ARIMA'
+    'ARIMA',
+    'ARMAtoMA'
 ]
 
 
+def ARMAtoMA(ar, ma, max_deg):
+    r"""
+    Convert ARMA coefficients to infinite MA coefficients.
+
+    Compute coefficients of MA model equivalent to given ARMA model.
+    MA coefficients are cut off at max_deg.
+    The same function as ARMAtoMA() in stats library of R
+
+    Parameters
+    ----------
+    ar : array-like, shape=(n_orders,)
+        The array of AR coefficients.
+
+    ma : array-like, shape=(n_orders,)
+        The array of MA coefficients.
+
+    max_deg : int
+        Coefficients are computed up to the order of max_deg.
+
+    Returns
+    -------
+    np.ndarray, shape=(max_deg,)
+        Equivalent MA coefficients.
+
+    Notes
+    -----
+    Here is the derivation. Suppose ARMA model is defined as
+    .. math::
+    x_t - ar_1*x_{t-1} - ar_2*x_{t-2} - ... - ar_p*x_{t-p}\\
+        = e_t + ma_1*e_{t-1} + ma_2*e_{t-2} + ... + ma_q*e_{t-q}
+    namely
+    .. math::
+    (1 - \sum_{i=1}^p[ar_i*B^i]) x_t = (1 + \sum_{i=1}^q[ma_i*B^i]) e_t
+    where :math:`B` is a backward operator.
+
+    Equivalent MA model is
+    .. math::
+        x_t = (1 - \sum_{i=1}^p[ar_i*B^i])^{-1}\\
+        * (1 + \sum_{i=1}^q[ma_i*B^i]) e_t\\
+        = (1 + \sum_{i=1}[ema_i*B^i]) e_t
+    where :math:``ema_i`` is a coefficient of equivalent MA model.
+    The :math:``ema_i`` satisfies
+    .. math::
+        (1 - \sum_{i=1}^p[ar_i*B^i]) * (1 + \sum_{i=1}[ema_i*B^i]) \\
+        = 1 + \sum_{i=1}^q[ma_i*B^i]
+    thus
+    .. math::
+        \sum_{i=1}[ema_i*B^i] = \sum_{i=1}^p[ar_i*B^i] \\
+        + \sum_{i=1}^p[ar_i*B^i] * \sum_{j=1}[ema_j*B^j] \\
+        + \Sum_{i=1}^q[ma_i*B^i]
+    therefore
+    .. math::
+        ema_i = ar_i (but 0 if i>p) \\
+        + \Sum_{j=1}^{min(i-1,p)}[ar_j*ema_{i-j}] + ma_i(but 0 if i>q) \\
+        = \sum_{j=1}{min(i,p)}[ar_j*ema_{i-j}(but 1 if j=i)] \\
+        + ma_i(but 0 if i>q)
+
+    Examples
+    --------
+    >>> ar = np.array([0.1])
+    >>> ma = np.empty(0)
+    >>> ARMAtoMA(ar, ma, 3)
+    array[0.1, 0.01, 0.001]
+    """
+    p = len(ar)
+    q = len(ma)
+    ema = np.empty(max_deg)
+    for i in range(0, max_deg):
+        temp = ma[i] if i < q else 0.0
+        for j in range(0, min(i + 1, p)):
+            temp += ar[j] * (ema[i - j - 1] if i - j - 1 >= 0 else 1.0)
+        ema[i] = temp
+    return ema
+
+
 def _aicc(model_results, nobs, add_constant):
     """Compute the corrected Akaike Information Criterion"""
     aic = model_results.aic
@@ -84,6 +162,42 @@ def _seasonal_prediction_with_confidence(arima_res,
 
     f = results.predicted_mean
     conf_int = results.conf_int(alpha=alpha)
+    if arima_res.specification['simple_differencing']:
+        # If simple_differencing == True, statsmodels.get_prediction returns
+        # mid and confidence intervals on differenced time series.
+        # We have to invert differencing the mid and confidence intervals
+        y_org = arima_res.model.orig_endog
+        d = arima_res.model.orig_k_diff
+        D = arima_res.model.orig_k_seasonal_diff
+        period = arima_res.model.seasonal_periods
+        # Forecast mid: undifferencing non-seasonal part
+        if d > 0:
+            y_sdiff = y_org if D == 0 else diff(y_org, period, D)
+            f_temp = np.append(y_sdiff[-d:], f)
+            f_temp = diff_inv(f_temp, 1, d)
+            f = f_temp[(2 * d):]
+        # Forecast mid: undifferencing seasonal part
+        if D > 0 and period > 1:
+            f_temp = np.append(y_org[-(D * period):], f)
+            f_temp = diff_inv(f_temp, period, D)
+            f = f_temp[(2 * D * period):]
+        # confidence interval
+        ar_poly = arima_res.polynomial_reduced_ar
+        poly_diff = np_polynomial.polypow(np.array([1., -1.]), d)
+        sdiff = np.zeros(period + 1)
+        sdiff[0] = 1.
+        sdiff[-1] = 1.
+        poly_sdiff = np_polynomial.polypow(sdiff, D)
+        ar = -np.polymul(ar_poly, np.polymul(poly_diff, poly_sdiff))[1:]
+        ma = arima_res.polynomial_reduced_ma[1:]
+        n_predMinus1 = end - start
+        ema = ARMAtoMA(ar, ma, n_predMinus1)
+        sigma2 = arima_res._params_variance[0]
+        var = np.cumsum(np.append(1., ema * ema)) * sigma2
+        q = results.dist.ppf(1. - alpha / 2, *results.dist_args)
+        conf_int[:, 0] = f - q * np.sqrt(var)
+        conf_int[:, 1] = f + q * np.sqrt(var)
+
     return check_endog(f, dtype=None, copy=False), \
         check_array(conf_int, copy=False, dtype=None)
 


Test patch:
diff --git a/pmdarima/arima/tests/test_arima.py b/pmdarima/arima/tests/test_arima.py
index 7d5ac5f76..0a45f5948 100644
--- a/pmdarima/arima/tests/test_arima.py
+++ b/pmdarima/arima/tests/test_arima.py
@@ -7,7 +7,7 @@
 import numpy as np
 import pandas as pd
 
-from pmdarima.arima import ARIMA, auto_arima, AutoARIMA
+from pmdarima.arima import ARIMA, auto_arima, AutoARIMA, ARMAtoMA
 from pmdarima.arima import _validation as val
 from pmdarima.compat.pytest import pytest_error_str
 from pmdarima.datasets import load_lynx, load_wineind, load_heartrate
@@ -712,3 +712,12 @@ def test_update_1_iter(model):
 
     # They should be close
     assert np.allclose(params1, params2, atol=0.05)
+
+
+def test_ARMAtoMA():
+    ar = np.array([0.5, 0.6])
+    ma = np.array([0.4, 0.3, 0.1, 0.05])
+    max_deg = 6
+    equivalent_ma = ARMAtoMA(ar, ma, max_deg)
+    ema_expected = np.array([0.9000, 1.3500, 1.3150, 1.5175, 1.5477, 1.6843])
+    assert_array_almost_equal(equivalent_ma, ema_expected, decimal=4)
diff --git a/pmdarima/arima/tests/test_auto.py b/pmdarima/arima/tests/test_auto.py
index 301ef37d7..eba3d5e71 100644
--- a/pmdarima/arima/tests/test_auto.py
+++ b/pmdarima/arima/tests/test_auto.py
@@ -13,6 +13,7 @@
 from pmdarima.warnings import ModelFitWarning
 from pmdarima.compat.pytest import pytest_error_str, pytest_warning_messages
 
+from numpy.testing import assert_allclose
 from numpy.testing import assert_array_almost_equal
 
 import os
@@ -364,6 +365,70 @@ def do_fit(simple_differencing):
     assert simple.seasonal_order == seasonal_fit.seasonal_order
 
 
+def test_stepwise_with_simple_differencing2():
+    def do_fit(simple_differencing):
+        return pm.auto_arima(austres, start_p=1, start_q=1, max_p=1,
+                             max_q=2, seasonal=False, d=1, stepwise=True,
+                             error_action='ignore',
+                             sarimax_kwargs={
+                                 'simple_differencing': simple_differencing
+                             },
+                             maxiter=2,
+                             trace=True)
+
+    # Without simple_differencing
+    fit = do_fit(False)
+    pred = fit.predict(n_periods=10, return_conf_int=True)
+    pred_mid = pred[0]
+    pred_ci = pred[1]
+
+    # With simple_differencing
+    fit_sd = do_fit(True)
+    pred_sd = fit_sd.predict(n_periods=10, return_conf_int=True)
+    pred_sd_mid = pred_sd[0]
+    pred_sd_ci = pred_sd[1]
+
+    # Expecting similar predictions with or without simple_differencing
+    assert_allclose(pred_mid, pred_sd_mid, rtol=0.01)
+    assert_allclose(pred_ci[:, 0], pred_sd_ci[:, 0], rtol=0.01)
+    assert_allclose(pred_ci[:, 1], pred_sd_ci[:, 1], rtol=0.01)
+
+
+# SARIMA with/without simple_differencing
+def test_stepwise_with_simple_differencing3():
+    def do_fit(simple_differencing):
+        return pm.auto_arima(wineind, start_p=1, start_q=1, max_p=1,
+                             max_q=2, m=12, start_P=0,
+                             seasonal=True,
+                             d=1, D=1, stepwise=True,
+                             error_action='ignore',
+                             sarimax_kwargs={
+                                 'simple_differencing': simple_differencing
+                             },
+                             maxiter=2,
+                             trace=True)
+
+    # Without simple_differencing
+    fit = do_fit(False)
+    pred = fit.predict(n_periods=24, return_conf_int=True)
+    pred_mid = pred[0]
+    pred_ci = pred[1]
+
+    # With simple_differencing
+    fit_sd = do_fit(True)
+    pred_sd = fit_sd.predict(n_periods=24, return_conf_int=True)
+    pred_sd_mid = pred_sd[0]
+    pred_sd_ci = pred_sd[1]
+
+    # Expecting similar predictions with or without simple_differencing
+    ave = np.average(pred_mid)
+    assert_allclose(pred_mid, pred_sd_mid, atol=ave * 0.15)
+    ave0 = np.average(pred_ci[:, 0])
+    ave1 = np.average(pred_ci[:, 1])
+    assert_allclose(pred_ci[:, 0], pred_sd_ci[:, 0], atol=0.35 * ave0)
+    assert_allclose(pred_ci[:, 1], pred_sd_ci[:, 1], atol=0.15 * ave1)
+
+
 def test_with_seasonality2():
     # show we can estimate D even when it's not there...
     pm.auto_arima(wineind, start_p=1, start_q=1, max_p=2, max_q=2, m=wineind_m,