Time series, decomposition and the empirical STL

1. Predictive analytics: forecast now, decide now

The deck opens with a compact definition: predictive analytics is an area of statistics that deals with extracting information from data and using it to predict trends and behavior patterns. And then the mechanism: its core relies on capturing relationships between explanatory variables and the predicted variables from past occurrences, and exploiting them to predict the unknown outcome.

The slides then split the name in two, and the split is the whole justification of the discipline:

A forecast nobody acts on is a hobby. The management examples given are all decisions with a lead time attached: forecast of goods or services demand; forecast of workforce availability or need; forecast of material requirements and warehouse availabilities; forecast of revenues, profits and losses in order to define future investments.

2. Horizons, and what a forecast actually is

Different decisions need different distances into the future, and the slides organise them into three bands. The band you are in determines both the model and the tolerance for error.

Four attributes of any forecast

The slides state four properties that a student is expected to repeat without hesitation, starting with the least flattering one.

3. A taxonomy of prediction models

The deck lays out the space of forecasting methods as a tree with two big branches — objective and subjective — and it is worth keeping the subjective branch in mind, because the rest of the course lives entirely on the objective one.

Subjective models

Very varied. The slides list: sales force composite (aggregation of estimates made by sales agents); polls among customers; opinions of managers; and Delphi — ask the question, gather individual opinions, share the opinions anonymously, re-ask the question, repeat the procedure until consensus. The anonymity is the mechanism: it is what stops the loudest voice from anchoring the room.

Objective models and their three phases

Objective models are based on mathematical models, and the slides break the workflow into three high-level phases which recur in every later chapter, under different names.

4. Time series: definition, examples, plots

A time series is simply a set of time-indexed values, and it represents the dynamics of a process over time. In the model, we assume to have n observations coming from as many dependent random variables. And the objective of studying them is stated with disarming honesty: to find patterns in past data which can be assumed to repeat in the future. The assumption is doing all the work; when it breaks, so does the forecast.

A time-series plot (timeplot, time series graph) is a bidimensional chart of the series: the vertical axis is the variable of interest, the horizontal axis the time units.

Series are everywhere

The three series this course keeps coming back to

df = pd.read_csv("BoxJenkins.csv")
plt.plot(df.Passengers, label="passengers")
plt.legend()

df = pd.read_csv("FilRouge.csv")
plt.plot(df.sales, label="TV sales")
plt.grid(); plt.legend(); plt.show()

The fil rouge itself, quarter by quarter, is worth having in front of you for the rest of the chapter:

5. Univariate, multivariate, exogenous — statistics or machine learning

Forecasting a time series means predicting the values it is going to take in the future. The slides distinguish three settings, and the third is the one most real problems fall into.

Two approaches, and an honest comparison

The slides put the two families side by side without pretending the choice is obvious.

6. Regression, prediction, forecast: residuals and forecast errors

Regression analysis tries to estimate the relationships among variables, usually between a dependent variable and one or more independent variables (or predictors). It estimates a function. It is not limited to prediction and forecasting: more generally, it wants to determine the value of the dependent variable when the independent variables are fixed.

Then comes a vocabulary distinction the slides are emphatic about, because ML literature blurs it. Regression is the general term often used in the machine learning literature to refer to prediction tasks. But:

7. Decomposition: trend, season, cycle, noise, accident

Time series decomposition deconstructs a time series into several components, each representing one of the underlying categories of patterns. The common decomposition is into trend, seasonal, cyclical and random components — plus, in the slides, an occasional component that deserves its own discussion.

The camper sales example from the slides shows all three of the structured components at once: a seasonality where Q3 > Q2 > Q1 > Q4 of the previous year, a long-term trend, and economic cycles visible as the recessions of the 1980s and 1990s. The GDP of Italy from 1861 to 2017, and a power-energy forecast, are shown as cases where an occasional component forces a model change: no amount of parameter fitting rescues a model across a structural break.

Additive, multiplicative, mixed

The hypothesis behind decomposition is that the observed values result from the combination of different components. How they combine is a modelling choice with a clear rule of thumb:

STL, and multiple seasonalities

A common model is STL — seasonal decomposition of time series by Loess, where Loess is short for local regression. It decomposes the series into seasonal, trend and irregular, with the cyclical component included in the trend plot.

from statsmodels.tsa.seasonal import seasonal_decompose
plt.rcParams['figure.figsize'] = (10.0, 6.0)
df = pd.read_csv('BoxJenkins.csv', usecols=["Passengers"])
ds = df[df.columns[0]]                  # converts to series
result = seasonal_decompose(ds, model='multiplicative', period=12)
result.plot()
plt.show()

Real series often have more than one period at once — a daily and a weekly rhythm, say. The slides list three packages that can identify multiple seasonalities: MSTL (in statsmodels), Facebook's Prophet, and BATS / TBATS.

8. The empirical STL, step by step on the fil rouge

Rather than call a library and admire the plot, the deck rebuilds STL by hand on the fil rouge, and projects it into the future. Three ingredients are defined first.

The trend by least squares

We look for the interpolation line y = a·x + b minimising the mean square error over all data pairs (x_i, y_i). Analytically, with mean values x̄ and ȳ:

