Tutorial

Let us consider chapter 7 of the excellent treatise on the subject of Exponential Smoothing By Hyndman and Athanasopoulos [1]. We will work through all the examples in the chapter as they unfold.

[1] Hyndman, Rob J., and George Athanasopoulos. Forecasting: principles and practice. OTexts, 2014.

Exponential smoothing

First we load some data. We have included the R data in the notebook for expedience.

In [1]:
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.api import ExponentialSmoothing, SimpleExpSmoothing, Holt

data = [446.6565,  454.4733,  455.663 ,  423.6322,  456.2713,  440.5881, 425.3325,  485.1494,  506.0482,  526.792 ,  514.2689,  494.211 ]
index= pd.DatetimeIndex(start='1996', end='2008', freq='A')
oildata = pd.Series(data, index)
oildata.index = pd.DatetimeIndex(oildata.index, freq=pd.infer_freq(oildata.index))

data = [17.5534,  21.86  ,  23.8866,  26.9293,  26.8885,  28.8314, 30.0751,  30.9535,  30.1857,  31.5797,  32.5776,  33.4774, 39.0216,  41.3864,  41.5966]
index= pd.DatetimeIndex(start='1990', end='2005', freq='A')
air = pd.Series(data, index)
air.index = pd.DatetimeIndex(air.index, freq=pd.infer_freq(air.index))

data = [263.9177,  268.3072,  260.6626,  266.6394,  277.5158,  283.834 , 290.309 ,  292.4742,  300.8307,  309.2867,  318.3311,  329.3724, 338.884 ,  339.2441,  328.6006,  314.2554,  314.4597,  321.4138, 329.7893,  346.3852,  352.2979,  348.3705,  417.5629,  417.1236, 417.7495,  412.2339,  411.9468,  394.6971,  401.4993,  408.2705, 414.2428]
index= pd.DatetimeIndex(start='1970', end='2001', freq='A')
livestock2 = pd.Series(data, index)
livestock2.index = pd.DatetimeIndex(livestock2.index, freq=pd.infer_freq(livestock2.index))

data = [407.9979 ,  403.4608,  413.8249,  428.105 ,  445.3387,  452.9942, 455.7402]
index= pd.DatetimeIndex(start='2001', end='2008', freq='A')
livestock3 = pd.Series(data, index)
livestock3.index = pd.DatetimeIndex(livestock3.index, freq=pd.infer_freq(livestock3.index))

data = [41.7275,  24.0418,  32.3281,  37.3287,  46.2132,  29.3463, 36.4829,  42.9777,  48.9015,  31.1802,  37.7179,  40.4202, 51.2069,  31.8872,  40.9783,  43.7725,  55.5586,  33.8509, 42.0764,  45.6423,  59.7668,  35.1919,  44.3197,  47.9137]
index= pd.DatetimeIndex(start='2005', end='2010-Q4', freq='QS')
aust = pd.Series(data, index)
aust.index = pd.DatetimeIndex(aust.index, freq=pd.infer_freq(aust.index))

Simple Exponential Smoothing

Lets use Simple Exponential Smoothing to forecast the below oil data.

In [2]:
ax=oildata.plot()
ax.set_xlabel("Year")
ax.set_ylabel("Oil (millions of tonnes)")
plt.show()
print("Figure 7.1: Oil production in Saudi Arabia from 1996 to 2007.")
Figure 7.1: Oil production in Saudi Arabia from 1996 to 2007.

Here we run three variants of simple exponential smoothing:

  1. In fit1 we do not use the auto optimization but instead choose to explicitly provide the model with the $\alpha=0.2$ parameter
  2. In fit2 as above we choose an $\alpha=0.6$
  3. In fit3 we allow statsmodels to automatically find an optimized $\alpha$ value for us. This is the recommended approach.
In [3]:
fit1 = SimpleExpSmoothing(oildata).fit(smoothing_level=0.2,optimized=False)
fcast1 = fit1.forecast(3).rename(r'$\alpha=0.2$')
fit2 = SimpleExpSmoothing(oildata).fit(smoothing_level=0.6,optimized=False)
fcast2 = fit2.forecast(3).rename(r'$\alpha=0.6$')
fit3 = SimpleExpSmoothing(oildata).fit()
fcast3 = fit3.forecast(3).rename(r'$\alpha=%s$'%fit3.model.params['smoothing_level'])

