Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
Implementation of Holt-Winters algorithms in Python 2
#The MIT License (MIT)
#
#Copyright (c) 2015 Andre Queiroz
#
#Permission is hereby granted, free of charge, to any person obtaining a copy
#of this software and associated documentation files (the "Software"), to deal
#in the Software without restriction, including without limitation the rights
#to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
#copies of the Software, and to permit persons to whom the Software is
#furnished to do so, subject to the following conditions:
#
#The above copyright notice and this permission notice shall be included in
#all copies or substantial portions of the Software.
#
#THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
#IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
#FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
#AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
#LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
#OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
#THE SOFTWARE.
#
# Holt-Winters algorithms to forecasting
# Coded in Python 2 by: Andre Queiroz
# Description: This module contains three exponential smoothing algorithms. They are Holt's linear trend method and Holt-Winters seasonal methods (additive and multiplicative).
# References:
# Hyndman, R. J.; Athanasopoulos, G. (2013) Forecasting: principles and practice. http://otexts.com/fpp/. Accessed on 07/03/2013.
# Byrd, R. H.; Lu, P.; Nocedal, J. A Limited Memory Algorithm for Bound Constrained Optimization, (1995), SIAM Journal on Scientific and Statistical Computing, 16, 5, pp. 1190-1208.
from __future__ import division
from sys import exit
from math import sqrt
from numpy import array
from scipy.optimize import fmin_l_bfgs_b
def RMSE(params, *args):
Y = args[0]
type = args[1]
rmse = 0
if type == 'linear':
alpha, beta = params
a = [Y[0]]
b = [Y[1] - Y[0]]
y = [a[0] + b[0]]
for i in range(len(Y)):
a.append(alpha * Y[i] + (1 - alpha) * (a[i] + b[i]))
b.append(beta * (a[i + 1] - a[i]) + (1 - beta) * b[i])
y.append(a[i + 1] + b[i + 1])
else:
alpha, beta, gamma = params
m = args[2]
a = [sum(Y[0:m]) / float(m)]
b = [(sum(Y[m:2 * m]) - sum(Y[0:m])) / m ** 2]
if type == 'additive':
s = [Y[i] - a[0] for i in range(m)]
y = [a[0] + b[0] + s[0]]
for i in range(len(Y)):
a.append(alpha * (Y[i] - s[i]) + (1 - alpha) * (a[i] + b[i]))
b.append(beta * (a[i + 1] - a[i]) + (1 - beta) * b[i])
s.append(gamma * (Y[i] - a[i] - b[i]) + (1 - gamma) * s[i])
y.append(a[i + 1] + b[i + 1] + s[i + 1])
elif type == 'multiplicative':
s = [Y[i] / a[0] for i in range(m)]
y = [(a[0] + b[0]) * s[0]]
for i in range(len(Y)):
a.append(alpha * (Y[i] / s[i]) + (1 - alpha) * (a[i] + b[i]))
b.append(beta * (a[i + 1] - a[i]) + (1 - beta) * b[i])
s.append(gamma * (Y[i] / (a[i] + b[i])) + (1 - gamma) * s[i])
y.append((a[i + 1] + b[i + 1]) * s[i + 1])
else:
exit('Type must be either linear, additive or multiplicative')
rmse = sqrt(sum([(m - n) ** 2 for m, n in zip(Y, y[:-1])]) / len(Y))
return rmse
def linear(x, fc, alpha = None, beta = None):
Y = x[:]
if (alpha == None or beta == None):
initial_values = array([0.3, 0.1])
boundaries = [(0, 1), (0, 1)]
type = 'linear'
parameters = fmin_l_bfgs_b(RMSE, x0 = initial_values, args = (Y, type), bounds = boundaries, approx_grad = True)
alpha, beta = parameters[0]
a = [Y[0]]
b = [Y[1] - Y[0]]
y = [a[0] + b[0]]
rmse = 0
for i in range(len(Y) + fc):
if i == len(Y):
Y.append(a[-1] + b[-1])
a.append(alpha * Y[i] + (1 - alpha) * (a[i] + b[i]))
b.append(beta * (a[i + 1] - a[i]) + (1 - beta) * b[i])
y.append(a[i + 1] + b[i + 1])
rmse = sqrt(sum([(m - n) ** 2 for m, n in zip(Y[:-fc], y[:-fc - 1])]) / len(Y[:-fc]))
return Y[-fc:], alpha, beta, rmse
def additive(x, m, fc, alpha = None, beta = None, gamma = None):
Y = x[:]
if (alpha == None or beta == None or gamma == None):
initial_values = array([0.3, 0.1, 0.1])
boundaries = [(0, 1), (0, 1), (0, 1)]
type = 'additive'
parameters = fmin_l_bfgs_b(RMSE, x0 = initial_values, args = (Y, type, m), bounds = boundaries, approx_grad = True)
alpha, beta, gamma = parameters[0]
a = [sum(Y[0:m]) / float(m)]
b = [(sum(Y[m:2 * m]) - sum(Y[0:m])) / m ** 2]
s = [Y[i] - a[0] for i in range(m)]
y = [a[0] + b[0] + s[0]]
rmse = 0
for i in range(len(Y) + fc):
if i == len(Y):
Y.append(a[-1] + b[-1] + s[-m])
a.append(alpha * (Y[i] - s[i]) + (1 - alpha) * (a[i] + b[i]))
b.append(beta * (a[i + 1] - a[i]) + (1 - beta) * b[i])
s.append(gamma * (Y[i] - a[i] - b[i]) + (1 - gamma) * s[i])
y.append(a[i + 1] + b[i + 1] + s[i + 1])
rmse = sqrt(sum([(m - n) ** 2 for m, n in zip(Y[:-fc], y[:-fc - 1])]) / len(Y[:-fc]))
return Y[-fc:], alpha, beta, gamma, rmse
def multiplicative(x, m, fc, alpha = None, beta = None, gamma = None):
Y = x[:]
if (alpha == None or beta == None or gamma == None):
initial_values = array([0.0, 1.0, 0.0])
boundaries = [(0, 1), (0, 1), (0, 1)]
type = 'multiplicative'
parameters = fmin_l_bfgs_b(RMSE, x0 = initial_values, args = (Y, type, m), bounds = boundaries, approx_grad = True)
alpha, beta, gamma = parameters[0]
a = [sum(Y[0:m]) / float(m)]
b = [(sum(Y[m:2 * m]) - sum(Y[0:m])) / m ** 2]
s = [Y[i] / a[0] for i in range(m)]
y = [(a[0] + b[0]) * s[0]]
rmse = 0
for i in range(len(Y) + fc):
if i == len(Y):
Y.append((a[-1] + b[-1]) * s[-m])
a.append(alpha * (Y[i] / s[i]) + (1 - alpha) * (a[i] + b[i]))
b.append(beta * (a[i + 1] - a[i]) + (1 - beta) * b[i])
s.append(gamma * (Y[i] / (a[i] + b[i])) + (1 - gamma) * s[i])
y.append((a[i + 1] + b[i + 1]) * s[i + 1])
rmse = sqrt(sum([(m - n) ** 2 for m, n in zip(Y[:-fc], y[:-fc - 1])]) / len(Y[:-fc]))
return Y[-fc:], alpha, beta, gamma, rmse
@surfer26th

