Chapter 47: Classical Forecasting: ARIMA & ETS

What Makes a Time Series Tick

A time series is a sequence of observations indexed in time order: $y_1, y_2, \ldots, y_T$. Unlike cross-sectional data, adjacent observations are statistically dependent — the value today carries information about the value tomorrow. Classical forecasting methods exploit this autocorrelation structure explicitly.

Before fitting any model, you must understand the structure of your series along four axes:

• Trend: a persistent long-run direction (upward, downward, or flat).
• Seasonality: a repeating pattern tied to the calendar (weekly, monthly, annual).
• Cyclicality: longer, irregular fluctuations not tied to fixed calendar periods.
• Irregular / noise: the unpredictable remainder.

Additive vs. Multiplicative Decomposition

An additive decomposition assumes the components sum:

A multiplicative decomposition assumes they multiply:

Multiplicative decomposition is appropriate when the amplitude of seasonal swings grows with the level of the series — a telltale sign is a cone-shaped seasonal pattern on a raw plot. Taking logs converts multiplicative to additive, which is why log-transforming is a ubiquitous first step.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import seasonal_decompose

# Load AirPassengers (classic benchmark series)
url = "https://raw.githubusercontent.com/jbrownlee/Datasets/master/airline-passengers.csv"
df = pd.read_csv(url, parse_dates=["Month"], index_col="Month")
df.index.freq = "MS"  # monthly start

# Multiplicative decomposition (seasonal amplitude grows with level)
result = seasonal_decompose(df["Passengers"], model="multiplicative", period=12)
result.plot()
plt.tight_layout()
plt.show()

Stationarity

Most classical models — in particular ARIMA — require weak stationarity: a constant mean, constant variance, and autocovariances that depend only on the lag, not on absolute time. A series with a trend or growing variance is non-stationary.

Two formal tests are standard:

• Augmented Dickey-Fuller (ADF): the null hypothesis is a unit root (non-stationary). Rejecting (low p-value) implies stationarity.
• KPSS: the null hypothesis is stationarity. Rejecting implies non-stationarity.

Using both together is more reliable because they have complementary power.

from statsmodels.tsa.stattools import adfuller, kpss

log_series = np.log(df["Passengers"])

adf_stat, adf_p, *_ = adfuller(log_series, autolag="AIC")
print(f"ADF  stat={adf_stat:.3f}, p={adf_p:.4f}")

kpss_stat, kpss_p, *_ = kpss(log_series, regression="c", nlags="auto")
print(f"KPSS stat={kpss_stat:.3f}, p={kpss_p:.4f}")
# ADF fails to reject -> non-stationary; KPSS rejects -> confirms non-stationary
# Solution: first-difference the log series

The standard remedies for non-stationarity are differencing (subtract consecutive values) and log transformation. The order of differencing $d$ is a key parameter of ARIMA.

The ACF and PACF: Your Diagnostic Compass

The autocorrelation function (ACF) at lag $k$ measures the linear correlation between $y_t$ and $y_{t-k}$, accounting for all intermediate lags:

The partial autocorrelation function (PACF) measures the correlation between $y_t$ and $y_{t-k}$ after removing the effect of $y_{t-1}, \ldots, y_{t-k+1}$. PACF isolates the direct linear contribution of each lag.

The classical pattern-matching rules for model identification:

• AR(p): PACF cuts off after lag <!--MATHBLOCK11-->; ACF decays geometrically.
• MA(q): ACF cuts off after lag <!--MATHBLOCK12-->; PACF decays geometrically.
• ARMA(p,q): both tail off gradually.

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

diff_log = np.log(df["Passengers"]).diff().dropna()

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
plot_acf(diff_log, lags=40, ax=axes[0], title="ACF — First-Differenced Log")
plot_pacf(diff_log, lags=40, ax=axes[1], title="PACF — First-Differenced Log")
plt.tight_layout()
plt.show()

Spikes at seasonal lags (12, 24, 36 for monthly data) in the ACF signal unresolved seasonal autocorrelation and indicate that a seasonal differencing step or SARIMA component is needed.

The AR Component: Learning from the Past

An autoregressive model of order $p$, written AR($p$), expresses the current value as a weighted sum of the $p$ most recent values plus white noise $\varepsilon_t \sim \mathcal{N}(0, \sigma^2)$:

For the process to be stationary, the roots of the characteristic polynomial $1 - \phi_1 B - \cdots - \phi_p B^p = 0$ must all lie outside the unit circle, where $B$ is the backshift operator ($B y_t = y_{t-1}$). An AR(1) with $|\phi_1| < 1$ is the simplest stationary AR model; $\phi_1 = 1$ is a random walk.

