Statistical models: from AR to SARIMAX, smoothing and Theta

1. The catalogue of classical methods

The deck opens with a list, and the list is worth reading as a family tree rather than an inventory. Every name in it is a combination of four ideas: regress on past values, regress on past errors, difference to remove trend, and repeat the whole thing at a seasonal lag.

2. Autoregressive models AR(p)

In autoregressive models of order p, the future value of a variable is assumed to be a linear combination of the last p observations, added to a random component and a constant:

yₜ = c + φ₁yₜ₋₁ + φ₂yₜ₋₂ + … + φₚyₜ₋ₚ + εₜ

Two restrictions are stated immediately and must be quoted in an exam: AR is feasible only for stationary processes, and it only models a linear dependency. Everything in chapter 3 exists so that this sentence can be satisfied. The slides repeat the point later: the method is suitable for univariate time series without trend and seasonal components.

Generating an AR(2)

The code that generates a synthetic AR(2) also documents its stability region:

np.random.seed(995)
n = 200
# for stationarity: phi2 + phi1 < 1, phi2 - phi1 < 1, |phi2| < 1
phi1, phi2 = 0.6, -0.3
errors = np.random.normal(0, 1, n)
y = np.zeros(n)
for t in range(2, n):
    y[t] = phi1 * y[t-1] + phi2 * y[t-2] + errors[t]
plt.plot(y); plt.title("AR(2) series")
plt.xlabel("time"); plt.ylabel("values"); plt.show()

AR(2) on the fil rouge

The deck fits an AR(2) on log-diff data to ensure stationarity, following a least-mean-squares analysis on data lagged up to p periods — that is, an ordinary regression of the log-differenced series on its own two previous values. The estimated model is:

Both coefficients are negative, which is exactly what an alternating quarterly series should produce: a high quarter predicts a low one. The fitted values track the log-diff column reasonably (0.16, 0.47, -0.69, -0.09, 0.66, …) but on twenty points nothing is being estimated with confidence — a point the slides make repeatedly.

from statsmodels.tsa.ar_model import AutoReg

data = [x + np.random.random() for x in range(1, 100)]   # artificial dataset
model = AutoReg(data, lags=1)
model_fit = model.fit()
pred = model_fit.predict(0, len(data))    # in-sample, both ends included
fore = model_fit.forecast(4)              # out of sample only
print(f"[AR]: {fore}")

A full AR(2) forecast, transforms included

The complete pipeline on the airline series is the template every later model reuses: preprocess, fit, predict, then undo every transform in reverse order.

3. How good is the model? AIC and AICc

AR gives us a first model of a time series. But how good is the model? The Akaike Information Criterion is an estimator of the prediction error, i.e. of the relative quality, of statistical models for a given set of data. The interpretation to remember: it estimates how much information is lost by using the model instead of the actual data.

AIC = -2 log(L) + 2k

where k is the number of estimated parameters in the model and L is the maximum likelihood of the model. Given a set of candidate models for the data, the preferred model is the one with the minimum AIC value.

The likelihood expression can be complicated, but if the errors are normally distributed (residuals Gaussian) and the sample size is large enough, the log-likelihood can be approximated by

log(L) ≈ -(n/2) · [ log(2π) + log(σ²) + 1 ]

with σ² the variance of the residuals. Two terms, two forces: the fit term rewards small residual variance, the 2k term punishes parameters. AIC is a bargain between them.

The worked comparison

Suppose there are three candidate AR(2) models. For AR(2), k = 3 (the parameters c, φ₁, φ₂). With n = 40 samples and residual variances 10, 12 and 15, the AIC values come out as 109, 112 and 117. The quantity

exp( (AICᵐᵢₙ - AICᵢ) / 2 )

is proportional to the probability that the i-th model is the one that minimises the information loss. The slides read off the conclusion: the second model is 0.161 times as probable as the first (the minimum) to minimise information loss, and the third is 0.017 times as probable.

