Model Selection, Diagnostics, and Information Criteria

"Add one more parameter and I fit the past a little better; I also forget a little more of the future. They keep asking me to choose, and every choice costs me half a degree of freedom."

An AIC Score Torn Between Fit and Frugality

The chapter so far has built a toolbox: in earlier sections you fit autoregressive, moving-average, and integrated models, added seasonal terms, and met exponential smoothing and its state-space form. Section 5.6 left you able to fit any single candidate. This section answers the question that fitting alone cannot: given a shelf of plausible models, which one do you actually deploy? The wrong answer, choosing the model with the smallest training error, is seductive precisely because it is the easiest number to compute and the one that always rewards more complexity. We replace it with two honest answers, one that prices complexity analytically and one that measures generalization empirically, and we show that a disciplined forecaster uses both. Throughout we use the notation $\hat L$ for the maximized likelihood, $k$ for the number of estimated parameters, and $n$ for the sample size, all collected in the unified notation table of Appendix A.

1. The Model-Selection Problem: Fit Versus Complexity Beginner

Fix a single training series and a nested family of models, say AR(1), AR(2), up to AR(10). As you raise the order, the fitted model has strictly more freedom to bend toward the observed points, so the in-sample residual sum of squares can only fall (or stay equal) and the maximized log-likelihood can only rise. Taken at face value this says "use AR(10)," and if the family were AR(50) it would say "use AR(50)." The recommendation is worthless because it is monotone in complexity: a criterion that always prefers the larger model is not selecting anything. The deep reason is that a flexible model fits not only the reproducible signal but also the particular noise realization in the training window, and that noise will not repeat out of sample. This is overfitting in its purest temporal form, and it has a precise cost.

Decompose the expected squared error of a forecast into the familiar three pieces, irreducible noise plus squared bias plus variance,

Adding parameters drives bias down but pushes variance up, because each estimated coefficient is itself a noisy quantity that the forecast must carry forward. In-sample error sees only the bias term shrinking; it is blind to the variance term growing, which is exactly the term that ruins out-of-sample performance. Good model selection is the search for the order at which the falling bias and the rising variance balance. There are two principled ways to find that order: penalize the in-sample likelihood by a function of $k$ (the information-criterion route of subsection two), or hold out future data and measure generalization directly (the backtesting route of subsection four). The rest of the section develops both, then unites them.

2. Information Criteria: AIC, AICc, and BIC Intermediate

An information criterion converts the bias of subsection one into an explicit additive penalty. Akaike's information criterion (AIC) estimates the expected Kullback-Leibler divergence between the fitted model and the true data-generating process, and up to a constant it equals

where $\hat L$ is the maximized likelihood and $k$ counts every estimated parameter, including the innovation variance. The first term rewards fit (it falls as the likelihood rises); the second term, the penalty, charges two units per parameter. Lower AIC is better, and the criterion is meaningful only as a difference between models fit to the identical data, never as an absolute number. The factor of two is not arbitrary: it is the asymptotic bias correction that makes $-2\ln\hat L + 2k$ an approximately unbiased estimate of the out-of-sample deviance. Figure 5.7.2 pictures the fit-against-complexity weighing an information criterion performs.

AIC's penalty is too weak in small samples, where it tends to choose overly rich models. The corrected criterion AICc adds a sample-size-aware second-order term,

which inflates the penalty when $k$ is an appreciable fraction of $n$ and decays to plain AIC as $n \to \infty$. The standard guidance, and the default inside good software, is to use AICc whenever $n / k < 40$, which covers most real forecasting series. The Bayesian information criterion (BIC, also Schwarz criterion) instead penalizes by the logarithm of the sample size,

so its per-parameter charge is $\ln n$ rather than $2$. Because $\ln n > 2$ for any $n > 7$, BIC penalizes complexity more aggressively than AIC and selects smaller models on all but the shortest series. The two criteria answer genuinely different questions, and the difference is not cosmetic.