The MA Component: Learning from Shocks

A moving average model of order $q$, MA($q$), expresses the current value as a linear combination of the $q$ most recent white noise shocks:

Note that the "moving average" here refers to shocks, not to a simple arithmetic average. For invertibility (so that the MA process has a unique representation), the roots of $1 + \theta_1 B + \cdots + \theta_q B^q = 0$ must lie outside the unit circle.

ARMA: Combining Both

An ARMA($p$, $q$) model combines both components:

Using the backshift operator, this is compactly written as $\Phi(B) y_t = c + \Theta(B) \varepsilon_t$. ARMA models are appropriate only for stationary series.

ARIMA: Adding Integration

Most real-world series are non-stationary due to trend. ARIMA($p$, $d$, $q$) — Integrated ARMA — applies $d$ rounds of differencing before fitting ARMA($p$, $q$):

Common choices: $d=1$ for series with a linear trend, $d=2$ for series with a quadratic trend (rare in practice and prone to overfitting).

import pandas as pd
import numpy as np
from statsmodels.tsa.arima.model import ARIMA
import warnings
warnings.filterwarnings("ignore")

url = "https://raw.githubusercontent.com/jbrownlee/Datasets/master/airline-passengers.csv"
df = pd.read_csv(url, parse_dates=["Month"], index_col="Month")
df.index.freq = "MS"
log_pass = np.log(df["Passengers"])

# Fit ARIMA(2,1,1) on log-transformed series (d=1 to handle trend)
model = ARIMA(log_pass, order=(2, 1, 1))
result = model.fit()
print(result.summary())

SARIMA: Extending to Seasonality

For seasonal data, the Seasonal ARIMA model — SARIMA($p$,$d$,$q$)($P$,$D$,$Q$)$_m$ — adds seasonal AR, differencing, and MA terms at lag multiples of the seasonal period $m$:

where $\Phi_s$, $\Theta_s$ are the seasonal polynomials of orders $P$ and $Q$.

For monthly AirPassengers data ($m=12$), a common starting model is SARIMA(1,1,1)(1,1,1)$_{12}$ — first difference for trend, seasonal difference for the $\times 12$ cycle, plus one AR/MA term at each timescale.

from statsmodels.tsa.statespace.sarimax import SARIMAX

# SARIMA(1,1,1)(1,1,1)_12 — the 'airline model'
sarima = SARIMAX(
    log_pass,
    order=(1, 1, 1),
    seasonal_order=(1, 1, 1, 12),
    enforce_stationarity=True,
    enforce_invertibility=True
)
sar_result = sarima.fit(disp=False)
print(f"AIC: {sar_result.aic:.2f}  BIC: {sar_result.bic:.2f}")

Information Criteria for Model Selection

With multiple candidate orders, use AIC (Akaike Information Criterion) or BIC (Bayesian Information Criterion):

where $k$ is the number of parameters and $\hat{L}$ is the maximized likelihood. Lower is better. BIC penalizes complexity more heavily and is preferred when parsimony matters. These criteria balance goodness-of-fit against model complexity, preventing overfitting.

Common Pitfalls

• Over-differencing: applying <!--MATHBLOCK47--> when <!--MATHBLOCK48--> suffices inflates variance and introduces spurious MA terms.
• Seasonal without non-seasonal differencing: if there is both a trend and seasonality, you usually need <!--MATHBLOCK49--> and <!--MATHBLOCK50-->.
• Large <!--MATHBLOCK51--> or <!--MATHBLOCK52-->: ARIMA models are typically parsimonious (orders 0–3). Large orders usually indicate a poor transformation choice.

The Intuition Behind Exponential Smoothing

All ARIMA models estimate their parameters by maximum likelihood on a fixed parametrization. Exponential smoothing takes a different route: it adapts forecasts by applying exponentially declining weights to past observations, giving more influence to recent data. The family ranges from simple single-parameter smoothing to the full Holt-Winters triple exponential smoothing model.

Simple Exponential Smoothing (SES)

For a series with no trend or seasonality, the level estimate $\ell_t$ is updated as:

The forecast for all future horizons is the constant $\hat{y}_{T+h|T} = \ell_T$. The smoothing parameter $\alpha$ controls adaptability: large $\alpha$ tracks the data closely (responsive, noisy); small $\alpha$ smooths heavily (slow, stable). The name "exponential" comes from expanding the recursion, which shows that the weight on observation $y_{t-k}$ is $\alpha(1-\alpha)^k$ — a geometric (exponential) decay.

Holt's Linear Trend Method

