Exponential smoothing with seasonality initial Level value is weird

[aelavender]

For Exponential Smoothing with seasonality, the initial Level (if not provided by the user) is set as follows:
y[np.arange(self.nobs) % m == 0].mean()

statsmodels/statsmodels/tsa/holtwinters.py

Line 769 in da9d7e9

l0 = y[np.arange(self.nobs) % m == 0].mean() if l0 is None else l0

In other words, take the mean of all the points corresponding to the first season.

This is an unusual way to initialize, and produces strange results for time series with a strongly non-zero slope: you'll start with a Level close to the mean of the series (which might not be anywhere close to the level at the start of the series). Then you'll have huge initial seasonality to compensate. The result is a forecast that looks fine in total but whose components are weird. Ideally, the mean of the seasonality time series should be zero (for addititive seasonality) or one (for multiplicative), but this method results in values far from that.

Hyndman (whose text is mentioned in the documentation) provides some suggestions. His simplest method is to take the mean of the first period (not the first season of each period, as is currently implemented above).

NIST's method is even simpler: just use y[0], identical to the non-seasonal version. (They don't explicitly mention it, but the example on the next page shows that result.)

I think we should use Hyndman's 1998 method, and change L769 to be

l0 = y[:m].mean() if l0 is None else l0

I can put up a PR if you wish.

Thanks!

[ChadFulton]

Thanks for this! I completely agree, and a PR would be very welcome.

[aelavender]

Here's an example. The source time series is generated from a random walk with drift, with weekly seasonality overlaid. I'm using additive trend and seasonality.

Note how, while the fitted data looks good, the trend is much higher than the data, and the seasonality is very negative to compensate.

[aelavender]

Thanks for the quick response, @ChadFulton! I'll put up a PR within a day or so.

[bashtage]

This seems like a change that is preferring aesthetics of the components over statistically sound choices. As @aelavender pointed out, the other components can compensate (fully) for this. It may also be necessary to change the initialization of other components when changing this setting or it will change the fit of the models in a bad way.

[ChadFulton]

I'm not sure I follow your reasoning - the seasonal term is not being normalized here, which I think is not just a case of aesthetics. For one thing, if the initial level and seasonals are estimated (as they are by default) then the model parameters are not identified.

Second, it means that a model including just a level term will differ from a model that also includes "seasonal" effects, but where the seasonal values are all equal. These models should produce identical results.

[aelavender]

@bashtage The initial seasonal components are rescaled (to be normalised) as a part of this change, since they are defined as s0 = list(y[:m] / l0) if seasonal == 'mul' else list(y[:m] - l0). b0 doesn't change, since it's calculated directly from the y's, but its evolution (which depends on l) does.

[josef-pkt]

l0 = y[:m].mean() if l0 is None else l0

Is it possible to change the name l0?
In the mail that I receive of issues, I was surprised why there is a hardcoded ten 10. Same for the font used in the github issue comments.

[bashtage]

Second, it means that a model including just a level term will differ from a model that also includes "seasonal" effects, but where the seasonal values are all equal. These models should produce identical results.

I don't see this -- it looks like the m-equal seasonals will pick up the average level and the level with look like a mean data.mean()-ish random walk. Something like the picture above but where the green curve is flat around -5000. Setting the seasonals to 0 will force the level to track the data.

The big question is are these two parameterizations identical. IOW do the same smoothing parameters but different initializations produce the same fitted values? If so, then it seems nearly free, even though it would break backward compatibility. The alternative would be to set it as an option (beyond allowing the initialization to directly user set).

FWIW in some testing this function tended to produce better fits than Hyndaman's choices in R, often with substantial differences in smoothing parameters.

[aelavender]

@josef-pkt lol, I agree. I'd prefer l0 to be written as either L0 or level0, and l[t] as either L[t] or level[t]. But changing it throughout would be a massive diff, so probably worth opening a separate issue for it.

[aelavender]

Something like the picture above but where the green curve is flat around -5000.

Agree, this is what would happen -- seasonality would start at the -1 * mean(y) (or 1 / mean(y) for multiplicative) and roughly stay there. The forecast would not change, only the components would.

The big question is are these two parameterizations identical. IOW do the same smoothing parameters but different initializations produce the same fitted values?

There are 8 variations (additive or multiplicative for trend and seasonality; damped or not), and I suspect that not all of them produce identical forecasts with this change. Agree that this would break backwards compatibility, since in current release, initial values cannot be set directly by the user.

FWIW in some testing this function tended to produce better fits than Hyndaman's choices in R, often with substantial differences in smoothing parameters.

This is very interesting; might you have research that you could share? Most sources that I've seen suggest that normalising is a good idea.

[bashtage]

I didn't keep anything -- just some quick experiments against some fo the canonical datasets like airlines. I was trying to verify sm against R until I realized that the versions were different seemingly by design. I suspect that averaging reduce the influence of some noise in the data.

…

On Fri, Jan 4, 2019 at 7:06 PM Lavender ***@***.***> wrote: Something like the picture above but where the green curve is flat around -5000. Agree, this is what would happen -- seasonality would start at the -1 * mean(y) (or 1 / mean(y) for multiplicative) and roughly stay there. The forecast would not change, only the components would. The big question is are these two parameterizations identical. IOW do the same smoothing parameters but different initializations produce the same fitted values? There are 8 variations (additive or multiplicative for trend and seasonality; damped or not), and I suspect that not all of them produce identical forecasts with this change. Agree that this would break backwards compatibility, since in current release, initial values cannot be set directly by the user. FWIW in some testing this function tended to produce better fits than Hyndaman's choices in R, often with substantial differences in smoothing parameters. This is very interesting; might you have research that you could share? Most sources that I've seen suggest that normalising is a good idea. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#5438 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AFU5RT8GnM6dCQ9Y7Ql59hpT9OgIv-5Mks5u_6ZFgaJpZM4ZjUyk> .