AICc, for small samples

When the sample size is small, AIC will select models that have too many parameters — AIC will prefer overfits. AICc is AIC with a correction for small sample sizes. If the model is univariate, linear in its parameters, and has normally distributed residuals:

AICc = AIC + 2k(k+1) / (n - k - 1)

where n is the sample size and k the number of parameters. It is AIC with an extra penalty term for the number of parameters, and the penalty vanishes as n grows.

4. Forecast accuracy metrics

Given actual data Y = (y₁, …, y_n) and forecast data F = (f₁, …, f_n), five error measures are defined. They differ in one respect only — how they treat the sign and the scale of the error — and that difference is the whole reason to know all five.

import numpy as np
from statsmodels.tsa.stattools import acf

def forecast_accuracy(forecast, actual):
    mape = np.mean(np.abs(forecast - actual)/np.abs(actual))   # MAPE
    me   = np.mean(forecast - actual)                          # ME
    mae  = np.mean(np.abs(forecast - actual))                  # MAE
    mpe  = np.mean((forecast - actual)/actual)                 # MPE
    rmse = np.mean((forecast - actual)**2)**.5                 # RMSE
    corr = np.corrcoef(forecast, actual)[0,1]                  # correlation coeff
    acf1 = acf(forecast-actual)[1]                             # ACF1
    return({'mape':mape, 'me':me, 'mae': mae, 'mpe': mpe, 'rmse':rmse,
            'acf1':acf1, 'corr':corr})

Two entries in that dictionary are worth pausing on. corr is the correlation between forecast and actual — a model can be perfectly correlated and still systematically wrong by a constant. acf1 is the lag-1 autocorrelation of the errors: if it is significantly non-zero, the residuals still contain structure the model failed to capture, which is the diagnostic criterion of section 9.

A worked evaluation: SARIMA on 2008

The slides evaluate a SARIMA model on the fil rouge, training on 2004-2007 and testing on the four quarters of 2008.

5. LASSO: selection and regularization

LASSO (L1 regularization) stands for Least Absolute Shrinkage and Selection Operator. It performs selection and regularization — constraining the values — of the coefficients of the linear model, and has been extended to other statistical models.

Given a linear regression model y = β₀ + β₁x, LASSO requires solving

min Σᵢ ( yᵢ - (β₀ + β₁xᵢ) )²  +  λ Σⱼ |βⱼ|

vector form:   min ‖ y - Xβ ‖²  +  λ ‖β‖₁      (λ a parameter)

The crucial consequence: LASSO constrains the absolute values of the coefficients. They can be zero, removing the corresponding variable. That is what makes it a selection operator and not merely a shrinkage one — an L2 penalty would shrink coefficients towards zero without ever reaching it.

It is mostly used for feature selection in multivariate forecasting, but can be adapted to the univariate case. And here the slides state the manoeuvre that unlocks the whole of chapters 5 and 6: to use LASSO for univariate forecast, one must transform the time series into a supervised learning problem by creating a matrix of lagged features.

# create lagged features for an AR(12)-style design matrix
p = 12                                              # AR order (number of lags)
X = np.column_stack([y[i:-p+i] for i in range(p)])  # y_{t-1} ... y_{t-p}
y_target = y[p:]                                    # target variable y_t

X_train, X_test, y_train, y_test = train_test_split(
    X, y_target, test_size=0.2, random_state=42)
lasso = Lasso(alpha=0.1)        # alpha is the regularization strength
lasso.fit(X_train, y_train)
print("Coefficients:", lasso.coef_)
predictions = lasso.predict(X_test)

6. Moving average models MA(q)

In moving average models of order q, the future value of a variable is assumed to be equal to the average of the observations, added to a linear combination of the last q errors:

yₜ = μ + θ₁εₜ₋₁ + θ₂εₜ₋₂ + … + θₔεₜ₋ₔ + εₜ