This comment has been minimized.

Show comment Hide comment
@surfer26th

surfer26th Mar 13, 2014

This looks great -- Holt Winters with Alpha, Beta, Gamma optimization taking advantage of Scipy optimzation. This seems to be the most complete Python-based Holt Winters I could find.

Have you or anyone had experience using this?

This looks great -- Holt Winters with Alpha, Beta, Gamma optimization taking advantage of Scipy optimzation. This seems to be the most complete Python-based Holt Winters I could find.

Have you or anyone had experience using this?

@gargvikram07

This comment has been minimized.

Show comment Hide comment
@gargvikram07

gargvikram07 Jun 23, 2014

Nice implementation. Exactly what I was looking for.
Thanks mate.

Nice implementation. Exactly what I was looking for.
Thanks mate.

@gargvikram07

This comment has been minimized.

Show comment Hide comment
@gargvikram07

gargvikram07 Jun 23, 2014

Dear andrequeiroz,
Actually need more help.
Is there any java implementation of the same. ?

Thanks

Dear andrequeiroz,
Actually need more help.
Is there any java implementation of the same. ?

Thanks

@jtprince

This comment has been minimized.

Show comment Hide comment
@jtprince

jtprince Nov 22, 2014

A. what exactly do "x", "m", "fc" represent?
B. Can this handle unevenly spaced time points? Thanks!

After digging into the code a bit
A. "x" is a list of data values (not a numpy array), "m" is the periodicity of the data (i.e., the number of points per period--or something of this nature), "fc" is the number of points to forecast into the future.
B. The above implementation only deals in evenly spaced data. Tomas Hanzak has done some work (for his dissertation) on implementing Holt Winters for unevenly spaced data. He has implemented it in the DMITS software (contact him directly about getting a copy).

