3  Lab 3: GAMs, Prophet, and Several Seasons at Once

Run this one live: the executable lab.

3.1 What we are doing

Lecture 3 made a claim: a GAM earns its keep when a series has several cycles running at once, because each one gets its own smooth term and you can still read the model.

Retail sales cannot show you that — they have essentially one season. So today we leave M5 and use electricity demand, which has three large cycles simultaneously, and an annual shape that will break the first model we try.

The business question: you are forecasting daily demand for a grid operator. Being wrong is expensive in both directions — under-forecast and you buy power at spot prices, over-forecast and you have paid to keep generation spinning for nothing.

Suggested time: ~40 minutes · You will need: pandas, numpy, matplotlib, prophet, pygam, and data/processed/electricity_daily.csv + electricity_hourly.csv (run python scripts/prep_electricity.py once).

---

3.2 Setup

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import warnings, logging
warnings.filterwarnings('ignore')
logging.getLogger('cmdstanpy').setLevel(logging.ERROR)   # Prophet is chatty; silence it

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_electricity.py` from the repository root first."
    )

daily = pd.read_csv(data_path("electricity_daily.csv"), parse_dates=["date"])
hourly = pd.read_csv(data_path("electricity_hourly.csv"), parse_dates=["datetime"])

# Five and a half years is plenty of yearly cycles and keeps every fit under a second.
daily = daily[daily.date >= "2013-01-01"].reset_index(drop=True)

SPLIT = "2017-08-01"
train, test = daily[daily.date < SPLIT], daily[daily.date >= SPLIT]

rmse = lambda a, b: float(np.sqrt(np.mean((np.asarray(a) - np.asarray(b)) ** 2)))
mae  = lambda a, b: float(np.mean(np.abs(np.asarray(a) - np.asarray(b))))
results = {}

print(f"{len(daily)} days — train {len(train)}, test {len(test)}")
print(f"test window: {test.date.min():%Y-%m-%d} to {test.date.max():%Y-%m-%d}")

---

Walkthrough

2041 days — train 1673, test 368
test window: 2017-08-01 to 2018-08-03

Two notes. We trim to 2013 onward — five and a half years, which is plenty of annual cycles and keeps every fit under a second, because you’re going to fit several. And we hold out just over a year, August 2017 through August 2018, so the test window contains a full summer and a full winter. That’s deliberate. Holding out only six months of an annual cycle tells you almost nothing.

The logging line silences Prophet, which is otherwise extremely chatty on every fit.

3.3 Step 1 — Three cycles, one series

Before any model, look at the three periodicities separately.

fig, axes = plt.subplots(1, 3, figsize=(12, 3.2))

hourly.groupby("hour").mw.mean().plot(ax=axes[0], marker="o", ms=3)
axes[0].set_title("hour of day"); axes[0].set_xlabel("hour")

daily.groupby("dow").mw_mean.mean().plot(ax=axes[1], marker="o", ms=4)
axes[1].set_title("day of week (0=Mon)"); axes[1].set_xlabel("")

daily.groupby("doy").mw_mean.mean().rolling(7, center=True).mean().plot(ax=axes[2])
axes[2].set_title("day of year (7d smoothed)"); axes[2].set_xlabel("")

for a in axes: a.set_ylabel("MW")
plt.tight_layout()

for label, g in [("hour of day", hourly.groupby("hour").mw.mean()),
                 ("day of week", daily.groupby("dow").mw_mean.mean()),
                 ("day of year", daily.groupby("doy").mw_mean.mean())]:
    swing = 100 * (g.max() - g.min()) / g.mean()
    print(f"{label:<14} swing {swing:5.1f}%   min {g.idxmin():>3}   max {g.idxmax():>3}")

Notice: all three cycles are large — none is a rounding error you could ignore. This is the situation Lecture 3 said GAMs are for.

---

Walkthrough

Before any model, look at the three periodicities separately.

Hour of day: demand bottoms out around four in the morning and peaks in the early evening, around hour nineteen. People wake up, go to work, come home, turn things on. Day of week: Monday highest, Sunday lowest. Commercial and industrial load.

Day of year: and this one is a strange shape, which is the whole point of Step 2.

hour of day    swing  33.3%   min   4   max  19
day of week    swing  11.6%   min   6   max   1
day of year    swing  44.7%   min 116   max 200

All three are large. The annual swing is forty-five percent of mean demand, the daily swing is thirty-three percent, and even the weakest — day of week — is nearly twelve percent. None of these is a rounding error you could ignore and still have a usable forecast.

That is the situation Lecture 3 said GAMs are for. Not one dominant season with noise around it. Three real cycles, all worth modeling, all on different time scales.

3.4 Step 2 — The annual cycle has two peaks

Look again at the third panel. It does not look like a wave.

doy = daily.groupby("doy").mw_mean.mean()
print(f"January  (days 1-31)    {doy.loc[1:31].mean():,.0f} MW   <- heating")
print(f"April    (days 95-115)  {doy.loc[95:115].mean():,.0f} MW   <- mild, nothing running")
print(f"July     (days 180-200) {doy.loc[180:200].mean():,.0f} MW   <- air conditioning")

Two peaks and two troughs per year, because heating and cooling are both electricity demand and they happen at opposite ends of the calendar.

Now the consequence. Fit the annual shape with a single sine wave, then with two.

t = np.arange(len(doy))
one = np.c_[np.ones_like(t), np.sin(2*np.pi*t/365), np.cos(2*np.pi*t/365)]
two = np.c_[one,             np.sin(4*np.pi*t/365), np.cos(4*np.pi*t/365)]

for label, X in [("one harmonic", one), ("two harmonics", two)]:
    beta, *_ = np.linalg.lstsq(X, doy.values, rcond=None)
    print(f"{label:<15} in-sample RMSE {rmse(doy.values, X @ beta):,.0f} MW")

This is the whole argument for a flexible term. One harmonic can put its peak in summer or in winter — not both. Adding the second harmonic (a cycle that repeats twice a year) cuts the error by about a factor of two and a half.

Prophet’s yearly_seasonality is exactly this idea with more terms: a Fourier series of order 10 by default, so 10 sine/cosine pairs. A pyGAM spline gets there differently — by bending freely — but solves the same problem.

---

Walkthrough

Look at that annual panel again. It does not look like a wave. It looks like a letter M.

January  (days 1-31)     34,403 MW   <- heating
April    (days 95-115)   26,945 MW   <- mild, nothing running
July     (days 180-200)  37,043 MW   <- air conditioning

There it is. High in January, low in April, high again in July. Two peaks and two troughs per year — because heating and cooling are both electricity demand, and they happen at opposite ends of the calendar. April is the month when nobody needs either.

Now the consequence, and this is the argument for the whole lecture. Fit that annual shape with a single sine wave. Then with two.

one harmonic     in-sample RMSE 3,428 MW
two harmonics    in-sample RMSE 1,348 MW

A single sine wave has one peak per cycle. It can put that peak in summer or in winter — it cannot do both. Whichever it picks, it is badly wrong for half the year. Three thousand four hundred megawatts of error on a series averaging around thirty thousand.

Add a second harmonic — a cycle that repeats twice a year — and the error drops to thirteen forty-eight. Two and a half times better, from one extra pair of terms, because that second harmonic is exactly the shape the physics needs.

That’s the mechanism. Prophet’s yearly_seasonality is this idea with more terms — a Fourier series of order ten by default, so ten sine-cosine pairs. A pyGAM spline gets there differently, by bending freely rather than by summing waves, but it’s solving the same problem: the annual shape is not a wave, so don’t fit it with one.

3.5 Step 3 — Prophet, and an honest benchmark

Prophet wants two columns named ds and y. That is most of the API.

from prophet import Prophet

# Benchmark first, always: same day last year (364 days keeps the weekday aligned).
lookup = daily.set_index("date").mw_mean
sn = [lookup.get(t - pd.Timedelta(days=364), np.nan) for t in test.date]
results["Seasonal naive (364d)"] = (rmse(test.mw_mean, sn), mae(test.mw_mean, sn))

m = Prophet(yearly_seasonality=True, weekly_seasonality=True, daily_seasonality=False)
m.fit(train.rename(columns={"date": "ds", "mw_mean": "y"})[["ds", "y"]])

fc = m.predict(test.rename(columns={"date": "ds"})[["ds"]])
results["Prophet"] = (rmse(test.mw_mean, fc.yhat), mae(test.mw_mean, fc.yhat))

for k, (r, a) in results.items():
    print(f"{k:<24} RMSE {r:7,.0f}   MAE {a:7,.0f}")

We turned daily_seasonality off deliberately — the data is one row per day, so there is no within-day variation left to model. Leave it on and Prophet will happily fit noise.

---

Walkthrough

Prophet wants two columns named ds and y. That is genuinely most of the API.

But benchmark first — always. Same day last year, and note it’s 364 days rather than 365, which keeps the weekday aligned. Comparing a Tuesday to a Wednesday would hand the benchmark an error that isn’t really about seasonality.

Seasonal naive (364d)    RMSE   4,608   MAE   3,333
Prophet                  RMSE   3,935   MAE   3,001

Prophet beats the benchmark — about fifteen percent better on RMSE.

One deliberate choice: daily_seasonality=False. The data is one row per day, so there is no within-day variation left to model. Leave it on and Prophet will happily fit a daily cycle to noise, and report it to you with a straight face.

3.6 Step 4 — The components plot

This is the payoff, and the reason Prophet gets used by people who will be questioned about their forecast.

m.plot_components(fc);

Three panels, one per additive term, each on the scale of the original series. The trend, the weekly cycle, the yearly cycle — separately readable.

Notice the yearly panel reproduces the double peak from Step 2. Nobody told Prophet about heating and cooling; the Fourier series had enough terms to find the shape.

---

Walkthrough

(Screen: run m.plot_components(fc).)

This is the payoff, and it’s why Prophet gets used by people who have to defend a forecast in a meeting.

Three panels, one per additive term, each on the scale of the original series. The trend. The weekly cycle. The yearly cycle. Separately readable, and each one is a sentence you can say to a non-technical person: “demand is drifting down slightly, Mondays run about a thousand megawatts above Sundays, and there are two annual peaks.”

And look at the yearly panel — it reproduces the double peak from Step 2. Nobody told Prophet about heating and cooling. Nobody coded a January term and a July term. The Fourier series simply had enough flexibility to find the shape that was in the data.

That combination — flexible enough to find it, structured enough to show you what it found — is the entire case for additive models.

3.7 Step 5 — How much did the extra harmonics buy?

Step 2 argued this in-sample. Test it properly, on the holdout.

yearly_seasonality accepts an integer — the Fourier order — instead of True.

Your turn: fit a Prophet with the yearly order forced to 1 (a single harmonic) and add it to the results table.

# ---- your turn: uncomment and fill in the order ----
# m1 = Prophet(yearly_seasonality=___, weekly_seasonality=True, daily_seasonality=False)
# m1.fit(train.rename(columns={"date": "ds", "mw_mean": "y"})[["ds", "y"]])
#
# fc1 = m1.predict(test.rename(columns={"date": "ds"})[["ds"]])
# results["Prophet (yearly order=1)"] = (rmse(test.mw_mean, fc1.yhat),
#                                        mae(test.mw_mean, fc1.yhat))
#
# print(f"single harmonic RMSE {rmse(test.mw_mean, fc1.yhat):,.0f}")

Report both numbers. Does the holdout agree with what Step 2 predicted in-sample?

---

Walkthrough

Step 2 made this argument in-sample, on the averaged annual profile. That’s suggestive, not evidence. Test it properly, on the holdout. yearly_seasonality takes an integer — the Fourier order — instead of True.

Try it

Pause here. Fit a Prophet with the yearly order forced to 1 — a single harmonic — and add it to the results table. Everything else stays the same.

single harmonic RMSE 4,269

Four thousand two hundred and sixty-nine, against three thousand nine hundred and thirty-five for the full-order model. The holdout agrees with Step 2 — crippling the annual term costs you real accuracy out of sample, not just in-sample fit.

But look where it lands relative to the benchmark. Seasonal naive was 4,608. The single-harmonic Prophet is 4,269 — it beats the benchmark by about seven percent, and that is all it beats it by.

And here’s the part worth stopping on. On MAE, it loses to the benchmark outright — 3,358 against 3,333. Same two models, same holdout, opposite verdict depending on which error measure you quote.

That is not a trick and it is not a bug. RMSE squares errors, so it punishes the occasional large miss; MAE treats all errors proportionally. A model can be better on one and worse on the other, and when a ranking flips between them the honest move is to say so rather than to quote the metric that flatters you.

3.8 Step 6 — pyGAM, and reading the smooth directly

Prophet is a GAM with the time-series decisions made for you. pyGAM makes you state the terms, which is more work and more control.

from pygam import LinearGAM, s, f

# s(0) = cyclic spline on day-of-year (basis="cp" makes Dec 31 join Jan 1)
# f(1) = day-of-week as a categorical factor
X_tr, X_te = train[["doy", "dow"]].values, test[["doy", "dow"]].values

gam = LinearGAM(s(0, n_splines=20, basis="cp") + f(1)).fit(X_tr, train.mw_mean.values)
results["pyGAM (doy + dow)"] = (rmse(test.mw_mean, gam.predict(X_te)),
                                mae(test.mw_mean, gam.predict(X_te)))

print(f"effective degrees of freedom: {gam.statistics_['edof']:.1f}")
print(pd.DataFrame(results, index=["RMSE", "MAE"]).T
        .sort_values("RMSE").round(0).astype(int).to_string())

Now the partial dependence — the GAM equivalent of the components plot.

fig, axes = plt.subplots(1, 2, figsize=(11, 3.2))
for i, (ax, name) in enumerate(zip(axes, ["day of year", "day of week"])):
    XX = gam.generate_X_grid(term=i)
    pdep, ci = gam.partial_dependence(term=i, X=XX, width=0.95)
    ax.plot(XX[:, i], pdep)
    ax.plot(XX[:, i], ci, c="grey", ls="--", lw=0.8)
    ax.set_title(name); ax.set_ylabel("effect on MW")
plt.tight_layout()

Read the left panel and say what it means in words, the way you would to a grid operator. The confidence band widens where the data is thin — say that too.

---

Walkthrough

Prophet is a GAM with the time-series decisions already made for you. pyGAM makes you state the terms yourself — more work, more control.

Two terms only. A cyclic spline on day-of-year — basis="cp" is what makes December 31st join back up with January 1st, so the smooth doesn’t have a discontinuity at the year boundary — and day-of-week as a categorical factor. That’s it. No trend term. No hourly anything.

effective degrees of freedom: 22.5

                          RMSE    MAE
pyGAM (doy + dow)         3444   2550
Prophet                   3935   3001
Prophet (yearly order=1)  4269   3358
Seasonal naive (364d)     4608   3333

Every model beat the benchmark. The single-harmonic Prophet is the worst of the three, as Step 2 predicted.

And pyGAM wins. Thirty-four forty-four against Prophet’s thirty-nine thirty-five — about twelve percent better, on both RMSE and MAE, using two terms against Prophet’s three, with no trend component at all. That should bother you slightly, and we’ll come back to it in the closing question.

Now the partial dependence — the GAM equivalent of the components plot.

Say what the left panel means in words, the way you’d say it to a grid operator: demand is high in January, falls through spring to a minimum around mid-April, climbs to an annual peak in mid-July, and falls again into autumn. The model learned the M shape directly from the data by bending, without being given a single Fourier term.

And say the other thing too: the confidence band widens where the data is thin. Look at the edges. That widening is the model telling you where it is guessing, and a forecast that comes with an honest statement of its own uncertainty is worth more than a point estimate that doesn’t.

3.9 Before you leave

You should have found that every method beat the seasonal-naive benchmark on RMSE, and that the single harmonic was the worst of the three models.

Look at the MAE column before you move on. The single-harmonic Prophet has a worse MAE than the benchmark it beat on RMSE. Same two models, same holdout, opposite verdict depending on which error measure you quote — because RMSE squares errors and so punishes the occasional large miss, while MAE treats every error proportionally.

Discuss: pyGAM was given only day-of-year and day-of-week — no trend term at all. Prophet was given a trend, a weekly cycle, and a yearly cycle with ten harmonics. Yet on this holdout the simpler model won, by roughly 12% on both RMSE and MAE.

Give one reason a more flexible model can lose on a held-out year, and name the thing you would check before concluding that pyGAM is genuinely the better choice here.

And the operational question: your forecast is for a grid operator who must buy generation in advance. Under-forecasting and over-forecasting cost different amounts. Which of these models would you actually ship, and what would you need to know first?

---

Solutions posted to Canvas after class.

Walkthrough

Try it

Post these in the Week 3 thread before Lecture 4.

First, the awkward result. pyGAM was given two things: day of year and day of week. Prophet was given a trend, a weekly cycle, and a yearly cycle with ten harmonics. Prophet has strictly more machinery — and it lost by twelve percent.

Give me one reason a more flexible model can lose on a held-out year. And name the thing you would check before concluding that pyGAM is genuinely the better choice here, rather than luckier on this particular holdout.

Second, the operational question. Your forecast is for a grid operator who must buy generation in advance, and the two directions of error cost different amounts. Under-forecast and you buy at spot during a peak. Over-forecast and you’ve paid for spinning reserve you didn’t need.

Which of these models would you actually ship — and what would you need to know before you could answer that at all?

Notice that nothing in the table we just built answers the second question. RMSE and MAE are both symmetric — they charge you the same for being a thousand megawatts high as a thousand low. If the costs aren’t symmetric, neither metric is measuring the thing you care about.

That’s not a flaw in this lab. It’s the gap between forecast accuracy and forecast value, and it’s worth carrying into the rest of the course.

Solutions for the coded parts go up on Canvas after the deadline. See you in Lecture 4, where we start on trees.

3.10 What your run should produce

Your numbers will differ where a model is stochastic; these are the values this lab produces on the pinned environment.

Reproducible; verified bit-identical across repeated runs. Prophet uses MAP estimation, not MCMC.

Quantity	Value
Data	2,041 days (train 1,673, test 368), 2017-08-01 → 2018-08-03
Swing — hour of day	33.3% (min hour 4, max hour 19)
Swing — day of week	11.6% (min Sunday, max Monday)
Swing — day of year	44.7% (min day 116, max day 200)
January / April / July	34,403 / 26,945 / 37,043 MW
One harmonic, in-sample	3,428 MW
Two harmonics, in-sample	1,348 MW (factor of 2.5, not 3)
pyGAM (doy + dow)	RMSE 3,444 · MAE 2,550
Prophet	RMSE 3,935 · MAE 3,001
Prophet (yearly order = 1)	RMSE 4,269 · MAE 3,358
Seasonal naive (364d)	RMSE 4,608 · MAE 3,333
pyGAM effective d.o.f.	22.5