with μ the average of the series and θ_j the model parameters. Random errors are usually assumed to have a normal distribution with mean 0 and variance σ². As with AR, the method is suitable for univariate time series without trend and seasonal components.

The slides summarise MA in one sentence that clears up most confusion: MA is just a regression on previous periods errors. Not on values — on the model's own past mistakes. Where AR says "yesterday was high, so today will be high", MA says "yesterday I under-predicted, so today I should aim higher".

MA(1) on the fil rouge

Fitted on the log-diff series, the MA(1) coefficients are μ = -0.06 and θ₁ = -0.21, giving fitted values -0.11, 0.12, -0.23, -0.17, 0.09, -0.02, … The slides append a diagnostic remark that matters more than the numbers: there is a hint of seasonality in the error. Structure left in the residuals means the specification is incomplete — which is precisely the argument for going seasonal in section 9.

from statsmodels.tsa.arima.model import ARIMA
from random import random

data = [x + random() for x in range(1, 100)]   # artificial dataset
model = ARIMA(data, order=(0, 0, 1))           # MA(1): p = 0, d = 0, q = 1
model_fit = model.fit()
yhat = model_fit.predict(len(data), len(data))
print(yhat)

7. ARMA(p,q)

Autoregressive and moving average models can be combined in a more general class, the ARMA models:

yₜ = c + Σᵢ₊₁ᵖ φᵢyₜ₋ᵢ  +  Σⱼ₊₁ₔ θⱼεₜ₋ⱼ  +  εₜ

ARMA models the next step in the sequence as a linear function of the observations and residual errors at preceding time steps. It remains suitable for univariate time series without trend and seasonal components — combining the two mechanisms does not lift the stationarity requirement.

n = 200
phi   = 0.6    # AR coefficient, for stationarity |phi| < 1
theta = 0.4    # MA coefficient
errors = np.random.normal(0, 1, n)
y = np.zeros(n)
for t in range(1, n):
    y[t] = phi * y[t-1] + theta * errors[t-1] + errors[t]

On the fil rouge, an ARMA(2,1) is fitted by downscaling ARIMA to order=(2, 0, 1). The verdict from the slides is candid: MA improves a little over AR. More data would be needed.

8. ARIMA(p,d,q) and the Box-Jenkins methodology

ARMA models can be used only on stationary processes, without trend or seasonality. ARIMA models are ARMA models that work on a diff time series — differencing as preprocessing folded into the model — and thus lift the trend application limit.

ARIMA(p,d,q):   ∇ᵈyₜ = c + Σᵢ₊₁ᵖ φᵢ∇ᵈyₜ₋ᵢ + Σⱼ₊₁ₔ θⱼεₜ₋ⱼ + εₜ

p   order of the autoregressive component
q   order of the moving average component
d   diff degree (usually d = 1; if d = 0, ARIMA = ARMA)

Or, as the slides put it in words:

Predicted yₜ = Constant
              + linear combination of lags of y      (up to p lags)
              + linear combination of lagged forecast errors (up to q lags)

Configuring it: the Box-Jenkins methodology

The classical approach for fitting the parameters of an ARIMA model is the Box-Jenkins methodology, first described in George Box and Gwilym Jenkins, Time Series Analysis: Forecasting and Control, 1970. Three steps:

1. Model identification. Identify trends, seasonality and autoregression elements to get an idea of the amount of differencing and the size of the lag that will be required.
2. Parameter estimation. Use a fitting procedure to find the coefficients of the regression model.
3. Model checking. Use plots and statistical tests of the residual errors to determine the amount and type of temporal structure not captured by the model.

The process is repeated until a desirable level of fit is achieved on the in-sample or out-of-sample observations. Note that this is precisely the specification / fitting / diagnosis triple of chapter 2, wearing different clothes.

Grid search and Hammersley sampling

Algorithms are very sensitive to their control parameters (aka hyperparameters), and optimizing parameter setting is an active area of research, with many methods from simple to very complicated. The simplest one is grid search: define an interval for each parameter and test all combinations, for integer parameters.

