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.
97,884 characters
Quantum Reservoir Computing for Realized Volatility Forecasting
\title{Quantum Reservoir Computing for Realized Volatility Forecasting}
\author{Qingyu Li}
\affiliation{Institute of Fundamental and Frontier Sciences, University of Electronic Sciences and Technology of China, Chengdu 611731, China}
\affiliation{Key Laboratory of Quantum Physics and Photonic Quantum Information, Ministry of Education, University of Electronic Science and Technology of China, Chengdu 611731, China}
\author{Chiranjib Mukhopadhyay}
\affiliation{Institute of Fundamental and Frontier Sciences, University of Electronic Sciences and Technology of China, Chengdu 611731, China}
\affiliation{Key Laboratory of Quantum Physics and Photonic Quantum Information, Ministry of Education, University of Electronic Science and Technology of China, Chengdu 611731, China}
\author{Abolfazl Bayat}
\affiliation{Institute of Fundamental and Frontier Sciences, University of Electronic Sciences and Technology of China, Chengdu 611731, China}
\affiliation{Key Laboratory of Quantum Physics and Photonic Quantum Information, Ministry of Education, University of Electronic Science and Technology of China, Chengdu 611731, China}
\affiliation{Shimmer Center, Tianfu Jiangxi Laboratory, Chengdu 641419, China}
\author{Ali Habibnia}
\affiliation{Department of Economics, Virginia Tech, U.S.A. }
\affiliation{Dataism Laboratory for Quantitative Finance, Virginia Tech, U.S.A.}
\begin{abstract}
Recent advances in quantum computing have demonstrated its potential to significantly enhance the analysis and forecasting of complex classical data. Among these, quantum reservoir computing has emerged as a particularly powerful approach, combining quantum computation with machine learning for modeling nonlinear temporal dependencies in high-dimensional time series. As with many data-driven disciplines, quantitative finance and econometrics can hugely benefit from emerging quantum technologies. In this work, we investigate the application of quantum reservoir computing for realized volatility forecasting. Our model employs a fully connected transverse-field Ising Hamiltonian as the reservoir with distinct input and memory qubits to capture temporal dependencies. The quantum reservoir computing approach is benchmarked against several econometric models and standard machine learning algorithms. The models are evaluated using multiple error metrics and the model confidence set procedures. To enhance interpretability and mitigate current quantum hardware limitations, we utilize wrapper-based forward selection for feature selection, identifying optimal subsets, and quantifying feature importance via Shapley values. Our results indicate that the proposed quantum reservoir approach consistently outperforms benchmark models across various metrics, highlighting its potential for financial forecasting despite existing quantum hardware constraints. This work serves as a proof-of-concept for the applicability of quantum computing in econometrics and financial analysis, paving the way for further research into quantum-enhanced predictive modeling as quantum hardware capabilities continue to advance.
\end{abstract}
\maketitle
\section{Introduction}
Quantum mechanics promises to enhance computing capacity in a fundamental way ~\cite{shor1994algorithms,steane1998quantum,kitaev2002classical}. In recent years, quantum computers, albeit with severe limitations, are rapidly emerging in various physical platforms including superconducting qubits \cite{arute2019quantum,wu2021strong,zhikun2024multilevel,ren2022experimental}, ion traps \cite{kielpinski2002architecture,zhang2017observation,ringbauer2022universal,monroe2021Programmable}, Rydberg atoms~\cite{bernien2017probing,ebadi2021quantum}, photonic setups~\cite{zhong2021phase,xiao2020nonHermitian}, nitrogen vacancy centers in diamond \cite{nemoto2014photonic} and topological qubits \cite{microsoft2025interferometric}.
While near-term quantum computers have demonstrated quantum supremacy over their classical counterparts, they suffer from various imperfections such as finite coherence time and a limited number of qubits~\cite{preskill2018quantum}.
Therefore, many algorithms with proven quantum advantages cannot be implemented on noisy near-term quantum computers.
A natural emergent question is - \emph{can near-term quantum computers be used for machine learning?}
The advent of variational quantum algorithms~\cite{cerezo2021variational} and related quantum approximate optimization algorithms~\cite{farhi2014quantum} are attempts to answer this question where they have
been adapted to solving problems in a truly diverse array of disciplines, including molecular simulations \cite{cao2019quantum,li2023fermionic}, biology \cite{baiardi2023quantum}, or particle physics \cite{paulson2021simulating}.
In parallel, methods of classical time series analysis like Crutchfield computational mechanics framework \cite{crutchfield1989inferring, shalizi2001computational}, have also been shown to be improved by quantum encodings which reduce model complexity \cite{gu2012quantum}. Quantum computing methods are also starting to become important for modern quantitative finance problems like asset pricing \cite{da2023quantum}, portfolio optimization \cite{mugel2021hybrid,mugel2022dynamic}, credit sales classification \cite{wisniewska2023variational}, risk management \cite{leclerc2023financial}, loan eligibility prediction \cite{innan2024lep} among others. See Refs ~\cite{orus2019quantum,herman2022survey,naik2025portfolio} for recent surveys.
The key approach behind these algorithms is to optimize a parametrized quantum circuit to optimize the loss function in an iterative way like a classical neural network~\cite{cerezo2021variational}. However, there is a different paradigm of classical machine learning, namely \textit{reservoir computing}~\cite{jaeger2004harnessing}, which does not seek to optimize the parameters of a neural network, and training takes place only at the final output layer.While this approach to machine learning has transformed time series modeling, it still operates in the realm of classical computing~\cite{LiDeep2024,kimComprehensiveSurveyDeep2025}.
Quantum versions of this approach are particularly promising for near-term quantum computers and have been proposed with various underlying platforms like disordered spin chains \cite{fujii2017harnessing}, quantum optical energy levels \cite{govia2021quantum}, spin-boson platforms \cite{das2025quantum}, or superconducting devices \cite{yasuda2023quantum} for some mathematical and physical problems. They are exceptionally useful for time-series forecasting, as demonstrated with canonical models like Mackey-Glass \cite{mackey1977oscillation} and autoregressive moving-average (\texttt{ARMA})~\cite{whittle1951hypothesis}. Thus, it is natural to ask whether quantum reservoir computing can be useful for accurate forecasting of real-world financial time series data.
Volatility modeling is a fundamental aspect of financial econometrics, essential to understand and manage the uncertainty or risks associated with financial markets. Accurate modeling and forecasting of volatility are crucial for various applications, including risk management, portfolio optimization, and derivative pricing \citep{black1973pricing, poon2003forecasting}. Given its critical role, developing robust volatility forecasting models has been a major focus of financial research. While traditional models such as Generalized Autoregressive Conditional Heteroskedasticity (\texttt{GARCH})~\citep{bollerslev1986generalized} and its standard extensions are widely used, their reliance on low-dimensional parametric recursions for conditional variance, typically linear in squared returns, imposes strong functional restrictions. These constraints may limit their ability to flexibly represent complex nonlinear and multiscale volatility dynamics governing realized volatility at medium and long horizons. Recognizing these limitations, subsequent research has emphasized the importance of understanding the distribution of realized stock return volatility for more effective financial analysis, as highlighted in \citet{andersen2001distribution}. Furthermore, econometric analysis of realized volatility, particularly its application in estimating stochastic volatility models, has provided crucial insights into the behavior of financial time series, as demonstrated by \citet{barndorff2002econometric}. To address some of the shortcomings of traditional models, the Heterogeneous Autoregressive (\texttt{HAR}) model \citep{corsi2009simple} was developed, offering a more nuanced approach by incorporating realized volatility over different time horizons, thus providing a more accurate and comprehensive measure of market risk. Its variations, including \texttt{HAR-J} (\texttt{HAR} with jumps) and \texttt{CHAR} (continuous \texttt{HAR}) \citep{andersen2007roughing}, \texttt{SHAR} (semivariance-\texttt{HAR}) \citep{patton2015good}, and \texttt{HARQ} (\texttt{HAR} with realized quarticity) \citep{bollerslev2016exploiting}, further enhanced its ability to capture different aspects of market volatility. The nonlinear and often non-Gaussian nature of financial data has driven researchers to explore more sophisticated approaches, such as machine learning techniques to capture complex relationships that traditional econometric models may miss \citep{kuan1994artificial, habibnia2016essays, gu2020empirical, bucci2020realized,
gu2021autoencoder, habibnia2021forecasting, zhu2023forecasting, jiang2023reimaging, chen2024deep}.
In the realm of realized volatility forecasting, machine learning models have shown promising results. While some studies found machine learning models to perform similarly to or slightly worse than traditional \texttt{HAR}-family models \citep{hillebrand2010bagging, fernandes2014modeling, audrino2016lassoing, branco2022forecasting}, others reported significant improvements using machine learning approaches \citep{bucci2020realized, christensen2023machine, zhu2023forecasting}. For instance, \citet{christensen2023machine} applied machine learning to volatility forecasting, demonstrating that these models can significantly improve prediction accuracy, especially when combined with macroeconomic variables. Similarly, \citet{zhang2024volatility} explored the integration of machine learning with intraday commonality, further enhancing the precision of volatility forecasts. \cite{bucci2020realized} also demonstrated the effectiveness of Long Short-Term Memory (\texttt{LSTM}) neural networks in forecasting monthly S\&P 500 realized volatility, and have found strong links between volatility and macroeconomic factors paving the way for more advanced machine learning applications in this field. The importance of capturing the multifaceted nature of volatility is underscored by the work of \citet{ghysels2006predicting}, who emphasized the value of high-frequency data for more accurate volatility estimation. By leveraging data sampled at different frequencies, their approach offers a richer understanding of market dynamics, which is crucial for effective risk management and portfolio allocation. Recent reviews, such as the one by \citet{gunnarsson2024prediction}, have highlighted the growing role of machine learning in volatility modeling, particularly in the prediction of realized and implied volatility indices. These advancements underscore the trend toward more data-driven approaches in finance, which offer superior adaptability to the complexities of modern financial markets. Among machine learning algorithms, recurrent neural networks and their variant, \texttt{LSTM} networks, have shown particular success in modeling time series data due to their inherent ability to capture temporal dependencies. Recurrent neural networks are designed to recognize patterns in sequences of data by maintaining a hidden state that is influenced by previous inputs, making them ideal for time-dependent data such as financial time series~\citep{hochreiter1997long,goodfellow2016deep,bucci2020realized}.
This work introduces a novel approach using quantum reservoir computing for realized volatility forecasting of the S\&P 500 index. Quantum reservoir computing, an extension of classical reservoir computing, has been chosen for this study due to its ability to efficiently process temporal data on near-term quantum computers~\citep{fujii2017harnessing}, while maintaining a low training cost and providing a state space that grows exponentially with the number of qubits to capture features. Quantum reservoir computing models harness quantum features to enhance the prediction power of time series data. They utilize the dynamics of fixed quantum systems known as the ``reservoir". Unlike conventional neural networks, in reservoir computing, the parameters of the reservoir are not trained, and the learning procedure takes place at the final output layer after performing measurements~\cite {lukovsevivcius2009reservoir}. This makes reservoir computing particularly advantageous for tasks requiring the capture of temporal dynamics without the heavy computational burden associated with training traditional neural networks. Quantum reservoir computing extends this concept into the quantum domain, leveraging quantum states to enhance computational power and efficiency. This approach can be promising for volatility forecasting, where capturing complex, time-dependent relationships is critical \citep{fujii2017harnessing, mujal2023time, garcia2023scalable,llodra2025quantum, garcia2024quantum, thakkar2023improved, rivera2022time}. Such quantum reservoir platforms have already been conceptualized in various atomic and spin lattices \cite{fujii2017harnessing,settino2024memory,llodra2025quantum} as well as photonic platforms \cite{garcia2023scalable,nerenberg2024photon, suprano2024experimental}, It is of particular interest to note that very recent results on reservoir computing already hint that quantum properties may lead to fast and reliable forecasts with smaller resources \citep{abbas2024reservoir}. We aim to showcase the proof of concept and capability of quantum machine learning in real-world time series analysis, particularly in financial forecasting. We benchmark our quantum model against several classical models, including the \texttt{HAR} and \texttt{HARX} models, as well as traditional neural network algorithms. By leveraging a comprehensive set of features, including market microstructure and macroeconomic variables, we explore whether quantum nonlinear models can more effectively capture the intricate relationships governing market volatility. Prior research, such as that by~\citet{alaminos2022forecasting} and~\citet{thakkar2023improved}, has demonstrated the potential of quantum machine learning in financial forecasting, suggesting that quantum models may eventually outperform classical machine learning algorithms, particularly when dealing with large datasets and complex patterns. \\
The paper is organized as follows. Section~\ref{sec: technical_background} outlines the classical models, including traditional regression methods and machine learning models, for realized volatility forecasting.
Sec.~\ref{sec: Proposed QRC} demonstrates our method of volatility forecasting with quantum reservoir computing. If you are not familiar with quantum computing, Sect. Appendix~\ref{sec: quantum computing} gives a brief discussion of the principle of quantum computing, bridging concepts from both financial econometrics and quantum computing to facilitate interdisciplinary understanding. Sec.~\ref{sec: Empirical Results} compares our results with other volatility forecasting strategies, before concluding discussions in Sec.~\ref{sec: conclusions and outlook}.
\section{Technical Background}
\label{sec: technical_background}
This paper lies at the intersection of three different subjects, namely econometrics, machine learning, and quantum computation. In this section, we provide a brief overview of the former two subjects for physicists, introducing the concepts and methodologies that we will use later. For non-physicists, we include an appendix on quantum computation preliminaries as well as a notation table. Depending on their expertise, the readers can skip the following subsections.
\subsection{Classical Models for Realized Volatility Forecasting}
Stock market volatility is an important economic factor that reflects the risk of investment at a given time. The concept of realized volatility, formally introduced by \citet{andersen1998answering}, marked a significant advancement in volatility measurement by providing an accurate, model-free estimate utilizing high-frequency financial data. Realized volatility, \(RV_t\), at time \(t\), is defined as the square root of the sum of squared returns within a given time interval:
\begin{equation}
RV_t = \sqrt{ \sum_{i=1}^{N_t} r_{i,t}^2 }, \quad \text{and we model} \; \log(RV_t)
\end{equation}
Where \(r_{i,t}\) represents the return on the day \(i\) within period \(t\), and \(N_t\) indicates the number of observations (e.g., trading days) in that period. Autoregressive (\texttt{AR}) models and their various extensions have since become essential tools for modeling and forecasting realized volatility due to their computational simplicity and effectiveness as baseline methodologies in time series analysis. Although traditional \texttt{AR} models effectively capture short-term volatility dependencies, they often fail to reflect the long-memory characteristics typically present in volatility dynamics. To overcome this limitation, the Heterogeneous Autoregressive (\texttt{HAR}) model proposed by \citet{corsi2009simple} incorporates realized volatility at multiple time horizons, thus effectively capturing both persistence and scaling features. The general form of the \texttt{HAR} model is expressed as the following:
\begin{equation}
RV_t = \beta_0 + \beta_d RV_{t-1}^{(d)} + \beta_w RV_{t-1}^{(w)} + \beta_m RV_{t-1}^{(m)} + \varepsilon_t,
\end{equation}
where \(RV_{t-1}^{(d)}\), \(RV_{t-1}^{(w)}\), and \(RV_{t-1}^{(m)}\) represent realized volatility averages computed over daily, weekly, and monthly horizons. The coefficients \(\beta_d\), \(\beta_w\), and \(\beta_m\) measure the contributions of these short-term, intermediate-term, and long-term volatility components to the current volatility level. The parameter $\varepsilon_t$ is the residual term assumed to follow a conditionally heteroskedastic process, with $\mathbb{E}[\varepsilon_t|\mathcal{F}_{t-1}] = 0$ and $\operatorname{Var}(\varepsilon_t|\mathcal{F}_{t-1}) = \sigma_t^2$, where \(\mathcal{F}_{t-1}\) denotes the information set available at time \(t-1\). This specification allows the \texttt{HAR} model to capture volatility at multiple frequencies, accommodating heterogeneity in market participants' investment horizons. The model assumes weak stationarity of the log-volatility series and serially uncorrelated errors. Estimation is typically performed via ordinary least squares on the log-transformed realized volatility series to stabilize variance and reduce the impact of heteroskedasticity. For statistical inference, heteroskedasticity and autocorrelation-consistent (HAC) standard errors are employed following \citet{newey1986simple}, ensuring robustness to serial correlation and time-varying error variance.
To incorporate exogenous macroeconomic and financial variables, we adopt a monthly version of the \texttt{HAR} model, resulting in a \texttt{HARX} specification that captures realized volatility dynamics at multiple time scales and permits additional predictors. Specifically, we define the monthly \texttt{HARX} model as:
\begin{equation}
\resizebox{\columnwidth}{!}{$
\begin{split}
RV_t &= \beta_0 + \beta_1 RV_{t-1} + \beta_2 \left( \frac{1}{3} \sum_{i=1}^{3} RV_{t-i} \right) + \beta_3 \left( \frac{1}{12} \sum_{i=1}^{12} RV_{t-i} \right) \\
&\quad + \sum_{k=1}^k\gamma_k^TX_{t-k} + \varepsilon_t,
\end{split}
$}
\end{equation}
where \(RV_{t-1}\) represents the monthly realized volatility lagged by one month, capturing short-term dynamics; \(\frac{1}{3}\sum_{i=1}^{3} RV_{t-i}\) and \(\frac{1}{12}\sum_{i=1}^{12} RV_{t-i}\) represent the quarterly and annual realized volatility averages, respectively; and \(X_{t-k}\) denotes the vector of macroeconomic and financial features available at lag \(k\), with corresponding coefficient vector \(\gamma_k\). Residual diagnostics and model stability tests are performed to ensure robustness. The \texttt{HARX} model serves as a high-performing linear benchmark, capturing persistence at different time horizons while allowing flexible enhancement through $X_t$. \\
However, linear models such as the \texttt{HAR} may still fall short in capturing the non-linear dynamics inherent in financial time series. This limitation has motivated the development of more advanced models, such as the Realized \texttt{GARCH} model by \citet{hansen2012realized}, which integrates realized measures into the GARCH framework, and extensions of the standard \texttt{HAR} framework designed to accommodate non-linearities and regime-switching behaviors \citep{mcaleer2008realized, hillebrand2010bagging}. More recently, hybrid approaches combining \texttt{HAR} structures with advanced methodologies, such as neural networks and regularization-based procedures (e.g., LASSO or elastic net), have further enhanced forecast performance and adaptability in high-dimensional settings \citep{audrino2016lassoing}. \\
Machine learning techniques, briefly discussed later in this section, offer promising alternatives due to their ability to effectively capture complex nonlinear relationships without requiring explicit parametric functional forms. Models like Long Short-Term Memory (\texttt{LSTM}) networks and Reservoir Computing (\texttt{RC}) have demonstrated effectiveness in financial forecasting tasks \citep{fischer2018deep, butcher2013reservoir}. Incorporating exogenous variables leads to extensions such as Long Short-Term Memory with Exogenous features (\texttt{LSTMX}) and Reservoir Computing with Exogenous features (\texttt{RCX}), enhancing forecast performance by utilizing additional information.
\begin{figure*}
\centering \includegraphics[width=0.8\textwidth]{Fig1_1.png}
\caption{\textbf{Schematic of various machine-learning paradigms.} (a) Traditional feedforward neural networks, where $W_{in}$, $\{W_l\}$ and $W_{out}$ are the trainable weight matrices for the input layer, hidden layer(s), and output layer, respectively. Noticeably, the information flows unidirectionally layer by layer from the input to the output.
(b) \texttt{LSTM} networks are a special form of the recurrent neural networks architecture. \texttt{LSTM} introduces two key states as cell state~$\boldsymbol{c}_i$ for long-term memory and hidden state~$\boldsymbol{h}_i$-for short-term memory. In addition, \texttt{LSTM} layers utilize additional gates to control which information in the hidden state is output and passed to the next hidden state. These additional gates overcome the common issues recurrent neural networks face in learning long-term dependencies.
(c) Reservoir Computing is a computational framework that inputs $\boldsymbol{x}_t$ at each step are mapped to a high-dimensional space by
a fixed random $W_{in}$. The hidden states $\boldsymbol{h}_t$ evolve dynamically through a fixed random reservoir $W_r$, which models the system's temporal denpendencies. Notably, both $W_{in}$ and $W_r$ are randomly initialized and reamin fixed throughout the process, while only the output weights matrix $W_{out}$ is trainable. This feature make reservoir computing efficient for handling time series and dynamic system tasks, requiring minimal training effort.}
\label{Fig: schematic_different_computational_models}
\end{figure*}
\subsection{Machine Learning Preliminaries}
Machine learning involves the development of algorithms and statistical models that allow systems to perform specific tasks effectively by analyzing data, identifying patterns, and making predictions.
The artificial neural network~\cite{zou2009overview}, as one of the most powerful algorithms in machine learning, is widely used in many domains, including image recognition, natural language processing, recommendation systems, predictive analytics, and time series processes.
\subsubsection{Feedforward Neural Networks}
Feedforward neural networks are the simplest type of artificial neural network, consisting of an input layer, an output layer, and one or more hidden layers that connect the input layer to output layer~(see Fig.~\ref{Fig: schematic_different_computational_models}(a)).
The term "feedforward`` indicates that the architecture of these networks relies on transforming inputs into outputs through a series of operations, where each operation involves multiplying the input by a weight matrix and applying an activation function.
Specifically, consider the output of the $l$-th layer, denoted as $\boldsymbol{h}_l$, such that the input of the network is represented as $\boldsymbol{h}_0=\boldsymbol{x}$, and the final output is $\boldsymbol{h}_L=\boldsymbol{\hat{y}}$.
For the $l$-th layer in the Feedforward neural networks, the output takes the form
\begin{equation}
\boldsymbol{h}_l = \mathcal{F}(W_l\boldsymbol{h}_{l-1}+\boldsymbol{b}_l),
\end{equation}
where $W_l$ is the weight matrix, $\boldsymbol{b}_l$ is the bias vector, and $\mathcal{F}(\cdot)$ is a certain activation function.
The overall mathematical model of an feedforward neural networks is a nested composition of these computations, where the output $\boldsymbol{h}_l$ of each layer serves as the input of the next layer.
Because there are no cycles or loops in the networks, the information flows strictly in a single direction: from the input node, through the hidden nodes, and to the output nodes~(see Fig.~\ref{Fig: schematic_different_computational_models}(a)).
Training for a neural network is based on a set of labeled data $\{ (\boldsymbol{x}_i,\boldsymbol{y}_i) \}$, where $\boldsymbol{x}_i$ and $\boldsymbol{y}_i$ represent input and output, respectively.
We train the weight matrices $\{W_l\}_l$ as well as the bias vectors $\{\boldsymbol{b}_l\}_l$ so that for any input data $\boldsymbol{x}_i$, the corresponding output of the neural network $\boldsymbol{\hat{y}}_i$ closely approximates the real label $\boldsymbol{y}_i$. This is accomplished by minimizing a loss function such as
\begin{equation}
L=\min_{\{W_l,\boldsymbol{b}_l\}} \frac{1}{n}\sum_i^n (\boldsymbol{y}_i-\boldsymbol{\hat{y}}_i)^2,
\label{cost_fun}
\end{equation}
where $n$ is the size of the training set.
Feedforward neural networks are well-suited for tasks where each input is independent of the others, such as image classification or regression. However, they are not ideal for tasks that require capturing temporal or sequential dependencies.
\subsubsection{Long Short-Term Memory Neural Networks}
For addressing tasks with spatio-temporal dependencies, it is often necessary to utilize information from past data to make accurate predictions about future results. For example, when reading an article, the meaning of each word is interpreted based on the understanding of the preceding words.
Long Short-Term Memory (\texttt{LSTM}) networks~\cite{hochreiter1997long} address this problem by introducing a specialized architecture designed to capture both short-term and long-term dependencies in sequential data. Unlike standard recurrent neural networks, which rely solely on a hidden state $\boldsymbol{h}_t$ to carry forward information, \texttt{LSTM} introduce an additional cell state $\boldsymbol{c}_t$, which serves as long-term memory. The cell state enables \texttt{LSTM} to selectively retain or discard information over time.
At each time step, as Fig.~\ref{Fig: schematic_different_computational_models}(b) shows, the \texttt{LSTM} cell receives the following inputs:
\begin{enumerate}
\item Hidden state $\boldsymbol{h}_{t-1}$: The short-term memory from the previous time step.
\item Cell state $\boldsymbol{c}_{t-1}$: The long-term memory from the previous time step.
\item Current input $\boldsymbol{x}_t$: The input in the current time step.
\end{enumerate}
The operation of the \texttt{LSTM} cell can be described in terms of three different gates: (i) forget gate; (ii) input gate; and (iii) output gate (see Fig.~\ref{Fig: schematic_different_computational_models}(b)).
The forget gate determines which parts of the previous cell state $\boldsymbol{c}_{t-1}$ should be ``forgotten" or retained. It uses a sigmoid activation function to produce values between 0~(completely forget) and 1~(completely keep) for each element, defined as
\begin{equation}
\boldsymbol{f}_t = \text{sigmoid}(W_f[\boldsymbol{h}_{t-1},\boldsymbol{x}_t]+\boldsymbol{b}_f),
\end{equation}
where $W_f$ is the weight matrix for the forget gate, $\boldsymbol{b}_f$ is the bias term, and $\boldsymbol{f}_t$ is a vector whose elements take values between $0$ and $1$.
The input gate determines which parts of the current input $\boldsymbol{x}_t$ should be added to the cell state $\boldsymbol{c}_t$.
It has two components, a sigmoid layer to decide which values to update, and a tanh layer to create new candidate values to potentially add to the cell state $\boldsymbol{c}_t$, which formula are defined as
\begin{equation}
\begin{split}
\boldsymbol{i}_t &= \text{sigmoid}(W_i[\boldsymbol{h}_{t-1},\boldsymbol{x}_t]+\boldsymbol{b}_i)\\
\tilde{\boldsymbol{c}}_t &= \text{tanh}(W_c[\boldsymbol{h}_{t-1},\boldsymbol{x}_t]+\boldsymbol{b}_c)
\end{split}
\end{equation}
where $\boldsymbol{i}_t$ is the input gate output, a vector of values between $0$ and $1$, and $\tilde{\boldsymbol{c}}_t$ is the candidate cell state.
The information from the forget gate and the input gate are used to update the cell state $\boldsymbol{c}_t$, as
\begin{equation}
\boldsymbol{c}_t = \boldsymbol{f}_t\odot \boldsymbol{c}_{t-1}+\boldsymbol{i}_t\odot \tilde{\boldsymbol{c}}_t,
\end{equation}
where $\odot$ is the element-wise multiplication operator.
As a result, $\boldsymbol{c}_t$ contains a nonlinear combination of the input state $\boldsymbol{x}_t$ and the previous input data.
The output gate determines the vector $\boldsymbol{o}_t$ which is a nonlinear function of the the hidden state $\boldsymbol{h}_{t-1}$ as well as the input $\boldsymbol{x}_t$ as
\begin{equation}
\boldsymbol{o}_t =\text{sigmoid}(W_o[\boldsymbol{h}_{t-1},\boldsymbol{x}_t]+\boldsymbol{b}_o).
\end{equation}
The next hidden state $\boldsymbol{h}_t$ is then determined through combination of $\boldsymbol{o}_t$ and the cell state $\boldsymbol{c}_t$ as
\begin{equation}
\boldsymbol{h}_t =\boldsymbol{o}_t \odot \text{tanh}(\boldsymbol{c}_t).
\end{equation}
In general, the output $\hat{\boldsymbol{y}}_t$ is a function of $h_t$ as $\hat{\boldsymbol{y}}_t=g(h_t)$, where $g(\cdot)$ depends on the specific problem.
Training is carried out by optimizing a loss function, of the form of Eq.~(\ref{cost_fun}), which aims to minimize the difference between $\hat{\boldsymbol{y}}_t$ and the real output $\boldsymbol{y}_t$.
During the training procedure, the weight matrices $\{W_f,W_i,W_c,W_o\}$ and bias vectors $\{\boldsymbol{b}_f,\boldsymbol{b}_i,\boldsymbol{b}_c,\boldsymbol{b}_o\}$ are all trained and updated iteratively, see Fig.~\ref{Fig: schematic_different_computational_models}(b).
The \texttt{LSTM} architecture uses gating mechanisms, see Fig.\ref{Fig: schematic_different_computational_models}(b), to regulate information flow and gradient propagation~\cite{hochreiter1997long}, effectively mitigating vanishing and exploding gradients common in standard recurrent neural networks and enabling the retention of long-term information.
These characteristics make \texttt{LSTM} widely used in fields such as natural language processing and time-series prediction.
\subsubsection{Classical Reservoir Computing}
In the early 2000s, echo state networks~\cite{jaeger2004harnessing} and liquid state machine~\cite{maass2002real} were independently proposed as the seminal approach of the time series model, which are grouped into the framework of classical reservoir computing.
Reservoir computing models share the common principle of using a ``reservoir'', comprising $W_r$~(the recurrent weight matrix) and $W_{in}$ (the input-to-reservoir weight matrix), to project inputs into a high-dimensional feature space.
This transformation enables the effective capture of complex patterns, relationships, and temporal dynamics in the data.
Unlike conventional neural networks, reservoir computing does not train weight matrices $W_r$ and $W_{in}$ and bias vectors $\boldsymbol{b}$.
Instead, these parameters are randomly initialized and remain fixed during training.
A simple and trainable read-out mechanism, i.e. linear regression, is then used to generate information for this high-dimensional representation.
In the whole reservoir computing process, only the weights $W_{out}$ in the linear regression layer are trained, significantly reducing training cost.
The mathematical representations of these models can be expressed as:
\begin{equation}
\begin{split}
h_t&=(1-\alpha)\boldsymbol{h}_{t-1}+ \mathcal{F}(W_r\boldsymbol{h}_{t-1}+W_{in}\boldsymbol{x}_t+\boldsymbol{b})\\
y_t&=W_{out}\boldsymbol{h}_t,
\end{split}\label{eq:classical_reservoir_learning}
\end{equation}
where $\boldsymbol{h}_t$ represents the hidden state at time $t$, and $\alpha$ is the leak rate, a hyperparameter controlling the update speed of hidden state, see Fig.~\ref{Fig: schematic_different_computational_models}(c).
Since only $W_{out}$ needs to be trained, reservoir computing is highly suitable for deployment in artificial or natural physical systems characterized by high-dimensionality and nonlinear transformations.
To date, reservoir computing has been widely implemented in various systems, such as cellular automata~\cite{yilmazSymbolicComputationUsing2015,mcdonaldReservoirComputingExtreme2017}, coupled oscillators~\cite{yamaneWaveBasedReservoirComputing2015}, analog circuits~\cite{royLiquidStateMachine2014,katumbaLowLossPhotonicReservoir2018}, optical node arrays~\cite{duportAllopticalReservoirComputing2012,mesaritakisMicroRingResonators2013,dejonckheereAllopticalReservoirComputer2014} and biological organization~\cite{goudarziDNAReservoirComputing2013}.
\section{Quantum Reservoir Computing for Realized Volatility Forecasting}
\label{sec: Proposed QRC}
\begin{figure*}
\centering
\includegraphics[width=0.8\textwidth]{Fig2.png}
\caption{\textbf{Quantum reservoir computing schematic.} The task is to forecast $\widehat{RV}_t$ based on the past information of features $\boldsymbol{x}$ at discrete time from $t{-}3$ to $t{-}1$. The Quantum Reservoir consists of two subsystems, the input qubits $\rho_I$ for inputting features and the hidden qubits $\rho_h$ for hidden state, composed of $n_1$ and $n_2$ qubits respectively. In the first, the feature variables $\boldsymbol{x}_{t-3}$ are encoded into the input qubits via single-qubit gates $R_Y(\boldsymbol{x}_{t-3})$ applied on $|0\rangle^{\otimes n_1}$. Therefore the input density matrix is given as $\rho_I(\boldsymbol{x}_{t-3})$. The hidden state at first is initialized as $\rho_h(t{-}3)=|0\rangle\langle 0|^{\otimes n_2}$.
Thus the quantum reservoir state is given by $\rho_I (\boldsymbol{x}_{t-3}) \otimes \rho_{h}(t{-}3)$. The state evolves under the action of the Hamiltonian $H$ for time $\tau$, described by $e^{-iH\tau}$. To update the hidden state for the next step, the input $n_1$ qubits are discarded such that the new hidden state is given by a partial trace as $\rho_h(t{-}2)$.
This process continues until all features from $t{-}3$ to $t{-}1$ have been encoded within the quantum reservoir. After processing all features, measurements on Pauli-Z basis of all qubits are performed to read out useful information from the quantum reservoir, yielding the expectation value $\langle Z_j\rangle_\tau$ for the $j$-th qubit. These expectation values collectively form the readout vector $\boldsymbol{m}_t^{(1)}=[\langle Z_1\rangle_\tau,\cdots,\langle Z_{n_1+n_2}\rangle_\tau]$. A linear regression model is used to predict $\widehat{RV}_t$ based on this readout vector as $\widehat{RV}_t=W_{out}\boldsymbol{m}_t^{(1)}$, where $W_{out}$ is a weight matrix trained using the ridge regression method.}\label{Fig: schematic_qrl}
\end{figure*}
By harnessing the unique properties of quantum mechanics, such as superposition and entanglement, quantum computing demonstrates a quantum advantage in solving specific problems, including integer factorization~\cite{shor1999polynomial}, random circuit sampling~\cite{arute2019quantum} and quantum simulation~\cite{sethUniversal1996}.
These achievements identify quantum computing as a powerful method for information processing and thus drive researchers to explore more fields that might benefit from quantum computing.
A natural extension of the quantum computing application is the investigation of reservoir computing with quantum systems.
Fujii and Nakajima took the first theoretical step by proposing disordered quantum spin ensembles as reservoirs \cite{fujii2017harnessing}, which has since been extended to photonic \cite{ghosh2019quantum,garcia2023scalable,nerenberg2024photon}, non-linear oscillator~\cite{govia2021quantum}, and neutral atomic Rydberg array~\cite{bravo2022quantum} platforms. In terms of physical implementation, nuclear-spin-based reservoirs~\cite{negoro2018machine} as well as superconducting qubit platforms such as IBM \cite{chen2020temporal} have already been successful in demonstrating various aspects of reservoir computing.
See Ref.~\cite{mujal2021opportunities} for a more detailed recent overview.
Here, we introduce a quantum reservoir computing approach, see Fig.~\ref{Fig: schematic_qrl}, based on a qubit system driven by a fully connected transverse-field Ising model as a reservoir.
The quantum reservoir consists of two subsystems: (i) input qubits, whose quantum state is shown as $\rho_I$, for encoding variables; and (ii) hidden qubits, whose quantum state is shown as $\rho_h$, to store past information to predict future outcomes.
The whole reservoir is controlled by full connected Hamiltonian $H$ which is given by
\begin{equation}\label{eq:Hamiltonian_Ising}
H=\sum_{ij}J_{ij}X_iX_j + v\sum_iZ_i
\end{equation}
where $X_i$ and $Z_i$ are the Pauli operators acting on $i$-th qubits (see Table~\ref{table:notation_table_adjusted} for the definition of Pauli matrices), $v$ is the strength of the magnetic field which is used as the unit of energy and thus set to be $v=1$ and finally
$J_{ij}$ are exchange couplings between qubits $i$ and $j$ which are randomly sampled from $J_{ij}/v\in [0,1]$.
After being randomly initialized, the Hamiltonian $H$ is fixed throughout the process, similar to a classical reservoir.
As shown in the leftmost part of Fig.~\ref{Fig: schematic_qrl}, the state of the ensemble quantum system is initialized as $|0\rangle$ at the beginning of the quantum reservoir computing process:
\begin{equation}
\rho_I\otimes\rho_h=|0\rangle\langle0|^{\otimes n_1}\otimes|0\rangle\langle0|^{\otimes n_2},
\end{equation}
where $n_1$ and $n_2$ are the respective sizes of the input and hidden subsystems, which can be adjusted for different tasks. In this paper, we consider $n_1+n_2=10$, which is more accessible size on current quantum computers.
The goal in this paper is to predict $RV_t$ by using past $k$ steps features, i.e., the time series $\boldsymbol{x}_{t-k}, \cdots,\boldsymbol{x}_{t-1}\}$.
For convenience of notation, we denote $\boldsymbol{x}_{t-k} = [x_{t-k,1}, x_{t-k,2},\cdots,x_{t-k,n_1}]^T$, where each element $x_{t-k, i}$ represents the $i$-th economic feature at time step $t{-}k$. For example the input vector $\boldsymbol{x}_{t-k}$ may contain three input data at time $t-k$ such as realized volatility $RV_{t-k}$, Dividend Yield ($DP_{t-k}$) ratio and Earning-Price ($EP_{t-k}$) ratio~(note that all external features used in this paper are defined in Table.~\ref{Table : Macro_features}). In this situation, $\boldsymbol{x}_{t-k}=[RV_{t-k},DP_{t-k}, EP_{t-k}]^T$ contains three input features.
In this work, we limit $k{=}3$, that is, we consider the memory depth for learning to be only three steps.
Although generalizing to other memory depths is straightforward and it may appear that considering large memories may yield better results, one must also consider the fact that larger memory depths require deeper quantum circuits.
Given the finite coherence times available in near-term quantum devices, this trade-off quickly becomes considerable.
As we shall see, our choice of $k{=}3$ already yields a satisfactory predictive accuracy.
Note that all features are scaled between $[-\pi,+\pi]$, fulfilling the requirement of the $R_Y$ gate. The input vectors $\{ \boldsymbol{x}_{t-3}, \boldsymbol{x}_{t-2}, \boldsymbol{x}_{t-1}\}$ are fed into the reservoir model sequentially.
The reservoir performs an iterative loop consisting of three steps: \\
\noindent \textbf{Step I: Encoding input data $\boldsymbol{x}_{t-3}$}. The classical variables $\boldsymbol{x}_{t-3}$, assumed to have $n_1$ features, are encoded to the input qubits of the quantum reservoir through phase encoding, i.e., single-qubit rotation gates around the $Y$ direction, such that the input state is given by
\begin{eqnarray}
|\psi_I (\boldsymbol{x}_{t-3}) \rangle &=&R_Y(\boldsymbol{x}_{t-3}) |0\rangle^{\otimes n_1}, \cr
\text{where } \quad R_Y(\boldsymbol{x}_{t-3})&=&\bigotimes_{j=1}^{n_1} e^{-i\frac{x_{t-3,j}}{2}Y_j}.
\end{eqnarray}
Rotation encoding is a common choice for near-term quantum hardware because it is easy to implement and embeds classical features into quantum states via parameterized single-qubit rotations, thereby effectively leveraging the continuous degrees of freedom of a qubit state.
Therefore, the input density matrix is given by $\rho_I (\boldsymbol{x}_{t-3}) = |\psi_I (\boldsymbol{x}_{t-3}) \rangle\langle \psi_I (\boldsymbol{x}_{t-3})|$. The hidden $n_2$ qubits (the total reservoir thus consists of $n_1+n_2$ qubits) at this stage are initialized as
\begin{equation}
\rho_{h}(t-3) = \big[|0\rangle\langle0|\big]^{\otimes n_2}.
\end{equation}
Thus, the collective input state of the quantum reservoir state is given by $\rho_I (\boldsymbol{x}_{t-3}) \otimes \rho_{h}(t{-}3)$. Then the whole reservoir evolves under the action of the Hamiltonian $H$ for a specific time $\tau$, described by $U = e^{-iH\tau}$. The evolution scrambles the input data $\boldsymbol{x}_{t-3}$ across the entire reservoir. Note that since the reservoir is a fully connected graph we set $\tau=1/v$ which is enough to guarantee that the information is scrambled across the whole system.
Following this, while the original $n_1$ qubits are discarded, the quantum state of the hidden $n_2$ qubits carry this information forward to the next step. The quantum state of the hidden qubits after discarding the input qubits now becomes
\begin{equation}
\rho_h (t-2) = \text{Tr}_{I} \left[ U(\tau) \left[ \rho_I (\boldsymbol{x}_{t-3}) \otimes \rho_{h}(t-3) \right] U^{\dagger}(\tau) \right]
\end{equation}
where $\text{Tr}_I[\cdot]$ represents partial trace over all the $n_1$ input qubits. \\
\noindent \textbf{Step II: Encoding input data $\boldsymbol{x}_{t-2}$}. After discarding the input $n_1$ qubits at the end of the last round, we replace them with fresh $n_1$ qubits with the encoding
\begin{eqnarray}
|\psi_I (\boldsymbol{x}_{t-2}) \rangle &=&R_Y(\boldsymbol{x}_{t-2}) |0\rangle^{\otimes n_1}, \cr
\text{where } \quad R_Y(\boldsymbol{x}_{t-2})&=&\bigotimes_{j=1}^{n_1} e^{-i\frac{x_{t-2,j}}{2}Y_j}.
\end{eqnarray}
Therefore the new state of the whole reservoir becomes $\rho_I (\boldsymbol{x}_{t-2}) \otimes \rho_{h}(t{-}2)$. Similar to the previous step the whole system undergoes another evolution for time $\tau$ which scrambles the input data $\boldsymbol{x}_{t-2}$ across the reservoir. Once again the $n_1$ input qubits are discarded and the quantum state of the hidden qubits is naturally updated to
\begin{equation}
\rho_h (t-1) = \text{Tr}_{I} \left[ U(\tau) \left[ \rho_I (\boldsymbol{x}_{t-2}) \otimes \rho_{h}(t-2) \right] U^{\dagger}(\tau) \right]
\end{equation}\\
\noindent \textbf{Step III: Encoding input data $\boldsymbol{x}_{t-1}$ and final measurement}. Similar to the previous step, the $n_1$ input qubits are replaced by fresh qubits to encode the input data $\boldsymbol{x}_{t-1}$ as
\begin{eqnarray}
|\psi_I (\boldsymbol{x}_{t-1}) \rangle &=&R_Y(\boldsymbol{x}_{t-1}) |0\rangle^{\otimes n_1}, \cr
\text{where } \quad R_Y(\boldsymbol{x}_{t-1})&=&\bigotimes_{j=1}^{n_1} e^{-i\frac{x_{t-1,j}}{2}Y_j}.
\end{eqnarray}
The new state of the whole reservoir thus becomes $\rho_I (\boldsymbol{x}_{t-1}) \otimes \rho_{h}(t{-}1)$. Again the whole system evolves freely for time $\tau$ which scrambles the input data $\boldsymbol{x}_{t-1}$ across the whole reservoir. At this stage, no qubit is discarded and all of them are measured in the Pauli $Z$ basis. The expectation of measurement outcomes of qubit $j$ is given by
\begin{equation}
\langle Z_j \rangle_\tau = \text{Tr} \left[ U(\tau) \left[ \rho_I (\boldsymbol{x}_{t-1}) \otimes \rho_{h}(t-1) \right] U^{\dagger}(\tau) Z_j \right] .
\end{equation}
At any given step $t$, the measurement outcomes on all the qubits form a vector $\boldsymbol{m}_t^{(1)} = [{\langle Z_1\rangle_\tau,\ldots,\langle Z_{n_1+n_2}\rangle_\tau}]^T$. The training is performed on these measured data, just as the classical reservoir learning. This is accomplished by a linear regression which is used to approximate $RV_t$ by
\begin{equation}
\widehat{RV}_t=W_{out}\boldsymbol{m}_t^{(1)},
\end{equation}
where $W_{out}$ is a weight matrix which has to be trained. If we consider the loss function as mean squared error, defined as:
\begin{equation}
\text{MSE}=\frac{1}{T}\sum_{t=1}^T(RV_t-\widehat{RV}_t)^2
\label{MSE_loss}
\end{equation}
Then the weight matrix takes an analytical form given by the ridge regression as
\begin{equation}\label{eq:training_Wout}
W_{out}=RV^\top M_1^\top {(M_1M_1^\top+\delta I)}^{-1},
\end{equation}
where $M_1=[\cdots,\boldsymbol{m}_{t}^{(1)},\cdots,\boldsymbol{m}_T^{(1)}]$, $I$ is the identity matrix and $\delta{=}10^{-8}$ is a small number which is used to guarantee that the matrix is non-singular.
We call this model as Quantum Reservoir 1 (\texttt{QR1}), since here we only use one quantum reservoir.{\color{black}~see Appendix~\ref{Numerical} for more information.}
\begin{figure}
\centering
\includegraphics[width=0.9\linewidth]{Fig3.png}
\caption{\textbf{Ensemble Reservoir}. The ensemble reservoir approach repeats the procedure twice, with the final step running the quantum dynamical evolution for times $\tau$ and $\tau/2$ in two separate ensembles respectively. This is referred to as \texttt{QR2} in this work, contrasted with \texttt{QR1} not employing ensemble learning. }
\label{Fig: virtual_node}
\end{figure}
Furthermore, to obtain richer information from the quantum reservoir, we adopt the ensemble reservoir approach~\cite{fujii2017harnessing}.
In this approach, at any given step $t$ apart from the measurement outcomes $\boldsymbol{m}^{(1)}_t$, one can use another reservoir which is almost identical to the previous setup except that the evolution at the step III is run for a time duration of $\tau/2$, instead of $\tau$, see Fig.~\ref{Fig: virtual_node}.
The ensemble of these two reservoir are combined to make a larger vector $\boldsymbol{m}_t^{(2)} = [\langle Z_1\rangle_\tau, \ldots,\langle Z_{n_1+n_2}\rangle_\tau,\langle Z_1\rangle_{\tau/2},\ldots,\langle Z_{n_1+n_2}\rangle_{\tau/2}]^T$, where $\langle Z_j\rangle_{\tau/2}$ represents the measurement of the Pauli operator $Z_j$ in the second reservoir setup in which the last evolution is run for time $\tau/2$, as schematically shown in Fig.~\ref{Fig: virtual_node}. We call this model Quantum Reservoir 2 (\texttt{QR2}) which is expected to be more precise than the \texttt{QR1} as it uses twice as much resources.
{\color{black}It is worth noting that, in the simulation stage, we directly compute the measurement outcomes of the quantum system. In an actual experiment, however, one must perform many repeated measurements (shots) and use their average as the output.
Recent studies~\cite{xiongFundamentalAspectsQuantum2025} have shown that, as the system dimension increases, QRC may suffer from an exponential concentration of measurement results toward a value that is independent of the input. Concretely, the variance of measurement outcomes across different input variables can decay exponentially to zero. As a result, an exponential number of shots would be required to reliably distinguish the outputs corresponding to different inputs, which undermines the potential advantage of QRC.
In our approach, by contrast, the QRC system is kept at a fixed size~(10 qubits), so this issue does not arise~(see Appendix~\ref{Exponential} for more information). The measurement outputs corresponding to different inputs still exhibit sufficiently large variance, indicating that our scheme does not require an excessive number of measurements to obtain accurate output estimates.}
\section{Empirical Analysis}
\label{sec: Empirical Results}
In this section, we present the empirical findings of our study, with a focus on the performance of quantum algorithms in forecasting the realized volatility of the S\&P 500 index. We begin by detailing the data set and the variables used, which include monthly realized volatility, as well as a set of macroeconomic and financial features. Following this, we outline the fitting procedure for competing models, both classical and quantum, and describe the training and evaluation methodologies applied.
A key technical contribution of this study lies in our approach to feature selection within the quantum framework. To optimize the set of features for the quantum reservoir computing models, we employ a forward selection method, a wrapper-based approach that incrementally identifies the most impactful features. This process enables the model to systematically build an optimal feature set by evaluating the performance impact of each addition of variables. Furthermore, to assess the relative importance of each selected feature, we utilize the Shapley value, a method grounded in game theory that fairly distributes the contribution of each feature across different model configurations. The use of Shapley values provides valuable information on the explanatory power of individual features, enhancing the interpretability of the forecast results of the quantum reservoir computing model.
Finally, we analyze the predictive performance of the models, comparing them based on accuracy metrics and statistical tests. This comprehensive evaluation demonstrates the effectiveness of the quantum reservoir computing approach in capturing the complex and nonlinear dynamics of financial market volatility, underscoring its potential advantages over classical models in volatility forecasting.
\begin{table*}[hbt]
\caption{\textbf{Economic variables used in this study.} Unless otherwise stated, equity data refers to S\&P 500 index.}\label{Table : Macro_features}
\centering
\begin{tabularx}{\textwidth}{l c l}
\hline
Variable & Symbol & Description \\
\hline
Realized Volatility & RV & \begin{tabular}[c]{@{}c@{}}The natural log of the square root of the sum of squared daily returns for a month.\\ RVq and RVa represent the quarterly and annual realized volatility averages.\end{tabular} \\
\begin{tabular}[c]{@{}c@{}}Dividend Yield Ratio \end{tabular} & DP & \begin{tabular}[c]{@{}c@{}}Dividends over the past year relative to current market prices\end{tabular} \\
\begin{tabular}[c]{@{}c@{}}Earning-Price Ratio \end{tabular} & EP & \begin{tabular}[c]{@{}c@{}}Earning over the past year relative to current market prices index\end{tabular} \\
Market Excess Return & MKT & Fama–French’s market factor: return of U.S. stock index minus the one-month T-bill rate\\
Value Factor & HML & Fama–French’s HML factor: average return on value stocks minus growth stocks\\
Size Premium Factor & SMB & Fama–French’s SMB factor: average return on small-cap stocks minus large-cap stocks\\
\begin{tabular}[c]{@{}c@{}}Short-Term Reversal Factor\end{tabular} & STR & \begin{tabular}[c]{@{}c@{}}Fama–French’s STR: average return on stocks with low prior minus high prior returns\\\end{tabular} \\
T-bill Rate & TB & Three-month T-bill rate\\
Monthly Inflation & INF & US monthly inflation rate\\
Default Spread & DEF & Estimate the credit risk\\
\begin{tabular}[c]{@{}c@{}}Monthly Industrial Production Growth Rate\end{tabular} & IP & Monthly growth rate of U.S. industrial production\\ \hline
\end{tabularx}
\end{table*}
\subsection{Data}
\label{Data_sec}
For this study, we utilize a dataset comprising monthly observations of Realized Volatility (RV) for the S\&P 500 index, covering the period from February 1950 to December 2017, resulting in a total of 815 data points. This approach aligns with that of \citet{bucci2020realized}, who employed a similar time frame to assess the performance of neural network models in forecasting realized volatility. Our primary variable of interest, realized volatility ($RV_t$), provides a foundational measure of market uncertainty. As illustrated in Figure~\ref{fig:RV}, realized volatility exhibits significant variability over time, reflecting alternating phases of market turbulence and stability. Capturing these dynamics accurately is critical for financial decision-making and risk management.
\begin{figure}[h]
\centering
\includegraphics[width=1\linewidth]{Fig4.png}
\caption{\textbf{Time series of the logarithm of monthly realized volatility data for the S\&P 500 Index.} Data is taken from Yahoo Finance's API~\cite{Yfinance} between January 1950 to December 2017.}
\label{fig:RV}
\end{figure}
To facilitate robust comparison, we align our study with the data and sample period used by \citet{bucci2020realized}. This alignment enables a direct and fair comparison between our quantum reservoir computing approach and established neural network models while maintaining consistency in data characteristics and market environments.
Beyond realized volatility, we include a set of macroeconomic and financial variables known to capture important economic forces and market behavior. Research has shown that adding these features can significantly enhance the performance of volatility forecasting models. For example, \citet{schwert1989} and \citet{engle2009economic} demonstrate that indicators such as Inflation Rates (INF) and Industrial Production Growth (IP) provide an essential context about the state of the broader economy, improving a model’s reliability in forecasting financial volatility. Similarly, valuation measures such as the Dividend-Price Ratio (DP) and the Earnings-Price ratio (EP) reflect market expectations and investor sentiment, providing valuable predictive power by signaling changes in expected returns and market risk \citep{welch2008comprehensive}. Additionally, market-based factors derived from the Fama-French model \citep{fama1993common,fama1996multifactor}, including Market Excess Returns (MKT), Value Factor (HML), Size Premium Factor (SMB) and Short-Term Reversal Factor (STR) capture critical aspects of equity market dynamics, offering insights into investor behavior and systematic market risk. The inclusion of financial indicators such as the Three-month Treasury Bill Rate (TB) and Default Spread (DEF) further enhances the predictive framework by incorporating measures of interest rate and credit risk conditions.
Including these diverse factors aligns with the findings of \citet{bucci2020realized}, \citet{christensen2023machine}, and \citet{zhang2024volatility}, who showed that combining macroeconomic and financial variables with machine learning models significantly enhances both explanatory power and forecasting accuracy. By incorporating this set of features into our quantum reservoir computing framework, we aim to capture the multifaceted drivers of market volatility. This approach is expected to not only improve forecast precision but also yield deeper insights into the complex interplay between financial markets and economic fundamentals.
Table~\ref{Table : Macro_features} provides a detailed summary of the features utilized in our study, consistent with the variables commonly employed in prior research.
\subsection{Training and Benchmarking}
We employ our quantum reservoir models \texttt{QR1} and \texttt{QR2} to predict realized volatility. To assess their performance, we conduct a comprehensive benchmarking study against established classical models. In particular, we compare quantum reservoir computing models \texttt{QR1} and \texttt{QR2} with the classical linear models \texttt{AR1} (which is \texttt{AR} with memory of one previous step), \texttt{AR3} (which is \texttt{AR} with memory of three previous steps), \texttt{ARMAX}, \texttt{HAR}, \texttt{HARX} as well as the non-linear machine learning based methods \texttt{LSTM}, \texttt{LSTMX}, \texttt{RC} and \texttt{RCX}, see the relevant parts of section \ref{sec: technical_background} for a short introduction to these models. For all models, including classical models and quantum reservoir computing, we employ a rolling-window approach to estimate and train the models, optimize their parameters, and perform one-step-ahead and five-steps-ahead out-of-sample predictions. Rolling re-estimation at each monthly step is consistent with standard out-of-sample evaluation practices in the volatility forecasting literature, ensuring alignment with recent information and accommodating structural change \citep{Feng2024,PattonSheppard2015,bucci2020realized,PesaranTimmermann2007}. This design promotes comparability across models by isolating learning effects from differences in estimation windows and is particularly important for monthly realized volatility, which reflects aggregated market activity and is influenced by slowly evolving macroeconomic regimes. Without regular updating, parameter estimates may become obsolete in the presence of structural breaks. Our approach, therefore, maintains regime adaptiveness while avoiding excessive sensitivity to high-frequency noise. While our study uses monthly updates to maximize predictive accuracy for benchmarking, a lower re-estimation frequency could be adopted in industrial settings to enhance stability and ease of debugging.
The initial training window spans February 1950 to June 1997 (approximately 571 months). After training and optimizing these models on the initial window, we generate a forecast for the next out-of-sample month (August 1997). Subsequently, the rolling window advances by one month, using data from March 1950 to August 1997 to predict September 1997. This process continues iteratively, rolling through all 245 out-of-sample observations, spanning from August 1997 to December 2017.
During this process, all models are re-estimated at each step to ensure optimal performance for the new rolling window.
Among the classical linear models that we use, \texttt{AR1}, \texttt{AR3}, \texttt{ARMAX}, \texttt{HAR}, \texttt{HARX} are estimated using ordinary least squares or maximum likelihood methods, depending on the structure of the model.
In addition, the machine learning based \texttt{LSTMX} model is implemented with two layers of \texttt{LSTM} cells. The hidden state size is set to 60 for the \texttt{LSTM} model and 50 for the \texttt{LSTMX} variant. After that, a linear regression is used to transform the output of LSTM(X) into the prediction of $\widehat{RV}_t$. The model parameters are optimized using the ADAM optimizer, with a learning rate of $0.001$, the batch size of $64$, and $100$ epochs.
For the classical reservoir computing models (\texttt{RC} and \texttt{RCX}), the reservoir consists of 50 hidden neurons for \texttt{RC} and 20 hidden neurons for \texttt{RCX}. The leak rate, which controls the speed at which the state of the reservoir evolves, is set at 0.6. The spectral radius, which influences the stability and non-linearity of the reservoir, is fixed at 0.9. The input scaling, which is a coefficient applied on $W_{in}$, is set at $0.1$. The readout matrix $W_{out}$, responsible for mapping the reservoir states to the output, is estimated using ridge regression to mitigate overfitting and improve generalization.
All machine learning hyperparameters, including architecture choices and training settings, are selected using time-series-aware cross-validation procedures within the training window to ensure generalization and robustness.{\color{black}~(see Appendix~\ref{Hyperparameters} for more informaiotn)}
For quantum reservoir computing, taking into account the current limitations of quantum devices, we used a quantum reservoir with 10 qubits (comprising both input and hidden qubits such that $n_1+n_2=10$) and 3 layers, which means that features with lagged as $t-1$, $t-2$, and $t-3$ will be used to predict $RV_t$.
As mentioned before, we implement two quantum reservoir computing models, \texttt{QR1} and \texttt{QR2}. \texttt{QR1} does not utilize the ensemble reservoir approach, while \texttt{QR2} incorporates an additional reservoir to enhance computational capacity (see Fig.~\ref{Fig: virtual_node}). The quantum reservoir outputs, corresponding to evolved quantum states driven by lagged input encodings, are used as nonlinear regressors in a regularized linear model estimated via ridge regression. In both quantum reservoir models, the weight matrix $W_{out}$ is estimated using ridge regression, as given in Eq.~(\ref{eq:training_Wout}). Implementation-wise, we consider 100 different quantum reservoir instances and train and evaluate each of them separately. We then select the best-performing quantum reservoir and report the results. Similarly for \texttt{LSTM}, we test multiple hyperparameter configurations and report the best-performing results for presentation. Notice that we are reporting
the best-performing results for both \texttt{LSTM} and QRC, hence QRC does not specifically gain any unfair advantage.
\subsection{Closed-Loop Prediction of Multi-Step Future Outcomes }
So far, we have focused on the prediction of outcomes which are only one step ahead of the input. In other words, all models are trained under a one-step open-loop setting: the model inputs are ground-truth values, and performance is evaluated against the ground truth. For example,
\[
\widehat{RV_t} = f\!\left(RV_{t-3}, Xs_{t-3},\, RV_{t-2}, Xs_{t-2},\, RV_{t-1}, Xs_{t-1}\right),
\]
where $Xs$ denotes other (exogenous) variables~(see Table~\ref{Table : Macro_features}). If a model name does not include the \texttt{X} suffix, then $Xs$ is empty (i.e., no exogenous variables are used). This is called open-loop strategy in which historic data is used for predicting one step ahead. However, one might be interested to predict the outcomes well ahead of the current data. For instance, by accessing the historic data one might be interested to know the outcome of $S$ steps ahead. This can be done through a closed-loop procedure which is explained below. To accomplish this, we rely on a conventional trained open-loop model which predict $S=1$ step ahead. For predicting $S>1$ steps in the future, we repeatedly call the model and feed its previous predictions back to the model as part of the next input. Specifically:
\begin{itemize}
\item For $S=1$,
\[
\widehat{RV_t} = f\!\left(RV_{t-3}, Xs_{t-3},\, RV_{t-2}, Xs_{t-2},\, RV_{t-1}, Xs_{t-1}\right).
\]
\item For $S=2$,
\[
\widehat{RV_{t+1}} = f\!\left(RV_{t-2}, Xs_{t-2},\, RV_{t-1}, Xs_{t-1},\, \widehat{RV_t}, Xs_t\right).
\]
At this step, $\widehat{RV_t}$ is taken from the previous step's output. Since our models predict only $RV$ (and not the other variables), $Xs_t$ is still provided as ground truth. However, if such data is not available then a new predictor has to be trained targeting $Xs_t$ and use its predicted value (i.e. $\widehat{Xs}_t$) as the input of the closed-loop model. Here, for simplicity we assume that the ground truth of $Xs_t$ is available.
\end{itemize}
Noted that when $S=1$, the closed-loop setting degenerates to the open-loop setting.
Using the closed-loop strategy allows us to assess a model's ability to forecast into the future, which is crucial in many forecasting scenarios. We consider $S=1$ and $S=5$ to evaluate both short-term and longer-term predictive performance. It is worth emphasizing that by increasing $S$ the accuracy is expected to go down as the inputs come from prediction rather than the actual ground truth data. In other words, at every steps, a small error can propagate to the next steps affecting the accuracy of the model for predicting distant future outcomes.
\subsection{Feature Selection for Quantum Reservoir Computing}
For the classical models that we selected, the features listed in Table~\ref{Table : Macro_features} can be easily applied to the models.
However, considering the performance of currently available quantum computers, we have limited the system size of the quantum reservoir model to $10$ qubits. Among these qubits, we use $n_1$ qubits for encoding the features and reserve the rest (i.e., $n_2=10-n_1$) for hidden qubits. Consequently, we cannot apply all the features listed in Table~\ref{Table : Macro_features} to the quantum model at the same time.
Therefore, we need to further select the most significant features for our quantum reservoir computing models \texttt{QR1} and \texttt{QR2}.
To refine the key features that best fit our model, we employ one of the wrapper methods, which are model-related feature selection algorithms. The wrapper method treats the machine learning model as a black box, optimizing it to evaluate different subsets of features. This approach simplifies feature selection by framing it as a search problem to identify the best subset of features. The wrapper method requires a strategy to search through possible feature subsets. Common search strategies include exhaustive search, forward selection, backward elimination, and various heuristic approaches, such as genetic algorithms or simulated annealing.
In this study, we use forward selection as the search strategy, which is an iterative method, see Fig.~\ref{Fig: schematic_feature_selection}.
In forward selection, the process begins with three key components: (i) a prediction model; (ii) an initially empty feature set; and (iii) a feature pool containing all candidate features available for selection.
At each iteration, the method evaluates the performance of the prediction model by temporarily adding one feature at a time from the candidate pool to the current feature set.
The feature that results in the best performance in the current iteration is permanently added to the optimal feature set and removed from the feature pool.
The process continues until no feature from the candidate pool can significantly improve the model's performance or until a predefined stopping criterion is met.
The algorithm is schematically shown in Fig.~\ref{Fig: schematic_feature_selection}.
Through this incremental feature selection process, the dimension of the feature space can be effectively reduced while retaining the most predictive features. This not only improves the explainability and stability of the model, but also avoids overfitting, thus enhancing the accuracy of predicting volatility.
\begin{figure}
\centering
\includegraphics[width=0.35\textwidth]{Fig5.png}
\caption{\textbf{Feature selection by wrapper method.} The process begins with an initial pool of features, which contains all candidate features without prior filtering. At each iteration, one feature from the pool is considered as a candidate for inclusion in the optimal feature set (initially empty). Each candidate feature is temporarily added to the current optimal feature set, and the resulting feature set is then used to train a quantum reservoir model. The performance of the trained model is evaluated using the mean squared error (MSE) loss function. Among the candidates, the feature that leads to the lowest MSE is selected and permanently added to the optimal feature set, while the remaining candidates return to the pool. This selection process is repeated iteratively until a predefined stopping criterion is met, such as reaching a target performance threshold or exhausting all candidate features.}\label{Fig: schematic_feature_selection}
\end{figure}
We use the forward selection method with the stop condiation as the maximum featurs in the optimal features is $10$, see Fig.~\ref{Fig: schematic_feature_selection}, for both the \texttt{QR1} and the \texttt{QR2} models to identify the best subsets of features.
For the \texttt{QR1}, the optimal subset is $F^*_{\texttt{QR1}}=\{RV, MKT, DP, IP, RVq, STR, DEF\}$, and for the \texttt{QR2}, it is $F^*_{\texttt{QR2}}=\{RV, MKT, STR, RVq, EP, INF, DEF\}$. Thus, we have selected $n_1=7$ qubits to input features for each model, allowing $n_2=3$ hidden qubits to carry information from the past data to be used for future prediction.
We measure the forecast Mean Squared Error (MSE) of \texttt{QR1} and \texttt{QR2} with different subsets of features and present the results in Figs.~\ref{Fig: Performance_MSE}(a) and (b). These figures clearly show that as the number of optimal features increases, the performance of the model initially improves but then deteriorates. One reason for this is that as the number of features increases, more qubits are required to input the values. Since the total size of the quantum system is fixed, the number of qubits for hidden states decreases, preventing the quantum model from capturing the long-memory property of past information. Based on the results in Fig.~\ref{Fig: Performance_MSE}, the optimal number of features for the \texttt{QR1} and \texttt{QR2} models is $n_1=7$. In addition, they share several features, $\{RV, MKT, STR, RVq, DEF\}$.
\begin{figure}
\centering
\includegraphics[width=0.8\linewidth]{Fig6.png}
\caption{\textbf{Performance (MSE) of \texttt{QR1} and \texttt{QR2} models with their optimal feature sets.} (a) \texttt{QR1} model (orange), the best performance is achieved with the optimal features $\{RV, MKT, DP, IP, RVq, STR, DEF\}$. (b) \texttt{QR2} model (green), the best performance is achieved with the optimal features $\{RV, MKT, STR, RVq, EP, INF, DEF\}$. Order of features in the optimal feature set is left to right along horizontal axis. Note that the features included in these two sets are not the same and have different orderings.}
\label{Fig: Performance_MSE}
\end{figure}
To see how increasing the features affects the prediction quality, in Fig.~\ref{Fig: performance with diff F}(a), we plot the realized volatility as well as the \texttt{QR1} forecast for the first $1$, $4$, $7$, and $9$ features, which are shown in Fig.~\ref{Fig: Performance_MSE}.
As the figure clearly shows by increasing the features, the prediction improves until the number of features reach $n_1=7$. By further increasing the number of features, the prediction quality cannot be further improved because the number of hidden qubits decreases; thus, the reservoir lacks sufficient memory to incorporate past information for future prediction. The corresponding results for the \texttt{QR2} model is also shown in Fig.~\ref{Fig: performance with diff F}(b). The same conclusion about the impact of features on prediction quality can be deduced from the \texttt{QR2} model.
In addition, by comparing the corresponding panels in Fig.~\ref{Fig: performance with diff F}(a) Fig.~\ref{Fig: performance with diff F}(b) one can clearly see that the \texttt{QR2} outperforms the \texttt{QR1} prediction. This is expected as the the \texttt{QR2} uses twice more measurement data for predicting the future outcomes.
\subsection{Model Interpretability }
Any machine learning based algorithm should not only minimize its loss function, but preferably also explains which features contribute the most towards this goal. In the context of this work, we want to understand which external features contribute the most towards realized volatility. In the previous section devoted to the forward feature selection method, we adopted a constructive approach of adding one feature at a time from scratch. However, this leaves an equally important question untouched - given a set of features already being encoded into our simulation, how do we separate out individual contributions coming from each feature towards the predictive success of our quantum reservoir computing based method?
To answer this we employ the Shapley value method, a model-agnostic technique grounded in cooperative game theory.
Originally introduced in Ref.~\cite{shapley1953value}, the Shapley value offers a principled framework for fairly distributing the total payoff among participants in a coalition based on their individual contributions.
This approach considers all possible combinations and permutations of features, helping to quantify the contribution of each feature to the model prediction. \citet{lundberg2017unified} adapted this concept to machine learning, proposing Shapley values as a unified and interpretable measure of feature importance that satisfies key axioms such as consistency and local accuracy.
In this spirit, our analysis applies Shapley values to the quantum reservoir computing framework to assess how macro-financial variables contribute to realized volatility forecasts, bridging interpretability with quantum-enhanced prediction.
In theory, exact computation of Shapley values is exponential in the number of features; therefore, we typically rely on approximation methods in practice. In this work, we estimate Shapley values for the \texttt{QR1} and \texttt{QR2} models using the Monte Carlo sampling method~\cite{strumbeljExplainingPredictionModels2014}, as implemented in the Julia package \texttt{ShapML}.
The choice of the features is decided according to one of following three strategies:
\begin{enumerate}
\item \textbf{Individual feature contribution:} In this strategy, we evaluate the contribution of each feature at each lag separately, e.g. $\{RV_{t-3},MKT_{t-3},\cdots,RV_{t-2},MKT_{t-2},\cdots,\newline RV_{t-1},MKT_{t-1},\cdots\}$. The results are shown in Figs.~\ref{Fig: Shapley_values}(a) and (b), for \texttt{QR1} and \texttt{QR2} models, respectively. As expected, in both models, the $RV_{t-1}$ has the most contribution in predicting the future outcome $RV_t$. \\
\item \textbf{Feature-family contribution:} In this strategy, we evaluate the contribution of each feature family irrespective of the time lag. Therefore, during the evaluation of the Shapley value when we evaluate the contribution of one feature type, e.g. for realized volatility $RV$, we consider contributions from all $\lbrace RV_{t-3}, RV_{t-2}, RV_{t-1}\rbrace$ time-lagged steps. The results are shown in Figs.~\ref{Fig: Shapley_values}(c) and (d) for the \texttt{QR1} and \texttt{QR2} models respectively. Interestingly, the Shapley value-based method furnishes a somewhat different ordering of feature importance compared to the forward selection method results depicted in Fig.~\ref{Fig: Performance_MSE}. Note that this is not inconsistent and stems from the fact that the Shapley value is computed when all the features are embedded in the model while the forward selection method works by considering each individual feature in isolation. \\
\item \textbf{Time lagged feature contribution:} In this strategy, we evaluate the contribution of all the features at a given time-lag, e.g. $F_{t-k}=\{RV_{t-k},MKT_{t-k},\cdots\}$. We expect that the more recent data contributes the most towards future prediction compared to historical data. The results are shown in Figs.~\ref{Fig: Shapley_values}(e) and (f) for \texttt{QR1} and \texttt{QR2} respectively. As expected, in both models our expectation is borne out and the more recent data are indeed more useful.
\end{enumerate}
\begin{figure*}[t]
\centering
\includegraphics[width=1\textwidth]{Fig7.png}
\caption{\textbf{Forecasting performance of \texttt{QR1} (a) and \texttt{QR2} (b) models on realized volatility with different feature subsets.} Each subfigure presents the actual realized volatility (black line) of the S\&P 500 index alongside predictions from the quantum reservoir computing models (\texttt{QR1} in orange and \texttt{QR2} in green) using varying feature configurations. \texttt{QR2} demonstrates higher accuracy, especially during volatility peaks, suggesting it has better responsiveness to market fluctuations.}
\label{Fig: performance with diff F}
\end{figure*}
\begin{figure*}
\centering
\includegraphics[width=1\textwidth]{Fig8.png}
\caption{\textbf{Shapley-value based feature importance test.} The Shapley values for the optimal feature set are presented for \texttt{QR1} and \texttt{QR2}. Panels (a), (c), and (e) correspond to \texttt{QR1}, while panels (b), (d), and (f) correspond to \texttt{QR2}.
(a) and (b) show the Shapley values for each feature at different time lags.
(c) and (d) show the Shapley values for each feature, aggregated across all time lags.
(e) and (f) show the Shapley values for each lagged input.}
\label{Fig: Shapley_values}
\end{figure*}
\subsection{Performance Metrics and Forecast Evaluation}
The selection of appropriate performance metrics is pivotal in the evaluation of volatility forecasting models~\cite{HansenLunde2005}. In all the models, either classical or quantum, we have used MSE, defined in Eq.~\eqref{MSE_loss}, as a loss function which is minimized. Although MSE is a popular figure of merit for evaluating the performance of a time series predictor, we stress that it is blind to the direction of error. In other words, overestimation as well as underestimation, both get equally penalized. However, in the concrete case of realized volatility prediction, underestimation is far more dangerous as it can lead to serious undermining of emerging risks in the market. Thus, we specifically need a figure of merit which quantifies the degree of underestimation. In order to accomplish this task, we employ the widely recognized Quasi-Likelihood (QLIKE) as our figure of merit which captures both symmetric and asymmetric aspects of forecast errors, thereby providing a comprehensive assessment of model performance. The QLIKE is defined as:
\[
\text{QLIKE} = \frac{1}{T} \sum_{t=1}^{T} \left( \log \widehat{RV}_t^2 + \Big(\frac{RV_t}{\widehat{RV}_t}\Big)^2 \right),
\]
As it is evident from the definition, the QLIKE function penalizes under-predictions more heavily, reflecting the higher costs of underestimating volatility in risk management and option pricing. Furthermore, \citet{Patton2011} demonstrates that QLIKE is robust to measurement errors in volatility proxies.
For statistical evaluation of predictive accuracy, we employ the Model Confidence Set (MCS) procedure proposed by \citet{HansenLundeNason2011} and the Diebold-Mariano (DM) test \citep{DieboldMariano1995}. We utilize MSE as our primary loss function due to its simplicity, while also reporting QLIKE values to assess asymmetry and penalize underestimation of volatility more strongly.\\
\emph{Model Confidence Set.--}
The MCS procedure is designed to identify a subset of forecasting models whose predictive performance is statistically indistinguishable from that of the best model at a specified confidence level. The procedure begins with an initial set of candidate models, denoted $\mathcal{M}_0$, and iteratively eliminates inferior models based on tests of equal predictive ability. Formally, let $L_{i,t}$ denote the loss (e.g., squared error) incurred by model $i$ at time $t$ and define the loss differential between models $i$ and $j$ as:
\[
d_{ij,t} = L_{i,t} - L_{j,t}, \quad i, j \in \mathcal{M}_0.
1\]
The null hypothesis of equal predictive ability (EPA) is:
\[
H_{0,\mathcal{M}} : \mu_{i,j} = \mathbb{E}[d_{ij,t}] = 0 \quad \forall i, j \in \mathcal{M},
\]
where $\mathcal{M} \subseteq \mathcal{M}_0$ denotes the set of models under consideration at a given iteration. Rejection of $H_{0,\mathcal{M}}$ implies that at least one model in $\mathcal{M}$ exhibits significantly inferior forecasting performance.
To test $H_{0,\mathcal{M}}$, a range-type test statistic is employed:
\[
T_{R,\mathcal{M}} = \max_{i \in \mathcal{M}} t_i,
\]
where $t_i$ denotes a studentized statistic comparing the performance of model $i$ to others (e.g., via average loss differences). If the null is rejected, the model with the most evidence against it is eliminated according to the elimination rule:
\[
e_{\mathcal{M}} = \arg\max_{i \in \mathcal{M}} \left( \sup_{j \in \mathcal{M}} \hat{d}_{i,j} \right),
\]
where $\hat{d}_{i,j}$ is the sample mean of the loss differential. This elimination process is repeated until the null hypothesis is no longer rejected.
The final surviving set of models is denoted by $\widehat{\mathcal{M}}^*_{1-\alpha}$, the $(1-\alpha)$ MCS, which is asymptotically guaranteed to contain the model(s) with the lowest expected loss with probability at least $1 - \alpha$. In this study, we set $\alpha = 0.05$, ensuring that the resulting confidence set contains the best-performing model(s) with at least 95\% probability in large samples. We complement the MCS analysis with pairwise Diebold-Mariano tests to assess the statistical significance of forecast performance differentials between models. This dual approach allows us to validate robustness in model rankings and mitigate the limitations associated with relying on a single evaluation metric.\\
\emph{Diebold-Mariano test (DM test).-- } The DM test is a particular adaptation of the above procedure for time-series data, and evaluates the null hypothesis of equal predictive accuracy between two competing models, formally stated as:
\[
H_0: \mathbb{E}[d_t] = 0,
\]
where \( d_t = L_{1,t} - L_{2,t} \) is the loss differential at time \( t \) between models 1 and 2. Here, \( L_{i,t} \) represents the loss at time \( t \) for model \( i \), and the choice of loss function \( L_{i,t} \) (e.g., MSE in this paper) reflects the predictive objective. The DM test statistic is computed as:
\[
\text{DM statistic} = \frac{\bar{d}}{\sqrt{\widehat{\sigma}^2 / T}},
\]
where \( \bar{d} = \frac{1}{T} \sum_{t=1}^T d_t \) is the sample mean of the loss differential, \( \widehat{\sigma}^2 \) is the Newey-West adjusted variance of \( d_t \), and \( T \) is the sample size. The test accounts for potential autocorrelation and heteroskedasticity in \( d_t \), making it suitable for time-series data with serially correlated forecast errors.
\begin{table*}[hbt]
\centering
\caption{\textbf{Model Confidence Set.} Performance of classical and quantum Models in predicting realized volatility with predicting $S$-ahead step. The MSE and QLIKE are presented alongside the MCS $p$-values ($P_{\text{MCS}}$). Asterisks ($^*$) indicate inclusion in the model confidence set. The quantum reservoir computing models outperform all classical models, with \texttt{QR2} exhibiting the best performance across all measures.}
\begin{tabularx}{\textwidth}{lcccc}
\multicolumn{5}{c}{\textbf{S=1 (entire sample: 1997.08-2017.12) (1950.01- 2017.12)}} \\
\hline
Model & \multicolumn{2}{l}{MSE} & \multicolumn{2}{l}{QLike} \\
& Loss of MSE & $P_{MCS}$ of MSE & Loss of QLike & $P_{MCS}$ of QLike \\
\hline
\multicolumn{5}{c}{\textbf{Classical Time Series Models}} \\
\hline
\texttt{HAR} & 0.1476 & 0.0004 & 2.0431 & 0.0008 \\
\texttt{HARX} & 0.1508 & 0.0004 & 2.2436 & 0.0008 \\
\texttt{AR1} & 0.1304 & 0.0065 & 1.7279 & 0.0050 \\
\texttt{AR3} & 0.1178 & 0.0936 & 1.5893 & 0.0861 \\
\texttt{ARMAX} & 0.1145 & $0.4406^*$ & 1.6196 & $0.4355^*$ \\
\texttt{LSTM} & 0.1295 & 0.0221 & 1.7909 & 0.0188 \\
\texttt{LSTMX} & 0.1185 & $0.4406^*$ & 1.7571 & $0.4355^*$ \\
\texttt{RC} & 0.1441 & 0.0084 & 2.1011 & 0.0061 \\
\texttt{RCX} & 0.1089 & $0.6086^*$ & 1.6480 & $0.6106^*$ \\
\hline
\multicolumn{5}{c}{\textbf{Quantum Reservoir Computing}} \\
\hline
\texttt{QR1} & 0.105 & $0.7603^*$ & 1.4427 & $0.7510^*$ \\
\texttt{QR2} & \textbf{0.103} & $1.0000^*$ & \textbf{1.4004} & $1.0000^*$ \\
\hline
\end{tabularx}
\vspace{1ex}
\begin{minipage}{\textwidth}
\footnotesize
\end{minipage}
\begin{tabularx}{\textwidth}{lcccc}
\multicolumn{5}{c}{\textbf{S=5 (entire sample: 1998.01-2017.12) (1950.01- 2017.12)}} \\
\hline
Model & \multicolumn{2}{l}{MSE} & \multicolumn{2}{l}{QLike} \\
& Loss of MSE & $P_{MCS}$ of MSE & Loss of QLike & $P_{MCS}$ of QLike \\
\hline
\multicolumn{5}{c}{\textbf{Classical Time Series Models}} \\
\hline
\texttt{HAR} & 0.2143 & $0.1291^*$ & 2.9041 & $0.1275^*$ \\
\texttt{HARX} & 0.2934 & 0.0044 & 4.5800 & 0.0040\\
\texttt{AR1} & 0.2642 & $0.0938^*$ & 3.4136 & $0.0937^*$ \\
\texttt{AR3} & 0.2134 & $0.1302^*$ & 2.8369 & $0.1297^*$ \\
\texttt{ARMAX}& 0.2134 & $0.1302^*$ & 3.0703 & $0.1297^*$ \\
\texttt{LSTM} & 0.1831 & $0.1291^*$ & 2.4600 & $0.1275^*$ \\
\texttt{LSTMX}& 0.2200 & $0.0925^*$ & 3.4512 & $0.0851^*$\\
\texttt{RC} & \textbf{0.1528} & $1.0000^*$ & \textbf{2.0551} & $1.0000^*$ \\
\texttt{RCX} & 0.1667 & $0.6333^*$ & 2.4605 & $0.6251^*$ \\
\hline
\multicolumn{5}{c}{\textbf{Quantum Reservoir Computing}} \\
\hline
\texttt{QR1} & 0.1556 & $0.7642^*$ & 2.1518 & $0.7703^*$\\
\texttt{QR2} & 0.1663 & $0.6333^*$ & 2.2332 & $0.6251^*$ \\
\hline
\end{tabularx}
\label{Table: performance_benchmarking}
\end{table*}
\begin{table*}
\centering
\caption{\textbf{Diebold-Mariano Test.} DM statistic values (lower triangular matrix) and $p$-values (upper triangular matrix) for model comparisons. The statistic values closer to zero indicate that the pair of models have similar predictive accuracy. Lower $p$-values suggest rejection of the null hypothesis of equal predictive ability.}
\begin{tabularx}{\textwidth}{lccccccccccc}
\multicolumn{12}{c}{\textbf{}} \\
\hline
& \texttt{HAR} & \texttt{HARX} & \texttt{AR1} & \texttt{AR3} & \texttt{ARMAX} & \texttt{LSTM} & \texttt{LSTMX} & \texttt{RC} & \texttt{RCX} & \texttt{QR1} & \texttt{QR2} \\
\hline
\texttt{HAR} & & 0.575 & 0.203 & 0.004 & 0.009 & 0.006 & 0.01 & 0.801 & 0.001 & 0 & 0 \\
\texttt{HARX} & -0.562 & & 0.116 & 0.001 & 0.003 & 0.004 & 0.002 & 0.597 & 0 & 0 & 0 \\
\texttt{AR1} & 1.275 & 1.577 & & 0.002 & 0.077 & 0.925 & 0.301 & 0.154 & 0.009 & 0.001 & 0 \\
\texttt{AR3} & 2.929 & 3.313 & 3.12 & & 0.652 & 0.065 & 0.937 & 0.003 & 0.17 & 0.036 & 0.014 \\
$\texttt{ARMAX}$ & 2.636 & 2.999 & 1.775 & 0.451 & & 0.134 & 0.657 & 0.006 & 0.393 & 0.12 & 0.111 \\
\texttt{LSTM} & 2.764 & 2.938 & 0.094 & -1.856 & -1.502 & & 0.268 & 0.168 & 0.017 & 0.006 & 0.002 \\
\texttt{LSTMX} & 2.583 & 3.168 & 1.037 & -0.079 & -0.445 & 1.109 & & 0.028 & 0.229 & 0.114 & 0.09 \\
\texttt{RC} & 0.252 & 0.53 & -1.431 & -3.032 & -2.758 & -1.383 & -2.206 & & 0 & 0 & 0 \\
\texttt{RCX} & 3.395 & 4.001 & 2.629 & 1.376 & 0.856 & 2.403 & 1.205 & 3.639 & & 0.501 & 0.296 \\
\texttt{QR1} & 3.695 & 4.149 & 3.322 & 2.104 & 1.561 & 2.788 & 1.584 & 3.82 & 0.674 & & 0.771 \\
\texttt{QR2} & 3.806 & 4.289 & 3.577 & 2.465 & 1.6 & 3.111 & 1.701 & 4.584 & 1.048 & 0.291 & \\
\hline
\end{tabularx}
\label{Table: DM_test}
\end{table*}
\emph{MCS results.--} The results of MCS are presented in Table~\ref{Table: performance_benchmarking}. The table reports the performance metrics (MSE and QLIKE) along with the MCS $p$-values for classical time series models and quantum reservoir computing models.
In particular, when S=1, namely one-step ahead prediction, quantum models, especially \texttt{QR2}, demonstrate superior performance across all measures. From the upper table in Table~\ref{Table: performance_benchmarking}, we further observe that the \texttt{QR2} model achieves the lowest MSE and QLIKE values, indicating its superior predictive accuracy. The MCS $p$-values further reinforce this finding, as \texttt{QR2} attains a $p$-value of 1.0, signifying its inclusion in the superior set of models at the 95\% significance level. In contrast, most classical models exhibit significantly higher loss values and lower MCS $p$-values, suggesting inferior performance. Examining the MCS $p$-values, models with asterisks are included in the Model Confidence Set $\mathcal{M}_{\text{MCS}}$, indicating that their predictive accuracy is statistically indistinguishable from the best-performing model. Among the classical models, only \texttt{ARMAX}, \texttt{LSTMX}, and \texttt{RCX} are included in the MCS based on their $p$-values, although their loss metrics are still higher than those of the quantum models. Additionally, \texttt{QR1} is also included in the MCS, further highlighting the strong performance of quantum models. Examining the results, we also observe that including additional features sometimes enhances the model performance. For example, \texttt{ARMAX} and \texttt{LSTMX} exhibit lower loss values and higher MCS $p$-values compared to their counterparts \texttt{AR3} and \texttt{LSTM}, suggesting that the inclusion of characteristics improves the accuracy of forecasting in these models. However, this improvement is not consistent across all models. Notably, \texttt{HARX} performs slightly worse than \texttt{HAR} in both MSE and QLIKE metrics, indicating that additional features do not always contribute to better performance. This suggests that the effectiveness of including features depends on the model structure and the relevance of the features.
When consider $S=5$, namely, the long-term prediction, both Quantum and Classical reservoir demonstrated outstanding performance compared with other models.
The performances of \texttt{QR1} and \texttt{QR2} are very similar to that of the classic reservoir computing. Considering that we only used a simple reservoir model, these indirectly demonstrate the capabilities of Quantum reservoir.
\\
\emph{Diebold-Mariano Test results.-- } Table~\ref{Table: DM_test} presents the Diebold-Mariano test statistic(s) and corresponding $p$-value(s) for pairwise model comparisons, providing further evidence of the quantum models' outperformance. Analyzing the DM test results, we find that the quantum models (\texttt{QR1} and \texttt{QR2}) significantly outperform most classical models. For instance, the $p$-values for comparisons between \texttt{QR2} and classical models like \texttt{HAR} and \texttt{AR1} are effectively zero, indicating strong evidence against the null hypothesis of equal predictive accuracy. Additionally, the test statistics are positive and relatively large, further supporting the superiority of the quantum models. Interestingly, the comparison between \texttt{QR1} and \texttt{QR2} yields a high $p$-value (0.771), suggesting no significant difference in predictive accuracy between the two quantum models. This implies that both quantum models perform comparably well, although \texttt{QR2} has a slight edge based on the loss metrics in Table~\ref{Table: performance_benchmarking}. The results also reveal that among classical models, \texttt{AR3}, and \texttt{ARMAX} exhibit relatively better performance, with higher $p$-values when compared to other classical models. However, they still fall short when compared to the quantum models. The results further illustrate the mixed impact of including additional features. For some model pairs, such as \texttt{LSTM} and \texttt{LSTMX}, the inclusion of features leads to significant improvements, as evidenced by lower test statistics and higher $p$-values. In contrast, comparisons between \texttt{AR3} and \texttt{ARMAX} yield $p$-values that indicate no significant difference, implying that the additional features in \texttt{ARMAX} do not provide substantial benefits over \texttt{AR3}. Overall, the results suggest that while incorporating additional features can enhance model performance in certain cases, it does not universally guarantee better forecasts. The effectiveness of additional features appears to be model-specific, highlighting the importance of selecting relevant variables that contribute meaningfully to volatility forecasting. \\
Overall, the combined evidence from the loss functions, MCS procedure, and DM test underscores the superior predictive performance of the quantum reservoir computing models over classical time series and machine learning models in forecasting realized volatility.
\section{Conclusions and Outlook}
\label{sec: conclusions and outlook}
Quantum reservoir computing presents a promising new paradigm for the modeling and forecasting of time series data. In this study, we have demonstrated how a specific instantiation of quantum reservoir computing that utilizes disordered spin systems as reservoirs can be effectively employed for forecasting realized volatility in financial markets. Our empirical evaluation reveals that even with a modest quantum architecture comprising only ten qubits and a moderated number of macroeconomic input features per instance, the proposed framework outperforms classical models in predictive accuracy. While we do not advance any claim to quantum supremacy in the rigorous sense of the term since this is a phenomenological study, the results obtained herein do indicate a potential future application of quantum reservoir computing towards applying them to other similar financial problems. Our protocol is especially suitable to use in real-life quantum computers based on trapped ion platforms, in particular, \cite{kielpinski2002architecture,erhard2021entangling,akhtar2023high}. These setups offer excellent all-to-all connectivity between spins, as required in our specific choice of the quantum reservoir. We anticipate that alternative quantum reservoir architectures with different connectivity constraints may yield comparable predictive performance, although empirical verification would be necessary to substantiate this speculation. \\
Quantum computing for machine learning still faces significant challenges, both in hardware and algorithmic development. On the hardware side, current quantum devices are still in their infancy, constrained by limited qubit counts, restricted connectivity, finite coherence times, and imperfect readout. Overcoming these limitations will require advances in error correction—a milestone likely still a decade away. Nevertheless, near-term quantum devices provide valuable opportunities to develop and test small-scale algorithms, paving the way for large-scale applications on future error-corrected quantum computers. From an algorithmic perspective, it remains unproven whether quantum computers can outperform classical machine learning methods on classical data. In this work, we demonstrate that quantum computers can indeed solve prediction tasks in stochastic time series. While our quantum reservoir learning approach shows advantages over several classical methods, we do not claim a definitive quantum advantage, as such a conclusion would require rigorous theoretical proof beyond the scope of this paper. Moreover, intrinsic limitations on QRC like the exponential concentration problem \cite{Xiong2025on}, whereby the predictions may become input-independent due to concentration of expectation values for many-qubit reservoirs, exist and have to be addressed in future work.\\
\noindent \textbf{Code Availability.} The Code used in this paper can be found \href{https://github.com/LeeQY1996/Quantum-Reservoir-computing-for-Realized-Volatility-Forecasting}{here} in the Github Repositories. \\
\noindent \textbf{Acknowledgments.} AB acknowledges support from the National Natural Science Foundation of China (grants No. 12274059,
No. 12574528, No. 1251101297 and No. W2541020).
\bibliography{rv_forecasting_PRR}