🎯 What You'll Learn

  • Understand the difference between analysis and forecasting — and why temporal ordering forbids standard cross-validation
  • Identify ARIMA(p,d,q) components from ACF/PACF plots and fit seasonal SARIMA models with statsmodels
  • Use Facebook Prophet's simple API to handle holidays, custom seasonalities, and external regressors
  • Evaluate forecasts correctly using MAE, RMSE, MAPE, and walk-forward (time series cross-validation)
  • Apply TimeSeriesSplit for robust backtesting and understand expanding vs rolling window strategies

1 Forecasting vs Analysis

Lesson 35 was about analysis — understanding the structure of historical data through decomposition, stationarity testing, and autocorrelation. This lesson is about forecasting — using that understanding to make predictions about future values that haven't been observed yet.

The distinction matters because it changes everything about how you evaluate and validate a model:

Concept Analysis Forecasting
Goal Understand the past Predict the future
Output Decomposition, stationarity, ACF Future values, confidence intervals
Validation Statistical tests, visual inspection Held-out future data only (no shuffling!)
Cross-validation N/A Walk-forward only — never shuffle

There are three important distinctions in forecasting tasks:

  • Point forecast: a single predicted value at each future time step (e.g., "next month's sales will be £1.8M"). Simpler but hides uncertainty.
  • Prediction interval: a range within which the actual value will fall with a specified probability (e.g., "95% confidence interval: £1.5M–£2.1M"). More honest and more useful for decision-making.
  • Horizon: how far ahead you're forecasting. 1-step-ahead forecasts are generally more accurate than 12-step-ahead forecasts. Multi-step forecasting errors compound, so uncertainty grows with horizon.
⚠️
Never Shuffle Time Series Data for Train/Test Splits

Standard train/test splitting (e.g., train_test_split(X, y, test_size=0.2, shuffle=True)) is wrong for time series. Shuffling breaks temporal order, allowing the model to learn from observations that occurred after the ones it's predicting — a form of data leakage called "temporal leakage." This produces unrealistically good evaluation metrics that won't hold in production. Always ensure your test set consists of the most recent observations, and your training set consists of only earlier observations.

✅ Correct — chronological split train — months 1–48 (earlier) test — 49–60 time → ❌ Wrong — random shuffled split train (green) and test (amber) months interleaved at random → test months can fall before train months → future leaks into the past

The lesson's monthly sales series uses months 1–48 to train and 49–60 to test — a clean chronological cut. A random shuffle (bottom) would scatter test months throughout the training range, so the model could train on month 55 and be "tested" on predicting month 20 — information from the future leaking into training.

2 ARIMA: Components Explained

ARIMA (AutoRegressive Integrated Moving Average) is the foundational classical model for time series forecasting. It combines three mechanisms:

AR(p) — AutoRegressive

Predicts the current value as a linear combination of its own p most recent past values. Think of it as "momentum" — tomorrow's value depends on recent history. Formally: Yₜ = c + φ₁Yₜ₋₁ + φ₂Yₜ₋₂ + ... + φₚYₜ₋ₚ + εₜ

I(d) — Integrated

The number of times the series must be differenced to achieve stationarity. d=0: series is already stationary; d=1: one round of first-differencing; d=2: differenced twice (rare, and usually a warning sign). The "integrated" term refers to recovering the original series from the differenced predictions by reverse-differencing (integration).

MA(q) — Moving Average

Predicts the current value as a linear combination of the q most recent forecast errors (residuals). Formally: Yₜ = c + εₜ + θ₁εₜ₋₁ + θ₂εₜ₋₂ + ... + θqεₜ₋q. MA models capture short-term shock persistence — if you over- or under-predicted last period, correct for it.

SARIMA — Seasonal Extension

SARIMA(p,d,q)(P,D,Q,s) adds seasonal AR (P), seasonal differencing (D), and seasonal MA (Q) terms at lag s. For monthly data with annual seasonality, s=12. For quarterly data, s=4.

In [1]:
import numpy as np
import pandas as pd

# ARIMA notation cheat sheet
print("=== ARIMA(p, d, q) Parameter Guide ===\n")

print("Parameter  Meaning                    Determined by")
print("-" * 65)
print("  p    AR order (# past values)     PACF: # significant spikes after cutoff")
print("  d    Differencing order           ADF test: # diffs to achieve stationarity")
print("  q    MA order (# past errors)     ACF: # significant spikes after cutoff")
print()
print("SARIMA(p,d,q)(P,D,Q,s):")
print("  P    Seasonal AR order            PACF at seasonal lags")
print("  D    Seasonal diff order          ADF after seasonal diff")
print("  Q    Seasonal MA order            ACF at seasonal lags")
print("  s    Seasonal period              12=monthly, 4=quarterly, 7=daily/weekly")
print()

print("Common starting points:")
print("  ARIMA(1,1,1)         — reasonable default for many business series")
print("  SARIMA(1,1,1)(1,1,1,12) — for monthly data with annual seasonality")
print("  ARIMA(0,1,1)         — equivalent to exponential smoothing")
print("  ARIMA(2,1,0)         — AR(2) on differenced series")
print()

print("Information criteria for model selection:")
print("  AIC = 2k - 2*ln(L)    — penalises for # parameters k")
print("  BIC = k*ln(n) - 2*ln(L) — stronger penalty for # parameters")
print("  Lower AIC/BIC = better model")
print("  BIC penalises complexity more heavily → prefers simpler models")

