Energy Demand Forecasting: Take It Further¶
Five extensions of the Chapter 138 model, each sharpening the multiple-seasonality forecast: choosing the number of Fourier harmonics, a rolling-origin backtest, measuring how much each seasonality contributes, adding an AR term to clean the residuals (dynamic regression), and checking the prediction-interval coverage. We start by rebuilding the series and the design matrix.
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.statespace.sarimax import SARIMAX
from statsmodels.stats.diagnostic import acorr_ljungbox
import warnings; warnings.filterwarnings('ignore')
def mape(a,f): return float(np.mean(np.abs((np.asarray(a)-np.asarray(f))/np.asarray(a)))*100)
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')
def design(idx, t0, K=3, weekly=True, yearly=True):
cols=[(np.arange(len(idx))+t0).reshape(-1,1)]
if weekly: cols.append(np.column_stack([(idx.dayofweek.values==k).astype(float) for k in range(1,7)]))
if yearly:
doy=idx.dayofyear.values; cols.append(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)]))
return sm.add_constant(np.hstack(cols))
H=60; train,test=s[:-H],s[-H:]; print('rebuilt:', len(s), 'days')
rebuilt: 1461 days
for K in [1,2,3,5,8]:
m=sm.OLS(train.values, design(train.index,0,K=K)).fit()
print('K=%d (%d yearly terms) AIC %.0f holdout MAPE %.2f%%' % (K, 2*K, m.aic, mape(test, m.predict(design(test.index,len(train),K=K)))))
print('add harmonics until AIC / holdout error stops improving -> a smooth annual curve, not a wiggly one')
K=1 (2 yearly terms) AIC 15466 holdout MAPE 3.31% K=2 (4 yearly terms) AIC 13614 holdout MAPE 2.16% K=3 (6 yearly terms) AIC 13617 holdout MAPE 2.15% K=5 (10 yearly terms) AIC 13625 holdout MAPE 2.15% K=8 (16 yearly terms) AIC 13633 holdout MAPE 2.15% add harmonics until AIC / holdout error stops improving -> a smooth annual curve, not a wiggly one
Each Fourier harmonic is a sine-cosine pair; more harmonics trace a wigglier annual curve. Too few cannot bend to the winter-and-summer shape; too many start fitting noise. Pick the order where AIC and holdout error stop improving, the Fourier equivalent of choosing a model's complexity.
scores=[]
for o in range(len(s)-180, len(s)-60, 30):
tr, te = s.iloc[:o], s.iloc[o:o+30]
m=sm.OLS(tr.values, design(tr.index,0)).fit()
scores.append(mape(te, m.predict(design(te.index,len(tr)))))
print('rolling-origin MAPE over %d windows: %s' % (len(scores), [round(x,2) for x in scores]))
print('average %.2f%% -> a more trustworthy estimate than a single holdout' % np.mean(scores))
rolling-origin MAPE over 4 windows: [2.79, 2.08, 1.99, 1.98] average 2.21% -> a more trustworthy estimate than a single holdout
A single 60-day holdout can flatter or punish the model by luck. Rolling-origin evaluation re-forecasts from several cutoffs and averages, the honest way to grade any forecaster, and it confirms the harmonic model's accuracy is stable, not a fluke of one window.
for name,wk,yr in [('trend only',False,False),('weekly only',True,False),('yearly only',False,True),('both',True,True)]:
m=sm.OLS(train.values, design(train.index,0,weekly=wk,yearly=yr)).fit()
print('%-12s MAPE %.2f%%' % (name, mape(test, m.predict(design(test.index,len(train),weekly=wk,yearly=yr)))))
print('each seasonality removes real error; you need BOTH to reach the best accuracy')
trend only MAPE 11.24% weekly only MAPE 8.21% yearly only MAPE 8.82% both MAPE 2.15% each seasonality removes real error; you need BOTH to reach the best accuracy
An ablation makes the case concretely: dropping either seasonality hurts. Weekly-only and yearly-only each leave a chunk of systematic error, and only the both model reaches the best accuracy. That is the whole argument for multiple-seasonality modeling in one table.
from statsmodels.tsa.ar_model import AutoReg
from statsmodels.tsa.stattools import acf
ols = sm.OLS(train.values, design(train.index,0)).fit()
resid = train.values - ols.predict(design(train.index,0))
lb_ols = acorr_ljungbox(resid, lags=[14])['lb_pvalue'].iloc[0]
ar = AutoReg(resid, lags=3).fit()
lb_ar = acorr_ljungbox(ar.resid, lags=[14])['lb_pvalue'].iloc[0]
print('plain harmonic-regression residual Ljung-Box p = %.3f (some autocorrelation left)' % lb_ols)
print('residual autocorrelation, lags 1-14:', np.round(acf(resid, nlags=14)[1:], 3))
print('after modeling the errors with an AR(3): p = %.3f (no significant autocorrelation left)' % lb_ar)
plain harmonic-regression residual Ljung-Box p = 0.031 (some autocorrelation left) residual autocorrelation, lags 1-14: [ 0.015 -0.053 -0.039 -0.019 -0.053 0.036 0.036 0.041 -0.01 -0.033 -0.005 -0.037 0.05 0.031] after modeling the errors with an AR(3): p = 0.179 (no significant autocorrelation left)
The harmonic regression captures the seasonality but leaves mild autocorrelation spread thinly across the first two weeks of lags (lag 1 is only 0.015, and none of the first fourteen exceeds about 0.05). Dynamic regression, the deterministic features plus an AR model on the residual errors, mops that up: the Ljung-Box p-value goes from about 0.03 to about 0.18, so no significant autocorrelation is left. It is the natural upgrade once the seasonality is handled.
hits=tot=0
for o in range(len(s)-150, len(s)-30, 30): # four separate 30-day windows, so no day counts twice
tr=s.iloc[:o]; m=sm.OLS(tr.values, design(tr.index,0)).fit()
pr=m.get_prediction(design(s.index[o:o+30],len(tr))).summary_frame(alpha=0.05)
act=s.iloc[o:o+30].values
hits+=int(((act>=pr['obs_ci_lower'].values)&(act<=pr['obs_ci_upper'].values)).sum()); tot+=len(act)
print('actuals inside the 95%% interval: %d of %d = %.0f%%' % (hits, tot, 100*hits/tot))
print('a 95%% band should hold about %.0f of %d' % (0.95*tot, tot))
actuals inside the 95% interval: 117 of 120 = 98% a 95% band should hold about 114 of 120
A prediction interval is a staffing promise: about 95% of actual days should land inside it. Backtesting the coverage checks that promise. If only 80% fell inside, the grid team would be caught short too often. Here 117 of 120 days landed inside, a little more than the 114 a 95% band would hold and within what chance allows over 120 days, so the band is honest and, if anything, errs slightly wide for reserve planning.
Take-it-further summary, in plain terms¶
- Choose the Fourier order by AIC / holdout: enough harmonics to trace the annual curve, no more.
- Rolling-origin backtesting confirms the accuracy is stable, not one lucky split.
- Each seasonality earns its place: dropping the weekly or yearly term measurably hurts.
- Dynamic regression adds an AR term on the errors, and the small leftover autocorrelation stops being significant.
- Interval coverage of 117 in 120 held-out days says the band is honest, erring slightly wide.