Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
131,152 characters · 22 sections · 87 citation commands
Reservoir Computing for Macroeconomic Forecasting with Mixed Frequency Data
Key Words: Reservoir Computing, Echo State Networks, Forecasting, US Output Growth, GDP, Mixed-frequency Data, Time Series, Multi-Frequency Echo State Network, MIDAS, DFM \\
JEL: C53, C45, E17
{ {\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. {[email removed] }} \addtocounter{footnote}{1} \footnotetext{ Department of Statistical Science, UCL, Gower Str., London WC1E 6BT, UK. {[email removed] }} \addtocounter{footnote}{1} \footnotetext{ Department of Statistics, Athens University of Economics and Business, 10434 Athens, Greece. {[email removed] }} \addtocounter{footnote}{1} \footnotetext{ Faculty of Mathematics and Statistics, University of St. Gallen, Bodanstrasse 6, CH-9000 St. Gallen, Switzerland. {[email removed] }} \addtocounter{footnote}{1} \footnotetext{ Honorary Associate Professor, Department of Statistics, University of Warwick, Coventry CV4 7AL, UK. {[email removed] }} \addtocounter{footnote}{1} \footnotetext{ Global Development Institute (GDI), University of Manchester, Manchester M13 9PL, UK. {[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
\setcounter{footnote}{0}
\setstretch{1.1}
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 (Hindrayanto2016) and the recent pandemic (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 (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 Borio2011, Borio2013,Morley2015; we also refer the reader to Fuleky2020 for more details regarding high-dimensional data and to Armesto2010 and 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 maass1, DynamicalSystemsMaass, Crutchfield2010, and lukosevicius, tanaka:review for reviews). Generally speaking, RC is a versatile class of recurrent neural network (RNN) models (see 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 (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 RC6, RC7, RC9, RC10, RC12, RC20, and references therein for more details), and excellent performance in forecasting, classification, and learning of dynamical systems (see Hart2021,RC21). While conventional RNNs have been adopted for macroeconomic forecasting in a few instances (see, for example, 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 Ghysels2004,Ghysels2007 has been adopted widely for macroeconomic forecasting with mixed-frequency data (see for instance 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 (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, 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 geweke1977dynamic and sargent1977business. DFMs have become the standard workhorse for macroeconomic nowcasting and prediction (for more details, we refer the reader to 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 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 Banbura2013 for a detailed review and 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 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) presents reservoir models and discusses their advantages, as well as estimation, hyperparameter tuning, penalization and nonlinear multistep forecasting. In Section (ref), 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) 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) 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}}.
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 lowest sampling frequency, which we call in what follows the 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 tempo index, written $t, {\ast}\vert{{\kappa}} \,$, where $\kappa$ is the frequency multiplier. Our tempo notation assumes that low- and high-frequency series are sampled with temporal alignment: this means that the reference time index $t$ and the tempo index ${\ast}\vert{{\kappa}}$ have the following properties.
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). 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 Clements2008,Clements2009,Chen2010,Jardet2020 and references therein.
In this section, we introduce {\it reservoir computing} models 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 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.
Reservoir computing (RC) models are nonlinear state-space systems that, in the forecasting setting, are defined by the following equations:
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 reservoir map, and the observation map $h_{\boldsymbol{\theta}}: \mathbb{R}^N \to \mathbb{R}^J$, $J\in \mathbb{N}^+$ is referred to as the 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 input and the output (target) of the system, respectively, and $(\boldsymbol{x}_t)_{t \in \mathbb{Z}}$, $\boldsymbol{x}_t \in \mathbb{R}^N$, are the associated reservoir states. In (ref), $(\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 RC8, and estimation and generalization error bounds have been established in RC10, RC12.
In the case of an ESN model, the state and observation equations (ref)-(ref) are given by
where $A \in \mathbb{M}_{N}$ is the reservoir matrix, $C \in \mathbb{M}_{N,K}$ is the input matrix, $\boldsymbol{\zeta} \in \mathbb{R}^N$ is the input shift, $\alpha \in [0,1)$ is the leak rate and $W \in \mathbb{M}_{N,J}$ are the 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 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} 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 (ref)-(ref) hold (see 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 (ref) holds. The result in Corollary 3.2 in RC7, which is also valid for the case of ESNs with the leak rate, shows that the sufficient condition of the ESP associated with (ref) 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 jaeger2001, Jaeger04, Buehner:ESN, zhang:echo, Yildiz2012, Wainrib2016, Manjunath:Jaeger for more details. The result in Corollary 3.2 in RC7 also shows that this condition implies the so-called {\it fading memory property} 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 RC16 proves that the condition $\left\|A\right\|_2<1$ guarantees variance stationarity of the states associated with variance stationary inputs. Moreover, 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 (ref) 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 RC23 and nonlinear ESNs RC15 as we discuss in Section (ref).
We now discuss in detail the estimation of coefficients $W$ in (ref). 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 (ref) as:
Collect the states and the targets into the state and the observation matrices, respectively, as
Consider the ridge regression estimator for $W$ given by
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 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 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 (ref). Finally, as mentioned in the properties of reservoir systems in Subsection (ref), 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 (ref) as the fixed-parameter estimator. In our empirical analyses, we also implement expanding and rolling window estimation strategies which update $\widehat{W}_\lambda$ as new observations become available (we refer the reader to Supplementary Appendix (ref) 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.
As discussed in Subsection (ref), 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 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
where $\rho(\widetilde{A})$ denotes the spectral radius of $\widetilde{A}$. As discussed in the properties of reservoir systems in Subsection (ref), 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 (ref) as
We refer to tuple $\boldsymbol{\varphi}:=(\alpha, \rho, \gamma, \omega)$ as the hyperparameters of the ESN. Specifically, $\alpha \in [0,1)$ is the leak rate and $\rho \in \mathbb{R}_+$ is called the spectral radius of the reservoir matrix, $\gamma \in \mathbb{R}_+$ is the input scaling, and $\omega \in \mathbb{R}_+$ is the shift scaling. The choice of the hyperparameters determines the properties of the state map. For simplicity, in Section (ref), we choose the hyperparameters based on the empirical ESN literature. In Supplementary Appendix (ref), we also propose a general though more computationally intensive procedure to select hyperparameters in a data-driven way that could be interesting to practitioners.
To apply ridge estimator (ref), 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 kockPenalizedTimeSeries2020,Ballarin2023a with their validity established in 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) for additional details.
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 RC12, RC24. Although it is customary in the empirical literature to take $N$ as large as possible Lukosevicius2012, some recent literature discusses both the statistical risk bounds and the approximation-risk trade-off bounds for various RC families (see RC10 and RC24 for details). Under simplified assumptions that $\alpha = 0$ and $\rho = 0$ in (ref), ESNs have a natural connection to random-weights neural networks caoReviewNeuralNetworks2018 and random projection regression 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 belloniNewAsymptoticTheory2015a. 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 hastie2009elements). Recent work has studied more in-depth the link between regularization and generalization: hastieSurprisesHighdimensionalRidgeless2022 show that “ridgeless", that is interpolation, solutions can be optimal in some scenarios. However, in our empirical evaluations in Section (ref), cross-validation consistently selects non-zero ridge penalties, confirming that ridge penalization plays an important role in ESN forecasting performance.
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
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$,
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
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), our proposed method shows superior performance with respect to the standard state-of-the-art benchmarks.
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 (ref) while simultaneously capitalizing on the available results using ESNs in the forecasting of dynamical systems. More explicitly, we add to the ESN specification (ref)-(ref) 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 (ref), while the target sequence is the same as the input sequence $(\boldsymbol{z}_t)_{t \in \mathbb{Z}}$,
Here, we use symbol $\mathcal{W}$ for the output coefficients to separate this case from the general ESN equations (ref)-(ref). In (ref), $(\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 (ref) is determined by (ref), and it is possible to re-feed the forecasted variables back into the state equation as inputs. This yields the following state recursion:
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” 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), the $h > 1$ steps autonomous state iteration is given by
and
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 (ref), we notice that for the conditional forecast of the states, it holds that
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 (ref) 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 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 (ref) with the help of (ref) 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 (ref) and the assumption of independence in the zero-mean noise $\boldsymbol{\epsilon}_{t+h}$, we write
and use the approximation
which originates from
where $\boldsymbol{u}_t$ is assumed to be zero-mean. The validity of (ref) 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.
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 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) 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 Single-Reservoir Multi-Frequency Echo State Networks (S-MFESNs). The second group, referred to as 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.
Recall that, in the temporal notation of Definition (ref), 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
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
where $\mathcal{W}\in \mathbb{M}_{N, \sum_{l=1}^{L}n_l}$ and $s >0$. We term this class of MFESN models the Single-Reservoir Multi-Frequency ESNs (S-MFESNs).
Notice that equations (ref)-(ref) 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 state alignment scheme. An 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
\paragraph{Estimation of aligned S-MFESN.} Both coefficient matrices ${W}$ and ${\mathcal{W}}$ can be estimated as explained in Subsection (ref) under appropriate choices of corresponding penalty strengths. In particular, in order to obtain $\widehat{\mathcal{W}}$, the state and the observation matrices in (ref) are given by
while
are used for the estimation of $\widehat{W}$. We note that the state equation (ref) 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)).
\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 (ref) is given by
which, composed with itself exactly $\kappa_{\max}$ times, yields the target-frequency-aligned autonomous state transition map
Finally, from (ref) the $h$-steps ahead low-frequency forecasts, $h\in \mathbb{N}$, can be computed as
Figure (ref) gives a graphical diagram of the 1-step forecasting procedure for an S-MFESN. Additionally, Figure (ref) in Supplementary Appendix (ref) provides a similar diagram for the case of multistep forecasts.
The following example illustrates this proposed forecasting strategy for the case of quarterly GDP forecasting using monthly and daily series inputs.
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 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
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) 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 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 aligned M-MFESN observation equation is given by
\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 (ref) are constructed as
while with the notation as in (ref)
are used for the estimation of $\widehat{W}$. Again, the state equations (ref) of M-MFESN can be started with $\boldsymbol{x}^{(l)}_{0,{0}\vert{\kappa_{l}}}=\mathbf{0}_{N_l}$ (see Subsection (ref) 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
The target-frequency-aligned autonomous state transition map associated with each frequency $l$ is hence defined as
Finally, from (ref) the $h$-steps ahead forecasts can be computed as
In Figure (ref) 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) in Supplementary Appendix (ref) provides a similar diagram for the case of multistep forecasting.
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.
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 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) in Appendix (ref).
The selection of predictors follows the seminal work by 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 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 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 Andreou2013,Jardet2020.
We follow instructions by 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 Borio2002,Gramlich2010,Qin2022). In addition to disequilibrium indicators, realized stock market volatility has been found to improve macroeconomic predictions 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 bollerslev:garch.\footnote{We include a control 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 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 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 ex ante the number of working days in 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).
In this section, we present the set of models that we use throughout our empirical exercises. For a general overview, Table (ref) 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).
\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 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 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 Ghysels2004, Ghysels2007. Our dynamic MIDAS specification includes autoregressive lags of the target series and uses an Almon weighting scheme. As shown in 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). To make optimization more efficient, we use explicit expressions for MIDAS loss gradients as in 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) 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 geweke1977dynamic and 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 forni2005generalized, doz2011two, stock2016dynamic). We generalize the standard two-frequency DFM modeling setup 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). Two distinct DFM specifications are used. The first one termed DFM [A] uses the standard linear aggregation scheme, as provided in Example (ref), while the second is a variation that implements an Almon weighting scheme as presented in Example (ref) (we name it DFM [B]). The latter is similar to a MIDAS-type aggregation scheme 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 stock2016dynamic), commonly used macroeconomic panels feature a number of challenges, such as weak factors onatskiAsymptoticsPrincipalComponents2012. Moreover, as mentioned in Supplementary Appendix (ref), 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).
The first set of ESNs we propose is given by two S-MFESN models, based on Example (ref). 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 (ref)). 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). 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 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]).
We start by commenting on the computational efficiency of competing models and report execution times in seconds in Table (ref). 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.
Competing forecasts are compared using the Model Confidence Set (MCS) test derived in 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 Aparicio2018. Given this fact, we also conduct pairwise model comparison tests with the Modified Diebold-Mariano (MDM) test for predictive accuracy 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 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). Finally, we do not report uMCS test outcomes for the expanding window setup, as quaedvliegMultiHorizonForecastComparison2021 argues that in such context the test is invalid.
We begin our discussion of the Small-MD forecasting results by reviewing Table (ref). 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) and (ref) in Appendix (ref); additional plots for cumulative SFE, cumulative RMSFEs and other metrics can be found in Supplementary Appendix (ref).
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) (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) (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) 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 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), Supplementary Appendix (ref). The MDM tests broadly agree with the results of Table (ref), 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) and (ref): 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) and (ref) in Appendix (ref) reproduce the RMSFE numbers of the aforementioned tables graphically.
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). 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) and (ref) 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.
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). Table (ref) 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) in Supplementary Appendix (ref) 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) and (ref) in Appendix (ref) 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) and (ref) present the relative RMSFE performance of multistep forecasts for all models, and we use Figures (ref) and (ref) of RMSFEs as references for our discussion, available in Appendix (ref). 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.
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 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.