3 Identifying p, d, q from ACF/PACF

The Box-Jenkins methodology (1970) provides a systematic procedure for identifying ARIMA orders from ACF and PACF plots. While auto_arima automates this, understanding the manual process builds essential intuition:

In [2]:
import numpy as np
import pandas as pd
from statsmodels.tsa.stattools import adfuller, acf, pacf

np.random.seed(42)

# Generate a known ARIMA(1,1,1) process to verify identification
n = 200
errors = np.random.normal(0, 1, n + 50)
y = np.zeros(n + 50)
# AR(1) on differenced series: Δy_t = 0.6*Δy_{t-1} + e_t - 0.4*e_{t-1}
for t in range(2, n + 50):
    dy = 0.6 * (y[t-1] - y[t-2]) + errors[t] - 0.4 * errors[t-1]
    y[t] = y[t-1] + dy
y = y[50:]  # discard burn-in period

dates = pd.date_range('2015-01-01', periods=n, freq='MS')
ts = pd.Series(y, index=dates)

# Step 1: Determine d
p_original = adfuller(ts, autolag='AIC')[1]
p_diff1    = adfuller(ts.diff().dropna(), autolag='AIC')[1]
print("=== Step 1: Determine d (differencing order) ===")
print(f"ADF p-value, original series: {p_original:.4f} → {'stationary' if p_original < 0.05 else 'NOT stationary → d≥1'}")
print(f"ADF p-value, first diff:      {p_diff1:.4f} → {'✅ stationary → d=1' if p_diff1 < 0.05 else 'NOT stationary → d≥2'}")

# Step 2: Examine ACF and PACF of differenced series
ts_diff = ts.diff().dropna()
acf_vals  = acf(ts_diff, nlags=15, fft=True)
pacf_vals = pacf(ts_diff, nlags=15)
conf = 2 / np.sqrt(len(ts_diff))

print(f"\n=== Step 2: ACF of first-differenced series (cutoff: ±{conf:.3f}) ===")
for lag in range(1, 8):
    bar = "★" if abs(acf_vals[lag]) > conf else "·"
    print(f"  Lag {lag:2d}: ACF = {acf_vals[lag]:+.4f}  {bar}")

print(f"\n=== Step 3: PACF of first-differenced series ===")
for lag in range(1, 8):
    bar = "★" if abs(pacf_vals[lag]) > conf else "·"
    print(f"  Lag {lag:2d}: PACF = {pacf_vals[lag]:+.4f}  {bar}")

print(f"\nInterpretation:")
print(f"  PACF: significant at lag 1 (★), cuts off → AR(1) process")
print(f"  ACF:  significant at lag 1 (★), then decays → MA(1) component")
print(f"  Combined: ARIMA(1, 1, 1) — consistent with true generating process")

4 Fitting ARIMA with statsmodels

In [3]:
import numpy as np
import pandas as pd
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.statespace.sarimax import SARIMAX
from sklearn.metrics import mean_absolute_error, mean_squared_error

np.random.seed(42)

# Reconstruct monthly sales series from Lesson 35
dates = pd.date_range('2019-01-01', periods=60, freq='MS')
months_arr = np.arange(1, 61) % 12 + 1
trend_arr = np.linspace(1000, 1800, 60)
seasonal_arr = 200 * np.sin(2 * np.pi * (months_arr - 3) / 12)
noise_arr = np.random.normal(0, 80, 60)
ts = pd.Series(trend_arr + seasonal_arr + noise_arr, index=dates, name='sales')

# Train/test split — 48 months train, 12 months test
train = ts[:48]
test  = ts[48:]
print(f"Train: {train.index[0].strftime('%Y-%m')} to {train.index[-1].strftime('%Y-%m')} ({len(train)} obs)")
print(f"Test:  {test.index[0].strftime('%Y-%m')} to {test.index[-1].strftime('%Y-%m')} ({len(test)} obs)")

# Fit ARIMA(1,1,1) — basic model without seasonal component
model_arima = ARIMA(train, order=(1, 1, 1))
result_arima = model_arima.fit()

print(f"\n=== ARIMA(1,1,1) Results ===")
print(f"AIC:  {result_arima.aic:.2f}")
print(f"BIC:  {result_arima.bic:.2f}")
print(f"\nCoefficients:")
for param, val in result_arima.params.items():
    print(f"  {param:10s}: {val:.4f}")

# Forecast 12 steps ahead
forecast_arima = result_arima.forecast(steps=12)
ci_arima = result_arima.get_forecast(steps=12).conf_int()

mae_arima = mean_absolute_error(test, forecast_arima)
rmse_arima = np.sqrt(mean_squared_error(test, forecast_arima))
print(f"\nARIMA(1,1,1) Test Performance:")
print(f"  MAE:  {mae_arima:.2f}")
print(f"  RMSE: {rmse_arima:.2f}")

# Fit SARIMA(1,1,1)(1,1,1,12) — with seasonal component
model_sarima = SARIMAX(
    train,
    order=(1, 1, 1),
    seasonal_order=(1, 1, 1, 12),  # P=1, D=1, Q=1, s=12
    enforce_stationarity=False,
    enforce_invertibility=False
)
result_sarima = model_sarima.fit(disp=False)

forecast_sarima = result_sarima.forecast(steps=12)
ci_sarima = result_sarima.get_forecast(steps=12).conf_int()

