←Back to the bookChapter 133 · Energy Demand Forecasting

Chapter 133 · Forecasting Case Study

Case Study: Energy Demand Forecasting

The multiple-seasonality forecasting problem on four years of daily electricity demand: separate the weekly and yearly cycles with MSTL, then beat a one-season Holt-Winters with a harmonic regression that captures both seasons using day-of-week dummies and yearly Fourier terms, backtested and graded with the full accuracy panel.

Statistics, Data Science and AI: A Visual Handbook · John Fisher · 2026

Case Study: Energy Demand Forecasting

Electricity demand is the classic multiple-seasonality forecasting problem: it rises and falls with the day of the week (weekday vs weekend) AND with the time of year (winter heating, summer cooling), both at once. That breaks the one-season tools like Holt-Winters and SARIMA. This case study takes four years of daily demand, decomposes the two cycles with MSTL, and forecasts with a harmonic regression that captures both, then grades it against a baseline and a one-season model.

import numpy as np, pandas as pd, matplotlib.pyplot as plt
import seaborn as sns   # seaborn = high-level statistical plots (heatmaps, pairplots, count/bar plots)
from matplotlib.colors import ListedColormap
EM="#0284c7"; DEEP="#075985"; LIGHT="#bae6fd"; INK="#1a2138"; GRID="#e6e9f2"; RED="#ef4444"; AMBER="#d97706"; GREEN="#059669"; BLUE="#2563eb"; PUR="#9333ea"; GREY="#94a3b8"; SLATE="#475569"; ORG="#0284c7"; CYAN="#0891b2"
plt.rcParams.update({"figure.facecolor":"white","axes.facecolor":"white","figure.dpi":110,"font.size":11,
   "axes.edgecolor":GRID,"axes.grid":True,"grid.color":GRID,"axes.axisbelow":True,"axes.spines.top":False,
   "axes.spines.right":False,"axes.titlesize":12,"axes.titleweight":"bold","legend.frameon":False})
sns.set_style("whitegrid")
BASE="https://raw.githubusercontent.com/johnfisher-ai/Statistics-Data-Science-AI-Visual-Book/main/data/"

import statsmodels.api as sm
from statsmodels.tsa.seasonal import MSTL
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.holtwinters import ExponentialSmoothing
from statsmodels.stats.diagnostic import acorr_ljungbox
import warnings; warnings.filterwarnings('ignore')
def mape(a,f): return float(np.mean(np.abs((a.values-np.asarray(f))/a.values))*100)
def mae(a,f): return float(np.mean(np.abs(a.values-np.asarray(f))))
def rmse(a,f): return float(np.sqrt(np.mean((a.values-np.asarray(f))**2)))

STEP 1 · DEFINE

What are we forecasting, and why?

Daily demand, to keep the grid balanced.

A grid operator must forecast daily electricity demand to schedule generation, buy power, and avoid both shortfalls and waste. Success: beat a seasonal-naive baseline, keep the error (MAPE) low, and, critically, capture both seasonal rhythms, missing either the weekend dip or the summer peak means mis-scheduling capacity.

STEP 2 · COLLECT

Four years of daily demand

One value per day.

try: raw = pd.read_excel('../../data/case-study-energy-demand-forecasting--energy_demand.xlsx', sheet_name='Data')
except FileNotFoundError: raw = pd.read_excel(BASE + 'case-study-energy-demand-forecasting--energy_demand.xlsx', sheet_name='Data')
raw['date']=pd.to_datetime(raw['date']); s = raw.set_index('date')['demand_mw'].asfreq('D')
print('days:', len(s), '| from', s.index.min().date(), 'to', s.index.max().date())

days: 1461 | from 2020-01-01 to 2023-12-31

STEP 3 · INSPECT

Two rhythms at once

Zoom out for the year, zoom in for the week.

