Time-Series ML III — Global / panel forecasting: where ML wins¶
One global model across many series vs a classical model fitted per series (the M5 lesson)¶
The two previous notebooks found classical econometrics winning on one clean series (ex1) and tying shrinkage-based ML on multivariate macro (ex2). This one reaches the regime where machine learning decisively wins — and it is the regime that dominates real quantitative practice: many related series forecast at once.
There are two ways to forecast a panel of series:
- Local — fit a separate classical model (ARIMA, ETS, an AR) to each series in isolation. Each model sees only its own history.
- Global — train one machine-learning model across all series pooled together, sharing parameters. It learns patterns that recur across the panel and borrows strength from every series to help each one.
Where the received wisdom comes from. Two competitions are usually cited, and they are worth separating because they reached different verdicts two years apart.
M4 (2018) ran 100,000 series across yearly, quarterly, monthly, weekly, daily and hourly frequencies, drawn from unrelated business and economic domains. Its headline finding was not a win for machine learning: several pure ML entries finished below a seasonal-naive benchmark. The winner, Slawek Smyl's ES-RNN, was a hybrid — exponential smoothing handling level and seasonality per series, a recurrent network learning the dynamics shared across them — and the runner-up, FFORMA, was a feature-based weighted combination of classical methods. The lesson drawn at the time was that statistical structure and learned components were complements, and that combinations beat individuals.
M5 (2020) ran Walmart's daily sales: roughly 42,840 series in a single coherent hierarchy, with exogenous drivers (prices, promotions, SNAP days, calendar events) and a great deal of intermittent demand. Here the verdict flipped hard — every one of the top 50 accuracy entries used machine learning, and the winner was a LightGBM. Global gradient boosting simply dominated.
Why the reversal? The two panels sit in different regimes, and the difference is precisely the one this notebook measures below. M4's series are individually longer and mutually unrelated, so there is limited common structure to pool and each series can largely fit itself. M5's are short, numerous, related by a shared calendar and shared shoppers, and accompanied by exogenous features — the wide-and-shallow regime where a model that borrows across the panel has an enormous advantage over one that does not.
So "global ML wins" is a claim about a regime, not about model classes in general. The experiment in section 4 identifies where the boundary falls. We reproduce that result on a panel of 48 stocks' volatility, and expose why it happens with a controlled experiment — the global advantage grows precisely as per-series data becomes scarce. Data: the 48-stock weekly returns from the asset-risk arc. Python-lead.
1. The panel and the task¶
For each of the 48 stocks we build a weekly volatility proxy — the log of a 4-week rolling standard deviation of returns — a persistent, forecastable series (volatility clusters). The goal is to forecast every stock's next-week volatility; we score by RMSE averaged across all 48. Crucially, the stocks share structure: volatility rises and falls together across the market, so a model that learns from the whole panel has information a single-stock model does not. That shared structure is what the global approach exploits.
import os
os.environ["KMP_DUPLICATE_LIB_OK"]="TRUE"
import numpy as np, pandas as pd, matplotlib.pyplot as plt, warnings, time
warnings.filterwarnings("ignore")
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; GREY="#a0aec0"; PURP="#6b46c1"
d=pd.read_csv("stocks_weekly.csv",index_col=0); tickers=list(d.columns); dates=pd.to_datetime(d.index)
R=d.values/100.0
V=np.log(pd.DataFrame(R).rolling(4).std().values+1e-6)[4:]; T,Nst=V.shape; vdates=dates[4:]
print(f"panel: {Nst} stocks x {T} weekly volatility observations")
fig,ax=plt.subplots(figsize=(12,4))
for s in range(0,Nst,6): ax.plot(vdates,V[:,s],lw=.6,alpha=.7)
ax.plot(vdates,V.mean(1),color="black",lw=2,label="panel average"); ax.set_ylabel("log rolling volatility")
ax.set_title("48 stock volatility series — they co-move (a global model can pool this)"); ax.legend()
plt.tight_layout(); plt.show()
print("Volatilities rise and fall together (2020 spike visible in all) -- shared structure the global model borrows across series.")
panel: 48 stocks x 308 weekly volatility observations
Volatilities rise and fall together (2020 spike visible in all) -- shared structure the global model borrows across series.
2. Local classical vs one global ML model¶
Local baselines fit one model per stock on that stock's history: a per-series AR($p$) and a per-series ETS (exponential smoothing). The global models train a single learner on the pooled data — every (lagged-features → next-volatility) example from all 48 stocks stacked together: a gradient booster and a global LSTM. All are scored 1-step-ahead over the last 40 weeks, averaged across the 48 stocks. A persistence baseline sets the floor — the forecast that repeats the last observed value, here "next week's log-volatility equals this week's", with nothing estimated at all. It is deliberately the crudest thing on the list: a model that cannot beat it has not extracted anything from the history. Note that it is a low bar only when a series is unpredictable — on a persistent series like this one (the property, not the baseline: volatility clusters, so this week resembles last week) repeating the last value is already a decent forecast, which is why the spread between top and bottom of the table below is narrow.
from statsmodels.tsa.holtwinters import ExponentialSmoothing
from numpy.linalg import lstsq
import xgboost as xgb, torch, torch.nn as nn; torch.set_num_threads(2)
p=8; Hte=40; te0=T-Hte
def feats(v,i): return [v[i-l] for l in range(1,p+1)]+[np.mean(v[i-4:i]),np.mean(v[i-8:i])]
def rmse(a,b): return float(np.sqrt(np.mean((np.array(a)-np.array(b))**2)))
# --- local models + assemble global training set ---
loc_ar=[]; loc_ets=[]; loc_ets_path=[]; persist=[]; GX=[]; GY=[]; tests=[]
for s in range(Nst):
v=V[:,s]
Xl=np.array([feats(v,i) for i in range(p,te0)]); Yl=v[p:te0]
b=lstsq(np.column_stack([np.ones(len(Xl)),Xl]),Yl,rcond=None)[0]
Xt=np.array([feats(v,i) for i in range(te0,T)]); Yt=v[te0:T]
loc_ar.append(rmse(np.column_stack([np.ones(len(Xt)),Xt])@b, Yt))
persist.append(rmse(v[te0-1:T-1], Yt))
# ETS must be scored the same way as everything else. Fitting once and calling .forecast(40) would be a
# single 40-step path, not 40 one-step forecasts -- and an additive trend extrapolated over 40 weeks of a
# mean-reverting log-volatility series is badly wrong. Refit on an expanding window, one step at a time.
_pr=[]
for _i in range(te0,T):
try: _pr.append(float(np.asarray(ExponentialSmoothing(v[:_i],trend=None).fit().forecast(1))[0]))
except Exception: _pr.append(v[_i-1])
loc_ets.append(rmse(_pr,Yt))
# and, for the comparison, the WRONG way: fit once and extrapolate a 40-step path with a trend.
try: _bad=np.asarray(ExponentialSmoothing(v[:te0],trend="add").fit().forecast(Hte),dtype=float)
except Exception: _bad=np.full(Hte,v[te0-1],dtype=float)
loc_ets_path.append(rmse(_bad,Yt))
for i in range(p,te0): GX.append(feats(v,i)); GY.append(v[i])
tests.append((Xt,Yt))
GX,GY=np.array(GX),np.array(GY)
gbm=xgb.XGBRegressor(n_estimators=400,learning_rate=0.05,max_depth=4,verbosity=0).fit(GX,GY)
glob_gbm=[rmse(gbm.predict(Xt),Yt) for Xt,Yt in tests]
# The comparison "local AR vs global GBM" changes TWO things at once: the data the model sees (one series vs
# all 48) and the model class (linear vs boosted trees). To attribute the gain we need the missing cell of the
# 2x2: one LINEAR AR fitted to the pooled panel. Global data, same model class as the local baseline.
bg=lstsq(np.column_stack([np.ones(len(GX)),GX]),GY,rcond=None)[0]
glob_ar=[rmse(np.column_stack([np.ones(len(Xt)),Xt])@bg,Yt) for Xt,Yt in tests]
# global LSTM on pooled windows
Wl=[]; Yl=[]
for s in range(Nst):
v=V[:,s]
for i in range(p,te0): Wl.append(v[i-p:i]); Yl.append(v[i])
mu,sd=np.mean(GY),np.std(GY)
Wt=torch.tensor(((np.array(Wl)-mu)/sd),dtype=torch.float32).unsqueeze(-1); Yt2=torch.tensor(((np.array(Yl)-mu)/sd),dtype=torch.float32).view(-1,1)
torch.manual_seed(0)
class G(nn.Module):
def __init__(s): super().__init__(); s.l=nn.LSTM(1,32,batch_first=True); s.f=nn.Linear(32,1)
def forward(s,x): o,_=s.l(x); return s.f(o[:,-1])
gl=G(); opt=torch.optim.Adam(gl.parameters(),3e-3); lf=nn.MSELoss()
t=time.time()
for e in range(12):
pm=torch.randperm(len(Wt))
for bi in range(0,len(Wt),512): j=pm[bi:bi+512]; opt.zero_grad(); lf(gl(Wt[j]),Yt2[j]).backward(); opt.step()
glob_lstm=[]
for s in range(Nst):
v=V[:,s]; Wtt=np.array([v[i-p:i] for i in range(te0,T)])
with torch.no_grad(): pr=gl(torch.tensor(((Wtt-mu)/sd),dtype=torch.float32).unsqueeze(-1)).numpy().ravel()*sd+mu
glob_lstm.append(rmse(pr,V[te0:T,s]))
res={"persistence":np.mean(persist),"local AR(p)":np.mean(loc_ar),"local ETS (1-step)":np.nanmean(loc_ets),
"GLOBAL AR (linear)":np.mean(glob_ar),"GLOBAL GBM":np.mean(glob_gbm),"GLOBAL LSTM":np.mean(glob_lstm)}
print(f"(global LSTM trained on {len(Wt):,} pooled examples, {time.time()-t:.0f}s)\n")
print("avg 1-step OOS RMSE across 48 stocks (lower=better):")
for k,v in sorted(res.items(),key=lambda z:z[1]): print(f" {k:20s} {v:.4f}")
_gp=np.mean(loc_ar)-np.mean(glob_ar) # pooling, holding the model class fixed
_gn=np.mean(glob_ar)-np.mean(glob_gbm) # flexibility, holding the data fixed
print(f"\nDecomposing the {np.mean(loc_ar)-np.mean(glob_gbm):+.4f} that global GBM gains over local AR:")
print(f" from POOLING (local AR -> global AR ): {_gp:+.4f}")
print(f" from NONLINEARITY (global AR -> global GBM): {_gn:+.4f}")
_E={"local AR":[], "global AR":[], "global GBM":[]}
for s,(Xt,Yt) in enumerate(tests):
v=V[:,s]
Xl=np.array([feats(v,i) for i in range(p,te0)]); Yl=v[p:te0]
b=lstsq(np.column_stack([np.ones(len(Xl)),Xl]),Yl,rcond=None)[0]
_E["local AR"].append((np.column_stack([np.ones(len(Xt)),Xt])@b-Yt)**2)
_E["global AR"].append((np.column_stack([np.ones(len(Xt)),Xt])@bg-Yt)**2)
_E["global GBM"].append((gbm.predict(Xt)-Yt)**2)
def _dm(a,b):
d=np.mean(np.array(_E[a]),axis=0)-np.mean(np.array(_E[b]),axis=0)
x=d-d.mean(); v_=np.mean(x*x)
for k in range(1,5): v_+=2*(1-k/5)*np.mean(x[k:]*x[:-k])
return d.mean()/np.sqrt(v_/len(d))
print("\nDiebold-Mariano (positive = row worse than column; |stat| > 1.96 to matter):")
for a,b in [("local AR","global AR"),("global AR","global GBM"),("local AR","global GBM")]:
_s=_dm(a,b)
print(f" {a:11s} vs {b:11s} {_s:+6.2f} {'significant' if abs(_s)>1.96 else 'indistinguishable'}")
print("\nSo on this panel the win is NOT the borrowing-strength story. Pooling on its own buys nothing measurable")
print("-- one linear AR fitted to all 48 series forecasts no better than 48 separate ones -- while swapping the")
print("linear map for a boosted one, on identical data, accounts for essentially the whole gain. The mechanism")
print("here is FLEXIBILITY, not the panel. Section 4 shows the conditions under which that reverses.")
print(f"\nSame model, same data, two scoring protocols: local ETS scored as 40 one-step\n"
f"forecasts on an expanding window gives {np.mean(loc_ets):.4f}; fitted once and extrapolated as a\n"
f"single {Hte}-step path with a trend it gives {np.mean(loc_ets_path):.4f}, "
f"{np.mean(loc_ets_path)/np.mean(loc_ets)-1:+.0%} worse.")
print("Nothing about the model changed. That gap is the protocol alone, and it is larger than every\n"
"difference between models on this page -- which is why the evaluation design has to be settled\n"
"before any of the comparisons below mean anything.")
(global LSTM trained on 12,480 pooled examples, 1s) avg 1-step OOS RMSE across 48 stocks (lower=better): GLOBAL GBM 0.4189 GLOBAL LSTM 0.4276 local AR(p) 0.4417 GLOBAL AR (linear) 0.4430 local ETS (1-step) 0.4779 persistence 0.4808 Decomposing the +0.0228 that global GBM gains over local AR: from POOLING (local AR -> global AR ): -0.0012 from NONLINEARITY (global AR -> global GBM): +0.0241 Diebold-Mariano (positive = row worse than column; |stat| > 1.96 to matter): local AR vs global AR -0.24 indistinguishable global AR vs global GBM +4.02 significant local AR vs global GBM +3.32 significant So on this panel the win is NOT the borrowing-strength story. Pooling on its own buys nothing measurable -- one linear AR fitted to all 48 series forecasts no better than 48 separate ones -- while swapping the linear map for a boosted one, on identical data, accounts for essentially the whole gain. The mechanism here is FLEXIBILITY, not the panel. Section 4 shows the conditions under which that reverses. Same model, same data, two scoring protocols: local ETS scored as 40 one-step forecasts on an expanding window gives 0.4779; fitted once and extrapolated as a single 40-step path with a trend it gives 0.7566, +58% worse. Nothing about the model changed. That gap is the protocol alone, and it is larger than every difference between models on this page -- which is why the evaluation design has to be settled before any of the comparisons below mean anything.
3. The result — and why global wins¶
The global models win: one gradient booster (and one LSTM) trained across all 48 stocks beats fitting an AR or ETS to each stock separately, averaged over the panel. The usual explanation is borrowing strength — the global model learning common volatility dynamics from 48× more data than any single-stock model sees. The decomposition in the cell above says that is not what happened here: pooling on its own is worth −0.0012, and the whole gain is the switch from a linear map to a boosted one on identical data. The left panel shows the scoreboard; the right splits the per-stock verdict by mechanism, and the pooling-only histogram sits centred on zero while the full global-GBM one does not.
tab=pd.Series(res).sort_values()
fig,ax=plt.subplots(1,2,figsize=(14,4.4))
col=lambda k:(GREY if "persist" in k else GREEN if "local" in k else BLUE)
ax[0].barh(tab.index,tab.values,color=[col(k) for k in tab.index]); ax[0].invert_yaxis(); ax[0].set_xlim(min(tab)*0.98,max(tab)*1.01)
ax[0].set_xlabel("avg 1-step OOS RMSE"); ax[0].set_title("Global ML (blue) beats per-series classical (green)")
for i,v in enumerate(tab.values): ax[0].text(v+0.001,i,f"{v:.3f}",va="center",fontsize=8)
diff=np.array(loc_ar)-np.array(glob_gbm) # >0 means global GBM better for that stock
diff_pool=np.array(loc_ar)-np.array(glob_ar) # >0 means POOLING alone better for that stock
ax[1].hist(diff,bins=20,color=BLUE,alpha=.7,edgecolor="white",label=f"vs global GBM ({int((diff>0).mean()*100)}% of stocks)")
ax[1].hist(diff_pool,bins=20,color=ORANGE,alpha=.7,edgecolor="white",label=f"vs global AR, pooling only ({int((diff_pool>0).mean()*100)}%)")
ax[1].axvline(0,color=RED,lw=2)
ax[1].set_xlabel("local AR RMSE - global RMSE (per stock)"); ax[1].set_ylabel("# stocks")
ax[1].set_title("Per-stock gain, split by mechanism"); ax[1].legend(fontsize=7)
plt.tight_layout(); plt.show()
print(f"Global GBM {res['GLOBAL GBM']:.3f} and LSTM {res['GLOBAL LSTM']:.3f} beat local AR {res['local AR(p)']:.3f} and 1-step ETS {res['local ETS (1-step)']:.3f}.")
print(f"The per-stock histogram shows the global GBM beating the local AR for {int((diff>0).mean()*100)}% of individual stocks, so this")
print("is not an average driven by a handful of series. What it does NOT show is WHY -- for that, compare the")
print(f"global GBM against the pooled LINEAR fit ({res['GLOBAL AR (linear)']:.3f}) rather than against the local one, which is what the")
print("decomposition in the previous cell does and what section 4 traces across data regimes.")
Global GBM 0.419 and LSTM 0.428 beat local AR 0.442 and 1-step ETS 0.478. The per-stock histogram shows the global GBM beating the local AR for 79% of individual stocks, so this is not an average driven by a handful of series. What it does NOT show is WHY -- for that, compare the global GBM against the pooled LINEAR fit (0.443) rather than against the local one, which is what the decomposition in the previous cell does and what section 4 traces across data regimes.
4. The mechanism — the global edge grows as data gets scarce¶
The clinching experiment. We vary how much per-series history each model gets and re-run local-AR vs global-GBM. When each series is data-rich, the local model has enough to fit itself and the gap is small; as per-series history shrinks, the local models starve while the global model — pooling across all 48 series — barely notices, and its advantage widens. This is exactly the M5 finding and the reason global ML dominates real forecasting problems, which are typically wide (many series) and shallow (short each): retail SKUs, sensor fleets, and cross-sectional finance all live in that regime.
def sweep(train_len):
"""Returns (local AR, global AR, global GBM) so pooling and flexibility can be read separately."""
loc=[]; GXt=[]; GYt=[]; tst=[]
for s in range(Nst):
v=V[:,s]; tr0=max(p,te0-train_len)
Xl=np.array([feats(v,i) for i in range(tr0,te0)]); Yl=v[tr0:te0]
b=lstsq(np.column_stack([np.ones(len(Xl)),Xl]),Yl,rcond=None)[0]
Xt=np.array([feats(v,i) for i in range(te0,T)]); Yt=v[te0:T]
loc.append(rmse(np.column_stack([np.ones(len(Xt)),Xt])@b,Yt))
for i in range(tr0,te0): GXt.append(feats(v,i)); GYt.append(v[i])
tst.append((Xt,Yt))
GXt,GYt=np.array(GXt),np.array(GYt)
bgl=lstsq(np.column_stack([np.ones(len(GXt)),GXt]),GYt,rcond=None)[0]
ga=np.mean([rmse(np.column_stack([np.ones(len(Xt)),Xt])@bgl,Yt) for Xt,Yt in tst])
g=xgb.XGBRegressor(n_estimators=300,learning_rate=0.05,max_depth=4,verbosity=0).fit(GXt,GYt)
return np.mean(loc), ga, np.mean([rmse(g.predict(Xt),Yt) for Xt,Yt in tst])
lens=[20,30,50,90,150,240]; L=[];GA=[];Gg=[]
for tl in lens: lo,ga,gg=sweep(tl); L.append(lo); GA.append(ga); Gg.append(gg)
fig,ax=plt.subplots(1,2,figsize=(13,4.2))
ax[0].plot(lens,L,"o-",color=GREEN,lw=2,label="local AR (48 separate fits)")
ax[0].plot(lens,GA,"o-",color=ORANGE,lw=2,label="global AR (one pooled linear fit)")
ax[0].plot(lens,Gg,"o-",color=BLUE,lw=2,label="global GBM")
ax[0].set_xlabel("per-series training length (weeks)"); ax[0].set_ylabel("avg OOS RMSE")
ax[0].set_title("Three models, two mechanisms"); ax[0].legend(fontsize=8)
ax[1].plot(lens,np.array(L)-np.array(GA),"o-",color=ORANGE,lw=2,label="pooling (local AR - global AR)")
ax[1].plot(lens,np.array(GA)-np.array(Gg),"o-",color=BLUE,lw=2,label="flexibility (global AR - global GBM)")
ax[1].axhline(0,color="k",lw=.5)
ax[1].set_xlabel("per-series training length (weeks)"); ax[1].set_ylabel("RMSE gain")
ax[1].set_title("The two effects trade places"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()
print(f" {'weeks':>7} {'local AR':>10} {'global AR':>11} {'global GBM':>12} {'pooling':>10} {'flexibility':>13}")
for i,tl in enumerate(lens):
print(f" {tl:>7} {L[i]:>10.4f} {GA[i]:>11.4f} {Gg[i]:>12.4f} {L[i]-GA[i]:>+10.4f} {GA[i]-Gg[i]:>+13.4f}")
print(f"\nThis is the experiment that rescues the M5 story, and it rescues it as a CONDITIONAL claim. With {lens[0]} weeks")
print(f"per series, pooling is worth {L[0]-GA[0]:+.4f} and flexibility {GA[0]-Gg[0]:+.4f} -- the panel is doing all the work and the")
print(f"boosted model is, if anything, a slight liability. With {lens[-1]} weeks the position is reversed: pooling is worth")
print(f"{L[-1]-GA[-1]:+.4f} and flexibility {GA[-1]-Gg[-1]:+.4f}. The two curves cross somewhere around 50-90 weeks.")
print("\nSo borrowing strength across a panel is real, and it is exactly what M5 found -- but it is a remedy for SHORT")
print("series. This panel has 260 weeks each, which puts it on the wrong side of the crossover, and that is why the")
print("headline comparison above is won by flexibility rather than by pooling. The right summary is not 'global beats")
print("local' but 'pooling pays when series are short, flexibility pays when they are long' -- and knowing which")
print("regime you are in is the whole decision.")
weeks local AR global AR global GBM pooling flexibility
20 0.5842 0.4472 0.4527 +0.1370 -0.0055
30 0.5071 0.4446 0.4519 +0.0625 -0.0072
50 0.4747 0.4439 0.4379 +0.0308 +0.0060
90 0.4518 0.4442 0.4236 +0.0075 +0.0207
150 0.4420 0.4431 0.4200 -0.0011 +0.0231
240 0.4420 0.4429 0.4172 -0.0009 +0.0257
This is the experiment that rescues the M5 story, and it rescues it as a CONDITIONAL claim. With 20 weeks
per series, pooling is worth +0.1370 and flexibility -0.0055 -- the panel is doing all the work and the
boosted model is, if anything, a slight liability. With 240 weeks the position is reversed: pooling is worth
-0.0009 and flexibility +0.0257. The two curves cross somewhere around 50-90 weeks.
So borrowing strength across a panel is real, and it is exactly what M5 found -- but it is a remedy for SHORT
series. This panel has 260 weeks each, which puts it on the wrong side of the crossover, and that is why the
headline comparison above is won by flexibility rather than by pooling. The right summary is not 'global beats
local' but 'pooling pays when series are short, flexibility pays when they are long' -- and knowing which
regime you are in is the whole decision.
5. Summary — and the subsection's arc¶
On a wide panel of series, one global machine-learning model beats a classical model fitted per series — the M4/M5 result, reproduced here on 48 stocks' volatility. Global gradient boosting and a global LSTM, trained across the whole panel, beat per-series AR and ETS on average and for the large majority of individual stocks. The reason is not the one usually given: on this panel, pooling buys nothing measurable and the gain is flexibility. The controlled experiment explains why, and rescues the M5 claim as a conditional one — the pooling edge is largest when per-series data is scarce and vanishes as each series grows, and 260 weeks each puts this panel past the crossover — so global ML dominates precisely the wide-and-shallow problems (many short series) that define real forecasting at scale.
This completes the through-line of the Time-Series ML subsection — the honest, conditional answer to "does ML beat traditional econometrics on time series?":
| setting | winner |
|---|---|
| ex1 one clean series | classical (SARIMA/ETS) dominates |
| ex2 multivariate macro | shrinkage wins; ML ties |
| ex3 wide panel of series | global ML wins |
The determinant is not "ML vs econometrics" in the abstract but the shape of the problem: single and clean favours purpose-built classical models; wide and shallow favours global ML; in between, shrinkage (Bayesian or L1) is what matters. That judgement — match the method to the structure of the forecasting problem — is the subsection's payoff, and the same lesson the whole ML arc has taught in every domain. Next (ex4): the hardest case — financial returns, where the signal is faint and honesty about predictability is the whole game.