When the series has a linear trend, Holt's method adds a trend component $b_t$:

Here $\beta^*$ is a second smoothing parameter for the trend. This produces linear extrapolation — fine for short horizons, but tends to over-forecast long-term. Damped trend (Gardner and McKenzie, 1985) multiplies the trend by a damping parameter $0 < \phi < 1$, yielding a more conservative long-horizon forecast:

Damped trend is empirically one of the best-performing forecasting methods across a wide range of series.

Holt-Winters Seasonal Model

The full Holt-Winters model adds a seasonal component $s_t$ with period $m$. In the additive variant:

In the multiplicative variant, $s_t$ scales the level rather than offsetting it. Choose multiplicative when seasonal amplitude is proportional to the series level.

from statsmodels.tsa.holtwinters import ExponentialSmoothing
import pandas as pd
import numpy as np

url = "https://raw.githubusercontent.com/jbrownlee/Datasets/master/airline-passengers.csv"
df = pd.read_csv(url, parse_dates=["Month"], index_col="Month")
df.index.freq = "MS"

# Multiplicative Holt-Winters (seasonal amplitude grows with level)
hw = ExponentialSmoothing(
    df["Passengers"],
    trend="add",
    seasonal="mul",
    seasonal_periods=12,
    damped_trend=True  # conservative long-horizon behaviour
)
hw_fit = hw.fit(optimized=True)
print(hw_fit.summary())
print("\nSmoothing params: alpha={:.3f}, beta={:.3f}, gamma={:.3f}, phi={:.3f}".format(
    hw_fit.params["smoothing_level"],
    hw_fit.params["smoothing_trend"],
    hw_fit.params["smoothing_seasonal"],
    hw_fit.params["damping_trend"],
))

The ETS Taxonomy

The modern ETS framework (Hyndman et al.) unifies all exponential smoothing variants under a state-space representation. ETS stands for Error / Trend / Seasonality, where each component can be:

• Error: Additive (A) or Multiplicative (M).
• Trend: None (N), Additive (A), Additive Damped (Ad).
• Seasonality: None (N), Additive (A), Multiplicative (M).

This gives 30 possible model combinations (some are ill-defined). The notation ETS(A,Ad,M) means additive errors, damped additive trend, multiplicative seasonality. This is equivalent to multiplicative Holt-Winters with damped trend — often the best default for business series.

Why ETS and ARIMA Are Related

Many ETS models have equivalent ARIMA representations. For example, SES is equivalent to ARIMA(0,1,1). Holt's linear trend is equivalent to ARIMA(0,2,2). This duality means the two families often produce similar forecasts, but ETS parameters are more interpretable for practitioners who think in terms of "how fast should the level adapt?" rather than lag polynomials.

# Forecast 24 months ahead with multiplicative Holt-Winters
forecast = hw_fit.forecast(24)
print(forecast)

When to Use ETS vs. ARIMA

• Use ETS when you want interpretable component-level parameters, when the series is short (ETS is sample-efficient), or when the decomposition framing (level/trend/season) aligns with your domain.
• Use ARIMA when you need to explicitly model autocorrelation structure, add exogenous regressors (SARIMAX), or when the autocorrelation signature is complex.
• In practice: run both, compare out-of-sample error, and ensemble if the improvement is material.

Fitting in statsmodels: MLE and State Space

statsmodels fits ARIMA models using Maximum Likelihood Estimation (MLE) via a state-space representation. The Kalman filter recursively computes the likelihood, handling missing values naturally. SARIMAX is the general class that supports seasonal models and exogenous regressors; plain ARIMA is a convenience wrapper for the non-seasonal case.

from statsmodels.tsa.statespace.sarimax import SARIMAX
import pandas as pd
import numpy as np
import warnings
warnings.filterwarnings("ignore")

url = "https://raw.githubusercontent.com/jbrownlee/Datasets/master/airline-passengers.csv"
df = pd.read_csv(url, parse_dates=["Month"], index_col="Month")
df.index.freq = "MS"
log_p = np.log(df["Passengers"])

fit = SARIMAX(log_p, order=(1,1,1), seasonal_order=(1,1,1,12)).fit(disp=False)
print(fit.summary())

The summary table reports: parameter estimates and their standard errors, $t$-statistics, $p$-values, log-likelihood, AIC, BIC, and Hannan-Quinn criterion. Coefficient $p$-values above 0.05 signal an order that may be reduced.

Residual Diagnostics: The Four-Panel Test