fig,ax=plt.subplots(2,1,figsize=(11,6))
ax[0].plot(s.index, s.values, color=EM, lw=0.7); ax[0].set(title='Full series: the YEARLY cycle (winter + summer peaks)', ylabel='MW')
wk=s['2022-03-01':'2022-03-28']; ax[1].plot(wk.index, wk.values, 'o-', color=EM, ms=4); ax[1].set(title='Four-week zoom: the WEEKLY cycle (weekday high, weekend low)', ylabel='MW')
for d in wk.index[wk.index.dayofweek>=5]: ax[1].axvline(d, color=AMBER, alpha=0.15)
plt.tight_layout(); plt.show()
print('weekday avg %d MW vs weekend avg %d MW' % (s[s.index.dayofweek<5].mean(), s[s.index.dayofweek>=5].mean()))

weekday avg 1113 MW vs weekend avg 911 MW

Two patterns are layered on top of each other. Zoomed out, a slow yearly cycle with winter and summer peaks; zoomed into a month, a sharp weekly cycle, weekdays run about 1,113 MW and weekends drop to about 911 MW (shaded). A good model must reproduce both.

STEP 4 · DECOMPOSE

Separate the two seasons with MSTL

One tool, multiple seasonal periods.

mstl = MSTL(s, periods=(7, 365)).fit()
fig = mstl.plot(); fig.set_size_inches(10,8); plt.tight_layout(); plt.show()
print('weekly component amplitude ~ %.0f MW | yearly component amplitude ~ %.0f MW' % (mstl.seasonal['seasonal_7'].std(), mstl.seasonal['seasonal_365'].std()))

weekly component amplitude ~ 94 MW | yearly component amplitude ~ 98 MW

Ordinary decomposition handles one season; MSTL (Multiple STL) peels off several. It splits the series into a trend, a weekly component, a yearly component, and a small residual. Both seasonal swings are sizable (about 94 and 98 MW), so a model that captures only one will leave real, systematic error on the table.

STEP 5 · STATIONARITY & THE CHALLENGE

No differencing needed, but two seasons

The real difficulty is the multiple seasonality.

print('ADF p = %.3f -> %s' % (adfuller(s)[1], 'trend-stationary (no unit root)' if adfuller(s)[1]<0.05 else 'non-stationary'))
print('the trend is mild, so we model trend + seasonality DIRECTLY rather than differencing.')
print('the genuine challenge: TWO seasonal periods (7 and 365) that one-season models cannot both capture.')

ADF p = 0.008 -> trend-stationary (no unit root)
the trend is mild, so we model trend + seasonality DIRECTLY rather than differencing.
the genuine challenge: TWO seasonal periods (7 and 365) that one-season models cannot both capture.

Unlike the retail series, this one is trend-stationary (ADF p < 0.05): the trend is gentle, so we model it directly instead of differencing. The real obstacle is the double seasonality. Holt-Winters and SARIMA each support only one seasonal period, so they can fit the weekly cycle or the yearly one, but not both. That is the modeling decision this chapter turns on.

STEP 6 · SPLIT & BASELINE

Hold out 60 days, set the bar

Chronological, and the number to beat.

H = 60; train, test = s[:-H], s[-H:]
snaive = np.array([train.iloc[-7 + (i % 7)] for i in range(H)])   # repeat the last week
print('train %d days | test %d days' % (len(train), len(test)))
print('weekly seasonal-naive baseline MAPE: %.2f%%  (captures the week, misses the year and trend)' % mape(test, snaive))

train 1401 days | test 60 days
weekly seasonal-naive baseline MAPE: 10.42%  (captures the week, misses the year and trend)

We hold out the last 60 days and split in time. The weekly seasonal-naive forecast (repeat the last week) captures the day-of-week pattern but nothing else, and lands at about 10.4% MAPE. That is the bar.

STEP 7 · ONE SEASONALITY

Holt-Winters (weekly)

Captures the week and the trend, blind to the year.

hw = ExponentialSmoothing(train, trend='add', seasonal='add', seasonal_periods=7).fit()
hw_fc = hw.forecast(H)
print('Holt-Winters (weekly) MAPE: %.2f%%  -> big improvement, but it has no yearly term' % mape(test, hw_fc))