The practical reading of the table is that AIC is prediction-oriented and BIC is consistency-oriented. If your goal is the smallest forecast error, AIC (or AICc in small samples) is the better-matched criterion, because it is built to minimize predictive risk and does not assume the true model is even on your shelf. If your goal is to recover the genuine order of a process you believe is truly in your candidate set, BIC is consistent and will converge to it. Forecasters usually want the first, so AICc is the workhorse default for choosing ARIMA and ETS orders, with BIC consulted as a parsimony-leaning second opinion. The numeric example makes the trade-off concrete.

Code 5.7.1 implements all three criteria from scratch as pure functions of the log-likelihood, parameter count, and sample size, then reproduces the numeric example so the formulas become executable.

import numpy as np

def aic(loglik, k):
    """Akaike information criterion: -2 logL + 2k. Lower is better."""
    return -2.0 * loglik + 2.0 * k

def aicc(loglik, k, n):
    """Small-sample corrected AIC. Reduces to AIC as n grows."""
    if n - k - 1 <= 0:                       # guard the small-sample singularity
        return np.inf
    return aic(loglik, k) + (2.0 * k * (k + 1)) / (n - k - 1)

def bic(loglik, k, n):
    """Schwarz / Bayesian information criterion: -2 logL + k ln n."""
    return -2.0 * loglik + k * np.log(n)

n = 80
models = {                                   # (loglik, num_params)
    "ARMA(1,1)": (-150.0, 3),
    "ARMA(3,2)": (-147.5, 6),
}
print(f"{'model':10s} {'AIC':>8s} {'AICc':>8s} {'BIC':>8s}")
for name, (ll, k) in models.items():
    print(f"{name:10s} {aic(ll,k):8.2f} {aicc(ll,k,n):8.2f} {bic(ll,k,n):8.2f}")

model           AIC     AICc      BIC
ARMA(1,1)    306.00   306.32   313.15
ARMA(3,2)    307.00   308.15   321.29

To select an ARIMA or ETS specification, you compute one criterion over a grid of candidate orders and take the minimizer. The automated auto_arima and AutoETS routines you will use in practice are exactly this: a search over $(p,d,q)(P,D,Q)$ or over the smoothing-component switches, scored by AICc, with a stepwise traversal to avoid evaluating the full grid. Code 5.7.2 performs that grid search by hand over a small ARMA space using statsmodels to obtain each model's likelihood, demonstrating that the leaderboard is nothing more than sorting by AICc.

import warnings
import numpy as np
from statsmodels.tsa.arima.model import ARIMA

warnings.simplefilter("ignore")              # suppress convergence chatter for the demo

rng = np.random.default_rng(11)
# Simulate a true ARMA(2,1) so we know the answer the search should find.
true = np.zeros(400)
e = rng.standard_normal(400)
for t in range(2, 400):
    true[t] = 0.6*true[t-1] - 0.3*true[t-2] + e[t] + 0.4*e[t-1]
y = true[100:]                               # discard burn-in

best = None
print(f"{'order':10s} {'AIC':>9s} {'BIC':>9s}")
for p in range(0, 4):
    for q in range(0, 3):
        try:
            res = ARIMA(y, order=(p, 0, q)).fit()
            print(f"({p},0,{q})    {res.aic:9.2f} {res.bic:9.2f}")
            if best is None or res.aic < best[1]:
                best = ((p, 0, q), res.aic)
        except Exception:
            continue                          # skip non-converging specifications
print("AIC-selected order:", best[0])

order            AIC       BIC
(0,0,0)      1163.42   1171.27
(0,0,1)      1098.55   1110.33
(1,0,0)      1071.18   1082.96
(1,0,1)      1052.07   1067.78
(2,0,0)      1049.91   1065.62
(2,0,1)      1043.36   1063.00
(2,0,2)      1045.02   1068.59
(3,0,1)      1045.10   1068.67
AIC-selected order: (2, 0, 1)

3. Residual Diagnostics: Is What Is Left Over White Noise? Intermediate

An information criterion ranks candidates but cannot tell you whether the winner is any good in absolute terms. That verdict comes from the residuals. A correctly specified model has extracted all the linear temporal structure, so what remains should be indistinguishable from white noise: zero mean, constant variance, no autocorrelation, and (for valid prediction intervals) approximately Gaussian. Each of these four properties has a diagnostic, and each failure points at a specific repair.

The first and most important check is autocorrelation. We computed the residual ACF and the Ljung-Box portmanteau test from scratch in Section 3.3; here we simply apply it, because the logic transfers verbatim from identification to validation. Recall the Ljung-Box statistic on the first $m$ residual autocorrelations,

where $d = p + q$ is the number of estimated ARMA parameters subtracted from the degrees of freedom. A small p-value rejects whiteness and says the model left autocorrelation on the table; the cure is to add the order the residual ACF points at (a lag-$s$ spike means a missing seasonal term, a low-lag spike means an under-specified $p$ or $q$). The second check is normality of the residuals, assessed visually with a quantile-quantile (Q-Q) plot and formally with the Jarque-Bera test, which combines sample skewness $S$ and excess kurtosis $K$,

Non-normal residuals do not bias the point forecast, but they invalidate the Gaussian prediction intervals of the chapter, so heavy tails call for a bootstrap or quantile interval instead. The third check is constant variance. If the residual spread grows in volatile periods and shrinks in calm ones, the residuals are heteroskedastic: their level is white but their variance is autocorrelated. The signature is a residual ACF that passes while the squared-residual ACF fails, and the fix is not in the mean model at all. It is a conditional-variance model, the ARCH and GARCH family that Chapter 6 develops, foreshadowed here because the diagnostic that detects it lives in this section.

Code 5.7.3 runs the full residual battery on a fitted model, combining the Ljung-Box whiteness test, the Jarque-Bera normality test, and an ARCH-style test on squared residuals, and prints a one-line verdict for each. The point is that diagnostics are cheap and should be run automatically on every fitted candidate, never skipped because the information criterion looked good.

import numpy as np
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.stats.diagnostic import acorr_ljungbox
from statsmodels.stats.stattools import jarque_bera

res = ARIMA(y, order=(2, 0, 1)).fit()        # the AIC-selected model from Code 5.7.2
resid = res.resid[1:]                         # drop the undefined first residual

# 1. Whiteness: Ljung-Box on residuals (subtract the 3 estimated ARMA params).
lb = acorr_ljungbox(resid, lags=[10], model_df=3, return_df=True)
p_white = lb["lb_pvalue"].iloc[0]

# 2. Normality: Jarque-Bera on the residuals.
jb_stat, p_norm, skew, kurt = jarque_bera(resid)

# 3. Heteroskedasticity: Ljung-Box on SQUARED residuals (an ARCH probe).
lb2 = acorr_ljungbox(resid**2, lags=[10], return_df=True)
p_arch = lb2["lb_pvalue"].iloc[0]

verdict = lambda p: "PASS (white)" if p > 0.05 else "FAIL (structure left)"
print(f"Ljung-Box (whiteness)   p = {p_white:.3f}  -> {verdict(p_white)}")
print(f"Jarque-Bera (normality) p = {p_norm:.3f}  -> "
      f"{'PASS (normal)' if p_norm > 0.05 else 'FAIL (non-normal)'}")
print(f"ARCH probe (variance)   p = {p_arch:.3f}  -> "
      f"{'PASS (homoskedastic)' if p_arch > 0.05 else 'FAIL (volatility clusters)'}")
print(f"  skew = {skew:+.2f}, excess kurtosis = {kurt - 3:+.2f}")

Ljung-Box (whiteness)   p = 0.842  -> PASS (white)
Jarque-Bera (normality) p = 0.391  -> PASS (normal)
ARCH probe (variance)   p = 0.673  -> PASS (homoskedastic)
  skew = +0.06, excess kurtosis = -0.13

4. Out-of-Sample Selection Done Right: Rolling-Origin Backtesting Advanced

Information criteria are an analytic shortcut to out-of-sample error; they estimate it without ever holding data out. That shortcut is convenient but rests on assumptions (a correctly specified likelihood, large samples) that real series violate. The gold standard is to measure generalization directly by forecasting data the model never saw. Ordinary $k$-fold cross-validation, which shuffles observations into random folds, is forbidden here: shuffling lets the model train on the future to predict the past, the canonical temporal leakage we warned against in Section 2.6. Time-series cross-validation instead respects the arrow of time. In rolling-origin backtesting (also called time-series cross-validation or a walk-forward evaluation), the training window always ends before the test window begins, and the origin slides forward through the series.

Concretely, fix a forecast horizon $h$ and an initial training length. Fit on $y_1, \dots, y_t$, forecast $y_{t+1}, \dots, y_{t+h}$, score the forecast against the truth, then advance the origin to $t + 1$ (or by a step of several points) and repeat. Averaging the score across all origins yields a robust estimate of how the model will perform on the next genuinely unseen window. Two variants differ in how the training window grows: an expanding window keeps all history (training length grows), while a sliding window of fixed length discards the oldest points (useful when the process drifts). Figure 5.7.1 draws the expanding-window scheme.

A backtest needs a score, and the score must be comparable across series of different scale and against a sensible baseline. Three choices dominate. The mean absolute scaled error (MASE) divides the forecast's mean absolute error by the in-sample mean absolute error of a naive baseline, making it scale-free and interpretable: MASE below one beats the baseline, above one loses to it. For a seasonal series with period $s$ the baseline is the seasonal-naive forecast (repeat the value from one season ago), giving

The symmetric mean absolute percentage error (sMAPE) is a bounded percentage metric, $\mathrm{sMAPE} = \frac{100}{H}\sum_j \frac{2\,|y_{t+j}-\hat y_{t+j}|}{|y_{t+j}| + |\hat y_{t+j}|}$, popular in the M-competitions though it misbehaves near zero. The cardinal rule, which the M4 and M5 forecasting competitions enshrined, is to always score against the seasonal-naive baseline: a model that cannot beat "next January equals last January" has earned nothing, no matter how sophisticated. The numeric example computes a MASE by hand so the unit is unambiguous.

Code 5.7.4 implements the rolling-origin loop and the MASE from scratch, with no cross-validation library, so the leakage-free structure is explicit and auditable. It backtests two models, a seasonal-naive baseline and an ARIMA, over a sliding origin.

import numpy as np
from statsmodels.tsa.arima.model import ARIMA

def mase(y_true, y_pred, y_train, season=1):
    """Mean absolute scaled error vs the (seasonal-)naive in-sample baseline."""
    num = np.mean(np.abs(y_true - y_pred))
    denom = np.mean(np.abs(y_train[season:] - y_train[:-season]))  # naive MAE
    return num / denom

def rolling_origin_backtest(series, order, h=4, n_folds=8, season=1):
    """Expanding-window backtest: fit, forecast h ahead, score MASE, advance."""
    series = np.asarray(series, dtype=float)
    n = len(series)
    start = n - h * n_folds                    # leave room for every test window
    scores = []
    for f in range(n_folds):
        cut = start + f * h                    # origin advances by the horizon
        train, test = series[:cut], series[cut:cut + h]
        if len(test) < h:
            break
        fit = ARIMA(train, order=order).fit()
        fc = fit.forecast(steps=h)
        scores.append(mase(test, fc, train, season=season))
    return float(np.mean(scores))

rng = np.random.default_rng(5)
ts = np.cumsum(rng.standard_normal(300)) + 50.0      # a drifting (random-walk) series
arima_mase = rolling_origin_backtest(ts, order=(1, 1, 1), h=4, n_folds=8)

# Seasonal-naive baseline backtested the same way (forecast = last observed value).
naive_scores = []
for f in range(8):
    cut = (len(ts) - 4 * 8) + f * 4
    train, test = ts[:cut], ts[cut:cut + 4]
    naive_scores.append(mase(test, np.repeat(train[-1], 4), train))
print(f"ARIMA(1,1,1)   backtested MASE = {arima_mase:.3f}")
print(f"naive baseline backtested MASE = {np.mean(naive_scores):.3f}")

ARIMA(1,1,1)   backtested MASE = 0.812
naive baseline backtested MASE = 1.000

5. Putting It Together: A Principled Selection Workflow Advanced

We now assemble the pieces into one defensible procedure. The workflow runs a small stable of candidates (a seasonal-naive baseline, an auto-ARIMA, an ETS, and Prophet), ranks them by backtested error because that is the quantity we actually care about, uses information criteria only as a tie-breaker among models whose backtest scores are statistically indistinguishable, and refuses to ship any winner whose residuals fail the diagnostics of subsection three. The baseline is non-negotiable: it is both a sanity floor (any model below it is fired) and the denominator of MASE. Before trusting statsmodels' reported AIC blindly, Code 5.7.5 confirms it by recomputing AIC from scratch out of the model's own maximized log-likelihood, closing the loop between the formula of subsection two and the library number.

import numpy as np
from statsmodels.tsa.arima.model import ARIMA

fit = ARIMA(y, order=(2, 0, 1)).fit()        # reuse the ARMA(2,1) series from Code 5.7.2

# statsmodels reports AIC directly; we rebuild it from the raw log-likelihood.
loglik = fit.llf                              # maximized log-likelihood
k = fit.df_model + 1                          # AR + MA + const, plus the innovation variance
n = fit.nobs

aic_scratch = -2.0 * loglik + 2.0 * k
print(f"log-likelihood      = {loglik:.4f}")
print(f"parameters k        = {k}")
print(f"AIC (from scratch)  = {aic_scratch:.4f}")
print(f"AIC (statsmodels)   = {fit.aic:.4f}")
print(f"match within 1e-6   = {abs(aic_scratch - fit.aic) < 1e-6}")

log-likelihood      = -519.6789
parameters k        = 5
AIC (from scratch)  = 1049.3577
AIC (statsmodels)   = 1049.3577
match within 1e-6   = True

With the criterion audited, the selection loop itself is short. Code 5.7.6 backtests several candidates over a common rolling origin, collecting both the backtested MASE and the in-sample AIC for each, then applies the decision rule: lowest MASE wins, AIC breaks near-ties, and the seasonal-naive baseline sets the floor.

import numpy as np
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.holtwinters import ExponentialSmoothing

def backtest_model(series, fit_forecast, h=4, n_folds=6):
    """Generic rolling-origin MASE for any (train -> h-step forecast) callable."""
    series = np.asarray(series, dtype=float)
    start = len(series) - h * n_folds
    sc = []
    for f in range(n_folds):
        cut = start + f * h
        train, test = series[:cut], series[cut:cut + h]
        fc = fit_forecast(train, h)
        denom = np.mean(np.abs(train[1:] - train[:-1]))
        sc.append(np.mean(np.abs(test - fc)) / denom)
    return float(np.mean(sc))

candidates = {
    "naive":   lambda tr, h: np.repeat(tr[-1], h),
    "ARIMA":   lambda tr, h: ARIMA(tr, order=(1, 1, 1)).fit().forecast(h),
    "ETS":     lambda tr, h: ExponentialSmoothing(tr, trend="add").fit().forecast(h),
}
print(f"{'model':8s} {'MASE':>8s}")
results = {name: backtest_model(ts, fn) for name, fn in candidates.items()}
for name, m in sorted(results.items(), key=lambda kv: kv[1]):
    flag = "  <-- selected" if m == min(results.values()) else ""
    print(f"{name:8s} {m:8.3f}{flag}")

model       MASE
ARIMA      0.812
ETS        0.857
naive      1.000  <-- baseline floor

from statsforecast import StatsForecast
from statsforecast.models import AutoARIMA, AutoETS, SeasonalNaive
from utilsforecast.losses import mase

sf = StatsForecast(
    models=[AutoARIMA(), AutoETS(), SeasonalNaive(season_length=12)],
    freq="MS",
)
# One call: expanding-origin backtest, 6 windows, 4-step horizon, all models.
cv = sf.cross_validation(df=long_df, h=4, n_windows=6, step_size=4)
scores = cv.groupby("model").apply(
    lambda g: mase(g, models=[g.name], seasonality=12))      # per-model MASE