# ARIMA grid search: p in [0,3], d in [0,2], q in [0,2]
for (p = 0; p <= 3; p++)
   for (d = 0; d <= 2; d++)
      for (q = 0; q <= 2; q++)
      { performance = ARIMA(p,d,q);
        if (performance > best_performance) ... etc

For continuous parameters it is possible to randomly generate values inside the interval, but independent generation along each dimension could lead to uneven sampling in multidimensional spaces. A possible solution is Hammersley sampling — the topic returns in chapter 9 with its low-discrepancy siblings Sobol and Halton.

ARIMA in practice, with the reconstruction

The slides give two utility versions of the same pipeline, one working on pandas Series and one on numpy arrays, and both make the same point: preprocessing must be undone in reverse order after forecasting.

import pandas as pd, numpy as np, matplotlib.pyplot as plt
from statsmodels.tsa.arima.model import ARIMA

df = pd.read_csv("M3C_monthly.csv")
rawdata    = df.iloc[490, 6:].values.astype(float)
train, test = rawdata[:-12], rawdata[-12:]
logdata    = pd.Series(np.log(train))          # log transform

model     = ARIMA(logdata, order=(1,1,1))      # let ARIMA do the differencing
model_fit = model.fit()

# with d=1 predictions come back in the input scale (here the log scale).
# in-sample predictions start from 1, to avoid the NaN from differencing
predlog = model_fit.predict(start=1, end=len(logdata)-1)
forelog = model_fit.forecast(12)               # out of sample, also in log scale

recpred = np.exp(predlog)                      # inverse transform
recfore = np.exp(forelog)

The array version spells the inversion out by hand, and it is the version worth studying because it shows what the library is doing for you:

logdata = np.log(train)
logdiff = np.array([logdata[i]-logdata[i-1] for i in range(1, len(logdata))])

fore = np.zeros(len(test))                     # forecast, you can do better

# Postprocessing, reconstruction
reclogdiff = np.insert(logdiff, 0, logdata[0]).cumsum()   # undo diff  = cumsum
recfore    = np.insert(fore, 0, reclogdiff[-1]).cumsum()
recpred    = np.exp(reclogdiff)                           # undo log   = exp
recfore    = np.exp(recfore)

9. SARIMA, SARIMAX and ARCH-GARCH

SARIMA is an extension, proposed by Box and Jenkins (1970), of ARIMA models, which can be applied also to data with seasonality. A seasonal diff is included, and the slides warn that model identification can be complex.

SARIMA(p, d, q)(P, D, Q)ₘ

Multiple seasonalities are addressed by an extension at the end of the deck — the BATS/TBATS and MSTL family already met in chapter 2.

Two ways to fit one

import pmdarima as pm            # pip install pmdarima
df = pd.read_csv('FilRouge.csv', names=['sales'], header=0)
ds = df.sales
model = pm.auto_arima(ds.values, start_p=1, start_q=1,
                      test='adf', max_p=3, max_q=3, m=4,
                      start_P=0, seasonal=True,
                      d=None, D=1, trace=True,
                      error_action='ignore',
                      suppress_warnings=True,
                      stepwise=True)      # stepwise=False for a full grid
morder     = model.order                  # p,d,q
mseasorder = model.seasonal_order         # P,D,Q,m
fitted = model.fit(ds)
yfore  = fitted.predict(n_periods=4)      # forecast
ypred  = fitted.predict_in_sample()

from statsmodels.tsa.statespace.sarimax import SARIMAX

sarima_model = SARIMAX(ds, order=(0,2,2), seasonal_order=(0,1,0,4))
sfit = sarima_model.fit()
sfit.plot_diagnostics(figsize=(10, 6)); plt.show()

ypred = sfit.predict(start=0, end=len(df))       # in-sample

forewrap     = sfit.get_forecast(steps=4)        # out-of-sample
forecast_ci  = forewrap.conf_int()
forecast_val = forewrap.predicted_mean

SARIMAX and exogenous variables

SARIMAX is an extension of SARIMA that also includes the modeling of exogenous variables. Exogenous variables — also called covariates, predictors or externals — are parallel input sequences that have observations at the same time steps as the original series. The primary series is referred to as endogenous data, to contrast it with the exogenous sequences.

The distinction that matters mechanically: the observations for exogenous variables are included in the model directly at each time step, and are not modeled in the same way as the primary endogenous sequence — they are not given an AR or MA process of their own. Consequence: to forecast h steps ahead you must supply the exogenous values for those h steps.

from statsmodels.tsa.statespace.sarimax import SARIMAX

m1 = SARIMAX(train_y,
             exog = train_x,
             order = (1,1,1),
             seasonal_order = (1,1,1,12),
             enforce_stationarity  = False,
             enforce_invertibility = False).fit(disp=False)
f1 = m1.forecast(12, exog=test_x)      # future exogenous values REQUIRED
print(m1.summary())

SARIMAX can also be used to model the subsumed models, adding exogenous variables: ARX, MAX, ARMAX and ARIMAX. It is suitable for univariate time series with trend and/or seasonal components and exogenous variables: it predicts a unique time series, using externals.

ARCH and GARCH

Autoregressive Conditional Heteroskedasticity (ARCH) is used to analyze volatility in time series in order to forecast future volatility. In the financial world, ARCH modeling is used to estimate risk by providing a model of volatility that more closely resembles real markets.

Var(yₜ | yₜ₋₁, …, yₜ₋ₘ) = σ²ₜ = α₀ + α₁y²ₜ₋₁ + … + αₘy²ₜ₋ₘ

The idea is a neat inversion of everything above: instead of an autoregression on the level, an autoregression on the variance. Large squared shocks predict large future variance — volatility clustering. The work earned R. Engle and C. Granger the Nobel prize in economics in 2003.

10. Prediction intervals

The topic is complex, and the slides give only a hint — but the framing is the same one chapter 2 opened with. Forecasts are usually wrong. We provide a specific value, meaning that the future value will be similar to the suggested one. Actually, we should provide a value distribution, which can be described by an interval on the expected values. A hypothesis on the distribution is needed — usually Student's t, here normality — and the forecast value has mean μ with prediction interval [μ - δ, μ + δ].

When the API does not give you one

A three-step recipe, which the slides note will be plausible also for neural models:

1. Associate a random variable with each value to forecast (time periods t, t+1, …).
2. Apply the model on sliding windows over historic data, computing the variance for each lag. The lag-1 forecast has its own error variance, the lag-2 forecast a larger one, and so on.
3. Apply the model for the forecast, assuming as prediction interval the forecast value ± the standard deviation computed for the corresponding lag.

This is why the interval fans out: the variance is measured separately per horizon, and empirically it grows. It is the mechanical explanation of the fourth forecast attribute in chapter 2 — reliability decreases with distance in the future.

11. Exponential smoothing: SES, Holt-Winters, Theta

Simple Exponential Smoothing

SES models the next time step as an exponentially weighted linear function of observations at prior time steps. It requires a single parameter, α, called the smoothing factor or smoothing coefficient. The equation is disarmingly simple:

fₜ₊₁ = αyₜ + (1 - α)fₜ

which, projected back in time, corresponds to

fₜ₊₁ = αyₜ + α(1-α)yₜ₋₁ + α(1-α)²yₜ₋₂ + …

The parameter controls the rate at which the influence of observations at prior time steps decays exponentially: values close to 1 mean the model pays attention mainly to the most recent observations, values close to 0 mean more of the history is taken into account. As with AR and MA, SES is suitable for univariate time series without trend and seasonal components — and the fil rouge demonstration in the slides is annotated exactly with what goes wrong: lost the trend, smoothed contributions.

def simple_exp_smooth(data, nforecasts=1, alpha=0.4):
   n = len(data)
   f = np.full(n + nforecasts, np.nan)              # forecast array
   data = np.append(data, [np.nan] * nforecasts)    # forecast placeholders
   f[1] = data[0]                                   # initialise first forecast
   for t in range(2, n+1):                          # predictions
      f[t] = alpha * data[t - 1] + (1 - alpha) * f[t - 1]
   for t in range(n+1, n+nforecasts):               # forecast
      f[t] = alpha * f[t - 1] + (1 - alpha) * f[t - 2]
   return pd.DataFrame.from_dict({"Data": data, "Forecast": f, "Error": data - f})

Holt-Winters Exponential Smoothing (HWES)

Also called Triple Exponential Smoothing — there would be a double in between — HWES models the next time step as an exponentially weighted linear function of observations at prior time steps, taking trends and seasonality into account. Three parameters, one per component:

Trend and seasonality may be modeled as either additive or multiplicative, for a linear or an exponential change respectively — the same choice chapter 2 made when decomposing. The method is suitable for univariate time series with trend and/or seasonal components: unlike everything else in this chapter so far, it needs no differencing.

from statsmodels.tsa.holtwinters import ExponentialSmoothing

model = ExponentialSmoothing(train, seasonal_periods=4, trend="add",
                             seasonal     = "mul",
                             damped_trend = True,
                             use_boxcox   = True,
                             initialization_method = "estimated")
hwfit = model.fit()
yfore = hwfit.predict(len(train), len(train)+3)
print(yfore)

The deck also gives a nuts-and-bolts implementation, and it repays reading because the initialisation reveals the structure: the level a[0] starts as the mean of the first season, the trend b[0] as the average difference between the second season and the first divided by m, and the seasonal factors s[i] as each of the first m observations divided by that initial level.

# x: series, m: seasonality, nfor: num forecast
def holtwinters(x, m, nfor, alpha=0.5, beta=0.5, gamma=0.5):
   Y = x[:]
   a = [sum(Y[0:m]) / float(m)]                          # level
   b = [(sum(Y[m:2 * m]) - sum(Y[0:m])) / m ** 2]        # trend
   s = [Y[i] / a[0] for i in range(m)]                   # seasonal (multiplicative)
   y = [(a[0] + b[0]) * s[0]]

   for i in range(len(Y) + nfor):
      if i == len(Y):
         Y.append((a[-1] + b[-1]) * s[-m])               # roll the forecast in
      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])
   return Y, rmse

