# Source code for statsmodels.tsa.holtwinters

"""
Notes
-----
Code written using below textbook as a reference.
Results are checked against the expected outcomes in the text book.

Properties:
Hyndman, Rob J., and George Athanasopoulos. Forecasting: principles and
practice. OTexts, 2014.

Author: Terence L van Zyl
Modified: Kevin Sheppard
"""
from statsmodels.compat.python import string_types

import numpy as np
import pandas as pd
from scipy.optimize import basinhopping, brute, minimize
from scipy.spatial.distance import sqeuclidean
from scipy.special import inv_boxcox
from scipy.stats import boxcox

from statsmodels.base.model import Results
from statsmodels.base.wrapper import populate_wrapper, union_dicts, ResultsWrapper
from statsmodels.tsa.base.tsa_model import TimeSeriesModel
from statsmodels.tsa.tsatools import freq_to_period
import statsmodels.tsa._exponential_smoothers as smoothers

def _holt_init(x, xi, p, y, l, b):
"""Initialization for the Holt Models"""
p[xi.astype(np.bool)] = x
alpha, beta, _, l0, b0, phi = p[:6]
alphac = 1 - alpha
betac = 1 - beta
y_alpha = alpha * y
l[:] = 0
b[:] = 0
l[0] = l0
b[0] = b0
return alpha, beta, phi, alphac, betac, y_alpha

def _holt__(x, xi, p, y, l, b, s, m, n, max_seen):
"""
Simple Exponential Smoothing
Minimization Function
(,)
"""
alpha, beta, phi, alphac, betac, y_alpha = _holt_init(x, xi, p, y, l, b)
for i in range(1, n):
l[i] = (y_alpha[i - 1]) + (alphac * (l[i - 1]))
return sqeuclidean(l, y)

def _holt_mul_dam(x, xi, p, y, l, b, s, m, n, max_seen):
"""
Multiplicative and Multiplicative Damped
Minimization Function
(M,) & (Md,)
"""
alpha, beta, phi, alphac, betac, y_alpha = _holt_init(x, xi, p, y, l, b)
if alpha == 0.0:
return max_seen
if beta > alpha:
return max_seen
for i in range(1, n):
l[i] = (y_alpha[i - 1]) + (alphac * (l[i - 1] * b[i - 1]**phi))
b[i] = (beta * (l[i] / l[i - 1])) + (betac * b[i - 1]**phi)
return sqeuclidean(l * b**phi, y)

def _holt_add_dam(x, xi, p, y, l, b, s, m, n, max_seen):
"""
Minimization Function
"""
alpha, beta, phi, alphac, betac, y_alpha = _holt_init(x, xi, p, y, l, b)
if alpha == 0.0:
return max_seen
if beta > alpha:
return max_seen
for i in range(1, n):
l[i] = (y_alpha[i - 1]) + (alphac * (l[i - 1] + phi * b[i - 1]))
b[i] = (beta * (l[i] - l[i - 1])) + (betac * phi * b[i - 1])
return sqeuclidean(l + phi * b, y)

def _holt_win_init(x, xi, p, y, l, b, s, m):
"""Initialization for the Holt Winters Seasonal Models"""
p[xi.astype(np.bool)] = x
alpha, beta, gamma, l0, b0, phi = p[:6]
s0 = p[6:]
alphac = 1 - alpha
betac = 1 - beta
gammac = 1 - gamma
y_alpha = alpha * y
y_gamma = gamma * y
l[:] = 0
b[:] = 0
s[:] = 0
l[0] = l0
b[0] = b0
s[:m] = s0
return alpha, beta, gamma, phi, alphac, betac, gammac, y_alpha, y_gamma

def _holt_win__mul(x, xi, p, y, l, b, s, m, n, max_seen):
"""
Multiplicative Seasonal
Minimization Function
(,M)
"""
alpha, beta, gamma, phi, alphac, betac, gammac, y_alpha, y_gamma = _holt_win_init(
x, xi, p, y, l, b, s, m)
if alpha == 0.0:
return max_seen
if gamma > 1 - alpha:
return max_seen
for i in range(1, n):
l[i] = (y_alpha[i - 1] / s[i - 1]) + (alphac * (l[i - 1]))
s[i + m - 1] = (y_gamma[i - 1] / (l[i - 1])) + (gammac * s[i - 1])
return sqeuclidean(l * s[:-(m - 1)], y)