a = Σ(xᵢ - x̄)(yᵢ - ȳ) / Σ(xᵢ - x̄)²
b = ȳ - a · x̄

which in the vectorised form the slides use for Python becomes:

On the twenty fil rouge points the arithmetic gives x̄ = 10.5, ȳ = 395.25, and

Sanity check on the first row: 25.888 × 1 + 123.426 = 149.314, against an observed 128. The gap is what the seasonal coefficient will absorb.

9. Correlograms and finding the period

Seasonality induces periodic fluctuations around the trend, and to model it we need a coefficient to multiply the trend by. But first we must find the periodicity (wavelength) — the number of seasons. The slides pick autocorrelograms, noting that the correlation of two time series can be quantified by Pearson's index (parametric) or Spearman's (nonparametric).

A correlogram, or autocorrelation plot, is a chart that represents the autocorrelation of a time series as a function of the lag used for computing the autocorrelation. Concretely: build a table where beside the column y_t you place the columns y_t-1, y_t-2, …, y_t-K, start from the (K+1)-th row so that all the compared series have equal length, and compute for each k the correlation

rₖ = Σₕ(yₜ - ȳ)(yₜ₋ₖ - ȳ) / Σₕ(yₜ - ȳ)²

Pairs (k, rₖ) are then plotted with lags on the x-axis and the corresponding correlation on the y-axis.

The fil rouge correlogram

Computed on the raw fil rouge with Pearson, the slides report:

The reading is immediate and is the one the slides write on the chart: peaks at lag 4 and lag 8, therefore seasonality of 4. Quarterly data, four seasons, as the calendar column of the original table promised.

Three typical shapes

Correlograms can be widely different, but the slides single out three cases you should be able to recognise on sight: the correlogram of the original data series; the correlogram of the residuals after detrending; and the correlogram of the residuals after detrending and deseasoning. A slowly decaying correlogram means a trend is still present; a comb of peaks at multiples of L means seasonality is still present; a flat correlogram near zero means there is nothing structured left to model — which is the goal.

import numpy as np, pandas as pd, matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

df = pd.read_csv('FilRouge.csv', usecols=[1], names=['value'], header=0)
plt.rcParams.update({'figure.figsize':(9,7), 'figure.dpi':120})
fig, axes = plt.subplots(2, 2, sharex=True)
axes[0, 0].plot(df.value); axes[0, 0].set_title('Original Series')
plot_acf(df.value, ax=axes[0, 1])
# 1st Differencing
axes[1, 0].plot(df.value.diff()); axes[1, 0].set_title('1st Order Differencing')
plot_acf(df.value.diff().dropna(), ax=axes[1, 1])
plt.show()

10. Seasonal coefficients, deseasoning and the projected forecast

With the period known, the seasonal coefficient is built in three moves.

1. Relative elongation. For each point, divide the observation by the trend value: y_t / T_t. For t = 1: 128 / 149.3143 = 0.857252.
2. Average by season index. Group the elongations by t mod 4 and average each group. The slides do it in Excel with a SUMPRODUCT over rows matching the modulus, divided by the number of complete seasons.
3. Look the coefficient up per period (a VLOOKUP on the modulus) and multiply the trend by it.

The four coefficients that come out of the fil rouge are:

Multiplying trend by coefficient reconstructs the seasonal model: at t = 1, 149.314 × 1.29234 = 192.96; at t = 5, 252.866 × 1.29234 = 326.79; at t = 20, 641.186 × 1.142187 = 732.35.

Deseasoning

To eliminate the effect of s_t we divide by s_t; what is left is trend plus random component. The slides tabulate it:

The projected forecast

Having the coefficients of the trend and of the seasonal components, we can project the model onto future periods. Extend the trend line past t = 20 and multiply by the coefficient of the corresponding season:

11. Sampling: train, test, and why k-fold does not apply

The slides are categorical: in any learning process, performance should be estimated with out-of-sample validation, withholding some of the sample data from model identification (the training phase), then using the model to make forecasts for the hold-out data (the test data). The testing dataset should never be used in the learning step.

Two definitions follow, and they close the loop with section 6: predictions on train values are called fitted values, and the corresponding errors are called residuals.

Overfitting

Errors on train data can be much smaller than those on test data, in case of overfitting. The slides name two circumstances in which it is likely:

1. a model with a large number of parameters has been fitted to a too small sample of data;
2. the model has been selected from a large set of potential models precisely by minimising the RMSE in the estimation period.

The second is subtler and more common: if you grid-search a hundred configurations and keep the one with the lowest in-sample error, the winning score is a selection artefact, not a performance estimate. Chapter 9 returns to this when it discusses grid search and Optuna.

Making the split

The mechanics are mundane: select an arbitrary split point in the ordered list of observations and create two datasets. Depending on how much data you have and how much the model needs, splits of 50-50, 70-30 or 90-10 are used.

ds = pd.read_csv('BoxJenkins.csv', header=0)
X = ds.Passengers.values
train_size = int(len(X) * 0.66)
train, test = X[0:train_size], X[train_size:len(X)]
print('Observations: %d' % (len(X)))
print('Training Observations: %d' % (len(train)))
print('Testing Observations: %d' % (len(test)))
plt.plot(train)
plt.plot([None for i in train] + [x for x in test])
plt.show()

Check your understanding