Lab 1: Benchmarks and Exponential Smoothing

ECON 8310: Business Forecasting — University of Nebraska at Omaha

What we are doing

Lecture 1 made a claim that sounded like housekeeping: always beat a benchmark first. In the next 35 minutes you are going to test that claim on a real Walmart series, and then watch it fail on a different one.

The business question: you are forecasting weekly HOUSEHOLD-goods demand at store TX_1 to set inventory a year out. Which method should you actually ship?

Suggested time: ~35 minutes · You will need: pandas, numpy, matplotlib, statsmodels, and data/processed/m5_weekly.csv (run python scripts/prep_m5.py if it is missing).

Setup

Nothing to fill in here — run it and move on.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.holtwinters import ExponentialSmoothing, SimpleExpSmoothing, Holt
import warnings
warnings.filterwarnings('ignore')

from pathlib import Path

def data_path(name):
    """Find the processed data file whether you render from the repo root or from
    this file's own folder."""
    for base in (Path.cwd(), *Path.cwd().parents):
        p = base / "data" / "processed" / name
        if p.exists():
            return p
    raise FileNotFoundError(
        f"Could not find data/processed/{name}. "
        "Run `python scripts/prep_m5.py` from the repository root first."
    )

weekly = pd.read_csv(data_path("m5_weekly.csv"), parse_dates=["week_start"])

def get_series(store, cat):
    return (weekly[(weekly.store_id == store) & (weekly.cat_id == cat)]
            .set_index("week_start")["units"].sort_index().astype(float))

y = get_series("TX_1", "HOUSEHOLD")

N_TEST = 52
train, test = y.iloc[:-N_TEST], y.iloc[-N_TEST:]
print(f"{len(y)} weeks total — train {len(train)}, test {len(test)}")

277 weeks total — train 225, test 52

Step 1 — Look at it first

Never fit anything before you have seen the series.

fig, ax = plt.subplots(figsize=(11, 3.5))
ax.plot(train.index, train.values, label="train", lw=1.4)
ax.plot(test.index, test.values, label="test (held out)", lw=1.4)
ax.set_ylabel("units / week"); ax.set_title("TX_1 HOUSEHOLD")
ax.legend(); plt.tight_layout()

Notice: the series roughly doubles over five years, and the level is still moving at the end. Hold that thought — it is the whole lab.

Step 2 — Establish the benchmark

Before any model, the thing to beat. Seasonal naive says: this week next year looks like this week last year.

def rmse(actual, pred):
    return float(np.sqrt(np.mean((np.asarray(actual) - np.asarray(pred)) ** 2)))

def mae(actual, pred):
    return float(np.mean(np.abs(np.asarray(actual) - np.asarray(pred))))

seasonal_naive = y.shift(52).iloc[-N_TEST:].values

results = {"Seasonal naive": (rmse(test, seasonal_naive), mae(test, seasonal_naive))}
print(f"Seasonal naive  RMSE {results['Seasonal naive'][0]:,.0f}")

Seasonal naive  RMSE 973

Your turn: before running anything else, write down a number — what RMSE do you think a proper exponential smoothing model will get? You will check it in Step 4.

Step 3 — Fit the three models

All three calls come straight from Lecture 1. Run this, then predict the ordering before you look at Step 4 — write your guess down.

# SES — level only.
ses_fit = SimpleExpSmoothing(train).fit()
results["SES"] = (rmse(test, ses_fit.forecast(N_TEST)),
                  mae(test, ses_fit.forecast(N_TEST)))

# Holt — level + trend. Undamped for now; we come back to damping in Step 5.
holt_fit = Holt(train).fit()
results["Holt"] = (rmse(test, holt_fit.forecast(N_TEST)),
                   mae(test, holt_fit.forecast(N_TEST)))

# Holt-Winters — level + trend + seasonality, additive, period 52.
hw_fit = ExponentialSmoothing(train, trend="add",
                              seasonal="add", seasonal_periods=52).fit()
results["Holt-Winters"] = (rmse(test, hw_fit.forecast(N_TEST)),
                           mae(test, hw_fit.forecast(N_TEST)))

One constraint to respect