ax = oildata.plot(marker='o', color='black', figsize=(12,8))
fcast1.plot(marker='o', ax=ax, color='blue', legend=True)
fit1.fittedvalues.plot(marker='o', ax=ax, color='blue')
fcast2.plot(marker='o', ax=ax, color='red', legend=True)

fit2.fittedvalues.plot(marker='o', ax=ax, color='red')
fcast3.plot(marker='o', ax=ax, color='green', legend=True)
fit3.fittedvalues.plot(marker='o', ax=ax, color='green')
plt.show()

Holt's Method

Lets take a look at another example. This time we use air pollution data and the Holt's Method. We will fit three examples again.

  1. In fit1 we again choose not to use the optimzer and provide explicit values for $\alpha=0.8$ and $\beta=0.2$
  2. In fit2 we do the same as in fit1 but choose to use an exponential model rather than a Holt's additive model.
  3. In fit3 we used a damped versions of the Holt's additive model but allow the dampening parameter $\phi$ to be optimized while fixing the values for $\alpha=0.8$ and $\beta=0.2$
In [4]:
fit1 = Holt(air).fit(smoothing_level=0.8, smoothing_slope=0.2, optimized=False)
fcast1 = fit1.forecast(5).rename("Holt's linear trend")
fit2 = Holt(air, exponential=True).fit(smoothing_level=0.8, smoothing_slope=0.2, optimized=False)
fcast2 = fit2.forecast(5).rename("Exponential trend")
fit3 = Holt(air, damped=True).fit(smoothing_level=0.8, smoothing_slope=0.2)
fcast3 = fit3.forecast(5).rename("Additive damped trend")

ax = air.plot(color="black", marker="o", figsize=(12,8))
fit1.fittedvalues.plot(ax=ax, color='blue')
fcast1.plot(ax=ax, color='blue', marker="o", legend=True)
fit2.fittedvalues.plot(ax=ax, color='red')
fcast2.plot(ax=ax, color='red', marker="o", legend=True)
fit3.fittedvalues.plot(ax=ax, color='green')
fcast3.plot(ax=ax, color='green', marker="o", legend=True)

plt.show()

Seasonally adjusted data

Lets look at some seasonally adjusted livestock data. We fit five Holt's models. The below table allows us to compare results when we use exponential versus additive and damped versus non-damped.

Note: fit4 does not allow the parameter $\phi$ to be otpimized by providing a fixed value of $\phi=0.98$

In [5]:
fit1 = SimpleExpSmoothing(livestock2).fit()
fit2 = Holt(livestock2).fit()
fit3 = Holt(livestock2,exponential=True).fit()
fit4 = Holt(livestock2,damped=True).fit(damping_slope=0.98)
fit5 = Holt(livestock2,exponential=True,damped=True).fit()
params = ['smoothing_level', 'smoothing_slope', 'damping_slope', 'initial_level', 'initial_slope']
results=pd.DataFrame(index=[r"$\alpha$",r"$\beta$",r"$\phi$",r"$l_0$","$b_0$","SSE"] ,columns=['SES', "Holt's","Exponential", "Additive", "Multiplicative"])
results["SES"] =            [fit1.params[p] for p in params] + [fit1.sse]
results["Holt's"] =         [fit2.params[p] for p in params] + [fit2.sse]
results["Exponential"] =    [fit3.params[p] for p in params] + [fit3.sse]
results["Additive"] =       [fit4.params[p] for p in params] + [fit4.sse]
results["Multiplicative"] = [fit5.params[p] for p in params] + [fit5.sse]
results
Out[5]:
SES Holt's Exponential Additive Multiplicative
$\alpha$ 1.000000 0.974306 0.977633 0.978848 0.974911
$\beta$ NaN 0.000000 0.000000 0.000000 0.000000
$\phi$ NaN NaN NaN 0.980000 0.981646
$l_0$ 263.917700 258.882590 260.341513 257.357600 258.951854
$b_0$ NaN 5.010792 1.013780 6.511648 1.019090
SSE 6761.350218 6004.138200 6104.194746 6036.555004 6081.995045

Plots of Seasonally Adjusted Data

The following plots allow us to evaluate the level and slope/trend components of the above table's fits.

In [6]:
for fit in [fit2,fit4]:
    pd.DataFrame(np.c_[fit.level,fit.slope]).rename(
        columns={0:'level',1:'slope'}).plot(subplots=True)
