In this tutorial we will study how to apply exponential smoothing methods to forecast seasonal data in Python.
Data: NSW retail turnover
Time series decomposition
Trend corrected exponential smoothing
Holt-Winters smoothing
Model diagnostics
Model validation
Forecast
This notebook relies on the following imports and settings
# Packages
import warnings
warnings.filterwarnings("ignore")
import numpy as np
from scipy import stats
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import statsmodels.api as sm
# Plot settings
sns.set_context('notebook')
sns.set_style('ticks')
red='#D62728'
blue='#1F77B4'
%matplotlib inline
In this tutorial we will keep working with the Australian retail turnover series. The next cells load the data and reproduces some of the steps from last week.
data = pd.read_csv('nswretail.csv', index_col='Month', parse_dates=True, dayfirst=True)
data.tail()
Turnover | |
---|---|
Month | |
2017-02-01 | 7298.9 |
2017-03-01 | 8085.8 |
2017-04-01 | 7883.7 |
2017-05-01 | 8132.0 |
2017-06-01 | 8130.1 |
y = data['Turnover'].copy()
y.index = y.index.to_period(freq='M')
ts = data['Turnover']
data.describe().round(2)
Turnover | |
---|---|
count | 138.00 |
mean | 6611.42 |
std | 1168.73 |
min | 4496.90 |
25% | 5769.78 |
50% | 6347.20 |
75% | 7450.70 |
max | 10783.10 |
fig, ax= plt.subplots(figsize=(8,5))
y.plot(color=red)
ax.set_xlabel('')
ax.set_ylabel('Turnover')
ax.set_title('NSW retail turnover (2005-2017)')
ax.set_xticks([], minor=True)
sns.despine()
plt.show()
fig, ax= plt.subplots(figsize=(8,5))
np.log(y).plot(color=red)
ax.set_xlabel('')
ax.set_ylabel('Log turnover')
ax.set_title('Log series')
ax.set_xticks([], minor=True)
sns.despine()
plt.show()
We use the X-13 ARIMA-SEATS method for time series decomposition. The X-13 method is a state-of-art time series decomposition approach devoloped and provided by the US Census Bureau as external software. You need to make the X-13 software available to Python by downloading it from the link and extracting the compressed folder to your computer. The statsmodels library has a function that interfaces with this software, and the easiest way to get it to work is to copy the x13as executable to the same folder as your notebook.
The X-13 is the current version of the X-12-ARIMA method described in the textbook. The only limitation is that the X-13 software makes calendar effect adjustments based on the US calendar by default, though it is possible to reconfigure it to other calendars. Refer to the ABS for the seasonal adjustment standards in Australia.
decomposition = sm.tsa.x13_arima_analysis(ts)
You can use the dir method to check the available output.
dir(decomposition)[-6:]
['observed', 'plot', 'results', 'seasadj', 'stdout', 'trend']
For example, if we want to retrieve and plot the seasonally adjusted component, we can do as follows.
decomposition.seasadj.tail()
Month 2017-02-01 8263.077126 2017-03-01 8321.825080 2017-04-01 8334.278109 2017-05-01 8425.606055 2017-06-01 8507.854903 Name: seasadj, dtype: float64
fig, ax= plt.subplots(figsize=(8,5))
decomposition.seasadj.plot(color=red)
ax.set_xlabel('')
ax.set_ylabel('Turnover')
ax.set_title('Seasonally NSW retail turnover 2005-2017)')
ax.set_xticks([], minor=True) # I prefer to remove the minor ticks for a cleaner plot
sns.despine()
plt.show()
The forecast module has a convenient function for plotting the components.
from forecast import plot_components_x13
plot_components_x13(decomposition)
plt.show()
The X-13-ARIMA software generates a very detailed analysis of the time series, which is stored in the results attribute. You can view it with the print(decomposition.results) command, but due to the length I recommend saving it to a text file and opening it elsewhere. You can do this as follows.
f=open('x13results.txt', mode='w+')
f.write(decomposition.results)
f.close()
Feel free to use information from the X-13-ARIMA report in your assignment, as long as you can justify and explain it.
In case you are not able to work with the X13 decomposition software, a more basic decomposition is available as follows.
basic_decomposition = sm.tsa.seasonal_decompose(ts, model='multiplicative')
Before we model and forecast the original series, the next cells demonstrate the application of the trend corrected exponential smoothing method, focusing on the seasonally adjusted series for coherence.
ts=decomposition.seasadj
ts.index=ts.index.to_period(freq='M')
import forecast
holt = forecast.holt(ts)
holt.fit()
holt.summary()
Holt (trend corrected) exponential smoothing Smoothing parameters: alpha (level) 0.584 (0.077) beta (trend) 0.042 (0.028) In-sample fit: MSE 6467.623 Log-likelihood -801.258 AIC 1608.517 BIC 1617.299
smoothed=pd.Series(holt.smooth(), index=y.index)
fig, ax= plt.subplots(figsize=(10,6))
ts.plot(color='#D62728', label='Seasonally adjusted retail turnover')
smoothed.plot(color='black', label='Smoothed', alpha=0.65)
plt.legend(loc='best')
ax.set_xlabel('')
ax.set_xticks([], minor=True)
ax.set_title(r'Holt exponential smoothing', fontsize=13.5, fontweight='bold')
sns.despine()
plt.show()
h=24
test=pd.period_range(start=ts.index[-1]+1, periods=h, freq='M')
pred=pd.Series(holt.forecast(h), index=test)
intv1=pd.DataFrame(holt.intervalforecast(h, level=.8), index=test)
intv2=pd.DataFrame(holt.intervalforecast(h, level=.9), index=test)
intv3=pd.DataFrame(holt.intervalforecast(h, level=.99), index=test)
fig, ax = forecast.fanchart(ts['01-2013':], pred, intv1, intv2, intv3)
ax.set_xlabel('')
ax.set_xticks([], minor=True)
plt.title('Seasonally adjusted retail turnover forecast (trend corrected)', fontsize=13.5, fontweight='bold')
sns.despine()
plt.show()
We estimate several candidate specifications below. The AIC suggests a multiplicative model without damping for the original series, but we would need to make a likelihood adjustment or a validation set comparison to selected between this method and an additive model for a log transformation.
ahw=forecast.holtwinters(y, additive=True, damped=False, m=12)
ahw.fit()
ahw.summary()
Additive Holt-winters exponential smoothing Smoothing parameters: alpha (level) 0.575 (0.072) beta (trend) 0.027 (0.017) delta (seasonal) 1.000 (0.173) In-sample fit: MSE 13865.368 RMSE 117.751 Log-likelihood -853.877 AIC 1715.754 BIC 1727.463
mhw=forecast.holtwinters(y, additive=False, damped=False, m=12)
mhw.fit()
mhw.summary()
Multiplicative Holt-winters exponential smoothing Smoothing parameters: alpha (level) 0.459 (0.071) beta (trend) 0.024 (0.012) delta (seasonal) 0.704 (0.144) In-sample fit: MSE 13016.536 RMSE 114.090 Log-likelihood -849.518 AIC 1707.036 BIC 1718.745
mhw_damped=forecast.holtwinters(y, additive=False, damped=True, m=12)
mhw_damped.fit()
mhw_damped.summary()
Multiplicative Holt-winters exponential smoothing (damped trend) Smoothing parameters: alpha (level) 0.330 (0.114) beta (trend) 1.000 (0.558) delta (seasonal) 0.552 (0.170) phi (damping) 0.270 (0.198) In-sample fit: MSE 13608.225 RMSE 116.654 Log-likelihood -852.585 AIC 1715.170 BIC 1729.807
ahw_log=forecast.holtwinters(np.log(y), additive=True, damped=False, m=12)
ahw_log.fit()
ahw_log.summary()
Additive Holt-winters exponential smoothing Smoothing parameters: alpha (level) 0.529 (0.074) beta (trend) 0.022 (0.014) delta (seasonal) 0.823 (0.162) In-sample fit: MSE 0.000 RMSE 0.018 Log-likelihood 358.566 AIC -709.132 BIC -697.423
The smoothed series based on the multiplicative method tracks the original series very closely.
smoothed=pd.Series(mhw.smooth(), index=y.index)
fig, ax= plt.subplots(figsize=(10,6))
y.plot(color='#D62728', label='NSW retail turnover')
smoothed.plot(color='black', label='Smoothed', alpha=0.6)
plt.legend(loc='best')
ax.set_xlabel('')
ax.set_xticks([], minor=True)
ax.set_title(r'Holt-winters exponential smoothing', fontsize=13.5, fontweight='bold')
sns.despine()
plt.show()
We base the diagnostics on the multiplicative Holt-Winters model. The model seems to adequately capture the time series patterns in the data, leading to small and insignicant residual correlations. The residuals do not follow the normal distribution, so that ideally we should use alternative assumptions for computing prediction intervals.
resid = (y-mhw.smooth())[12:] # we remove the first 12 observations as they are for initialisitation only
fig, ax= plt.subplots(figsize=(8,5))
resid.plot(color=blue)
ax.set_xlabel('')
ax.set_xticks([], minor=True)
ax.set_title('Residual plot')
sns.despine()
plt.show()
fig, ax = plt.subplots(figsize=(8,5))
sm.graphics.tsa.plot_acf(resid, lags=40, ax=ax)
sns.despine()
plt.show()
from forecast import histogram, qq_plot
histogram(resid)
plt.show()
qq_plot(resid)
plt.show()
print('Residual skewness: {:.3f}'.format(resid.skew()))
print('Residual kurtosis: {:.3f}'.format(resid.kurt()))
Residual skewness: 0.163 Residual kurtosis: 0.788
The real time forecasting analysis suggests that the log additive Holt Winters model for the log series is the most accurate method for one step ahead forecasts. The Holt-Winters forecasts perform significantly better than a seasonal random walk.
# Real time forecasting
ts=np.log(y) # log series for the log-additive model
validation=y['2012-01':].index
start = y.index.get_loc('2012-01')
pred1 = [] # seasonal random walk
pred2 = [] # additive holt winters
pred3 = [] # mutiplicative holt winters
pred4 = [] # log additive holt winters
actual= []
for i in range(start, len(y)):
actual.append(y.iloc[i]) # actual value
pred1.append(y.iloc[i-12]) # seasonal random walk forecast
model = forecast.holtwinters(y.iloc[:i], additive=True, damped=False, m=12)
model.fit()
pred2.append(model.forecast(1)[0]) # additive holt winters forecast
model = forecast.holtwinters(y.iloc[:i], additive=False, damped=False, m=12)
model.fit()
pred3.append(model.forecast(1)[0]) # multiplicate holt winters forecast
model = forecast.holtwinters(ts.iloc[:i], additive=True, damped=False, m=12)
model.fit()
resid = (ts.iloc[:i]-model.smooth())[12:] # residuals
y_pred = np.exp(model.forecast(1)[0])*np.mean(np.exp(resid)) # forecast with a retransformation adjustment
pred4.append(y_pred) # long additive holt winters forecast
columns=['Seasonal RW', 'Additive', 'Multiplicative', 'Log additive', 'Actual']
results = np.vstack([pred1,pred2,pred3,pred4,actual]).T
results = pd.DataFrame(results, columns=columns, index=validation)
from statlearning import rmse_jack
table = pd.DataFrame(0.0, index=results.columns[:-1], columns=['RMSE','SE'])
for i in range(4):
table.iloc[i,0], table.iloc[i,1] = rmse_jack(results.iloc[:,i], results.iloc[:,-1])
table.round(2)
RMSE | SE | |
---|---|---|
Seasonal RW | 385.73 | 21.12 |
Additive | 146.49 | 17.78 |
Multiplicative | 134.39 | 10.07 |
Log additive | 129.95 | 10.75 |
Finally, we generate a two year forecast based on the selected model.
h=24
test=pd.period_range(start=y.index[-1]+1, periods=h, freq='M')
pred=pd.Series(ahw_log.forecast(h), index=test)
var = ahw_log.forecastvariance(h)
pred=np.exp(pred+var/2) # retransformation under the assumption of normality, for simplicity
intv1=pd.DataFrame(ahw_log.intervalforecast(h, level=.8), index=test)
intv2=pd.DataFrame(ahw_log.intervalforecast(h, level=.9), index=test)
intv3=pd.DataFrame(ahw_log.intervalforecast(h, level=.99), index=test)
intv1=np.exp(intv1)
intv2=np.exp(intv2)
intv3=np.exp(intv3)
fig, ax = forecast.fanchart(y['01-2014':], pred, intv1, intv2, intv3)
ax.set_xlabel('')
ax.set_xticks([], minor=True)
plt.title('NSW retail turnover forecast (log additive Holt-winters model)', fontsize=13.5, fontweight='bold')
sns.despine()
plt.show()