Holt-Winters (weekly) MAPE: 3.51%  -> big improvement, but it has no yearly term

Holt-Winters with a weekly period captures the weekday/weekend swing and the trend, cutting the error to about 3.5%. But it has no yearly component, so as the 60-day horizon drifts across the seasonal calendar, it slowly loses the yearly signal. To do better we need both seasons.

STEP 8 · BOTH SEASONALITIES

Harmonic regression with Fourier terms

Day-of-week dummies + yearly sine/cosine.

def design(idx, t0, K=3):
    trend = (np.arange(len(idx)) + t0).reshape(-1,1)
    dow = np.column_stack([(idx.dayofweek.values == k).astype(float) for k in range(1,7)])   # weekly dummies
    doy = idx.dayofyear.values
    fourier = np.column_stack([f(2*np.pi*k*doy/365.25) for k in range(1,K+1) for f in (np.sin, np.cos)])   # yearly
    return sm.add_constant(np.hstack([trend, dow, fourier]))
Xtr, Xte = design(train.index, 0), design(test.index, len(train))
hr = sm.OLS(train.values, Xtr).fit(); hr_fc = hr.predict(Xte)
print('harmonic regression (weekly + yearly) MAPE: %.2f%%  <- captures BOTH seasons' % mape(test, hr_fc))

harmonic regression (weekly + yearly) MAPE: 2.15%  <- captures BOTH seasons

The fix for multiple seasonality is a harmonic regression: a linear model with a trend, day-of-week dummies for the weekly cycle, and Fourier terms (pairs of sines and cosines) for the yearly cycle. Fourier terms are how you bolt a long seasonal period onto any model, a few harmonics trace the smooth annual curve. This captures both rhythms and drops the error to about 2.2%.

STEP 9 · VALIDATE

Grade all three on the holdout

One picture, one table.

nmae = mae(test, snaive)
rows=[]
for name,f in [('weekly seasonal-naive',snaive),('Holt-Winters (1 season)',hw_fc),('harmonic reg (2 seasons)',hr_fc)]:
    rows.append({'model':name,'MAE':round(mae(test,f)),'RMSE':round(rmse(test,f)),'MAPE_%':round(mape(test,f),2),'MASE':round(mae(test,f)/nmae,2)})
print(pd.DataFrame(rows).to_string(index=False))
fig,ax=plt.subplots(figsize=(11,4)); ax.plot(s.index[-90:], s.values[-90:], color='k', lw=1.6, label='actual')
ax.plot(test.index, hw_fc, '--', color=AMBER, label='Holt-Winters (1 season)'); ax.plot(test.index, hr_fc, '--', color=EM, lw=2, label='harmonic reg (2 seasons)')
ax.set(title='Backtest: capturing BOTH seasons wins', ylabel='MW'); ax.legend(); plt.tight_layout(); plt.show()

                   model  MAE  RMSE  MAPE_%  MASE
   weekly seasonal-naive  127   147   10.42  1.00
 Holt-Winters (1 season)   44    54    3.51  0.34
harmonic reg (2 seasons)   25    33    2.15  0.20

The verdict panel is unambiguous: the two-season harmonic regression wins (about 2.2% MAPE, MASE well below 1), beating the one-season Holt-Winters and cutting the baseline error by nearly 80%. The plot shows why: Holt-Winters gets the weekly zig-zag right but drifts off the yearly level, while the harmonic model tracks both.

STEP 10 · INTERPRET

Are the residuals clean?

And what does the model say?

resid = train.values - hr.predict(Xtr)
lb = acorr_ljungbox(resid, lags=[14])['lb_pvalue'].iloc[0]
print('harmonic-regression residual Ljung-Box p = %.2f -> %s' % (lb, 'little structure left' if lb>0.05 else 'some autocorrelation remains'))
dow_names=['Tue','Wed','Thu','Fri','Sat','Sun']
print('\nweekday effects vs Monday (MW):')
for name,coef in zip(dow_names, hr.params[2:8]): print('  %-4s %+6.0f' % (name, coef))