Holt-Winters must estimate a seasonal index for each of the 52 weeks, and it needs two full cycles — 104 observations — to initialize. You have 225 training weeks, so you are fine here. You will not always be.

Step 4 — Compare honestly

tbl = (pd.DataFrame(results, index=["RMSE", "MAE"]).T
         .sort_values("RMSE").round(0).astype(int))
print(tbl.to_string())

                RMSE   MAE
SES              427   352
Holt-Winters     885   814
Seasonal naive   973   848
Holt            1166  1029

Stop before reading on: did any model beat the benchmark? By how much? And was it the model you expected in Step 3?

In the room: say your answer to the person next to you, and make them say theirs.
On your own: score yourself against what you already wrote down. Scroll back to the RMSE number from Step 2 and the ordering from Step 3, and write two lines in a comment:

# guess ____, actual ____, off by ____
# ordering: right / wrong — first model I misplaced was ____

Grading your own prediction in writing is the part that matters. It is what makes the next result inform you rather than just agree with you.

Step 5 — Break it, then fix it

Look at where Holt’s forecast goes.

fig, ax = plt.subplots(figsize=(11, 3.5))
ax.plot(test.index, test.values, label="actual", lw=1.6)
ax.plot(test.index, ses_fit.forecast(N_TEST), label="SES", ls="--")
ax.plot(test.index, holt_fit.forecast(N_TEST), label="Holt (undamped)", ls="--")
ax.set_ylabel("units / week"); ax.legend(); plt.tight_layout()

Lecture 1 warned that undamped linear extrapolation gets unrealistic past about six periods. You are forecasting 52.

Your turn: refit Holt with a damped trend and add it to the table. One keyword argument.

# ---- your turn: uncomment and fill in the one keyword argument ----
# holt_damped = Holt(train, ___).fit()             # <- damped_trend=True
#
# results["Holt (damped)"] = (rmse(test, holt_damped.forecast(N_TEST)),
#                             mae(test, holt_damped.forecast(N_TEST)))
#
# print(pd.DataFrame(results, index=["RMSE", "MAE"]).T
#         .sort_values("RMSE").round(0).astype(int).to_string())
# print(f"estimated phi: {holt_damped.params['damping_trend']:.3f}")

Report the value of \(\phi\) — then check it against the range statsmodels searches, which is \(0.8 \le \phi \le 0.995\). If yours reads exactly 0.800, the optimizer stopped at the edge of that range: it wanted more damping than it was allowed. That is the bound talking, not an estimate. Now compare the damped model’s RMSE with SES’s. What does that suggest the damped trend is actually doing?

Step 6 — Same store. Different aisle. Opposite answer.

Change one word — HOUSEHOLD becomes FOODS. Same store, same code, same horizon.

y2 = get_series("TX_1", "FOODS")
train2, test2 = y2.iloc[:-N_TEST], y2.iloc[-N_TEST:]

sn2 = y2.shift(52).iloc[-N_TEST:].values
ses2 = SimpleExpSmoothing(train2).fit()

print(f"TX_1 FOODS   seasonal naive RMSE {rmse(test2, sn2):,.0f}")
print(f"TX_1 FOODS   SES            RMSE {rmse(test2, ses2.forecast(N_TEST)):,.0f}")

TX_1 FOODS   seasonal naive RMSE 1,308
TX_1 FOODS   SES            RMSE 1,474

Notice: the ranking flips. On the household aisle an exponential smoother won comfortably. On the food aisle of the same store, the benchmark wins and SES is about 13% worse.

The reason is visible in one number. For TX_1 HOUSEHOLD the test year sits about 15% above the year before it, so seasonal naive carries a stale level forward and loses badly. For TX_1 FOODS the level barely moves — about 1.5% down — and last year is an excellent guess.

So it is not the store, and it is not the method. It is whether the level moved.

Before you leave

Discuss: You now have two series where the best method is different. You are responsible for all 30 series in this dataset and you have to pick one approach to run in production every month.

What do you actually do — and what would you need to measure to decide? “Fit everything and pick the winner per series” is an answer, but say what it costs you.

Solutions posted to Canvas after class.