Three lines, three components: each of a, b and s is updated as a blend between what the newest observation implies and what the previous estimate said — the same convex combination as SES, applied three times.

Theta

The Theta model of Assimakopoulos and Nikolopoulos (2000) is a simple method for forecasting that involves three moves:

1. Fitting two lines. The original time series (Theta = 1): the raw data without modification, retaining the trend and seasonal patterns present in the data. And the damped trend line (Theta = 2): it emphasises the long-term trend, derived by applying a double exponential smoothing (Holt's method with no dampening) to estimate a more stable trend.
2. Forecasting the lines using a Simple Exponential Smoother.
3. Combining the forecasts from the two lines to produce the final forecast.

A seasonality test examines the ACF at the seasonal lag m; if that lag is significantly different from zero, the data is deseasoned with either a multiplicative or an additive method — and the forecasts are reseasoned at the end if needed. The parameters are θ and α, the first estimated from an OLS regression and the second the SES smoothing parameter. The slides note the limiting behaviour: ultimately θ only plays a role in determining how much the trend is damped; if θ is very large, the forecast of the model is identical to that from an Integrated Moving Average with a drift.

from statsmodels.tsa.forecasting.theta import ThetaModel
from statsmodels.datasets import get_rdataset

airpass = get_rdataset('AirPassengers').data.value.values
theta_model = ThetaModel(airpass, period=12)
fit = theta_model.fit()
print(fit.summary())
forecast = fit.forecast(steps=12)      # 12 months ahead

Check your understanding