plt.show()
print('Figure 7.4: Level and slope components for Holt’s linear trend method and the additive damped trend method.')
Figure 7.4: Level and slope components for Holt’s linear trend method and the additive damped trend method.

Comparison

Here we plot a comparison Simple Exponential Smoothing and Holt's Methods for various additive, exponential and damped combinations. All of the models parameters will be optimized by statsmodels.

In [7]:
fit1 = SimpleExpSmoothing(livestock2).fit()
fcast1 = fit1.forecast(9).rename("SES")
fit2 = Holt(livestock2).fit()
fcast2 = fit2.forecast(9).rename("Holt's")
fit3 = Holt(livestock2, exponential=True).fit()
fcast3 = fit3.forecast(9).rename("Exponential")
fit4 = Holt(livestock2, damped=True).fit(damping_slope=0.98)
fcast4 = fit4.forecast(9).rename("Additive Damped")
fit5 = Holt(livestock2, exponential=True, damped=True).fit()
fcast5 = fit5.forecast(9).rename("Multiplicative Damped")

ax = livestock2.plot(color="black", marker="o", figsize=(12,8))
livestock3.plot(ax=ax, color="black", marker="o", legend=False)
fcast1.plot(ax=ax, color='red', legend=True)
fcast2.plot(ax=ax, color='green', legend=True)
fcast3.plot(ax=ax, color='blue', legend=True)
fcast4.plot(ax=ax, color='cyan', legend=True)
fcast5.plot(ax=ax, color='magenta', legend=True)
ax.set_ylabel('Livestock, sheep in Asia (millions)')
plt.show()
print('Figure 7.5: Forecasting livestock, sheep in Asia: comparing forecasting performance of non-seasonal methods.')
Figure 7.5: Forecasting livestock, sheep in Asia: comparing forecasting performance of non-seasonal methods.

Holt's Winters Seasonal

Finally we are able to run full Holt's Winters Seasonal Exponential Smoothing including a trend component and a seasonal component. Statsmodels allows for all the combinations including as shown in the examples below:

  1. fit1 additive trend, additive seasonal of period season_length=4 and the use of a Boxcox transformation.
  2. fit2 additive trend, multiplicative seasonal of period season_length=4 and the use of a Boxcox transformation..
  3. fit3 additive damped trend, additive seasonal of period season_length=4 and the use of a boxcox transformation.
  4. fit4 additive damped trend, multiplicative seasonal of period season_length=4 and the use of a boxcox transformation.

The plot shows the results and forecast for fit1 and fit2. The table allows us to compare the results and parameterizations.

In [8]:
fit1 = ExponentialSmoothing(aust, seasonal_periods=4, trend='add', seasonal='add').fit(use_boxcox=True)
fit2 = ExponentialSmoothing(aust, seasonal_periods=4, trend='add', seasonal='mul').fit(use_boxcox=True)
fit3 = ExponentialSmoothing(aust, seasonal_periods=4, trend='add', seasonal='add', damped=True).fit(use_boxcox=True)
fit4 = ExponentialSmoothing(aust, seasonal_periods=4, trend='add', seasonal='mul', damped=True).fit(use_boxcox=True)
results=pd.DataFrame(index=[r"$\alpha$",r"$\beta$",r"$\phi$",r"$\gamma$",r"$l_0$","$b_0$","SSE"])
params = ['smoothing_level', 'smoothing_slope', 'damping_slope', 'smoothing_seasonal', 'initial_level', 'initial_slope']
results["Additive"]       = [fit1.params[p] for p in params] + [fit1.sse]
results["Multiplicative"] = [fit2.params[p] for p in params] + [fit2.sse]
results["Additive Dam"]   = [fit3.params[p] for p in params] + [fit3.sse]
results["Multiplica Dam"] = [fit4.params[p] for p in params] + [fit4.sse]

ax = aust.plot(figsize=(10,6), marker='o', color='black', title="Forecasts from Holt-Winters' multiplicative method" )
ax.set_ylabel("International visitor night in Australia (millions)")
ax.set_xlabel("Year")
fit1.fittedvalues.plot(ax=ax, style='--', color='red')
fit2.fittedvalues.plot(ax=ax, style='--', color='green')

fit1.forecast(8).plot(ax=ax, style='--', marker='o', color='red', legend=True)
fit2.forecast(8).plot(ax=ax, style='--', marker='o', color='green', legend=True)