harmonic-regression residual Ljung-Box p = 0.03 -> some autocorrelation remains

weekday effects vs Monday (MW):
  Tue     +19
  Wed     +20
  Thu      +8
  Fri     -12
  Sat    -174
  Sun    -213

The residuals are small, though a Ljung-Box test flags a little leftover day-to-day autocorrelation, a short AR term (dynamic regression, in Take It Further) would mop it up. The fitted coefficients are readable: the weekend effect is a large negative shift (Saturday about 174 MW and Sunday about 213 MW below a Monday), and the Fourier terms trace the annual curve. A model you can both trust and explain is exactly what an operations team needs.

STEP 11 · DEPLOY

Refit on all data, forecast ahead

With a prediction interval.

Xall = design(s.index, 0); final = sm.OLS(s.values, Xall).fit()
fut_idx = pd.date_range(s.index[-1] + pd.Timedelta(days=1), periods=30, freq='D')
pred = final.get_prediction(design(fut_idx, len(s))).summary_frame(alpha=0.05)
print('next 30 days: average %d MW | day 1 = %d MW  (95%% interval %d to %d)' % (pred['mean'].mean(), pred['mean'].iloc[0], pred['obs_ci_lower'].iloc[0], pred['obs_ci_upper'].iloc[0]))
fig,ax=plt.subplots(figsize=(11,4)); ax.plot(s.index[-120:], s.values[-120:], color=GREY, lw=1.2, label='history')
ax.plot(fut_idx, pred['mean'], 'o-', color=EM, ms=3, label='forecast')
ax.fill_between(fut_idx, pred['obs_ci_lower'], pred['obs_ci_upper'], color=EM, alpha=0.18, label='95% interval')
ax.set(title='Deployed forecast: the next 30 days with uncertainty', ylabel='MW'); ax.legend(); plt.tight_layout(); plt.show()

next 30 days: average 1312 MW | day 1 = 1343 MW  (95% interval 1281 to 1404)

Before forecasting the real future we refit on all four years. The model projects the next month day by day, weekend dips and all, around an average near 1,310 MW, with a 95% interval on each day. In production this reruns every morning as yesterday's actuals land, so the grid team always has a fresh, calendar-aware forecast with its uncertainty. One honesty note on the band: because the residuals still carry a little day-to-day autocorrelation (the Ljung-Box flag above), each day's 95% interval is a touch narrower than the true uncertainty; the AR error model in Take It Further widens it back to honest coverage.

STEP 12 · COMMUNICATE

One paragraph for the control room

Plain English, with the two rhythms.

For the operations team: demand follows two clocks at once, a weekly one (weekdays run about 200 MW above weekends) and a yearly one (winter and summer peaks). Our forecast captures both and was accurate to about 2.2% over the last two months, far better than repeating last week. Expect the next month to average around 1,310 MW; schedule extra capacity for weekdays and the seasonal peak, and pull back on weekends. Plan against the daily range, and the forecast refreshes every morning.

The forecasting method, in one view

1. Define the decision (daily demand for grid scheduling) and the baseline.
2. Collect the daily series.
3. Inspect: spot both the weekly and yearly rhythms.
4. Decompose with MSTL to separate the two seasonal cycles.
5. Stationarity check: trend-stationary, so model directly; the real challenge is multiple seasonality.
6. Split in time and set a weekly seasonal-naive baseline.
7. One season: Holt-Winters captures the week but misses the year.
8. Both seasons: harmonic regression with day-of-week dummies + yearly Fourier terms.
9. Validate: the two-season model wins clearly on the holdout.
10. Interpret: clean residuals and readable weekday/weekend effects.
11. Deploy: refit on all data, forecast ahead with an interval.
12. Communicate both rhythms and the range in plain English.

The one idea to keep: when a series has more than one seasonal cycle, a one-season model (Holt-Winters, SARIMA) can only capture one, so reach for MSTL to see them and Fourier terms to model them all at once.

---

Statistics, Data Science and AI: A Visual Handbook · © 2026 John Fisher

←Back to the bookChapter 133 · Energy Demand Forecasting