A. what exactly do "x", "m", "fc" represent?
B. Can this handle unevenly spaced time points? Thanks!

After digging into the code a bit
A. "x" is a list of data values (not a numpy array), "m" is the periodicity of the data (i.e., the number of points per period--or something of this nature), "fc" is the number of points to forecast into the future.
B. The above implementation only deals in evenly spaced data. Tomas Hanzak has done some work (for his dissertation) on implementing Holt Winters for unevenly spaced data. He has implemented it in the DMITS software (contact him directly about getting a copy).

@jtprince

This comment has been minimized.

Show comment Hide comment
@jtprince

jtprince Nov 25, 2014

Note: there is a minor bug in the trend initialization code (which appears in multiple places) when running in python 2 -- it's doing integer division and should be doing float division. For example:

b = [(sum(Y[m:2 * m]) - sum(Y[0:m])) / m ** 2]
# easy fix, should be something like:
b = [(sum(Y[m:2 * m]) - sum(Y[0:m])) / float(m ** 2)]

Note: there is a minor bug in the trend initialization code (which appears in multiple places) when running in python 2 -- it's doing integer division and should be doing float division. For example:

b = [(sum(Y[m:2 * m]) - sum(Y[0:m])) / m ** 2]
# easy fix, should be something like:
b = [(sum(Y[m:2 * m]) - sum(Y[0:m])) / float(m ** 2)]
@danuker

This comment has been minimized.

Show comment Hide comment
@danuker

danuker Jan 30, 2015

Note: in the multiplicative RMSE, the formula for y is copy-paste-mistaken (line 61).
It should be

 y.append((a[i + 1] + b[i + 1])* s[i + 1])

to match the one in the actual model (line 151).
I forked it.

danuker commented Jan 30, 2015

Note: in the multiplicative RMSE, the formula for y is copy-paste-mistaken (line 61).
It should be

 y.append((a[i + 1] + b[i + 1])* s[i + 1])

to match the one in the actual model (line 151).
I forked it.

@muguangyuze

This comment has been minimized.

Show comment Hide comment
@muguangyuze

muguangyuze Mar 27, 2015

nice! But i want to know weather can you give a example for the newer.

nice! But i want to know weather can you give a example for the newer.

@arturoramos

This comment has been minimized.

Show comment Hide comment
@arturoramos

arturoramos Mar 31, 2015

Can someone explain the logic of this? This is the initialization of the trend in the cases of additive or multiplicative but the division by m squared doesn't make sense. Should it not be m_2 (m times two) instead of m_*2 (m squared)?

b = [(sum(Y[m:2 * m]) - sum(Y[0:m])) / float(m ** 2)]

Can someone explain the logic of this? This is the initialization of the trend in the cases of additive or multiplicative but the division by m squared doesn't make sense. Should it not be m_2 (m times two) instead of m_*2 (m squared)?

b = [(sum(Y[m:2 * m]) - sum(Y[0:m])) / float(m ** 2)]

@denis-bz

This comment has been minimized.

Show comment Hide comment
@denis-bz

denis-bz Apr 8, 2015

For a small self-contained example of Holt-Winters on synthetic ECG data,
see https://gist.github.com/denis-bz/85795665a261c14e5c05
It's NOT based on Andre's holtwinters.py, no optimization at all, but may help understanding.

(Bytheway, from __future__ import division makes e.g. 1/2 == 0.5, not 0, so you hardly ever need float() .)

denis-bz commented Apr 8, 2015

For a small self-contained example of Holt-Winters on synthetic ECG data,
see https://gist.github.com/denis-bz/85795665a261c14e5c05
It's NOT based on Andre's holtwinters.py, no optimization at all, but may help understanding.

(Bytheway, from __future__ import division makes e.g. 1/2 == 0.5, not 0, so you hardly ever need float() .)

@andrequeiroz

This comment has been minimized.

Show comment Hide comment
@andrequeiroz

andrequeiroz Jun 12, 2015

Thanks everyone for feedback. I did all the suggested corrections. :)

Owner

andrequeiroz commented Jun 12, 2015

Thanks everyone for feedback. I did all the suggested corrections. :)

@flipdazed

This comment has been minimized.

Show comment Hide comment
@flipdazed

flipdazed Jun 28, 2015

Nice, saved me a load of effort! I like the use of L_BFGS too.

Nice, saved me a load of effort! I like the use of L_BFGS too.

@rsnape

This comment has been minimized.

Show comment Hide comment
@rsnape

rsnape Jul 7, 2015

Hi, I mentioned this implementation in a stackoverflow answer here - thanks for the code! I notice as well that there has been interest in incorporating it into statsmodels here. I was wondering about the license on your code - can you (do you want to?) license it under a free software license?