plt.show()
print("Figure 7.6: Forecasting international visitor nights in Australia using Holt-Winters method with both additive and multiplicative seasonality.")

results
Figure 7.6: Forecasting international visitor nights in Australia using Holt-Winters method with both additive and multiplicative seasonality.
Out[8]:
Additive Multiplicative Additive Dam Multiplica Dam
$\alpha$ 4.546206e-01 3.658961e-01 9.514088e-09 0.000181
$\beta$ 8.173643e-09 1.116099e-21 3.276534e-73 0.000181
$\phi$ NaN NaN 9.435278e-01 0.912994
$\gamma$ 5.244301e-01 1.614569e-15 1.139715e-06 0.000000
$l_0$ 1.421750e+01 1.454957e+01 1.416110e+01 14.534939
$b_0$ 1.307797e-01 1.661715e-01 2.305158e-01 0.444003
SSE 5.001650e+01 4.307077e+01 3.525563e+01 39.839183

The Internals

It is possible to get at the internals of the Exponential Smoothing models.

Here we show some tables that allow you to view side by side the original values $y_t$, the level $l_t$, the trend $b_t$, the season $s_t$ and the fitted values $\hat{y}_t$.

In [9]:
pd.DataFrame(np.c_[aust, fit1.level, fit1.slope, fit1.season, fit1.fittedvalues], 
columns=[r'$y_t$',r'$l_t$',r'$b_t$',r'$s_t$',r'$\hat{y}_t$'],index=aust.index). \
append(fit1.forecast(8).rename(r'$\hat{y}_t$').to_frame())
/usr/lib/python3/dist-packages/pandas/core/frame.py:6211: FutureWarning: Sorting because non-concatenation axis is not aligned. A future version
of pandas will change to not sort by default.

To accept the future behavior, pass 'sort=False'.

To retain the current behavior and silence the warning, pass 'sort=True'.

  sort=sort)
Out[9]:
$\hat{y}_t$ $b_t$ $l_t$ $s_t$ $y_t$
2005-01-01 41.721183 -34.969345 49.317629 -7.593360 41.7275
2005-04-01 24.190165 -35.452844 49.932481 -25.834735 24.0418
2005-07-01 31.460535 -36.532820 51.126193 -19.182999 32.3281
2005-10-01 36.634766 -37.397670 52.210144 -15.209749 37.3287
2006-01-01 45.097765 -38.467245 53.476765 -7.837217 46.2132
2006-04-01 27.191742 -40.276302 55.513849 -27.017380 29.3463
2006-07-01 36.544223 -40.625049 56.224488 -19.713866 36.4829
2006-10-01 41.449476 -42.041771 57.766101 -15.523340 42.9777
2007-01-01 50.934531 -41.543814 57.536753 -7.588665 48.9015
2007-04-01 31.418274 -42.197843 58.151018 -26.874278 31.1802
2007-07-01 38.718335 -42.303505 58.363026 -20.191946 37.7179
2007-10-01 44.140672 -41.085809 57.181922 -14.982807 40.4202
2008-01-01 49.315795 -42.961527 58.853048 -8.619915 51.2069
2008-04-01 32.307038 -43.186572 59.367040 -27.309207 31.8872
2008-07-01 39.207470 -44.827035 61.095659 -20.925630 40.9783
2008-10-01 44.551251 -44.899536 61.462249 -17.319996 43.7725
2009-01-01 54.358160 -46.192273 62.816889 -7.881677 55.5586
2009-04-01 35.153876 -45.980227 62.832220 -28.447936 33.8509
2009-07-01 43.066536 -46.228158 63.082712 -20.550689 42.0764
2009-10-01 45.871169 -46.852570 63.748932 -17.997895 45.6423
2010-01-01 57.166668 -48.770444 65.777651 -7.372181 59.7668
2010-04-01 36.761385 -48.308396 65.650095 -29.816914 35.1919
2010-07-01 44.932522 -48.798715 66.119516 -21.517551 44.3197
2010-10-01 48.399596 -49.269879 66.667503 -18.522357 47.9137
2011-01-01 61.338048 NaN NaN NaN NaN
2011-04-01 37.242899 NaN NaN NaN NaN
2011-07-01 46.842699 NaN NaN NaN NaN
2011-10-01 51.005341 NaN NaN NaN NaN
2012-01-01 64.470996 NaN NaN NaN NaN
2012-04-01 39.777024 NaN NaN NaN NaN
2012-07-01 49.636014 NaN NaN NaN NaN
2012-10-01 53.901641 NaN NaN NaN NaN
In [10]:
pd.DataFrame(np.c_[aust, fit2.level, fit2.slope, fit2.season, fit2.fittedvalues], 
columns=[r'$y_t$',r'$l_t$',r'$b_t$',r'$s_t$',r'$\hat{y}_t$'],index=aust.index). \
append(fit2.forecast(8).rename(r'$\hat{y}_t$').to_frame())
/usr/lib/python3/dist-packages/pandas/core/frame.py:6211: FutureWarning: Sorting because non-concatenation axis is not aligned. A future version
of pandas will change to not sort by default.

