9.2 Forecasting

A forecast extends what a series has shown so far: its level, trend, season and memory. Good forecasting is mostly modest: compare against simple benchmarks, check that the residuals hold no more memory, and judge on a future the model never saw.

import warnings

import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import numpy as np
from scipy import stats
import pandas as pd
import statsmodels.api as sm
from matplotlib.axes import Axes
from statsmodels.stats.diagnostic import acorr_ljungbox
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.arima_process import ArmaProcess
from statsmodels.tsa.holtwinters import ExponentialSmoothing
from statsmodels.tsa.statespace.mlemodel import MLEResults
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.tsa.stattools import acf, pacf

from statly import style

style.use()
rng = np.random.default_rng(1)
warnings.filterwarnings("ignore")  # the fitting routines warn about starting values

co2 = sm.datasets.co2.load_pandas().data.interpolate()  # weekly CO2 at Mauna Loa, 1958 to 2001
monthly = co2.co2.resample("MS").mean()
train, actual = monthly[:"1997"], monthly["1998":]
HORIZON = len(actual)
sun = sm.datasets.sunspots.load_pandas().data  # sunspots per year, 1700 to 2008
spots = pd.Series(
    sun.SUNACTIVITY.to_numpy(), index=pd.period_range("1700", periods=len(sun), freq="Y")
)


def plot_acf_bars(
    ax: Axes, values: np.ndarray, lags: int, partial: bool = False, color: str = style.BLUE
) -> None:
    """Autocorrelation (or partial autocorrelation) bars with the band chance alone stays inside."""
    correlations = (pacf if partial else acf)(values, nlags=lags)[1:]
    ax.vlines(range(1, lags + 1), 0, correlations, color=color, lw=2)
    band = 2 / np.sqrt(len(values))
    ax.axhspan(-band, band, color=style.GREY, alpha=0.2, lw=0)
    ax.axhline(0, color=style.INK, lw=0.8)
    ax.set(ylim=(-1, 1), xlabel="lag")

9.2.1 Start with benchmarks

Before any model, four simple forecasts set the bar. With $y_T$ the last value seen and h the number of steps ahead:

• mean: the average of all the past, $\hat y_{T+h} = \bar y$;
• naive: the last value, $\hat y_{T+h} = y_T$, the best forecast for a random walk;
• seasonal naive: the value from the same season last year, $\hat y_{T+h} = y_{T+h-12}$;
• drift: the last value plus the average change per step, the line from the first to the last point continued.