rsnape commented Jul 7, 2015

Hi, I mentioned this implementation in a stackoverflow answer here - thanks for the code! I notice as well that there has been interest in incorporating it into statsmodels here. I was wondering about the license on your code - can you (do you want to?) license it under a free software license?

@andrequeiroz

This comment has been minimized.

Show comment Hide comment
@andrequeiroz

andrequeiroz Aug 5, 2015

Hi @rsnape I'm glad this code has been proved itself useful. I licensed it under the MIT license, so anyone can use it or modify it according to their needs. :)

Owner

andrequeiroz commented Aug 5, 2015

Hi @rsnape I'm glad this code has been proved itself useful. I licensed it under the MIT license, so anyone can use it or modify it according to their needs. :)

@kannbaba

This comment has been minimized.

Show comment Hide comment
@kannbaba

kannbaba Aug 22, 2015

hi denis-bz

your link https://gist.github.com/denis-bz/85795665a261c14e5c05 is broken. appreciate if you can provide another. thnx in advance.

hi denis-bz

your link https://gist.github.com/denis-bz/85795665a261c14e5c05 is broken. appreciate if you can provide another. thnx in advance.

@rsvp

This comment has been minimized.

Show comment Hide comment
@rsvp

rsvp Sep 21, 2015

Here's Holt-Winters implemented for IPython / Jupyter notebook / pandas: https://github.com/rsvp/fecon235/blob/master/lib/yi_timeseries.py -- which uses numpy and also handles DataFrame argument. It's Python 2 and 3 compatible, so thus it will work under both Python kernels for Jupyter.

We use it for financial economics, see https://git.io/fecon235 for topics which use Holt-Winters for smoothing and forecasting. Its license is BSD.

Robust parameters are discussed in: Sarah Gelper, 2007, Robust Forecasting with Exponential and Holt-Winters smoothing -- which are adopted as default parameters. However, the alpha and beta can be optimized given specific data, so we use scipy with a modifiable loss function (currently, not RMSE, but rather: median absolute error for robustness), see https://github.com/rsvp/fecon235/blob/master/lib/ys_opt_holt.py

The tests for the modules show numerically how the various functions are validated. Error bounds are included in h-step ahead forecast().

rsvp commented Sep 21, 2015

Here's Holt-Winters implemented for IPython / Jupyter notebook / pandas: https://github.com/rsvp/fecon235/blob/master/lib/yi_timeseries.py -- which uses numpy and also handles DataFrame argument. It's Python 2 and 3 compatible, so thus it will work under both Python kernels for Jupyter.

We use it for financial economics, see https://git.io/fecon235 for topics which use Holt-Winters for smoothing and forecasting. Its license is BSD.

Robust parameters are discussed in: Sarah Gelper, 2007, Robust Forecasting with Exponential and Holt-Winters smoothing -- which are adopted as default parameters. However, the alpha and beta can be optimized given specific data, so we use scipy with a modifiable loss function (currently, not RMSE, but rather: median absolute error for robustness), see https://github.com/rsvp/fecon235/blob/master/lib/ys_opt_holt.py

The tests for the modules show numerically how the various functions are validated. Error bounds are included in h-step ahead forecast().

@surfer26th

This comment has been minimized.

Show comment Hide comment
@surfer26th

surfer26th Jan 28, 2016

Hi Andre,

I believe that s[i] should be changed to s[i-m+1] in lines 69, 71, 81, 83, 149, 151, 182, 184 refer to Hyndman otext https://www.otexts.org/fpp/7/5

Hi Andre,

I believe that s[i] should be changed to s[i-m+1] in lines 69, 71, 81, 83, 149, 151, 182, 184 refer to Hyndman otext https://www.otexts.org/fpp/7/5

@welch

This comment has been minimized.

Show comment Hide comment
@welch

welch Mar 8, 2016

for a more robust initialization of the seasonal and trend components, there's seasonal

welch commented Mar 8, 2016

for a more robust initialization of the seasonal and trend components, there's seasonal

@cici7941

This comment has been minimized.

Show comment Hide comment
@cici7941

cici7941 Jul 8, 2016

Hi Andre,

I just tried using your optimizing approach to initialize my params for additive models and unfortunately it always gives me params like alpha=1, beta=0, gamma=0. I'm guessing it happens because in this way, it fitted your training dataset perfectly. But for forecasting, since you will no longer have actual data points, you will definitely need beta and gamma to be nozeros. Have you ever thought about it?

Thanks

cici7941 commented Jul 8, 2016