To accept the future behavior, pass 'sort=False'.

To retain the current behavior and silence the warning, pass 'sort=True'.

  sort=sort)
Out[10]:
$\hat{y}_t$ $b_t$ $l_t$ $s_t$ $y_t$
2005-01-01 41.859741 -36.536816 51.252556 0.815776 41.7275
2005-04-01 25.838629 -35.873460 50.744389 0.495306 24.0418
2005-07-01 31.659691 -37.290590 52.068620 0.612908 32.3281
2005-10-01 35.189842 -39.176092 54.195378 0.664094 37.3287
2006-01-01 44.928996 -40.313350 55.714621 0.814935 46.2132
2006-04-01 27.933829 -42.095035 57.765220 0.492986 29.3463
2006-07-01 35.824846 -43.108031 59.136209 0.610014 36.4829
2006-10-01 39.768910 -45.651969 61.916436 0.661617 42.9777
2007-01-01 51.174368 -45.129560 61.866043 0.813475 48.9015
2007-04-01 30.814710 -46.417035 63.145045 0.490231 31.1802
2007-07-01 39.009742 -46.395029 63.337175 0.608146 37.7179
2007-10-01 42.486372 -46.179425 63.153575 0.660354 40.4202
2008-01-01 52.174331 -46.768283 63.711851 0.813268 51.2069
2008-04-01 31.677624 -47.844624 64.881067 0.489487 31.8872
2008-07-01 40.036382 -49.248297 66.478119 0.607596 40.9783
2008-10-01 44.516283 -49.584255 67.075800 0.659493 43.7725
2009-01-01 55.343426 -50.611703 68.200495 0.812650 55.5586
2009-04-01 33.773408 -51.524884 69.295663 0.487812 33.8509
2009-07-01 42.644872 -52.034901 69.981645 0.606295 42.0764
2009-10-01 46.778801 -52.320367 70.376720 0.658604 45.6423
2010-01-01 58.009207 -54.103871 72.223140 0.812172 59.7668
2010-04-01 35.648665 -54.510016 72.921325 0.486452 35.1919
2010-07-01 44.785013 -55.173877 73.694772 0.605321 44.3197
2010-10-01 49.174847 -55.399718 74.041485 0.657736 47.9137
2011-01-01 60.967531 NaN NaN NaN NaN
2011-04-01 36.994137 NaN NaN NaN NaN
2011-07-01 46.713232 NaN NaN NaN NaN
2011-10-01 51.483555 NaN NaN NaN NaN
2012-01-01 64.457006 NaN NaN NaN NaN
2012-04-01 39.017904 NaN NaN NaN NaN
2012-07-01 49.293066 NaN NaN NaN NaN
2012-10-01 54.321401 NaN NaN NaN NaN

Finally lets look at the levels, slopes/trends and seasonal components of the models.

In [11]:
states1 = pd.DataFrame(np.c_[fit1.level, fit1.slope, fit1.season], columns=['level','slope','seasonal'], index=aust.index)
states2 = pd.DataFrame(np.c_[fit2.level, fit2.slope, fit2.season], columns=['level','slope','seasonal'], index=aust.index)
fig, [[ax1, ax4],[ax2, ax5], [ax3, ax6]] = plt.subplots(3, 2, figsize=(12,8))
states1[['level']].plot(ax=ax1)
states1[['slope']].plot(ax=ax2)
states1[['seasonal']].plot(ax=ax3)
states2[['level']].plot(ax=ax4)
states2[['slope']].plot(ax=ax5)
states2[['seasonal']].plot(ax=ax6)
plt.show()
In [ ]: