Model Comparison & Selection

Model Comparison & Selection#

Navigation: ← Previous: State-Space Models

We compare all classical methods head-to-head using time-series cross-validation, proper forecast metrics, and the Diebold-Mariano test for statistical significance.

from pathlib import Path
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')
plt.style.use('seaborn-v0_8-whitegrid')

DATA_DIR = Path('data')
if not DATA_DIR.exists():
    DATA_DIR = Path('projects/time-series-analysis/data')

def load_series(name):
    """Load a bundled CSV as a DatetimeIndex Series."""
    df = pd.read_csv(DATA_DIR / name, parse_dates=['date'])
    return df.set_index('date')['value'].sort_index()

def load_air_passengers():
    path = DATA_DIR / 'air_passengers.csv'
    if path.exists():
        return load_series('air_passengers.csv')
    from statsmodels.datasets import get_rdataset
    air = get_rdataset('AirPassengers', 'datasets').data
    s = air['value']
    s.index = pd.date_range('1949-01', periods=len(s), freq='MS')
    return s

def load_milk():
    path = DATA_DIR / 'monthly_milk.csv'
    if path.exists():
        return load_series('monthly_milk.csv')
    url = 'https://raw.githubusercontent.com/jbrownlee/Datasets/master/monthly-milk-production-pounds.csv'
    return pd.read_csv(url, index_col=0, parse_dates=True).squeeze()

def load_co2():
    from statsmodels.datasets import co2
    return co2.load().data.resample('MS').mean().ffill().squeeze()

def load_sunspots():
    from statsmodels.datasets import sunspots
    d = sunspots.load_pandas().data
    return d.set_index('YEAR')['SUNACTIVITY']

def load_retail():
    path = DATA_DIR / 'uk_retail_sales.csv'
    if path.exists():
        return load_series('uk_retail_sales.csv')
    raise FileNotFoundError('uk_retail_sales.csv not found in data/')
from sklearn.model_selection import TimeSeriesSplit
from statsmodels.tsa.holtwinters import ExponentialSmoothing, SimpleExpSmoothing
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.tsa.statespace.structural import UnobservedComponents
from scipy import stats

def mae(y, yhat): return np.mean(np.abs(y - yhat))
def rmse(y, yhat): return np.sqrt(np.mean((y - yhat) ** 2))
def mape(y, yhat): return np.mean(np.abs((y - yhat) / y)) * 100
def smape(y, yhat): return np.mean(2 * np.abs(y - yhat) / (np.abs(y) + np.abs(yhat))) * 100

def crps_normal(y, mu, sigma):
    """CRPS for Gaussian predictive distribution."""
    z = (y - mu) / sigma
    return sigma * (z * (2 * stats.norm.cdf(z) - 1) + 2 * stats.norm.pdf(z) - 1 / np.sqrt(np.pi))

def diebold_mariano(e1, e2, h=1):
    """DM test: H0 equal forecast accuracy."""
    d = e1**2 - e2**2
    d_mean = d.mean()
    n = len(d)
    gamma0 = np.var(d, ddof=1)
    var_d = gamma0 / n
    if var_d <= 0:
        return 0.0, 1.0
    dm_stat = d_mean / np.sqrt(var_d)
    p = 2 * (1 - stats.norm.cdf(np.abs(dm_stat)))
    return dm_stat, p

Time Series Cross-Validation#

Expanding-window CV respects temporal ordering — no random shuffling.

def evaluate_models(series, test_size=12):
    train, test = series.iloc[:-test_size], series.iloc[-test_size:]
    results = {}

    # Holt-Winters multiplicative
    hw = ExponentialSmoothing(
        train, trend='add', seasonal='mul', seasonal_periods=12,
        initialization_method='estimated'
    ).fit(optimized=True)
    fc = hw.forecast(len(test))
    results['HW Multiplicative'] = {'MAE': mae(test, fc), 'RMSE': rmse(test, fc), 'MAPE': mape(test, fc)}

    # SARIMA airline
    sar = SARIMAX(train, order=(0,1,1), seasonal_order=(0,1,1,12)).fit(disp=False)
    fc = sar.forecast(len(test))
    results['SARIMA Airline'] = {'MAE': mae(test, fc), 'RMSE': rmse(test, fc), 'MAPE': mape(test, fc)}

    # Structural
    uc = UnobservedComponents(train, level='local linear trend', seasonal=12).fit(disp=False)
    fc = uc.forecast(len(test))
    results['Structural'] = {'MAE': mae(test, fc), 'RMSE': rmse(test, fc), 'MAPE': mape(test, fc)}

    return pd.DataFrame(results).T