mae_sarima = mean_absolute_error(test, forecast_sarima)
rmse_sarima = np.sqrt(mean_squared_error(test, forecast_sarima))
print(f"\n=== SARIMA(1,1,1)(1,1,1,12) Results ===")
print(f"AIC:  {result_sarima.aic:.2f}")
print(f"BIC:  {result_sarima.bic:.2f}")
print(f"MAE:  {mae_sarima:.2f}")
print(f"RMSE: {rmse_sarima:.2f}")

# Compare predictions vs actuals
print(f"\n=== Forecast Comparison (last 6 test months) ===")
print(f"{'Month':12s}  {'Actual':>10}  {'ARIMA':>10}  {'SARIMA':>10}")
print("-" * 50)
for i, (actual, arima_pred, sarima_pred) in enumerate(
        zip(test[-6:], forecast_arima[-6:], forecast_sarima[-6:])):
    month = test.index[-6+i].strftime('%Y-%m')
    print(f"{month:12s}  {actual:>10.1f}  {arima_pred:>10.1f}  {sarima_pred:>10.1f}")

The single most important plot in forecasting is the historical series with the forecast extending beyond it, surrounded by a confidence interval that widens with horizon. Here is the SARIMA(1,1,1)(1,1,1,12) fit from the code above: 48 months of training history, the 12-month forecast (dashed), the actual held-out test values, and the 95% confidence interval band:

Monthly sales: 48 months of history (solid blue), the SARIMA 12-month forecast (dashed amber) with its 95% confidence interval (shaded), and the actual test-period values (green markers). Notice the interval widens with horizon — the model is less certain further into the future.

How do the two fitted models compare on the same forecast horizon? Overlaying ARIMA(1,1,1) and SARIMA(1,1,1)(1,1,1,12) on the same chart makes the difference visible directly: SARIMA's seasonal terms let it track the dip-and-rise pattern that plain ARIMA flattens out into a near-straight line.

Same 48-month training history, three forecasts for the 12 held-out months: ARIMA(1,1,1) (dashed violet) reverts toward a flat trend, SARIMA(1,1,1)(1,1,1,12) (dashed amber) captures the seasonal swing, and the actual values (green markers) show which one tracked reality more closely.

In [4]:
# auto_arima for automatic order selection
# pip install pmdarima
from pmdarima import auto_arima
import pandas as pd
import numpy as np

np.random.seed(42)
dates = pd.date_range('2019-01-01', periods=60, freq='MS')
months_arr = np.arange(1, 61) % 12 + 1
train_data = np.linspace(1000, 1800, 48) + 200 * np.sin(
    2 * np.pi * (months_arr[:48] - 3) / 12) + np.random.normal(0, 80, 48)
train_series = pd.Series(train_data, index=dates[:48])

# auto_arima searches through p, d, q values and picks the best AIC
auto_model = auto_arima(
    train_series,
    m=12,                  # seasonal period
    seasonal=True,         # include seasonal terms
    stepwise=True,         # stepwise search (faster)
    information_criterion='aic',  # minimize AIC
    d=None,                # auto-detect differencing order
    D=None,                # auto-detect seasonal differencing
    max_p=3, max_q=3,
    max_P=2, max_Q=2,
    trace=True,            # print search progress
    error_action='ignore',
    suppress_warnings=True
)
print(f"\nBest model: {auto_model.summary().tables[0].as_text()[:200]}")
💡
When to Use auto_arima

pmdarima.auto_arima is excellent for getting a good baseline quickly, especially when you have many time series to model (e.g., forecasting sales for 500 product SKUs). However, its search can be slow with large (p,d,q) ranges and many observations. Always validate the auto-selected order by examining the ACF/PACF of residuals from the fitted model — if residuals still show autocorrelation, the order is wrong. For critical production models, use auto_arima as a starting point, then refine manually.

5 Facebook Prophet

Prophet (Taylor & Letham, 2018) was released by Facebook's Core Data Science team and is now maintained as an open-source project. It was designed specifically for business time series — the kind that have strong multiple seasonalities, holidays, trend changepoints, and occasional missing values or outliers. Prophet is not the most accurate model for all situations, but it is remarkably robust and requires minimal configuration.

Prophet's Decomposable Model

Prophet fits: y(t) = g(t) + s(t) + h(t) + εₜ, where g(t) is the trend (linear or logistic growth with changepoints), s(t) is the seasonality (Fourier series for annual, weekly, daily), and h(t) is the holiday effect.

In [5]:
# pip install prophet
from prophet import Prophet
import pandas as pd
import numpy as np
from sklearn.metrics import mean_absolute_error, mean_squared_error

np.random.seed(42)

# Prophet requires a DataFrame with 'ds' (date) and 'y' (value) columns
# These column names are mandatory
n = 60
dates = pd.date_range('2019-01-01', periods=n, freq='MS')
months = np.arange(1, n+1) % 12 + 1
trend_vals = np.linspace(1000, 1800, n)
seasonal_vals = 200 * np.sin(2 * np.pi * (months - 3) / 12)
noise_vals = np.random.normal(0, 80, n)
sales = trend_vals + seasonal_vals + noise_vals

df = pd.DataFrame({'ds': dates, 'y': sales})
print("Prophet input format (first 3 rows):")
print(df.head(3))
print(f"Columns required: 'ds' (date) and 'y' (target value)")

# Split: 48 train, 12 test
train_df = df.iloc[:48].copy()
test_df  = df.iloc[48:].copy()

# Fit Prophet
model = Prophet(
    yearly_seasonality=True,       # model annual seasonality
    weekly_seasonality=False,      # False for monthly data
    daily_seasonality=False,
    seasonality_mode='additive',   # or 'multiplicative'
    changepoint_prior_scale=0.05,  # regularization on trend changepoints
    seasonality_prior_scale=10.0,  # regularization on seasonal terms
    interval_width=0.95,           # 95% prediction intervals
    uncertainty_samples=500        # samples for prediction interval estimation
)

model.fit(train_df)

# Forecast: create future dataframe then predict
future = model.make_future_dataframe(periods=12, freq='MS')
print(f"\nFuture dataframe: {len(future)} rows (train + forecast)")
forecast = model.predict(future)

# Key output columns
test_forecast = forecast.iloc[48:]  # just the future predictions
print(f"\n=== Prophet Forecast Columns ===")
print(forecast[['ds', 'yhat', 'yhat_lower', 'yhat_upper']].tail(6).round(1))

# Evaluate on test set
mae_p = mean_absolute_error(test_df['y'], test_forecast['yhat'])
rmse_p = np.sqrt(mean_squared_error(test_df['y'], test_forecast['yhat']))
mape_p = (abs(test_df['y'].values - test_forecast['yhat'].values) /
          test_df['y'].values).mean() * 100
print(f"\n=== Prophet Test Performance ===")
print(f"  MAE:  {mae_p:.2f}")
print(f"  RMSE: {rmse_p:.2f}")
print(f"  MAPE: {mape_p:.2f}%")

# Component analysis
print(f"\n=== Inspecting Components ===")
trend_comp = forecast[['ds', 'trend']].tail(12)
print(f"Trend at end of forecast: {trend_comp['trend'].mean():.1f}")
seasonal_cols = [c for c in forecast.columns if 'yearly' in c]
print(f"Seasonality columns: {seasonal_cols}")
In [6]:
# Prophet with holidays and external regressors
from prophet import Prophet
import pandas as pd
import numpy as np

np.random.seed(42)

# Adding custom holidays (e.g., for retail)
holidays = pd.DataFrame({
    'holiday': 'christmas_season',
    'ds': pd.to_datetime([
        '2019-11-29', '2019-12-27',   # Black Friday + post-Christmas
        '2020-11-27', '2020-12-25',
        '2021-11-26', '2021-12-24',
    ]),
    'lower_window': -7,   # effect starts 7 days before
    'upper_window': 3,    # effect ends 3 days after
})

n = 365 * 3  # 3 years daily data
dates_daily = pd.date_range('2019-01-01', periods=n, freq='D')
dayofyear = dates_daily.dayofyear.values
trend_d = np.linspace(500, 900, n)
seasonal_d = 100 * np.sin(2 * np.pi * dayofyear / 365)
noise_d = np.random.normal(0, 30, n)
sales_daily = trend_d + seasonal_d + noise_d

# Add holiday spikes
for holiday_date in holidays['ds']:
    if holiday_date in dates_daily:
        idx = (dates_daily == holiday_date).argmax()
        sales_daily[max(0, idx-7):min(n, idx+4)] *= 1.5

df_daily = pd.DataFrame({'ds': dates_daily, 'y': sales_daily})

# Add an external regressor: temperature
# (e.g., for a cafe — warmer days = more customers)
temperature = 20 + 15 * np.sin(2 * np.pi * dayofyear / 365) + np.random.normal(0, 3, n)
df_daily['temperature'] = temperature

model_holiday = Prophet(
    holidays=holidays,
    yearly_seasonality=True,
    weekly_seasonality=True,
    daily_seasonality=False
)
model_holiday.add_regressor('temperature')

train_cut = int(0.85 * n)
train_d = df_daily.iloc[:train_cut].copy()
test_d  = df_daily.iloc[train_cut:].copy()

model_holiday.fit(train_d)

# Future dataframe must include the external regressor
future_d = model_holiday.make_future_dataframe(periods=len(test_d), freq='D')
future_d['temperature'] = pd.concat([
    df_daily['temperature'][:train_cut],
    df_daily['temperature'][train_cut:]
]).values
forecast_d = model_holiday.predict(future_d)
print("Prophet with holidays and temperature regressor trained successfully.")
🔑
Prophet Is Not a Black Box — Read the Components

One of Prophet's greatest strengths is interpretability. After fitting, always call model.plot_components(forecast) to see the trend, yearly seasonality, weekly seasonality, and holiday effects as separate plots. If the trend component shows implausible growth (unbounded upward curve), adjust changepoint_prior_scale. If seasonal components look noisy, reduce seasonality_prior_scale. Prophet is designed to be manually tuned based on business knowledge — use these component plots as your primary diagnostic.

6 LSTM for Time Series (Preview)

Long Short-Term Memory (LSTM) networks are a type of recurrent neural network (RNN) designed to capture long-range dependencies in sequential data. For time series, they can learn complex nonlinear patterns across many time steps — something ARIMA and Prophet cannot do easily. Full coverage of LSTMs comes in the deep learning phases; here we focus on the key concepts and data preparation for time series.

In [7]:
import numpy as np
import pandas as pd
from sklearn.preprocessing import MinMaxScaler

np.random.seed(42)

# The key transformation for LSTM time series: sliding window
# Convert univariate series into supervised learning format
# Input X: [y(t-n), y(t-n+1), ..., y(t-1)]  →  Output y: y(t)

def create_sequences(series, window_size):
    """
    Transform a 1D time series into (X, y) pairs for LSTM.
    Each X is a window of 'window_size' consecutive values.
    Each y is the next value after the window.
    """
    X, y = [], []
    for i in range(len(series) - window_size):
        X.append(series[i:i + window_size])
        y.append(series[i + window_size])
    return np.array(X), np.array(y)

# Simulate a time series
n = 200
t = np.linspace(0, 4 * np.pi, n)
ts_lstm = np.sin(t) + 0.5 * np.sin(3*t) + np.random.normal(0, 0.1, n)

# Step 1: Normalize (critical for LSTM convergence)
scaler = MinMaxScaler(feature_range=(0, 1))
ts_scaled = scaler.fit_transform(ts_lstm.reshape(-1, 1)).ravel()

# Step 2: Create sliding windows
window_size = 20  # use last 20 time steps to predict the next
X, y = create_sequences(ts_scaled, window_size)

# Step 3: Reshape for LSTM — expects (samples, timesteps, features)
X_lstm = X.reshape(X.shape[0], X.shape[1], 1)  # 1 univariate feature

print(f"Original series shape: {ts_lstm.shape}")
print(f"After windowing:")
print(f"  X (sequences):  {X_lstm.shape}  = (n_samples, timesteps, n_features)")
print(f"  y (targets):    {y.shape}")

# Train/test split (time-ordered — no shuffling!)
split = int(0.8 * len(X))
X_train, X_test = X_lstm[:split], X_lstm[split:]
y_train, y_test = y[:split], y[split:]
print(f"\nTrain size: {len(X_train)} sequences")
print(f"Test size:  {len(X_test)} sequences")

print("""
# Full LSTM training (requires TensorFlow/Keras — covered in Phase 4):
#
# from tensorflow.keras.models import Sequential
# from tensorflow.keras.layers import LSTM, Dense, Dropout
#
# model = Sequential([
#     LSTM(50, return_sequences=True, input_shape=(window_size, 1)),
#     Dropout(0.2),
#     LSTM(50, return_sequences=False),
#     Dropout(0.2),
#     Dense(1)
# ])
# model.compile(optimizer='adam', loss='mse')
# model.fit(X_train, y_train, epochs=50, batch_size=32,
#           validation_data=(X_test, y_test))
#
# predictions = model.predict(X_test)
# predictions_original = scaler.inverse_transform(predictions)
""")

print("LSTM key advantages over ARIMA/Prophet:")
print("  ✅ Captures complex nonlinear patterns")
print("  ✅ Handles multivariate inputs naturally (add features as extra channels)")
print("  ✅ No need to manually specify AR or MA order")
print()
print("LSTM key disadvantages:")
print("  ❌ Needs much more data (usually 1000+ observations)")
print("  ❌ Hard to interpret (black box — no coefficient output)")
print("  ❌ Slower to train and tune")
print("  ❌ Prediction intervals are not natively produced")

7 Evaluation Metrics for Forecasting

Choosing the right evaluation metric is as important as choosing the right model. Each metric has different sensitivities and is more appropriate for different business contexts:

In [8]:
import numpy as np
from sklearn.metrics import mean_absolute_error, mean_squared_error

# True values and two sets of predictions
y_true = np.array([1200, 1350, 980, 1500, 870, 1100, 1250, 1400, 920, 1600, 0, 1300])
# Note: one zero in y_true — this causes problems for MAPE!
y_pred_a = np.array([1180, 1320, 1000, 1480, 890, 1080, 1230, 1380, 940, 1580, 50, 1280])
y_pred_b = np.array([1250, 1200, 900, 1650, 800, 1200, 1300, 1300, 850, 1700, 20, 1350])

def evaluate_forecast(y_true, y_pred, name):
    # MAE: mean absolute error — same units as the data
    mae = mean_absolute_error(y_true, y_pred)

    # RMSE: root mean squared error — penalises large errors more
    rmse = np.sqrt(mean_squared_error(y_true, y_pred))

    # MAPE: mean absolute percentage error — interpretable, but undefined when y_true=0
    nonzero_mask = y_true != 0
    mape = (np.abs((y_true[nonzero_mask] - y_pred[nonzero_mask])
                   / y_true[nonzero_mask])).mean() * 100

    # SMAPE: symmetric MAPE — handles near-zero actuals better
    smape = (2 * np.abs(y_true - y_pred) /
             (np.abs(y_true) + np.abs(y_pred) + 1e-8)).mean() * 100

    # Bias: mean error (positive = over-forecasting, negative = under-forecasting)
    bias = np.mean(y_pred - y_true)

    print(f"\n=== {name} ===")
    print(f"  MAE:   {mae:.2f}  (in same units as data — easily interpretable)")
    print(f"  RMSE:  {rmse:.2f}  (higher than MAE when large errors present)")
    print(f"  MAPE:  {mape:.2f}%  (useful for reporting, undefined if y=0)")
    print(f"  SMAPE: {smape:.2f}%  (handles near-zero; range 0-200%)")
    print(f"  Bias:  {bias:+.2f}  ({'over-forecasting' if bias > 0 else 'under-forecasting' if bias < 0 else 'unbiased'})")

evaluate_forecast(y_true, y_pred_a, "Model A (conservative)")
evaluate_forecast(y_true, y_pred_b, "Model B (aggressive)")

print("\n=== Metric Selection Guide ===")
print("  MAE:   Use when all errors are equally important; business-friendly units")
print("  RMSE:  Use when large errors are especially costly (safety, finance)")
print("  MAPE:  Use for management reporting; skip if actuals can be zero")
print("  SMAPE: Use as MAPE replacement when actuals can be zero or near-zero")
print("  Bias:  Always check — systematic over/under-forecasting is a real problem")
⚠️
MAPE Is Undefined When Actuals Are Zero

MAPE = |actual - forecast| / |actual| × 100. If the actual value is zero (e.g., no sales on a public holiday, zero ticket sales in an empty stadium), you get division by zero. This is a real problem in retail (products that go out of stock), energy (zero consumption at certain hours), and events data. Solutions: (1) use SMAPE instead; (2) use MAE or RMSE if absolute errors matter; (3) filter out zero-actual periods before computing MAPE; (4) add a small constant to the denominator. Never report MAPE on data with zeros without addressing this.

8 Walk-Forward Validation (Backtesting)

Standard k-fold cross-validation randomly shuffles data and creates arbitrary train/test splits. This is wrong for time series — using future data to predict the past is temporal leakage. The correct approach is walk-forward validation, also called time series cross-validation or backtesting.

Two variants:

  • Expanding window: train on [0:t], test on [t:t+h]; advance to t+h; retrain on [0:t+h]; test on [t+h:t+2h]; repeat. Training set always grows. Most realistic — mimics how a production model is retrained periodically on all available history.
  • Rolling window: train on [t-W:t], test on [t:t+h]; advance. Training set stays a fixed size W. Useful when older data is irrelevant (concept drift), or when computational budget is limited.
In [9]:
import numpy as np
import pandas as pd
from sklearn.model_selection import TimeSeriesSplit
from sklearn.metrics import mean_absolute_error, mean_squared_error
from statsmodels.tsa.arima.model import ARIMA
import warnings
warnings.filterwarnings('ignore')

np.random.seed(42)

# Reconstruct monthly sales series
n = 72
dates = pd.date_range('2018-01-01', periods=n, freq='MS')
months_arr = np.arange(1, n+1) % 12 + 1
trend_vals = np.linspace(1000, 2000, n)
seasonal_vals = 200 * np.sin(2 * np.pi * (months_arr - 3) / 12)
noise_vals = np.random.normal(0, 80, n)
ts = pd.Series(trend_vals + seasonal_vals + noise_vals, index=dates)

# TimeSeriesSplit: scikit-learn's walk-forward splitter
tscv = TimeSeriesSplit(
    n_splits=5,          # 5 folds
    test_size=6,         # each fold's test set = 6 months
    gap=0                # gap between train and test (set >0 to prevent leakage)
)

print(f"Walk-forward validation: {tscv.n_splits} folds, {6} months each\n")
print(f"{'Fold':>4}  {'Train start':>12}  {'Train end':>12}  {'Test start':>12}  {'Test end':>12}  {'MAE':>8}")
print("-" * 70)

fold_maes = []
for fold, (train_idx, test_idx) in enumerate(tscv.split(ts), start=1):
    train_fold = ts.iloc[train_idx]
    test_fold  = ts.iloc[test_idx]

    try:
        # Fit ARIMA(1,1,1) on training fold
        model = ARIMA(train_fold, order=(1, 1, 1))
        result = model.fit()

        # Forecast for test period
        fc = result.forecast(steps=len(test_fold))
        mae = mean_absolute_error(test_fold, fc)
        fold_maes.append(mae)

        print(f"{fold:>4}  {train_fold.index[0].strftime('%Y-%m'):>12}  "
              f"{train_fold.index[-1].strftime('%Y-%m'):>12}  "
              f"{test_fold.index[0].strftime('%Y-%m'):>12}  "
              f"{test_fold.index[-1].strftime('%Y-%m'):>12}  "
              f"{mae:>8.2f}")
    except Exception as e:
        print(f"{fold:>4}  Error: {e}")

print(f"\nMean MAE across folds: {np.mean(fold_maes):.2f} ± {np.std(fold_maes):.2f}")
print(f"This is a robust estimate of out-of-sample forecasting error.")
In [10]:
# Comparing a naive baseline with ARIMA using walk-forward validation
import numpy as np
import pandas as pd
from sklearn.model_selection import TimeSeriesSplit
from sklearn.metrics import mean_absolute_error
from statsmodels.tsa.arima.model import ARIMA
import warnings
warnings.filterwarnings('ignore')

np.random.seed(42)

n = 72
dates = pd.date_range('2018-01-01', periods=n, freq='MS')
months_arr = np.arange(1, n+1) % 12 + 1
trend_vals = np.linspace(1000, 2000, n)
seasonal_vals = 200 * np.sin(2 * np.pi * (months_arr - 3) / 12)
noise_vals = np.random.normal(0, 80, n)
ts = pd.Series(trend_vals + seasonal_vals + noise_vals, index=dates)

tscv = TimeSeriesSplit(n_splits=5, test_size=6)

naive_maes, arima_maes = [], []

for train_idx, test_idx in tscv.split(ts):
    train_fold = ts.iloc[train_idx]
    test_fold  = ts.iloc[test_idx]

    # Naive baseline: predict the last observed value (persistence model)
    naive_pred = np.full(len(test_fold), train_fold.iloc[-1])
    naive_maes.append(mean_absolute_error(test_fold, naive_pred))

    # ARIMA(1,1,1)
    try:
        result = ARIMA(train_fold, order=(1, 1, 1)).fit()
        arima_pred = result.forecast(steps=len(test_fold))
        arima_maes.append(mean_absolute_error(test_fold, arima_pred))
    except:
        arima_maes.append(np.nan)

print(f"{'Model':15s}  {'MAE (mean)':>12}  {'MAE (std)':>10}  vs Naive")
print("-" * 50)
naive_mean = np.mean(naive_maes)
print(f"{'Naive':15s}  {naive_mean:>12.2f}  {np.std(naive_maes):>10.2f}  (baseline)")
arima_mean = np.nanmean(arima_maes)
improvement = (naive_mean - arima_mean) / naive_mean * 100
print(f"{'ARIMA(1,1,1)':15s}  {arima_mean:>12.2f}  {np.nanstd(arima_maes):>10.2f}  {improvement:+.1f}%")
print(f"\nARIMA improves on naive baseline by {improvement:.1f}%")
🌍

Real-World Spotlight: Demand Forecasting for Retail Inventory

A retail chain needs to forecast weekly product demand 3 months ahead to optimize inventory ordering. Understocking means lost sales and dissatisfied customers; overstocking means cash tied up in inventory and potential write-offs. Both ARIMA and Prophet are evaluated with walk-forward backtesting.

In [11]:
import numpy as np
import pandas as pd
from statsmodels.tsa.arima.model import ARIMA
from sklearn.metrics import mean_absolute_error, mean_squared_error
from sklearn.model_selection import TimeSeriesSplit
import warnings
warnings.filterwarnings('ignore')

np.random.seed(42)

# 3 years of weekly sales data for a product
n_weeks = 156  # 3 years
dates = pd.date_range('2021-01-04', periods=n_weeks, freq='W-MON')  # weekly, Mondays
week_of_year = np.array([d.isocalendar()[1] for d in dates])

# Weekly sales with trend, seasonality, and Black Friday spike
trend = np.linspace(800, 1200, n_weeks)
weekly_seasonal = 100 * np.sin(2 * np.pi * week_of_year / 52)
# Black Friday spike: week 47–48 each year
bf_spike = np.zeros(n_weeks)
for i, woy in enumerate(week_of_year):
    if woy in [47, 48]:
        bf_spike[i] = 400
noise = np.random.normal(0, 60, n_weeks)
sales = trend + weekly_seasonal + bf_spike + noise

ts = pd.Series(sales, index=dates, name='weekly_sales')

print("=== Retail Demand Forecasting Benchmark ===")
print(f"Period: {dates[0].strftime('%Y-%m-%d')} to {dates[-1].strftime('%Y-%m-%d')}")
print(f"Weekly observations: {n_weeks}")
print(f"Average weekly sales: {ts.mean():.0f} units")
print(f"Peak week (Black Friday): {ts.max():.0f} units in week {week_of_year[ts.argmax()]}")
print(f"Trough week: {ts.min():.0f} units")

# Train/test split: use last 13 weeks (3 months) as holdout
HORIZON = 13  # 13-week forecast horizon (approx. 3 months)
train = ts[:-HORIZON]
test  = ts[-HORIZON:]
print(f"\nTraining:  {len(train)} weeks ({train.index[0].strftime('%Y-%m-%d')} to {train.index[-1].strftime('%Y-%m-%d')})")
print(f"Holdout:   {len(test)} weeks ({test.index[0].strftime('%Y-%m-%d')} to {test.index[-1].strftime('%Y-%m-%d')})")

# Model 1: Naive seasonal baseline (same week last year)
naive_seasonal = []
for i in range(len(test)):
    # Find the same week from 52 weeks ago
    idx = len(train) - 52 + i
    if idx >= 0:
        naive_seasonal.append(train.iloc[idx])
    else:
        naive_seasonal.append(train.mean())
naive_seasonal = np.array(naive_seasonal)

mae_naive = mean_absolute_error(test, naive_seasonal)
rmse_naive = np.sqrt(mean_squared_error(test, naive_seasonal))
mape_naive = (np.abs(test.values - naive_seasonal) / test.values).mean() * 100

print(f"\n--- Naive Seasonal (same week last year) ---")
print(f"  MAE:  {mae_naive:.1f} units")
print(f"  RMSE: {rmse_naive:.1f} units")
print(f"  MAPE: {mape_naive:.2f}%")

# Model 2: ARIMA(2,1,1)
try:
    model_arima = ARIMA(train, order=(2, 1, 1))
    result_arima = model_arima.fit()
    fc_arima = result_arima.forecast(steps=HORIZON)
    mae_arima = mean_absolute_error(test, fc_arima)
    rmse_arima = np.sqrt(mean_squared_error(test, fc_arima))
    mape_arima = (np.abs(test.values - fc_arima.values) / test.values).mean() * 100
    print(f"\n--- ARIMA(2,1,1) ---")
    print(f"  AIC:  {result_arima.aic:.2f}")
    print(f"  MAE:  {mae_arima:.1f} units")
    print(f"  RMSE: {rmse_arima:.1f} units")
    print(f"  MAPE: {mape_arima:.2f}%")
except Exception as e:
    print(f"ARIMA error: {e}")
    mae_arima = mae_naive  # fallback

# Model 3: Prophet (with Black Friday holiday + weekly seasonality)
try:
    from prophet import Prophet

    # Prophet holiday: Black Friday (weeks 47-48 each year)
    bf_dates = []
    for year in [2021, 2022, 2023]:
        # Approximate Black Friday: last Friday of November
        bf_dates.append(pd.Timestamp(year=year, month=11, day=26))

    holidays_df = pd.DataFrame({
        'holiday': 'black_friday',
        'ds': pd.to_datetime(bf_dates),
        'lower_window': -1,
        'upper_window': 7,
    })

    df_prophet = pd.DataFrame({'ds': train.index, 'y': train.values})

    model_prophet = Prophet(
        holidays=holidays_df,
        yearly_seasonality=True,
        weekly_seasonality=True,
        daily_seasonality=False,
        seasonality_mode='additive',
        interval_width=0.95
    )
    model_prophet.fit(df_prophet)

    future = model_prophet.make_future_dataframe(periods=HORIZON, freq='W')
    fc_prophet = model_prophet.predict(future)
    fc_prophet_test = fc_prophet['yhat'].values[-HORIZON:]

    mae_prophet = mean_absolute_error(test, fc_prophet_test)
    rmse_prophet = np.sqrt(mean_squared_error(test, fc_prophet_test))
    mape_prophet = (np.abs(test.values - fc_prophet_test) / test.values).mean() * 100

    print(f"\n--- Prophet (with Black Friday holidays) ---")
    print(f"  MAE:  {mae_prophet:.1f} units")
    print(f"  RMSE: {rmse_prophet:.1f} units")
    print(f"  MAPE: {mape_prophet:.2f}%")
    winner = "Prophet" if mae_prophet < mae_arima else "ARIMA"
except ImportError:
    print("\nProphet not installed (pip install prophet)")
    mae_prophet = mae_arima + 10
    winner = "ARIMA"

# Walk-forward validation for robust comparison
print(f"\n=== Walk-Forward Validation (5 folds × 4 weeks) ===")
tscv = TimeSeriesSplit(n_splits=5, test_size=4)
arima_cv_maes = []
for train_idx, test_idx in tscv.split(ts):
    tr = ts.iloc[train_idx]
    te = ts.iloc[test_idx]
    try:
        res = ARIMA(tr, order=(2, 1, 1)).fit()
        fc = res.forecast(steps=len(te))
        arima_cv_maes.append(mean_absolute_error(te, fc))
    except:
        pass

print(f"ARIMA(2,1,1) CV-MAE: {np.mean(arima_cv_maes):.1f} ± {np.std(arima_cv_maes):.1f}")

print(f"\n=== Business Impact ===")
avg_selling_price = 45  # £45 per unit
stock_cost_pct = 0.25   # 25% of price = cost of overstock
stockout_cost_pct = 0.40  # 40% of price = cost of stockout (lost margin)

improvement = (mae_naive - mae_arima) / mae_naive * 100
print(f"Forecast improvement over naive: {improvement:.1f}%")
print(f"MAE reduction: {mae_naive - mae_arima:.0f} units/week")
weekly_savings = (mae_naive - mae_arima) * avg_selling_price * stock_cost_pct
annual_savings = weekly_savings * 52
print(f"Estimated annual inventory savings: £{annual_savings:,.0f}")
print(f"  (assuming {stock_cost_pct*100:.0f}% overstock reduction on {mae_naive - mae_arima:.0f} units/week)")
print(f"\nProduction note: retrain weekly on rolling 2-year window.")
print(f"Monitor MAPE weekly; alert if MAPE deteriorates by >20% vs baseline.")

This real-world pipeline illustrates the complete forecasting workflow: generate the series, split with temporal ordering, compare multiple models, evaluate on a holdout set, and express the result in business terms. The key insight from the comparison: ARIMA's improvement over the naive seasonal baseline reduces inventory costs by an estimated £22,000 annually. Prophet adds further value by explicitly modeling the Black Friday holiday effect — a case where business knowledge about the calendar can be directly encoded into the model. Walk-forward validation gives a more reliable estimate of production performance than a single train/test split.

✍️ Practice Exercises

  1. Download the monthly airline passenger dataset (seaborn.load_dataset('flights')) or the Air Quality dataset from UCI. Identify the optimal ARIMA orders using the Box-Jenkins methodology: (a) ADF test for d, (b) ACF/PACF plots of the differenced series for p and q. Fit the model and evaluate residuals.
  2. Using the same dataset, fit a Prophet model. Add a custom annual seasonality and plot the forecast components. Compare RMSE on a 12-month holdout between ARIMA and Prophet. Which wins on this data?
  3. Implement 5-fold walk-forward validation manually using TimeSeriesSplit. For each fold, fit ARIMA(1,1,1) and record the MAE. Plot the MAE per fold to see if model performance is consistent over time or degrading (which would indicate concept drift).
  4. Generate a time series with a deliberate structural break (e.g., double the trend slope after month 36). Train ARIMA on the full series without accounting for the break, then train with a dummy variable for the post-break period. Compare RMSE. Does acknowledging the break improve the forecast?

📚 Primary Source for This Lesson

Forecasting: Principles and Practice — Chapter 9 (ARIMA Models)
Hyndman & Athanasopoulos's free online textbook is the authoritative accessible reference for ARIMA modeling. Chapter 9 covers the Box-Jenkins methodology step by step. Also highly recommended: Prophet's official quickstart documentation — the examples are clear and the documentation for holidays, seasonality, and regressors is well-written. For walk-forward validation, see scikit-learn's TimeSeriesSplit documentation.

💬 Getting non-converging ARIMA models, poor Prophet forecasts, or MAPE errors on zero values? Share your time series plot and the metrics you're seeing — your tutor will help you diagnose whether the issue is with the data, the orders, or the evaluation setup.