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)))
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.
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
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('whole series: weekday avg %d MW vs weekend avg %d MW' % (s[s.index.dayofweek<5].mean(), s[s.index.dayofweek>=5].mean()))
whole series: 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, with weekends shaded. Across all four years weekdays average about 1,113 MW and weekends about 911 MW. A good model must reproduce both.
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.
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.
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.
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.
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%.
nmae = mae(test, snaive) # rel_MAE = MAE / the baseline's MAE on this holdout (MASE would scale by the in-sample naive error)
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),'rel_MAE':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_% rel_MAE 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, with an MAE a fifth of the baseline's), 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.
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 over 14 lags flags a little leftover autocorrelation (p = 0.03). It is spread thinly across the lags, with no single one above about 0.05 and almost none at lag 1, and an AR(3) on the errors in Take It Further leaves none that is significant. 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.
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. Take It Further backtests the band in four separate 30-day windows, and 117 of 120 actual days land inside it, a little more than the 95% promised, so for planning it errs slightly wide.
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¶
- Define the decision (daily demand for grid scheduling) and the baseline.
- Collect the daily series.
- Inspect: spot both the weekly and yearly rhythms.
- Decompose with MSTL to separate the two seasonal cycles.
- Stationarity check: trend-stationary, so model directly; the real challenge is multiple seasonality.
- Split in time and set a weekly seasonal-naive baseline.
- One season: Holt-Winters captures the week but misses the year.
- Both seasons: harmonic regression with day-of-week dummies + yearly Fourier terms.
- Validate: the two-season model wins clearly on the holdout.
- Interpret: clean residuals and readable weekday/weekend effects.
- Deploy: refit on all data, forecast ahead with an interval.
- 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.