The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
131,048 characters
Reservoir Computing for Macroeconomic Forecasting with Mixed Frequency Data
\title{\vspace{-25pt}\huge Reservoir Computing for Macroeconomic Forecasting\\ with Mixed Frequency Data}
\author{Giovanni Ballarin$^{1}$, Petros Dellaportas$^{2,3}$, Lyudmila Grigoryeva$^{4,5}$, \\ Marcel Hirt$^{8}$, Sophie van Huellen$^{6,7}$, Juan-Pablo Ortega$^{8}$}
\date{ \ }
\maketitle
\vspace{-20pt}
\begin{abstract}
Macroeconomic forecasting has recently started embracing techniques that can deal with large-scale datasets and series with unequal release periods. MIxed-DAta Sampling (MIDAS) and Dynamic Factor Models (DFM) are the two main state-of-the-art approaches that allow modeling series with non-homogeneous frequencies. We introduce a new framework called the Multi-Frequency Echo State Network (MFESN) based on a relatively novel machine learning paradigm called reservoir computing. Echo State Networks (ESN) are recurrent neural networks formulated as nonlinear state-space systems with random state coefficients where only the observation map is subject to estimation. MFESNs are considerably more efficient than DFMs and allow for incorporating many series, as opposed to MIDAS models, which are prone to the curse of dimensionality. All methods are compared in extensive multistep forecasting exercises targeting US GDP growth. We find that our MFESN models achieve superior or comparable performance over MIDAS and DFMs at a much lower computational cost.
\end{abstract}
\bigskip
\textbf{Key Words:} Reservoir Computing, Echo State Networks, Forecasting, US Output Growth, GDP, Mixed-frequency Data, Time Series, Multi-Frequency Echo State Network, MIDAS, DFM \\
\textbf{JEL:} C53, C45, E17
\vspace{1em}
{\footnotesize\noindent {\bf Acknowledgments:} The authors acknowledge financial support from UK Research and Innovation (grant number ES/V006347/1). GB thanks the hospitality and the generosity of the University of Warwick without which this paper would not have happened. GB is grateful to the Center for Doctoral Studies in Economics (CDSE) at the University of Mannheim for providing travel funding through the 2020 CDSE Teaching Award. The authors acknowledge support by the state of Baden-Württemberg through bwHPC. The authors acknowledge the use of the UCL Myriad High Performance Computing Facility (Myriad@UCL), and associated support services, in completing this work. MH acknowledges that the majority of his work was done when he was affiliated with UCL and hence thanks the support of UCL in this project. The authors acknowledge the support of the Swiss National Science Foundation (grant number 200021\_175801/1). The authors also thank Thanos Moraitis and Alfonso Silva Ruiz for excellent research assistance and help in data collection and curation. }
\addtocounter{footnote}{1} \footnotetext{
Department of Economics, University of Mannheim, L7, 3-5, Mannheim, 68131, Germany. {\texttt{[email removed]} }}
\addtocounter{footnote}{1} \footnotetext{
Department of Statistical Science, UCL, Gower Str., London WC1E 6BT, UK. {\texttt{[email removed]} }}
\addtocounter{footnote}{1} \footnotetext{
Department of Statistics, Athens University of Economics and Business, 10434 Athens, Greece. {\texttt{[email removed]} }}
\addtocounter{footnote}{1} \footnotetext{
Faculty of Mathematics and Statistics, University of St. Gallen, Bodanstrasse 6, CH-9000 St.~Gallen, Switzerland. {\texttt{[email removed]} }}
\addtocounter{footnote}{1} \footnotetext{
Honorary Associate Professor, Department of Statistics, University of Warwick, Coventry CV4 7AL, UK. {\texttt{[email removed]} }}
\addtocounter{footnote}{1} \footnotetext{
Global Development Institute (GDI), University of Manchester, Manchester M13 9PL, UK.
{\texttt{[email removed]} }}
\addtocounter{footnote}{1} \footnotetext{
Department of Economics, SOAS University of London, WC1H 0XG, London, UK.
{\texttt{[email removed]} }}
\addtocounter{footnote}{1} \footnotetext{
Division of Mathematical Sciences,
Nanyang Technological University,
21 Nanyang Link,
Singapore 637371.
\texttt{[email removed]} (some of the work done while at UCL), {\texttt{[email removed]}}}
\makeatother
\newpage
\setcounter{footnote}{0}
\newpage
\setstretch{1.1}
\section{Introduction}
\label{Introduction}
The availability of timely and accurate forecasts of key macroeconomic variables is of crucial importance to economic policymakers, businesses, and the banking sector alike. Fundamental macroeconomic figures, such as GDP growth, become available at low frequency with a considerable time lag and are subject to various rounds of revisions after their release. This is particularly problematic in a fast-changing and uncertain economic environment, as experienced during the Great Recession of 2007-2008 (\cite{Hindrayanto2016}) and the recent pandemic (\cite{Buell2021, Huber2021}). However, a large number of the potentially predictive financial market (and other macroeconomic) indicators are available at a daily or even higher frequency (\cite{Andreou2013}). The desire to utilize such high-frequency data for macroeconomic forecasting has led to the exploration of techniques that can deal with large-scale datasets and series with unequal release periods (see \cite{Borio2011, Borio2013,Morley2015}; we also refer the reader to \cite{Fuleky2020} for more details regarding high-dimensional data and to \cite{Armesto2010} and \cite{Banbura2013} for a review on mixed-frequency data).
We contribute to the existing literature by proposing a new macroeconomic forecasting framework that utilizes high-dimensional and mixed-frequency input data, the Multi-Frequency Echo State Network (MFESN). The MFESN originates from a machine learning paradigm called Reservoir Computing (RC). RC is a family of learning models that take advantage of the information processing capabilities of complex dynamical systems~(see \cite{maass1, DynamicalSystemsMaass, Crutchfield2010}, and \cite{lukosevicius, tanaka:review} for reviews). Generally speaking, RC is a versatile class of recurrent neural network (RNN) models (see \cite{Salehinejad2017} for a detailed survey). Although conventional RNNs are well-suited for handling sequence data and dynamic problems, estimating their weights during the training phase is inherently difficult~(\cite{Pascanu2013a, Doya92}). Reservoir networks stand out due to the fact that their inner weights can be {\it randomly generated} and {\it fixed}, and only the output (readout) layer weights are subject to estimation (supervised training). Echo State Network (ESN) is one of the most popular instances of RC models with provable universality, generalization properties (see \cite{RC6, RC7, RC9, RC10, RC12, RC20}, and references therein for more details), and excellent performance in forecasting, classification, and learning of dynamical systems (see \cite{Hart2021,RC21}). While conventional RNNs have been adopted for macroeconomic forecasting in a few instances (see, for example, \cite{Paranhos2021}), to the best of our knowledge, we are the first to explore easily-trainable reservoir models in this context.
Our main contribution is three-fold. First, inspired by the remarkable empirical success of ESNs in prediction tasks, we propose the so-called Multi-Frequency Echo State Network (MFESN) framework, which allows multistep forecasting of the target variable at lower or the same frequencies as those of the input series. Second, we introduce two different approaches to predicting within the MFESN framework, namely {\it Single-Reservoir MFESN} (S-MFESN) and {\it Multi-Reservoir MFESN} (M-MFESN). S-MFESN is determined by modifying the ESN architecture to accommodate input and target variables of mixed frequencies. In M-MFESN, several Echo State Networks are adopted to handle input time series, each ESN corresponding to a group of input variables quoted at one given frequency. Finally, our third contribution consists of an extensive empirical comparative analysis of the forecasting capability of the proposed approaches in a concrete task of predicting the quarterly U.S. output growth. We inspect the forecasting capabilities of the MFESN framework compared to two well-established benchmarks widely used in the macroeconomic literature and among practitioners and show its empirical superiority in several thoroughly conducted forecasting exercises. Moreover, as a bi-product, we propose a new data aggregation scheme that allows bridging these two standard forecasting approaches, which is not available in the literature.
In our empirical study, we evaluate the multistep forecasting performance of the MFESN framework targeting quarterly U.S. output growth (Gross Domestic Product (GDP) growth) and utilizing a small- and medium-sized set of monthly and daily financial and macroeconomic variables. We compare the MFESN approach against two state-of-the-art methods, MIDAS and DFM, known for their ability to incorporate data of heterogeneous frequencies and utilize high-dimensional data inputs. The MIxed DAta Sampling (MIDAS) model developed in \cite{Ghysels2004,Ghysels2007} has been adopted widely for macroeconomic forecasting with mixed-frequency data (see for instance \cite{Clements2008, Clements2009, Ghysels2009, Francis2011, Monteforte2012, Galvao2010, Galvao2013, Andreou2013, Ghysels2016, Jardet2020}). However, MIDAS is prone to curse-of-dimensionality problems and performs poorly when the set of predictors is of even moderate size (\cite{Clements2009,Kostrov2021essays}) due to optimization-related issues. Recently, some attempts have been made in the literature to overcome these issues by employing variable selection techniques under some additional assumptions. For instance, \cite{Babii2022} proposes the MIDAS projection approach, which is more amenable to high-dimensional data environments under the assumption of sparsity. Even with these improvements, practical high-dimensional implementations of MIDAS remain challenging. This is in part caused by the ragged edges of the ``raw'' macroeconomic data, incomplete observations, and uneven sampling frequencies. The relative inflexibility of MIDAS regression lag specifications makes integrating daily and weekly data at true calendar frequencies (that is, without interpolation or aggregation) very complex. State-space models effectively mitigate these issues.
A strong state-of-the-art state-space competitor for our MFESN framework is the Dynamic Factor Model (DFM), which has been first introduced in \cite{geweke1977dynamic} and \cite{sargent1977business}. DFMs have become the standard workhorse for macroeconomic nowcasting and prediction (for more details, we refer the reader to \cite{Stock1996, Stock2002, stock2016dynamic, Giannone2008, Banbura2011, Chauvet2015, Hindrayanto2016}). Conventional DFMs for data of multiple sampling frequencies are linear state-space models with a latent low-frequency process of interest and high-dimensional input time series. Although their linear structure lends itself to inference with likelihood-based methods and Kalman filtering, using DFMs in the high-dimensional setting is limited by the associated computational effort. For Gaussian state-space models, some of these issues are proposed to be handled with a more compact matrix representation as in \cite{DELLEMONACHE201922}. Still, in the particular settings of nowcasting and forecasting of GDP growth, the computational complexity is one of the main reasons why DFMs are rarely used with daily input series, see \cite{Banbura2013} for a detailed review and \cite{Aruoba2009} for a mixed-frequency DFM wherein the latent factor process is updated daily, with the highest input frequency being weekly. We address these numerical difficulties using novel Python libraries for auto-differentiation and using GPUs for parallel computing, which allow the estimation of DFMs even in instances of high-frequency input observations. Further, to adapt the DFM to mixed frequency tasks, we propose a new DFM aggregation scheme with Almon polynomial structure that bridges MIDAS and the DFM for our forecasting comparison. To our knowledge, we are the first to present this aggregation scheme which reduces the number of parameters subject to estimation. In contrast, previous DFM such as in \cite{mariano2003new, Banbura2011, Camacho2010a, Frale2011a} commonly assume a fixed aggregation scheme a-priori depending on whether the macroeconomic variable is a flow or stock variable.
To carry out a fair comparison of our MFESN framework with the state-of-the-art MIDAS and DFM models, we designed two model evaluation settings that differ regarding whether the financial crisis of 2007-2008 is included in the estimation period or not. In the first forecasting setting, all the competing models are estimated using the data from January 1st, 1990, until December 31st, 2007. Their performance in the forecasting into and after the financial crisis period is assessed. In the second evaluation setting, fitting is done with data largely encompassing the crisis period, again from January 1st, 1990 but now up to December 31st, 2011. In both cases, the forecasting (testing) period spans time up to the COVID-19 pandemic events, namely the fourth quarter of 2019.
Along with the two state-of-the-art DFM and MIDAS models, we use the unconditional mean of the sample as a baseline benchmark against the reservoir models. We find that our ESN-inspired models attain comparable or much better performance than DFMs at a much lower computational cost, even for a relatively long forecasting horizon of four quarters. Additionally, ESNs do not suffer from curse-of-dimensionality problems, which are known to be pervasive for MIDAS models and hence consistently outperform them in a number of forecasting exercises.
The remainder of the paper is structured as follows. Section~\ref{Reservoir models} presents reservoir models and discusses their advantages, as well as estimation, hyperparameter tuning, penalization and nonlinear multistep forecasting. In Section~\ref{Multi-Frequency Echo State Models}, we introduce the Multi-Frequency Echo State Network (MFESN) framework, propose the single-reservoir and multi-reservoir MFESN models, and spell out their defining features. Section~\ref{sec:empirical_study} contains the empirical study of the comparative GDP forecasting performance of MFESNs with respect to the set of benchmark models. We assess one-step and multistep forecasting results in several setups, with a small and a medium-sized set of regressors. We fit models with data before and after the 2007-08 financial crisis, and with different estimation windows. Section~\ref{Conclusions} concludes and discusses future research avenues and applications. Finally, the Appendix contains information regarding data sources, forecasting figures and formal details regarding our forecasting setups. The Supplementary Appendix gives detailed information on the implementation of all models, robustness checks and provides additional figures.
\paragraph{Code.} Our code, the data, and all results presented in the paper are made available in the GitHub repository at {\href{https://github.com/RCEconModelling/Reservoir-Computing-for-Macroeconomic-Modelling}{github.com/RCEconModelling/Reservoir-Computing-for-Macroeconomic-Modelling}}.
\subsection{Notation}
\label{Notation}
We use the symbol $\mathbb{N}$ (respectively, $\mathbb{N}^+$) to denote the set of natural numbers with the zero element included (respectively, excluded). $\mathbb{Z}$ denotes the set of all integers. We use $\mathbb{R}$ (respectively, $\mathbb{R}_+$) to denote the set of all (respectively, positive excluding zero element) reals. We abbreviate the set $[n]=\{1, \ldots, n\}$, with $ n \in \mathbb{N}^+$.
\paragraph{Vector notation.}
A column vector is denoted by a bold lowercase symbol like $\boldsymbol{r}$ and $\boldsymbol{r} ^\top $ indicates its transpose. Given a vector $\boldsymbol{v} \in \mathbb{R} ^n $, we denote its entries by $v_i$, with $i \in \left\{ 1, \dots, n
\right\} $; we also write $\boldsymbol{v}=(v _i)_{i \in \left\{ 1, \dots, n\right\} }$. The symbols $\boldsymbol{i} _n, \boldsymbol{0} _n \in \mathbb{R} ^n $ stand for the vectors of length $n$ consisting of ones and of zeros, respectively. Additionally, given $n \in
\mathbb{N}^+ $, $\boldsymbol{e}_n ^{(i)} \in \mathbb{R} ^n$, $i \in \left\{ 1, \dots, n\right\} $ denotes the
canonical unit vector of length $n$ determined by $\boldsymbol{e} _n ^{(i)} = (\delta _{ij} )_{j\in \left\{ 1, \dots, n\right\} }$. For any $\boldsymbol{v} \in \mathbb{R} ^{n} $, $\|\boldsymbol{v}\| $ denotes its Euclidean norm.
\paragraph{Matrix notation.}
We denote by $\mathbb{M}_{n , m }$ the space of real $n\times m$ matrices with $m, n \in \mathbb{N}^+ $. When $n=m$, we use the symbols $\mathbb{M}_n $ and $\mathbb{D} _n $ to refer to the space of square and diagonal matrices of order
$n$, respectively. Given a matrix $A \in \mathbb{M} _{n , m} $, we denote its components by $A _{ij} $ and we write $A=(A_{ij})$, with $i \in \left\{ 1, \dots, n\right\} $, $j \in \left\{ 1, \dots m\right\} $. The symbol $\mathbb{I}_n \in \mathbb{D} _n $ denotes the identity matrix, and the symbol $\mathbb{O}_n$ stands for the zero matrix of dimension $n$. For any $A \in \mathbb{M} _{n , m} $, $\|A\| _2 $ denotes its matrix norm induced by the Euclidean norms in $\mathbb{R}^m $ and $\mathbb{R} ^n $, and $\|A\| _2=\sigma_{{\rm max}}(A)$, with $\sigma_{{\rm max}}(A)$ the largest singular value of $A$.
\paragraph{Input and target stochastic processes.}
We fix a probability space $(\Omega,\mathcal{A},\mathbb{P})$
on which all random variables are defined.
The input and target signals are modeled by discrete-time stochastic processes $\boldsymbol{ z} = (\boldsymbol{ z}_t)_{t \in \mathbb{Z}}$ and $\boldsymbol{ y} = (\boldsymbol{ y}_t)_{t \in \mathbb{Z}}$ taking values in $\mathbb{R}^K$ and $\mathbb{R}^J$, respectively. Moreover, we write $\boldsymbol{ z}(\omega)=(\boldsymbol{ z}_t( \omega))_{t \in \mathbb{Z}}$ and $\boldsymbol{ y}(\omega)=(\boldsymbol{ y}_t( \omega))_{t \in \mathbb{Z}}$ for each outcome $\omega \in \Omega$ to denote the realizations or sample paths of $\boldsymbol{ z}$ and $\boldsymbol{ y}$, respectively. Since $\boldsymbol{ z}$ can be seen as a random sequence in $\mathbb{R}^K$, we write interchangeably $\boldsymbol{ z}:{\mathbb{Z}} \times \Omega \longrightarrow \mathbb{R}^K$ and $\boldsymbol{ z}: \Omega\longrightarrow (\mathbb{R}^K)^{\mathbb{Z}}$.
The same applies to the analogous assignments involving $\boldsymbol{y}$.
\paragraph{Temporal notation.}
Let $(u_t)_{t \in I}$, $u_t \in \mathbb{R}$ be a (scalar) time series with $I$ some index set (in this paper it will always be discrete). Time series $(u _t)_{t \in I}$ will be denoted just as $(u _t)$ when the index set $I$ is specified by the context. We write $u_{s_1:s_2}=(u_t)_{t\in \{s_1, \ldots, s_2\}}$ for integers $s_1<s_2$ and time series $(u_t)$. To define the concept of the sampling frequency, we must introduce an additional series, call it $(v_s)_{s\in J} $. The time index $J$ is not the same as $I$. We assume that $u_t$ is sampled at the coarsest rate; equivalently, it has the \textit{lowest} sampling frequency, which we call in what follows the \textit{reference frequency}. In practice, this means that in the same window of time, $u_t$ will be observed at most as frequently as $v_s$. The case when the sampling frequency of $v_s$ is strictly higher than that of $u_t$ is of primary interest.
We assume that all sampling happens in instants that are evenly spaced in time. Series other than the reference one and with higher sampling frequencies are given an additional time index, the \textit{tempo index}, written $t, {\ast}\vert{{\kappa}} \,$, where $\kappa$ is the \textit{frequency multiplier}. Our tempo notation assumes that low- and high-frequency series are sampled with temporal \textit{alignment}: this means that the reference time index $t$ and the tempo index ${\ast}\vert{{\kappa}}$ have the following properties.
\begin{definition}
\label{Temporal notation}
A reference time index $t \in \mathbb{N}$ and a tempo index ${\ast}\vert{{\kappa}}$ for a given high-frequency ${\kappa} \in \mathbb{N}^+$ are such that the following relations hold
\begin{description}
\item[(i)] $t, {0}\vert{{\kappa}} \equiv t $
\item[(ii)] $t, {{\kappa}}\vert{{\kappa}} \equiv t+1$
\item[(iii)] $t, {s}\vert{{\kappa}} \equiv t + \lfloor{s / {\kappa}}\rfloor, {(s{\mskip 3mu}\textnormal{mod}{\mskip 3mu}{\kappa})}\vert{{\kappa}} \quad {\rm for} \enspace \forall s \in \mathbb{N}$
\item[(iv)] $t, {-s}\vert{{\kappa}} \equiv (t - 1) - \lfloor{s / {\kappa}}\rfloor, {{\kappa} - (s{\mskip 3mu}\textnormal{mod}{\mskip 3mu}{\kappa})}\vert{{\kappa}} \quad {\rm for} \enspace \forall s \in \mathbb{N}$,
\end{description}
where ${\mskip 3mu}\textnormal{mod}{\mskip 3mu}$ is the modulo operation and for any $x\in \mathbb{R}$ the floor operator $\lfloor x \rfloor$ outputs the greatest $z \in \mathbb{N}$ such that $z\leq x$.
\end{definition}
Since we can exchange ``frequency'' and ``frequency multiplier'' in the tempo notation, we will make no distinction between the two terms in what follows.
\paragraph{Forecasting schemes.}
The theoretical setup and design of the forecasting exercises conducted in this paper are carefully discussed in Appendix~\ref{Forecasting Schemes}. There, we formally distinguish between the so-called high-frequency and low-frequecy forecasting in the presence of mixed-frequency data. For more details regarding time series forecasting with economic data, we also refer the reader to \cite{Clements2008,Clements2009,Chen2010,Jardet2020} and references therein.
\section{Reservoir Models}
\label{Reservoir models}
In this section, we introduce {\it reservoir computing} models \citep{Jaeger04} for forecasting of stochastic time series of a single frequency. We focus on a family of RC systems called {\it Echo State Networks} (ESNs), which have been successfully applied to forecasting of deterministic dynamical systems \citep{pathak:chaos, Pathak:PRL, wikner2021using, arcomano2022hybrid}. In the following, we discuss the linear estimation of ESN model parameters, the hyperparameters tuning, the loss penalty selection, and how to carry out nonlinear forecasting.
\subsection{Reservoir Models}
\label{subsec:reservoir_models}
Reservoir computing (RC) models are nonlinear state-space systems that, in the forecasting setting, are defined by the following equations:
\begin{align}
\boldsymbol{x}_t & = F(\boldsymbol{x}_{t-1}, \boldsymbol{z}_t), \label{state_map}\\
\boldsymbol{y}_{t+1} & = h_{\boldsymbol{\theta}}(\boldsymbol{x}_t) + \boldsymbol{\epsilon}_t, \label{output_map}
\end{align}
for all $t \in \mathbb{Z}$, where the {\it state map} $F: \mathbb{R}^N \times \mathbb{R}^K \to \mathbb{R}^N$, $N, K \in \mathbb{N}^+$ is called also the \textit{reservoir map}, and the \textit{observation map} $h_{\boldsymbol{\theta}}: \mathbb{R}^N \to \mathbb{R}^J$, $J\in \mathbb{N}^+$ is referred to as the \textit{readout} layer, parametrized by $\boldsymbol{\theta} \in \Theta$.
Sequences $(\boldsymbol{z}_t)_{t \in \mathbb{Z}}$, $\boldsymbol{z}_t \in \mathbb{R}^K$, and $(\boldsymbol{y}_t)_{t \in \mathbb{Z}}$, $\boldsymbol{y}_t \in \mathbb{R}^J$, stand for the \textit{input} and the \textit{output (target)} of the system, respectively, and $(\boldsymbol{x}_t)_{t \in \mathbb{Z}}$, $\boldsymbol{x}_t \in \mathbb{R}^N$, are the associated \textit{reservoir states}.
In \eqref{output_map}, $(\boldsymbol{\epsilon}_t)_{t \in \mathbb{Z}}$ are $J$-dimensional independent zero-mean innovations with variance $\sigma_\epsilon^2\mathbb{I}_J$ that are also independent of $\boldsymbol{x}_t$ across all $t$.
Importantly, many families of RC systems have been proven to have universal approximation properties for $L^p$-integrable stochastic processes \citep{RC8}, and estimation and generalization error bounds have been established in \cite{RC10, RC12}.
In the case of an ESN model, the state and observation equations \eqref{state_map}-\eqref{output_map} are given by
\begin{align}
\boldsymbol{x}_t & = \alpha \boldsymbol{x}_{t-1} + (1-\alpha) \sigma(A \boldsymbol{x}_{t-1} + C \boldsymbol{z}_t + \boldsymbol{\zeta}) \label{eq:esn_model_states} \\
\boldsymbol{y}_{t+1} & =\boldsymbol{a} + W^\top \boldsymbol{x}_t + \boldsymbol{\epsilon}_t, \label{eq:esn_model_output}
\end{align}
where $A \in \mathbb{M}_{N}$ is the \textit{reservoir matrix}, $C \in \mathbb{M}_{N,K}$ is the \textit{input matrix}, $\boldsymbol{\zeta} \in \mathbb{R}^N$ is the \textit{input shift}, $\alpha \in [0,1)$ is the \textit{leak rate} and $W \in \mathbb{M}_{N,J}$ are the \textit{readout coefficients}. The map $\sigma: \mathbb{R} \to \mathbb{R}$ is an activation function applied elementwise, which in what follows we take to be the hyperbolic tangent. We refer to $A$, $C$, $\boldsymbol{\zeta}$ as \textit{state parameters} that are randomly generated. Notice that if $A = 0$ and $\alpha = 0$ the state equation reduces to a nonlinear regression model with random coefficients (or a feedforward neural network with random weights) which is usually referred to as an {\it Extreme Learning Machine} \citep{caoReviewNeuralNetworks2018, RC12}.
\paragraph{Properties of ESN models.}
We focus on ESNs with the so-called {\it echo state property (ESP)}, that is, when for any $ \boldsymbol{z} \in (\mathbb{R}^K)^{\mathbb{Z}}$ there exists a unique $ \boldsymbol{y} \in (\mathbb{R}^J)^{\mathbb{Z}}$ such that \eqref{eq:esn_model_states}-\eqref{eq:esn_model_output} hold (see \cite{RC6, RC7, RC9} and references therein). One can require that the ESP holds only on the level of the state equation, that is for any input sequence $ \boldsymbol{z} \in (\mathbb{R}^K)^{\mathbb{Z}}$ there exists a unique state sequence $ \boldsymbol{x} \in (\mathbb{R}^N)^{\mathbb{Z}}$ such that \eqref{eq:esn_model_states} holds. The result in Corollary 3.2 in \citet{RC7}, which is also valid for the case of ESNs with the leak rate, shows that the sufficient condition of the ESP associated with \eqref{eq:esn_model_states} to hold is $\left\|A\right\|_2 L _\sigma<1 $ where $L _\sigma$ is the Lipschitz constant of the activation function $\sigma$ (in our setting, $L _{tanh}=1$). This sufficient ESP condition has been extensively studied in the ESN literature; see \cite{jaeger2001, Jaeger04, Buehner:ESN, zhang:echo, Yildiz2012, Wainrib2016, Manjunath:Jaeger} for more details. The result in Corollary 3.2 in \citet{RC7} also shows that this condition implies the so-called {\it fading memory property} \citep{Boyd1985}, which from the practical point of view means that the impact of initial $\boldsymbol{x}_0$ is negligible for sufficiently long samples.
In the stochastic setting, part (i) of Proposition 4.2 in \citet{RC16} proves that the condition $\left\|A\right\|_2<1$ guarantees variance stationarity of the states associated with variance stationary inputs. Moreover, \cite{RC27} show that this condition implies the so-called stochastic state contractivity ensuring a stochastic analog of the ESP. Notably, violations of $\left\|A\right\|_2 <1$ do not have detrimental implications for the performance of ESNs in various learning tasks, as reported in multiple empirical studies.
\paragraph{Computational advantages of ESNs.}
We emphasize that the core computational advantage of ESNs is that state parameters $A$, $C$, and $\boldsymbol{\zeta}$ are randomly sampled and need not be estimated. Additionally, since observation equation \eqref{eq:esn_model_output} is linear in $\boldsymbol{x}_t$, coefficients $W$ can be estimated via (penalized) least squares regression, as we explain in the following subsection. The choice of properties of state parameters determines memory properties and forecasting performance of linear \citep{RC23} and nonlinear ESNs \citep{RC15} as we discuss in Section~\ref{subsec:Hyperparameter Tuning}.
\subsection{Estimation}
\label{Estimation}
We now discuss in detail the estimation of coefficients $W$ in \eqref{eq:esn_model_output}. Let a sample $(\boldsymbol{z}_t, \boldsymbol{y}_t)_{t=1}^T$ of input and target pairs be available. Given an initial state $\boldsymbol{x}_0$, the reservoir states can be computed iteratively according to state equation \eqref{eq:esn_model_states} as:
\begin{equation*}
\boldsymbol{x}_1
= \alpha \boldsymbol{x}_{0} + (1-\alpha) \sigma(A \boldsymbol{x}_{0} + C \boldsymbol{z}_1 + \boldsymbol{\zeta}) ,
\quad \ldots, \quad
\boldsymbol{x}_T
= \alpha \boldsymbol{x}_{T-1} + (1-\alpha) \sigma(A \boldsymbol{x}_{T-1} + C \boldsymbol{z}_T + \boldsymbol{\zeta}) .
\end{equation*}
Collect the states and the targets into the state and the observation matrices, respectively, as
\begin{equation*}
X = (\boldsymbol{x}_1, \boldsymbol{x}_2, \ldots, \boldsymbol{x}_{T-1})^\top \in \mathbb{M}_{T-1,N},
\qquad
Y = (\boldsymbol{y}_2, \boldsymbol{y}_3, \ldots, \boldsymbol{y}_T)^\top \in \mathbb{M}_{T-1,J} .
\end{equation*}
Consider the ridge regression estimator for $W$ given by
\begin{equation}\label{eq:ridge_regression}
\widehat{W}_\lambda
\: := \:
\operatorname*{arg\,min}_{W \in \mathbb{R}^N} \sum_{t=1}^{T-1} \left\lVert \boldsymbol{y}_{t+1} - {W}^\top \boldsymbol{x}_t \right\rVert_2^2 + \lambda \lVert {W} \rVert_2^2
\: = \:
\left( X^\top X + \lambda ((T-1)\, \mathbb{I}_N) \right)^{-1} X^\top Y ,
\end{equation}
where $\lambda \in \mathbb{R}_+$ is the ridge penalty strength. When $\lambda \to 0$, the estimator $\widehat{W} _\lambda$ converges to the minimum-norm least squares solution \citep{ishwaranGeometryPropertiesGeneralized2014b}.
In applications, ridge regression is the most commonly used estimation method applied to ESNs, as it provides a straightforward regularization scheme both when $N < T$ and $N \geq T$. This is especially important since in practice the ESN state dimension is often chosen to be $10^3$--$10^4$ (see for example \cite{pathak:chaos}). Additionally, a virtue of the ridge regression problem is the fact that the associated objective function is convex and, hence, it can be efficiently solved using stochastic gradient descent even when ${\rm min}\{N, T\}$ is large and one decides against the closed-form solution~\eqref{eq:ridge_regression}. Finally, as mentioned in the properties of reservoir systems in Subsection~\ref{subsec:reservoir_models}, we notice that in the presence of the fading memory property, the estimation does not depend significantly on the choice of $\boldsymbol{x}_0$ as sample size $T$ increases.
We refer to \eqref{eq:ridge_regression} as the \textit{fixed-parameter} estimator. In our empirical analyses, we also implement \textit{expanding} and \textit{rolling window} estimation strategies which update $\widehat{W}_\lambda$ as new observations become available (we refer the reader to Supplementary Appendix \ref{sec:estimation_setup} for details). In the rest of the paper, for brevity, we use $\widehat{W}$ to denote the ridge estimator of coefficients $W$ assuming that the appropriate choice of the penalty strength $\lambda$ is made for each concrete situation.
\subsubsection{Hyperparameter Tuning}
\label{subsec:Hyperparameter Tuning}
As discussed in Subsection~\ref{subsec:reservoir_models}, the performance of ESNs depends on the choice of randomly drawn state parameters $A$, $C$, $\boldsymbol{\zeta}$. Much work has been put into determining optimal specifications (see for example \cite{Rodan2011, Goudarzi2016, farkas:bosak:2016, GHLO2014_capacity, RC3, RC15}). We construct these parameters by first sampling $\widetilde{A}$, $\widetilde{C}$ and $\widetilde{\boldsymbol{\zeta}}$ from appropriately chosen laws. Then, we normalize each element of the tuple such that
\begin{equation}
\label{normalizing}
\overline{A} = {\widetilde{A}}/{\rho(\widetilde{A})}
, \quad
\overline{C} = {\widetilde{C}}/{\lVert\widetilde{C}\rVert}
, \quad
\overline{\boldsymbol{\zeta}} = {\widetilde{\boldsymbol{\zeta}}}/{\lVert\widetilde{\boldsymbol{\zeta}}\rVert}
,
\end{equation}
where $\rho(\widetilde{A})$ denotes the spectral radius of $\widetilde{A}$. As discussed in the properties of reservoir systems in Subsection~\ref{subsec:reservoir_models}, the sufficient condition of the ESP is $\|A\|_2<1$. By this normalizing choice, we allow for some more flexibility in terms of marginal violations of the non-sharp ESP constraint. Finally, defining $A= \rho \overline{A}$, $C=\gamma \overline{C}$, and $\boldsymbol{\zeta}=\omega \overline{\boldsymbol{\zeta}}$, we can rewrite state equation \eqref{eq:esn_model_states} as
\begin{equation}\label{eq:esn_hyperparams}
\boldsymbol{x}_t = \alpha \boldsymbol{x}_{t-1} + (1-\alpha) \sigma(\rho \overline{A} \boldsymbol{x}_{t-1} + \gamma \overline{C} \boldsymbol{z}_t + \omega \overline{\boldsymbol{\zeta}}) .
\end{equation}
We refer to tuple $\boldsymbol{\varphi}:=(\alpha, \rho, \gamma, \omega)$ as the \textit{hyperparameters} of the ESN. Specifically, $\alpha \in [0,1)$ is the leak rate and $\rho \in \mathbb{R}_+$ is called the \textit{spectral radius} of the reservoir matrix, $\gamma \in \mathbb{R}_+$ is the \textit{input scaling}, and $\omega \in \mathbb{R}_+$ is the \textit{shift scaling}.
The choice of the hyperparameters determines the properties of the state map. For simplicity, in Section~\ref{sec:empirical_study}, we choose the hyperparameters based on the empirical ESN literature. In Supplementary Appendix \ref{Hyperparameter Tuning}, we also propose a general though more computationally intensive procedure to select hyperparameters in a data-driven way that could be interesting to practitioners.
\subsubsection{Penalty Selection}
To apply ridge estimator \eqref{eq:ridge_regression}, it is necessary to first select a penalty $\lambda$. Cross-validation (CV) is a common selection procedure for regularization strength in penalized methods such as ridge, LASSO, and Elastic Net. CV techniques have also been applied in the time series context \citep{kockPenalizedTimeSeries2020,Ballarin2023a} with their validity established in \cite{bergmeir2018}.
In our empirical study, to account for temporal dependence, we use a sequential CV strategy with ten validation folds.
More precisely, we reserve the last 50 observations for validation and all other previous data points for training. The first fold consists of the first five observations out of the validation set, and the model is fitted using all training data. The following validation fold comprises the next five subsequent validation observations while the training set is expanded by five data points (from the previous fold). This procedure is repeated ten times and the CV loss is the average of the one-step-ahead forecast MSE on each fold.
In expanding or rolling window setups, we rerun the CV penalty selection to ensure that estimated ESN coefficients do not induce oversmoothing. We refer the reader to Supplementary Appendix~\ref{sec:cross_validation} for additional details.
\subsection{Relation to Nonparametric Regression}
Together with hyperparameters and penalty strength selection, the choice of the state dimension $N$ is a key ingredient of an ESN model. A large state space generally implies better approximation bounds \citep{RC12, RC24}. Although it is customary in the empirical literature to take $N$ as large as possible \citep{Lukosevicius2012}, some recent literature discusses both the statistical risk bounds and the approximation-risk trade-off bounds for various RC families (see \cite{RC10} and \cite{RC24} for details). Under simplified assumptions that $\alpha = 0$ and $\rho = 0$ in \eqref{eq:esn_hyperparams}, ESNs have a natural connection to random-weights neural networks \citep{caoReviewNeuralNetworks2018} and random projection regression \citep{maillardLinearRegressionRandom2012}, and are thus comparable to nonparametric sieve methods. If the data were independently sampled, known results on sieve estimation would require that at most $N/T = o(1)$ up to logarithmic factors for consistency \citep{belloniNewAsymptoticTheory2015a}. \cite{chenOptimalUniformConvergence2015} have extended this result to $\beta$-mixing data with B-spline and wavelet sieves.
Sieve rates appear to suggest that choosing $N = O(T)$ in echo state networks could lead to nontrivial forecasting bias owing to poor approximation properties. Unfortunately, this comparison relies on neglecting the dynamic component of the ESN model, and as such it is only qualitative. It is, therefore, an important topic for future research.
A different but related problem is the potential degradation of forecasting performance when a model is at the interpolation threshold in the overparametrized regime, $N \geq T$. Ridge regression is also commonly applied to address generalization concerns in statistical learning (see \cite{hastie2009elements}). Recent work has studied more in-depth the link between regularization and generalization: \cite{hastieSurprisesHighdimensionalRidgeless2022} show that ``ridgeless", that is interpolation, solutions can be optimal in some scenarios. However, in our empirical evaluations in Section~\ref{sec:empirical_study}, cross-validation consistently selects non-zero ridge penalties, confirming that ridge penalization plays an important role in ESN forecasting performance.
\subsection{ESN Forecasting}\label{subsec:esn_forecasting}
We are primarily interested in using ESN models to construct conditional forecasts of target variables. Given that the conditional mean is the best mean square error estimator for $h$-step-ahead target $\boldsymbol{y}_{t+h}$, $h \geq 1$, our main focus is approximating
\begin{equation*}
\widehat{\boldsymbol{y}}_{t+h|t} := \mathbb{E}\left[\boldsymbol{y}_{t+h} | \boldsymbol{x}_{0:t}, \boldsymbol{z}_{0:t}\right].
\end{equation*}
The case $h=1$ is trivial, since the ESN model is estimated by regressing $\boldsymbol{y}_{t+1}$ on state $\boldsymbol{x}_t$, and thus we can set $\widetilde{\boldsymbol{y}}_{t+1|t} = \widehat{W}^\top \boldsymbol{x}_t $. However, when $h > 1$ the nonlinear state dynamics precludes a direct computation of the conditional mean. This is in contrast to linear models like VARMAs or DFMs, where the assumption of linearity implies that conditional expectations reduce to simple matrix-vector operations. In particular, linear models are such that the variance (and any other higher-order moments) of the noise term do not impact the conditional mean forecast.
Let $p_{{\theta}}(\boldsymbol{x}_t \vert \boldsymbol{x}_{t-1}, \boldsymbol{z}_t)$ and $g_{\theta}(\boldsymbol{y}_{t+1} \vert \boldsymbol{x}_t)$ be the state transition and observation densities, respectively. Then, for $h > 1$,
\begin{equation}\label{eq:esn_condexp_integral}
\widehat{\boldsymbol{y}}_{t+h|t} =
\int \boldsymbol{y}_{t+h} \; g_{{\theta}}(\boldsymbol{y}_{t+h} \vert \boldsymbol{x}_{t+h-1}) \prod_{j=1}^{h-1} p_{{\theta}}(\boldsymbol{x}_{t+j} \vert \boldsymbol{x}_{t+j-1}, \boldsymbol{z}_{t+j})
\nu(\boldsymbol{z}_{t+j} \vert \boldsymbol{x}_{t+j-1}) \mathrm{d} \boldsymbol{z}_{t+j} \mathrm{d} \boldsymbol{x}_{t+j} \mathrm{d} \boldsymbol{y}_{t+h}
,
\end{equation}
where $\nu(\boldsymbol{z}_{t+j} \vert \boldsymbol{x}_{t+j-1})$ is the conditional density of inputs. Here, we introduce the additional assumption that $\boldsymbol{x}_{t+j-1}$ is sufficient to condition on past states and inputs, that is
\begin{equation}
\nu(\boldsymbol{z}_{t+j} \vert \boldsymbol{x}_{t+j-1}) \equiv \nu(\boldsymbol{z}_{t+j} \vert \boldsymbol{x}_{0:t+j-1}, \boldsymbol{z}_{0:t+j-1}).
\label{eq:state_conditioning}
\end{equation}
Some elements in the expectation integral are not directly available. Specifically, while an ESN explicitly models both $p_{{\theta}}(\boldsymbol{x}_t \vert \boldsymbol{x}_{t-1}, \boldsymbol{z}_t)$ and $g_{{\theta}}(\boldsymbol{y}_{t+1} \vert \boldsymbol{x}_t)$, the density $\nu(\boldsymbol{z}_{t+j} \vert \boldsymbol{x}_{t+j-1})$ is unavailable.
In the remaining part of this subsection, we present a novel ESN-based approach to forecasting the target variable. Our idea is to enrich the ESN model with an auxiliary observation equation for the input covariates. As we demonstrate in Section~\ref{sec:empirical_study}, our proposed method shows superior performance with respect to the standard state-of-the-art benchmarks.
\subsubsection{Multi-step Forecasting of Targets via Iterative Forecasting of Inputs}
In general, we are interested in constructing forecasts of target variables that are not the same as the model inputs. To do so, we resolve the issue of the intractability of \eqref{eq:esn_condexp_integral} while simultaneously capitalizing on the available results using ESNs in the forecasting of dynamical systems. More explicitly, we add to the ESN specification \eqref{eq:esn_model_states}-\eqref{eq:esn_model_output} an equation that allows sidestepping modeling the density $\nu$ directly, thus making the computation of $\widehat{\boldsymbol{y}}_{t+h|t}$ feasible even when $h > 1$.
Consider the ESN where the reservoir states $(\boldsymbol{x}_t)_{t \in \mathbb{Z}}$ follow \eqref{eq:esn_model_states}, while the target sequence is the same as the input sequence $(\boldsymbol{z}_t)_{t \in \mathbb{Z}}$,
\begin{align}
\boldsymbol{x}_t & = \alpha \boldsymbol{x}_{t-1} + (1-\alpha) \sigma(A \boldsymbol{x}_{t-1} + C \boldsymbol{z}_t + \boldsymbol{\zeta}) \label{eq:autonomous_esn_states} \\
\boldsymbol{z}_{t+1} & = \mathcal{W}^\top \boldsymbol{x}_t + \boldsymbol{u}_{t+1}. \label{eq:autonomous_esn_obs}
\end{align}
Here, we use symbol $\mathcal{W}$ for the output coefficients to separate this case from the general ESN equations \eqref{eq:esn_model_states}-\eqref{eq:esn_model_output}. In \eqref{eq:autonomous_esn_obs}, $(\boldsymbol{u}_t)_{t \in \mathbb{Z}}$ are $K$-dimensional independent zero-mean innovations with variance $\sigma_u^2\mathbb{I}_K$ that are also independent of $\boldsymbol{x}_t$ across all $t$.
In this case, the reservoir map $F(\boldsymbol{x}_{t-1}, \boldsymbol{z}_{t})$ in \eqref{state_map} is determined by \eqref{eq:autonomous_esn_states}, and it is possible to re-feed the forecasted variables back into the state equation as inputs. This yields the following state recursion:
\begin{equation*}
\boldsymbol{x}_t = F(\boldsymbol{x}_{t-1}, \mathcal{W}^\top \boldsymbol{x}_{t-1} + \boldsymbol{u}_{t}) =: G_{\theta}(\boldsymbol{x}_{t-1}, \boldsymbol{u}_{t}),
\end{equation*}
where the subscript $\theta$ denotes the dependence on the model coefficients. In the reservoir computing literature, regimes, where the ESN state equation is iteratively fed with the model outputs, are called ``autonomous'' \citep{RC15}. They are widely and successfully utilized for the prediction of deterministic dynamical systems. Indeed, in those instances, provided that the ridge estimate $\widehat{\mathcal{W}}$ is available from data according to Subsection~\ref{Estimation}, the $h > 1$ steps autonomous state iteration is given by
\begin{equation*}
{F}^{*}_{\theta}(\boldsymbol{x}_{t}) := \alpha \boldsymbol{x}_{t} + (1-\alpha) \sigma((A + C \widehat{\mathcal{W}}^\top) \boldsymbol{x}_t + \boldsymbol{\zeta})
\end{equation*}
and
\begin{equation*}
\boldsymbol{x}_{t+h} = \underbrace{{F}^{*}_{\theta} \circ {F}^{*}_{\theta} \circ \cdots \circ {F}^{*}_{\theta} }_{h \textnormal{ times}} \; (\boldsymbol{x}_t ).
\end{equation*}
Hence one can directly obtain the $h$-steps ahead predictions of the input time series as $\boldsymbol{z}_{t+h} = \widehat{\mathcal{W}}^\top \boldsymbol{x}_{t+h-1}$.
In the case of stochastic target variables, assuming \eqref{eq:state_conditioning}, we notice that for the conditional forecast of the states, it holds that
\begin{equation}\label{eq:esn_condexp_nu_int}
\widehat{\boldsymbol{x}}_{t+1|t} = \mathbb{E}\left[\boldsymbol{x}_{t+1} | \boldsymbol{x}_{0:t}, \boldsymbol{z}_{0:t}\right] = \int \boldsymbol{x}_t \, p_{{\theta}}(\boldsymbol{x}_t \vert \boldsymbol{x}_{t-1}, \boldsymbol{z}_t) \nu(\boldsymbol{z}_t \vert \boldsymbol{x}_{t-1}) \mathrm{d} \boldsymbol{z}_t
=
\int G_{{\theta}}(\boldsymbol{x}_{t-1}, \boldsymbol{u}_t) \phi(\boldsymbol{u}_t) \mathrm{d} \boldsymbol{u}_t,
\end{equation}
where density $\phi$ of $\boldsymbol{u}_t$ is, again, unavailable. Note that, even under the assumption $\boldsymbol{u}_t \sim \mathcal{N}(\boldsymbol{0}, \Sigma_{\boldsymbol{u}})$, which is standard in the filtering literature, the presence of nonlinear map $G_{{\theta}}$ makes the computation of the forecasts of $\boldsymbol{z}_{t+h}$ a non-straightforward exercise. Nevertheless, this forecast construction can be readily used when one is interested exclusively in predicting the time series $\boldsymbol{z}_t$.
Whenever the final goal of the exercise is forecasting some other explained variable $\boldsymbol{y}_{t+h}$ $h$-steps ahead, additional issues arise. In this case, one needs to compute the conditional expectation in \eqref{eq:esn_condexp_integral} which is intractable even under Gaussian assumptions on the innovations. One option is to apply particle filtering techniques such as bootstrap sampling or sequential importance sampling (SIS) to evaluate the expectation \citep{doucet2001sequential}. We emphasize that the state dimension is usually chosen to be large, and hence implementing filtering techniques requires some care.
Our approach is to avoid dealing with the nonlinear densities involved in \eqref{eq:esn_condexp_integral} with the help of \eqref{eq:esn_condexp_nu_int} and, instead, to reduce the computation of the conditional expectation $\widehat{\boldsymbol{y}}_{t+h|t}$ to a composition of functions. By the linearity of observation equation \eqref{eq:esn_model_output} and the assumption of independence in the zero-mean noise $\boldsymbol{\epsilon}_{t+h}$, we write
\begin{align}
\widehat{\boldsymbol{y}}_{t+h|t}
& = W^\top \widehat{\boldsymbol{x}}_{t+h-1|t} =
\bigintssss W^\top \boldsymbol{x}_{t+h-1} \prod_{j=1}^{h-1} p_{{\theta}}(\boldsymbol{x}_{t+j} \vert \boldsymbol{x}_{t+j-1}, \boldsymbol{z}_{t+j})
\nu(\boldsymbol{z}_{t+j} \vert \boldsymbol{x}_{t+j-1}) \mathrm{d} \boldsymbol{x}_{t+j} \mathrm{d} \boldsymbol{z}_{t+j} \nonumber
\end{align}
and use the approximation
\begin{align}
\widehat{\boldsymbol{y}}_{t+h|t} \approx \widetilde{\boldsymbol{y}}_{t+h} = W^\top \underbrace{{F}^{*}_{{\theta}} \circ {F}^{*}_{{\theta}} \circ \cdots \circ {F}^{*}_{{\theta}} }_{h-1 \textnormal{ times}} \; (\boldsymbol{x}_{t}), \label{eq:esn_forecast_integral_approx}
\end{align}
which originates from
\begin{equation}\label{eq:esn_condexp_state_auto_approx}
\widehat{\boldsymbol{x}}_{t|t-1} = \int G_{{\theta}}(\boldsymbol{x}_{t-1}, \boldsymbol{u}_t) \phi(\boldsymbol{u}_t) \mathrm{d} \boldsymbol{u}_t
\approx
G_{{\theta}}(\boldsymbol{x}_{t-1}, \mathbb{E}[\boldsymbol{u}_t])
=
F( \boldsymbol{x}_{t-1}, \mathcal{W}^\top \boldsymbol{x}_{t-1} )
\equiv
{F}^{*}_{{\theta}}( \boldsymbol{x}_{t-1} ),
\end{equation}
where $\boldsymbol{u}_t$ is assumed to be zero-mean. The validity of \eqref{eq:esn_condexp_state_auto_approx} itself requires implicit assumptions on the nature of the distribution of $\boldsymbol{u}_t$, but here we want to keep the analysis of $\widehat{\boldsymbol{y}}_{t+h|t}$ to a minimum, and just use the insights from the dynamical systems ESN literature. We are hence not delving deeper into alternative approaches to estimate forecasts or, more generally, to compute conditional expectations of ESN models with stochastic inputs.
\section{Multi-Frequency Echo State Models}
\label{Multi-Frequency Echo State Models}
In this subsection, we construct a broad class of ESN models that can accommodate input and target time series sampled at distinct sampling frequencies. We call this family of reservoir models the \textit{Multi-Frequency Echo State Networks} (MFESNs). The state-space structure of MFESNs is naturally amenable to the setting of time series with mixed frequencies. Additionally, the prediction strategy discussed in Section~\ref{subsec:esn_forecasting} is straightforward to extend to MFESNs.
We present two groups of MFESN architectures. The first family is based on a single echo state network architecture and we call these models \textit{Single-Reservoir Multi-Frequency Echo State Networks} (S-MFESNs). The second group, referred to as \textit{Multi-Reservoir Multi-Frequency Echo State Networks} (M-MFESNs), allows for as many state equations as the number of distinct sampling frequencies present in the input data.
\subsection{Single-Reservoir MFESN}\label{subsec:Single-reservoir MFESN}
Recall that, in the temporal notation of Definition~\ref{Temporal notation}, we reserve $t$ to be the reference time index, which is also used for the target variable, and all other frequencies will be measured with respect to the reference frequency.
Consider $L$ collections of different time series. We assume that the $l$th collection, $l\in [L]$, consists of $n_l$ time series that are sampled at a common frequency $\kappa_l$ and contain observations $(\boldsymbol{z}^{(l)}_{t,{s}\vert{\kappa_l}})_{t,s}$ with $\boldsymbol{z}^{(l)}_{t,{s}\vert{\kappa_l}} \in \mathbb{R}^{n_l}$ for all $t \in \mathbb{Z}$ and $s \in \{0, \ldots, \kappa_l-1\}$. Let $\kappa_{\max} = \max_l \kappa_l$ be the highest sampling frequency among the $L$ time series groups and let $q_l := \kappa_{\max} / \kappa_l$ indicate how low each $\kappa_l$ sampling frequency is with respect to $\kappa_{\max}$. We can now stack together and repeat the observations in a way that is consistent with the high-frequency index by defining
\begin{equation*}
\boldsymbol{z}_{t,{s}\vert{\kappa_{\max}}} :=
\left( \boldsymbol{z}^{(1)\,\top}_{t,{\lfloor s / q_1 \rfloor}\vert{\kappa_1}} ,
\: \boldsymbol{z}^{(2)\,\top}_{t,{\lfloor s / q_2 \rfloor}\vert{\kappa_2}} ,
\: \ldots ,
\: \boldsymbol{z}^{(L)\,\top}_{t,{ \lfloor s / q_L \rfloor}\vert{\kappa_L}} \right)^\top\in \mathbb{R}^{\sum_{l=1}^{L}n_l}, \enspace s\in \{0, \ldots, \kappa_{\max}-1\},
\end{equation*}
where for all $l\in [L]$, $\boldsymbol{z}^{(l)}_{0,{0}\vert{\kappa_l}} = \mathbf{0}_{n_l}$.
Thus, it is possible to write a single high-frequency ESN as
\begin{align}
\boldsymbol{x}_{t,{s}\vert{\kappa_{\max}}} &= \alpha \boldsymbol{x}_{t,{s-1}\vert{\kappa_{\max}}} + (1-\alpha) \sigma(A \boldsymbol{x}_{t,{s-1}\vert{\kappa_{\max}}} + C \boldsymbol{z}_{t,{s}\vert{\kappa_{\max}}} + \boldsymbol{\zeta}),\label{single ESN state}\\
\boldsymbol{z}_{t,{s+1}\vert{\kappa_{\max}}} & = \mathcal{W}^\top \boldsymbol{x}_{t,{s}\vert{\kappa_{\max}}} + \boldsymbol{u}_{t,{s+1}\vert{\kappa_{\max}}},\label{single ESN input obs}
\end{align}
where $\mathcal{W}\in \mathbb{M}_{N, \sum_{l=1}^{L}n_l}$ and $s >0$.
We term this class of MFESN models the \textit{Single-Reservoir Multi-Frequency ESNs} (S-MFESNs).
Notice that equations \eqref{single ESN state}-\eqref{single ESN input obs} of the S-MFESN model prescribe the dynamics at the highest frequency, $\kappa_{\max}$.
In order to forecast a lower frequency target, we map high-frequency states $\boldsymbol{x}_{t,{s}\vert{\kappa_{\max}}}$ to low-frequency targets $\boldsymbol{y}_{t+1} \in \mathbb{R}^J$ by introducing a \textit{state alignment} scheme.
An \textit{aligned} S-MFESN uses the most recent state with respect to the reference time index $t$ to construct the forecast. More precisely, the state equation of an S-MFESN is iterated $\kappa_{\max}$ times until the state $\boldsymbol{x}_{t-1,{\kappa_{\max}}\vert{\kappa_{\max}}} = \boldsymbol{x}_{t,{0}\vert{\kappa_{\max}}}$ is obtained and then target $\boldsymbol{y}_{t+1}$ is forecast with observation equation
\begin{equation}
\label{single ESN obs}
\boldsymbol{y}_{t+1} = W^\top \boldsymbol{x}_{t,{0}\vert{\kappa_{\max}}} + \boldsymbol{\epsilon}_{t+1}, \quad W\in \mathbb{M}_{N,J}.
\end{equation}
\paragraph{Estimation of aligned S-MFESN.}
Both coefficient matrices ${W}$ and ${\mathcal{W}}$ can be estimated as explained in Subsection~\ref{Estimation} under appropriate choices of corresponding penalty strengths. In particular, in order to obtain $\widehat{\mathcal{W}}$, the state and the observation matrices in \eqref{eq:ridge_regression} are given by
\begin{align*}
X_{\kappa_{\max}} &= (\boldsymbol{x}_{1,{0}\vert{\kappa_{\max}}}, \ldots, \boldsymbol{x}_{1,{\kappa_{\max}-1}\vert{\kappa_{\max}}}, \ldots, \boldsymbol{x}_{T-1,{0}\vert{\kappa_{\max}}}, \ldots, \boldsymbol{x}_{T-1,{\kappa_{\max}-1}\vert{\kappa_{\max}}})^\top \in \mathbb{M}_{(T-1)\kappa_{\max}-1,N},\\
Y_{\kappa_{\max}} &= (\boldsymbol{z}_{1,{1}\vert{\kappa_{\max}}}, \ldots, \boldsymbol{z}_{1,{\kappa_{\max}}\vert{\kappa_{\max}}}, \ldots, \boldsymbol{z}_{T-1,{1}\vert{\kappa_{\max}}}, \ldots, \boldsymbol{z}_{T-1,{\kappa_{\max}}\vert{\kappa_{\max}}})^\top \in \mathbb{M}_{(T-1)\kappa_{\max}-1,\sum_{l=1}^{L}n_l},
\end{align*}
while
\begin{align*}
X &= \left(\boldsymbol{x}_{1,{0}\vert{\kappa_{\max}}}, \boldsymbol{x}_{2,{0}\vert{\kappa_{\max}}}, \ldots, \boldsymbol{x}_{T-1,{0}\vert{\kappa_{\max}}}\right)^\top \in \mathbb{M}_{T-1,N},\\
Y &= \left(\boldsymbol{y}_{2}, \ldots, \boldsymbol{y}_{T}\right)^\top \in \mathbb{M}_{T-1,J},
\end{align*}
are used for the estimation of $\widehat{W}$. We note that the state equation \eqref{single ESN state} of S-MFESN can be initialized by $\boldsymbol{x}_{0,{0}\vert{\kappa_{\max}}}$, which under the fading memory property is inconsequential for long enough samples (see the discussion in Subsection~\ref{Reservoir models}).
\paragraph{Forecasting with aligned S-MFESN.}
Let $\widehat{W}$ and $\widehat{\mathcal{W}}$ be the sample estimates of the readout matrices as explained above. The fitted high-frequency autonomous state transition map associated with \eqref{single ESN state} is given by
\begin{equation}
\label{single ESN state map}
F_{\kappa_{\max}}(\boldsymbol{x}_{t,{s-1}\vert{\kappa_{\max}}}) := \alpha \boldsymbol{x}_{t,{s-1}\vert{\kappa_{\max}}} + (1-\alpha) \sigma\left( (A + C \widehat{\mathcal{W}}^\top) \boldsymbol{x}_{t,{s-1}\vert{\kappa_{\max}}} + \boldsymbol{\zeta} \right),
\end{equation}
which, composed with itself exactly $\kappa_{\max}$ times, yields the target-frequency-aligned autonomous state transition map
\begin{equation}
\label{single ESN state aligned}
F(\boldsymbol{x}_{t,{0}\vert{\kappa_{\max}}}) := \, \underbrace{F_{\kappa_{\max}} \circ F_{\kappa_{\max}} \circ \cdots \circ F_{\kappa_{\max}}}_{\kappa_{\max} \textnormal{ times}}\, ( \boldsymbol{x}_{t, {0}\vert{\kappa_{\max}}} ).
\end{equation}
Finally, from \eqref{eq:esn_forecast_integral_approx} the $h$-steps ahead low-frequency forecasts, $h\in \mathbb{N}$, can be computed as
\begin{align}
\label{single ESN forecast}
\widetilde{y}_{T+h \vert T} & = \widehat{W}^\top \big( \underbrace{F_{ } \circ F \circ \cdots \circ F}_{h-1 \textnormal{ times}}\, ( \boldsymbol{x}_{T, {0}\vert{\kappa_{\max}}} ) \big).
\end{align}
Figure~\ref{fig:single_reservoir_sfesn_diagram} gives a graphical diagram of the 1-step forecasting procedure for an S-MFESN. Additionally, Figure \ref{fig:single_reservoir_sfesn_multistep_diagram} in Supplementary Appendix \ref{subsec:additional_figures} provides a similar diagram for the case of multistep forecasts.
\begin{figure}[!t]
\centering
\includegraphics[width=0.97\linewidth]{figures/smfesn_1step_for_diagram}
\caption{Scheme of a Single-Reservoir MFESN (S-MFESN) model combining input data sampled at two frequencies with state alignment and estimation for one-step ahead forecasting of the target series.}
\label{fig:single_reservoir_sfesn_diagram}
\end{figure}
The following example illustrates this proposed forecasting strategy for the case of quarterly GDP forecasting using monthly and daily series inputs.
\begin{example}
\normalfont
\label{ex:s-mfesn}
Suppose that we wish to use an aligned {S-MFESN} model to forecast a quarterly one-dimensional target $(y_t)$ using $n_{(\textnormal{\texttt{m}})}$ monthly and $n_{(\textnormal{\texttt{d}})}$ daily series, $(\boldsymbol{z}^{(\textnormal{\texttt{m}})}_{t,{s}\vert{\kappa_1}})$ and $(\boldsymbol{z}^{(\textnormal{\texttt{d}})}_{t,{s}\vert{\kappa_2}})$, respectively.
We adopt the assumption that daily data is released 24 days over each calendar month and hence $\kappa_1 = 3$, $\kappa_2 = 72$ and $\kappa_{\max} = 72$, while $q_1 = 24$ and $q_2=1$. Let $t, {\ast}\vert{72}$ be the temporal index with a quarterly reference frequency. The input vector for the S-MFESN state equation consistent with the daily frequency is given by
\begin{equation*}
\boldsymbol{z}^{(\textnormal{\texttt{m,d}})}_{t, {s}\vert{72}} := (
{\boldsymbol{z}^{(\textnormal{\texttt{m}})}_{t,{\lfloor s / 24 \rfloor}\vert{3}}}^\top,\:
{\boldsymbol{z}^{(\textnormal{\texttt{d}})}_{t, {s}\vert{72}}}^\top
)^\top\in \mathbb{R}^{n_{(\textnormal{\texttt{m}})}+n_{(\textnormal{\texttt{d}})}} \enspace {\rm with}\enspace \boldsymbol{z}^{(\textnormal{\texttt{d}})}_{0,{0}\vert{3}} = \mathbf{0}_{n_{(\textnormal{\texttt{d}})}}\enspace {\rm and}\enspace\boldsymbol{z}^{(\textnormal{\texttt{m}})}_{0,{0}\vert{24}} = \mathbf{0}_{n_{(\textnormal{\texttt{m}})}}.
\end{equation*}
The complete S-MFESN model with the state space dimension $N$ can be written as:
\begin{align}
\boldsymbol{x}^{(\textnormal{\texttt{m,d}})}_{t,{s}\vert{72}} & = \alpha \boldsymbol{x}^{(\textnormal{\texttt{m,d}})}_{t,{s-1}\vert{72}} + (1-\alpha) \sigma(A \boldsymbol{x}^{(\textnormal{\texttt{m,d}})}_{t,{s-1}\vert{72}} + C \boldsymbol{z}^{(\textnormal{\texttt{m,d}})}_{t,{s}\vert{72}} + \boldsymbol{\zeta}), \label{eq:s-mfesn_state} \\
\boldsymbol{z}^{(\textnormal{\texttt{m,d}})}_{t,{s+1}\vert{72}} & = \mathcal{W}^\top \boldsymbol{x}^{(\textnormal{\texttt{m,d}})}_{t,{s}\vert{72}} + \boldsymbol{u}_{t,{s+1}\vert{72}}, \label{eq:s-mfesn_input_obs} \\
y_{t+1} & = W^\top \boldsymbol{x}^{(\textnormal{\texttt{m,d}})}_{t,{0}\vert{\kappa_{\max}}} + \boldsymbol{\epsilon}_{t+1}, \label{eq:s-mfesn_target_obs}
\end{align}
where the state equations \eqref{eq:s-mfesn_state}-\eqref{eq:s-mfesn_input_obs} are run in their own maximum frequency temporal index $s>0$, and only the states $\boldsymbol{x}_{t-1,{\kappa_{\max}}\vert{\kappa_{\max}}} = \boldsymbol{x}_{t,{0}\vert{\kappa_{\max}}}$ are used in the observation equation \eqref{eq:s-mfesn_target_obs}.
Provided the input-target pairs sample of length $T$, the coefficient matrices $\mathcal{W}\in \mathbb{M}_{N, n_{(\textnormal{\texttt{m}})}+n_{(\textnormal{\texttt{d}})}}$ in \eqref{eq:s-mfesn_input_obs} and $W\in \mathbb{R}^{N}$ in \eqref{eq:s-mfesn_target_obs} can be estimated via ridge regression as explained above.
From \eqref{single ESN state map} the high-frequency autonomous state transition map is given by
\begin{equation*}
F_{72}^{(\textnormal{\texttt{m,d}})}(\boldsymbol{x}^{(\textnormal{\texttt{m,d}})}_{t,{s-1}\vert{72}}) := \alpha \boldsymbol{x}^{(\textnormal{\texttt{m,d}})}_{t,{s-1}\vert{72}} + (1-\alpha) \sigma\left( (A + C \widehat{\mathcal{W}}^\top) \boldsymbol{x}^{(\textnormal{\texttt{m,d}})}_{t,{s-1}\vert{72}} + \boldsymbol{\zeta} \right),
\end{equation*}
which, composed with itself exactly $72$ times, by \eqref{single ESN state aligned} yields the target-frequency-aligned autonomous state transition map
\begin{equation*}
F^{(\textnormal{\texttt{m,d}})}(\boldsymbol{x}^{(\textnormal{\texttt{m,d}})}_{t,{0}\vert{72}}) := \, \underbrace{F_{72}^{(\textnormal{\texttt{m,d}})} \circ F_{72}^{(\textnormal{\texttt{m,d}})} \cdots \circ F_{72}^{(\textnormal{\texttt{m,d}})}}_{72 \textnormal{ times}}\, ( \boldsymbol{x}^{(\textnormal{\texttt{m,d}})}_{t, {0}\vert{72}} ).
\end{equation*}
By applying $F^{(\textnormal{\texttt{m,d}})}$ to state $\boldsymbol{x}^{(\textnormal{\texttt{m,d}})}_{t,{0}\vert{72}}$ we iterate the S-MFESN forward in time to provide an estimate for $\boldsymbol{x}^{(\textnormal{\texttt{m,d}})}_{t+1,{0}\vert{72}}$, which can then be linearly projected using $\widehat{W}$ to yield a forecast for $y_{t+2}$. For the target variable, as well as forecasts, we do not use our temporal notation for the sake of compactness and clarity of exposition. Finally, the quarterly forecasts for $h\in \mathbb{N}$ can be computed using \eqref{single ESN forecast} as
\begin{align*}
\widetilde{y}_{T+h \vert T} & = \widehat{W}^\top \big( \underbrace{F^{(\textnormal{\texttt{m,d}})}_{\phantom{.}} \circ F^{(\textnormal{\texttt{m,d}})} \circ \cdots \circ F^{(\textnormal{\texttt{m,d}})}}_{h-1 \textnormal{ times}}\, ( \boldsymbol{x}^{(\textnormal{\texttt{m,d}})}_{T, {0}\vert{72}} ) \big).
\end{align*}
\end{example}
\subsection{Multi-Reservoir MFESN}\label{subsec:Multi-reservoir MFESN}
Constructing an MFESN with a single reservoir is not necessarily the most effective modeling strategy. Having more than one reservoir allows more flexible modeling of state dynamics for different subsets of input variables sampled at common frequencies. For example, suppose quarterly and monthly data are used as regressors. Our presentation is general enough to accommodate other types of partitioning of series into the corresponding reservoir models. We leave it to future research to test other approaches based, for instance, on markets or data types as done in \cite{Huellen2020}.
Assume again $L$ groups of series with input observations $(\boldsymbol{z}^{(l)}_{t,{s}\vert{\kappa_l}})_{t,s}$ with $\boldsymbol{z}^{(l)}_{t,{s}\vert{\kappa_l}} \in \mathbb{R}^{n_l}$, $l\in [L]$, for all $t \in \mathbb{Z}$ and $s \in \{0, \ldots, \kappa_l-1\}$ sampled at common frequencies $\{\kappa_1, \ldots, \kappa_L\}$, respectively.
For each of the $L$ groups of input series we define the corresponding ESN model as
\begin{align}
\boldsymbol{x}^{(l)}_{t,{s}\vert{\kappa_l}} & = \alpha_l \boldsymbol{x}^{(l)}_{t,{s-1}\vert{\kappa_l}} + (1-\alpha_l) \sigma(A_l \boldsymbol{x}^{(l)}_{t,{s-1}\vert{\kappa_l}} + C_l \boldsymbol{z}^{(l)}_{t,{s}\vert{\kappa_l}} + \boldsymbol{\zeta}_l),\label{mult ESN state}\\
\boldsymbol{z}^{(l)}_{t,{s+1}\vert{\kappa_l}} & = \mathcal{W}_{l}^\top \boldsymbol{x}^{(l)}_{t,{s}\vert{\kappa_l}} + \boldsymbol{u}^{(l)}_{t,{s+1}\vert{\kappa_l}}, \label{mult ESN input obs} \quad l \in [L],
\end{align}
with $s>0$, $\mathcal{W}_{l}\in \mathbb{M}_{N_l, n_l}$ with $N_l$ the dimension of the state space. Notice that the time index $s$ is different for each $l$ according to our temporal notation introduced in Definition~\ref{Temporal notation} and each state equation runs at its own frequency $\kappa_l$. The dimensions $\{N_1, N_2, \ldots, N_L\}$ of the state spaces can be chosen for the $L$ reservoir models individually. Additionally, multiple reservoirs have the associated hyperparameter tuples $\{\boldsymbol{\varphi}_1, \ldots, \boldsymbol{\varphi}_L\}$ to be tuned. This requires some care whenever one wants to optimize all hyperparameters jointly. Since there are $L$ reservoir state equations, we call this class of MFESN models \textit{Multi-Reservoir Multi-Frequency ESN} (M-MFESN).
Similar to S-MFESN, all $L$ state equations are iterated each $\kappa_l$ times respectively until the states $\boldsymbol{x}^{(l)}_{t-1,{\kappa_l}\vert{\kappa_l}}=\boldsymbol{x}^{(l)}_{t,{0}\vert{\kappa_l}}$ are obtained. The \textit{aligned} M-MFESN observation equation is given by
\begin{align}
\label{mult ESN obs}
\boldsymbol{y}_{t+1} = W^\top \boldsymbol{x}_{t,L} + \boldsymbol{\epsilon}_{t+1}, \enspace {\rm with} \enspace \boldsymbol{x}_{t,L} = \left(\begin{array}{c}
\boldsymbol{x}^{(1)}_{t,{0}\vert{\kappa_1}} \\
\vdots \\
\boldsymbol{x}^{(L)}_{t,{0}\vert{\kappa_L}}
\end{array}\right) \in \mathbb{R}^{\sum_{l=1}^L N_l}, \enspace W\in \mathbb{M}_{\sum_{l=1}^L N_l, J}.
\end{align}
\paragraph{Estimation of aligned M-MFESN.}
The coefficient matrices ${W}_l$, $l\in [L]$, and ${\mathcal{W}}$ can be estimated similarly to the case of S-MFESN. The state and observation matrices for the estimation of $\widehat{\mathcal{W}}_l$, $l\in [L]$, in \eqref{eq:ridge_regression} are constructed as
\begin{align*}
X^{(l)} &= (\boldsymbol{x}^{(l)}_{1,{0}\vert{\kappa_{l}}}, \ldots, \boldsymbol{x}^{(l)}_{1,{\kappa_{l}-1}\vert{\kappa_{l}}}, \ldots, \boldsymbol{x}^{(l)}_{T-1,{0}\vert{\kappa_{l}}}, \ldots, \boldsymbol{x}^{(l)}_{T-1,{\kappa_{l}-1}\vert{\kappa_{l}}})^\top \in \mathbb{M}_{(T-1)\kappa_{l}-1,N_l},\\
Y^{(l)} &= (\boldsymbol{z}^{(l)}_{1,{1}\vert{\kappa_{l}}}, \ldots, \boldsymbol{z}^{(l)}_{1,{\kappa_{l}}\vert{\kappa_{l}}}, \ldots, \boldsymbol{z}^{(l)}_{T-1,{1}\vert{\kappa_{l}}}, \ldots, \boldsymbol{z}^{(l)}_{T-1,{\kappa_{l}}\vert{\kappa_{l}}})^\top \in \mathbb{M}_{(T-1)\kappa_{l}-1,n_l},
\end{align*}
while with the notation as in \eqref{mult ESN obs}
\begin{align*}
X &= \left(\boldsymbol{x}_{1,L},
\boldsymbol{x}_{2,L},
\ldots,
\boldsymbol{x}_{T-1,L}
\right)^\top \in \mathbb{M}_{T-1,\sum_{l=1}^L N_l},\\
Y &= \left(\boldsymbol{y}_{2}, \ldots, \boldsymbol{y}_{T}\right)^\top \in \mathbb{M}_{T-1,J},
\end{align*}
are used for the estimation of $\widehat{W}$. Again, the state equations \eqref{mult ESN state} of M-MFESN can be started with $\boldsymbol{x}^{(l)}_{0,{0}\vert{\kappa_{l}}}=\mathbf{0}_{N_l}$ (see Subsection~\ref{Reservoir models} for more details).
\paragraph{Forecasting with aligned M-MFESN.}
Let $\widehat{W}$ and $\widehat{\mathcal{W}}_{l}$, $l\in [L]$, be the sample estimates of the readout matrices. For any $l\in[L]$ the $\kappa_l$-frequency autonomous state transition map is given by
\begin{equation}
\label{mult ESN l state}
F^{(l)}_{\kappa_l}(\boldsymbol{x}^{(l)}_{t,{s-1}\vert{\kappa_{l}}}) := \alpha_l \boldsymbol{x}^{(l)}_{t,{s-1}\vert{\kappa_{l}}} + (1-\alpha_l) \sigma\left( (A_l + C_l \widehat{\mathcal{W}}_l^\top) \boldsymbol{x}^{(l)}_{t,{s-1}\vert{\kappa_{l}}} + \boldsymbol{\zeta}_l \right).
\end{equation}
The target-frequency-aligned autonomous state transition map associated with each frequency $l$ is hence defined as
\begin{equation}
\label{mult ESN state aligned}
F^{(l)}(\boldsymbol{x}_{t,{0}\vert{\kappa_{l}}}) := \, \underbrace{F^{(l)}_{\kappa_l} \circ F^{(l)}_{\kappa_l} \circ \cdots \circ F^{(l)}_{\kappa_l}}_{\kappa_{l} \textnormal{ times}}\, ( \boldsymbol{x}^{(l)}_{t, {0}\vert{\kappa_{l}}} ).
\end{equation}
Finally, from \eqref{eq:esn_forecast_integral_approx} the $h$-steps ahead forecasts can be computed as
\begin{align}
\label{mult ESN forecast}
\widetilde{y}_{T+h \vert T} & = \widehat{W}^\top \left(\begin{array}{c}
\underbrace{ F^{(1)} \circ F^{(1)} \circ \cdots \circ F^{(1)}}_{h-1 \textnormal{ times}}\, ( \boldsymbol{x}^{(1)}_{T,{0}\vert{\kappa_1}} ) \\
\vdots \\
\underbrace{ F^{(L)} \circ F^{(L)} \circ \cdots \circ F^{(L)}}_{h-1 \textnormal{ times}}\, ( \boldsymbol{x}^{(L)}_{T,{0}\vert{\kappa_L}} )
\end{array}\right).
\end{align}
In Figure~\ref{fig:multi_reservoir_mfesn_diagram} we provide a diagram for the case of 1-step ahead forecasting with an aligned M-MFESN involving regressors of only two frequencies. Figure~\ref{fig:multi_reservoir_mfesn_multistep_diagram} in Supplementary Appendix \ref{subsec:additional_figures} provides a similar diagram for the case of multistep forecasting.
\begin{figure}[!t]
\centering
\includegraphics[width=0.97\linewidth]{figures/mmfesn_1step_for_diagram}
\caption{Scheme of a Multi-Reservoir MFESN (M-MFESN) model combining input data sampled at two frequencies with state alignment and estimation for one-step ahead forecasting of the target series.}
\label{fig:multi_reservoir_mfesn_diagram}
\end{figure}
\begin{example}
\label{ex:m-mfesn}
\normalfont
Similar to Example \ref{ex:s-mfesn}, we aim to forecast a quarterly target with monthly and daily series, but this time we use an M-MFESN model. We have to define two independent state equations, one for monthly and one for daily series; in the observation equations, two states must be aligned temporally and stacked to form the full set of regressors. The data consists again of quarterly $(y_t)$, $n_{(\textnormal{\texttt{m}})}$ monthly series $(\boldsymbol{z}^{(\textnormal{\texttt{m}})}_{t,{s}\vert{3}})$ and $n_{(\textnormal{\texttt{d}})}$ daily series $(\boldsymbol{z}^{(\textnormal{\texttt{d}})}_{t,{s}\vert{72}})$.
The aligned M-MFESN model with two reservoirs of dimensions $N_{(\textnormal{\texttt{m}})}$ and $N_{(\textnormal{\texttt{d}})}$, respectively, is given by
\begin{align}
\boldsymbol{x}^{(\textnormal{\texttt{m}})}_{t,{s}\vert{3}} & = \alpha_1 \boldsymbol{x}^{(\textnormal{\texttt{m}})}_{t,{s-1}\vert{3}} + (1-\alpha_1) \sigma(A_1 \boldsymbol{x}^{(\textnormal{\texttt{m}})}_{t,{s-1}\vert{3}} + C_1 \boldsymbol{z}^{(\textnormal{\texttt{m}})}_{t,{s}\vert{3}} + \boldsymbol{\zeta}_1) \label{eq:m-mfesn_state_1}, \\
\boldsymbol{z}^{(\textnormal{\texttt{m}})}_{t,{s+1}\vert{3}} & = \mathcal{W}_{(\textnormal{\texttt{m}})}^\top \boldsymbol{x}^{(\textnormal{\texttt{m}})}_{t,{s}\vert{3}} + \boldsymbol{u}^{(\textnormal{\texttt{m}})}_{t,{s+1}\vert{3}} \label{eq:m-mfesn_input_obs_1}, \\[3pt]
\boldsymbol{x}^{(\textnormal{\texttt{d}})}_{t,{s}\vert{72}} & = \alpha_2 \boldsymbol{x}^{(\textnormal{\texttt{d}})}_{t,{s-1}\vert{72}} + (1-\alpha_2) \sigma(A_2 \boldsymbol{x}^{(\textnormal{\texttt{d}})}_{t,{s-1}\vert{72}} + C_2 \boldsymbol{z}^{(\textnormal{\texttt{d}})}_{t,{s}\vert{72}} + \boldsymbol{\zeta}_2) \label{eq:m-mfesn_state_2}, \\
\boldsymbol{z}^{(\textnormal{\texttt{d}})}_{t,{s+1}\vert{72}} & = \mathcal{W}_{(\textnormal{\texttt{d}})}^\top \boldsymbol{x}^{(\textnormal{\texttt{d}})}_{t,{s}\vert{72}} + \boldsymbol{u}^{(\textnormal{\texttt{d}})}_{t,{s+1}\vert{72}} \label{eq:m-mfesn_input_obs_2}, \\[3pt]
y_{t+1} & = W^\top \left(\begin{array}{c}
\boldsymbol{x}^{(\textnormal{\texttt{m}})}_{t,{0}\vert{3}} \\[3pt]
\boldsymbol{x}^{(\textnormal{\texttt{d}})}_{t,{0}\vert{72}}
\end{array}\right) + \boldsymbol{\epsilon}_{t+1}, \label{eq:m-mfesn_target_obs}
\end{align}
where $s>0$, $\mathcal{W}_{(\textnormal{\texttt{m}})} \in \mathbb{M}_{N_{(\textnormal{\texttt{m}})},n_{(\textnormal{\texttt{m}})}}$, $\mathcal{W}_{(\textnormal{\texttt{d}})} \in \mathbb{M}_{N_{(\textnormal{\texttt{d}})},n_{(\textnormal{\texttt{d}})}}$ and $W\in \mathbb{R}^{N_{(\textnormal{\texttt{m}})}+N_{(\textnormal{\texttt{d}})}}$.
Here, the monthly reservoir $(\boldsymbol{x}^{(\textnormal{\texttt{m}})}_{t,{s}\vert{3}})$ has the temporal index of frequency $3$, while the daily reservoir $(\boldsymbol{x}^{(\textnormal{\texttt{d}})}_{t,{s}\vert{72}})$ of $72$; the high-frequency index $s$ is different for the two models. Notice that in an M-MFESN model it is necessary to introduce $2$ additional observation equations for the states, that is \eqref{eq:m-mfesn_input_obs_1} and \eqref{eq:m-mfesn_input_obs_2}). Notice that the state equations are iterated each $\kappa_l$ times to collect the states to be aligned in the observation equation \eqref{eq:m-mfesn_target_obs}. Again, the sample-based estimates of coefficient matrices $\widehat{\mathcal{W}}_{(\textnormal{\texttt{m}})}$, $\widehat{\mathcal{W}}_{(\textnormal{\texttt{d}})}$ and $\widehat{W}$ in \eqref{eq:m-mfesn_input_obs_1}, \eqref{eq:m-mfesn_state_2}, and in \eqref{eq:m-mfesn_target_obs}, respectively, can be obtained via the ridge regression as discussed above.
Exactly as in Example~\ref{ex:s-mfesn}, using \eqref{mult ESN l state} we can introduce high-frequency autonomous state maps $F_{3}^{(\textnormal{\texttt{m}})}$ and $F_{72}^{(\textnormal{\texttt{d}})}$ as
\begin{align*}
F_{3}^{(\textnormal{\texttt{m}})}(\boldsymbol{x}^{(\textnormal{\texttt{m}})}_{t,{s-1}\vert{3}}) & := \alpha_1 \boldsymbol{x}^{(\textnormal{\texttt{m}})}_{t,{s-1}\vert{3}} + (1-\alpha_1) \sigma\left( (A_1 + C_1 \widehat{\mathcal{W}}_{(\textnormal{\texttt{m}})}^\top) \boldsymbol{x}^{(\textnormal{\texttt{m}})}_{t,{s-1}\vert{3}} + \boldsymbol{\zeta}_1 \right), \\
F_{72}^{(\textnormal{\texttt{d}})}(\boldsymbol{x}^{(\textnormal{\texttt{d}})}_{t,{s-1}\vert{72}}) & := \alpha_2 \boldsymbol{x}^{(\textnormal{\texttt{d}})}_{t,{s-1}\vert{72}} + (1-\alpha_2) \sigma\left( (A_2 + C _2\widehat{\mathcal{W}}_{(\textnormal{\texttt{d}})}^\top) \boldsymbol{x}^{(\textnormal{\texttt{d}})}_{t,{s-1}\vert{72}} + \boldsymbol{\zeta}_2 \right),
\end{align*}
as well as their target-frequency aligned counterparts $F^{(\textnormal{\texttt{m}})}$ and $F^{(\textnormal{\texttt{d}})}$, by \eqref{mult ESN state aligned}, as
\begin{align*}
F^{(\textnormal{\texttt{m}})}(\boldsymbol{x}^{(\textnormal{\texttt{m}})}_{t,{0}\vert{3}}) &:= \, \underbrace{F_{3}^{(\textnormal{\texttt{m}})} \circ F_{3}^{(\textnormal{\texttt{m}})} \circ F_{3}^{(\textnormal{\texttt{m}})}}_{3 \textnormal{ times}}\, ( \boldsymbol{x}^{(\textnormal{\texttt{m}})}_{t, {0}\vert{3}} ),\\
F^{(\textnormal{\texttt{d}})}(\boldsymbol{x}^{(\textnormal{\texttt{d}})}_{t,{0}\vert{72}}) &:= \, \underbrace{F_{72}^{(\textnormal{\texttt{d}})} \circ F_{72}^{(\textnormal{\texttt{d}})} \circ \cdots \circ F_{72}^{(\textnormal{\texttt{d}})}}_{72 \textnormal{ times}}\, ( \boldsymbol{x}^{(\textnormal{\texttt{d}})}_{t, {0}\vert{72}} ).
\end{align*}
The $h$-step ahead forecasts can be computed using the approximation in \eqref{mult ESN forecast} as
\begin{align*}
\widetilde{y}_{T+h \vert T} & = \widehat{W}^\top \left(\begin{array}{c}
\underbrace{ F^{(\textnormal{\texttt{m}})}_{\phantom{.}} \circ F^{(\textnormal{\texttt{m}})} \circ \cdots \circ F^{(\textnormal{\texttt{m}})} }_{h-1 \textnormal{ times}}\, ( \boldsymbol{x}^{(\textnormal{\texttt{m}})}_{T,{0}\vert{3}} ) \\[13pt]
\underbrace{ F^{(\textnormal{\texttt{d}})}_{\phantom{.}} \circ F^{(\textnormal{\texttt{d}})} \circ \cdots \circ F^{(\textnormal{\texttt{d}})} }_{h-1 \textnormal{ times}}\, ( \boldsymbol{x}^{(\textnormal{\texttt{d}})}_{T,{0}\vert{72}} )
\end{array}\right).
\end{align*}
In this case, it is important to note that while both $F^{(\textnormal{\texttt{m}})}$ and $F^{(\textnormal{\texttt{d}})}$ are composed $h-1$ times at step $h$, the underlying number of autonomous reservoir iterations is different for the monthly and daily reservoirs, namely $3$ and $72$, and depends on their own frequencies. This also suggests that one should take into account the different time dynamics when, for example, tuning M-MFESN hyperparameters $\boldsymbol{\varphi}^{(\textnormal{\texttt{m}})}$ and $\boldsymbol{\varphi}^{(\textnormal{\texttt{d}})}$, as proposed in Subsection~\ref{Hyperparameter Tuning}.
\end{example}
\section{Empirical Study}
\label{sec:empirical_study}
In this section, we compare the forecasting performance of our proposed MFESN to state-of-the-art benchmarks. We use a combination of macroeconomic and financial data sampled at low and high-frequency intervals, respectively. Our empirical exercises encompass several setups, with a small and a medium-sized set of regressors, fitting models with data before and after the 2007-08 crisis, and with fixed, rolling, and expanding estimation windows.
\subsection{Data}
\label{subsec:data}
Two sets of predictors of different sizes are compiled: Small-MD with 9 predictors and Medium-MD with 33 predictors in monthly and daily frequency. The reference frequency is quarterly: this is the frequency at which the target variable, US GDP growth, is available. Seasonally adjusted quarterly and monthly data is obtained from the Federal Reserve Bank of St. Louis Monthly (FRED-MD) and Quarterly (FRED-QD) Databases for Macroeconomic Research (see \cite{McCracken2016,McCracken2020} for detail). Daily data is obtained from Refinitiv Datastream, a subscription-based data service. All data is the last revised vintage data. The macroeconomic target and predictors, their transformations, and availability are provided in full detail in Table \ref{tab:dataMD} in Appendix~\ref{subsec:data_table}.
The selection of predictors follows the seminal work by \cite{Stock1996,Stock2006} in which the FRED-MD and FRED-QD data are proposed. Variations of their dataset have been used profusely in the literature (for example, see \cite{Boivin2005,Marcellino2006,Hatzius2010}). Indicators from ten macroeconomic and financial categories are considered: (1) output and income, (2) labor market, (3) housing, (4) orders and inventories, (5) price indices, (6) money and credit, (7) interest rates, (8) exchange rates, (9) equity, and (10) derivatives. The latter five categories represent financial market conditions and are sourced at daily frequency. The exception is interest rates, which move relatively slowly and enter as monthly aggregates, available in the FRED-MD data. We refer to this dataset as Medium-MD. A subset of predictors is selected for the Small-MD dataset by choosing variables that have been identified as leading indicators in the empirical literature \citep{Ingenito1996, Clements2008, Andreou2013, marsilli:thesis, clement:midas, Carriero2019, Jardet2020}. Data availability is an additional criterion, and predictors unavailable before 1990 are not considered. This excludes the VIX volatility index, which has been identified as a leading indicator in some studies, for example in \cite{Andreou2013,Jardet2020}.
We follow instructions by \cite{McCracken2016,McCracken2020} on pre-processing macroeconomic predictors before they are used as input for forecasting. These are mainly differenced for detrending. We further transform financial predictors to capture market disequilibrium and volatility. Disequilibrium indicators, such as interest rate spreads, have been found to be more relevant for macroeconomic prediction than routine changes captured by differencing (see \cite{Borio2002,Gramlich2010,Qin2022}). In addition to disequilibrium indicators, realized stock market volatility has been found to improve macroeconomic predictions \citep{Chauvet2015}. In the absence of intraday trading data from the 1990s onward, which prevents us from utilizing conventional daily realized volatility indicators, we extract volatility indicators from daily price series by fitting a GARCH(1,1) by \citet{bollerslev:garch}.\footnote{We include a control \texttt{scale = 1} to ensure convergence of the optimization algorithm and only include a constant mean term in the return process for simplicity.} In addition to volatility of stock and commodity prices, term structure indicators are used. The term structure is forward-looking, capturing information about future demand and supply, and has been found to be a leading predictor of GDP growth (see for example \citet{Hong2012,Kang2020}).
The data spans the period January 1st, 1990 to December 31st, 2019.\footnote{In the Small-MD dataset experiments we make a small variation and instead include data starting from 1st January 1975, but \textit{only} for the initial CV selection of ridge penalties for MFESN models. Our aim is to make sure that at least for the fixed window estimation strategy -- where $\lambda$ is cross-validated once and only one $\widehat{W}$ is estimated -- the ridge estimator is robust. In practice, when we compare to expanding and rolling window estimators, where $\lambda$ is re-selected at each window, we find that extending the initial CV data window has little impact on out-of-sample performance.} We are interested in evaluating model performance under two stylized settings. First, a researcher fits all models up until the Great Recession, including data from Q1 1990 to Q4 2007. Second, fitting is done with data largely encompassing the crisis period, again from Q1 1990 but now up to Q4 2011. In both cases, the testing sample ranges from the next GDP growth observation after fitting up to Q4 2019. All exercises exclude the global COVID-19 economic depression, as we consider it as an extreme, unpredictable event that induces significant structural changes in the underlying macroeconomic dynamics.\footnote{In the macroeconomic literature this falls under the category of ``natural disaster'' events, and should not be na\"ively modeled together with previous observations. In this section, we therefore avoid dealing with post-COVID-19 macroeconomic data altogether.}
To avoid having to handle the many edge cases that daily data in its ``raw'' calendar releases involves, we use an interpolation approach. We set \textit{ex ante} the number of working days in \textit{any} month to be exactly 24: given that in forecasting the most recent information sets are more relevant, when interpolating daily data over months with less than calendar 24 observations, we linearly interpolate the ``missing'' data starting from a months' beginning (using the previous months' last observation). The choice of 24 as a daily frequency is transparent by noting that this is the closest number to actual commonly observed data releases, whilst also being a multiple of both 4 (approximate number of weeks per month) and 6 (upper bound on the number of working days per week).
\subsection{Models}
In this section, we present the set of models that we use throughout our empirical exercises. For a general overview, Table \ref{table:models_empirical} summarizes all models, including hyperparameters. In our analysis, we compare the competing models based on several performance measures, which we introduce in Supplementary Appendix~\ref{Performance measures}.
\begin{table}[!t]
\centering
\setlength\tabcolsep{6pt}
\setlength\extrarowheight{0pt}
\linespread{1}\selectfont\centering
\begin{tabular}{m{0.15\textwidth}<{\centering}m{0.45\textwidth}m{0.32\textwidth}}
\hline
Model Name & Description & Specification \\
\hline
Mean &
Unconditional mean of target series over estimation sample. &
None \\[4pt]
AR(1) &
Autoregressive model of target series estimated using OLS. &
None \\[8pt]
MIDAS &
Almon-weighted MIDAS regression, linear (unconstrained) autoregressive component. &
\shortstack[l]{
Autoregressive lags: 3 \\
Monthly freq. lags: 9 \\
Daily freq. lags: 30
} \\[15pt]
DFM [A] &
Stock aggregation, VAR(1) factor process. &
\shortstack[l]{
Factors: 5 for Small-MD \\
\phantom{Factors:} 10 for Medium-MD
} \\[10pt]
DFM [B] &
Almon aggregation, VAR(1) factor process &
\shortstack[l]{
Factors: 5 for Small-MD \\
\phantom{Factors:} 10 for Medium-MD
} \\[10pt]
singleESN [A] &
\shortstack[l]{
S-MFESN model: \\
$\quad$ Sparse-normal $\widetilde{A}$, sparse-uniform $\widetilde{C}$, $\widetilde{\boldsymbol{\zeta}} = \mathbf{0}$. \\
$\quad$ Isotropic ridge regression fit.\vspace{10pt}
} & \shortstack[l]{
Reservoir dim: 30 \\
Sparsity: $33.3\%$ \\
$\rho = 0.5$, $\gamma = 1$, $\alpha = 0.1$
} \\[20pt]
singleESN [B] &
\shortstack[l]{
S-MFESN model: \\
$\quad$ Sparse-normal $\widetilde{A}$, sparse-uniform $\widetilde{C}$, $\widetilde{\boldsymbol{\zeta}} = \mathbf{0}$. \\
$\quad$ Isotropic ridge regression fit.
} & \shortstack[l]{
Reservoir dim: 120 \\
Sparsity: $8.3\%$ \\
$\rho = 0.5$, $\gamma = 1$, $\alpha = 0.1$
} \\[20pt]
multiESN [A] &
\shortstack[l]{
M-MFESN model: \\
$\quad$ Monthly and daily frequency reservoirs. \\
$\quad$ Sparse-normal $\widetilde{A}_1$, $\widetilde{A}_2$,\\
$\quad$ sparse-uniform $\widetilde{C}_1$, $\widetilde{C}_2$, $\widetilde{\boldsymbol{\zeta}}_1 = \mathbf{0}$, $\widetilde{\boldsymbol{\zeta}}_2 = \mathbf{0}$. \\
$\quad$ Isotropic ridge regression fit.\vspace{10pt}
} & \shortstack[l]{
Reservoir dims: M=100, D=20 \\
Sparsity: M=$10\%$, D=$50\%$ \\
M: $\rho = 0.5$, $\gamma = 1.5$, $\alpha = 0$ \\
D: $\rho = 0.5$, $\gamma = 0.5$, $\alpha = 0.1$
} \\[20pt]
multiESN [B] &
\shortstack[l]{
M-MFESN model: \\
$\quad$ Monthly and daily frequency reservoirs. \\
$\quad$ Sparse-normal $\widetilde{A}_1$, $\widetilde{A}_2$, \\
$\quad$ sparse-uniform $\widetilde{C}_1$, $\widetilde{C}_2$, $\widetilde{\boldsymbol{\zeta}}_1 = \mathbf{0}$, $\widetilde{\boldsymbol{\zeta}}_2 = \mathbf{0}$. \\ $\quad$ Isotropic ridge regression fit.\vspace{10pt}
} & \shortstack[l]{
Reservoir dims: M=100, D=20 \\
Sparsity: M=$10\%$, D=$50\%$ \\
M: $\rho = 0.08$, $\gamma = 0.25$, $\alpha = 0.3$ \\
D: $\rho = 0.01$, $\gamma = 0.01$, $\alpha = 0.99$
} \\[15pt]
\hline
\end{tabular}
\caption{Table of models used in applied forecasting exercises. MFESN hyperparameters are defined with respect to normalized state parameters c.f. \eqref{normalizing}.}
\label{table:models_empirical}
\end{table}
\subsubsection{Benchmarks}
\label{Benchmarks}
\paragraph{Unconditional mean.}
We use the unconditional mean of the sample used for fitting as a baseline benchmark. For GDP growth forecasting, there is evidence that the unconditional mean produces forecasts that are competitive with linear models such as VARs in terms of mean square forecasting errors (MSFE), even at relatively short horizons \citep{aroraNonlinearNonparametricModeling2013}. It is therefore an important reference for the performance of all other models and we report relative MSFE with respect to the unconditional mean in the tables below.
\paragraph{AR(1) model.}
A simple autoregressive process of order one on the target variable is included as a benchmark model.\footnote{Suggested by an anonymous referee.} This is also a common benchmark in the literature, as AR(1) models are often able to capture key dynamics and produce meaningful forecasts for macroeconomic variables \citep{Stock2002,Bai2008}. We emphasize that since AR(1) model is fit to the series of quarterly GDP targets and does not use any additional information, its forecasts are identical for both the Small-MD and Medium-MD samples.
\paragraph{MIxed DAta Sampling (MIDAS).}
The first mixed-frequency model benchmark is given by a MIDAS model \citep{Ghysels2004, Ghysels2007}. Our dynamic MIDAS specification includes autoregressive lags of the target series and uses an Almon weighting scheme. As shown in \cite{Bai2010}, exponential Almon MIDAS regressions are related to dynamic factor models, which we also consider as benchmarks. The MIDAS model includes three lags of quarterly GDP target variable, and 30 daily and 9 monthly lags for all daily and monthly series, respectively. This model prescription allows for some parsimony as the Almon polynomial weighing reduces the number of daily and monthly lag coefficients.
A thorough description of our MIDAS implementation can be found in Supplementary Appendix~\ref{subsec_midas}. To make optimization more efficient, we use explicit expressions for MIDAS loss gradients as in \cite{Kostrov2021essays}. The MIDAS estimation can be hard to perform in practice due to the complexity of nonlinear optimization. First, exponential weighting schemes might require computing floating-point numbers that exceed numerical precision. Therefore, it is a better choice to start the gradient descent close to the origin of the parameter space. Second, even with this choice of starting points, one may encounter issues with optimization results since the Almon-scheme MIDAS loss can have a large number of distinct local minima. In Supplementary Appendix \ref{sec:midas_robustness} we document, using a simple replication experiment, that even small changes in the initial conditions can result in different local minima picked by the numerical optimization algorithm.\footnote{We set the initial coefficient values to zero in all empirical exercises.} These important robustness issues are present even when using closed-form gradients and multi-start optimization routines for the MIDAS models. The computational issues become more pronounced as the number of MIDAS parameters increases unless a careful model/variable selection step is performed. We, therefore, do not include any MIDAS model specifications in the Medium-MD setup.
\paragraph{Dynamic Factor Model (DFM).}
The dynamic factor model framework has been extensively applied in macroeconometrics, starting with \cite{geweke1977dynamic} and \cite{sargent1977business}. A DFM specification assumes that predictable dynamics of a large set of time series can be explained by a small number of factors with an autoregressive dependence (see for example \cite{forni2005generalized, doz2011two, stock2016dynamic}). We generalize the standard two-frequency DFM modeling setup \citep{mariano2003new,banbura2014maximum} to a flexible mixed-frequency DFM that encompasses any number of data frequencies. Moreover, we derive a novel weighting scheme that effectively links the MIDAS and DFM approaches. For a detailed discussion of our factor model setup, we refer the reader to Supplementary Appendix \ref{subsec_dfm}. Two distinct DFM specifications are used. The first one termed DFM [A] uses the standard linear aggregation scheme, as provided in Example \ref{ex:stock}, while the second is a variation that implements an Almon weighting scheme as presented in Example \ref{ex:almon} (we name it DFM [B]). The latter is similar to a MIDAS-type aggregation scheme \citep{Marcellino2010}: the factor structure effectively mitigates the parameter proliferation.
A key choice for a DFM model is the dimension of the factor process. While a number of methods have been developed over the years to systematically derive the number of factors (see, for example, the review of \cite{stock2016dynamic}), commonly used macroeconomic panels feature a number of challenges, such as weak factors \citep{onatskiAsymptoticsPrincipalComponents2012}. Moreover, as mentioned in Supplementary Appendix~\ref{sec:DFM_implementation}, factor number selection in the mixed-frequency setting has not been sufficiently addressed in the literature. To sidestep these issues, we construct both DFM models with 5 unobserved factors for Small-MD and 10 for Medium-MD, respectively, and assume that they follow a VAR(1) process.
One extant issue with integrating daily data is its very high release frequency compared to monthly and especially quarterly releases: computationally this can be extremely taxing, which might be one of the reasons why to our knowledge we are {\it the first to provide DFM forecasts that include daily data}. Our solution is to reduce aggregate daily data every 6 days by averaging, thus leaving 4 observations per month. This eases the computational burden to estimate coefficients and latent states considerably (12 versus 72 daily observations per quarter).
\subsubsection{Multi-Frequency ESNs}
The first set of ESNs we propose is given by two S-MFESN models, based on Example \ref{ex:s-mfesn}. One model uses a reservoir of 30 neurons (we call it singleESN [A]); the other has a larger reservoir of dimension 120 (named singleESN [B]). The sparsity degree of state parameters for both models is set to be $10/N$, where $N$ is the reservoir size. Both MFESNs share the same hyperparameters, $\rho = 0.5$, $\gamma = 1$, $\alpha = 0.1$ (see \eqref{eq:esn_hyperparams}). These values have not been tuned but are presumed credible given other ESN implementations in the literature. To make a fair comparison with DFMs, we fit the S-MFESN models using 6-day-averaged daily data. Note here that for MFESN models the computational gains of averaging are negligible, and are most apparent when tuning the ridge penalty via cross-validation.
Our second set of proposed models consists of two M-MFESNs according to Example \ref{ex:m-mfesn}. Both models have two reservoirs, one for each data frequency -- monthly and daily -- with 100 and 20 neurons, respectively; sparsity degrees are again adjusted to be $10/N$, where $N$ is the reservoir state dimension. The first M-MFESN has hyperparameters that are hand-selected among reasonable values: we note that the monthly-frequency reservoir has no state leak and a larger input scaling, while the daily frequency reservoir features smaller scaling than usual (to avoid compressing high volatility events with the activation function) and the same leak rate as in S-MFESN models (we call this specification multiESN~[A]). For the second M-MFESN, we change hyperparameters more radically: we aim to set up a model that has a very high input memory \citep{RC23}, and that also features long-term smoothing of states. Note that here input scaling values are small, spectral radii are an order of magnitude smaller than in previous models, and leak rates are large (we term this model multiESN~[B]).
\subsection{Results}
\label{subsec:results}
We start by commenting on the computational efficiency of competing models and report execution times in seconds in Table~\ref{tab:execution_times}. Firstly, DFM models appear to be the most computationally effortful models among all specifications. For the Small-MD dataset, the simplest MFESN models, that is, singleESN [A] and [B], have execution times which are at most $3.5$ times higher than the MIDAS model, while still being at least $15.6$ times computationally cheaper than any of the DFM models. The more resource-demanding models MFESN, multiESN [A] and [B], are nevertheless at least $2.6$ times faster to run than the best DFM model (DFM [A]). When moving to the Medium-MD dataset, where the MIDAS model is not, as explained earlier, a feasible choice, the most inefficient MFESN model (singleESN [B]) still outperforms the best DFM model, DFM [A], by $8.4$ times, while the same holds for multiESN [A] model versus DFM [A] model by $2.7$ times. We can conclude that our proposed MFESN architectures provide an attractive and computationally efficient framework for GDP forecasting in the multifrequency framework which is feasible for computations on low-cost machine configurations available to practitioners.
\begin{table}[!t]
\centering
\setlength\extrarowheight{2pt}
\linespread{1.2}\selectfont\centering
Execution Time (Seconds) for Model Estimation \\[5pt]
\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}}*{10}{c}}
\toprule
& & & & \multicolumn{2}{c}{DFM} & \multicolumn{2}{c}{singleESN} & \multicolumn{2}{c}{multiESN} \\
\cline{5-6} \cline{7-8} \cline{9-10}
Dataset & Mean & AR(1) & MIDAS & [A] & [B] & [A] & [B] & [A] & [B] \\
\midrule
Small-MD & 0.1 & 0.7 & 1.3 & $40.5$ & $85.5$ & $2.6$ & $4.5$ & $15.3$ & $14.6$ \\
Medium-MD & 0.1 & 0.8 & $-$ & $48.0$ & $226.5$ & $2.5$ & $5.7$ & $17.7$ & $14.7$ \\
\bottomrule
\end{tabular*}
\caption{\small Execution time in seconds for model estimation measured over a single run on a quad-core computer. MFESN models timing includes ridge penalty cross-validation. MIDAS estimation time refers to optimization from a single initial value. DFM models were estimated on a single-core server and times are adjusted by a factor of $1/4$ for comparison. }
\label{tab:execution_times}
\end{table}
Competing forecasts are compared using the Model Confidence Set (MCS) test derived in \cite{Hansen2011}. One should note that due to the intrinsic nature of data availability of macroeconomic time series and panels, our sample sizes are modest. This implies that the small sample sensitivities of the MCS test need to be taken into account when evaluating our comparisons. Recent analyses of the finite sample properties of the MCS methodology have shown that it requires signal-to-noise ratios which are unattainable in most empirical settings, an issue that undermines its applicability \citep{Aparicio2018}. Given this fact, we also conduct pairwise model comparison tests with the Modified Diebold-Mariano (MDM) test for predictive accuracy \citep{Diebold2002,Harvey1997}.
As we also provide multiple-steps-ahead forecasts, we test for the best subset of models uniformly across all horizons using the Uniform Multi-Horizon MCS (uMCS) test proposed by \cite{quaedvliegMultiHorizonForecastComparison2021}. Since there is relatively little systematic knowledge regarding the power properties of the uMCS test in small samples, our inclusion of this procedure is meant as a statistical counterpoint to simple relative forecasting error comparisons, which provide limited information about the significance of performance differences. We provide more details on our implementation of the test in Supplementary Appendix~\ref{uMCS}.
Finally, we do not report uMCS test outcomes for the expanding window setup, as \cite{quaedvliegMultiHorizonForecastComparison2021} argues that in such context the test is invalid.
\subsubsection{Small Dataset}
We begin our discussion of the Small-MD forecasting results by reviewing Table \ref{table:1sa_GDP_smallMD}. For both sample setups (2007 and 2011) and all three estimation strategies (fixed, expanding, and rolling windows) we provide relative MSFE metrics, with the unconditional mean being used as a reference. Plots of each of the model's forecasts are given in Figures~\ref{fig:1sa_forecasts_2007_GDP_smallMD} and \ref{fig:1sa_forecasts_2011_GDP_smallMD} in Appendix~\ref{subsec:forecast_figures}; additional plots for cumulative SFE, cumulative RMSFEs and other metrics can be found in Supplementary Appendix~\ref{subsec:additional_figures}.
\begin{table}[!t]
\centering
\setlength\tabcolsep{3.3pt}
\setlength\extrarowheight{2pt}
\linespread{1}\selectfont\centering
1-Step-ahead GDP Forecasting - Small-MD Dataset \\[5pt]
\begin{tabular*}{\textwidth}{@{}*{13}{c}}
\toprule
& \multicolumn{4}{c}{Fixed Parameters} & \multicolumn{4}{c}{Expanding Window} & \multicolumn{4}{c}{Rolling Window} \\
\cline{2-5} \cline{6-9} \cline{10-13}
Model & \multicolumn{2}{c}{2007} & \multicolumn{2}{c}{2011} & \multicolumn{2}{c}{2007} & \multicolumn{2}{c}{2011} & \multicolumn{2}{c}{2007} & \multicolumn{2}{c}{2011} \\
\cline{2-3} \cline{4-5} \cline{6-7} \cline{8-9} \cline{10-11} \cline{12-13}
& MSFE & MCS & MSFE & MCS & MSFE & MCS & MSFE & MCS & MSFE & MCS & MSFE & MCS \\
\midrule
Mean & 1.000 & * & 1.000 & ** & 1.000 & ** & 1.000 & ** & 1.000 & ** & 1.000 & ** \\
AR(1) & 0.758 & * & 1.230 & ** & 0.789 & ** & 1.226 & ** & 0.775 & ** & 1.209 & ** \\
MIDAS & \cellcolor{black!15}\textbf{0.533} & \cellcolor{black!15}** & 1.300 & & 0.596 & ** & 1.129 & * & 0.709 & ** & 1.170 & * \\
DFM [A] & 0.799 & * & 1.337 & & 0.980 & * & 1.320 & & 0.919 & * & 1.226 & \\
DFM [B] & 0.885 & & 1.221 & ** & 0.982 & * & 1.022 & ** & 0.948 & & 1.028 & ** \\
singleESN [A] & 0.721 & ** & 1.015 & ** & 0.597 & ** & \cellcolor{black!8}0.867 & \cellcolor{black!8}** & \cellcolor{black!15}\textbf{0.529} & \cellcolor{black!15}** & \cellcolor{black!8}0.863 & \cellcolor{black!8}** \\
singleESN [B] & 0.758 & * & \cellcolor{black!15}\textbf{0.921} & \cellcolor{black!15}** & 0.602 & ** & \cellcolor{black!15}\textbf{0.844} & \cellcolor{black!15}** & 0.561 & ** & 0.930 & ** \\
multiESN [A] & 0.802 & * & 1.250 & & 0.635 & ** & 0.874 & ** & 0.621 & ** & \cellcolor{black!15}\textbf{0.859} & \cellcolor{black!15}** \\
multiESN [B] & \cellcolor{black!8}0.590 & \cellcolor{black!8}** & \cellcolor{black!8}0.969 & \cellcolor{black!8}** & \cellcolor{black!15}\textbf{0.552} & \cellcolor{black!15}** & 0.895 & ** & \cellcolor{black!8}0.530 & \cellcolor{black!8}** & 0.921 & ** \\
\bottomrule
\end{tabular*}
\caption{\small Relative MSFE and Model Confidence Set (MCS) comparison between models in 1-step-ahead forecasting exercises. Unconditional mean MSFE is used as a reference. MCS columns show inclusion among best models: $\ast$ indicates inclusion at 90\%, $\ast\ast$ indicates inclusion at 75\% confidence.}
\label{table:1sa_GDP_smallMD}
\end{table}
The overall finding is that MFESN models perform excellent, and, when we exclude the 2007 fixed parameters setup, they perform the best. It is easy to see from Figure \ref{fig:1sa_forecasts_2007_GDP_smallMD} (a) why the 2007 fixed window estimation case is different from other cases: the 2008 Financial Crisis induced a deep drop in quarter-to-quarter GDP growth that was in stark contrast with previous business cycle fluctuations. By keeping model parameters fixed, and using only information from 1990 to 2007 -- periods where systematic fluctuations are small -- DFM and MFESN models are fit to produce smooth, low-volatility forecasts. MIDAS, on the other hand, yields an exponential smoothing which can be more responsive to changes in monthly and daily series. From Figure \ref{fig:1sa_forecasts_2007_GDP_smallMD} (b) and (c) it is possible to see that expanding and rolling window estimation resolves this weakness of state-space models. At the same time, the AR(1) model outperforms the unconditional mean only in the 2007 sample with fixed parameters, losing to the MIDAS model in all but one scenario.
Table \ref{table:1sa_GDP_smallMD} shows that MFESN models always perform better than the mean in terms of MSFE, something which no other model class achieves across all setups. In both expanding and rolling window setups they also always outperform the AR(1) model. Furthermore, at least one MFESN model for each subclass (single or multi-reservoir) is always included in the model confidence set at the highest confidence level. We remind again that the MCS test of \cite{Hansen2011} might be distorted due to the modest sample sizes considered, even more so in the 2011 test sample. To complement the MCS, we provide graphical tables for pairwise Modified Diebold-Mariano (MDM) tests, with 10\% level rejections highlighted in Figure \ref{fig:1sa_GDP_smallMD_diebold_mariano}, Supplementary Appendix \ref{subsec:additional_figures}. The MDM tests broadly agree with the results of Table \ref{table:1sa_GDP_smallMD}, although they do not account for multiple testing, and therefore cannot be interpreted as yielding subsets of the most accurate forecasting models in a statistical sense.
For multiple-steps-ahead forecasts, relative RMSFE and uMCS are reported in Tables \ref{table:msa_2007_GDP_smallMD} and \ref{table:msa_2011_GDP_smallMD}: we constrain our exercise to $h \in \{ 1, \ldots, 8 \}$ steps, since we are interested in GDP growth forecasts within 2 years. Note that for $h = 1$ our results are similar, but do not reduce to the one-step-ahead results. To make correct multistep RMSFE evaluations and execute the uMCS procedure one must select $h$ different vectors of residuals of the same length: this implies that residuals at the end of the forecasting sample must be trimmed off to compute short-term multistep RMSFEs that are comparable to the long-term ones. Generally, we notice that MIDAS, as well as S-MFESNs, provide the worst-performing multistep forecasts, with RMSFEs considerably exceeding the unconditional mean baseline after horizon 1. Figures \ref{fig:msa_2007_GDP_smallMD} and \ref{fig:msa_2011_GDP_smallMD} in Appendix~\ref{subsec:forecast_figures} reproduce the RMSFE numbers of the aforementioned tables graphically.
\begin{table}[!p]
\vspace{5em}
\centering
\setlength\tabcolsep{6.5pt}
\setlength\extrarowheight{2pt}
\linespread{1}\selectfont\centering
Multistep-ahead GDP Forecasting - Small-MD Dataset - 2007 Sample \\[5pt]
\begin{tabular*}{\textwidth}{@{}*{11}{c}}
\toprule
Setup & Model & \multicolumn{8}{c}{Horizon} & uMCS \\
\cline{3-10}
& & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & \\
\midrule
FIX & Mean & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & \\
FIX & AR(1) & 0.870 & \cellcolor{black!15}\textbf{0.950} & \cellcolor{black!15}\textbf{0.982} & \cellcolor{black!8}0.991 & 0.992 & 0.991 & 0.992 & 0.992 & ** \\
FIX & MIDAS & \cellcolor{black!8}0.823 & 1.672 & 2.737 & 1.816 & 2.213 & 2.791 & 1.888 & 1.921 & \\
FIX & DFM [A] & 0.890 & 0.969 & 1.014 & 1.077 & 1.341 & 1.701 & 2.001 & 2.180 & * \\
FIX & DFM [B] & 0.937 & 1.069 & 1.202 & 1.344 & 1.799 & 2.310 & 2.638 & 2.801 & \\
FIX & singleESN [A] & 0.852 & 0.994 & 0.995 & 0.995 & 0.993 & 0.991 & 0.991 & 0.991 & * \\
FIX & singleESN [B] & 0.871 & 0.986 & 0.989 & \cellcolor{black!15}\textbf{0.989} & \cellcolor{black!15}\textbf{0.985} & \cellcolor{black!15}\textbf{0.981} & \cellcolor{black!15}\textbf{0.981} & \cellcolor{black!15}\textbf{0.981} & ** \\
FIX & multiESN [A] & 0.898 & 0.980 & 0.990 & \cellcolor{black!8}0.991 & \cellcolor{black!8}0.988 & \cellcolor{black!8}0.985 & \cellcolor{black!8}0.985 & \cellcolor{black!8}0.985 & ** \\
FIX & multiESN [B] & \cellcolor{black!15}\textbf{0.767} & \cellcolor{black!8}0.954 & \cellcolor{black!8}0.983 & \cellcolor{black!8}0.991 & 0.991 & 0.990 & 0.991 & 0.991 & ** \\
\midrule
EW & Mean & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & \cellcolor{black!15}\textbf{1.000} & \cellcolor{black!15}\textbf{1.000} & - \\
EW & AR(1) & 0.887 & 0.922 & 0.951 & \cellcolor{black!8}0.962 & \cellcolor{black!8}0.957 & \cellcolor{black!15}\textbf{0.981} & \cellcolor{black!8}1.001 & \cellcolor{black!8}1.008 & - \\
EW & MIDAS & 0.814 & 1.283 & 1.518 & 1.596 & 1.697 & 1.391 & 1.951 & 1.800 & - \\
EW & DFM [A] & 0.985 & 1.109 & 1.123 & 1.114 & 1.217 & 1.226 & 1.241 & 1.539 & - \\
EW & DFM [B] & 0.989 & 1.082 & 1.149 & 1.199 & 1.315 & 1.412 & 1.373 & 1.425 & - \\
EW & singleESN [A] & \cellcolor{black!8}0.771 & 1.260 & 1.485 & 1.564 & 2.070 & 2.728 & 2.550 & 2.834 & - \\
EW & singleESN [B] & 0.772 & 1.031 & 1.135 & 1.319 & 1.831 & 2.279 & 2.449 & 2.556 & - \\
EW & multiESN [A] & 0.792 & \cellcolor{black!8}0.897 & \cellcolor{black!8}0.941 & 0.976 & 1.015 & 1.240 & 1.377 & 1.227 & - \\
EW & multiESN [B] & \cellcolor{black!15}\textbf{0.740} & \cellcolor{black!15}\textbf{0.853} & \cellcolor{black!15}\textbf{0.894} & \cellcolor{black!15}\textbf{0.911} & \cellcolor{black!15}\textbf{0.873} & \cellcolor{black!8}0.993 & 1.020 & 1.020 & - \\
\midrule
RW & Mean & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & * \\
RW & AR(1) & 0.898 & \cellcolor{black!8}0.943 & \cellcolor{black!8}0.968 & \cellcolor{black!8}0.974 & \cellcolor{black!8}0.963 & \cellcolor{black!15}\textbf{0.968} & \cellcolor{black!8}0.970 & \cellcolor{black!15}\textbf{0.962} & ** \\
RW & MIDAS & 0.933 & 1.438 & 1.642 & 1.993 & 1.794 & 1.661 & 1.816 & 1.973 & * \\
RW & DFM [A] & 0.931 & 1.017 & 1.033 & 1.020 & 1.024 & \cellcolor{black!8}1.003 & \cellcolor{black!15}\textbf{0.918} & \cellcolor{black!8}1.062 & * \\
RW & DFM [B] & 0.942 & 0.973 & 0.970 & 1.045 & 1.059 & 1.203 & 1.225 & 1.263 & * \\
RW & singleESN [A] & \cellcolor{black!15}\textbf{0.714} & 1.320 & 1.693 & 1.972 & 2.733 & 3.669 & 3.391 & 3.719 & * \\
RW & singleESN [B] & 0.737 & 1.100 & 1.248 & 1.667 & 2.327 & 2.765 & 2.842 & 2.792 & * \\
RW & multiESN [A] & 0.773 & 0.972 & 1.053 & 1.111 & 1.187 & 1.293 & 1.505 & 1.131 & * \\
RW & multiESN [B] & \cellcolor{black!8}0.716 & \cellcolor{black!15}\textbf{0.895} & \cellcolor{black!15}\textbf{0.916} & \cellcolor{black!15}\textbf{0.926} & \cellcolor{black!15}\textbf{0.890} & 1.041 & 1.102 & 1.105 & ** \\
\bottomrule
\end{tabular*}
\caption{Relative RMSFE and Uniform Multi-Horizon Model Confidence Set (uMCS) comparison between models in multiple-steps-ahead forecasting exercises. Unconditional mean RMSFE used as reference. FIX: Fixed parameters, EW: Expanding window, and RW: Rolling window. uMCS columns show inclusion among best models: $\ast$ indicates inclusion at 90\% confidence, $\ast\ast$ indicates inclusion at 75\% confidence.}
\label{table:msa_2007_GDP_smallMD}
\end{table}
\begin{table}[!p]
\vspace{5em}
\centering
\setlength\tabcolsep{6.5pt}
\setlength\extrarowheight{2pt}
\linespread{1}\selectfont\centering
Multistep-ahead GDP Forecasting - Small-MD Dataset - 2011 Sample \\[5pt]
\begin{tabular*}{\textwidth}{@{}*{11}{c}}
\toprule
Setup & Model & \multicolumn{8}{c}{Horizon} & uMCS \\
\cline{3-10}
& & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & \\
\midrule
FIX & Mean & 1.000 & \cellcolor{black!8}1.000 & 1.000 & \cellcolor{black!8}1.000 & \cellcolor{black!15}\textbf{1.000} & 1.000 & \cellcolor{black!8}1.000 & \cellcolor{black!8}1.000 & ** \\
FIX & AR(1) & 1.119 & 1.031 & 1.008 & 1.001 & \cellcolor{black!8}1.001 & \cellcolor{black!8}0.999 & \cellcolor{black!15}\textbf{0.999} & \cellcolor{black!15}\textbf{0.998} & * \\
FIX & MIDAS & 1.090 & 1.721 & 1.793 & 2.203 & 2.363 & 1.997 & 2.846 & 2.328 & \\
FIX & DFM [A] & 1.112 & 1.051 & 0.999 & 1.079 & 1.084 & 1.025 & 1.020 & 1.061 & * \\
FIX & DFM [B] & 1.058 & \cellcolor{black!15}\textbf{0.945} & \cellcolor{black!15}\textbf{0.916} & 1.003 & 1.012 & \cellcolor{black!15}\textbf{0.970} & 1.038 & 1.033 & ** \\
FIX & singleESN [A] & \cellcolor{black!8}0.978 & 1.705 & 2.561 & 2.704 & 3.314 & 3.151 & 2.999 & 3.316 & \\
FIX & singleESN [B] & \cellcolor{black!15}\textbf{0.930} & 1.095 & 1.885 & 2.356 & 2.650 & 2.704 & 2.880 & 2.844 & ** \\
FIX & multiESN [A] & 1.059 & 1.148 & 1.262 & 1.312 & 1.339 & 1.409 & 1.424 & 1.162 & \\
FIX & multiESN [B] & 0.981 & 1.007 & \cellcolor{black!8}0.985 & \cellcolor{black!15}\textbf{0.994} & 1.008 & \cellcolor{black!8}0.999 & \cellcolor{black!15}\textbf{0.999} & \cellcolor{black!15}\textbf{0.998} & ** \\
\midrule
EW & Mean & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & \cellcolor{black!8}1.000 & 1.000 & - \\
EW & AR(1) & 1.117 & 1.033 & 1.011 & 1.002 & 1.007 & 1.003 & 1.004 & 1.003 & - \\
EW & MIDAS & 1.005 & 1.382 & 1.339 & 1.354 & 1.609 & 1.444 & 1.803 & 1.263 & - \\
EW & DFM [A] & 1.144 & 1.132 & 1.057 & 1.093 & 1.076 & 1.067 & 1.038 & 1.016 & - \\
EW & DFM [B] & 0.985 & \cellcolor{black!15}\textbf{0.940} & \cellcolor{black!15}\textbf{0.918} & 0.995 & 1.010 & \cellcolor{black!15}\textbf{0.980} & 1.050 & \cellcolor{black!15}\textbf{0.971} & - \\
EW & singleESN [A] & 0.935 & 1.645 & 2.184 & 1.929 & 2.388 & 1.959 & 1.810 & 2.266 & - \\
EW & singleESN [B] & \cellcolor{black!15}\textbf{0.911} & 1.092 & 1.101 & 1.529 & 2.195 & 1.843 & 1.847 & 2.060 & - \\
EW & multiESN [A] & \cellcolor{black!8}0.922 & \cellcolor{black!8}0.965 & 1.089 & \cellcolor{black!8}0.978 & \cellcolor{black!15}\textbf{0.977} & 1.043 & 1.278 & \cellcolor{black!8}0.995 & - \\
EW & multiESN [B] & 0.944 & 0.992 & \cellcolor{black!8}0.978 & \cellcolor{black!15}\textbf{0.977} & \cellcolor{black!8}0.991 & \cellcolor{black!8}0.985 & \cellcolor{black!15}\textbf{0.990} & 0.996 & - \\
\midrule
RW & Mean & 1.000 & \cellcolor{black!8}1.000 & 1.000 & \cellcolor{black!8}1.000 & 1.000 & 1.000 & 1.000 & 1.000 & \\
RW & AR(1) & 1.080 & \cellcolor{black!8}1.000 & 0.984 & \cellcolor{black!15}\textbf{0.989} & \cellcolor{black!8}0.982 & 0.976 & \cellcolor{black!15}\textbf{0.963} & 0.968 & \\
RW & MIDAS & 1.051 & 1.303 & 1.310 & 1.674 & 1.762 & 1.467 & 1.643 & 1.463 & \\
RW & DFM [A] & 1.061 & 1.033 & 1.012 & 1.088 & 1.077 & 1.015 & 1.040 & 1.069 & \\
RW & DFM [B] & 0.947 & \cellcolor{black!15}\textbf{0.893} & \cellcolor{black!15}\textbf{0.901} & 1.009 & 1.040 & \cellcolor{black!15}\textbf{0.966} & 1.030 & \cellcolor{black!15}\textbf{0.949} & ** \\
RW & singleESN [A] & \cellcolor{black!8}0.919 & 1.788 & 2.359 & 2.483 & 2.981 & 2.401 & 2.234 & 2.690 & \\
RW & singleESN [B] & 0.944 & 1.132 & 1.214 & 1.762 & 2.608 & 2.552 & 2.517 & 2.541 & \\
RW & multiESN [A] & \cellcolor{black!15}\textbf{0.896} & 1.047 & 1.222 & 1.124 & 1.122 & 1.410 & 1.666 & 1.316 & \\
RW & multiESN [B] & 0.940 & 1.003 & \cellcolor{black!8}0.969 & \cellcolor{black!15}\textbf{0.989} & \cellcolor{black!15}\textbf{0.979} & \cellcolor{black!8}0.972 & \cellcolor{black!8}0.967 & \cellcolor{black!8}0.961 & ** \\
\bottomrule
\end{tabular*}
\caption{Relative RMSFE and Uniform Multi-Horizon Model Confidence Set (uMCS) comparison between models in multiple-steps-ahead forecasting exercises. Unconditional mean RMSFE used as reference. FIX: Fixed parameters, EW: Expanding window, and RW: Rolling window. uMCS columns show inclusion among best models: $\ast$ indicates inclusion at 90\%, $\ast\ast$ indicates inclusion at 75\% confidence.}
\label{table:msa_2011_GDP_smallMD}
\end{table}
For MIDAS, we have already discussed how the existence of multiple loss minima can generate numerical instabilities. Model re-fitting at each horizon can amplify this problem, as the loss landscape itself changes as new observations are added to the fitting sample. We provide more discussions in Supplementary Appendix~\ref{sec:midas_robustness}.
In the case of S-MFESN models, the reason is structural: we have discussed how in our framework multistep MFESN forecasting entails iterating the state map, which can have multiple attraction (stable) points. If the hyperparameters and estimated full model $\widehat{W}$s jointly do not define a contraction, the limit of the multistep forecast does not have to be the estimated MFSEN model intercept. However, Figures \ref{fig:msa_2007_GDP_smallMD} and \ref{fig:msa_2011_GDP_smallMD} show that our M-MFESN models, multiESN [A] and multiESN [B], both perform on par or better than DFM models even after horizon $h = 4$. For example, in the 2007 expanding and rolling window experiments, multiESN [B] is able to outperform both DFMs and an unconditional mean forecast by meaningful margins for forecasts up to a year into the future.
\subsubsection{Medium Dataset}
\begin{table}[t!]
\centering
\setlength\tabcolsep{3.3pt}
\setlength\extrarowheight{2pt}
\linespread{1}\selectfont\centering
1-Step-ahead GDP Forecasting - Medium-MD Dataset \\[5pt]
\begin{tabular*}{\textwidth}{@{}*{13}{c}}
\toprule
& \multicolumn{4}{c}{Fixed Parameters} & \multicolumn{4}{c}{Expanding Window} & \multicolumn{4}{c}{Rolling Window} \\
\cline{2-5} \cline{6-9} \cline{10-13}
Model & \multicolumn{2}{c}{2007} & \multicolumn{2}{c}{2011} & \multicolumn{2}{c}{2007} & \multicolumn{2}{c}{2011} & \multicolumn{2}{c}{2007} & \multicolumn{2}{c}{2011} \\
\cline{2-3} \cline{4-5} \cline{6-7} \cline{8-9} \cline{10-11} \cline{12-13}
& MSFE & MCS & MSFE & MCS & MSFE & MCS & MSFE & MCS & MSFE & MCS & MSFE & MCS \\
\midrule
Mean & 1.000 & & \cellcolor{black!8}1.000 & \cellcolor{black!8}** & 1.000 & ** & 1.000 & ** & 1.000 & ** & 1.000 & ** \\
AR(1) & \cellcolor{black!8}0.758 & \cellcolor{black!8}* & 1.230 & ** & 0.789 & ** & 1.226 & ** & 0.775 & * & 1.209 & ** \\
DFM [A] & 0.841 & * & 1.325 & * & 0.682 & ** & 1.272 & ** & 0.747 & * & 1.517 & ** \\
DFM [B] & 1.118 & * & 1.408 & ** & 0.821 & * & 1.117 & ** & 0.926 & & 1.186 & ** \\
singleESN [A] & 0.967 & * & 1.717 & * & 0.775 & ** & 1.072 & ** & 0.791 & * & 1.493 & * \\
singleESN [B] & 0.826 & * & 1.278 & ** & 0.655 & ** & 1.028 & ** & 0.561 & ** & 0.944 & ** \\
multiESN [A] & 0.901 & * & 1.080 & ** & \cellcolor{black!8}0.618 & \cellcolor{black!8}** & \cellcolor{black!8}0.913 & \cellcolor{black!8}** & \cellcolor{black!8}0.556 & \cellcolor{black!8}** & \cellcolor{black!8}0.884 & \cellcolor{black!8}** \\
multiESN [B] & \cellcolor{black!15}\textbf{0.682} & \cellcolor{black!15}** & \cellcolor{black!15}\textbf{0.748} & \cellcolor{black!15}** & \cellcolor{black!15}\textbf{0.587} & \cellcolor{black!15}** & \cellcolor{black!15}\textbf{0.774} & \cellcolor{black!15}** & \cellcolor{black!15}\textbf{0.547} & \cellcolor{black!15}** & \cellcolor{black!15}\textbf{0.728} & \cellcolor{black!15}** \\
\bottomrule
\end{tabular*}
\caption{\small Relative MSFE and Model Confidence Set (MCS) comparison between models in 1-step-ahead forecasting exercises. Unconditional mean MSFE is used as a reference. MCS columns show inclusion among best models: $\ast$ indicates inclusion at 90\% confidence, $\ast\ast$ indicates inclusion at 75\% confidence.}
\label{table:1sa_GDP_mediumMD}
\end{table}
We now present the results for the Medium-MD dataset, which includes more than 30 regressors and many high-frequency daily series. The same metrics as in the previous subsection are used for this dataset to evaluate the relative performance of different methods.
The main difference in our empirical exercises is that now we a priori exclude MIDAS from the set of forecasting methods as explained in detail in Subsection~\ref{Benchmarks}. Table \ref{table:1sa_GDP_mediumMD} showcases the relative performance of DFM and MFESN models in the Medium-MD forecast setup. We find that the MFESN model multiESN [B] performs best in all setups, particularly under fixed parameters, where MCS testing reveals that it is the only model included at a 75\% confidence level. Of course, for the MCS results we must again take into account the relatively small sample size, which could distort the selection of best model subsets. MDM tests of Figure \ref{fig:1sa_GDP_mediumMD_diebold_mariano} in Supplementary Appendix \ref{subsec:additional_figures} largely agree with the MCS results: in the fixed parameter setup any pairwise comparison of an alternative model against MFESN multiESN [B] is rejected in favor of the latter.
A visual inspection of one-step-ahead forecasts in Figures \ref{fig:1sa_forecasts_2007_GDP_mediumMD} and \ref{fig:1sa_forecasts_2011_GDP_mediumMD} in Appendix~\ref{subsec:forecast_figures} also shows that DFM models estimated over the Medium-MD datasets produce forecasts with larger variability than MFESN methods, which is likely the key driver of the difference in performance.
The multistep-ahead experiments are run as for the Small-MD dataset, with a maximum horizon of 8 quarters. Tables \ref{table:msa_2007_GDP_mediumMD} and \ref{table:msa_2011_GDP_mediumMD} present the relative RMSFE performance of multistep forecasts for all models, and we use Figures \ref{fig:msa_2007_GDP_mediumMD} and \ref{fig:msa_2011_GDP_mediumMD} of RMSFEs as references for our discussion, available in Appendix~\ref{subsec:forecast_figures}. What can be seen visually -- and is also reproduced in the Tables -- is that multi-reservoir MFESN models and DFM model [A] have the best performance up to 4 quarters ahead; overall, taking into account also the longer term, expanding or rolling window estimation of model multiESN [B] yields the best forecasting results in the 2007 sample setup. The post-crisis 2011 sample setup makes comparison harder, as DFM and M-MFESN models largely produce results in line with the unconditional sample mean. This evaluation is confirmed by uMCS tests, consistently with the multistep results obtained with the Small-MD dataset.
\begin{table}[!p]
\vspace{5em}
\centering
\setlength\tabcolsep{6.5pt}
\setlength\extrarowheight{2pt}
\linespread{1}\selectfont\centering
Multistep-ahead GDP Forecasting - Medium-MD Dataset - 2007 Sample \\[5pt]
\begin{tabular*}{\textwidth}{@{}*{11}{c}}
\toprule
Setup & Model & \multicolumn{8}{c}{Horizon} & uMCS \\
\cline{3-10}
& & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & \\
\midrule
FIX & Mean & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & * \\
FIX & AR(1) & \cellcolor{black!8}0.870 & \cellcolor{black!8}0.950 & 0.982 & 0.991 & 0.992 & 0.991 & 0.992 & 0.992 & \\
FIX & DFM [A] & 0.914 & \cellcolor{black!15}\textbf{0.947} & \cellcolor{black!15}\textbf{0.955} & 0.988 & 1.015 & 1.027 & 1.034 & 0.995 & ** \\
FIX & DFM [B] & 1.046 & 1.204 & 1.293 & 1.341 & 1.649 & 1.984 & 2.101 & 2.070 & * \\
FIX & singleESN [A] & 0.985 & 0.995 & 0.995 & 0.995 & 0.994 & 0.992 & 0.992 & 0.992 & * \\
FIX & singleESN [B] & 0.912 & 0.985 & 0.985 & \cellcolor{black!15}\textbf{0.985} & \cellcolor{black!15}\textbf{0.980} & \cellcolor{black!15}\textbf{0.976} & \cellcolor{black!15}\textbf{0.976} & \cellcolor{black!15}\textbf{0.976} & * \\
FIX & multiESN [A] & 0.950 & 0.993 & 0.994 & 0.994 & 0.992 & 0.990 & 0.990 & 0.990 & * \\
FIX & multiESN [B] & \cellcolor{black!15}\textbf{0.826} & 0.972 & \cellcolor{black!8}0.988 & \cellcolor{black!8}0.990 & \cellcolor{black!8}0.989 & \cellcolor{black!8}0.986 & \cellcolor{black!8}0.985 & \cellcolor{black!8}0.985 & * \\
\midrule
EW & Mean & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & \cellcolor{black!15}\textbf{1.000} & \cellcolor{black!15}\textbf{1.000} & - \\
EW & AR(1) & 0.887 & 0.922 & \cellcolor{black!8}0.951 & \cellcolor{black!8}0.962 & \cellcolor{black!8}0.957 & \cellcolor{black!8}0.981 & \cellcolor{black!8}1.001 & \cellcolor{black!8}1.008 & - \\
EW & DFM [A] & 0.805 & \cellcolor{black!8}0.916 & 0.978 & 1.038 & 1.077 & 1.126 & 1.077 & 1.073 & - \\
EW & DFM [B] & 0.893 & 1.134 & 1.418 & 1.567 & 2.238 & 2.964 & 3.375 & 3.629 & - \\
EW & singleESN [A] & 0.879 & 1.125 & 1.305 & 1.442 & 1.860 & 2.166 & 2.361 & 2.443 & - \\
EW & singleESN [B] & 0.802 & 1.174 & 1.439 & 1.744 & 2.305 & 2.869 & 2.935 & 3.167 & - \\
EW & multiESN [A] & \cellcolor{black!8}0.780 & 0.935 & 1.012 & 1.005 & 1.093 & 1.337 & 1.328 & 1.313 & - \\
EW & multiESN [B] & \cellcolor{black!15}\textbf{0.760} & \cellcolor{black!15}\textbf{0.874} & \cellcolor{black!15}\textbf{0.911} & \cellcolor{black!15}\textbf{0.891} & \cellcolor{black!15}\textbf{0.863} & \cellcolor{black!15}\textbf{0.971} & 1.030 & 1.051 & - \\
\midrule
RW & Mean & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & \cellcolor{black!8}1.000 & \cellcolor{black!8}1.000 & \\
RW & AR(1) & 0.898 & 0.943 & 0.968 & 0.974 & \cellcolor{black!8}0.963 & \cellcolor{black!15}\textbf{0.968} & \cellcolor{black!15}\textbf{0.970} & \cellcolor{black!15}\textbf{0.962} & \\
RW & DFM [A] & 0.837 & \cellcolor{black!8}0.913 & \cellcolor{black!8}0.924 & \cellcolor{black!8}0.954 & 1.012 & 0.997 & 1.018 & 1.005 & \\
RW & DFM [B] & 0.932 & 1.116 & 1.232 & 1.414 & 1.952 & 2.704 & 3.183 & 3.294 & \\
RW & singleESN [A] & 0.873 & 1.274 & 1.530 & 1.652 & 2.095 & 2.575 & 2.786 & 3.014 & \\
RW & singleESN [B] & \cellcolor{black!8}0.732 & 1.190 & 1.490 & 1.712 & 2.218 & 2.861 & 2.967 & 3.094 & \\
RW & multiESN [A] & \cellcolor{black!8}0.732 & 0.914 & 0.960 & 1.011 & 1.202 & 1.618 & 1.683 & 1.572 & \\
RW & multiESN [B] & \cellcolor{black!15}\textbf{0.731} & \cellcolor{black!15}\textbf{0.871} & \cellcolor{black!15}\textbf{0.875} & \cellcolor{black!15}\textbf{0.844} & \cellcolor{black!15}\textbf{0.771} & \cellcolor{black!8}0.971 & 1.014 & 1.014 & ** \\
\bottomrule
\end{tabular*}
\caption{Relative RMSFE and Uniform Multi-Horizon Model Confidence Set (uMCS) comparison between models in multiple-steps-ahead forecasting exercises. Unconditional mean RMSFE used as reference. FIX: Fixed parameters, EW: Expanding window, and RW: Rolling window. uMCS columns show inclusion among best models: $\ast$ indicates inclusion at 90\% confidence, $\ast\ast$ indicates inclusion at 75\% confidence.}
\label{table:msa_2007_GDP_mediumMD}
\end{table}
\begin{table}[!p]
\vspace{5em}
\centering
\setlength\tabcolsep{6.5pt}
\setlength\extrarowheight{2pt}
\linespread{1}\selectfont\centering
Multistep-ahead GDP Forecasting - Medium-MD Dataset - 2011 Sample \\[5pt]
\begin{tabular*}{\textwidth}{@{}*{11}{c}}
\toprule
Setup & Model & \multicolumn{8}{c}{Horizon} & uMCS \\
\cline{3-10}
& & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & \\
\midrule
FIX & Mean & 1.000 & 1.000 & 1.000 & \cellcolor{black!8}1.000 & \cellcolor{black!15}\textbf{1.000} & 1.000 & \cellcolor{black!8}1.000 & \cellcolor{black!8}1.000 & * \\
FIX & AR(1) & 1.119 & 1.031 & 1.008 & 1.001 & \cellcolor{black!8}1.001 & \cellcolor{black!8}0.999 & \cellcolor{black!15}\textbf{0.999} & \cellcolor{black!15}\textbf{0.998} & ** \\
FIX & DFM [A] & 1.126 & 0.987 & \cellcolor{black!8}0.962 & 1.054 & 1.031 & \cellcolor{black!15}\textbf{0.988} & 1.001 & 1.002 & ** \\
FIX & DFM [B] & 1.149 & 0.987 & \cellcolor{black!15}\textbf{0.885} & 1.064 & 1.142 & 1.134 & 1.273 & 1.296 & \\
FIX & singleESN [A] & 1.283 & 1.921 & 2.527 & 3.038 & 3.285 & 3.154 & 3.193 & 3.655 & \\
FIX & singleESN [B] & 1.059 & 1.523 & 1.918 & 2.417 & 2.812 & 2.683 & 2.703 & 2.970 & \\
FIX & multiESN [A] & 1.011 & 1.061 & 1.434 & 1.477 & 1.748 & 2.030 & 2.023 & 1.994 & \\
FIX & multiESN [B] & \cellcolor{black!15}\textbf{0.841} & \cellcolor{black!15}\textbf{0.945} & 0.997 & \cellcolor{black!15}\textbf{0.978} & 1.004 & 1.015 & 1.013 & 1.014 & ** \\
\midrule
EW & Mean & 1.000 & 1.000 & \cellcolor{black!8}1.000 & \cellcolor{black!8}1.000 & \cellcolor{black!15}\textbf{1.000} & \cellcolor{black!8}1.000 & \cellcolor{black!8}1.000 & \cellcolor{black!8}1.000 & - \\
EW & AR(1) & 1.117 & 1.033 & 1.011 & 1.002 & 1.007 & 1.003 & 1.004 & 1.003 & - \\
EW & DFM [A] & 1.092 & \cellcolor{black!8}0.942 & \cellcolor{black!15}\textbf{0.944} & 1.049 & 1.026 & \cellcolor{black!15}\textbf{0.994} & \cellcolor{black!15}\textbf{0.996} & \cellcolor{black!15}\textbf{0.999} & - \\
EW & DFM [B] & 0.971 & 1.046 & 1.031 & 1.114 & 1.238 & 1.116 & 1.223 & 1.310 & - \\
EW & singleESN [A] & 1.039 & 1.451 & 1.980 & 2.385 & 2.699 & 2.353 & 2.506 & 2.608 & - \\
EW & singleESN [B] & 0.992 & 1.828 & 2.465 & 3.072 & 3.547 & 3.357 & 3.368 & 3.610 & - \\
EW & multiESN [A] & \cellcolor{black!8}0.934 & 1.014 & 1.391 & 1.252 & 1.371 & 1.369 & 1.228 & 1.279 & - \\
EW & multiESN [B] & \cellcolor{black!15}\textbf{0.857} & \cellcolor{black!15}\textbf{0.931} & 1.003 & \cellcolor{black!15}\textbf{0.973} & \cellcolor{black!8}1.002 & 1.009 & 1.025 & 1.029 & - \\
\midrule
RW & Mean & 1.000 & 1.000 & 1.000 & \cellcolor{black!8}1.000 & \cellcolor{black!8}1.000 & 1.000 & 1.000 & 1.000 & ** \\
RW & AR(1) & 1.080 & 1.000 & \cellcolor{black!8}0.984 & \cellcolor{black!15}\textbf{0.989} & \cellcolor{black!15}\textbf{0.982} & \cellcolor{black!15}\textbf{0.976} & \cellcolor{black!15}\textbf{0.963} & \cellcolor{black!15}\textbf{0.968} & ** \\
RW & DFM [A] & 1.113 & 0.982 & \cellcolor{black!15}\textbf{0.927} & 1.038 & 1.030 & 0.997 & 1.016 & 1.028 & * \\
RW & DFM [B] & 0.881 & 0.996 & 1.021 & 1.098 & 1.150 & 1.114 & 1.114 & 1.212 & ** \\
RW & singleESN [A] & 1.193 & 2.267 & 3.265 & 3.580 & 4.090 & 3.790 & 4.015 & 4.562 & \\
RW & singleESN [B] & 0.927 & 1.933 & 2.612 & 3.265 & 3.753 & 3.567 & 3.556 & 3.792 & \\
RW & multiESN [A] & 0.900 & 1.049 & 1.500 & 1.465 & 1.789 & 1.707 & 1.505 & 1.462 & \\
RW & multiESN [B] & \cellcolor{black!15}\textbf{0.816} & \cellcolor{black!15}\textbf{0.916} & 0.977 & 1.009 & \cellcolor{black!15}\textbf{0.982} & \cellcolor{black!8}0.988 & \cellcolor{black!8}0.974 & \cellcolor{black!8}0.981 & ** \\
\bottomrule
\end{tabular*}
\caption{Relative RMSFE and Uniform Multi-Horizon Model Confidence Set (uMCS) comparison between models in multiple-steps-ahead forecasting exercises. Unconditional mean RMSFE used as reference. FIX: Fixed parameters, EW: Expanding window, and RW: Rolling window. uMCS columns show inclusion among best models: $\ast$ indicates inclusion at 90\% confidence, $\ast\ast$ indicates inclusion at 75\% confidence.}
\label{table:msa_2011_GDP_mediumMD}
\end{table}
\section{Conclusions}
\label{Conclusions}
Macroeconomic forecasting -- especially long-term forecasting of macroeconomic aggregates -- is a topic of crucial importance for institutional policymakers, private companies, and economic researchers. Given the modern-day availability of ``big data'' resources, methods capable of integrating heterogeneous data sources are increasingly sought to provide more precise and robust forecasts.
This paper presents a new methodological framework inspired by the Reservoir Computing literature to deal with data sampled at multiple frequencies and with multiple-step-ahead forecasts. We have then taken Echo State Networks -- a type of RC models -- and formally extended them to allow the modeling of data with multiple release frequencies. Our discussion encompasses model fitting, hyperparameter tuning, and forecast computation. As a result, we provide two classes of models, single- and multiple reservoir multi-frequency ESNs, that can be effectively applied to our empirical setup: forecasting US GDP growth using monthly and daily data series. Along with the unconditional mean and AR(1) model, we considered two well-known methods, MIDAS and DFMs, as the current benchmarks available in the literature. In our applications, we find that MFESN models are computationally more efficient and easier to implement than DFMs and MIDAS, respectively, and perform better than or as well as benchmarks in terms of MSFE. These improvements are statistically significant in a number of setups, as shown by our MCS and MDM tests. Thus, we argue that our machine learning-based methodology can be a useful addition to the toolbox of contemporary macroeconomic forecasters.
Lastly, we wish to highlight the many potential areas of research that we believe would be interesting to explore in the future. We have not discussed the role of the distribution from which we sample the entries of the reservoir matrices. While it is known that these can have significant effects on the forecasting capacity of an ESN model, the literature lacks definitive theoretical results (even for dynamical systems applications) or systematic studies with stochastic inputs and targets. The hyperparameter tuning routine we have developed neither allows separating individual hyperparameters nor does it tackle the identification problem. Moreover, we assume that the ridge regression penalty strength, $\lambda$, is tuned \textit{ex ante}: it would be interesting and desirable to understand if it is possible to jointly tune $\lambda$ and $\boldsymbol{\varphi}$, or rather if one can fully separate their selection. In our preliminary experiments, we have noticed that the roles of the ridge penalty and the input scaling, for example, cannot be trivially disentangled -- thus prompting the $\psi$-form normalization. Model selection for the dimension of MFESN models is another question that would be key to exploring and designing more efficient and effective ESN models, especially when dealing with multiple frequencies and reservoirs. Finally, practitioners may be interested in identifying the combination of frequencies in the regressor series that would lead to the most accurate GDP forecasts produced by MFESN models.
\newpage
{\footnotesize
\bibliography{GOLibrary.bib}
}