def _holt_win__add(x, xi, p, y, l, b, s, m, n, max_seen):
"""
Minimization Function
(,A)
"""
alpha, beta, gamma, phi, alphac, betac, gammac, y_alpha, y_gamma = _holt_win_init(
x, xi, p, y, l, b, s, m)
if alpha == 0.0:
return max_seen
if gamma > 1 - alpha:
return max_seen
for i in range(1, n):
l[i] = (y_alpha[i - 1]) - (alpha * s[i - 1]) + (alphac * (l[i - 1]))
s[i + m - 1] = y_gamma[i - 1] - (gamma * (l[i - 1])) + (gammac * s[i - 1])
return sqeuclidean(l + s[:-(m - 1)], y)

def _holt_win_add_mul_dam(x, xi, p, y, l, b, s, m, n, max_seen):
"""
Minimization Function
"""
alpha, beta, gamma, phi, alphac, betac, gammac, y_alpha, y_gamma = _holt_win_init(
x, xi, p, y, l, b, s, m)
if alpha * beta == 0.0:
return max_seen
if beta > alpha or gamma > 1 - alpha:
return max_seen
for i in range(1, n):
l[i] = (y_alpha[i - 1] / s[i - 1]) + \
(alphac * (l[i - 1] + phi * b[i - 1]))
b[i] = (beta * (l[i] - l[i - 1])) + (betac * phi * b[i - 1])
s[i + m - 1] = (y_gamma[i - 1] / (l[i - 1] + phi *
b[i - 1])) + (gammac * s[i - 1])
return sqeuclidean((l + phi * b) * s[:-(m - 1)], y)

def _holt_win_mul_mul_dam(x, xi, p, y, l, b, s, m, n, max_seen):
"""
Multiplicative and Multiplicative Damped with Multiplicative Seasonal
Minimization Function
(M,M) & (Md,M)
"""
alpha, beta, gamma, phi, alphac, betac, gammac, y_alpha, y_gamma = _holt_win_init(
x, xi, p, y, l, b, s, m)
if alpha * beta == 0.0:
return max_seen
if beta > alpha or gamma > 1 - alpha:
return max_seen
for i in range(1, n):
l[i] = (y_alpha[i - 1] / s[i - 1]) + \
(alphac * (l[i - 1] * b[i - 1]**phi))
b[i] = (beta * (l[i] / l[i - 1])) + (betac * b[i - 1]**phi)
s[i + m - 1] = (y_gamma[i - 1] / (l[i - 1] *
b[i - 1]**phi)) + (gammac * s[i - 1])
return sqeuclidean((l * b**phi) * s[:-(m - 1)], y)

"""
Minimization Function
"""
alpha, beta, gamma, phi, alphac, betac, gammac, y_alpha, y_gamma = _holt_win_init(
x, xi, p, y, l, b, s, m)
if alpha * beta == 0.0:
return max_seen
if beta > alpha or gamma > 1 - alpha:
return max_seen
for i in range(1, n):
l[i] = (y_alpha[i - 1]) - (alpha * s[i - 1]) + \
(alphac * (l[i - 1] + phi * b[i - 1]))
b[i] = (beta * (l[i] - l[i - 1])) + (betac * phi * b[i - 1])
s[i + m - 1] = y_gamma[i - 1] - (gamma * (l[i - 1] + phi * b[i - 1])) + (gammac * s[i - 1])
return sqeuclidean((l + phi * b) + s[:-(m - 1)], y)