history = train.to_numpy()
months = np.arange(1, HORIZON + 1)
BENCHMARKS = {
    "mean": np.full(HORIZON, history.mean()),
    "naive": np.full(HORIZON, history[-1]),
    "seasonal naive": np.tile(history[-12:], HORIZON // 12 + 1)[:HORIZON],
    "drift": history[-1] + months * (history[-1] - history[0]) / (len(history) - 1),
}
COLORS = {
    "mean": style.GOLD,
    "naive": style.PURPLE,
    "seasonal naive": style.GREEN,
    "drift": style.CORAL,
}


def plot_benchmarks(ax: Axes) -> None:
    """Recent history, the four benchmarks, and what happened."""
    recent = train["1993":]
    ax.plot(recent.index, recent, color=style.BLUE, lw=1.2)
    ax.plot(actual.index, actual, color=style.GREY, lw=1.6)
    for name, forecast in BENCHMARKS.items():
        mae = np.mean(abs(actual.to_numpy() - forecast))
        ax.plot(
            actual.index, forecast, color=COLORS[name], lw=1.4, label=f"{name}: MAE {mae:.1f}"
        )
    ax.legend(fontsize=8, loc="upper left")
    ax.set(xlabel="year", ylabel="CO₂ (ppm)", ylim=(350, 375))


fig, ax = plt.subplots(figsize=(8, 3.6))
plot_benchmarks(ax)
plt.show()

Figure 9.2.1: Four benchmark forecasts of monthly CO₂ for 1998 to 2001, made from the data up to 1997 (blue). What actually happened is grey. The mean (far below, off the chart) and the naive forecast miss the trend or the season; drift catches the trend; seasonal naive catches the season but lags a year behind the trend. Their mean absolute errors (MAE, ppm) are in the legend.

A model that cannot beat the best of these is not worth its complexity. Here each benchmark catches only part of the pattern, which suggests the model must handle trend and season.

9.2.2 Autoregressive models

An autoregressive model, AR(p), predicts each value from the p values before it:

$$y_t = c + \phi_1 y_{t-1} + \phi_2 y_{t-2} + \dots + \phi_p y_{t-p} + \varepsilon_t$$

a regression of the series on its own past. The coefficients shape the memory. AR(1) with φ = 0.9 drifts slowly and returns to its mean; with φ = −0.6 it zigzags; AR(2) can produce cycles. Each has a signature in its ACF and PACF (Chapter 9.1): the PACF cuts off after lag p, because once the last p values are known, older ones add nothing.

PROCESSES = [
    ("AR(1), φ = 0.9", [1, -0.9], [1], style.BLUE),
    ("AR(1), φ = −0.6", [1, 0.6], [1], style.GOLD),
    ("AR(2), cycles", [1, -1.3, 0.7], [1], style.PURPLE),
]

fig, axes = plt.subplots(3, 3, figsize=(10.5, 6.2))
for col, (name, ar, ma, color) in enumerate(PROCESSES):
    values = ArmaProcess(ar, ma).generate_sample(300, distrvs=rng.standard_normal)
    axes[0, col].plot(values, color=color, lw=0.9)
    axes[0, col].set(title=name, xlabel="time")
    plot_acf_bars(axes[1, col], values, 20, color=color)
    plot_acf_bars(axes[2, col], values, 20, partial=True, color=color)
axes[1, 0].set_ylabel("ACF")
axes[2, 0].set_ylabel("PACF")
fig.tight_layout()
plt.show()

Figure 9.2.2: Three simulated autoregressive series (top), their ACF (middle) and PACF (bottom). AR(1) with φ = 0.9: slowly fading ACF, one PACF spike. AR(1) with φ = −0.6: alternating ACF, one negative PACF spike. AR(2) with φ₁ = 1.3, φ₂ = −0.7: a damped wave in the ACF, like the sunspots, and two PACF spikes.

An AR model is stationary only if its memory fades: for AR(1), |φ| < 1. At φ = 1 it is the random walk, the unit root that the Dickey-Fuller test looks for (Chapter 9.1).

9.2.3 Moving-average and ARMA models

A moving-average model, MA(q), lets the value depend on the last q surprises instead of the last values:

$$y_t = c + \varepsilon_t + \theta_1 \varepsilon_{t-1} + \dots + \theta_q \varepsilon_{t-q}$$

A shock affects exactly q future values and is then forgotten, so the ACF cuts off after lag q. An ARMA(p, q) model combines both, often with fewer coefficients than either alone. The signatures guide the choice:

ma_values = ArmaProcess([1], [1, 0.8]).generate_sample(300, distrvs=rng.standard_normal)

fig, axes = style.subplots(3, height=2.8, width=10.5)
axes[0].plot(ma_values, color=style.GREEN, lw=0.9)
axes[0].set(title="MA(1), θ = 0.8", xlabel="time")
plot_acf_bars(axes[1], ma_values, 20, color=style.GREEN)
axes[1].set_title("ACF")
plot_acf_bars(axes[2], ma_values, 20, partial=True, color=style.GREEN)
axes[2].set_title("PACF")
fig.tight_layout()
plt.show()

Figure 9.2.3: A simulated MA(1) series with θ = 0.8 (left), its ACF with a single spike at lag 1 (middle), and its PACF, which fades with alternating signs (right): the mirror image of an AR(1).

Model	ACF	PACF
AR(p)	fades (exponentially or as a damped wave)	cuts off after lag p
MA(q)	cuts off after lag q	fades
ARMA(p, q)	fades after lag q	fades after lag p

Real series rarely show textbook patterns, so in practice a few candidate orders are fitted and compared by the AIC (Chapter 2.3), then checked.

9.2.4 Fitting and forecasting an AR model

The sunspots’ PACF (Chapter 9.1) had large spikes at lags 1 and 2 and smaller ones up to 9. Fitting AR models of every order from 1 to 12 and comparing their AIC:

spots_train, spots_test = spots[:"1959"], spots["1960":"1990"]
orders = range(1, 13)
aics = [ARIMA(spots_train, order=(p, 0, 0)).fit().aic for p in orders]
best_order = orders[int(np.argmin(aics))]
fits = {p: ARIMA(spots_train, order=(p, 0, 0)).fit() for p in (2, best_order)}


def plot_aic(ax: Axes) -> None:
    """AIC against the AR order."""
    ax.plot(list(orders), aics, "o-", color=style.BLUE, ms=4)
    ax.plot(best_order, min(aics), "o", ms=10, color=style.CORAL, fillstyle="none", mew=2)
    ax.set(title="Choosing the order", xlabel="AR order p", ylabel="AIC", xticks=list(orders))


def plot_spot_forecasts(ax: Axes) -> None:
    """Forecasts from two AR models, with the interval of the better one."""
    years = np.arange(1960, 1991)
    recent = spots_train["1900":]
    ax.plot(recent.index.year, recent, color=style.BLUE, lw=1)
    ax.plot(years, spots_test, color=style.GREY, lw=1.6)
    for p, color in ((2, style.GOLD), (best_order, style.CORAL)):
        frame = fits[p].get_forecast(len(years)).summary_frame()
        if p == best_order:
            ax.fill_between(
                years, frame.mean_ci_lower, frame.mean_ci_upper, color=color, alpha=0.15, lw=0
            )
        ax.plot(years, frame["mean"], color=color, lw=1.6, label=f"AR({p})")
    ax.legend(fontsize=8, loc="upper left")
    ax.set(title="Forecasts", xlabel="year", ylabel="sunspots")


fig = plt.figure(figsize=(10.5, 3.3))
grid = fig.add_gridspec(1, 2, width_ratios=[1, 1.6])
plot_aic(fig.add_subplot(grid[0]))
plot_spot_forecasts(fig.add_subplot(grid[1]))
fig.tight_layout()
plt.show()

Figure 9.2.4: Left: AIC of AR models of order 1 to 12 for the sunspots (up to 1959), lowest at order 9. Right: forecasts for 1960 to 1990 from AR(2) and AR(9), with the 95% interval of AR(9), against what happened (grey). Both forecasts are damped waves that fade towards the long-run mean; the uncertainty grows with the horizon until it equals the whole range of the series.

For a stationary series, every forecast eventually returns to the mean, and its interval to the spread of the series: far enough ahead, the past has nothing more to say. The useful horizon is the stretch before that happens.

9.2.5 ARIMA and seasonal ARIMA

Non-stationary series are differenced first (Chapter 9.1). An ARIMA(p, d, q) model is an ARMA(p, q) model for the series differenced d times; the “I” stands for integrated. Seasonal series add a second set of terms at the seasonal lag: ARIMA(p, d, q)(P, D, Q)ₛ, where D is the number of seasonal differences and P and Q are AR and MA terms at lags s, 2s, … For monthly CO₂ a classic choice is ARIMA(1, 1, 1)(0, 1, 1)₁₂: one ordinary and one seasonal difference, as in Chapter 9.1, plus a little short-term and yearly memory.

seasonal = SARIMAX(train, order=(1, 1, 1), seasonal_order=(0, 1, 1, 12)).fit(disp=False)
forecast = seasonal.get_forecast(HORIZON).summary_frame(alpha=0.05)


def plot_forecast(ax: Axes) -> None:
    """The recent data, the forecast with its band, and what actually happened."""
    recent = train["1992":]
    ax.plot(recent.index, recent, color=style.BLUE, lw=1.2, label="data used")
    ax.fill_between(
        forecast.index,
        forecast.mean_ci_lower,
        forecast.mean_ci_upper,
        color=style.CORAL,
        alpha=0.2,
        lw=0,
    )
    ax.plot(
        forecast.index,
        forecast["mean"],
        color=style.CORAL,
        lw=1.6,
        label=f"forecast: MAE {np.mean(abs(actual - forecast['mean'])):.1f}",
    )
    ax.plot(actual.index, actual, color=style.GREY, lw=1.2, label="what happened")
    ax.legend(fontsize=8, loc="upper left")
    ax.set(xlabel="year", ylabel="CO₂ (ppm)")


fig, ax = plt.subplots(figsize=(8, 3.4))
plot_forecast(ax)
plt.show()

Figure 9.2.5: A seasonal ARIMA(1,1,1)(0,1,1)₁₂ model fitted to monthly CO₂ up to 1997, forecasting four years ahead (red) with a 95% band. The actual values (grey) stay close to the forecast. The band widens with the horizon: the further ahead, the less certain.

The model’s mean absolute error is far below that of every benchmark in Section 9.2.1.

9.2.6 Checking the residuals

If a model has captured the memory of a series, its one-step-ahead residuals are white noise: no autocorrelation left. The ACF of the residuals shows it, and the Ljung-Box test sums up the first h autocorrelations $r_k$ of the residuals:

$$Q = n(n + 2) \sum_{k=1}^{h} \frac{r_k^2}{n - k}$$

which follows a chi-square distribution if the residuals are white noise. A small p-value says there is memory the model missed.

too_simple = ARIMA(train, order=(0, 1, 0)).fit()


def plot_residual_acf(ax: Axes, model: MLEResults, title: str, color: str) -> None:
    """ACF of a model's residuals with the Ljung-Box p-value."""
    residuals = model.resid[13:]  # skip the start, where differencing has no history
    plot_acf_bars(ax, residuals, 36, color=color)
    p = acorr_ljungbox(residuals, lags=[24]).lb_pvalue.iloc[0]
    style.note(
        ax,
        f"Ljung-Box p = {p:.2f}" if p >= 0.001 else "Ljung-Box p < 0.001",
        x=0.55,
        y=0.98,
        color=style.INK,
    )
    ax.set_title(title)


fig, axes = style.subplots(2, height=3, width=9.5, sharey=True)
plot_residual_acf(axes[0], too_simple, "ARIMA(0,1,0)", style.GOLD)
plot_residual_acf(axes[1], seasonal, "ARIMA(1,1,1)(0,1,1)₁₂", style.BLUE)
axes[0].set_ylabel("autocorrelation of residuals")
fig.tight_layout()
plt.show()

Figure 9.2.6: ACF of the residuals of two models for monthly CO₂, with the Ljung-Box p-value over the first 24 lags. Left: a model with only an ordinary difference, ARIMA(0, 1, 0), leaves the whole yearly season in its residuals. Right: the seasonal ARIMA leaves only small, scattered autocorrelations, consistent with white noise.

Where do those p-values come from? For residuals that really are white noise, Q varies around its number of lags, 24, following a chi-square distribution with 24 degrees of freedom. Simulating white noise of the same length shows it:

good_residuals = seasonal.resid[13:]
bad_residuals = too_simple.resid[13:]
q_good = float(acorr_ljungbox(good_residuals, lags=[24]).lb_stat.iloc[0])
q_bad = float(acorr_ljungbox(bad_residuals, lags=[24]).lb_stat.iloc[0])
null_q = np.array(
    [
        acorr_ljungbox(rng.normal(size=len(good_residuals)), lags=[24]).lb_stat.iloc[0]
        for _ in range(2000)
    ]
)

fig, ax = plt.subplots(figsize=(8, 3.2))
style.null_panel(
    ax,
    q_good,
    samples=null_q,
    dist=stats.chi2(24),
    tail="right",
    xlabel="Ljung-Box Q over 24 lags",
)
ax.annotate(
    f"ARIMA(0,1,0): Q = {q_bad:,.0f} →",
    xy=(0.97, 0.5),
    xycoords="axes fraction",
    ha="right",
    color=style.GOLD,
    fontsize=9,
)
plt.show()

Figure 9.2.7: The Ljung-Box statistic Q over 24 lags for 2,000 simulated white-noise series as long as the residuals (grey), with the chi-square curve with 24 degrees of freedom (blue). The seasonal ARIMA’s Q (dashed) lies in the bulk: nothing left to find. The simple model’s Q, from its large yearly autocorrelations, is far beyond the right edge of this plot.

With long series, the Ljung-Box test rejects for tiny leftover correlations; as always, look at the size of what is left, not only at the p-value.

9.2.7 Exponential smoothing

A second family of forecasting methods weighs the past directly. Simple exponential smoothing keeps a running level,

$$\ell_t = \alpha\, y_t + (1 - \alpha)\, \ell_{t-1}$$

and forecasts that level. Unrolled, the level is a weighted average of all past values with weights $\alpha, \alpha(1-\alpha), \alpha(1-\alpha)^2, \dots$ that fade exponentially. Holt’s method adds a smoothed trend, and Holt-Winters adds a smoothed season, each with its own smoothing constant, all fitted by minimising the one-step-ahead errors.

holt_winters = ExponentialSmoothing(
    train, trend="add", seasonal="add", seasonal_periods=12
).fit()
hw_forecast = holt_winters.forecast(HORIZON)


def plot_weights(ax: Axes) -> None:
    """Weights on past values for three smoothing constants."""
    lags = np.arange(0, 20)
    for alpha, color in ((0.2, style.BLUE), (0.5, style.GREEN), (0.8, style.CORAL)):
        ax.plot(
            -lags, alpha * (1 - alpha) ** lags, "o-", ms=3, color=color, label=f"α = {alpha}"
        )
    ax.legend(fontsize=8)
    ax.set(title="Weights on the past", xlabel="steps back", ylabel="weight")


def plot_holt_winters(ax: Axes) -> None:
    """Holt-Winters and seasonal ARIMA forecasts against what happened."""
    ax.plot(actual.index, actual, color=style.GREY, lw=1.6, label="what happened")
    ax.plot(
        actual.index,
        hw_forecast,
        color=style.CORAL,
        lw=1.4,
        label=f"Holt-Winters: MAE {np.mean(abs(actual - hw_forecast)):.2f}",
    )
    ax.plot(
        actual.index,
        forecast["mean"],
        color=style.BLUE,
        lw=1.4,
        ls="--",
        label=f"seasonal ARIMA: MAE {np.mean(abs(actual - forecast['mean'])):.2f}",
    )
    ax.legend(fontsize=8, loc="upper left")
    ax.xaxis.set_major_locator(mdates.YearLocator())
    ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y"))
    ax.set(title="Two methods", xlabel="year", ylabel="CO₂ (ppm)")


fig = plt.figure(figsize=(10.5, 3.3))
grid = fig.add_gridspec(1, 2, width_ratios=[1, 1.5])
plot_weights(fig.add_subplot(grid[0]))
plot_holt_winters(fig.add_subplot(grid[1]))
fig.tight_layout()
plt.show()

Figure 9.2.8: Left: the weights exponential smoothing gives to past values, for three smoothing constants. A large α follows the latest values; a small one averages over a long past. Right: Holt-Winters forecasts (additive trend and season, red) for CO₂ from 1998, next to the seasonal ARIMA (blue) and what happened (grey). The two methods agree closely.

9.2.8 Judging forecasts honestly

A forecast must be judged on a future the model did not see, and one test period is a single lucky or unlucky draw. The rolling origin of Chapter 6.3 repeats the exercise: fit on the data up to a date, forecast the next two years, move the date forward six months, and repeat. Then compare the errors at each horizon (one month ahead, two months ahead, …). Common error measures for forecasts $\hat y$ of values y:

$$\text{MAE} = \text{mean}\,|y - \hat y|, \qquad \text{RMSE} = \sqrt{\text{mean}\,(y - \hat y)^2}, \qquad \text{MAPE} = \text{mean}\,\left|\frac{y - \hat y}{y}\right| \times 100\%$$

ORIGINS = pd.date_range("1990-01-01", "1995-07-01", freq="6MS")
AHEAD = 24
errors = {"seasonal naive": [], "Holt-Winters": [], "seasonal ARIMA": []}
for origin in ORIGINS:
    past = monthly[: origin - pd.offsets.MonthBegin()]
    future = monthly[origin:].iloc[:AHEAD].to_numpy()
    naive = np.tile(past.to_numpy()[-12:], 3)[:AHEAD]
    hw = (
        ExponentialSmoothing(past, trend="add", seasonal="add", seasonal_periods=12)
        .fit()
        .forecast(AHEAD)
        .to_numpy()
    )
    sa = (
        SARIMAX(past, order=(1, 1, 1), seasonal_order=(0, 1, 1, 12))
        .fit(disp=False)
        .forecast(AHEAD)
        .to_numpy()
    )
    for name, prediction in (
        ("seasonal naive", naive),
        ("Holt-Winters", hw),
        ("seasonal ARIMA", sa),
    ):
        errors[name].append(abs(future - prediction))

fig, ax = plt.subplots(figsize=(8, 3.3))
for (name, values), color in zip(errors.items(), (style.GREEN, style.CORAL, style.BLUE)):
    ax.plot(range(1, AHEAD + 1), np.mean(values, axis=0), "o-", ms=3, color=color, label=name)
ax.legend(fontsize=8)
ax.set(xlabel="months ahead", ylabel="mean absolute error (ppm)")
plt.show()

Figure 9.2.9: Rolling-origin evaluation for monthly CO₂: 12 forecast origins from 1990 to 1995, each forecasting 24 months ahead from all data before it. Mean absolute error by horizon for seasonal naive (green), Holt-Winters (red) and seasonal ARIMA (blue). Errors grow with the horizon for all three; the two models beat the benchmark at every horizon, and are close to each other.

MAPE is popular because it has no units, but it explodes for values near zero and punishes over-forecasts and under-forecasts differently. Whatever the measure, compare models on the same origins and horizons, and always include a benchmark.

9.2.9 In one picture

• Start with benchmarks: mean, naive, seasonal naive, drift.
• AR(p): $y_t = c + \sum \phi_j y_{t-j} + \varepsilon_t$, PACF cuts off after p. MA(q): ACF cuts off after q. ARMA combines them; choose orders by AIC.
• ARIMA(p, d, q)(P, D, Q)ₛ adds differencing and seasonal terms.
• Residuals should be white noise: check their ACF and the Ljung-Box test.
• Exponential smoothing (Holt-Winters) is a strong, simple alternative.
• Judge forecasts by rolling origin, at each horizon, with MAE or RMSE; intervals widen with the horizon.

Next: cycles, several series at once, changing volatility and change points.