After fitting, the residuals $\hat{\varepsilon}_t = y_t - \hat{y}_{t|t-1}$ should look like white noise — zero mean, constant variance, and no autocorrelation. Diagnostic violations mean the model has not captured all the structure in the data.

fig = fit.plot_diagnostics(figsize=(14, 8))
fig.suptitle("SARIMA Residual Diagnostics", y=1.01)
plt.tight_layout()
plt.show()

The four panels produced by plot<em>diagnostics are:

1. Standardized residuals over time: should be centred at zero with no obvious trends, shifts, or growing variance.
2. Histogram + KDE vs. Normal: residuals should be approximately Gaussian.
3. Normal Q-Q plot: points should lie on the 45-degree line; heavy tails signal leptokurtosis.
4. Correlogram (ACF of residuals): all bars should be within the 95% confidence bands — any significant spike indicates unmodelled autocorrelation.

Formal Tests on Residuals

from statsmodels.stats.diagnostic import acorr_ljungbox
from scipy.stats import normaltest

# Ljung-Box test: H0 = no autocorrelation up to lag K
lb = acorr_ljungbox(fit.resid, lags=[10, 20, 30], return_df=True)
print("Ljung-Box p-values:")
print(lb[["lb_stat", "lb_pvalue"]].round(4))
# All p-values should be > 0.05 for a well-specified model

# Normality test (D'Agostino-Pearson)
stat, p = normaltest(fit.resid.dropna())
print(f"\nNormality test: stat={stat:.3f}, p={p:.4f}")

Common violations and remedies:

• ACF spike at seasonal lag: add <!--MATHBLOCK4--> or <!--MATHBLOCK5--> seasonal component.
• Growing residual variance: apply Box-Cox / log transformation.
• Non-normal residuals with heavy tails: consider Student-<!--MATHBLOCK6--> innovations (available via SARIMAX(innovations</em>order=...)), or accept the non-normality and use bootstrap intervals.
• Structural breaks in residuals: add outlier dummies or switch to a model with regime changes.

Automated Order Selection with pmdarima

Manual ACF/PACF identification is an art; for large-scale or automated pipelines, pmdarima.auto<em>arima implements a stepwise AIC-guided search similar to R's auto.arima:

import pmdarima as pm

log_p = np.log(df["Passengers"])

auto_model = pm.auto_arima(
    log_p,
    start_p=0, max_p=3,
    start_q=0, max_q=3,
    d=None,           # auto-select d via KPSS test
    seasonal=True,
    m=12,
    D=None,           # auto-select D
    start_P=0, max_P=2,
    start_Q=0, max_Q=2,
    information_criterion="aic",
    stepwise=True,    # stepwise search (faster than exhaustive)
    trace=True,
    error_action="ignore",
    suppress_warnings=True
)
print(auto_model.summary())

auto</em>arima uses the following strategy: fit an initial model, then explore neighbouring orders by incrementing/decrementing $p$, $q$, $P$, $Q$ one at a time, keeping whichever has lower AIC. Stepwise search is $O(pq)$ rather than $O(2^{pq})$ exhaustive search — important when the series is long or the pipeline processes many series.

Pitfalls of Automated Selection

• The selected order is optimal in-sample. Always validate out-of-sample.
• auto_arima may select <!--MATHBLOCK13--> on noisy series — check that this is genuinely needed, not over-differencing.
• AIC rewards likelihood; RMSE rewards point forecast accuracy. These are not identical. Use out-of-sample CV as the final arbiter.
• The stepwise heuristic can miss global optima. Use stepwise=False for small series where you can afford exhaustive search.

Why In-Sample Fit Is Not Enough

A model can achieve a perfect in-sample fit by memorizing noise, but such a model will forecast poorly. Honest evaluation requires held-out data — either a single holdout set or, preferably, a systematic walk-forward procedure that simulates how the model would have performed as it was deployed over time.

Point Forecast Error Metrics

Let $y_t$ be the actual value and $\hat{y}_t$ the forecast at horizon $h$. The standard metrics are:

Mean Absolute Error (MAE):

MAE is robust to outliers and has the same units as the series — easy to communicate to stakeholders.

Root Mean Squared Error (RMSE):

RMSE penalises large errors more heavily than MAE. Use when large errors are disproportionately costly.

Mean Absolute Percentage Error (MAPE):

Scale-independent — useful for comparing across series of different magnitudes. Pitfall: undefined when $y_t = 0$ and biased (penalises over-forecasts more than under-forecasts). The symmetric MAPE (sMAPE) partially addresses the asymmetry.

Mean Absolute Scaled Error (MASE):

The denominator is the MAE of the seasonal naive forecast on the training set. MASE < 1 means the model beats the seasonal naive baseline. MASE is the recommended scale-independent metric for forecast competitions.