def _holt_win_mul_add_dam(x, xi, p, y, l, b, s, m, n, max_seen):
"""
Multiplicative and Multiplicative Damped with Additive Seasonal
Minimization Function
"""
alpha, beta, gamma, phi, alphac, betac, gammac, y_alpha, y_gamma = _holt_win_init(
x, xi, p, y, l, b, s, m)
if alpha * beta == 0.0:
return max_seen
if beta > alpha or gamma > 1 - alpha:
return max_seen
for i in range(1, n):
l[i] = (y_alpha[i - 1]) - (alpha * s[i - 1]) + \
(alphac * (l[i - 1] * b[i - 1]**phi))
b[i] = (beta * (l[i] / l[i - 1])) + (betac * b[i - 1]**phi)
s[i + m - 1] = y_gamma[i - 1] - \
(gamma * (l[i - 1] * b[i - 1]**phi)) + (gammac * s[i - 1])
return sqeuclidean((l * phi * b) + s[:-(m - 1)], y)

('mul', 'mul'): smoothers._holt_win_mul_mul_dam,
('mul', None): smoothers._holt_win__mul,
(None, 'mul'): smoothers._holt_mul_dam,
(None, None): smoothers._holt__}

('mul', 'mul'): _holt_win_mul_mul_dam,
('mul', None): _holt_win__mul,
(None, 'mul'): _holt_mul_dam,
(None, None): _holt__}

