ExponentionalSmoothing forecasts bug?

[mloning]

Describe the bug

Some forecasts from ExponentialSmoothing(y_train, trend='add', damped=True) seem to be way off, see code and graph below where "replicated" comes from the code below and "original" is the forecast reported in the M4 competition, only part of the training series is shown. Data comes from the M4 competition.

• Seems to happen only on longer series (below one example, D703 series from the M4)
• Returned forecasts are constant and equal to level attribute.
• With damped=False seems to work as expected.

Code Sample, a copy-pastable example if possible

To run code, clone M4 methods Github repo and change paths accordingly.

# load packages
import pandas as pd
import numpy as np
import os
from statsmodels.tsa.holtwinters import ExponentialSmoothing
import matplotlib.pyplot as plt

# define paths
repodir = ".../m4-methods/"  # repo root dir
datadir = os.path.join(repodir, "Dataset")
traindir = os.path.join(datadir, 'Train')
testdir = os.path.join(datadir, 'Test')

# select series from M4 dataset
dataset = 'Daily'
series_id = dataset[0] + str(703)

# load data
alltrain = pd.read_csv(os.path.join(traindir, f'{dataset}-train.csv'), index_col=0)
alltest = pd.read_csv(os.path.join(testdir, f'{dataset}-test.csv'), index_col=0)

y_train = alltrain.loc[series_id].dropna().reset_index(drop=True)
y_test = alltest.loc[series_id].dropna().reset_index(drop=True)
y_test.index = y_test.index + y_train.shape[0]

y_preds_original = pd.read_csv(os.path.join(repodir, 'Point Forecasts', 'submission-Damped.csv'), index_col=0)
y_pred_original = y_preds_original.loc[series_id].reset_index(drop=True).dropna()
y_pred_original.index = y_pred_original.index + y_train.shape[0]

# fit/forecast
m = ExponentialSmoothing(y_train, trend='add', damped=True)
mf = m.fit()
y_pred = mf.forecast(14)

# plot series/forecasts
fig, ax = plt.subplots(1)
y_train.iloc[-100:].plot(ax=ax, label='train')
y_test.plot(ax=ax, label='test')
y_pred.plot(ax=ax, label='replicated')
y_pred_original.plot(ax=ax, label='M4 forecasts');
plt.legend();

Expected Output

Forecasts closer to forecasts reported in M4 competition based on R forecast package or error/warning message.

Output of import statsmodels.api as sm; sm.show_versions()

INSTALLED VERSIONS

Python: 3.7.3.final.0
OS: Darwin 18.6.0 Darwin Kernel Version 18.6.0: Thu Apr 25 23:16:27 PDT 2019; root:xnu-4903.261.4~2/RELEASE_X86_64 x86_64
byteorder: little
LC_ALL: None
LANG: en_GB.UTF-8

Statsmodels

Installed: 0.9.0 (/Users/mloning/.conda/envs/sktime/lib/python3.7/site-packages/statsmodels)

Required Dependencies

cython: 0.29.7 (/Users/mloning/.conda/envs/sktime/lib/python3.7/site-packages/Cython)
numpy: 1.16.4 (/Users/mloning/.conda/envs/sktime/lib/python3.7/site-packages/numpy)
scipy: 1.2.1 (/Users/mloning/.conda/envs/sktime/lib/python3.7/site-packages/scipy)
pandas: 0.24.2 (/Users/mloning/.conda/envs/sktime/lib/python3.7/site-packages/pandas)
dateutil: 2.8.0 (/Users/mloning/.conda/envs/sktime/lib/python3.7/site-packages/dateutil)
patsy: 0.5.1 (/Users/mloning/.conda/envs/sktime/lib/python3.7/site-packages/patsy)

Optional Dependencies

matplotlib: 3.0.3 (/Users/mloning/.conda/envs/sktime/lib/python3.7/site-packages/matplotlib)
backend: module://ipykernel.pylab.backend_inline
cvxopt: Not installed
joblib: 0.12.5 (/Users/mloning/.conda/envs/sktime/lib/python3.7/site-packages/joblib)

Developer Tools

IPython: 7.4.0 (/Users/mloning/.conda/envs/sktime/lib/python3.7/site-packages/IPython)
jinja2: 2.7.3 (/Users/mloning/.conda/envs/sktime/lib/python3.7/site-packages/jinja2)
sphinx: 2.0.1 (/Users/mloning/.conda/envs/sktime/lib/python3.7/site-packages/sphinx)
pygments: 2.3.1 (/Users/mloning/.conda/envs/sktime/lib/python3.7/site-packages/pygments)
nose: Not installed
pytest: 4.4.2 (/Users/mloning/.conda/envs/sktime/lib/python3.7/site-packages)
virtualenv: Not installed

[bashtage]

Any chance you could make the example more self-contained? For example, using simulated data that is like the data in the competition?

[mloning]

I tried, but couldn't replicate the error with any of the simulated datasets that I tried, but it does happen on quite a few of the hourly and daily M4 datasets.

[mloning]

Damping in general does not seem to work properly for me as it returns a constant line in many cases, but that may be intended, not sure:

# generate data
rng = np.random.RandomState(1)
n = 100
y = np.zeros(n)
y[0] = 3
alpha = .75
beta = 0.01
for i in range(1, n):
    y[i] = (alpha * y[i - 1]) + (beta * i) + rng.normal(loc=1, scale=0.5)
y = pd.Series(y)