air = load_air_passengers()
milk = load_milk()
air_results = evaluate_models(air)
milk_results = evaluate_models(milk)
print('Air Passengers (12-month holdout):')
display(air_results.style.format('{:.2f}'))
print('Milk Production (12-month holdout):')
display(milk_results.style.format('{:.2f}'))

/opt/hostedtoolcache/Python/3.11.16/x64/lib/python3.11/site-packages/statsmodels/tsa/base/tsa_model.py:480: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/opt/hostedtoolcache/Python/3.11.16/x64/lib/python3.11/site-packages/statsmodels/tsa/base/tsa_model.py:480: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/opt/hostedtoolcache/Python/3.11.16/x64/lib/python3.11/site-packages/statsmodels/tsa/base/tsa_model.py:480: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/opt/hostedtoolcache/Python/3.11.16/x64/lib/python3.11/site-packages/statsmodels/tsa/base/tsa_model.py:480: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)

/opt/hostedtoolcache/Python/3.11.16/x64/lib/python3.11/site-packages/statsmodels/tsa/base/tsa_model.py:480: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/opt/hostedtoolcache/Python/3.11.16/x64/lib/python3.11/site-packages/statsmodels/tsa/base/tsa_model.py:480: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/opt/hostedtoolcache/Python/3.11.16/x64/lib/python3.11/site-packages/statsmodels/tsa/base/tsa_model.py:480: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)

/opt/hostedtoolcache/Python/3.11.16/x64/lib/python3.11/site-packages/statsmodels/tsa/base/tsa_model.py:480: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)

Air Passengers (12-month holdout):
  MAE RMSE MAPE
HW Multiplicative 10.30 15.81 2.21
SARIMA Airline 16.22 21.09 3.65
Structural 78.78 91.29 15.77 
Milk Production (12-month holdout):
  MAE RMSE MAPE
HW Multiplicative 5.24 6.49 0.63
SARIMA Airline 4.78 5.51 0.57
Structural 4.57 5.47 0.55 
# Diebold-Mariano: HW vs SARIMA on Air Passengers
train = air.iloc[:-12]
test = air.iloc[-12:]
hw = ExponentialSmoothing(train, trend='add', seasonal='mul', seasonal_periods=12,
                        initialization_method='estimated').fit(optimized=True)
sar = SARIMAX(train, order=(0,1,1), seasonal_order=(0,1,1,12)).fit(disp=False)
e_hw = test.values - hw.forecast(12).values
e_sar = test.values - sar.forecast(12).values
dm, p = diebold_mariano(e_hw, e_sar)
print(f'Diebold-Mariano (HW vs SARIMA): stat={dm:.3f}, p={p:.4f}')

Diebold-Mariano (HW vs SARIMA): stat=-1.575, p=0.1154

/opt/hostedtoolcache/Python/3.11.16/x64/lib/python3.11/site-packages/statsmodels/tsa/base/tsa_model.py:480: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/opt/hostedtoolcache/Python/3.11.16/x64/lib/python3.11/site-packages/statsmodels/tsa/base/tsa_model.py:480: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)
/opt/hostedtoolcache/Python/3.11.16/x64/lib/python3.11/site-packages/statsmodels/tsa/base/tsa_model.py:480: ValueWarning: No frequency information was provided, so inferred frequency MS will be used.
  self._init_dates(dates, freq)

When to Use Which#

  • SES: No trend/season; maximum simplicity

  • Holt-Winters: Clear trend+season; transparent parameters

  • ARIMA/SARIMA: Complex autocorrelation; formal inference

  • Structural: Interpretable components; missing data; calendar effects

M-competitions (M1–M5) showed exponential smoothing often matches ARIMA; combinations frequently win.

Key Takeaways#

  • Out-of-sample matters: In-sample AIC can mislead; always hold out recent data.

  • No universal winner: Best model depends on DGP, horizon, and series length.

  • DM test: Formal comparison of forecast accuracy; use with caution at small sample sizes.

  • Looking ahead: Neural forecasters (N-BEATS, TFT), Prophet, and gradient boosting extend this toolkit in future projects.

Navigation: ← Previous: State-Space Models

Contents