[docs]class HoltWintersResults(Results): """ Holt Winter's Exponential Smoothing Results Parameters ---------- model : ExponentialSmoothing instance The fitted model instance params : dict All the parameters for the Exponential Smoothing model. Attributes ---------- params: dict All the parameters for the Exponential Smoothing model. params_formatted: pd.DataFrame DataFrame containing all parameters, their short names and a flag indicating whether the parameter's value was optimized to fit the data. fittedfcast: array An array of both the fitted values and forecast values. fittedvalues: array An array of the fitted values. Fitted by the Exponential Smoothing model. fcastvalues: array An array of the forecast values forecast by the Exponential Smoothing model. sse: float The sum of squared errors level: array An array of the levels values that make up the fitted values. slope: array An array of the slope values that make up the fitted values. season: array An array of the seasonal values that make up the fitted values. aic: float The Akaike information criterion. bic: float The Bayesian information criterion. aicc: float AIC with a correction for finite sample sizes. resid: array An array of the residuals of the fittedvalues and actual values. k: int the k parameter used to remove the bias in AIC, BIC etc. optimized: bool Flag indicating whether the model parameters were optimized to fit the data. mle_retvals: {None, scipy.optimize.optimize.OptimizeResult} Optimization results if the parameters were optimized to fit the data. """ def __init__(self, model, params, **kwargs): self.data = model.data super(HoltWintersResults, self).__init__(model, params, **kwargs)
[docs] def predict(self, start=None, end=None): """ In-sample prediction and out-of-sample forecasting Parameters ---------- start : int, str, or datetime, optional Zero-indexed observation number at which to start forecasting, ie., the first forecast is start. Can also be a date string to parse or a datetime type. Default is the the zeroth observation. end : int, str, or datetime, optional Zero-indexed observation number at which to end forecasting, ie., the first forecast is start. Can also be a date string to parse or a datetime type. However, if the dates index does not have a fixed frequency, end must be an integer index if you want out of sample prediction. Default is the last observation in the sample. Returns ------- forecast : array Array of out of sample forecasts. """ return self.model.predict(self.params, start, end)
[docs] def forecast(self, steps=1): """ Out-of-sample forecasts Parameters ---------- steps : int The number of out of sample forecasts from the end of the sample. Returns ------- forecast : array Array of out of sample forecasts """ try: freq = getattr(self.model._index, 'freq', 1) start = self.model._index[-1] + freq end = self.model._index[-1] + steps * freq return self.model.predict(self.params, start=start, end=end) except (AttributeError, ValueError): # May occur when the index doesn't have a freq return self.model._predict(h=steps, **self.params).fcastvalues
[docs] def summary(self): """ Summarize the fitted Model Returns ------- smry : Summary instance This holds the summary table and text, which can be printed or converted to various output formats. See Also -------- statsmodels.iolib.summary.Summary """ from statsmodels.iolib.summary import Summary from statsmodels.iolib.table import SimpleTable model = self.model title = model.__class__.__name__ + ' Model Results' dep_variable = 'endog' if isinstance(self.model.endog, pd.DataFrame): dep_variable = self.model.endog.columns[0] elif isinstance(self.model.endog, pd.Series): dep_variable = self.model.endog.name seasonal_periods = None if self.model.seasonal is None else self.model.seasonal_periods lookup = {'add': 'Additive', 'additive': 'Additive', 'mul': 'Multiplicative', 'multiplicative': 'Multiplicative', None: 'None'} transform = self.params['use_boxcox'] box_cox_transform = True if transform else False box_cox_coeff = transform if isinstance(transform, string_types) else self.params['lamda'] if isinstance(box_cox_coeff, float): box_cox_coeff = '{:>10.5f}'.format(box_cox_coeff) top_left = [('Dep. Variable:', [dep_variable]), ('Model:', [model.__class__.__name__]), ('Optimized:', [str(np.any(self.optimized))]), ('Trend:', [lookup[self.model.trend]]), ('Seasonal:', [lookup[self.model.seasonal]]), ('Seasonal Periods:', [str(seasonal_periods)]), ('Box-Cox:', [str(box_cox_transform)]), ('Box-Cox Coeff.:', [str(box_cox_coeff)])] top_right = [ ('No. Observations:', [str(len(self.model.endog))]), ('SSE', ['{:5.3f}'.format(self.sse)]), ('AIC', ['{:5.3f}'.format(self.aic)]), ('BIC', ['{:5.3f}'.format(self.bic)]), ('AICC', ['{:5.3f}'.format(self.aicc)]), ('Date:', None), ('Time:', None)] smry = Summary() smry.add_table_2cols(self, gleft=top_left, gright=top_right, title=title) formatted = self.params_formatted # type: pd.DataFrame def _fmt(x): abs_x = np.abs(x) scale = 1 if abs_x != 0: scale = int(np.log10(abs_x)) if scale > 4 or scale < -3: return '{:>20.5g}'.format(x) dec = min(7 - scale, 7) fmt = '{{:>20.{0}f}}'.format(dec) return fmt.format(x) tab = [] for _, vals in formatted.iterrows(): tab.append([_fmt(vals.iloc[1]), '{0:>20}'.format(vals.iloc[0]), '{0:>20}'.format(str(bool(vals.iloc[2])))]) params_table = SimpleTable(tab, headers=['coeff', 'code', 'optimized'], title="", stubs=list(formatted.index)) smry.tables.append(params_table) return smry
class HoltWintersResultsWrapper(ResultsWrapper): _attrs = {'fittedvalues': 'rows', 'level': 'rows', 'resid': 'rows', 'season': 'rows', 'slope': 'rows'} _wrap_attrs = union_dicts(ResultsWrapper._wrap_attrs, _attrs) _methods = {'predict': 'dates', 'forecast': 'dates'} _wrap_methods = union_dicts(ResultsWrapper._wrap_methods, _methods) populate_wrapper(HoltWintersResultsWrapper, HoltWintersResults)
[docs]class SimpleExpSmoothing(ExponentialSmoothing): """ Simple Exponential Smoothing Parameters ---------- endog : array-like Time series Returns ------- results : SimpleExpSmoothing class Notes ----- This is a full implementation of the simple exponential smoothing as per [1]_. SimpleExpSmoothing is a restricted version of :class:ExponentialSmoothing. See Also -------- ExponentialSmoothing Holt References ---------- .. [1] Hyndman, Rob J., and George Athanasopoulos. Forecasting: principles and practice. OTexts, 2014. """ def __init__(self, endog): super(SimpleExpSmoothing, self).__init__(endog)
[docs]class Holt(ExponentialSmoothing): """ Holt's Exponential Smoothing Parameters ---------- endog : array-like Time series exponential : bool, optional Type of trend component. damped : bool, optional Should the trend component be damped. Returns ------- results : Holt class Notes ----- This is a full implementation of the Holt's exponential smoothing as per [1]_. Holt is a restricted version of :class:ExponentialSmoothing. See Also -------- ExponentialSmoothing SimpleExpSmoothing References ---------- .. [1] Hyndman, Rob J., and George Athanasopoulos. Forecasting: principles and practice. OTexts, 2014. """ def __init__(self, endog, exponential=False, damped=False): trend = 'mul' if exponential else 'add' super(Holt, self).__init__(endog, trend=trend, damped=damped)