# split into train/test
fh = 20 # forecast horizon
y_train = pd.Series(y[:-fh])
y_test = pd.Series(y[-fh:])

# fit/forecast
m = ExponentialSmoothing(y_train, trend='add', damped=False)
mf = m.fit()
y_undamped = mf.forecast(fh)

m = ExponentialSmoothing(y_train, trend='add', damped=True)
mf = m.fit()
y_damped = mf.forecast(fh)

# plot
fig, ax = plt.subplots(1)
y_train.iloc[-100:].plot(ax=ax, label='train')
y_test.plot(ax=ax, label='test')
y_damped.plot(ax=ax, label='damped')
y_undamped.plot(ax=ax, label='undamped')
plt.legend()

• I'm not sure why, maybe the damping slope parameter should be bounded to [0.8, 0.98] as described in Forecasting: Principles and Practice?

• Was Holt Winters tested against R forecast package? As said, there seem to be quite a few datasets in the M4 competition where forecasts from statsmodels deviate from the reported results?

[bashtage]

Fixed in #5893

[bashtage]

Appears that this was not fixed by #5893

[timdhondt1]

I can confirm this bug. Fitted values are fine but forecasts seem to be way off for any model that uses additive trend, although I haven't found the bug in the actual code yet. My guess is that something goes wrong using the trended(lvls, b) function when calculating the fitted values

[ChadFulton]

An update on this issue, although I haven't had time to fully look into it.

I'm not sure that this is a "bug" in the sense of incorrect results from wrong code. Instead, what happens is that the trend smoothing parameter is getting fitted to be exactly zero. When there is no damping, this means that the trend is always fixed at the initial trend, which in e.g. the latest example here looks about right for the series as a whole (but in general wouldn't necessarily work out so well). When there is damping (even a relatively small amount), the initial trend gets damped to zero pretty quickly, and you get the flat line.

But it seems this shouldn't happen (and doesn't happen in R), so more investigation is definitely necessary.

[bashtage]

According to @mloning sounds like you have to restrict the parameter to be high (>0.8) when included.

[ChadFulton]

According to @mloning sounds like you have to restrict the parameter to be high (>0.8) when included.

This wasn't enough to make it work for me, since even relatively large damping parameters will damp the initial trend to zero relatively quickly. But I know that R has additional parameter restrictions that we don't use, e.g. on the trend smoothing parameter, and maybe in general this is a problem of parameter restrictions. I haven't had a chance to look into it yet though.

[ChadFulton]

I finally looked at the results from ets in R, and it is giving the same type of forecast when the trend is damped here - i.e. it estimates the smoothing trend to be zero and then the damping implies that the trend at the end of the sample about 0.

We do have different estimates because we don't restrict the damping parameter, which we estimate as about 0.2, while ets gives 0.98. Given the length of the sample, this doesn't have a practical effect on the forecasts.

So I'm not sure there's a reproducible bug here. Clearly the graph in the first issue comment is concerning, but I think that we haven't been able to reproduce the problem?

[mloning]

@ChadFulton Do you get a different forecast when you run the code on the selected series from the M4 datasets (using my code above)?

[ChadFulton]

@ChadFulton Do you get a different forecast when you run the code on the selected series from the M4 datasets (using my code above)?

Thanks for pointing this out. Yes, you're right, I get the same forecast and that's clearly a bug.

[ChadFulton]

I'm not quite sure what this bug is. It seems like a combination of the following:

• Really bad starting parameters (smoothing_level = 0 and smoothing_slope = 0, initial_slope < 0 when the series is actually trending upwards).
• Starting parameters on the bounds of the constraints.
• Long time series means that these starting parameters yield a huge sum of squared errors.

[ChadFulton]

Although something that I don't understand is going on with l-bfgs-b:

reload(holtwinters)
mod = holtwinters.ExponentialSmoothing(endog, trend='add', damped=True)
res = mod.fit()
print(res.sse)
print('--------------')
print(res.mle_retvals)

yields:

121510574038.16998
--------------
      fun: 30893297.84558988
 hess_inv: <5x5 LbfgsInvHessProduct with dtype=float64>
      jac: array([1.21478039e+19, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00])
  message: b'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'
     nfev: 12
      nit: 1
   status: 0
  success: True
        x: array([0.00000000e+00, 0.00000000e+00, 4.12770000e+03, 0.00000000e+00,
       7.36842105e-01])

but the x that is given in those results does not yield an SSE that corresponds with the fun result from minimize; i.e. those x values lead to the res.sse value of 121510574038.16998 and do not produce the SSE from fun of 30893297.84558988.

[ChadFulton]

Solutions for this particular problem include:

• Forcing starting parameters away from the bounds
• Making the bounds (1e-5, 1 - 1e-5) instead of (0, 1)
• Using either BFGS or SLSQP instead of L-BFGS-B
• Using naive starting parameters for the smoothing level and slope (i.e. both set to 0.5) rather than the values found by brute (which are 1.0 and 0.0, respectively).

Hard to know how this might generalize to other problems.

[ChadFulton]

I think that what might be best to implement, just to get a fix in, will be:

Making the bounds (1e-5, 1 - 1e-5) instead of (0, 1)

This is what the R forecast package does, at least for some of the parameters.

If problems still persist, we can re-evaluate some of the other solutions.

[mloning]

I think it comes from the seasonality adjustment in ExponentialSmoothing. I've managed to reproduce the published results when using my own de-seasonalization wrapper based on statsmodels' seasonal_decompose function.

[bashtage]

fixed in #6870