Hedging Memory Horizons for Non-Stationary Prediction via Online Aggregation
arXiv:2605.06541v1 [cs.LG] 7 May 2026
Yutong Wang Department of Statistics, LSE London, UK [email protected]
Yannig Goude Électricité de France R&D Palaiseau, France [email protected]
Qiwei Yao Department of Statistics, LSE London, UK [email protected]
Abstract We study online prediction under distribution shift, where inputs arrive chronologically and outcomes are revealed only after prediction. In this setting, predictors must remain stable in quiet regimes yet adapt when regimes shift, and the right adaptation memory is unknown in advance. We propose MELO (Memory-hedged Exponentially Weighted Least-Squares Online aggregation), a model-agnostic method that hedges across adaptation scales: it wraps any non-anticipating base-predictor pool with exponentially weighted least-squares (EWLS) adaptation experts at multiple forgetting factors, and aggregates raw and EWLS-adapted forecasts with MLpol which is a parameter-free online aggregation rule. Under boundedness conditions, we establish deterministic oracle inequalities showing that it competes with both the best raw predictor and the best bounded, time-varying affine combinations of the base predictions, up to a path-length-dependent tracking cost and a sublinear aggregation overhead. We evaluate MELO on French national electricity-load forecasting through the COVID-19 lockdown using no regime indicators, lockdown dates, or policy covariates. MELO reduces overall RMSE by 34.7% relative to base-only MLpol and achieves lower overall RMSE than a TabICL reference supplied with an external COVID policy-response covariate. MELO requires only lightweight per-step recursive updates without model retraining.
1
Introduction
Operational forecasters can fail abruptly when reality changes faster than their update mechanism. The COVID-19 lockdown sharply disrupted French national electricity load within days. Models calibrated on years of pre-pandemic data incurred large short-run errors, while reliable retraining required enough post-shift labels and therefore could not provide an immediate correction. This pattern recurs across deployed forecasting settings, including short-term energy load, retail demand, financial quantities, and user behaviour. In these settings, inputs arrive chronologically and models are evaluated on future periods; the deployment distribution may also shift after training. Throughout, causal means non-anticipating in time: the prediction at time t may use only inputs and labels revealed before t, and not future outcomes. The appropriate adaptation memory is therefore both unknown in advance and itself non-stationary: a horizon well tuned to one regime can be badly miscalibrated for the next. Existing approaches address only part of this trade-off. Offline-trained models, such as tree ensembles [Chen and Guestrin, 2016, Ke et al., 2017], GAMs [Hastie and Tibshirani, 1986, Wood, 2017], deep tabular networks [Gorishniy et al., 2021], and tabular foundation models [Hollmann et al., 2023, Qu et al., 2025], encode rich historical structure but cannot adapt without explicit retraining or reconditioning, which requires labels not yet available immediately after a regime shift. Adaptive filters update online, including recursive least-squares methods [Hayes, 1996, Sayed, 2003], and Preprint.
Kalman-style forecast corrections for load forecasting [Obst et al., 2021, De Vilmarest and Goude, 2022]. However, a single adaptation memory typically governs their behaviour: long memory is stable but slow at breaks, whereas short memory adapts quickly but is noisy in quiet periods. Adaptive Kalman variants [Huang et al., 2018, de Vilmarest and Wintenberger, 2024] attempt to infer adaptation scales online, but each does so within its own state-space parametrisation and update mechanism. Classical forecast combination and online expert aggregation [Bates and Granger, 1969, Stock and Watson, 2004, Cesa-Bianchi and Lugosi, 2006, Gaillard et al., 2014, Devaine et al., 2013] weight heterogeneous forecasters by realised loss. This gives robustness to the available expert pool, but the aggregate cannot move outside their convex hull with simplex weights over raw forecasts. A shared post-shift bias is therefore inherited rather than corrected. We propose MELO (Memory-hedged Exponentially Weighted Least-Squares Online aggregation), a causal adaptation layer for online prediction under deployment-time shift. Rather than committing to a single memory horizon, MELO exposes several horizons as competing experts and lets realised forecasting loss decide how much to trust each. Concretely, the entire base-prediction vector is fed into a time-varying affine combination layer. Unlike convex aggregation over the raw forecasts, this layer can move outside their convex hull when post-shift bias must be corrected. It maintains online corrections at multiple adaptation speeds in parallel, while keeping the raw forecasts as conservative fallbacks. The choice of which correction speed to trust is then deferred to a loss-driven outer aggregation rule. The resulting procedure requires no regime labels, change-point annotations, or base-model retraining. Contributions. • We introduce MELO, a model-agnostic causal adaptation layer motivated by residual-aware correction. It treats base forecasts as nonlinear representations with potentially complementary errors, wraps them with multi-scale EWLS experts, and aggregates raw and corrected forecasts with MLpol using lightweight recursive updates. • We prove deterministic oracle inequalities, under explicit boundedness/stability conditions, showing that MELO competes with both the best raw predictor and the best correction memory scale in the EWLS grid, balancing prediction robustness with the cost of adapting to a changing affine comparator. • On COVID-era French electricity-load forecasting, MELO outperforms base-only aggregation and matched adaptive-filter baselines, and also beats a stronger-information TabICL+GRI reference on overall RMSE.
2
Related Work
Online aggregation and prediction with expert advice. Prediction with expert advice gives an online-learning perspective on classical forecast combination [Bates and Granger, 1969, Timmermann, 2006, Cesa-Bianchi and Lugosi, 2006]: a forecaster combines a finite pool of experts sequentially and is evaluated by regret against hindsight benchmarks, from the best fixed expert [Cesa-Bianchi and Lugosi, 2006], to shifting experts [Herbster and Warmuth, 1998, Bousquet and Warmuth, 2002] and, more broadly, a time-varying comparator sequence in dynamic regret [Zinkevich, 2003, Cesa-Bianchi et al., 2012, Zhang et al., 2018]. Second-order aggregation methods such as BOA [Wintenberger, 2017] and MLpol [Gaillard et al., 2014] use data-dependent updates that adapt to the realised loss sequence; we use MLpol as a parameter-free aggregation rule. Conceptually closer to our multiscale design, Bousquet and Warmuth [2002] introduce mixing past posteriors as an implicit form of memory-scale hedging, and Cesa-Bianchi et al. [2012] provide a unified mirror-descent analysis of fixed-share, weight-sharing, and related tracking schemes; MELO makes the memory hedging explicit by exposing a grid of EWLS adaptation horizons as competing experts. Online expert aggregation has also been applied to financial nonstationarity [Miao and Polak, 2023], where a multiplicative-weights aggregator over a heterogeneous model pool (linear, tree, neural) achieves robust performance through COVID; MELO differs in hedging across memory horizons of the same learner family rather than across model classes, and in maintaining recursive expert state rather than refitting on rolling windows. Within aggregation theory, Q-aggregation [Dai et al., 2012, 2014] and ϕ-divergence FTRL [Alquier, 2021] sharpen the deviation behaviour and unboundedness handling of EWA in the batch and unbounded-loss settings respectively; MLpol’s polynomial potential plays 2
an analogous role online. Our contribution lies at the level of the expert pool: we augment raw predictors with a multi-scale family of online residual-correction experts, so the aggregate can choose, by realised squared loss, between stable historical predictors and rapidly adapting combinations. Tabular prediction under temporal distribution shift. Modern tabular prediction is dominated by strong tree ensembles [Chen and Guestrin, 2016, Ke et al., 2017], with growing interest in deep tabular models and tabular foundation models [Gorishniy et al., 2021, Hollmann et al., 2023, Qu et al., 2025]. Much of this literature evaluates on random splits, where training and test samples are approximately exchangeable; the recent chronological benchmark TabReD [Rubachev et al., 2025] shows that this can hide large performance gaps under temporal distribution shift. Our goal is orthogonal to tabular architecture design: any causal tabular predictor can be used as a base expert, while the EWLS expert learns online combinations from realised forecast errors. Concept drift and deployment-time adaptation. Streaming and concept-drift methods address non-stationarity through sliding windows, drift detection, online ensembles, or adaptive retraining [Gama et al., 2014, Bifet and Gavaldà, 2007, Gomes et al., 2017]. Test-time and deployment-time adaptation methods modify model statistics, representations, or predictions during inference, with or without labels [Kim et al., 2025, Liu et al., 2025, Lee et al., 2025]. We study the label-revealedafter-prediction setting, where adaptation can be driven directly by realised forecast errors. The proposed EWLS layer therefore provides a model-agnostic supervised adaptation mechanism, without committing to a single drift detector, retraining rule, or hand-chosen window schedule. Unlike testtime adaptation methods that modify base-model internals, MELO leaves base predictors untouched and operates entirely as an outer aggregation layer. Adaptive filtering and adaptive Kalman methods. EWLS/RLS adaptive filters track time-varying linear relationships through least-squares recursions [Hayes, 1996, Sayed, 2003], while Kalman filtering gives the classical recursive prediction–correction framework for linear state-space models [Kalman, 1960]. Adaptive Kalman-style variants can estimate process uncertainty or adaptation scales online [Huang et al., 2018, de Vilmarest and Wintenberger, 2024], and related state-space corrections have also been developed for electricity-load forecasting through Kalman-style updates of forecasting-model coefficients [Obst et al., 2021, De Vilmarest and Goude, 2022]. These approaches are usually tied to a particular state-space or model-update parametrisation. We instead combine the entire base-prediction pool, run several fixed-scale filters as experts in parallel, and use an outer loss-driven aggregator to hedge the unknown adaptation time scale.
3
Methodology Base forecasts
Expert pool
Online aggregation
Raw forecast experts (1)
(M )
zt = ft , . . . , ft
XGB, GAM, . . . static in main experiment
M pass-through base forecasts
EWLS correction experts
available information
MLpol weights over expert pool
ŷt
EWLS (γ1 ) EWLS (γ2 )
Ht−1
updated online
Raw forecast experts
update MLpol weights
...
yt
EWLS (γK ) updated online
observed target update EWLS states
Figure 1: Overview. Causal base predictions zt enter the expert pool directly and through K EWLS experts with different forgetting factors. MLpol aggregates the M + K raw and adapted candidates into ŷt ; after yt is observed, both EWLS states and MLpol weights are updated online. We study sequential one-step-ahead prediction. Let Ht−1 denote all information available to the learner before yt is revealed. By convention, Ht−1 subsumes the entire past (previous covariates, 3
predictions, and revealed outcomes) as well as any current period inputs known prior to predicting yt , such as calendar variables, lagged outcomes, or exogenous covariates whose values for time t are already observed or forecast. At round t, the learner uses Ht−1 to form a prediction ŷt , then observes yt and incurs the one-step squared loss ℓt (ŷt ) = (ŷt − yt )2 . For the prediction sequence ŷ = (ŷ1 , . . . , ŷT ), the objective is to minimize the cumulative loss T T X X LT (ŷ) = ℓt (ŷt ) = (ŷt − yt )2 . (1) t=1
t=1
(m)
Let zt,m = ft (Ht−1 ) be the prediction produced by base forecaster m ∈ [M ], and collect the base predictions as zt ≡ (zt,1 , . . . , zt,M )⊤ ∈ RM . This formulation is deliberately general: a (m) rolling-retrain model corresponds to an ft that is refitted each round using training data contained (m) in Ht−1 ; a static model corresponds to a fixed map ft ≡ f (m) evaluated on the current-period inputs encoded in Ht−1 ; and physics-based or rule-based forecasters arise as deterministic mappings of Ht−1 . The base forecasters are treated as black-box causal prediction streams. We only require their forecasts to be available before yt . They may be fixed or causally updated, provided the expert identities remain coherent across rounds. In the main experiments they are trained once and kept fixed, isolating MELO’s adaptation from any benefit of base retraining. Step 1. EWLS experts. Put z̃t = (zt⊤ , 1)⊤ . For a forgetting factor γ ∈ [1/2, 1), the EWLS expert (γ) (γ) gives a prediction y̌t = z̃t⊤ wt , where the coefficient vector is the exponentially weighted ridge solution based only on past observations: t−1 nX o (γ) 2 wt ∈ arg min γ (t−1)−s (ys − w⊤ z̃s )2 + γ t−1 δ∥w∥2 . (2) w
s=1
Here δ > 0 is a small ridge ensuring the discounted Gram matrix is invertible. Thus each EWLS expert learns a time-varying affine combination of the base predictions, with an intercept. The constant 1 appended to z gives each EWLS expert a time-varying intercept, allowing it to absorb level shifts in yt , such as a sudden drop in baseline demand, without forcing the adaptation through the base coefficients. The forgetting factor controls the adaptation scale via nominal memory length h(γ) = 1/(1 − γ): smaller h(γ) reacts faster but is noisier, whereas larger h(γ) is more stable but adapts more slowly. Since the appropriate adaptation scale is unknown and may itself change over time, MELO maintains K EWLS experts with forgetting factors Γ = {γ1 , . . . , γK }. Step 2. MLpol aggregation over base and EWLS experts. At time t, MELO forms a pool of N = M + K candidate predictions by concatenating the raw base predictors and the EWLS experts: (γ ) (γ ) ⊤ ỹt ≡ (ỹt,1 , · · · , ỹt,N )⊤ = zt,1 , · · · , zt,M , y̌t 1 , · · · , y̌t K .
The final predictor is MLpol aggregation [Gaillard et al., 2014] ŷt = p⊤ t ỹt . MLpol assigns weights according to positive cumulative pseudo-regrets: t−1 X [Rt−1,j ]+ pt,j = PN , Rt−1,j = 2(ŷs − ys )(ŷs − ỹs,j ), (3) s=1 i=1 [Rt−1,i ]+
where [x]+ = max{x, 0}. If all positive parts vanish, uniform weights apply. This is the gradienttrick form for squared loss: an expert receives positive pseudo-regret when it would have moved the aggregate prediction toward the observed target. Raw base predictors are included directly in the pool, so the aggregate can fall back to them when online correction is unnecessary. Each candidate prediction is an affine function of z̃t , so the aggregate can equivalently be analysed against time-varying affine combinations. Grid design. We place forgetting factors on a geometric grid in nominal memory length, hk = 1/(1 − γk ), over [hmin , hmax ]. This uses only K = O(log(hmax /hmin )) correction experts while ensuring that, if the oracle memory length lies in the range, some grid point approximates it up to a constant factor. Section 4 makes this precise. 4
Implementation and Complexity. EWLS solutions are updated by the standard recursive leastsquares (RLS) recursion, with causal initialization and a small covariance inflation for numerical stability, see Algorithm 1 in Appendix B.1. Let d = M + 1. Each EWLS expert costs O(d2 ) per step via the RLS recursion, giving total per-step cost O(Kd2 + N ). This is independent of T and of the base predictors’ own training cost.
4
Theoretical guarantee
We summarise the main deterministic, sequence-wise guarantees that motivate the two design choices: maintaining several EWLS time scales, and aggregating raw base forecasts together with their online corrections. The statements in this section are informal main-text versions; the formal oracle inequalities, including constants, clipping conventions, and proofs, are given in Appendix A. PT Setup. Write LT (a1:T ) = t=1 (at − yt )2 for any prediction sequence a1:T = (at )Tt=1 . For a time-varying affine combination path u1:T = (u1 , . . . , uT ), define LT (u1:T ) =
T X t=1
2 z̃t⊤ ut − yt ,
PT (u1:T ) =
T X t=2
∥ut − ut−1 ∥2 .
The path length PT measures the total movement of the comparator path and accommodates both gradual drift and abrupt changes. Write UT (R) = {u1:T : ∥ut ∥2 ≤ R for all t}. Assumptions. We assume bounded data, ∥z̃t ∥2 ≤ Bz and |yt | ≤ By , and bounded EWLS iterates along the realised sequence: for the EWLS expert with forgetting factor γ, the unprojected forgotten(γ) RLS coefficients satisfy ∥wt ∥2 ≤ R for t ∈ [T ], applied to the relevant EWLS expert(s) in each statement below. This is an explicit stability condition on the realised recursion, not a consequence of forgetting alone. (γ)
(γ)
Single-scale tracking. Let y̌1:T = (y̌t )Tt=1 denote the EWLS prediction sequence generated by the idealised recursion in (2). Denote hγ = 1/(1 − γ). Theorem 1 (Single-scale tracking; informal, see Theorem A.3). Under the assumptions above, for any γ ∈ [1/2, 1) and any comparator path u1:T ∈ UT (R), C1 dT (γ) LT y̌1:T ≤ LT (u1:T ) + + C2 hγ PT (u1:T ) + l.o.t., hγ where C1 , C2 hide only boundedness and initialisation constants, and l.o.t. collects δ-initialisation and logarithmic-in-h terms. The two leading terms display a stability–adaptivity trade-off: short memory reduces the tracking term hγ PT but inflates the finite-memoryp term dT /hγ ; long memory does the opposite. When PT (u1:T ) > 0, the balancing scale is h⋆ ≍ dT /PT (u1:T ), which depends on the unknown amount of comparator movement; when PT (u1:T ) = 0, the bound reduces to LT (u1:T ) + C1 dT /hγ , which favours the longest available memory. Multi-scale aggregation. MELO sidesteps the choice of h by exposing a grid of EWLS scales as competing experts and letting MLpol decide by realised loss. Let Γ ⊂ [1/2, 1) be a finite-forgetting (j) grid with |Γ| = K, and let EΓ = {f (1) , . . . , f (M ) } ∪ {y̌ (γ) : γ ∈ Γ}. Let f1:T denote the raw prediction sequence of base expert j and let ŷ1:T denote the resulting MLpol aggregate over EΓ . Theorem 2 (MELO oracle inequality; informal, see Theorem A.6). Under the same assumptions and the bounded-loss MLpol convention formalised in Appendix A, suppressing logarithmic and initialisation terms, the MLpol aggregate over EΓ satisfies, for every comparator path u1:T ∈ UT (R), p (j) C1 dT LT (ŷ1:T ) ≤ min min LT (f1:T ), LT (u1:T ) + min hγ + C2 hγ PT (u1:T ) + C3 (M + K)T , γ∈Γ 1≤j≤M {z } {z } | | best EWLS scale tracking u1:T
best raw predictor
where C1 , C2 , C3 hide only boundedness and initialisation constants. 5
The aggregate competes with two benchmarks simultaneously. The first is a safety comparison: the bounded-loss MLpol aggregate over the enlarged pool remains competitive with the best individual raw predictor up to the sublinear aggregation overhead. The second captures distribution shifts through the best online affine correction at the best memory scale in the grid. Geometric grids. A geometric grid over [hmin , hmax ] with ratio ρ contains a point within factor ρ of h⋆ whenever h⋆ lies in the range. Therefore, up to the p suppressed logarithmic and initialisation terms, the inner minimum in Theorem 2 is of order Cρ dT PT (u1:T ) when PT (u1:T ) > 0, using only K = O(log(hmax /hmin )) experts. The precise corollary is in Appendix A.
5
Forecasting French daily electricity load
Task and protocol. We evaluate one-step-ahead forecasting of French national electricity load, a standard operational short-term load forecasting task [Devaine et al., 2013, Obst et al., 2021, De Vilmarest and Goude, 2022]. We use the public RTE (French TSO) daily load dataset from 2012-01-01 to 2021-01-15 and refer to this benchmark as RTE-FR Load. Base models are trained on 2012–2018. Online methods are then evaluated sequentially from 2019-01-01 to 2021-01-15 (T = 746 days), with each update using only true values available after the corresponding prediction. We report RMSE overall and over three regimes: pre-lockdown (2019-01-01 to 2020-03-16, 441 days), lockdown (2020-03-17 to 2020-05-11, 56 days), and the post-first-lockdown period (2020-05-12 to 2021-01-15, 249 days). Models and tuning. The base pool contains M = 7 heterogeneous forecasters: Lag-1, Linear(Ridge) [Hoerl and Kennard, 1970], XGBoost [Chen and Guestrin, 2016], LightGBM-GBDT [Ke et al., 2017], GAM [Hastie and Tibshirani, 1986, Wood, 2017], ResNet [He et al., 2016, Gorishniy et al., 2021], and FT-Transformer [Gorishniy et al., 2021]. All are trained once on 2012–2018 data with a shared feature set including lagged load, temperature, calendar, and holiday variables. To isolate algorithmic adaptation from weather-forecast quality, all methods use realised-temperature covariates in this controlled benchmark; an operational deployment would replace them by temperature forecasts, see Appendix B.2. Hyperparameters use published defaults or tail-of-training validation, with no tuning on the test period, see Appendix B.3. The EWLS layer uses K = 16 correction experts: 15 finite-forgetting factors are placed on a geometric grid over γ ∈ [0.950, 0.9998], equivalently nominal scales h(γ) = 1/(1 − γ) ∈ [20, 5000], and one no-forgetting endpoint γ = 1 is appended. The theory applies to the finite-forgetting experts; the γ = 1 endpoint is included empirically as a static long-memory reference. Numericalstability parameters, including the covariance-inflation scale ε0 = 10−8 , are selected by walk-forward validation on 2018 only, see Appendices B.4 and B.5. Comparisons. Our main ablations are MLpol on the base pool only, MLpol on the EWLS pool only, and MELO, which is MLpol on the combined Base+EWLS pool. We also report the best single EWLS expert selected in hindsight on overall test RMSE as a within-method oracle baseline. External references include individual static baselines, tabular foundation-model baselines run daily with an expanding in-context training window (TabPFN and TabICL), stronger-information variants augmented with the Oxford Government Response Index (GRI; Hale et al. 2021), and adaptive Kalman-filter alternatives tuned under the same 2018 walk-forward protocol. 5.1
Main performance and external baselines
Table 1 shows the main empirical message: adding the EWLS experts to the raw base-forecast pool substantially improves online aggregation under distribution shift. Relative to MLpol on the base pool alone, MELO reduces RMSE before lockdown (+9.8%), by the largest margin during the COVID-19 regime break (+55.7%), and still substantially after the first lockdown period (+36.1%), yielding a +34.7% overall improvement. A paired moving-block bootstrap on the daily squared-loss sequences confirms that the overall margins over both base-only MLpol and TabICL(+GRI) are significant at the 5% level; details are in Appendix C.2. The ablations separate the two roles of the Base+EWLS pool. Base-only aggregation is reasonably competitive before lockdown, but deteriorates sharply under the break. EWLS-only aggregation 6
Table 1: Per-regime and overall RMSE (MW) on the French electricity-load test period. Bold: best per column among non-hindsight methods; underline: second best among non-hindsight methods. Hindsight rows (italicised) are non-deployable EWLS reference benchmarks and do not compete for highlighting. “Best single EWLS (hindsight)” selects one forgetting factor using overall test RMSE. “Per-regime best EWLS (hindsight)” selects the in-segment best forgetting factor separately for each regime (γ = 0.9990, 0.9773, 0.9997), using regime boundaries and future losses, and is therefore not realisable a priori. The expanded table with all individual base baselines, all single EWLS scales, and additional external references is in Appendix C. “+GRI” denotes an additional regime-informed covariate; MELO uses no such signal. Method Regime length (days)
Pre-lockdown 441
Lockdown 56
Post-lockdown 249
Overall 746
GAM (best static overall) TabICL (online) TabICL (online, +GRI)
946.6 670.3 668.0
3201.3 2943.6 1132.6
1076.2 687.4 699.5
1298.3 1036.2 723.4
MLpol on Base only Best single EWLS (γ = 0.9986, hindsight) Per-regime best EWLS (hindsight) MLpol on EWLS only MELO (Base+EWLS)
690.5 657.9 656.2 653.9 623.1
2452.7 1168.3 998.8 1073.6 1086.1
907.0 585.4 573.6 584.1 579.3
1004.0 687.5 662.8 673.2 655.8
RMSE reduction vs. Base only (MW) RMSE reduction vs. Base only (%)
+67.4 +9.8%
+1366.6 +55.7%
+327.8 +36.1%
+348.3 +34.7%
Table 2: Comparison with adaptive Kalman-filter alternatives under the same tuning protocol. All nuisance hyperparameters are selected on 2018 out-of-sample data; the 2019–2021 test period is never used for selection. Each adaptive-KF layer is combined with the same direct base experts by the same outer MLpol aggregator. Full diagnostics and filter-only ablations are in Appendix D. Method MELO (Base+EWLS) MLpol on Base+VIKING MLpol on Base+IMM-KF MLpol on Base+VBAKF
Pre-lockdown
Lockdown
Post-lockdown
Overall
623.1 634.0 636.5 665.1
1086.1 1185.3 1276.7 1227.4
579.3 599.7 648.6 665.4
655.8 680.5 708.7 722.8
adapts especially well during the short lockdown segment, where it slightly outperforms the full pool, but the raw base forecasts remain useful in quieter periods. The full Base+EWLS aggregation therefore gives the best overall balance: it is the best non-hindsight method before lockdown, after lockdown, and overall. Within the EWLS family alone, no single forgetting factor matches MELO’s overall RMSE. Even a non-deployable per-regime EWLS oracle remains slightly worse overall. The comparison with TabICL is a stress test rather than a claim that the EWLS correction layer replaces foundation-model forecasts. The fair same-information comparison is TabICL (online), which uses the same causal information available to MELO. TabICL (online, +GRI) is instead a stronger-information reference, since the Oxford Government Response Index changes contemporaneously with the policy regime. TabICL (online) is substantially worse than MELO, whereas TabICL (online, +GRI) narrows the gap by adding this regime-informed covariate. MELO still achieves lower overall RMSE than this stronger-information reference while using no such exogenous regime signal. Finally, the framework is extensible rather than tied to a particular base pool. Adding the TabICL (online, +GRI) forecast as an additional causal expert improves MELO from 655.8 to 621.1 RMSE, and including the same forecast in the EWLS design vector further improves RMSE to 601.9, see Appendix H. Adaptive-filter alternatives. A natural alternative to a fixed grid of correction scales is a filter that estimates its own adaptation scale online. We compare against three adaptive Kalman-filter baselines with complementary process-uncertainty parameterisations: VIKING [de Vilmarest and Wintenberger, 2024] (isotropic log-variance random walk), VBAKF [Huang et al., 2018] (variational 7
covariance/noise adaptation), and IMM-KF [Blom and Bar-Shalom, 1988] (Bayesian mixing over a fixed grid of Q-scales). VIKING and VBAKF are run as small pools over responsiveness settings and aggregated with the base experts by MLpol; IMM-KF uses its internal Markov mixture over fixed-Q modes and is then also exposed to the same outer MLpol comparison with the direct base experts. All nuisance hyperparameters are selected by the same 2018 walk-forward protocol used for ε0 ; the 2019–2021 test period is never used for selection. Among the matched Base+adaptive-filter pools in Table 2, MELO achieves the lowest RMSE in each of the four reported regimes. The pairwise differences in overall-test-period RMSE are significant under the paired moving-block bootstrap, see Appendix D, which also reports the full bootstrap table, filter-only ablations, and diagnostics helping explain why these adaptive-Q mechanisms do not match the explicit multi-scale EWLS + regret-aggregation separation on this benchmark. 5.2
Adaptation diagnostics Forgetting factor γ
Excess RMSE vs. best fixed scale (%)
0.95
40
0.97
0.98
0.99
0.995
0.999
0.9995
Pre-lockdown Post-lockdown Lockdown
30
20
10
0 102
103
Nominal EWLS scale h(γ) = 1/(1 − γ)
Figure 2: No single forgetting scale is optimal across regimes. Each curve shows the excess RMSE of a standalone EWLS correction with fixed forgetting factor γ, relative to the best fixed scale within that regime. The horizontal axis is the nominal scale h(γ) = 1/(1 − γ), used as an ordered forgetting-scale index because the implemented experts also include small covariance inflation. Empty markers indicate the best fixed scale in each regime: lockdown favours a much shorter scale than the non-lockdown regimes. Figure 2 isolates this scale dependence: lockdown is minimized at much shorter nominal scales than the surrounding regimes, consistent with the trade-off in Theorem 1. Figure 3 shows that MLpol exploits this online: at lockdown onset, weight reallocates sharply into the Fast EWLS bucket, and away from the base experts, so MELO ends at cumulative regret −28 × 106 MW2 relative to the best static convex combination fitted in hindsight. Base-only MLpol diverges to +403 × 106 MW2 over the same window. On two TabReD datasets with much higher base-residual correlation the same framework yields gains below 0.5%, see Appendix E, suggesting a boundary condition: when base residuals are nearly collinear, affine combination has little additional signal to exploit.
8
Stacked MLpol weight
Cumulative regret (×106 MW2)
400
+403
COVID-19 lockdown
MELO (ours) Base-only MLpol Best single EWLS in hindsight (γ=0.99856)
300
Hindsight static optimum (below this line = better)
200
100
+3 -28
0
1.0 0.8
Fast EWLS (h < 100)
Slow EWLS (h ≥ 1000)
Medium EWLS (100 ≤ h < 1000)
Base experts (total)
0.6 0.4 0.2 0.0 Jan 2019
Apr 2019
Jul 2019
Oct 2019
Jan 2020
Apr 2020
Jul 2020
Oct 2020
Jan 2021
Date
Figure 3: Cumulative regret and MLpol weights over the test period. Top: Cumulative excess squared loss of three aggregators relative to the best static convex combination of all N = M + K = 23 experts, fitted in hindsight. Values below zero mean that the online aggregate beats this static hindsight benchmark. Endpoint labels are in 106 MW2 . Bottom: MLpol weights stack with EWLS experts aggregated into three buckets by nominal scale hγ : Fast (h < 100, red), Medium (100 ≤ h ≤ 1000, gold), and Slow (h ≥ 1000, blue); grey on top is the total weight on the M = 7 raw base experts.
6
Conclusion
We introduced MELO, a causal, model-agnostic adaptation layer for time-indexed tabular prediction under distribution shift. The method keeps raw base predictors as conservative competitors, augments them with multi-scale EWLS correction experts, and uses MLpol to select among raw and corrected forecasts by realised loss. Our deterministic oracle inequalities formalise the two intended safeguards: competitiveness with both the best raw predictor protects against unnecessary correction and the best EWLS memory scale gives a path-length-dependent guarantee for bounded affine comparator paths when residual drift is exploitable. Empirically, on the RTE-FR load stress test, MELO improves substantially across the COVID-19 lockdown, outperforming base-only aggregation and the adaptive Kalman-filter alternatives considered under matched validation tuning. The limitation is that large gains require both exploitable non-stationarity and a base pool with residual diversity. When base predictors make highly correlated errors, the EWLS combination space has limited room to help, as seen in the TabReD experiments, motivating a simple residual-correlation diagnostic on the validation period. Future work includes broader evaluations with actively retrained base models and extensions to multivariate targets, structured losses, and online expansion of the memory-scale grid.
9
References T. Akiba, S. Sano, T. Yanase, T. Ohta, and M. Koyama. Optuna: A next-generation hyperparameter optimization framework. In Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pages 2623–2631, 2019. P. Alquier. Non-exponentially weighted aggregation: Regret bounds for unbounded loss functions. In Proceedings of the 38th International Conference on Machine Learning (ICML), volume 139 of Proceedings of Machine Learning Research, pages 207–218, 2021. J. M. Bates and C. W. J. Granger. The combination of forecasts. Journal of the Operational Research Society, 20(4):451–468, 1969. doi: 10.1057/jors.1969.103. A. Bifet and R. Gavaldà. Learning from time-changing data with adaptive windowing. In Proceedings of the 2007 SIAM International Conference on Data Mining, pages 443–448. SIAM, 2007. doi: 10.1137/1. 9781611972771.42. H. A. P. Blom and Y. Bar-Shalom. The interacting multiple model algorithm for systems with Markovian switching coefficients. IEEE Transactions on Automatic Control, 33(8):780–783, 1988. doi: 10.1109/9.1299. O. Bousquet and M. K. Warmuth. Tracking a small set of experts by mixing past posteriors. Journal of Machine Learning Research, 3:363–396, 2002. N. Cesa-Bianchi and G. Lugosi. Prediction, learning, and games. Cambridge university press, 2006. N. Cesa-Bianchi, P. Gaillard, G. Lugosi, and G. Stoltz. Mirror descent meets fixed share (and feels no regret). In Advances in Neural Information Processing Systems 25 (NeurIPS), 2012. T. Chen and C. Guestrin. Xgboost: A scalable tree boosting system. In KDD, 2016. D. Dai, P. Rigollet, and T. Zhang. Deviation optimal learning using greedy Q-aggregation. The Annals of Statistics, 40(3):1878–1905, 2012. D. Dai, P. Rigollet, L. Xia, and T. Zhang. Aggregation of affine estimators. Electronic Journal of Statistics, 8(1): 302–327, 2014. J. De Vilmarest and Y. Goude. State-space models for online post-covid electricity load forecasting competition. IEEE Open Access Journal of Power and Energy, 9:192–201, 2022. J. de Vilmarest and O. Wintenberger. Viking: Variational Bayesian variance tracking. Statistical Inference for Stochastic Processes, 27(3):839–860, 2024. doi: 10.1007/s11203-024-09312-7. M. Devaine, P. Gaillard, Y. Goude, and G. Stoltz. Forecasting electricity consumption by aggregating specialized experts. Machine Learning, 90(2):231–260, 2013. doi: 10.1007/s10994-012-5314-7. P. Gaillard, G. Stoltz, and T. van Erven. A second-order bound with excess losses. COLT, 2014. P. Gaillard, Y. Goude, and R. Nedellec. Additive models and robust aggregation for GEFCom2014 probabilistic electric load and electricity price forecasting. International Journal of Forecasting, 32(3):1038–1050, 2016. doi: 10.1016/j.ijforecast.2015.12.001. J. Gama, I. Žliobaitė, A. Bifet, M. Pechenizkiy, and A. Bouchachia. A survey on concept drift adaptation. ACM Computing Surveys, 46(4):1–37, 2014. doi: 10.1145/2523813. H. M. Gomes, A. Bifet, J. Read, J. P. Barddal, F. Enembreck, B. Pfharinger, G. Holmes, and T. Abdessalem. Adaptive random forests for evolving data stream classification. In Machine Learning and Knowledge Discovery in Databases, pages 711–725. Springer, 2017. doi: 10.1007/978-3-319-71273-4_45. Y. Gorishniy, I. Rubachev, V. Khrulkov, and A. Babenko. Revisiting deep learning models for tabular data. In Advances in Neural Information Processing Systems, volume 34, 2021. Y. Gorishniy, I. Rubachev, and A. Babenko. On embeddings for numerical features in tabular deep learning. In Advances in Neural Information Processing Systems, 2022. Y. Goude, R. Nedellec, and N. Kong. Local short and middle term electricity load forecasting with semiparametric additive models. IEEE Transactions on Smart Grid, 5(1):440–446, 2014. doi: 10.1109/TSG.2013. 2278425.
10
T. Hale, N. Angrist, R. Goldszmidt, B. Kira, A. Petherick, T. Phillips, S. Webster, E. Cameron-Blake, L. Hallas, S. Majumdar, et al. A global panel database of pandemic policies (oxford covid-19 government response tracker). Nature human behaviour, 5(4):529–538, 2021. T. Hastie and R. Tibshirani. Generalized additive models. Statistical Science, 1(3):297–318, 1986. doi: 10.1214/ss/1177013604. M. H. Hayes. Statistical digital signal processing and modeling. John Wiley & Sons, 1996. K. He, X. Zhang, S. Ren, and J. Sun. Deep residual learning for image recognition. In IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 770–778, 2016. doi: 10.1109/CVPR.2016.90. M. Herbster and M. K. Warmuth. Tracking the best expert. Machine learning, 32(2):151–178, 1998. A. E. Hoerl and R. W. Kennard. Ridge regression: Biased estimation for nonorthogonal problems. Technometrics, 12(1):55–67, 1970. doi: 10.1080/00401706.1970.10488634. N. Hollmann et al. Tabpfn: A transformer that solves small tabular classification problems in a second. In ICLR, 2023. Y. Huang, Y. Zhang, Z. Wu, N. Li, and J. Chambers. A novel adaptive Kalman filter with inaccurate process and measurement noise covariance matrices. IEEE Transactions on Automatic Control, 63(2):594–601, 2018. doi: 10.1109/TAC.2017.2730480. R. Kalman. A new approach to linear filtering and prediction problems. Journal of Basic Engineering, 1960. G. Ke et al. Lightgbm: A highly efficient gradient boosting decision tree. In NeurIPS, 2017. H. Kim, S. Kim, J. Mok, and S. Yoon. Battling the non-stationarity in time series forecasting via test-time adaptation. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 39, pages 17868–17876, 2025. doi: 10.1609/aaai.v39i17.33965. H. R. Künsch. The jackknife and the bootstrap for general stationary observations. The Annals of Statistics, 17 (3):1217–1241, 1989. T. L. Lee, W. Toner, R. Singh, A. Joosen, and M. Asenov. Lightweight online adaption for time series foundation model forecasts. In A. Singh, M. Fazel, D. Hsu, S. Lacoste-Julien, F. Berkenkamp, T. Maharaj, K. Wagstaff, and J. Zhu, editors, Proceedings of the 42nd International Conference on Machine Learning, volume 267 of Proceedings of Machine Learning Research, pages 33736–33764. PMLR, 2025. Z. Liu, M. Cheng, G. Zhao, J. Yang, Q. Liu, and E. Chen. Improving time series forecasting via instance-aware post-hoc revision. In Advances in Neural Information Processing Systems, volume 39, 2025. J. Miao and P. Polak. Online ensemble learning for sector rotation: A gradient-free framework. arXiv preprint arXiv:2304.09947, 2023. D. Obst, J. de Vilmarest, and Y. Goude. Adaptive methods for short-term electricity load forecasting during COVID-19 lockdown in France. IEEE Transactions on Power Systems, 36(5):4754–4763, 2021. doi: 10.1109/TPWRS.2021.3067551. J. Qu, D. Holzmüller, G. Varoquaux, and M. L. Morvan. TabICL: A tabular foundation model for in-context learning on large data. In Forty-second International Conference on Machine Learning, 2025. I. Rubachev, N. Kartashev, Y. Gorishniy, and A. Babenko. Tabred: Analyzing pitfalls and filling the gaps in tabular deep learning benchmarks. In International Conference on Learning Representations, 2025. A. H. Sayed. Fundamentals of adaptive filtering. John Wiley & Sons, 2003. J. H. Stock and M. W. Watson. Combination forecasts of output growth in a seven-country data set. Journal of Forecasting, 23(6):405–430, 2004. doi: 10.1002/for.928. A. Timmermann. Forecast combinations. In G. Elliott, C. W. J. Granger, and A. Timmermann, editors, Handbook of Economic Forecasting, Volume 1, chapter 4, pages 135–196. Elsevier, 2006. O. Wintenberger. Optimal learning with bernstein online aggregation. Machine Learning, 106(1):119–141, 2017. S. N. Wood. Generalized additive models: an introduction with R. chapman and hall/CRC, 2017. L. Zhang, S. Lu, and Z.-H. Zhou. Adaptive online learning in dynamic environments. In Advances in Neural Information Processing Systems 31 (NeurIPS), pages 1323–1333, 2018. M. Zinkevich. Online convex programming and generalized infinitesimal gradient ascent. In Proceedings of the 20th International Conference on Machine Learning (ICML), pages 928–936, 2003.
11
A
Proofs for the pathwise online guarantees
This appendix proves the deterministic guarantees stated in Section 4: the dynamic-regret bound for a single forgotten-RLS correction expert, the MLpol oracle inequality over the enlarged Base+EWLS pool, and the geometric-grid corollary. The analysis is sequence-wise: apart from the explicitly stated bounded-iterate (or, equivalently, uniform residual/separation) condition, we impose no distributional, stationarity, or persistence-of-excitation assumptions on the realised sequence. The proof has five ingredients. First, a one-step potential identity for forgotten RLS under squared loss. Second, a log-determinant inequality adapted to the non-monotone recursion At = γAt−1 + z̃t z̃t⊤ . Third, a dynamic-comparator telescoping argument that converts the fixed-comparator identity into a path-length bound. Fourth, an MLpol potential bound for the linearised excess losses. Fifth, a local-band refinement showing that, for dense EWLS grids, the effective aggregation cost can depend on the inverse mass of a useful time-scale band rather than on the raw number of grid points. A.1
Setup and notation
We use the same notation as in the main text. At round t, the EWLS correction expert observes the augmented base-prediction vector z̃t = (zt⊤ , 1)⊤ ∈ Rd , where d = M + 1, predicts z̃t⊤ wt , then observes yt and incurs squared loss 2 ℓt (w) = z̃t⊤ w − yt . Here wt denotes the pre-update coefficient used to predict at round t. The proof below uses only boundedness of z̃t , not any property of the raw covariates or base models. A.2
Map to the informal statements in Section 4
The main text states Theorems 1 and 2 in informal form, suppressing constants and lower-order terms for readability. The corresponding formal statements are Theorem A.3 (single-scale tracking) and Theorem A.6 (MELO oracle inequality), proved in Appendices A.5 and A.8. The geometric-grid statement of Section 4 is formalised as Corollary A.7. We restate the correspondence here so that the notation in the formal statements is unambiguous. The function Ψh . The single-scale tracking penalty appearing implicitly in the informal theorems is, with all constants made explicit, 2D2 dT B2 h Ψh (T, P ) = |{z} δR2 + D2 d log 1 + δz + + 4R(δ + Bz2 h) P , (4) {z } | h | {z } | {z } initialisation log term
finite-memory term
tracking term
where D = By + Bz R and d = M + 1. The "log and initialisation terms" referenced in the informal Theorem 1 are the first two summands; the constants C1 , C2 in the informal statement correspond to C1 = 2D2 ,
C2 = 8RBz2
(under the mild regime δ ≤ Bz2 h).
The unsimplified C2 is 4R(δ + Bz2 h)/h, which reduces to the above once δ ≤ Bz2 h; see (15). Equivalent affine view of the aggregate. The informal Theorem 2 compares ŷt to time-varying affine combinations of the base-prediction vector. The justification is that every candidate in the MELO pool is itself affine in z̃t : each raw base expert is a coordinate of zt , and each EWLS expert (γ) is z̃t⊤ wt . Hence the MLpol aggregate satisfies ŷt = z̃t⊤ w̄t ,
w̄t =
N X
pt,j wt,j ,
j=1
for coefficients w̄t determined causally by past losses and EWLS states, where wt,j is the affine (γ ) representation of expert j (a canonical basis vector for raw base experts, wt k for EWLS experts). The comparator class UT (R) in Theorem A.3 is therefore the natural benchmark: it is the same affine class that ŷt itself ranges over. 12
Clipping convention for the named base forecasts. In the formal oracle inequality (24), LT (fj ) denotes the loss of the clipped raw prediction sequence, since clipping is applied to all experts before they enter MLpol. When the unclipped raw predictions are the deployed base forecasts and outcomes lie in [−B, B], Lemma A.4 shows that clipping cannot increase squared loss against the target. Hence the bound stated against the clipped LT (fj ) also holds against the unclipped raw-forecast loss in any practical reading where |yt | ≤ B. This is the convention used throughout Section 5. Forgotten-RLS recursion. Fix a forgetting factor γ ∈ (0, 1) and define the nominal scale h = (1 − γ)−1 . Let A0 = δI with δ > 0, b0 = 0, and w1 = 0. The forgotten-RLS recursion is At = γAt−1 + z̃t z̃t⊤ ,
bt = γbt−1 + yt z̃t ,
wt+1 = A−1 t bt .
(5)
Equivalently, At = γ t δI +
t X
γ t−s z̃s z̃s⊤ .
s=1
Since At ≻ 0, the update is well-defined. Moreover, because At−1 wt = bt−1 , a direct calculation gives 1 wt+1 = wt − A−1 gt , gt := ∇ℓt (wt ) = 2(z̃t⊤ wt − yt )z̃t . (6) 2 t Boundedness assumptions.
We use the following deterministic boundedness conditions:
(B1) ∥z̃t ∥2 ≤ Bz for all t. (B2) |yt | ≤ By for all t. (B3) The comparator path satisfies ∥ut ∥2 ≤ R for all t. (B4) The unprojected forgotten-RLS iterates are stable along the realised sequence: ∥wt ∥2 ≤ R for all t. Condition (B4) is an explicit stability assumption on the unprojected RLS recursion, not a consequence of forgetting alone. Without a conditioning or persistence-of-excitation condition, At may be illconditioned in directions weakly observed by the realised sequence, and wt = A−1 t−1 bt−1 need not be uniformly bounded. We state the pathwise tracking bound conditionally on this stability event, thereby avoiding any stochastic excitation assumption. Empirical diagnostic for (B4). Appendix B.7 provides a sanity check for this condition on the RTE-FR Load experiment. Since the augmented EWLS coefficient contains an intercept on the response scale, the diagnostic separates the full augmented norm from the slope-only norm. Along the realised test trajectory, the slope component remains uniformly bounded across the K − 1 = 15 finite-forgetting EWLS experts and all T = 746 test days, with (γ)
sup ∥βt ∥2 = 1.69. t,γ
Alternatively, the proof can be read under the direct bounds |z̃t⊤ wt − yt | ≤ D
and
∥wt − ut ∥2 ≤ 2R
for the comparator class under consideration. We do not claim that Euclidean projection preserves the exact recursion (5); the proof below is for the unprojected recursion. In the implementation, predictions are clipped before being passed to MLpol, which ensures bounded losses for the aggregation layer; see Lemma A.4 and Remark A.4. Set D := By + Bz R. Under (B1), (B2), and (B4), the realised RLS residuals satisfy |z̃t⊤ wt − yt | ≤ D. Under (B1), (B2), and (B3), any comparator ut also satisfies |z̃t⊤ ut − yt | ≤ D. 13
For a comparator path u1:T = (u1 , . . . , uT ), define LT (u1:T ) =
T X (z̃t⊤ ut − yt )2 ,
PT (u1:T ) =
t=1
t=2
We write (γ)
y̌t
= z̃t⊤ wt ,
T X
LT (y̌ (γ) ) =
T X t=1
(γ)
y̌t
∥ut − ut−1 ∥2 .
− yt
(7)
2
for the prediction sequence and cumulative loss of the forgotten-RLS expert with nominal scale h = (1 − γ)−1 . A.3
A one-step potential identity
Lemma A.1 (One-step identity for forgotten RLS). For every u ∈ Rd and every t ≥ 1, 1 ℓt (wt ) − ℓt (u) = γ∥wt − u∥2At−1 − ∥wt+1 − u∥2At + gt⊤ A−1 t gt . 4
(8)
Proof. For squared loss, expanding around wt gives the exact identity 2 ℓt (wt ) − ℓt (u) = gt⊤ (wt − u) − z̃t⊤ (wt − u) .
(9)
From (6), gt = 2At (wt − wt+1 ).
Therefore
gt⊤ (wt − u) = 2(wt − wt+1 )⊤ At (wt − u)
= ∥wt − u∥2At − ∥wt+1 − u∥2At + ∥wt − wt+1 ∥2At 1 = ∥wt − u∥2At − ∥wt+1 − u∥2At + gt⊤ A−1 t gt . 4
(10)
The recursion At = γAt−1 + z̃t z̃t⊤ implies 2 ∥wt − u∥2At = γ∥wt − u∥2At−1 + z̃t⊤ (wt − u) . 2 Substituting this expression and (10) into (9) cancels the curvature term z̃t⊤ (wt − u) , yielding (8). A.4
A forgetting-adapted log-determinant inequality
The classical online-Newton log-determinant argument relies on monotonicity of At = At−1 + z̃t z̃t⊤ . With forgetting, At = γAt−1 + z̃t z̃t⊤ is no longer monotone in the PSD order. The next lemma is the corresponding log-determinant bound for forgotten RLS. Lemma A.2 (Forgetting-adapted log-determinant inequality). Let At = γAt−1 + z̃t z̃t⊤ with A0 = δI, δ > 0, γ = 1 − 1/h, and h ≥ 2. Under (B1), T X Bz2 h 2dT z̃t⊤ A−1 z̃ ≤ d log 1 + + . (11) t t δ h t=1 Proof. Define qt := z̃t⊤ (γAt−1 )−1 z̃t = γ −1 z̃t⊤ A−1 t−1 z̃t ≥ 0. By the Sherman–Morrison formula, −1 A−1 − t = (γAt−1 )
(γAt−1 )−1 z̃t z̃t⊤ (γAt−1 )−1 . 1 + z̃t⊤ (γAt−1 )−1 z̃t
Contracting with z̃t gives z̃t⊤ A−1 t z̃t = 14
qt . 1 + qt
(12)
The matrix-determinant lemma gives det At = det(γAt−1 ) 1 + qt = γ d det(At−1 )(1 + qt ). Thus log(1 + qt ) = log det At − log det At−1 + d log(1/γ).
Summing over t and using log(1 + q) ≥ q/(1 + q) for q ≥ 0, T X
z̃t⊤ A−1 t z̃t =
t=1
T X qt ≤ log(1 + qt ) 1 + qt t=1 t=1
T X
= log
det AT + dT log(1/γ). det A0
(13)
Next, AT = γ T δI +
T X s=1
γ T −s z̃s z̃s⊤ ⪯ δ + Bz2
Therefore log
T −1 X j=0
γ j I ⪯ (δ + Bz2 h)I.
det AT B2h . ≤ d log 1 + z det A0 δ
Finally, since h ≥ 2, γ = 1 − 1/h ≥ 1/2, and log(1/γ) = − log(1 − 1/h) ≤
2 . h
Substituting into (13) proves (11). Remark A.1 (Origin of the T /h term). When γ = 1, the second term in (13) vanishes and one recovers the classical log-determinant control. Forgetting contributes dT log(1/γ) ≈ dT /h. This is the explicit cost of using a finite memory scale rather than retaining all past curvature. A.5
Dynamic regret of a single forgotten-RLS expert
Theorem A.3 (Dynamic regret of forgotten RLS). Assume (B1)–(B4), let h = (1 − γ)−1 ≥ 2, and let u1:T be any comparator path satisfying (B3). Then Bz2 h 2D2 dT (h) 2 2 LT (y̌ ) − LT (u1:T ) ≤ δR + D d log 1 + + + 4R(δ + Bz2 h)PT (u1:T ). (14) δ h In particular, if δ ≤ Bz2 h, then LT (y̌
(h)
B2h ) − LT (u1:T ) ≤ δR + D d log 1 + z δ 2
2
+
2D2 dT + 8RBz2 h PT (u1:T ). h
(15)
Ignoring logarithmic and initialization terms, the h-dependent part is dT O + hPT (u1:T ) . h Proof. Apply Lemma A.1 with u = ut and sum over t: LT (y̌ (h) ) − LT (u1:T ) =
T X t=1
ℓt (wt ) − ℓt (ut )
T h T i 1X X ≤ γ∥wt − ut ∥2At−1 − ∥wt+1 − ut ∥2At + gt⊤ A−1 t gt . 4 t=1 t=1 {z } | {z } | S1
15
S2
(16)
Bounding S2 . Hence
By (B1), (B2), and (B4), |z̃t⊤ wt − yt | ≤ D,
2 ⊤ −1 gt⊤ A−1 t gt ≤ 4D z̃t At z̃t .
Using Lemma A.2, S2 ≤ D2 Bounding S1 .
gt = 2(z̃t⊤ wt − yt )z̃t .
T X
2D2 dT Bz2 h 2 + . z̃t⊤ A−1 z̃ ≤ D d log 1 + t t δ h t=1
Shift the second sum in S1 : S1 = γ∥w1 − u1 ∥2A0 − ∥wT +1 − uT ∥2AT T h i X + γ∥wt − ut ∥2At−1 − ∥wt − ut−1 ∥2At−1 .
(17)
(18)
t=2
The terminal term is non-positive and can be dropped. Since w1 = 0, A0 = δI, γ ≤ 1, and ∥u1 ∥ ≤ R, γ∥w1 − u1 ∥2A0 ≤ δR2 . For the bracketed term, let xt := wt − ut , ∆t := ut − ut−1 . Then wt − ut−1 = xt + ∆t , and 2 γ∥xt ∥2At−1 − ∥xt + ∆t ∥2At−1 = −(1 − γ)∥xt ∥2At−1 − 2∆⊤ t At−1 xt − ∥∆t ∥At−1 ≤ 2|∆⊤ t At−1 xt |.
Moreover, At−1 = γ
t−1
δI +
t−1 X s=1
(19)
γ t−1−s z̃s z̃s⊤ ⪯ (δ + Bz2 h)I.
By (B3)–(B4), ∥xt ∥ ≤ 2R. Therefore 2 2|∆⊤ t At−1 xt | ≤ 2∥∆t ∥ ∥At−1 ∥op ∥xt ∥ ≤ 4R(δ + Bz h)∥ut − ut−1 ∥. Summing (19) over t = 2, . . . , T , S1 ≤ δR2 + 4R(δ + Bz2 h)PT (u1:T ). Combining this with (17) proves (14). If δ ≤ Bz2 h, then δ + Bz2 h ≤ 2Bz2 h, giving (15). Remark A.2 (About bounded iterates and projection). The exact identity (6) is what makes Lemma A.1 possible. A Euclidean projection after the RLS update would generally destroy this identity. One can instead analyse an At -metric projection, as in online Newton methods, but the implementation in this paper is closer to the unprojected RLS recursion with prediction clipping. We therefore state the deterministic bound under the bounded-iterate condition (B4). The same proof also goes through under direct bounds |z̃t⊤ wt − yt | ≤ D and ∥wt − ut ∥ ≤ 2R for the comparator class under consideration. A.6
Clipping and bounded expert losses
The MLpol bound below requires bounded outcomes and bounded expert predictions. We use the standard clipping device. Lemma A.4 (Clipping does not increase squared loss). Let |y| ≤ B, and define clipB (a) := max{−B, min{a, B}}. Then 2 clipB (a) − y ≤ (a − y)2 for all a ∈ R. Proof. If a ∈ [−B, B], the claim is equality. If a > B, then clipB (a) = B ≥ y, so |B − y| ≤ |a − y|. The case a < −B is symmetric. Thus the clipped forgotten-RLS expert has cumulative squared loss no larger than the un-clipped expert analysed in Theorem A.3, provided the target lies in [−B, B]. 16
A.7
MLpol regret over a finite pool
We next record the MLpol potential bound used for the aggregation layer. We state a distributional form against any fixed convex combination of experts: the point-expert oracle inequality is the special case where the distribution is a point mass, and the local-band refinement below uses the uniform distribution over a band of nearby forgetting factors. Let N be the number of experts. At round t, write ỹt,j for the prediction of expert j, and let ŷt =
N X
pt,j ỹt,j
j=1
be the aggregate prediction. We use the linearised excess loss rt,j = 2(ŷt − yt )(ŷt − ỹt,j ),
Rt,j =
t X
rs,j .
s=1
The MLpol weights are [Rt−1,j ]+ pt,j = PN , i=1 [Rt−1,i ]+ with the uniform distribution used when all positive parts vanish. Lemma A.5 (MLpol finite-pool regret, distributional form). Assume that all expert predictions and outcomes lie in [−B, B]. For any fixed distribution q ∈ ∆N , define the corresponding convex-mixture prediction N X q yt = qj ỹt,j . j=1
Then there exists a universal constant C > 0 such that T X t=1
(ŷt − yt )2 ≤
T X t=1
(ytq − yt )2 + CB 2
q T N ∥q∥22 .
(20)
In particular, taking q to be a point mass gives T T X X √ (ỹt,j − yt )2 + CB 2 T N . (ŷt − yt )2 ≤ min 1≤j≤N
t=1
(21)
t=1
Moreover, if S ⊆ {1, . . . , N } is non-empty and q S is uniform on S, then s T T X X N 1 X 2 S 2 2 (ŷt − yt ) ≤ (ȳt − yt ) + CB T , ȳtS = ỹt,j . |S| |S| t=1 t=1
(22)
j∈S
p Remark A.3. A tighter, second-order bound of the form O(B (N · LT )) can be obtained via the MLpol potential analysis of [Gaillard et al., 2014]. We use the simpler first-order form here because the leading rate-determining term in Theorem 2 is the EWLS tracking term Ψh (T, PT ), not the MLpol overhead. Proof. Define the potential N
Φt =
1 X [Rt,j ]2+ . 2N j=1
For any scalar x, a, the inequality [x + a]2+ − [x]2+ ≤ 2[x]+ a + a2 holds. Therefore Φt − Φt−1 ≤
N N 1 X 1 X 2 [Rt−1,j ]+ rt,j + r . N j=1 2N j=1 t,j
17
The first term is zero. Indeed, if some positive part is non-zero, then N N N X X X [Rt−1,j ]+ rt,j = [Rt−1,j ]+ pt,j rt,j j=1
j=1
j=1
N N X X = [Rt−1,j ]+ 2(ŷt − yt ) ŷt − pt,j ỹt,j = 0. j=1
j=1
If all positive parts vanish, the same term is trivially zero. Since ŷt , yt , ỹt,j ∈ [−B, B], |rt,j | = 2|ŷt − yt | |ŷt − ỹt,j | ≤ 8B 2 . Thus Φt − Φt−1 ≤ CB 4 , and hence ΦT ≤ CB 4 T.
Now fix q ∈ ∆N . For each round t, by the elementary identity a2 − b2 = 2a(a − b) − (a − b)2 , with a = ŷt − yt and b = ytq − yt , we have N X (ŷt − yt )2 − (ytq − yt )2 ≤ 2(ŷt − yt )(ŷt − ytq ) = qj rt,j . j=1
Summing over t, LT (ŷ) − LT (y q ) ≤
N X j=1
qj RT,j ≤
N X
qj [RT,j ]+ .
j=1
By Cauchy–Schwarz, 1/2 N X 1 2 qj [RT,j ]+ ≤ N [RT,j ]+ N j=1 j=1 j=1 q q p = N ∥q∥22 2ΦT ≤ CB 2 T N ∥q∥22 .
N X
N X
1/2
qj2
This proves (20). The point-mass and uniform-on-S special cases follow from ∥ej ∥22 = 1 and ∥q S ∥22 = 1/|S|, respectively. A.8
Oracle inequality for the Base+EWLS pool
Let H = {h1 , . . . , hK } be a finite set of nominal EWLS scales. For each h ∈ H, let y̌ (h) denote the corresponding forgotten-RLS correction sequence, with predictions clipped to [−B, B] before aggregation. The enlarged pool is IH = {f1 , . . . , fM } ∪ {y̌ (h) : h ∈ H}. Let N = M + K. Define
2D2 dT B2h + + 4R(δ + Bz2 h)P. (23) Ψh (T, P ) := δR2 + D2 d log 1 + z δ h Theorem A.6 (Oracle inequality for multi-scale online correction). Assume the conditions of Theorem A.3 hold for each h ∈ H, and assume outcomes lie in [−B, B]. All predictions entering MLpol are clipped to [−B, B]. Then, writing fj and y̌ (h) for the clipped raw and EWLS prediction sequences, for every comparator path u1:T satisfying (B3), p LT (ŷ) ≤ min min LT (fj ), LT (u1:T ) + min Ψh (T, PT (u1:T )) +CB 2 T (M + K). (24) h∈H
1≤j≤M
Proof. By the point-mass special case of Lemma A.5 applied to the enlarged pool, p LT (ŷ) ≤ min min LT (fj ), min LT (y̌ (h) ) + CB 2 T (M + K). 1≤j≤M
h∈H
By Lemma A.4, clipping cannot increase squared loss. Therefore Theorem A.3 gives, for every h ∈ H, LT (y̌ (h) ) ≤ LT (u1:T ) + Ψh (T, PT (u1:T )). Taking the minimum over h ∈ H proves (24). 18
A.9
Geometric grids and unknown correction scale
The leading h-dependent terms in Ψh are C1 dT + C2 hP, h where P = PT (u1:T ). If P = 0, the leading expression is minimised by the largest available memory scale, corresponding to the no-drift case. For P > 0, the continuous minimizer is r C1 dT ⋆ h = , C2 P √ with minimum value 2 C1 C2 dT P . Corollary A.7 (Geometric-grid adaptation). Let P = PT (u1:T ) > 0, and let H be a geometric grid over [hmin , hmax ] with ratio ρ > 1. Suppose h⋆ ∈ [hmin , hmax ]. Then there exists hk ∈ H such that √ C1 dT + C2 hk P ≤ Cρ dT P , hk
(25)
where Cρ depends only on ρ, C1 , C2 . Consequently, ignoring logarithmic and initialization terms, p p LT (ŷ) ≤ min min LT (fj ), LT (u1:T ) + Cρ dT PT (u1:T ) + O B 2 T (M + K) . (26) j
Proof. Since H is geometric and contains h⋆ within its range, there is a grid point hk with h⋆ ≤ hk ≤ ρh⋆
h⋆ /ρ ≤ hk ≤ h⋆ .
or
In either case, h⋆ ≤ ρ. hk
hk ≤ ρ, h⋆ Therefore C1 dT + C2 hk P ≤ ρ hk
C1 dT + C 2 h⋆ P h⋆ p = 2ρ C1 C2 dT P .
√ Absorbing 2ρ C1 C2 into Cρ proves (25). Substituting this bound into Theorem A.6 gives (26). A.10
Local-band refinement for dense EWLS grids
The oracle inequality above compares MLpol with the pbest single EWLS scale in a finite grid, which gives a worst-case aggregation overhead of order T (M + K). This point-oracle view can be pessimistic for dense EWLS grids. Nearby forgetting factors often generate very similar prediction sequences, and the useful comparator may be a local band of time scales rather than a single isolated expert. The distributional form of Lemma A.5 makes this precise. Let the enlarged pool contain M base experts and K EWLS experts, so that N = M + K. For a non-empty subset S ⊆ {1, . . . , K} of EWLS scales, define the band-averaged EWLS prediction ȳtS =
1 X (hk ) y̌t . |S| k∈S
Corollary A.8 (Band-oracle bound for EWLS grids). Assume all base predictions, EWLS predictions after clipping, and outcomes lie in [−B, B]. Then for every non-empty EWLS band S ⊆ {1, . . . , K}, s M +K LT (ŷ) ≤ LT (ȳ S ) + CB 2 T . (27) |S| 19
Proof. Apply Lemma A.5 with N = M + K to the distribution q S that assigns mass 1/|S| to the EWLS experts in S and zero mass to all other experts, including the M base experts. Then ∥q S ∥22 =
1 , |S|
and the convex-mixture prediction induced by q S is precisely ȳtS . Substituting this into (20) gives (27). We next relate the band-average comparator to a single central forgetting scale. The following continuity statement is for the idealised forgotten-RLS recursion analysed above, without covariance inflation. It is used only to justify why nearby forgetting factors in a dense grid should not be viewed as unrelated experts. It is convenient to parametrize forgetting factors by h = (1 − γ)−1 .
θ = log h,
Write y̌ (θ) for the clipped EWLS prediction sequence with h = eθ . Lemma A.9 (Continuity of finite-horizon EWLS predictions). Fix T < ∞, δ > 0, and a compact (θ) (θ) interval Θ ⊂ R. For every t ≤ T , the idealised EWLS coefficient wt and prediction y̌t are −θ continuous functions of θ ∈ Θ. If Θ is compact and γ(θ) = 1 − e stays in (0, 1), then the map is Lipschitz: there exists Lθ < ∞, depending on the realised finite sequence and on Θ, such that (θ)
max |y̌t
1≤t≤T
(θ ′ )
− y̌t
| ≤ Lθ |θ − θ′ |
for all θ, θ′ ∈ Θ.
(28)
Proof. For the idealised EWLS recursion, At (γ) = γ t δI +
t X
γ t−s z̃s z̃s⊤ ,
bt (γ) =
s=1
t X
γ t−s ys z̃s .
s=1
Both At (γ) and bt (γ) are polynomial functions of γ. Moreover At (γ) ⪰ γ t δI ≻ 0 for γ ∈ (0, 1], so At (γ)−1 is continuous wherever γ > 0. Hence (γ)
wt+1 = At (γ)−1 bt (γ) is continuous in γ, and therefore in θ through the continuous map γ(θ) = 1 − e−θ . On a compact interval contained in (0, 1), the derivative is bounded, so the map is Lipschitz. The prediction (θ) (θ) y̌t = z̃t⊤ wt , followed by clipping, preserves continuity and Lipschitzness. Corollary A.10 (Local-band control for dense EWLS grids). Assume the Lipschitz condition (28). Fix a target scale θ⋆ ∈ Θ and a radius r > 0. Let Sr (θ⋆ ) = {k ∈ {1, . . . , K} : |θk − θ⋆ | ≤ r} ,
θk = log hk ,
and assume Sr (θ⋆ ) ̸= ∅. Then
s ⋆
LT (ŷ) ≤ LT (y̌ (θ ) ) + 4BLθ T r + CB 2
T
M +K . |Sr (θ⋆ )|
(29)
⋆
Consequently, combining with Theorem A.3, for h⋆ = eθ , s LT (ŷ) ≤ LT (u1:T ) + Ψh⋆ (T, PT (u1:T )) + 4BLθ T r + CB 2 Proof. For every k ∈ Sr (θ⋆ ), the Lipschitz condition gives (θ ⋆ )
(h )
max |y̌t k − y̌t
1≤t≤T
| ≤ Lθ r.
Therefore the band average satisfies S (θ ⋆ )
max |ȳt r
1≤t≤T
(θ ⋆ )
− y̌t
20
| ≤ Lθ r.
T
M +K . |Sr (θ⋆ )|
(30)
For clipped predictions and outcomes in [−B, B], squared loss is 4B-Lipschitz in the prediction: |(a − y)2 − (b − y)2 | ≤ 4B|a − b|, Hence
⋆
a, b, y ∈ [−B, B].
⋆
LT (ȳ Sr (θ ) ) ≤ LT (y̌ (θ ) ) + 4BLθ T r. Combining this inequality with Corollary A.8 proves (29). The final display follows by applying ⋆ Theorem A.3 to the EWLS scale h⋆ = eθ . A.11
Implementation remarks: covariance inflation and nominal scale
Remark A.4 (Covariance inflation). The deterministic proof above analyses the idealised forgottenRLS recursion (5). The implementation used in the experiments adds a small covariance-inflation term to the Pt = A−1 recursion for numerical stability and responsiveness, with the inflation t scale selected on the 2018 walk-forward validation period. This modification is an implementation device rather than the source of the theoretical stability condition: adding εI to the covariance caps overconfident precision directions but does not by itself impose a uniform lower bound on λmin (At ). We therefore treat h(γ) = 1/(1 − γ) in the experiments as a nominal forgetting scale rather than a calibrated number of days of memory.
21
Algorithm 1 Multi-Time-Scale EWLS with MLpol Aggregation M
Require: Base experts {f (j) }j=1 ; forgetting factors Γ = {γ1 , . . . , γK }; inflation levels {εk }K k=1 ; diffuse-prior scale δ0 > 0; optional clipping radius B 1: Set N ← M + K (k) (k) −1 2: Initialise: w1 = 0 and P0 = δ0 IM +1 for k = 1, . . . , K; p1 = 1/N ; R0 = 0 3: for t = 1, 2, . . . , T do ⊤ 4: Observe xt and form zt = f (1) (xt ), . . . , f (M ) (xt ) ; set z̃t = (zt⊤ , 1)⊤ 5: // Expert predictions 6: ỹt,j = zt,j for j = 1, . . . , M {raw base experts} (k) 7: ỹt,M +k = z̃t⊤ wt for k = 1, . . . , K {EWLS correction experts} 8: Optionally clip each ỹt,j to [−B, B] 9: // MLpol aggregate prediction PN 10: ŷt = j=1 pt,j ỹt,j 11: Observe yt and set gt = 2(ŷt − yt ) 12: // EWLS/RLS updates 13: for k = 1, . . . , K do (k) (k) 14: st ← γk + z̃t⊤ Pt−1 z̃t (k) (k) (k) 15: Kt ← Pt−1 z̃t /st {RLS/Kalman gain} (k) (k) (k) (k) 16: wt+1 ← wt + Kt yt − z̃t⊤ wt 17:
(k)
(k)
(k)
(k)
(k)
Pt|t ← γk−1 Pt−1 − Pt−1 z̃t z̃t⊤ Pt−1 /st (k)
(k)
{small inflation carried to the next round} Pt ← Pt|t + εk I 19: end for 20: // MLpol weight update 21: r̃t,j ← gt (ŷt − ỹt,j ) for all j = 1, . . . , N 22: Rt ← Rt−1 + r̃t PN P 23: pt+1,j ← [Rt,j ]+ / i=1 [Rt,i ]+ for all j = 1, . . . , N , using pt+1,j = 1/N if i [Rt,i ]+ = 0 24: end for 18:
B
Experimental protocol and implementation details
B.1
The algorithm
Algorithm 1 gives the recursive update after the cold-start initialization used in the experiments. The EWLS state update always uses the unprojected RLS recursion. The clipping line in the algorithm is optional and is included only to align with the bounded-loss version of the MLpol oracle inequality in Appendix A; it is not activated in the reported experiments. Thus the empirical results use the unclipped expert predictions, while the theoretical aggregation guarantee applies to the standard clipped variant. The single-scale EWLS tracking bound is stated separately for the idealised unprojected recursion under the explicit stability condition. B.2
RTE-FR Load data and splits
Raw series and daily target. We use the public French national electricity-load data described in Section 5. The native short-term load-forecasting problem is intraday: load is observed at half-hourly resolution, and production systems typically forecast a vector of future half-hourly loads. In this paper we focus on a scalar daily-load benchmark. For each calendar day t, the target is the average load over the intraday half-hourly observations, nt 1 X yt = Loadt,q , nt q=1 where Loadt,q denotes the q-th half-hourly load observation on day t. On ordinary days nt = 48. Thus the target remains an average power quantity, and RMSE is reported in MW rather than in daily energy units. This aggregation removes the need to model the intraday load shape and focuses the evaluation on non-stationarity in the daily demand level. 22
Weather covariates. All methods are given the same temperature covariates constructed from realised temperature observations rather than from day-ahead weather forecasts. This is a controlledcovariate design: it prevents differences in weather-forecast quality from confounding the comparison between online adaptation methods. The experiment should therefore be read as an ex-post benchmark of load forecasting conditional on common weather information, not as a claim that realised target-day temperature would be available in a real deployment. In deployment, the same pipeline would instead receive forecast-temperature covariates. Feature availability and leakage control. The feature vector contains calendar variables, holiday indicators, temperature variables, and lagged load features. Calendar and holiday variables are known before the target day. Lagged load features are constructed only from strictly previous days, so no future load values enter the feature vector. During the test period, all online updates are sequential: each method first predicts yt , then observes yt , and only then updates its EWLS/RLS states and MLpol weights. Thus the online adaptation step is causal with respect to load labels. The only ex-post covariate convention is the realised-temperature choice above, which is shared by all compared methods. Train, validation, and test periods. The data span 2012-01-01 to 2021-01-15. Hyperparameters are selected only within the pre-test period: 2018 is used as a walk-forward validation window when needed, and selected static base models are then refit on the full 2012–2018 window. All reported online results are evaluated sequentially on 2019-01-01 to 2021-01-15. The test period is split into pre-lockdown, lockdown, and post-lockdown recovery regimes as described in Section 5. No observation from the 2019–2021 test period is used for base-model tuning, EWLS-grid selection, or nuisance-parameter selection. B.3
Base-model features and tuning
This appendix documents how base-model hyperparameters are selected for the main experiment. The guiding rule is that no base-model hyperparameter is selected using any part of the test period beginning on 2019-01-01. Scope.
Of the M = 7 base predictors, hyperparameters are set as follows:
• Lag-1: no hyperparameters. • GAM: structural specification taken from Gaillard et al. [2016]; fitted with library-level defaults of mgcv, with no test-period tuning. The exact formula is given below. • Ridge, XGBoost, LightGBM, ResNet, FT-Transformer: tuned on a held-out 2018 validation window within the training period, using the Optuna-based protocol described below. For reproducibility, the released code includes both the selected configurations used in the main experiment and the tuning routine used to obtain them. Re-running the tuning stage is therefore optional for reproducing the reported results. GAM specification. The additive structure follows the semi-parametric load-forecasting template established by Goude et al. [2014] and refined in Gaillard et al. [2016], which couples thin-plate smooths of exponentially smoothed temperatures with calendar effects and lagged load. The exponentially smoothed temperature features (Temps95 and its daily min/max) and the toy time-of-year covariate with cyclic spline basis follow the conventions of the opera R package [Gaillard et al., 2016]. Writing ηt for the linear predictor of daily load at day t, ηt = te(ℓ(Loadt−1 ), f(WeekDayt )) + ℓ(Loadt−7 ) s95,min + s(Tempt ) + s(Temps95 ) + s(Temps95,max ) t t ) + s(Tempt + f(WeekDayt ) + f(BHt ) + scp (toyt ) + s(Montht ),
where ℓ(·) is a linear term, s(·) a thin-plate regression smooth, scp (·) a cyclic cubic smooth (used on toy, time-of-year), f(·) a categorical factor, and te(·, ·) a tensor-product interaction. The covariates Temps95 and its daily min/max variants are exponentially smoothed temperatures (smoothing factor 0.95); BH flags French bank holidays. All smoothing parameters are estimated by REML; no covariate, basis dimension, or smoothness penalty is tuned on the test period. 23
Train/validation split. Each tuned model is fit on 2012-01-01 to 2017-12-31 and scored on a 2018 validation window (2018-01-01 to 2018-12-31) by RMSE on the one-step-ahead prediction task. After tuning, the selected configuration is refit on the full 2012–2018 window for evaluation on the 2019–2021 test period. No test-period information enters either the search or the final fit. Separation of convergence and modelling hyperparameters. Convergence controls—optimiser max_iter/tol for linear models, training epochs/patience for deep models—are fixed a priori rather than searched. This avoids conflating optimisation budget with modelling capacity: if epochs is searched as a free parameter, the search tends to select the largest budget available. We fix the deep-model training loop at 5000 epochs with early-stopping patience of 200 validation steps, which is sufficient for every retained configuration to converge on this data. Search algorithm and trial budget. We use Optuna [Akiba et al., 2019] with the Tree-structured Parzen Estimator (TPE) sampler and seed 42. For the two deep models, we additionally enable Optuna’s MedianPruner (5 startup trials, 20 warmup epochs, interval 1), which terminates a trial early if its validation RMSE at the current epoch lies above the median of prior trials at the same epoch. Linear and tree models are run without pruning. Trial budgets reflect per-trial cost: linear models receive 2000 trials, tree ensembles 500, and deep models 200 completed or pruned trials. Extending the search to 2× the budget did not improve any of the retained configurations. Search spaces. Table 3 reports the per-model search spaces used in the reported runs. Ranges are deliberately generous rather than tight around prior beliefs. Log-scale ranges are indicated by “log”; step-discretised integer ranges by “step s”. Selected configurations and reproducibility. The selected configuration for each model— including tuned modelling hyperparameters and fixed convergence controls—is included in the released supplementary code. Reproducing Section 5 therefore requires only the released configurations and does not require re-running the hyperparameter search. Re-running the tuning protocol with the same validation split and random seed reproduces the selected configurations up to the usual library- and hardware-level numerical variation. Total tuning cost is approximately 3 hours for the non-deep models and 6 hours for ResNet and FT-Transformer combined; hardware details are given in the Compute paragraph below. Foundation-model baselines. TabPFN and TabICL are applied with their realised library defaults. The expanding-window online variants are rerun at each test step using the available historical sample, with only data observed before the prediction. We do not tune task-specific hyperparameters for these models, preserving their role as strong out-of-the-box baselines. B.4
EWLS implementation and covariance inflation
Algorithm 1 differs from the idealised EWLS recursion in Section 3 in one implementation detail: after the standard forgotten-RLS covariance update, we add a small γ-dependent covariance inflation term εk I. This term is used for numerical stability and responsiveness in finite samples. It is not used to justify the bounded-iterate condition in Appendix A. Cold start for EWLS experts. The diffuse-prior initialisation in Algorithm 1 describes the recursive update after the EWLS state has been initialised. In the experiments we use a short cold-start period of M + 5 test observations for the EWLS experts. During this period, each EWLS correction expert outputs the uniform average of the raw base predictions, M
(γ )
y̌t k =
1 X (j) f (xt ), M j=1
t ≤ M + 5,
while MLpol treats these predictions like ordinary expert predictions and updates its cumulative pseudo-regrets normally. After observing the first M + 5 labels, for each γk we initialise each EWLS state by solving the corresponding batch EWLS problem on those cold-start observations, obtaining (k) (k) the pre-update coefficient wM +6 and covariance PM +5 . From the round M + 6 onward, the EWLS experts use the RLS recursion in Algorithm 1. Thus the cold start affects only the internal EWLS states; the MLpol layer is active from the first test day and all updates remain causal. 24
Table 3: Per-model Optuna search spaces used for base-model hyperparameter tuning. “log”: loguniform sampling. “step s”: discretised with stride s. Convergence controls are fixed a priori and not searched. Model
Hyperparameter
Range −4
Notes 4
Ridge
alpha
[10
XGBoost
n_estimators learning_rate max_depth min_child_weight subsample colsample_bytree reg_alpha reg_lambda
[100, 2000] [10−3 , 0.3] [2, 10] [10−2 , 20] [0.3, 1.0] [0.3, 1.0] [10−8 , 10] [10−8 , 10]
step 100 log integer log uniform uniform log log
LightGBM
n_estimators learning_rate num_leaves max_depth min_child_samples subsample colsample_bytree reg_alpha reg_lambda
[100, 2000] [10−3 , 0.3] [7, 255] [−1, 12] [5, 100] [0.5, 1.0] [0.5, 1.0] [10−8 , 10] [10−8 , 10]
step 100 log log integer (−1 = unlimited) step 5 uniform uniform log log
ResNet
n_layers d d_hidden_factor hidden_dropout residual_dropout d_embedding lr batch_size
[1, 8] [64, 512] [1.0, 4.0] [0, 0.5] [0, 0.5] [4, 64] [10−5 , 10−2 ] {128, 256, 512}
integer step 64 uniform uniform uniform step 4 log categorical
FT-Transformer
n_layers d_token n_heads d_ffn_factor attention_dropout ffn_dropout residual_dropout lr
[1, 6] {64, 128, 192, 256} {4, 8} [2/3, 8/3] [0, 0.5] [0, 0.5] [0, 0.2] [10−4 , 10−2 ]
integer categorical categorical uniform uniform uniform uniform log
Covariance inflation with γ-dependent scaling. (k)
Pt
← γk−1
(k) (k) P z̃t z̃ ⊤ P (k) Pt−1 − t−1 (k)t t−1 st
, 10 ]
log
The implemented covariance recursion is ! + εk I,
(k)
st
(k)
= γk + z̃t⊤ Pt−1 z̃t .
(31)
Under the Kalman-filter interpretation of RLS, adding εk I to the covariance corresponds to a small process-noise covariance in the implicit state equation wt+1 = wt + η t , with Cov(η t ) = εk I. It therefore encodes the expert’s prior scale for how much the combination coefficients may move from one step to the next. A single shared εk ≡ ε is unsatisfactory at both ends of the γ-grid. For large γ, the EWLS expert is intended to represent a long-memory, stable correction; too much injected covariance would dominate the slow forgetting and effectively shorten its memory. For smaller γ, past-data contributions are rapidly discounted, and the recursive covariance can become numerically unstable along directions not recently excited. The inflation level should therefore be coupled to the forgetting scale. We use the power-law schedule εk = ε0 (1 − γk )α , α = 1. (32) The exponent α = 1 is a scale-matching choice. P For an EWLS expert with forgetting factor γ, the total discounted mass of past observations is ℓ≥0 γ ℓ = 1/(1 − γ). Consequently, when the feature 25
directions are of comparable scale, the weighted Gram matrix has nominal magnitude (1 − γ)−1 , so its inverse has nominal magnitude (1 − γ). Setting εk = ε0 (1 − γk )α therefore keeps the added covariance on the same nominal inverse-Gram scale across the forgetting grid. Long-memory experts receive less inflation, so that their slow scale is not overwhelmed; short-memory experts receive more inflation, which helps numerical stability when recent data provide limited curvature. This argument is only a scale heuristic for the implementation. The pathwise guarantees in Appendix A are proved for the idealised recursion without covariance inflation, and ε0 is chosen by walk-forward validation on 2018 (Appendix B.5). Relation to the theoretical analysis. The deterministic analysis in Section 4 and Appendix A studies the idealised forgotten-RLS recursion without the εk I term. The inflation in (31) changes the exact recursion and therefore the calibrated mapping between γ and a physical effective sample size. It should not be read as a proof of the bounded-iterate condition (B4): adding εk I to the covariance caps overconfident precision directions but does not by itself impose a uniform lower bound on λmin (At ). Accordingly, h(γ) = 1/(1 − γ) is treated throughout the experiments as a nominal forgetting-scale index: smaller h corresponds to more responsive corrections and larger h to more stable corrections, but the absolute value of h should not be interpreted as an exact number of days of memory. B.5
Selection of ε0
The covariance-inflation scale ε0 is the only EWLS implementation parameter selected by validation. The exponent α = 1 and the forgetting-factor grid are fixed a priori. We select ε0 by a walk-forward validation protocol on a window preceding the main test period; the main test period (2019-01-01 onwards) is never inspected during selection. Protocol We run a separate walk-forward selection pipeline with all data boundaries shifted one year earlier than in the main experiment. All training and initialisation stages of this pipeline are isolated from 2018; 2018 is used only as the final out-of-sample validation window: • Base-model hyperparameters. The base-model hyperparameters are re-tuned from scratch under the shifted protocol, using 2012–2016 as the training window and 2017 as the tuning validation window. They are not reused from the main run of Appendix B.3, preventing 2018 information from entering through hyperparameter choices. • Base-model fits. The selected configurations are refit on 2012–2017. Neural networks use a tail-of-2017 split for early stopping, while static models are refit on the full 2012–2017 window. The resulting base predictors have never seen 2018. • Aggregation initialisation. The EWLS and MLpol states are initialised using predictions and labels from 2012–2017 only. • Validation window for ε0 . Candidate values are scored on 2018-01-01 to 2018-12-31, a 365-day window of genuine out-of-sample predictions for every base model. During 2018 the online updates follow the same predict-then-update protocol as in the main evaluation. Thus, no 2018 information enters the base-hyperparameter search, base-model fit, or aggregation initialisation. The 2018 RMSE is used only to choose ε0 . For each ε0 in the logarithmic grid {10−13 , 10−12 , . . . , 10−5 }, with α = 1 fixed, we run the full EWLS + MLpol pipeline and record the RMSE on the 2018 validation window. The minimiser is then fixed for the main run and is not revised using the 2019+ test period. The selection script enforces that the validation window cannot end after 2018-12-31, thereby preventing 2019+ dates from entering the selection. Sweep curve
Figure 4 shows the resulting validation RMSE as a function of ε0 .
Diagnostics • Interior minimum, both sides bracketed. ε0 = 10−8 is not at a grid boundary. Performance worsens when ε0 is made much smaller than the useful range (ε0 ≤ 10−11 ), and the curve rises again for larger values, with a clear deterioration by ε0 = 10−5 . Thus the grid brackets the useful scale rather than selecting an endpoint. 26
RMSE on 2018 out-of-sample (MW)
700
680
660 10−13 10−12 10−11 10−10 10−9 10−8 10−7 10−6 10−5 ε0 (log scale) Base-only (705.3899)
MLpol on Base+EWLS
Selected ε0 = 10−8
Figure 4: Walk-forward ε0 sweep on 2018 out-of-sample. The validation curve has an interior minimum at ε0 = 10−8 (RMSE 677.3). The shaded region marks a conservative two-decade plateau [10−9 , 10−7 ], where the RMSE remains within 1% of the best value. The neighbouring point at 10−6 is also close, but the curve has already started to rise. The left tail worsens when the inflation is too small relative to the conditioning of the EWLS covariance recursion, while the right tail worsens when excessive inflation prevents high-γ experts from retaining their intended long-memory behaviour. The main test period, starting on 2019-01-01, is not used in this selection. • Wide plateau, robust selection. The three grid points {10−9 , 10−8 , 10−7 }, spanning two orders of magnitude, are all within 1% of the best RMSE. The neighbouring value 10−6 is also close, but the upward trend has already begun. The validation procedure therefore identifies a stable scale rather than a finely tuned scalar, mitigating concerns that ε0 is overfit to the particular 2018 validation year. • EWLS layer helps even in a non-crisis validation year. MLpol on Base only, an ε0 -independent reference, achieves RMSE 705.4 on 2018. At the selected ε0 = 10−8 , the full Base+EWLS pool reaches 677.3, a 28.0 MW (4.0%) improvement even in a year without a COVID-scale regime shift. This suggests that the EWLS layer can add value beyond the base experts alone, even outside extreme regime shifts. Why we do not use in-sample validation An alternative validation protocol is to sweep ε0 using in-sample base predictions on the training period (2012–2017), evaluated on the last year of training. We report such a sweep for completeness: it prefers the smallest ε0 in the grid and does not produce an interior minimum over [10−13 , 10−5 ]. This behaviour is consistent with the theoretical role of ε0 : in-sample residuals can under-represent the out-of-sample variation that an EWLS expert must track, biasing the sweep toward insufficient injection. The walk-forward protocol above avoids this failure mode by evaluating on genuine out-of-sample predictions. Agreement with an unrestricted test-period sweep As a post-hoc robustness check, not used for selection, we also swept ε0 directly on the main test period under the same diagnostic setup. The absolute RMSE values from this diagnostic run are not meant to replace the final main-experiment numbers; they are reported only to compare the selected and test-optimal ε0 under a common setup. The test-optimal value ε0 = 10−9 attains overall RMSE 660.57, while the walk-forward-selected value ε0 = 10−8 attains 661.61. The 1.04 MW (0.16%) gap is within the plateau already identified on 2018, and quantifies the small protocol cost of selecting ε0 without touching the test period. B.6
Compute and reproducibility
Compute. All experiments were run on a single Apple M4 Pro workstation (14-core CPU, 20-core integrated GPU, 48 GB unified memory). Deep-model training (ResNet, FT-Transformer) used the 27
PyTorch MPS backend; all other components ran on CPU. Base-model hyperparameter tuning took approximately 3 hours for the linear and tree ensembles combined and approximately 6 hours for ResNet and FT-Transformer combined. The main online-aggregation pipeline over the 746-day test period with N = 23 experts runs end-to-end in under 10 s on CPU—the recursion is dominated by M × M linear-algebra updates and is not the computational bottleneck. The paired moving block bootstrap (Appendix C.2, R = 10,000) runs in under 30 s. B.7
Empirical bounded-iterate diagnostic
The pathwise tracking guarantee in Theorem 1 (restated as Theorem A.3 in Appendix A) is conditional (γ) on the bounded-iterate condition (B4): ∥wt ∥2 ≤ R uniformly along the realised sequence. This appendix provides an empirical diagnostic for this condition on RTE-FR Load. Decomposing slope and intercept. The EWLS expert is fitted on the augmented base-prediction vector zet = (zt⊤ , 1)⊤ , so that (γ)
wt
(γ)
(γ)
= βt , αt
(γ)
,
(γ)
where βt ∈ RM combines the base forecasts and αt ∈ R is an intercept. These two components have different scales. The slope vector is the online combination component: large oscillations in (γ) ∥βt ∥2 would indicate unstable adaptation of the base-forecast weights. The intercept, in contrast, absorbs the level of the response. Since daily French electricity load is measured in megawatts (γ) (γ) and is of order 104 –105 , the augmented norm ∥wt ∥2 can be dominated by |αt | even when the combination weights remain stable. Relation to the bounded-iterate assumption. Condition (B4) is a sufficient condition used to make the pathwise bound finite. In the augmented parametrisation, its numerical value mixes two effects: the scale-free stability of the combination slope and the scale-dependent magnitude of the intercept. The diagnostic below therefore reports both the full augmented norm and the slope-only norm. The former checks the literal quantity appearing in (B4), while the latter isolates the part of the iterate that governs online aggregation of base forecasts. Empirical result. Figure 5 shows the trajectories for the K − 1 = 15 finite-forgetting EWLS experts over the 2019–2021 test period. The no-forgetting endpoint γ = 1 is excluded because Theorem 1 is stated for finite forgetting factors. Two observations are useful. First, the full augmented norm is mainly driven by the intercept. Its running supremum is 1.50 × 104 , which is on the same scale as the response level in megawatts. This behaviour reflects the scale of the intercept. Second, after removing the intercept coordinate, the slope norm remains small: (γ)
sup ∥βt ∥2 = 1.69 t,γ
over all finite-forgetting experts and all T = 746 test days. The fastest-forgetting experts respond most strongly around the onset of the COVID-19 lockdown, where the running supremum is attained, while the long-memory experts change more gradually. After the lockdown period, the slope norm relaxes to a range around 0.5–0.7 across most of the grid. Overall, this diagnostic supports the plausibility of the bounded-iterate condition along the realised RTE-FR Load trajectory: the literal augmented iterate remains finite, and the scale-free combination component remains uniformly bounded by a small constant. We therefore use this figure as a sanity check for the conditional pathwise guarantee, rather than as an independent proof of (B4). 28
With intercept 15000
running sup ‖w‖2 = 14983.37 s≤t, γ
5000
10000 7500 5000 2500
1000
0
500 Slope only
1.75 1.50
running sup ‖β‖2 = 1.69 s≤t, γ
100
1.25 ‖βγ(t)‖2
h(γ) = 1/(1 − γ)
‖wγ(t)‖2
12500
50
1.00 0.75
20
0.50 0.25 0.00 Dec 2018
Apr 2019
Aug 2019
Dec 2019
Apr 2020
Aug 2020
Dec 2020
Figure 5: Empirical bounded-iterate diagnostic on RTE-FR Load. Trajectory of the EWLS coefficient norm for the K − 1 = 15 finite-forgetting experts over the 2019–2021 test period. The no-forgetting endpoint γ = 1 is outside the finite-forgetting setting of Theorem 1 and is excluded. Curves are coloured by effective memory h(γ) = 1/(1 − γ), with short-memory experts in warm colours and long-memory experts in cool colours; the shaded band marks the COVID-19 lockdown (γ) period. Top: full augmented norm ∥wt ∥2 . The running supremum is 1.50 × 104 , reflecting the (γ) scale of the intercept. Bottom: slope-only norm ∥βt ∥2 , with the intercept coordinate removed. The running supremum is 1.69, attained near the lockdown onset. Thus the EWLS slope iterate responds to the structural break while remaining uniformly bounded across forgetting factors and test days.
C
Detailed RTE-FR results and statistical significance
C.1
Full per-regime RMSE table
Table 4 reports the full per-regime RMSE table behind Table 1. The rows include the final aggregates, all standalone EWLS/RLS correction experts, external foundation-model references, and all individual base-pool models. The column “Excl. warm-up” reports the same overall test metric after removing the initial M + 5 cold-start observations used to initialise the EWLS states; all methods are evaluated on the same remaining dates in that column. C.2
Bootstrap uncertainty for the main comparison
Because MLpol is deterministic given its inputs and the base experts are fit once on a fixed training window, Table 1 contains no seed-level variability to report. The uncertainty that remains is sampling variability of the finite test trajectory itself: had the test window been a different but comparable trajectory of comparable length, how would the headline numbers move? We quantify this via a paired moving block bootstrap on the per-step squared-loss sequences. Protocol. For each method we treat the sequence of daily squared errors {(ŷt −yt )2 }Tt=1 as the basic observation and resample with a moving block bootstrap [Künsch, 1989]. Within each regime (overall, pre-lockdown, lockdown, post-lockdown) we draw ⌈Tregime /B⌉ block start positions uniformly from the valid start set, concatenate the resulting length-B blocks, and truncate to the regime length. The replicate RMSE is computed from the resampled sequence. We use B = 14 days on overall, 29
Table 4: Full per-regime RMSE (MW) on RTE-FR Load. Rows are sorted by overall RMSE. “Excl. warm-up” removes the initial M + 5 EWLS cold-start observations from the overall evaluation window. External online TabICL/TabPFN variants are rerun daily on an expanding context window; “+GRI” denotes the additional regime-informed Oxford Government Response Index covariate. Group
Method
Pre-lockdown Lockdown
Post-lockdown Excl. warm-up Overall
Final agg. Final agg.
MELO (Base+EWLS) MLpol on EWLS only
623.1 653.9
1086.1 1073.6
579.3 584.1
651.5 656.9
655.8 673.2
EWLS expert EWLS expert EWLS expert EWLS expert EWLS expert EWLS expert EWLS expert EWLS expert EWLS expert EWLS expert EWLS expert EWLS expert
EWLS, γ = 0.998563 EWLS, γ = 0.997868 EWLS, γ = 0.999031 EWLS, γ = 0.996838 EWLS, γ = 0.995309 EWLS, γ = 0.999347 EWLS, γ = 0.993041 EWLS, γ = 0.999560 EWLS, γ = 0.989676 EWLS, γ = 0.984685 EWLS, γ = 0.999703 EWLS, γ = 0.977280
657.9 661.5 656.2 666.9 673.9 656.6 682.2 659.0 691.2 700.6 663.3 710.9
1168.3 1138.2 1204.0 1110.9 1084.1 1247.8 1056.1 1301.9 1028.3 1006.9 1367.9 998.8
585.4 592.8 579.9 602.1 613.2 576.1 625.3 574.0 637.3 648.2 573.6 659.2
673.0 673.0 675.6 675.2 679.1 681.1 684.0 689.9 689.7 696.1 702.2 704.7
687.5 687.9 689.7 690.3 694.4 694.7 699.5 702.8 705.2 711.7 714.4 720.2
External EWLS expert EWLS expert EWLS expert External Final agg. External External EWLS expert
TabICL (online, +GRI) EWLS, γ = 0.999800 EWLS, γ = 0.966295 EWLS, γ = 0.950000 TabPFN (online, +GRI) MLpol on Base only TabICL (online) TabPFN (online) EWLS, γ = 1.000000
668.0 669.2 723.0 738.4 670.9 690.5 670.3 669.3 736.7
1132.6 1446.3 1009.1 1037.9 1165.5 2452.7 2943.6 2943.3 3212.4
699.5 575.2 673.6 695.5 788.4 907.0 687.4 719.7 878.5
710.2 718.2 717.8 737.4 749.3 1013.0 1038.6 1047.8 1170.2
723.4 729.5 732.7 751.4 758.6 1004.0 1036.2 1043.1 1163.2
Base pool External Base pool External Base pool Base pool Base pool Base pool Base pool
GAM TabPFN (offline) ResNet TabICL (offline) XGBoost LightGBM FT-Transformer Ridge Lag-1
946.6 677.7 771.2 683.5 842.8 821.9 1060.4 1999.0 3918.6
3201.3 4103.1 3582.5 4453.8 3972.8 4410.9 4131.1 3116.8 2626.3
1076.2 1055.8 1336.7 1137.9 1347.4 1544.3 1763.3 1918.7 3286.4
1314.2 1396.5 1401.4 1500.2 1510.7 1656.9 1744.8 2067.4 3574.4
1298.3 1381.1 1382.5 1482.3 1486.8 1629.7 1727.3 2078.5 3633.2
pre-lockdown, and post-lockdown (two weekly cycles), and B = 7 on the 56-day lockdown regime so that the set of valid block starts remains non-trivial. All methods share the same block-start draws within a replicate (paired bootstrap), so the distribution of pairwise RMSE differences is well-defined. Results below use R = 10,000 replicates and seed 0; Monte Carlo variation in the reported CIs is under 1 MW (verified against independent runs at R = 1,000). The bootstrap consumes the same per-step squared-loss sequences that produce Table 1; the point estimates therefore match the main results exactly. Given a method m and anchor method m0 (our full Base+EWLS aggregate), we report: • the point estimate RMSEm (identical to Table 1) and the 2.5–97.5% percentile bootstrap CI, as a marginal uncertainty on m’s performance; • the paired difference ∆m = RMSEm − RMSEm0 , computed on the real data, and the 2.5–97.5% percentile CI of the bootstrap distribution of the difference. A CI excluding zero corresponds to a significant difference at the 5% level under this resampling scheme.
Per-method RMSE with 95% bootstrap CIs. Table 5 reports, for each method and regime, the point estimate and the 95% bootstrap CI. The CI widths scale as expected with regime length (narrow on the 441-day pre-lockdown, very wide on the 56-day lockdown) and with RMSE level (wider for the offline baselines that are dominated by a few catastrophic lockdown days). 30
Table 5: Per-method RMSE (MW) with 95% moving-block bootstrap CIs (R = 10,000, B = 14 for overall/pre/post, B = 7 for lockdown). Point estimates match Table 1 exactly. CIs are rounded to the nearest MW for readability. The same bootstrap resamples are shared across methods and are used for the paired differences in Table 6. Pre-lockdown
Lockdown
Post-lockdown
Method
RMSE 95% CI
RMSE 95% CI
RMSE 95% CI
MLpol on Base+EWLS (ours) MLpol on EWLS only MLpol on Base only Best EWLS (γ ≈ 0.9986, hindsight) TabICL (online, +GRI) TabPFN (online, +GRI) TabICL (online) TabPFN (online) TabICL (offline) TabPFN (offline)
623.1 653.9 690.5 657.9 668.0 670.9 670.3 669.3 683.5 677.7
1086.1 1073.6 2452.7 1168.3 1132.6 1165.5 2943.6 2943.3 4453.8 4103.1
[532, 680] [546, 704] [595, 752] [552, 711] [564, 718] [568, 723] [566, 715] [566, 711] [581, 727] [575, 721]
Overall RMSE 95% CI
[600, 1323] 579.3 [466, 661] [673, 1315] 584.1 [473, 666] [1654, 2984] 907.0 [740, 1026] [867, 1367] 585.4 [475, 665] [579, 1425] 699.5 [544, 854] [522, 1452] 788.4 [637, 936] [1914, 3742] 687.4 [556, 790] [1854, 3714] 719.7 [559, 828] [3847, 5176] 1137.9 [810, 1333] [3244, 4953] 1055.8 [744, 1213]
655.8 673.2 1004.0 687.5 723.4 758.6 1036.2 1043.1 1482.3 1381.1
[565, 745] [576, 749] [780, 1245] [591, 761] [618, 816] [644, 865] [681, 1390] [683, 1388] [972, 1954] [899, 1842]
Pairwise differences vs. our full method. Table 6 is the more informative summary: each row reports ∆ = RMSEm − RMSEours on the real data, and the 95% CI of its bootstrap distribution. A checkmark in the sig. column indicates that the CI excludes zero. The paired structure concentrates the uncertainty on the difference, which is considerably narrower than the individual marginal CIs in Table 5 because the methods share common noise on each day. Table 6: Paired ∆RMSE = RMSEm − RMSEours (MW) with 95% moving-block bootstrap CIs. Positive ∆ favours our full Base+EWLS aggregate. ✓ marks CIs excluding zero. R = 10,000, with block lengths as in Table 5. Panel A: Pre-lockdown and lockdown. Pre-lockdown Comparator m
∆
MLpol on EWLS only MLpol on Base only Best single EWLS (γ ≈ 0.9986, hindsight) TabICL (online, +GRI) TabPFN (online, +GRI) TabICL (online) TabPFN (online) TabICL (offline) TabPFN (offline)
30.8 [ 67.4 [ 34.8 [ 44.9 [ 47.8 [ 47.2 [ 46.2 [ 60.4 [ 54.6 [
95% CI
sig.
−3.5, 50.4] 27.5, 112.1] ✓ 1.5, 55.3] ✓ −31.0, 95.0] −21.7, 92.6] −31.2, 95.1] −28.3, 87.4] −18.4, 108.5] −19.3, 96.6]
Lockdown ∆
Panel B: Post-lockdown and overall. Post-lockdown Comparator m
∆
MLpol on EWLS only MLpol on Base only Best single EWLS (γ ≈ 0.9986, hindsight) TabICL (online, +GRI) TabPFN (online, +GRI) TabICL (online) TabPFN (online) TabICL (offline) TabPFN (offline)
4.8 [ 327.8 [ 6.2 [ 120.2 [ 209.1 [ 108.2 [ 140.4 [ 558.6 [ 476.6 [
95% CI
sig.
−4.8, 13.8] 197.2, 442.1] ✓ −14.1, 25.5] 7.7, 269.6] ✓ 118.3, 340.5] ✓ 25.7, 208.3] ✓ 43.7, 219.3] ✓ 250.1, 774.5] ✓ 198.0, 649.6] ✓
95% CI
sig.
−12.5 [ −30.5, 112.2] 1366.6 [ 993.3, 1772.5] ✓ 82.2 [ 16.5, 316.7] ✓ 46.5 [ −154.0, 268.4] 79.4 [ −225.1, 216.1] 1857.5 [ 1237.3, 2545.8] ✓ 1857.2 [ 1191.0, 2485.5] ✓ 3367.7 [ 3064.1, 4042.8] ✓ 3017.1 [ 2543.0, 3774.0] ✓ Overall ∆ 17.4 [ 348.3 [ 31.7 [ 67.7 [ 102.8 [ 380.4 [ 387.3 [ 826.6 [ 725.3 [
95% CI
sig.
−7.0, 32.9] 185.9, 527.6] ✓ −2.5, 56.7] 9.4, 125.3] ✓ 39.7, 160.9] ✓ 86.4, 675.2] ✓ 90.8, 667.1] ✓ 363.5, 1254.6] ✓ 288.5, 1137.8] ✓
What is and is not statistically established. We summarise the findings of Table 6; all inferences are at the 5% level under the paired moving block bootstrap. 1. The EWLS correction layer adds significant value in every regime. Against MLpol on Base only, Ours wins by +67.4 MW on pre-lockdown, +1366.6 MW on lockdown, +327.8 MW on post-lockdown, and +348.3 MW overall, with every CI strictly positive. This is the cleanest empirical support of Section 4 and the statistical basis for the headline finding in Section 5.1. 31
2. Base+EWLS beats the hindsight-optimal single-γ EWLS expert in pre-lockdown and lockdown; the overall improvement is positive but marginally inside bootstrap noise. Against the best single EWLS expert chosen in hindsight on overall RMSE (γ ≈ 0.9986, overall RMSE 687.5), our full pool wins by +34.8 MW on pre-lockdown (CI [1.5, 55.3], ✓) and +82.2 MW on lockdown (CI [16.5, 316.7], ✓); the overall improvement is +31.7 MW (CI [−2.5, 56.7], marginally inside bootstrap noise) and the post-lockdown difference (+6.2 MW, CI [−14.1, 25.5]) is essentially zero. This is the expected shape: no single memory scale dominates the whole trajectory. The overall-selected single γ remains locally competitive post-lockdown, but pays in regimes where the optimal effective memory differs. Per-regime hindsight optima make this point quantitatively: the best γ is ≈ 0.9990 on prelockdown, ≈ 0.977 on lockdown, and ≈ 0.9997 on post-lockdown—an order-of-magnitude shift in optimal effective memory across the three regimes. In pre-lockdown, the full Base+EWLS aggregate (623.1) in fact beats every single EWLS expert in hindsight, including the per-regime oracle, because the base-expert channel carries predictive structure orthogonal to any memorylength choice. 3. Base experts contribute a small orthogonal signal, most visibly in pre-lockdown. Adding base experts to MLpol on EWLS-only yields a positive pre-lockdown gain (∆ = +30.8 MW, CI [−3.5, 50.4], marginally inside bootstrap noise); the overall gap (+17.4 MW, CI [−7.0, 32.9]) and the post-lockdown gap (+4.8 MW, CI [−4.8, 13.8]) are favourable but within bootstrap noise, while the lockdown comparison marginally favours EWLS-only (−12.5 MW, CI [−30.5, 112.2], not significant). The picture is coherent with the design: stable historical structure—which base experts encode and EWLS cannot recover from a fixed memory grid alone—contributes most where the base models’ training distribution applies, i.e. pre-lockdown. Under the structural break, EWLS-only and Base+EWLS are statistically indistinguishable, confirming that the lockdown gain in Table 1 is delivered by the EWLS layer rather than by the base pool. 4. Significantly beats TabICL+GRI overall and post-lockdown, despite using no external regime-response signal. Against TabICL (online, +GRI)—the strongest competitor, supplied with a regime-informed covariate—Ours wins significantly on overall RMSE (∆ = +67.7 MW, CI [9.4, 125.3], ✓, a 9.4% relative improvement) and in post-lockdown (∆ = +120.2 MW, CI [7.7, 269.6], ✓, a 17.2% relative improvement). The pre-lockdown and lockdown margins are within bootstrap noise. This is consistent with the role of the GRI covariate: before the lockdown shock it provides little regime-shift information, while during lockdown it gives TabICL+GRI an immediate external response signal. The recovery regime is where this coincident signal is less directly informative and our method’s memory-scale tracking takes over. 5. Significantly beats every no-GRI online adaptive baseline in the regimes where adaptation is stress-tested. Against TabICL/TabPFN online (without GRI), every difference is significant in lockdown, post-lockdown, and overall; pre-lockdown is the lone null cell, where no method actually exercises adaptation and all reasonable predictors lie within ∼ 50 MW of each other. Sensitivity to block length. We verified the robustness of the overall-regime CIs to block length by repeating the bootstrap at B ∈ {7, 14, 28}. Point estimates are identical across B by construction. The width of each CI increases with B (as fewer independent blocks fit inside the trajectory), but the qualitative sig. column of Table 6 is stable—no entry flips between ✓ and blank over this range. Implementation. The bootstrap procedure is implemented as scripts/bootstrap_ci.py in the released code base and runs in under 30 s for R = 10,000 on a laptop; it consumes only the saved aggregation outputs and requires no pipeline re-execution.
D
Adaptive-filter and state-space baselines
This section clarifies the relationship between the EWLS/RLS update used in Algorithm 1 and Kalman filtering. The purpose is interpretive: our method does not attempt to estimate a latent process-noise covariance. Instead, each forgetting factor defines a fixed time-scale objective, and the outer MLpol layer selects among these objectives according to realised forecasting loss. Consider first the ideal EWLS recursion without the inflation term. We write wt|t−1 and Pt|t−1 for the state mean and covariance used before observing yt , and wt|t and Pt|t for the corresponding 32
quantities after observing yt . In the forgotten-RLS interpretation, Pt|t−1 = γ −1 Pt−1|t−1 ,
wt|t−1 = wt−1|t−1 .
With unit observation-noise scale, the measurement update is Kt = Pt|t−1 z̃t St−1 , wt|t = wt|t−1 + Kt yt − z̃t⊤ wt|t−1 ,
St = 1 + z̃t⊤ Pt|t−1 z̃t ,
Pt|t = Pt|t−1 − Kt z̃t⊤ Pt|t−1 .
Equivalently, since Pt|t−1 = γ −1 Pt−1|t−1 , this can be written in the standard RLS form Pt−1|t−1 z̃t , st wt|t = wt−1|t−1 + Kt yt − z̃t⊤ wt−1|t−1 , ! Pt−1|t−1 z̃t z̃t⊤ Pt−1|t−1 −1 Pt|t = γ Pt−1|t−1 − . st st = γ + z̃t⊤ Pt−1|t−1 z̃t ,
Kt =
The through-t coefficient wt|t is the minimiser of t X s=1
γ t−s (ys − w⊤ z̃s )2 + γ t δ0 ∥w∥22 .
Thus the coefficient used to predict at round t is wt|t−1 , while wt|t is available only after yt has been observed. This recursion can be read as a Kalman filter for the linear observation model yt = z̃t⊤ wt + εt , with random-walk state equation wt = wt−1 + ηt . In this interpretation, forgetting replaces the previous posterior covariance Pt−1|t−1 by the prior covariance Pt|t−1 = γ −1 Pt−1|t−1 . Hence the implicit process-noise covariance is (γ)
Qt
= Pt|t−1 − Pt−1|t−1 = (γ −1 − 1)Pt−1|t−1 .
Thus smaller γ corresponds to a larger prior covariance inflation and therefore to a filter that is willing to move its coefficients more rapidly. Conversely, γ close to one imposes a long-memory, slowly moving coefficient path. The implementation adds a small isotropic inflation term after the measurement update, Pt = Pt|t + εk I. Under the Kalman interpretation, this can be viewed as carrying a small additional process-noise component into the next round. Its role in our experiments is numerical stabilisation, especially for fast-forgetting experts whose effective design matrices may become poorly conditioned. This interpretation is different from adaptive Kalman filtering. Methods such as adaptive Kalman filters try to infer or update a process-noise covariance from innovations, likelihoods, or variational criteria. Our construction does not estimate Qt . Instead, it defines a finite family of fixed EWLS objectives indexed by γ, equivalently a finite family of implicit time-scale assumptions. The choice among these assumptions is deferred to the outer MLpol aggregation layer, which is driven directly by realised squared forecasting loss. In this sense, the forgetting factors are not latent parameters to be estimated inside a single state-space model; they are competing online objectives exposed to expert aggregation. 33
D.1
RLS, Kalman filtering, and time-scale modelling
A natural question is whether the separation of roles in our framework—a fixed γ-grid of EWLS time scales at the expert layer, followed by loss-driven MLpol aggregation—could instead be replaced by a Kalman-filter-style adaptive layer that updates process-noise uncertainty online. We therefore compare against three representative adaptive-filter alternatives on the French electricity benchmark. The goal of this appendix is not to rule out all possible adaptive state-space designs, but to test whether standard adaptive-Q mechanisms, tuned under the same validation protocol, provide the same empirical behaviour as the γ-grid + MLpol construction. We benchmark γ-grid + MLpol against three adaptive-KF baselines that cover complementary modelling choices: • VIKING [de Vilmarest and Wintenberger, 2024]: a variational Kalman filter with isotropic process noise Qt = ebt I. A Gaussian posterior over the log-variance bt is updated online, and a scalar parameter ρb controls the variability of this log-variance process. • VBAKF [Huang et al., 2018]: a variational Bayesian adaptive Kalman filter. An inverseWishart prior is placed on the predicted state covariance and an inverse-gamma prior on the scalar observation noise; a forgetting factor ρ ∈ (0, 1] controls how quickly the covariance posterior discounts past information. • IMM-KF [Blom and Bar-Shalom, 1988]: an interacting multiple-models filter. It maintains several Kalman filters on a fixed grid of Qi = qi I, mixes their state through a Markov transition matrix, and updates mode probabilities by Gaussian innovation likelihoods. IMM is the closest state-space analogue of a multi-scale construction, but its combination rule is Bayesian Markov mixing rather than regret minimisation over predictions exposed as separate experts. These baselines vary along two axes: the parameterisation of process uncertainty (isotropic scalar, full covariance, or discrete grid) and the adaptation rule used to update and combine that uncertainty. This makes them a useful stress test for whether the explicit expert-layer / aggregation-layer separation is empirically useful on this benchmark. Setup and tuning protocol. For a like-for-like comparison, each adaptive-KF baseline is paired with the same M = 7 direct base experts and combined by the same outer MLpol aggregator as our main method. Only the adaptive expert-layer component varies. For VIKING and VBAKF, the responsiveness axis is swept as a small pool for MLpol to aggregate over. IMM-KF already contains an internal grid of Q scales, so we use a single IMM instance with that grid: • VIKING pool: three instances at ρb ∈ {10−4 , 10−3 , 10−2 }, spanning conservative → responsive log-variance random-walk rates. • VBAKF pool: three instances at ρ ∈ {0.95, 0.98, 0.995}, spanning responsive → conservative IW-concentration forgetting. • IMM: a single instance with a five-mode grid q ∈ {10−8 , 10−6 , 10−4 , 10−2 , 100 }. The IMM architecture already integrates a grid of Q scales internally, so a pool of multiple IMMs would be redundant. All remaining nuisance hyperparameters (VIKING’s b̂0 , Σ0 , VB iteration count, and initial covariance scale; VBAKF’s τP , τR , VB iterations, and δ; IMM’s πstay and δ) are tuned on the 2018 out-of-sample window via the same walk-forward protocol used for our method’s ε0 selection (Appendix B.5). In this protocol, base models are fit on 2012–2016, the aggregator runs online over 2017–2018, and the 2018 RMSE of the corresponding MLpol pipeline is minimised over a grid of nuisance-parameter combinations. The real test period (2019-01-01 onward) is never used for selection. Thus all adaptivefilter baselines are tuned under the same information constraints as our method. Observation-noise priors are initialised from the training residual variance of the uniform base-expert mean and are not tuned. Selected hyperparameters are listed in Table 7. Main comparison. Table 8 reports RMSE by regime. Under the matched walk-forward tuning protocol, MELO has the lowest overall RMSE among all compared variants. Relative to the adaptiveKF pools that are also exposed to the same direct base experts and the same outer MLpol aggregator, 34
Table 7: Nuisance-parameter search grids and selected values for the three adaptive-KF baselines. All values selected by minimising the 2018-out-of-sample RMSE of MLpol on Base+pool, using the walk-forward protocol of Appendix B.5. The responsiveness-axis grids (VIKING ρb , VBAKF ρ, and IMM’s internal q-grid) are listed in the Setup of §D and are not tuned: they play the role of the pool that MLpol aggregates over. Observation-noise priors (â0 for VIKING, R0 for VBAKF/IMM) are initialised from train-residual variance of the uniform base-expert mean and are also not tuned. Selected values at grid boundaries were confirmed by boundary-extension sweeps. Parameter
Role
Search grid
Selected
VIKING (fixed: learn_sigma = learn_Q = True, ρa = 0) b̂0 prior mean of log Q {−14, −12, −10, −8, −6, −4, −2, 0} Σ0 prior variance of log Q {10−5 , 10−4 , 10−3 , 10−2 , 10−1 , 1} NVB Newton steps per observation {1, 2, 3, 5, 8} δ initial covariance P0 = δI {10−3 , 10−2 , 10−1 , 1, 10, 102 }
−8 10−1 8 10−1
VBAKF τP τR NVB δ
prior IW dof on P prior IG dof on R variational iterations/step initial P0 = δI
{3, 4, 5, 7, 10, 15, 25, 50} {3, 4, 5, 7, 10, 15, 25, 50} {3, 5} {10−2 , 10−1 , 1, 10, 102 , 103 , 104 }
3 50 3 10−2
IMM-KF πstay
Markov self-transition prob.
δ
initial P0 = δI per mode
{0.30, 0.50, 0.60, 0.70, 0.80, 0.90, 0.95, 0.98, 0.99, 0.995, 0.999} {10−3 , 10−2 , 10−1 , 1, 10, 102 }
0.5 10−2
MELO reduces overall RMSE by 3.6% against Base+VIKING, 7.5% against Base+IMM, and 9.3% against Base+VBAKF. These overall gains are also confirmed by the paired moving-block bootstrap reported below. The regime decomposition gives a more nuanced picture. MELO is the best point estimate before lockdown and after lockdown, and remains close to the best adaptive variants during the short lockdown window. In that window, EWLS-only and VBAKF-only have slightly lower point estimates, while IMM-only is essentially tied with MELO. Thus the advantage of MELO in this comparison is not that it uniformly dominates every adaptive filter in every regime, but that it gives the best overall balance across the full 2019–2021 online test period. The filter-only rows provide a diagnostic on how the adaptive-filter outputs interact with the direct base experts. Adding direct experts improves overall RMSE for VIKING, VBAKF, and our EWLS pool, while it worsens the overall RMSE of IMM. The regime decomposition shows that this interaction is not uniform: during lockdown, adding direct experts helps VIKING, slightly hurts the EWLS pool, and hurts IMM and VBAKF more substantially. This suggests that the benefit of exposing direct base forecasts to the outer aggregator depends on whether the adaptive component remains complementary to the base pool during the break. We return to this point in the IMM diagnostic below. Bootstrap significance. We extend the paired moving-block bootstrap of Appendix C.2 to the three adaptive-KF baselines with direct base experts exposed to the same outer MLpol aggregator. We use R = 10,000 bootstrap replicates, with block length B = 14 for the overall, pre-lockdown, and post-lockdown evaluations, and B = 7 for the short lockdown regime. Table 9 reports the paired ∆RMSE = RMSEm − RMSEMELO , where positive values favour MELO, together with the 95% percentile CI of the bootstrap distribution. The bootstrap results support the same overall conclusion as the point estimates. Each of the three adaptive-KF baselines with direct experts is significantly worse than MELO on overall RMSE: the paired difference is +24.7 MW for Base+VIKING (CI [1.9, 51.1]), +52.9 MW for Base+IMM-KF (CI [24.8, 86.5]), and +67.0 MW for Base+VBAKF (CI [32.0, 113.6]). The regime-level results are more mixed. MELO is significantly better than all three baselines after lockdown and significantly better than Base+VBAKF before lockdown, but the shorter 56-day lockdown window has wide uncertainty intervals and none of the three pool-plus-direct lockdown comparisons excludes zero. We therefore interpret the per-regime diagnostics as qualitative evidence about behaviour, and the overall paired bootstrap as the main significance statement. 35
Table 8: Per-regime and overall RMSE (MW) of MELO against three adaptive Kalman-filter baselines on the 2019-01-01 to 2021-01-15 test period. All nuisance hyperparameters of each baseline were selected on 2018 out-of-sample data with the same walk-forward protocol used for our own ε0 selection (Appendix B.5); the real test period was never touched during selection. filter-only rows drop direct base experts from the outer MLpol pool. Bold: best per column; underline: second best. Method
Pre-lockdown
Lockdown
Post-lockdown
Overall
Our method MELO (MLpol on Base+EWLS) MLpol on EWLS only
623.1 653.9
1086.1 1073.6
579.3 584.1
655.8 673.2
Adaptive KF baselines MLpol on Base+VIKING MLpol on VIKING only MLpol on Base+IMM MLpol on IMM only MLpol on Base+VBAKF MLpol on VBAKF only
634.0 676.8 636.5 670.8 665.1 869.3
1185.3 1293.0 1276.7 1086.5 1227.4 1066.9
599.7 595.0 648.6 605.0 665.4 633.6
680.5 717.3 708.7 690.5 722.8 816.2
Reference MLpol on Base only
690.5
2452.7
907.0
1004.0
Table 9: Paired ∆RMSE = RMSEm − RMSEMELO (MW) for the three adaptive-KF baselines with direct base experts included in the outer MLpol pool, with 95% paired moving-block bootstrap CIs. Positive ∆ favours MELO. ✓ indicates that the CI excludes zero. The corresponding adaptive-filteronly ablations are discussed in the text. Pre-lockdown Comparator m
Lockdown
∆ / 95% CI sig.
Base+VIKING +10.9 [−2.2, 25.6] Base+IMM-KF +13.4 [−5.5, 33.2] Base+VBAKF +42.0 [12.0, 76.3]
✓
∆ / 95% CI sig. +99.2 [−144.0, 180.8] +190.6 [−35.5, 315.4] +141.4 [−89.5, 282.3]
Post-lockdown
Overall
∆ / 95% CI sig. +20.5 [5.0, 43.3] +69.4 [16.9, 108.8] +86.2 [25.3, 119.6]
✓ ✓ ✓
∆ / 95% CI sig. +24.7 [1.9, 51.1] +52.9 [24.8, 86.5] +67.0 [32.0, 113.6]
✓ ✓ ✓
The adaptive-filter-only ablations reach the same qualitative verdict on overall RMSE for VIKING and VBAKF: VIKING-only and VBAKF-only are significantly worse than MELO, with paired differences of +61.5 MW (CI [18.9, 92.6]) and +160.4 MW (CI [55.6, 253.1]), respectively. IMMonly is closer: the point estimate still favours MELO (+34.7 MW), but the overall CI crosses zero (CI [−4.0, 52.6]). This is consistent with the diagnostic in §D.2: IMM’s Bayesian mixture remains competitive when isolated, but appears less complementary to the direct experts under the outer MLpol aggregator. The diagnostics below examine why these adaptive-filter alternatives behave differently from the fixed-scale EWLS grid on this benchmark. They should be read as empirical interpretations of this comparison, rather than as general negative results about adaptive Kalman filtering. D.2
IMM diagnostics
IMM is informative because its internal mode probabilities respond visibly to the lockdown transition. Figure 6 plots the mode probabilities µi,t over the full test period and zooms in on the lockdown onset. Mode probabilities respond to the break. Outside the lockdown window, the IMM posterior places most of its mass on the two slowest modes (q = 10−8 and q = 10−6 ), consistent with relatively stable pre- and post-lockdown behaviour at the daily frequency. Around the lockdown onset, the posterior shifts toward the intermediate mode q = 10−4 : its mass rises to 53% on 2020-03-17 and peaks at 57% on 2020-03-18. Thus IMM does react to the change in the innovation sequence, even though it receives no lockdown-aware covariate. Detection does not fully translate into aggregate performance. Despite this visible mode shift, Base+IMM trails MELO by 52.9 MW overall, a 7.5% reduction in RMSE relative to Base+IMM. The lockdown period is more nuanced: IMM-only is essentially tied with MELO in point estimate, 36
(a) full test period
mode probability μi
1.0
1.0
0.8
(b) lockdown transition
0.8 q = 100
0.6
0.6
0.4
0.4
q = 10−2 q = 10−4 q = 10−6 q = 10−8
0.2
0.2
0.0
lockdown start
0.0 2019
2020
2021
Mar 16
Mar 23
Mar 30
Figure 6: IMM-KF mode probabilities µi,t on a five-mode grid q ∈ {10−8 , . . . , 100 }. Colours encode modes by Q-scale, from slow to fast. (a) Full test period. Dashed vertical lines mark the COVID-19 lockdown window. (b) Zoom on the lockdown transition. With the tuned πstay = 0.5, IMM is relatively willing to change mode under the selected validation protocol. Around the lockdown onset, posterior mass shifts away from the two slowest modes and toward the intermediate mode q = 10−4 : from µ(2020-03-14) = [0.39, 0.39, 0.20, 0.02, 0.00] to µ(2020-03-17) = [0.19, 0.20, 0.53, 0.07, 0.01], ordered from q = 10−8 to q = 100 . The q = 10−4 mode absorbs about 53% of the posterior mass on 2020-03-17 and peaks at 57% on 2020-03-18.
whereas Base+IMM performs worse in that short window, with wide bootstrap uncertainty. The underperformance of Base+IMM is therefore not simply due to an absence of mode response at the break. Rather, the way IMM compresses its adaptive behaviour into a single Bayesian mixture prediction appears less favourable for the outer MLpol aggregator on this benchmark. One interpretation is that IMM performs its time-scale selection internally: it mixes states, updates mode probabilities, and then exposes only the resulting Bayesian mixture forecast to the outer aggregator. Once posterior mass concentrates around an intermediate operating point, the output behaves much like a single adaptive filter at that effective scale. By contrast, the EWLS construction exposes all K = 16 fixed-scale predictions separately to MLpol, so the outer aggregator can place weight directly on short-memory corrections during the break while downweighting both slow EWLS experts and raw base forecasts that are no longer useful. Internal mixing may reduce the diversity visible to MLpol. This distinction is apparent in the filter-only ablations. IMM-only is competitive during lockdown (RMSE 1086.5), but adding direct base experts worsens lockdown RMSE to 1276.7. In MELO, adding the base experts to the EWLS pool also slightly hurts lockdown RMSE, from 1073.6 to 1086.1, but the deterioration is much smaller (+12.5 MW), while the improvement before and after lockdown is enough to give the best overall RMSE. This suggests that IMM’s internal Markov mixing leaves fewer distinct adaptive signals for the outer MLpol layer to combine. The issue is not that IMM fails to react to the break; rather, its reaction is largely compressed into a single state-space mixture prediction before the regret aggregator observes it. D.3
Adaptive-Q pool diagnostics
A second diagnostic concerns the pools of VIKING and VBAKF filters with different responsiveness settings. Figure 7 reports the average MLpol weight allocated to each member of these pools, broken down by regime. One might expect the most responsive configurations to receive more weight under the lockdown break. Empirically, however, the outer MLpol layer does not favour the most responsive setting. In the VIKING pool, the most conservative configuration ρb = 10−4 rises from 17.2% pre-lockdown to 40.3% during lockdown, while the most responsive ρb = 10−2 receives only 9.0% in lockdown. The intermediate setting ρb = 10−3 also gains weight, reaching 27.0%. Thus the lockdown allocation within VIKING is concentrated toward the more conservative end of the pool. 37
VIKING pool (varying ρb) 0.40
mean MLpol weight
0.35
VBAKF pool (varying ρ)
pre-lockdown
← most responsive
lockdown
← most responsive
post-lockdown exclude-warm
0.30 0.25 0.20 0.15 0.10 0.05 0.00 ρb = 10−2
ρb = 10−3
ρ = 0.95
ρb = 10−4
ρ = 0.98
ρ = 0.995
Figure 7: Mean MLpol weight allocated to each member of the VIKING pool (left) and the VBAKF pool (right), broken down by regime. For VIKING, smaller ρb is more conservative; for VBAKF, larger ρ is more conservative. Under lockdown, both pools shift weight away from the most responsive configuration toward more conservative ones, although the effect is sharper for VIKING than for VBAKF. The two pools use different internal adaptation mechanisms—an isotropic log-variance random walk for VIKING and a variational covariance/noise update for VBAKF—yet both exhibit this same qualitative reallocation. The VBAKF pool shows the same tendency, though less sharply. Before lockdown, the most responsive configuration ρ = 0.95 receives the largest weight among the three VBAKF members (13.2%, versus 7.2% for ρ = 0.98 and 4.6% for ρ = 0.995). During lockdown, however, the weight shifts toward the slower settings: ρ = 0.98 and ρ = 0.995 rise to 33.4% and 37.5%, respectively, while the more responsive ρ = 0.95 remains non-negligible at 24.8% but is no longer dominant. At the same time, the direct-expert share collapses from 74.9% before lockdown to only 4.3% during lockdown, so the outer aggregator is indeed relying almost entirely on the VBAKF pool in that regime—just not primarily on its most responsive member. This pattern should not be read as a failure of MLpol. Rather, it suggests that the more aggressive adaptive-Q configurations did not consistently deliver lower realised squared loss on this trajectory. A plausible explanation is a bias–variance trade-off: aggressive process-noise adaptation can reduce bias after a break, but it can also increase prediction variability. Since MLpol is driven by realised squared loss, rt,j = 2(ŷt − yt )(ŷt − ỹt,j ), it downweights experts whose additional responsiveness does not translate into lower loss. The fact that this qualitative pattern appears for both VIKING and VBAKF suggests that, on this benchmark, filter-internal adaptive-Q objectives need not align with the outer squared-loss objective.
The contrast with the EWLS grid is that each EWLS expert has a fixed update scale. The outer aggregator therefore compares predictions generated by stable, pre-specified adaptation rules, rather than predictions whose internal process-noise scale is itself being re-estimated online. In our main method, the lockdown period assigns 87.9% of the total MLpol weight to the EWLS pool, with 43.9% on fast experts (γ < 0.99), 30.1% on intermediate experts (γ ∈ [0.99, 0.999)), and only 13.9% on the slowest experts (γ ≥ 0.999). Thus, in our implementation, the EWLS grid leaves the candidate time scales explicit and separate, allowing MLpol to place substantial weight directly on short-memory corrections during the break. D.4
VBAKF diagnostics
VBAKF shows an additional behaviour that helps interpret its weaker overall performance in this experiment. Figure 8 plots the estimated state-uncertainty scale tr(E[P ])/d and the estimated observation-noise scale E[R] over the test period. A natural response to the lockdown break would be a robust increase in state-evolution uncertainty across responsiveness settings. The fitted VBAKF variants do not show such a common pattern. The most responsive setting ρ = 0.95 rises sharply 38
(a) VBAKF process-noise estimate
(b) VBAKF observation-noise estimate
ρ = 0.95
10−2
ρ = 0.98
[R] (MW²)
tr([P]) / d
ρ = 0.995
9 × 10−3
105 8 × 10−3
2019
2020
2021
2019
2020
2021
Figure 8: VBAKF’s estimated state-uncertainty scale tr(E[P ])/d (left) and observation-noise scale E[R] (right) over the test period, for three values of the IW-concentration forgetting factor ρ. The lockdown window is shaded. The most responsive setting ρ = 0.95 shows a pronounced rise in state uncertainty around the transition, whereas the two slower settings remain nearly flat or slightly lower relative to their pre-lockdown levels. The observation-noise estimates also evolve in a strongly configuration-dependent way, especially for smaller ρ. Thus VBAKF does not expose a common, configuration-consistent signal of increased state-evolution uncertainty at the break. around the transition and remains volatile during and after lockdown, but the two slower settings ρ = 0.98 and ρ = 0.995 remain essentially flat or slightly lower relative to the pre-lockdown regime. Quantitatively, the lockdown-to-pre-lockdown ratio of tr(E[P ])/d is 1.161 for ρ = 0.95, but only 0.987 for ρ = 0.98 and 0.988 for ρ = 0.995. The corresponding post-lockdown-to-pre-lockdown ratios are 1.124, 0.979, and 0.982. Thus the VBAKF state-uncertainty posterior provides a response for the most responsive configuration, but not a consistent cross-configuration signal of increased state volatility around the break. At the same time, the estimated observation-noise scale is far from constant: the coefficient of variation of E[R] is 1.188, 0.670, and 0.266 for ρ = 0.95, 0.98, 0.995, respectively. This suggests that the adaptive update partly reallocates uncertainty between state and observation noise in a configuration-dependent way, which need not align with the squared-loss objective used by the outer MLpol layer. One possible mechanism is that the variational covariance update can become path-dependent. In stable periods, the state updates ∆θ are small and the posterior covariance after the Kalman update is also small. The sufficient statistic added to the inverse-Wishart update, AIW = Pnew + (∆θ)(∆θ)⊤ . is therefore small, which can reinforce a tight state posterior. When a break arrives, the filter may respond through a mixture of state-uncertainty and observation-noise adaptation, rather than by exposing a clean set of alternative time-scale predictions to the outer aggregator. This behaviour is consistent with the empirical results: VBAKF-only is competitive during the short lockdown window in point estimate, but the Base+VBAKF pool is significantly worse than MELO overall, before lockdown, and after lockdown. D.5
Synthesis
The adaptive-filter comparison gives a consistent empirical message on this benchmark. We consider three representative adaptive-KF alternatives: isotropic variational adaptation (VIKING), inverseWishart covariance/noise adaptation (VBAKF), and Bayesian multi-model mixing (IMM). Under the same walk-forward tuning protocol, all three trail γ-grid + MLpol on overall RMSE. The paired bootstrap confirms that the overall differences are statistically significant for all three adaptive-KF baselines when their outputs are combined with the direct base experts by the same outer MLpol layer (Table 9). The diagnostics suggest three complementary interpretations. First, IMM reacts to the lockdown transition in its mode probabilities, but its internal Bayesian mixing exposes only a single combined 39
prediction to the outer aggregator. This may reduce the diversity that MLpol can exploit. Second, in the VIKING and VBAKF pools, the outer aggregator does not simply put all lockdown mass on the most responsive adaptive-Q configuration. VIKING shifts strongly toward the more conservative end of its pool, while VBAKF shifts toward slower settings and away from the direct base experts. This suggests that greater internal process-noise responsiveness need not translate into lower realised squared loss. Third, the tuned VBAKF variants do not expose a common, configuration-consistent signal of increased state-evolution uncertainty at the break: the most responsive setting reacts strongly, whereas the slower settings remain nearly flat or slightly lower relative to their pre-lockdown levels. Part of the adaptation is instead absorbed by the estimated observation-noise scale, in a configuration-dependent way. Together, these diagnostics support the design choice made in MELO: keep the adaptation scales fixed and explicit at the expert layer, and use the outer MLpol layer to select among the resulting predictions according to realised squared loss. The comparison does not imply that adaptive Kalman filters are unsuitable in general, nor that no alternative state-space design could close the gap. It shows that, on this deployment-shift benchmark and under the same walk-forward tuning protocol, the adaptive-Q mechanisms considered here do not match the overall performance of the explicit multi-scale EWLS + regret-aggregation separation.
E
Cross-dataset behavior on TabReD industrial datasets
The empirical analysis in Section 5 and Appendix D is centered on RTE-FR Load, a French national electricity-load benchmark spanning the COVID-19 lockdown. This appendix evaluates the framework on two additional industrial regression datasets from the TabReD benchmark [Rubachev et al., 2025]1 : Weather (predicting station temperature, Ttest = 40,840 instances over July 2023) and Delivery-ETA (predicting food-delivery arrival time, Ttest = 36,927 instances). The goal of this appendix is not to claim uniformly large gains across all chronological tabular tasks. Rather, it clarifies the boundary of the framework: MELO is most useful when two conditions hold simultaneously: (1) there is exploitable non-stationarity, and (2) the base pool contains residual diversity that a covariance-aware combination layer can use. The TabReD datasets provide informative negative-control settings in which the base-model residuals are highly correlated on early held-out test segments; consistently with this diagnostic, the EWLS layer remains mildly beneficial but produces much smaller gains than on RTE-FR Load.
E.1
Setup
Base pool. For each TabReD dataset we use four base models: LightGBM, XGBoost, MLP, and MLP-PLR (MLP with periodic-linear numerical embeddings; Gorishniy et al., 2022). Hyperparameters are taken directly from the TabReD authors’ tuned configurations committed in their public repository, without re-tuning. Each model is trained at three random seeds for Weather and one seed for Delivery-ETA; predictions are averaged across seeds.
Online evaluation. We use TabReD’s default time-based split. Base models are trained on the train split with early stopping on the val split; the test split is then traversed once in chronological order, with EWLS and MLpol updated after each observation. Forgetting factors and δ0 are inherited from the RTE-FR Load main experiment without modification (K = 16, γ ∈ [0.95, 0.9998], equivalently h ∈ [20, 5000], δ0 = 10−3 ). The covariance-inflation scale ε0 is selected by walk-forward validation on the val split using the protocol of Appendix B.5.
Patch to the TabReD pipeline. The TabReD code base saves test-set predictions to disk by default, but does not distinguish them by random seed in the prediction filename; we extracted the per-seed predictions from each evaluation directory’s predictions.npz. No model code or hyperparameters were modified. 40
Table 10: Cross-dataset summary. RMSE gains are measured on the full test period. “Mean offdiag corr” is the average off-diagonal entry of the residual correlation matrix among base models, computed only on an initial diagnostic segment of the chronological test stream: the first 300 test days for RTE-FR Load and the first 20,000 test instances for the two TabReD datasets. These correlations are reported only as a post-hoc diversity diagnostic; they are not used for model selection, hyperparameter tuning, expert-pool construction, or online updating. RTE-FR Load uses the 7-base pool of Section 5; TabReD entries use the 4-base pool described in Section E.1. Dataset RTE-FR Load TabReD Weather TabReD Delivery-ETA
E.2
Mean off-diag corr
MELO gain
Ttest
0.43 0.97 0.98
+34.7% +0.30% +0.24%
746 40,840 36,927
Cross-dataset results
Table 10 reports the headline finding. MELO improves over MLpol on Base alone in all three datasets we evaluate, but the magnitude of the improvement varies by two orders of magnitude. This pattern is not a contradiction of the main RTE-FR Load result; it identifies the regime in which the correction layer is most useful. In our experiments, large gains appear to require both a non-stationary test trajectory and a base pool whose residuals contain multiple directions of error. When diagnosticsegment residuals are nearly collinear, as in the two TabReD datasets, the affine combination space has little additional signal to exploit, and the observed gain is correspondingly small. Per-method test-set RMSE for the two TabReD datasets is given in Table 11. Both datasets show the same qualitative pattern as RTE-FR Load: MLpol on Base only improves marginally over a uniform simple average, and Base+EWLS adds a further small improvement. The relative size of these improvements is two orders of magnitude smaller than on RTE-FR Load, consistent with the high diagnostic-segment residual correlations reported in Table 10. Table 11: Per-method RMSE on the two TabReD test sets. Lower is better. Bold: best per column. Method
E.3
Weather (RMSE)
Delivery-ETA (RMSE)
Simple average over base pool MLpol on Base only MELO
1.4668 1.4577 1.4533
0.5449 0.5447 0.5434
Best base predictor (single) ∆ MELO vs Base only
1.4586 (LightGBM) +0.30%
0.5467 (XGBoost) +0.24%
Residual correlation as a diversity diagnostic
The two-orders-of-magnitude gap between the RTE-FR Load and TabReD gains is consistent with the amount of base-pool diversity visible in an initial segment of each chronological test stream. Figure 9 measures this diversity through pairwise residual correlations of the base models. This diagnostic is computed after the online protocols, hyperparameters, and expert pools are fixed. It is not used to select hyperparameters, construct the expert pool, or update the online aggregation pipeline for the reported comparisons; it is reported only to explain whether affine combination had a rich enough signal space to exploit. In a prospective deployment, the analogous diagnostic would be computed on a held-out pre-deployment validation period or on an initial calibration stream not used for model selection. RTE-FR Load’s 7-base pool has substantially more residual diversity than the TabReD pools. In the early-window diagnostic, Lag-1 has low correlation with most feature-based models, and the feature-based models themselves separate into weaker and stronger correlation blocks. This structural asymmetry is important: Lag-1 uses the previous-day target as its information source, whereas the other models map the contemporaneous feature vector to load. The EWLS layer can exploit such 1 TabReD comprises eight tabular datasets curated for studying temporal robustness with realistic industrial feature pipelines
and time-based train/val/test splits. See https://github.com/yandex-research/tabred for code and data.
41
Residual correlation across base models, by dataset RTE-FR Load (first 300/746) mean off-diag = 0.43
1.0
1.00
0.62
0.56
0.54
0.50
XGB
0.11
0.36
0.62
1.00
0.84
0.56
0.57
LGB
0.11
0.34
0.56
0.84
1.00
0.57
0.57
ResNet
0.23
0.37
0.54
0.56
0.57
1.00
0.56
FT-Tr
0.14
0.34
0.50
0.57
0.57
0.56
1.00
0.95
0.97
LGB
1.00
1.00
0.98
0.98
XGB
1.00
1.00
0.96
0.97
XGB
1.00
1.00
0.98
0.98
MLP
0.95
0.96
1.00
0.97
MLP
0.98
0.98
1.00
0.98
MLP-PLR
0.97
0.97
0.97
1.00
MLP-PLR
0.98
0.98
0.98
1.00
FT -T r
B
et N Re s
B
LG
XG
e
AM G
1
dg
g-
Ri
La
0.8
0.6
0.4
Residual correlation
0.54
1.00
0.2
LR
0.17
1.00
LP
GAM
LGB
LP -P
0.34
M
0.37
B
0.34
M
0.36
XG
0.54
B
1.00
LG
0.39
LR
Ridge
LP
0.14
M
0.23
LP -P
0.11
TabReD Delivery-ETA (first 20000/36927) mean off-diag = 0.98
M
0.11
B
0.17
B
0.39
XG
1.00
LG
Lag-1
TabReD Weather (first 20000/40840) mean off-diag = 0.97
0.0
Figure 9: Early-window residual correlation matrices across the three datasets, sorted left-toright by base-pool homogeneity. Correlations are computed on an initial diagnostic segment of each chronological test stream: the first 300 test days for RTE-FR Load and the first 20,000 test instances for the two TabReD datasets. Each cell reports the Pearson correlation between two base models’ residuals (prediction minus truth), and panel titles report mean off-diagonal correlation. These correlations are post-hoc diagnostics only and are not used by MELO or by any compared online method. RTE-FR Load exhibits visibly lower and more structured residual correlation, whereas TabReD Weather and Delivery-ETA remain almost uniformly correlated across base models. (γ)
inductive-bias diversity because the time-varying affine combination z̃t⊤ wt base forecasts make meaningfully different errors.
is useful only when the
This interpretation matches the role of the comparator in Theorem 2. The theorem bounds the loss of the aggregate against a time-varying affine comparator path, but the comparator class is only empirically useful when the base-prediction vector spans several error directions. The residual correlation matrix is therefore a simple empirical proxy for the spanning quality of the base pool. By contrast, the two TabReD datasets show near-uniformly high residual correlation. The four base models consist of two GBDT families and two MLP-based families. They appear to converge to nearly the same mapping from engineered features to target, leaving residuals dominated by shared per-instance errors. In this setting, the effective affine combination space is narrow, so the EWLS layer has little room to improve over base-only aggregation. These results clarify a boundary of the method rather than undermining the main RTE-FR Load result. Large gains require both deployment-time non-stationarity and a base pool with non-redundant error directions. Ruling out tuning convergence as the cause. A natural concern is that the high TabReD correlations are an artifact of all four models being hyperparameter-tuned toward the same loss-landscape optimum, and that under-tuned models would expose hidden diversity. We test this by training a default-hyperparameter XGBoost model with 100 trees, depth 6, and learning rate 0.3 on each TabReD dataset, then measuring its validation residual correlation with the four tuned models. Results are summarised in Table 12. The default-XGBoost probe argues against the tuning-convergence hypothesis. An explicitly undertuned model whose test RMSE is 7% worse on Weather and 1.2% worse on Delivery-ETA still produces validation residuals close to those of the tuned models. The high TabReD correlations therefore appear to reflect a structural property of the datasets and feature pipelines: the heavily engineered features leave little room for inductive-bias diversity, regardless of how each model is configured. E.4
A diagnostic for expected aggregation gain
The empirical regularity in Figure 9 and Table 10 suggests a practical diagnostic, but the figures reported above should be read as post-hoc diagnostics because they use initial segments of the test streams. The online evaluation itself is unchanged: the full test stream is traversed once 42
Table 12: Default-XGBoost probe on the TabReD validation splits. Despite materially worse single-model test RMSE, the default model’s validation residuals remain highly correlated with the four tuned models, and adding it to the pool lowers mean off-diagonal validation correlation by less than 0.01 in both datasets. Weather
Delivery-ETA
1.5643 1.4614 (LightGBM)
0.5534 0.5467 (XGBoost)
Default-XGB vs LightGBM corr Default-XGB vs XGBoost corr Default-XGB vs MLP corr Default-XGB vs MLP-PLR corr
0.944 0.944 0.905 0.912
0.986 0.986 0.973 0.969
Mean off-diag, 4 tuned bases Mean off-diag, +default-XGB (5) ∆
0.946 0.938 −0.008
0.983 0.981 −0.002
Default-XGB test RMSE Best tuned-model test RMSE
chronologically, and no residual-correlation statistic enters the online updates. In a prospective deployment, the analogous quantity would be computed on a held-out pre-deployment validation period or on an initial calibration stream not used for model selection. Low residual correlation indicates that the pool contains multiple error directions, making affine combination more promising. High residual correlation indicates that the base models make nearly the same errors, so MELO has limited room to improve over base-only aggregation. These results suggest a qualitative diagnostic. Mean off-diagonal correlation above 0.95 suggests that large gains from affine combination are unlikely unless the deployment period introduces a new source of residual diversity. In this regime, the framework may still be useful as a conservative adaptation layer, but the expected improvement over base-only aggregation should be modest. Mean off-diagonal correlation below 0.7 suggests that the base pool contains several distinct error directions, so larger gains are plausible, particularly when the deployment period also involves non-stationarity. This diagnostic complements the path-length analysis in Theorem 1. Empirical gains depend on both the variability of the useful comparator path and the spanning quality of the base-prediction vector. Either condition can be the binding constraint. RTE-FR Load combines residual diversity with a severe deployment-time shift; the two TabReD datasets have highly correlated residuals, and the corresponding gains are small. Limitations of the diagnostic. The diagnostic is directional, not a calibrated estimator. The gain magnitude depends on the joint distribution of base errors, the realised deployment path, and the losses used by the online aggregation layer, not only on marginal pairwise correlations. Three datasets span a useful range but do not identify a precise functional relationship. A more refined diagnostic would use the spectrum of the residual covariance matrix, i.e., the effective number of independent error directions; on RTE-FR Load this spectrum has approximately four non-trivial directions, whereas on TabReD it concentrates near one. E.5
Discussion
The cross-dataset comparison is consistent with the intended deployment profile of the framework. In the RTE-FR Load stress test, where residual drift is exploitable and the base pool contains diverse error directions, MELO delivers large gains through the COVID-19 lockdown. On the two TabReD datasets considered here, where residuals are highly correlated, the same mechanism has limited room to improve over base-only aggregation. This should be read as a boundary condition rather than a failure mode. The oracle inequality gives a conservative comparison with the original base experts and the added EWLS experts up to the aggregation overhead; it does not imply finite-sample monotonic improvement whenever experts are added. Empirically, Base+EWLS is positive relative to Base-only on all three datasets considered here, but the TabReD results show that the magnitude of the benefit is governed by detectable residual 43
diversity. The diagnostic of Section E.4 therefore helps set realistic expectations before deployment when computed on a genuine pre-deployment validation or calibration period.
F
Robustness to base-pool diversity: a heterogeneity sweep
The main experiment (Section 5) uses F = 7 base predictors spanning linear, tree, additive, and deep families—an intentionally heterogeneous pool. A natural question is how the aggregation gain reported there depends on this heterogeneity. We answer it here by sweeping the number of distinct families in the pool and enumerating, at each level, every possible subset. 7 Setup. For each heterogeneity level H ∈ {1, 2, . . . , 7}, we enumerate all H size-H subsets of the seven families (127 subsets in total). For each subset we rerun the complete aggregation pipeline: MLpol on the base subset only, MLpol on the EWLS pool constructed from that subset, and MLpol on the combined Base+EWLS pool. All other hyperparameters (γ grid, ε0 , α) and the test period (2019-01-01 to 2021-01-15) are identical to the main experiment. For each (H, method) we report 7 the mean and the 2.5%–97.5% percentile band of overall RMSE across the H subsets at that level, so the shown spread captures the sensitivity of the result to which families happen to be in the pool. RMSE vs. heterogeneity. Figure 10 reports the three MLpol variants (Base only, EWLS only, Base+EWLS) and the best standalone single-γ EWLS (hindsight reference) as a function of H. The headline pattern is sharp: • MLpol on Base only needs pool heterogeneity to work. At H = 1 (one family, one predictor) it is trivially the best single base and sits at mean RMSE ≈ 1900 MW; adding families drives it down monotonically, reaching 1004 MW at H = 7 (the main experiment). • Both EWLS-augmented methods are already within 110 MW of their H = 7 value by H = 2, and within 35 MW by H = 3. The marginal return to further heterogeneity is small. • Already at H=2, a Base+EWLS pool outperforms the full seven-family Base-only pool. • The p2.5–p97.5 spread shrinks steeply with H: at H = 1, the choice of family determines whether overall RMSE is 750 or 3200 MW; by H = 4 every subset sits within 100 MW of the mean. Thus small pools are sensitive to which model family happens to be chosen, whereas larger pools are much more stable. Since increasing H only adds base forecasts to the MLpol pool and has negligible computational cost relative to fitting the base models, this variability favours using as broad a base pool as is available. The two sources of diversity are substitutable, not additive. Figure 11 decomposes the aggregation gain into three pairwise contrasts. Two observations: (i) EWLS substitutes for family heterogeneity, monotonically. The gain Base+EWLS − Base-only is 568 MW at H = 1 and falls monotonically to 325 MW at H = 7 (the main-experiment value). The EWLS-only vs. Base-only gain traces the same curve, offset by a near-constant small amount. This is the core substitutability signal: as family heterogeneity accumulates in the base pool, the marginal value of adding the EWLS layer declines—precisely because the base pool itself is now contributing the diversity EWLS would otherwise provide. (ii) Base experts contribute an orthogonal signal of roughly constant size. The gain Base+EWLS − EWLS-only is ≈ 12 MW for all H ≥ 2, with p2.5–p97.5 bands tightly concentrated around the mean. The base-expert channel contributes a small but persistent signal that is not captured by any convex combination of single-memory EWLS experts, and whose magnitude is essentially independent of how many families are in the pool. Together (i) and (ii) characterise a clear trade-off: most of the value of the two-layer framework at low H comes from EWLS (it can rescue a homogeneous pool), and a roughly constant residual value at every H comes from the complementary base channel (orthogonal to memory-length aggregation). Family heterogeneity does not replace adaptation under drift. Figure 12 splits the results by regime. The pre-lockdown and post-lockdown panels show the same diminishing-returns pattern as the overall metric: Base-only catches up as H grows, while Base+EWLS is close to asymptote by H = 3. The lockdown panel is qualitatively different—in the structural break, Base+EWLS 44
Overall RMSE vs. number of base-expert families (mean ± p2.5--p97.5 across all (HF ) subsets) MLpol on Base only MLpol on EWLS only MELO (MLpol on Base+EWLS) Best single EWLS (hindsight)
Overall RMSE (MW)
3 × 103
2 × 103
103
1
2
3 4 5 Number of base-expert families H
6
7
Figure 10: Overall RMSE as afunction of the number of base-expert families H (log scale). Each 7 curve is the mean over all H size-H subsets at that level; the shaded band is the 2.5%–97.5% percentile range across subsets. MLpol on Base only (vermilion) needs pool heterogeneity to become competitive; both EWLS-augmented methods (blue, green) plateau within 35 MW of their H = 7 value by H = 3. The grey dash-dot line is the hindsight-best standalone single γ per subset, shown as a reference for what a “best single memory length” would achieve without aggregation. (a) Gain decomposition: mean ± p2.5--p97.5 band
600
400
200
Base+EWLS gain over Base-only (MW)
800 Overall RMSE gain (MW)
(b) Per-subset gain with medians (horizontal bars)
Base+EWLS vs. Base-only EWLS-only vs. Base-only Base+EWLS vs. EWLS-only
0
800
600
400
200
0 1
2
3 4 5 Number of base-expert families H
6
7
1
2
3 4 5 6 Number of base-expert families H
7
Figure 11: (a) Decomposition of the aggregation gain into three pairwise contrasts as a function of H. Base+EWLS vs. Base-only (green) and EWLS-only vs. Base-only (blue) trace a near-identical monotone decline from ∼ 570 MW at H = 1 to ∼ 320 MW at H = 7, the signature of substitutability between EWLS and family heterogeneity. Base+EWLS vs. EWLS-only (orange) is flat at ≈ 12 MW for all H ≥ 2: base experts contribute a small orthogonal signal whose magnitude is independent of pool diversity. (b) Per-subset values of the Base+EWLS gain (green scatter) with per-H medians (horizontal bars). Positive at every subset we evaluate.
RMSE is essentially flat in H (median 1151 MW at H = 2; 1099 MW at H = 7, a 5% reduction), while Base-only plateaus at around 2400 MW even at H = 7. This matches the theoretical picture of Section 4: under heavy drift, static base experts tend to fail together even when diverse, so adding families contributes no tracking capability—only the EWLS layer provides explicit online adaptation. Family heterogeneity is a partial remedy for poor stable-regime predictions; it is not a substitute for multi-time-scale adaptation under drift. Summary. Sweeping the F = 7 base families yields three findings that together characterise how the two-layer framework interacts with base-pool construction: 45
pre-lockdown
lockdown
post-lockdown
MLpol on Base only MLpol on EWLS only MELO (MLpol on Base+EWLS)
3 × 103
2 × 103
3 × 103 RMSE (MW)
exclude-warm
3 × 103
4 × 103
3 × 103
2 × 103
2 × 103 2 × 103 103 103
103
6 × 10
103
6 × 102 1
2
3
4 H
5
6
7
1
2
3
4 H
5
6
7
2
1
2
3
4 H
5
6
7
1
2
3
4 H
5
6
7
Figure 12: Per-regime RMSE as a function of H (log scale). In stable regimes (pre- and postlockdown), Base-only catches up as heterogeneity grows. Under the lockdown structural break, the Base+EWLS curve is nearly flat in H while Base-only plateaus at ∼ 2400 MW: the drift-tracking contribution of the framework comes from EWLS and cannot be replaced by pool heterogeneity. The main-experiment 1. EWLS and family heterogeneity are partial substitutes. Base+EWLS − Base-only gain of 325 MW grows monotonically to 568 MW as the pool narrows to a single family. EWLS provides an alternative source of diversity, expressed along the time-scale axis rather than the model-family axis. 2. A homogeneous EWLS-augmented pool outperforms a heterogeneous static one. A twofamily Base+EWLS pool already beats the seven-family Base-only pool overall. Practitioners without access to a rich model zoo can still obtain most of the benefit reported in the main experiment from a small pool, as long as the EWLS layer is present. 3. Base experts contribute a small, approximately H-independent complementary signal. Adding base experts to an EWLS-only pool improves RMSE by ≈ 12 MW at any pool diversity, so the two layers are complementary rather than redundant—even when diversity is abundant. The headline aggregation gain reported in the main experiment is therefore not driven by the structural richness of the heterogeneous base pool: it is present at every pool composition tested, and is especially large in the low-diversity configurations and drift regimes where static predictors struggle most. Because the identity of the chosen family matters substantially at small H, and because admitting additional available families is cheap at aggregation time, the sweep also supports the practical default of using the broadest causal base pool available.
G
Weight dynamics: per-expert details
Section 5.2 presented the MLpol weight trajectory with base experts aggregated into a single band and EWLS experts shown as coloured strata, which emphasises the across-pool base/EWLS reallocation. Two further phenomena are important for interpreting the regret dynamics in Figure 3 but cannot be read off that view: (i) the within-base reallocation, where Lag-1 absorbs much of the remaining base-expert mass under the break, and (ii) the within-EWLS reshuffling, where the weight mass moves across time scales from slow to fast and back. Fully disaggregated view. Figure 13 shows the individual weight trajectory of all N = M +K = 23 experts, with base experts coloured individually and EWLS experts coloured by γ on the RdYlBu ramp. The three macro-level claims of Section 5.2—base dominance in pre-lockdown, abrupt lockdown reallocation, and medium-memory dominance in the recovery—are all visible, but two finer patterns become apparent as well. First, the within-base mass at lockdown onset is not simply spread uniformly across the original leaders: Lag-1, which had negligible weight before the break, becomes the dominant member of the remaining base allocation. Second, within the EWLS strata the reddest bands (smallest γ) widen abruptly at lockdown onset at the expense of the bluest bands (near-static γ → 1), consistent with the optimal memory length h∗ shrinking under drift (Theorem 1). We unpack each of these two axes next. Within-base reallocation: Lag-1 takes over during the break. Figure 14 isolates the M = 7 base experts. The top panel shows their raw MLpol weight; the bottom panel normalises to the base sub-simplex (each expert’s share of the total base weight), which removes the confounding effect of base-pool shrinkage during the lockdown and isolates internal reallocation. 46
1.0
Base experts
Lockdown
Lag-1 Ridge GAM XGBoost LightGBM ResNet FT-Tr.
MLpol weight
0.8
0.6
0.4 EWLS (γ) γ = 0.95000
0.2
γ = 0.98968 γ = 0.99787 γ = 0.99956 γ = 1.00000
0.0
Apr 2019
Jul 2019
Oct 2019
Jan 2020
Apr 2020
Jul 2020
Oct 2020
Jan 2021
Date
Figure 13: Fully disaggregated MLpol weight trajectory. All N = M +K = 23 experts are shown as individual strata. Base experts (bottom of the stack) are coloured individually; EWLS experts (top) are coloured by forgetting factor γ on the RdYlBu ramp (red = fast forgetting, γ = 0.95; blue = slow / near-static, γ → 1). Dashed red vertical lines mark the COVID-19 lockdown window (2020-03-17 to 2020-05-11). The plot begins on 2019-02-15 for readability, omitting the initial uniform-weight transient. The EWLS legend shows five representative γ’s; all K = 16 values appear in the plot. Three regularities stand out: • In the stable pre-lockdown regime, the base allocation is dominated by the structurally trained models, roughly in the order LightGBM > ResNet > XGBoost > GAM, each receiving 6–13% of the total MLpol weight. Lag-1 and Ridge sit at essentially zero. • During lockdown, the structurally trained base experts collapse jointly: GAM, XGBoost, LightGBM, ResNet, and FT-Transformer fall from a combined 41.3% pre-lockdown weight to 2.4% during lockdown. Lag-1 moves in the opposite direction, rising from 0.5% to 9.7%. Thus the base channel does not disappear entirely; rather, most of the remaining base-expert mass is reassigned to the one-step autoregressive baseline. On the base sub-simplex, Lag-1 accounts for roughly 80% of the lockdown base allocation. This is consistent with the role of Lag-1 as an extremely short-memory raw expert, analogous in the base channel to the fast-forgetting EWLS experts in the correction channel. • Post-lockdown, the structurally trained experts partially recover, but not to their pre-lockdown levels: their combined weight rises to 12.0%, compared with 41.3% before lockdown. GAM and ResNet recover most visibly, while XGBoost and LightGBM remain near zero. The base pool has learned, in effect, that its offline-trained members carry less reliable structure after the regime shift.
H
Extensibility: signal-informed forecasts as additional experts
The main experiments deliberately use frozen pre-lockdown base models and no external regime signal, in order to isolate the effect of the online correction layer. The framework is not tied to this particular base pool. Any causal forecast whose loss can be evaluated online can be added as an expert in the final MLpol layer, and can optionally be included in the EWLS design vector. This gives a simple oracle-inequality interpretation of extensibility. Suppose that a new causal predictor mX is added to an existing pool of N experts. Applying Lemma A.5 to the enlarged pool gives an aggregate ŷ N +1 satisfying T X t=1
ℓt (ŷtN +1 ) ≤
T X
ℓt (mX (t)) + G
p
(N + 1)T ,
t=1
where G depends on the clipping radius, under the same boundedness conditions as in the MLpol oracle inequality. Thus adding a signal-informed forecast gives the aggregate a new comparator, at 47
MLpol weight
Lockdown
Lag-1 Ridge GAM XGBoost LightGBM ResNet FT-Tr.
0.6 0.4 0.2 0.0
Share of base pool
1.0 0.8 0.6 0.4 0.2 0.0
Apr 2019
Jul 2019
Oct 2019
Jan 2020
Apr 2020
Jul 2020
Oct 2020
Jan 2021
Date
Figure 14: Within-base reallocation. Top: raw MLpol weight of each of the M = 7 base experts. Bottom: the same weights normalised to sum to one across the base pool, so one reads each expert’s share of the base allocation independently of the base-pool total. The first six weeks covering the uniform-weight warm-up transient are clipped. Dashed red vertical lines mark the COVID-19 lockdown window. Under the structural break the structurally-trained base experts (GAM/XGBoost/LightGBM/ResNet/FT-Transformer) collapse jointly while Lag-1, which had been receiving essentially zero weight in pre-lockdown, takes over the majority of the base sub-simplex. p √ the cost of increasing the finite-pool overhead from order N T to order (N + 1)T . This does not imply finite-sample monotonic improvement over the previous pool, but it means that useful auxiliary signals can be incorporated without changing the algorithmic structure. Table 13 tests this extensibility by adding a strong daily-retrained TabICL forecast augmented with the Oxford Government Response Index (TabICL+GRI). This is a stronger-information setting than the main experiment, so we report it as a diagnostic rather than as the primary comparison. Three increasingly rich integrations all reduce overall RMSE: adding TabICL+GRI to a base-only MLpol pool reduces RMSE from 1004.03 to 669.14; adding it only as an extra raw expert to the Base+EWLS pool reduces RMSE from 655.78 to 621.06; and additionally including it in the EWLS design vector reaches 601.86, with the largest gain during lockdown. Table 13: Extensibility experiment with a regime-informed foundation-model forecast. TabICL+GRI denotes the daily-retrained TabICL forecast augmented with the Oxford Government Response Index. “Raw only” adds TabICL+GRI only as an expert in the final MLpol pool. “Via EWLS” also includes TabICL+GRI among the inputs to the EWLS correction layer. Lower is better. Method
Pre-lockdown
Lockdown
Post-lockdown
Overall
690.54 667.96 607.81 623.10 594.90 582.98
2452.70 1132.62 1192.50 1086.07 980.20 895.87
907.04 699.46 606.08 579.25 559.19 550.27
1004.03 723.43 669.14 655.78 621.06 601.86
MLpol on Base only TabICL+GRI MLpol on Base + TabICL+GRI MLpol on Base+EWLS MLpol on Base+EWLS + TabICL+GRI raw only MLpol on Base+EWLS + TabICL+GRI via EWLS
This supports the intended deployment interpretation of MELO: it is not a closed ensemble of weak static models, but an online adaptation layer around any available causal forecasts. Stronger base forecasts can be incorporated directly, and the EWLS layer can still exploit residual structure in their joint prediction vector. 48