import numpy as np

def forecast_metrics(actual, predicted):
    e = np.array(actual) - np.array(predicted)
    mae  = np.mean(np.abs(e))
    rmse = np.sqrt(np.mean(e**2))
    mape = np.mean(np.abs(e / np.array(actual))) * 100
    return {"MAE": round(mae, 3), "RMSE": round(rmse, 3), "MAPE%": round(mape, 3)}

Walk-Forward (Time-Series Cross-Validation) Backtesting

A single train/test split is sensitive to the particular test window chosen. Walk-forward validation produces multiple test windows:

• Expanding window: training set grows by one step each iteration; test is the next <!--MATHBLOCK9--> observations. Simulates a model that is periodically refitted as new data arrives.
• Rolling (sliding) window: training set size is fixed; both start and end roll forward. Useful when older data is less relevant (concept drift).

import pandas as pd
import numpy as np
from statsmodels.tsa.statespace.sarimax import SARIMAX
import warnings
warnings.filterwarnings("ignore")

url = "https://raw.githubusercontent.com/jbrownlee/Datasets/master/airline-passengers.csv"
df = pd.read_csv(url, parse_dates=["Month"], index_col="Month")
df.index.freq = "MS"
log_p = np.log(df["Passengers"])

def walk_forward_arima(series, train_size, horizon=1, order=(1,1,1),
                       seasonal_order=(1,1,1,12)):
    """Expanding-window walk-forward CV."""
    predictions, actuals = [], []
    T = len(series)
    for t in range(train_size, T - horizon + 1):
        train = series.iloc[:t]
        try:
            fit = SARIMAX(train, order=order,
                          seasonal_order=seasonal_order).fit(disp=False)
            fc = fit.forecast(horizon)
            predictions.append(float(fc.iloc[-1]))
            actuals.append(float(series.iloc[t + horizon - 1]))
        except Exception:
            pass
    return np.array(actuals), np.array(predictions)

actuals, preds = walk_forward_arima(log_p, train_size=100, horizon=1)
errors = actuals - preds
print(f"Walk-forward MAE (log scale): {np.mean(np.abs(errors)):.4f}")
print(f"Walk-forward RMSE (log scale): {np.sqrt(np.mean(errors**2)):.4f}")
# Back-transform to original scale
print(f"Walk-forward MAPE (original): {np.mean(np.abs(np.exp(actuals)-np.exp(preds))/np.exp(actuals))*100:.2f}%")

Always Beat the Naive Baseline

Before reporting any metric, compare against the seasonal naive forecast: $\hat{y}_{T+h} = y_{T+h-m}$ (repeat the value from the same season last year). Any serious model should beat this benchmark — if it does not, something is wrong with the modelling or the data.

# Seasonal naive forecast as baseline
def seasonal_naive_errors(series, train_size, m=12):
    preds = [float(series.iloc[t - m]) for t in range(train_size, len(series))]
    actuals = [float(series.iloc[t]) for t in range(train_size, len(series))]
    e = np.array(actuals) - np.array(preds)
    return {"MAE": np.mean(np.abs(e)):.4f, "RMSE": np.sqrt(np.mean(e**2)):.4f}

A Note on Prophet

Prophet (Taylor & Letham, Meta, 2018) models a time series as:

where $g(t)$ is a piecewise-linear or logistic growth trend with automatic change-point detection, $s(t)$ is a Fourier-series seasonal component, and $h(t)$ models holidays/special events. Prophet is fitted via Stan (Bayesian MAP), requires no stationarity transformation, handles missing data, and exposes intuitive parameters (changepoint<em>prior</em>scale, seasonality<em>prior</em>scale). It excels for daily business time series with strong human-driven seasonality (weekly + yearly) and is less appropriate for series where fine-grained autocorrelation structure matters (high-frequency or highly irregular data). When evaluating classical vs. modern methods, always include Prophet as a strong baseline.

from prophet import Prophet
import pandas as pd

url = "https://raw.githubusercontent.com/jbrownlee/Datasets/master/airline-passengers.csv"
df = pd.read_csv(url, parse_dates=["Month"], index_col="Month").reset_index()
df.columns = ["ds", "y"]

m = Prophet(yearly_seasonality=True, weekly_seasonality=False,
            daily_seasonality=False, changepoint_prior_scale=0.05)
m.fit(df)
future = m.make_future_dataframe(periods=24, freq="MS")
forecast = m.predict(future)
print(forecast[["ds", "yhat", "yhat_lower", "yhat_upper"]].tail(6))