Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
76,807 characters · 21 sections · 81 citation commands
From Many Models, One: Macroeconomic Forecasting with Reservoir Ensembles
\makeatletter \let\ofootnote\@ifstar\footnote@star\ofootnote \def\footnote@star#1{{\let\thefootnote\relax\footnotetext{#1}}} \@ifstar\footnote@star\ofootnote*{ YCL and LG acknowledge the financial support of the FoKo of the University of St. Gallen (Project Nr. 1022324, “Machine Learning Techniques for Macroeconomic Forecasting and Pricing of Financial Derivatives”). GB thanks the Great Minds Postdoctoral Fellowship program of the University of St. Gallen, which made this research possible. Python code accompanying this manuscript is available at the online repository \url{https://github.com/RCEconModelling/Ensemble-MFESN} } \makeatother
\noindentKeywords: echo state networks, model combination, online learning, time series prediction
Many mathematical, statistical, and econometric tools have been developed to solve the problem of predicting stochastic and deterministic processes. In empirical applications, however, it is often difficult to identify which model, or even which model class, is most appropriate for a given task. Models that perform well on average may fail during particular periods. This situation is common in forecasting settings, for example, involving macroeconomic time series. Structural breaks, financial crises, and policy regime shifts can significantly affect predictive accuracy. Since forecasting performance can vary substantially over time, it is natural, when model updating is infeasible or unsuitable, to revise the choice of model(s) as new information becomes available.
The model combination approach formalizes the use of multiple prediction schemes. The forecaster relies on a finite set, referred to as an ensemble, of distinct models or algorithms, termed experts, to predict an outcome of interest. Experts in an ensemble are evaluated and reweighted over time based on their predictive losses. Ensemble methods can hence adapt to changing time series dynamics by favoring different models across periods. The objective is to construct data-driven weights that aggregate individual expert forecasts to minimize cumulative loss over time. A well-designed combination scheme is expected to select the best-performing models, achieving performance close to the minimum loss in hindsight.
Forecast combination is a long-standing idea in econometrics and time series analysis Bates1969,Winkler1983,clemen1989combining,grangerInvitedReviewCombining1989\@ifstar\footnote@star\ofootnote{Boosting, bagging, and random forest methods, which explicitly combine many learners into one predictor to improve performance, have also become increasingly popular Athey2019,atheyMachineLearningMethods2019.}, with a large literature documenting its empirical effectiveness in a variety of macroeconomic and financial forecasting contexts. Subsequent contributions have studied both linear and nonlinear combination schemes, as well as the role of time variation in forecast weights (see timmermann2006forecast, for a comprehensive survey). More recent work has emphasized adaptive and data-driven approaches to forecast combination in high-dimensional and potentially unstable environments, including settings with structural change and model misspecification (e.g., hansenLeastsquaresForecastAveraging2008,chengForecastingFactoraugmentedRegression2015,kimMiningBigData2018,DenReijer2019,bolhuisDeusExMachina2020,fulekyMacroeconomicForecastingEra2020,masiniMachineLearningAdvances2023), as well as considered Bayesian WANG20231518 and density combination schemes wallisCombiningForecastsForty2011.
The literature on online combination schemes and ensembles is rich. In the standard framework of prediction with expert advice Littlestone1994,Freund1997,vovk1995game,Cesa-Bianchi2006, or its variant known as the Hedge setting, much of the work focuses on developing combination algorithms that enjoy guarantees on regret, defined as the difference between the cumulative loss of the forecaster and that of the best among $K$ experts. Prominent examples include {\it Follow-the-Leader} (FTL), which assigns uniform weights at each round to the model(s) with the smallest cumulative loss so far; it is known to achieve an $O(\log T)$ regret for strongly convex losses in stochastic settings (with $T$ a time horizon), while potentially performing poorly under worst-case data vanervenAdaptiveHedge2011. The {\it Hedge} algorithm Freund1997, also called the exponentially weighted forecaster, provides a smoothed alternative to FTL and enjoys an anytime worst-case $O(\sqrt{T \log K})$ regret bound under an optimal learning rate. In stochastic environments, allowing the learning rate of Hedge to be data-adaptive yields regret comparable to FTL, motivating schemes such as Adaptive Hedge vanervenAdaptiveHedge2011,rooijFollowLeaderIf2014. Significant effort has been devoted to develop strategies that combine the best properties of these approaches (see, for example, koolenLearningLearningRate2014,wintenbergerOptimalLearningBernstein2017,mourtadaOptimalityHedgeAlgorithm2019,itoBestofBestWorlds2024,wintenbergerStochasticOnlineConvex2024).
In this paper, we bridge the literature on online ensemble learning (see, e.g., shalev_online) with econometric forecast combination timmermann2006forecast in a general time series setup. From a theoretical perspective, we study the properties of Follow-the-Leader and {\it Decreasing Hedge} (DecHedge) algorithms in the stochastic setting for i.i.d. (in time) and $\varphi$-mixing losses, building upon the developments in mourtadaOptimalityHedgeAlgorithm2019 and rooijFollowLeaderIf2014, respectively. For both algorithms, we derive Hoeffding- and Bernstein-type regret bounds that characterize their cumulative performance relative to the best expert in hindsight. While the dependence on the number of experts is mild (logarithmic), the gap between the least expected loss and the runner-up plays a key role in our bounds.
On the empirical side, we construct our ensembles using Multi-Frequency Echo State Network (MFESN) models introduced in ballarin2022reservoir, which have been shown to outperform state-of-the-art mixed-frequency forecasting methods such as MIxed DAta Sampling (MIDAS) Ghysels2007 and Dynamic Factor Models (DFMs) geweke1977dynamic,sargent1977business. The MFESN models are based on Echo State Networks (ESNs), a class of recurrent-type neural networks with randomly sampled state parameters jaeger2001,maass2,lukosevicius. This model family has been actively studied in terms of approximation, generalization, and memory properties RC6,RC7,RC8,RC10,RC12,RC23. Using quarterly U.S. GDP growth and a set of $33$ monthly and daily financial and macroeconomic variables in ballarin2022reservoir, we show that MFESN ensembles substantially reduce mean-squared forecasting errors. Additionally, allowing for heterogeneity in random parameter draws and state leak rates across the ensemble yields up to $40\%$ improvements relative to baseline MFESNs. These results set new benchmarks in the prediction of U.S. GDP growth with multi-frequency regressors. Finally, we report the experts' weight evolution in our forecasting exercise, showing that a handful of models receive most of the weight.
The remainder of the paper is organized as follows. Section (ref) introduces the ensemble forecasting framework and formalizes the online learning setting for combination strategies. Section (ref) reviews popular combination strategies and presents theoretical results, deriving regret bounds for Follow-the-Leader and decreasing Hedge under i.i.d. and $\varphi$-mixing losses. Section (ref) describes the Multi-Frequency Echo State Network architecture and the construction of our ensemble variants, EN-MFESN-RP and EN-MFESN-$\alpha$RP. Section (ref) presents an application of these ensembles in the mixed-frequency U.S. GDP growth forecasting setting of ballarin2022reservoir and comparing them to state-of-the-art methods. Section (ref) concludes.
\paragraph{Notation.} For $n \in \mathbb{N}^+$, we define $[n] := \{1, \ldots, n\}$. Given a vector ${v} \in \mathbb{R}^n$ (or $\bm{v} \in \mathbb{R}^n$), we denote its entries by $v_i$, $i \in [n]$. For a vector ${v} \in \mathbb{R}^n$, we write ${v} \geq 0$ whenever ${v}_i\geq 0$ for all $i\in [n]$. The symbol $\bm{1}_n \in \mathbb{R}^n$ stands for the vector of ones. All random variables are defined on a fixed probability space $(\Omega, \mathcal{F}, \mathbb{P})$ and follow context-specific notation: Quantities traditionally denoted by Latin or Greek lowercase letters (for example, $\ell$, $\varepsilon$, $\zeta$) remain lowercase even when random. For all other random variables, we use uppercase letters.
In this section, we introduce the general econometric setting of interest, which involves constructing time series predictions by combining (or selecting) models from a set of experts.
Let $\{Y_t\}_{t \in \mathbb{Z}}$, $Y_t \in \mathbb{R}$, be the target time series, and $\{Z_t\}_{t \in \mathbb{Z}}$, $Z_t \in \mathbb{R}^{d}$, a vector of regressors, which may include lagged values of $Y_t$. We consider the problem of predicting the target series $Y_{t+h}$ at the forecasting horizon $h \in \mathbb{N}^+$. Let $(\Omega,\mathcal F,\mathbb P)$ be a fixed probability space on which all the random variables are defined. For each $t \in[T]$, let $Z_t$ be $\mathcal{F}_{t}$-measurable with $\mathcal{F}_t:=\sigma\left(Z_t\right)$. Under relatively weak conditions of the joint law of $Y_t$ and $Z_t$ the $h$-step-ahead minimum mean square error (MSE) predictor of $Y_{t+h}$ at time $t$, with $h \in[H]$, is the conditional expectation function (CEF) given by
where $\mathcal{G}$ is the class of $\mathcal{F}_t$-measurable functions.
The core challenge arising when solving (ref), even in the one-step-ahead prediction setting ($h=1$), is that minimizing MSE over (possibly infinite-dimensional) $\mathcal G$ class of functions can be infeasible or computationally and statistically demanding. This is typically addressed by structural assumptions on the joint data-generating process $(Y_t,Z_t)$. For example, standard choices include linear models BrocDavisYellowBook, semi- or nonparametric classes of sufficient regularity stoneOptimalGlobalRates1982,Li2009,Tsybakov2009, or a reproducing kernel Hilbert space setting Steinwart2008,scholkopf2002learning,berlinet2011reproducing. In this work, we study a standard online aggregation setting in which multiple models (“experts”) are combined, each potentially excelling under different finite sample realizations, but always assuming that a best-performing model exists (see Section (ref)). The focus of the paper is on studying schemes to weight experts in order to produce a combined forecast that over time aligns with the optimal expert.
The concept of prediction with experts goes back to the works of Littlestone1994,foster1991prediction,foster1993randomization,vovk1995game,cesa1997use. A forecaster is interested in obtaining, at each time step $t \in [T]$, a prediction $\widehat{Y}_{t+h}$, for a horizon $h\in [H]$, based on $K \geq 1$ experts indexed by $k \in [K]$, each providing their own forecast $\widehat{Y}_{t+h}^{(k)}$. The forecaster constructs a convex combination of $K$ expert predictions using a weight vector $\bm{\omega}_{t,h} := (\omega^{(1)}_{t,h}, \ldots, \omega^{(K)}_{t,h})^\top \in \mathbb{R}^K$ such that
and producing a forecast
We will also use $\Delta^{K-1}:=\{\bm{\omega} \in \mathbb{R}_{+}^K: \bm{\omega}^\top \bm{1}_K=1\}$ to indicate the real simplex.
The accuracy of each forecast is assessed using a loss function $\ell: \mathbb{R} \times \mathbb{R} \rightarrow \mathbb{R}_{+}$. We define, for each expert $k \in [K]$ and each $t \in [T]$, the associated loss at $h$-horizon forecasting task of predicting ${Y}_{t+h}$, $h \in [H]$, as
the vector containing the losses of all $K$ experts as $\bm{\ell}_{t,h} := (\ell^{(1)}_{t,h}, \ldots, \ell^{(K)}_{t,h})^\top$ and the forecaster's loss $\overline{\ell}_{t,h} := \bm{\omega}_{t,h}^\top \bm{\ell}_{t,h}$. We also define the cumulative loss of expert $k$ and of the forecaster, incurred by time $t$ when forecasting $h$ steps ahead, as
respectively, with the convention that $L^{(k)}_{t,h}=0$ and $\overline{L}_{t,h}=0$ for all $t\leq h$, and construct a vector $\bm{L}_{t,h} := (L^{(1)}_{t,h}, \ldots, L^{(K)}_{t,h})^\top$ of cumulative losses of all $K$ experts.
The natural objective of online learning procedures is to minimize the loss $\overline{L}_{t,h}$ experienced by the forecaster when combining $K$ experts' predictions. Assuming that $[K]$ does not change over time, one can introduce cumulative regret at $t \in [T]$ for horizon $h\in [H]$ as
which, for a given forecasting horizon $h$, is the difference between the forecaster's cumulative loss and that of the best expert in hindsight. Equivalently, set $\ell^{*}_{t,h} := \min_{k \in [K]} \ell^{(k)}_{t,h}$, so that the instantaneous prediction regret is $r_{t,h} := \overline{\ell}_{t,h} - \ell^{*}_{t,h}$, while the $k$th model regret is $r^{(k)}_{t,h} := \ell^{(k)}_{t,h} - \ell^{*}_{t,h}$. Notice that the optimal expert at time $t$, indexed by $k^\star_{t,h} \in \arg\min_{k \in [K]} \ell^{(k)}_{t,h}$, by definition yields $r^{(k^\star_{t,h})}_{t,h} = 0$. Therefore, the goal of an online combination scheme consists in minimizing cumulative regret $\overline{R}_{t,h}$ at time $t\in [T]$ for horizon $h\in [H]$, or, equivalently, reliably selecting, over the prediction window, the best possible candidate among the available $K$ experts. One may also consider an adversarial setting, where losses are defined via $\ell(y, \hat{y})$ and experts' $\widehat{Y}_{t+h}^{(k)}$ are fixed before the environment acts. Equivalently, the environment chooses $Y_{t+h}$, to induce $ \boldsymbol{\ell}_{t, h} $, such that the regret is maximized, leading to a worst-case scenario.
Throughout the paper, unless stated otherwise, we impose the following assumptions on expert and forecaster losses for all $h\in [H]$.
The definition of $\varphi$-mixing dependence in (ii) and (iii) is in the sense of Definition (ref) in the Appendix, following rioAsymptoticTheoryWeakly2017. The dependence conditions are imposed on losses (hence, data and experts) over time, but no independence is assumed across the $K$ losses.
Assumption (ref) is standard in online learning and, up to a fixed rescaling, requires losses to be bounded. Assumption (ref) imposes conditions on the temporal dependence of loss vectors, as the prediction algorithm evolves. While much of the literature focuses on the i.i.d. case, mixing conditions such as Assumption (ref)(ii) are empirically more relevant in time-series applications. Notice also that for any $k\in [K]$, since $0 \leq \varphi_n^{(k)} \leq 1$, Assumption (ref)(i) implies Assumption (ref)(i) and Assumption (ref)(ii) implies Assumption (ref)(ii). Moreover, in each of these assumptions, part (ii) implies part (i). All assumptions are imposed on the loss processes $\{\ell^{(k)}_{t,h}\}_t$, rather than on the joint data process $(Y_t,Z_t)$. If experts are fixed and $(Y_t,Z_t)$ is mixing, the induced loss processes inherit this property under measurable maps. Even under mild nonstationarity (e.g., local stationarity), these losses often remain weakly dependent, making uniform mixing a plausible assumption in practice.
In the remainder of the paper, we focus on a single forecasting horizon $h\in \mathbb{N}$ and hence drop the $h$ subscript where convenient. Our results extend directly to the multi-horizon case by analyzing cumulative regret separately for each horizon $h \in[H]$.
We now turn to an in-depth discussion of the Follow-the-Leader and Hedge online learning methods. Table (ref) provides an overview of these schemes, as well as classical combination approaches based on averaging and rolling mean-squared-error weighting (see Appendix (ref) for more details). In Section (ref), we theoretically study the stochastic setting with and without temporal dependence, and derive regret guarantees for both FTL and decreasing Hedge.
An intuitive first step to construct expert weights is the so-called {\it Follow-the-Leader} (FTL) strategy: At each round, uniform weights are assigned to the (subset of) model(s) with the smallest cumulative loss so far. Formally, let $\textnormal{\otextsc{ftl}}_t \in \arg\min_{k \in [K]} L^{(k)}_{t-1}$ be the subset of leader-experts at time $t$, that is, the experts with the least cumulative loss up to $t-1$. The Follow-the-Leader weights are defined as
with the convention that $\omega^{(k)}_{\textnormal{\otextsc{ftl}},1} = 1/K$. FTL is known to perform well in several settings, particularly in the case of i.i.d. losses. However, it has the downside of being sensitive to adversarially-generated losses, as it overreacts to short-term fluctuations. rooijFollowLeaderIf2014 derive an upper bound on the cumulative regret incurred by a forecaster using the FTL strategy due to the variability of the losses and the number of leader changes.
Let $s_t := \max_{k \in [K]} \ell^{(k)}_{t} - \min_{k \in [K]} \ell^{(k)}_{t}$ and $S_t := \max\{s_1, \ldots, s_t\}$ denote the instantaneous maximal loss differential and maximum loss differential over $t$ periods, respectively. Define $c_t=1$ if there exists an expert $k'$ with $k' \in \mathrm{FTL}_{t-1}$ but $k' \notin \mathrm{FTL}_{t}$ (i.e., the leader set changes at $t$), and $c_t=0$ otherwise. We can write $c_t=\mathds{1}\{ \textnormal{\otextsc{ftl}}_{t-1} \not= \textnormal{\otextsc{ftl}}_{t} \}$. Let then $C_t := \sum_{\tau=1}^t c_\tau$ be the total number of times the leader set changes up to time $t\in [T]$.
Lemma (ref) shows that FTL is effective whenever leader changes are infrequent, and the loss range is small. Under i.i.d. losses whose expectations are well separated, frequent leader changes are unlikely rooijFollowLeaderIf2014. In Section (ref), we provide a sharper result in a stochastic setting, while also allowing for dependence.
Even though FTL relies on identifying a single best expert for prediction, the underlying setup can be generalized to non-sparse weighting vectors. The overarching framework is that of the {\it Hedge} family of algorithms Littlestone1994,Freund1997,Freund1999,cesa1997use, also known as exponentially weighted average forecasters. In Hedge, experts are weighted according to
where $\{\eta_t\}_{t=1}^T$, with $\eta_t > 0$ for all $t\in [T]$, is the sequence of the so-called learning rates Cesa-Bianchi2006. We use the standard convention $\omega^{(k)}_{\textnormal{\otextsc{HDG}},1} = 1/K$. Variations of Hedge propose different sequences of learning rates aimed at minimizing the cumulative regret of the forecaster. A natural way to interpret (ref) is to note that it corresponds to the gradient of a scaled LogSumExp function of the negative cumulative losses
where $\eta_t>0$ acts as a smoothing parameter. Large values of $\eta_t$ make the distribution over experts concentrate around the best expert(s) (recovering FTL), while small values yield nearly uniform weights. Equivalently, Hedge arises as the solution to an entropy-regularized linear optimization problem over the simplex, with Shannon entropy acting as a regularizer (see boyd2004convex). From this perspective, Hedge can be viewed as a smoothed version of FTL, with $\eta_t$ controlling the per-round smoothing of the loss.
\paragraph{Constant Hedge.} Setting $\eta_t:=\eta>0$ yields a constant-rate forecaster with a regret bound
where $S_T$ is the maximal per-round loss range, which follows from Hoeffding's lemma. Optimizing (ref) with respect to $\eta$ gives
leading to $\overline{R}_T \le S_T \sqrt{T\log K /2}$. If losses take values in $[0,1]$, then $S_T\le 1$ and optimal $\eta^\star=\sqrt{8\log K/T}$ Littlestone1994. A constant rate allows for the standard closed-form online update via instantaneous losses:
In practice, $\eta$ is a learning hyperparameter: $\eta\to\infty$ recovers FTL; $\eta=0$ yields no learning. The downside of a constant $\eta$ is its lack of adaptivity over arbitrarily long learning periods.
\paragraph{Decreasing Hedge (DecHedge).} {\it Decreasing Hedge} is obtained by setting $\eta_t = c_0 \sqrt{\log K/t}$ for a constant $c_0 > 0$ auer2016anAlgorithm. While this version of Hedge generally requires tuning $c_0$ for a specific learning run (as it is also the case for constant Hedge), mourtadaOptimalityHedgeAlgorithm2019 prove that the choice $c_0 = 2$ is worst-case optimal, without requiring further assumptions on the loss sequences. Decreasing Hedge, too, achieves a worst-case regret bound of $\mathcal{O}(\sqrt{T\log K })$ chernovPredictionExpertAdvice2010. However, in stochastic settings where losses are well-behaved, more refined guarantees are possible. In fact, the generic constant $c_0 = 2$ may be conservative when the environment is not adversarial.
\paragraph{Adaptive Hedge (AdaHedge).} {\it AdaHedge} dynamically tunes $\eta_t$ based on past performance, and neither horizon nor additional parameters are used vanervenAdaptiveHedge2011.
Let $\bm{\omega}_t\in\Delta^{K-1}$ be the weight vector at time $t$, let $\bm{\ell}_t$ be the loss vector of $K$ experts, and let the forecaster's loss be $\overline{\ell}_t$. The mix loss at time $t$ for a learning rate $\eta_t>0$ is
(see, e.g., Cesa-Bianchi2006, Chapter 2). The approximation error $\delta_t := \overline{\ell}_t-\overline{m}_t$ is the mixability gap, with cumulative gap $\nabla_t := \sum_{s=1}^t \delta_s$. Consequently, the regret of the forecaster can be decomposed as
where \(\overline{M}_{t} - \min_{k \in [K]} L^{(k)}_{t}\) is the regret incurred under the mix loss rooijFollowLeaderIf2014. Then, by (ref),
AdaHedge chooses the learning rate at time $t$ by balancing the complexity term ${\log K }/{\eta}$ with the observed cumulative gap $\nabla_{t-1}$,\@ifstar\footnote@star\ofootnote{At time $t$, while the forecaster makes their decision and the current loss is not realized yet, only the $t-1$ cumulative mixability gap is available.} leading to:
with the convention that $\eta_{\textnormal{\otextsc{ah}},t}=\infty$ if $\nabla_{t-1}=0$. The scheme is thus parameter-free and independent of quantities that are available only at terminal time $T$.
If $s_t\le S$ for all $t$, this gives $\overline{R}_{\textnormal{\otextsc{ah}},T}\le S\sqrt{T \log K}+S \left(\frac{4}{3}\ln K+2\right)$, matching the $\mathcal{O}(S\sqrt{T\ln K})$ worst-case rate. Additionally, the leading term depends on the empirical variability $\sum_{t=1}^T s_t^2$ rather than directly on $T$, so the bound tightens in low-variance regimes.
Finally, for completeness, Appendix (ref) also discusses the doubling-trick Hedge algorithm.
Following luoAchievingAllNo2015 and mourtadaOptimalityHedgeAlgorithm2019, one can see that the complexity of the online forecasting problem for a forecaster following a certain scheme can be generally understood in terms of the gaps between the expected loss of the best expert and the expected losses of other experts. Let the per-expert gap be defined as
where $\mu_k = \mathbb{E}[ \ell^{(k)}_t ]$ and $k^\star \in \arg\min_{k \in [K]} \mu_k$, and let
denote the so-called sub-optimality gap between the expected loss of the best expert with respect to the runner-up.
Our first result is a set of regret guarantees for FTL online combination under both independence and mixing assumptions. We apply Hoeffding-type and Bernstein-type bounds to showcase the different kinds of guarantees that can be obtained, which depend on two constants that quantify the dependence of losses over time, and that are necessary ingredients in our proofs. For the sake of space, we do not include here an additional bound obtained using a Bernstein inequality which can be found in Appendix (ref). We emphasize that, even though the presented bounds assume a constant-in-time second moment for losses, they remain valid if one uses time-dependent variances or their uniform-in-time bound.
mourtadaOptimalityHedgeAlgorithm2019 showed that, in the stochastic setting with i.i.d. losses, decreasing Hedge naturally adapts to the hardness of the combination problem, whereas constant and doubling trick Hedge do not.
Let $\chi^{(k)}_t := (\ell^{(k)}_t - \ell^{(k^\star)}_t )$ for each $k \in [K]$, $k \not= k^\star$. The process $\{ \chi^{(k)}_t \}_{t \in \mathbb{Z}}$ is the excess loss of expert $k$ relative to the best expert $k^\star$.
An additional but more complex bound can be constructed using a different Bernstein inequality, which, for clarity of exposition, we derive in Appendix (ref). Again, our Hedge regret bounds are valid for time-dependent variance or for time-independent bounds thereof.
As shown in Section (ref), both Follow-the-Leader and decreasing Hedge algorithms admit constant regret upper bounds, for any terminal time $T$ and number of experts $K$. These guarantees are valid for an arbitrary choice of models-experts. This naturally leads to the question: What, in general, would be the best method for a time series forecasting setting?
The choice of a combination method is non-trivial. Theoretical upper bounds for cumulative regret depend on the sub-optimality gap $\Delta$, which depends on the unknown population mean of losses. Moreover, some methods, such as non-adaptive Hedge, typically require tuning of their learning rate, which could be problematic in practice. However, since in our non-adversarial stochastic setting where mixing dependence conditions are mildly different from independence, we can rely on numerical evidence from the literature to inform our selection. According to the numerical results in mourtadaOptimalityHedgeAlgorithm2019, FTL, decreasing Hedge, and AdaHedge all demonstrate strong performance, both for large $\Delta$-gap scenarios and “hard” stochastic zero-gap ($\Delta =0$) settings, considerably outperforming the constant and doubling-trick Hedge. However, the simulations by vanervenAdaptiveHedge2011 show that constant and decreasing Hedge perform similarly in practice. In Section (ref), we implement constant Hedge as a computationally simple baseline exponential combination method, decreasing Hedge with theoretically optimal scaling ($c_0 = 2$) and AdaHedge, which does not depend on any tuning parameter.
Echo State Networks (ESN) are a particular family of recurrent neural networks with randomly sampled connectivity weights. ESN models are, in general, nonlinear state-space systems that, in the forecasting setting, are defined by the following equations:
where $A \in \mathbb{R}^{D \times D}$ is the reservoir matrix, $C \in \mathbb{R}^{D \times d}$ is the input matrix, ${\zeta} \in \mathbb{R}^D$ is the input shift, $\alpha \in [0, 1)$ is the leak rate, and ${\theta} := (b, W^\top)^\top \in \mathbb{R}^{D+1}$ with $W \in \mathbb{R}^D$ denotes the readout coefficients. The map $\sigma : \mathbb{R} \rightarrow \mathbb{R}$ is an activation function applied elementwise and we assume that $\sigma$ is the hyperbolic tangent in the rest of the paper. We refer to $A$, $C$, and ${\zeta}$ as state parameters. These parameters are randomly drawn and kept fixed throughout, while $\theta$ is the only set of parameters that need to be estimated.
The randomly drawn parameters $A$, $C$, and ${\zeta}$ in (ref) are constructed by first sampling $\widetilde{A}$, $\widetilde{C}$, and $\widetilde{{\zeta}}$ from appropriately chosen distributions, and then normalizing as follows:
where $\rho(\widetilde{A})$ denotes the spectral radius of $\widetilde{A}$. We then define $A = \rho \overline{A}$, $C = \gamma \overline{C}$, and ${\zeta} = \varsigma \overline{{\zeta}}$, which allows to write the state equation (ref) as
The tuple ${\varphi}:= (\alpha, \rho, \gamma, \varsigma)$ is referred to as the hyperparameters of the ESN. Specifically, $\alpha \in [0, 1)$ is the so-called leak rate, $\rho \in \mathbb{R}^+$ denotes the spectral radius of the reservoir matrix $A$, $\gamma \in \mathbb{R}^+$ is the input scaling, and $\varsigma \in \mathbb{R}^+$ corresponds to the shift scaling. These hyperparameters play a crucial role in defining the properties of the state map in (ref). Although there is no clear consensus on how to best select $\varphi$, in this work we construct ensembles of ESN models that not only differ in the draws of the state parameters $A$, $C$, and $\zeta$, but also in the state hyperparameters. This allows us to explore a wide range of different reservoir architectures, and thus sidestep the problem of model tuning.\@ifstar\footnote@star\ofootnote{ballarin2022reservoir propose a general data-driven algorithm to tune $\varphi$. However, the resulting hyperparameter optimization problem is computationally complex and highly non-convex.}
As is common in the ESN literature, we consider the ridge regression estimator for $\theta$. More precisely, let $X := (X_1, \ldots, X_{T-1})^\top \in \mathbb{R}^{(T-1)\times D}$ and $Y := (Y_2, \ldots, Y_T)^\top \in \mathbb{R}^{T-1}$. Define $ M= \big(\mathbb{I}_{T-1} - \frac{1}{T-1} \mathbf{1}_{T-1} \mathbf{1}_{T-1}^\top\big)$, $ \widetilde{Y} = M Y$, and $ \widetilde{X} = M X. $ Then the ridge regression estimator for $\theta$ is given by
where $\lambda \in \mathbb{R}^+$ is the regularization hyperparameter. In practice, to select $\lambda$ we rely on time-series-adapted cross-validation (CV) (see, e.g., bergmeir2012use,bergmeir2018,Hyndman2013).
ballarin2022reservoir introduced a class of Echo State Network models designed for time series observed at multiple sampling frequencies. These models are collectively referred to as {\it Multi-Frequency ESNs} (MFESNs).
The first category of MFESNs, known as {\it Single-reservoir MFESNs} (S-MFESNs), enables the forecasting of lower-frequency targets by mapping high-frequency reservoir states to low-frequency outputs through a state alignment scheme. In an aligned S-MFESN, the forecast is generated based on the most recent state corresponding to the reference low-frequency time index $t$. The second category, termed {\it Multi-reservoir MFESNs} (M-MFESNs), introduces multiple state equations, each associated with a distinct sampling frequency in the input data. For example, quarterly and monthly variables can be incorporated jointly as regressors. By employing multiple reservoirs evolving at their respective frequencies, these models flexibly capture the dynamics of different groups of input variables sampled at common frequencies.
We now briefly introduce the M-MFESN framework following ballarin2022reservoir. For brevity, we do not discuss S-MFESN models explicitly, as they are empirically dominated by multi-reservoir MFESNs (see Section (ref)). Suppose there are $Q$ groups of input time series with observations $\{ Z^{(q)}_{t,{s}\vert{\kappa_q}} \}_{t,s}$, where $Z^{(q)}_{t,{s}\vert{\kappa_q}} \in \mathbb{R}^{d_q}$, $\kappa_q\in \mathbb{N}$, $q \in [Q]$, $t \in [T]$, and $s \in \{0, \ldots, \kappa_q-1\}$, sampled at common frequencies $\{\kappa_1, \ldots, \kappa_Q\}$. Using the tempo indexing of ballarin2022reservoir, the high-frequency observations occur at fixed fractional sub-intervals relative to the reference low-frequency index $t$. For each frequency $\kappa_q $, the fractional index $s=0$ indicates alignment with $t$, so that $(t,{0}\vert{\kappa_q}) \equiv t$, and advancing $\kappa_q$ sub-intervals corresponds to one unit of low-frequency time, that is $(t,{\kappa_q}\vert{\kappa_q}) \equiv (t+1,{0}\vert{\kappa_q}) \equiv t+1$. For each input series group $q$, we define a frequency-specific $D_q$-dimensional reservoir state equation
where $A_q$, $C_q$, and $\zeta_q$ are group-specific normalized reservoir parameters, and ${\varphi}_q:= (\alpha_q, \rho_q, \gamma_q, \varsigma_q)$ is the corresponding hyperparameter vector of the reservoir. To construct the forecast at the reference time $t$, the states are aligned by stacking the frequency-specific states evaluated at $s=0$, that is, $X_{t,Q} := (X^{(1)\top}_{t,{0}\vert{\kappa_1}}, \ldots, X^{(Q)\top}_{t,{0}\vert{\kappa_Q}} )^\top \in \mathbb{R}^{D_Q}$, $D_Q := \sum_{q=1}^Q D_q$. The M-MFESN prediction equation is given by
The readout $W_Q$ defines a linear map from the nonlinear, frequency-specific states to $Y_{t+1}$.
An important concern when working with models initialized or inherently constructed with random weights is the impact that such randomness may have on model performance. A recent and growing body of work on the “lottery ticket” hypothesis frankleLotteryTicketHypothesis2018,malachProvingLotteryTicket2020, argues, informally, that a key component of the empirical success of complex neural networks is “being lucky”, that is, drawing an initial configuration of weights that is favorable maSanityChecksLottery2021,sreenivasanRareGemsFinding2022. Related research has focused explicitly on designing weight distributions that are guaranteed to be a priori advantageous for initialization (see, e.g., zhaoZerOInitializationInitializing2022 and bolagerSamplingWeightsDeep2023). Echo State Networks, although proven to satisfy universal approximation and generalization guarantees RC7,RC10,RC12, also depend critically on the properties of randomly drawn reservoir parameters. Once used for forecasting such systems may offer very different forecasting accuracy. It is hence suggestive to construct {\it ensembles of ESNs} (EN-ESN). In our multifrequency setting we directly extend this idea to {\it ensembles of MFESNs} (EN-MFESN), with single- or multi-reservoir architectures.
Following the notation introduced in Section (ref), consider an ensemble of $K \geq 1$ distinct MFESN models-experts of multiple reservoirs (M-MFESNs), each used to construct the corresponding expert prediction out of $Q$ groups of input series sampled at distinct frequencies as $\widehat{Y}^{(k)}_{t+1} = \widehat{b}^{(k)}_{Q, \lambda} + \widehat{W}^{(k)\top}_Q X^{(k)}_{t,Q}$, $k\in [K]$, with $\widehat{W}^{(k)\top}_Q$ and $\widehat{b}^{(k)}_{Q, \lambda}$ estimated as in (ref)-(ref). The ensemble forecast is constructed according to (ref) as
\paragraph{ESN Ensembles over Random Parameters (EN-MFESN-RP).} Since, as discussed above, the primary source of diversity MFESN instances is the randomness of the state parameters, the first ensemble class we consider is obtained by re-drawing the state parameters for each $k \in [K]$, as in (ref), while keeping all other reservoir hyperparameters $\varphi_q=(\alpha_q, \rho_q, \gamma_q, \varsigma_q)$ fixed (for single- and multi-reservoir cases). This type of reservoir ensemble, which we call {\it MFESN-RP ensemble} (EN-MFESN-RP), aims to span the family of models that differ only in terms of the sampled state parameters. As a result, only moderate variation can be expected between members in terms of autoregressive state dynamics, input scaling, and possible activation function saturation.
\paragraph{ESN Ensembles over Leak Rate and Random Parameters (EN-MFESN-$\alpha$RP).} To further explore the potential gains of the ensemble reservoir approach, we consider varying both the random state parameters and the leak rate $\alpha$. We term this approach as {\it MFESN-$\alpha$RP ensemble} (EN-MFESN-$\alpha$RP). Different leak rates induce different memory characteristics (see RC15), which can strongly affect predictive performance. Combining random reservoir weights with varying leak rates increases ensemble diversity and enhances the ability to capture temporal dependencies at different time scales, hence improving accuracy and robustness. Finally, since there is no known theoretically grounded method of tuning $\alpha$, an EN-MFESN-$\alpha$RP provides a practical alternative: Rather than selecting a single leak rate, the ensemble adaptively identifies effective values in an online, data-driven manner. Our empirical analysis confirms that the leak rate can have a substantial impact on the MFESN forecasting accuracy even when all other architectural choices remain unchanged.
In this section, we evaluate and compare the one-step-ahead forecasting performance of our proposed EN-MFESN against individual MFESN models and state-of-the-art approaches. To ensure a fair comparison, we use the same empirical experimental setup as in ballarin2022reservoir: We use the same data, preprocessing procedures, and benchmark specifications.
We conduct our analysis for the original dataset used in ballarin2022reservoir, which contains two groups of predictors: small-MD, comprising 9 predictors, and medium-MD, consisting of 33 predictors. The data is collected at different sampling frequencies, while the target variable is U.S. quarterly GDP growth. In this paper, we focus on a more challenging medium-MD data sample, while the results for the small-MD dataset can be made available by request. For more in-depth information on specific data sources, data preprocessing and availability, we refer the reader to Appendix (ref) and Section 4.1 in ballarin2022reservoir.
\paragraph{Benchmarks.} The first benchmark we consider is the in-sample mean of U.S. GDP growth computed over the estimation sample, which serves a natural baseline for comparing mean squared forecasting errors. We also include an AR(1) model, standard in univariate time series analysis and macroeconomic forecasting Stock2002, baiForecastingEconomicTime2008. Our analysis demonstrates that AR(1) can outperform na\"ive forecast and poorly calibrated ensembles, making it a non-trivial reference for evaluating multi-frequency nonlinear state-space models. Additionally, the two dynamic factor models (DFMs) from ballarin2022reservoir are also included. DFMs are widely applied in macroeconomic analysis to recover low-dimensional common components from high-dimensional time series stock2016dynamic,Fuleky2020,barigozziDynamicFactorModels2024. \paragraph{MFESNs.} Two classes of ESN-based forecasting models are considered. The first set comprises single-frequency ESN models (S-MFESNs) with reservoir sizes $D=30$ and $D=120$, respectively. The sparsity degree of state matrices is set to $10/D$. The second group of models includes multi-frequency ESNs (M-MFESNs) with separate reservoirs for monthly and daily data input series, allowing to capture the dynamics across both temporal resolutions. For each model class, two variants of MFESNs are considered: “A” and “B” types, which differ, e.g., in the state hyperparameters used. Our implementation of MFESNs follows ballarin2022reservoir in the choice of architecture, reservoir sizes, activation functions, hyperparameters, time-series cross validation for the regularization parameter, and other relevant model-specific details. A summary of our model designs is provided in Appendix (ref). \paragraph{Ensembles.} As outlined in Section (ref), we consider MFESN ensembles with state coefficients resampling (EN-MFESN-RP), as well as ensembles that additionally vary the leak rate value (EN-MFESN-$\alpha$RP). For each class and type of multi-frequency ESN model (single- or multi-reservoir, type A or B), we generate $K$ variants of the “baseline” specification by independently resampling the random reservoir weight matrices. When varying the leak rate, we select $\alpha$ from a finite grid, $\alpha \in \{0.1, 0.3, 0.5, 0.7, 0.9\}$. To ensure comparability between ensemble specifications, a total ensemble size of $K=1000$ models is maintained, allocating an equal number of models to each $\alpha$ value. In the specific case of M-MFESN ensembles, we use the same leak rate for both monthly and daily frequency state equations. We apply the direct combination strategies of simple averaging (SA) and rolling MSE (rollMSE) (see Appendix (ref)) as well as Constant Hedge (Hedge), Decreasing Hedge (DecHedge), Follow-the-Leader (FTL), and the Adaptive Hedge (AdaHedge) of rooijFollowLeaderIf2014 introduced in Section (ref).\@ifstar\footnote@star\ofootnote{We provide pseudo-code algorithms for FTL and Hedge schemes in Appendix (ref).}
We begin by analyzing forecasting performance separately for ensembles without and with leak rate variation, namely EN-MFESN-RP and EN-MFESN-$\alpha$RP ensembles. Our evaluation metric is relative mean-squared forecasting error (MSFE) with respect to AR(1). For each ensemble, we report the empirical cumulative distribution function (ECDF) of relative MSE across all MFESN experts in the ensemble, together with vertical lines indicating the median of the relative errors, and the relative MSFEs of the baseline models and combination schemes.
Figure (ref) presents the results for EN-MFESN-RP ensembles and reveals substantial sensitivity to the random draw of reservoir state coefficients. In particular, individual S-MFESN A/B and M-MFESN A/B model draws taken from ballarin2022reservoir display heterogeneous performance: While B-type single- and multi-reservoir MFESNs perform favorably, A-type models often exhibit higher forecasting errors, in some cases exceeding the median relative MSE with respect to the AR(1) benchmark. Simple averaging and constant Hedge follow each other closely across panels, primarily due to the choice of a small learning rate. In contrast, both FTL and AdaHedge deliver consistent, sizable, and robust performance gains. It is important to notice that Follow-the-Leader often to models in the extreme left tail of the ECDF, indicating near-optimal performance within the ensemble. Although FTL dominates AdaHedge throughout, the difference is higher for S-MFESN models than for M-MFESNs, consistent with their greater flexibility and baseline accuracy of the latter.
The results for EN-MFESN-$\alpha$RP ensembles, displayed in Figure (ref), broadly align with those obtained with EN-MFESN-RP ensembles, while exhibiting further performance gains due to the increased model richness induced by leak-rate variation. Once again, one can see that individual MFESN draws depend significantly on properties of one incidental sampling instance: For example, M-MFESN B has relative MSFE below the median, whereas S-MFESN A remains noticeably above it. Importantly, the left tail of the ECDF changes substantially in some cases (see, e.g., Figure (ref)(c)), indicating that high-performing models (previously unavailable in EN-MFESN-RP) may emerge when leak rate $\alpha$ varies. A further distinction is that no single aggregation method dominates across all EN-MFESN-$\alpha$RP ensembles: FTL and AdaHedge alternate in delivering superior results depending on the MFESN type, as illustrated by the contrast between panels (a)–(c) and (b)–(d) in Figure (ref). This suggests that introducing leak-rate variability fundamentally changes ensemble composition and, in turn, the relative performance of combination schemes.
We summarize our results in Table (ref), where benchmark forecasts are included for reference. For both the median ensemble error and each combination method, we report the relative MSFE (with respect to the in-sample mean) together with percentage improvements relative to the original MFESN model instance used in ballarin2022reservoir. Overall, both FTL and AdaHedge are robust in achieving significant gains over the original MFESNs, with FTL appearing as the most effective scheme. These results agree well with Section (ref) which drew on the theoretical and simulation-based literature on online learning. In the EN-MFESN-RP setting, online ensemble combination yields forecasting error reduction up to 31% within a given MFESN class; for the best model type, M-MFESN B, FTL achieves reductions exceeding 41% and 32% MSFE over the in-sample mean and AR(1), respectively. Even higher gains arise in the EN-MFESN-$\alpha$RP setting. For S-MFESN A models, AdaHedge reduces MFSE by more than 42%, while for M-MFESN B, FTL achieves MSFE reductions exceeding 51% and 45% relative to the in-sample mean and AR(1) benchmarks, respectively. We discuss the weight evolution of ensembles for FTL amd AdaHedge algorithms in Appendix (ref).
In this subsection we examine the effectiveness of different model groups at a more granular level than the aggregate analysis in Section (ref). For the $\alpha$RP ensembles, the total budget of $K=1000$ models is partitioned into subsets of $200$ MFESN specifications, each corresponding to a distinct leak rate. This allows us to explore whether there is substantial variability in forecasting performance across these subsets.
Figure (ref) displays the overlaid ECDFs of EN-MFESN-$\alpha$RP ensemble models grouped by leak rate. The plots reveal pronounced variability across both the leak rates $\alpha$ and model types. For S-MFESNs (types A and B), a leak rate of $\alpha = 0.7$ consistently dominates, whereas both low and high values ($\alpha = 0.1$ and $\alpha = 0.9$) perform considerably worse than the ensemble as a whole. The behavior of M-MFESN ensembles is more nuanced: For both A- and B-type models, $\alpha=0.9$ produces the worst results by a wide margin, while low leak rates ($\alpha \in \{0.1, 0.3\}$) degrade the high-performing experts and shrink the upper tail, compressing the distribution around the median. In contrast, $\alpha = 0.5$ and $\alpha = 0.7$ yield the best overall results. The first value uniformly dominates the total $\alpha$RP ECDF, while the second produces a pronounced rise up to approximately the 3$^\text{rd}$ decile, indicating that a substantial fraction of models populate the lower MSFE tail of the error distribution with this leak rate value.
These findings suggest revisiting the leak-rate choices in ballarin2022reservoir. Across MFESN classes, the most effective leak rates differ substantially from the original specifications. For instance, baseline S-MFESNs are constructed with $\alpha = 0.1$ (Table (ref)), while Figure (ref)(a)-(b) suggest the values around $0.7$ are a uniformly more robust choice. This insight, however, comes at a cost. Attempting to study multiple hyperparameters with a comparable level of precision, such as distinct monthly and daily frequency reservoir leak rates in M-MFESNs, would require the ensemble budget to grow exponentially. Model combination therefore offers a parsimonious and adaptive alternative, relying on a fixed-budget set of models rather than the construction of a single optimal specification. This adaptivity is particularly valuable for MFESNs, where randomness induces persistent model heterogeneity that cannot be corrected ex-post, a feature shared by many modern neural network architectures.
The ensemble combination approach is a well-studied framework that enables a forecaster to incorporate the predictions of a pool of models into a single, adaptive forecast. However, these techniques remain under-explored in the fields of applied statistics and econometrics.
In this setting, our work makes contributions in both theoretical and empirical directions. First, we provide concentration bounds for stochastic ensemble losses under i.i.d. and $\varphi$-mixing assumptions, exploiting Hoeffding- and Bernstein-type results. Then, we apply these results to Follow-the-Leader and decreasing Hedge methods, for which we derive finite regret bounds. We then consider an application where model combination allows for robust and substantial performance improvements against previous state-of-the-art results: Macroeconomic forecasting with mixed-frequency data using multi-frequency ESN models. We find that the ensemble approach is extremely effective at sharpening the performance of all MFESN models, with Follow-the-Leader and AdaHedge schemes being the overall most effective schemes. Our empirical analysis also reveals that combinations can efficiently sidestep the hyperparameter tuning, which for (MF)ESN models is highly non-trivial.