Hi Andre,

I just tried using your optimizing approach to initialize my params for additive models and unfortunately it always gives me params like alpha=1, beta=0, gamma=0. I'm guessing it happens because in this way, it fitted your training dataset perfectly. But for forecasting, since you will no longer have actual data points, you will definitely need beta and gamma to be nozeros. Have you ever thought about it?

Thanks

@anshumanksharma

This comment has been minimized.

Show comment Hide comment
@anshumanksharma

anshumanksharma May 2, 2017

Y = args[0] type = args[1] rmse = 0

What does Y stand for?

Y = args[0] type = args[1] rmse = 0

What does Y stand for?

@sinmmd

This comment has been minimized.

Show comment Hide comment
@sinmmd

sinmmd Jun 2, 2017

I tried using your implementation for multiplicative models and passed Y=(2,3,4,5,6,3,4,5,76,7,4,3), tupl=(Y,'multiplicative',1) and multiplicative(y,4,1) however there is an indexerror= tuple index out of range for a.append(alpha * (Y[i] / s[i-m+1]) + (1 - alpha) * (a[i] + b[i]))

Am I missing something?

sinmmd commented Jun 2, 2017

I tried using your implementation for multiplicative models and passed Y=(2,3,4,5,6,3,4,5,76,7,4,3), tupl=(Y,'multiplicative',1) and multiplicative(y,4,1) however there is an indexerror= tuple index out of range for a.append(alpha * (Y[i] / s[i-m+1]) + (1 - alpha) * (a[i] + b[i]))

Am I missing something?

@Lammatian

This comment has been minimized.

Show comment Hide comment
@Lammatian

Lammatian Aug 25, 2017

Hello,

I haven't tried using the code yet, but I have problem understanding why do the methods return Y[-fc]. Isn't Y just the data we have passed (with one appended value for whatever reason?) and y the prediction we are making? If that is the case, shouldn't y[-fc:] be the returned array? Any help appreciated.

Hello,

I haven't tried using the code yet, but I have problem understanding why do the methods return Y[-fc]. Isn't Y just the data we have passed (with one appended value for whatever reason?) and y the prediction we are making? If that is the case, shouldn't y[-fc:] be the returned array? Any help appreciated.

@ValantisTsiakos

This comment has been minimized.

Show comment Hide comment
@ValantisTsiakos

ValantisTsiakos Sep 26, 2017

@andrequeiroz Hi, i tried to use the code but i was not succesful. Could you please explain further the variables needed for each function ?
Use case: I am trying to predict the values of a raster based on the values of raster files of previous time periods.

Thanks in advance!!!

@andrequeiroz Hi, i tried to use the code but i was not succesful. Could you please explain further the variables needed for each function ?
Use case: I am trying to predict the values of a raster based on the values of raster files of previous time periods.

Thanks in advance!!!

@Daijunxu

This comment has been minimized.

Show comment Hide comment
@Daijunxu

Daijunxu Nov 25, 2017

Also fmin_l_bfgs_b is deprecated according to Scipy doc

The specific optimization method interfaces below in this subsection are not recommended for use in new scripts; all of these methods are accessible via a newer, more consistent interface provided by the functions above.

Also fmin_l_bfgs_b is deprecated according to Scipy doc

The specific optimization method interfaces below in this subsection are not recommended for use in new scripts; all of these methods are accessible via a newer, more consistent interface provided by the functions above.

@Vniex

This comment has been minimized.

Show comment Hide comment
@Vniex

Vniex Nov 27, 2017

Hello,
This code is great! But I have some questions, the form in code " s.append(gamma * (Y[i] / (a[i] + b[i])) + (1 - gamma) * s[i]) " , in some paper and wikipedia the form is " gamma * (Y[i] / a[i] ) + (1 - gamma) * s[i] ) ", are both forms all ok? Am I missing something?

Vniex commented Nov 27, 2017

Hello,
This code is great! But I have some questions, the form in code " s.append(gamma * (Y[i] / (a[i] + b[i])) + (1 - gamma) * s[i]) " , in some paper and wikipedia the form is " gamma * (Y[i] / a[i] ) + (1 - gamma) * s[i] ) ", are both forms all ok? Am I missing something?

@qx0731

This comment has been minimized.

Show comment Hide comment
@qx0731

qx0731 Mar 22, 2018

@arturoramos, that is the average of the slope , there is y(m+1) -y(1)/m for slope, then average those need divide another m

qx0731 commented Mar 22, 2018

@arturoramos, that is the average of the slope , there is y(m+1) -y(1)/m for slope, then average those need divide another m

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment