Keywords: volatility spillover; realized volatility; graph signal processing; spectral analysis.

1 Introduction↩︎

In the complex ecosystem of global financial markets, understanding and forecasting volatility is critical for investors, risk managers and policymakers. To capture the dynamics of financial volatility, numerous measures have been proposed [1]. Among these, the high-frequency-data-based Realized Volatility (RV) [2] has become a standard tool in both research and practice. It is widely used for forecasting, risk management and asset pricing.

Volatility is rarely confined to individual assets or markets-it propagates across financial systems via spillovers, contagion and co-movements, especially during crises. Recent studies on financial volatility analysis have highlighted volatility spillover phenomena, with interest intensifying in the wake of the 2008 Global Financial Crisis (GFC) and the COVID-19 pandemic. This phenomenon refers to the transmission of volatility shocks from one market (or asset) to others [1], [3]–[7]. The significance of the volatility spillover effect necessitates a comprehensive volatility forecasting framework that jointly captures the cross-sectional volatility interconnections across different financial markets and the temporal autocorrelation of each market’s volatility. Recent work has applied graph-based models to RV forecasting by exploiting spillover networks [8], demonstrating that explicitly specifying the structure of inter-market linkages can enhance the RV forecast accuracy.

Yet, an important aspect remains largely unexplored in the graph-based RV forecasting studies: the lack of a systematic way to evaluate how well a given network construction represents the true interdependence structure. Since different graph-building methods (e.g., correlation-based, variance decomposition–based, or directional measures) can yield vastly different network topologies, the absence of an evaluation framework, which is agnostic to volatility forecasting models, makes it difficult to judge which constructions meaningfully represent volatility dynamics. Without such diagnostics, forecasts may rely on networks whose structure is potentially misaligned with the underlying economic reality, reducing interpretability and performance. To bridge this gap, we introduce graph signal energy (GSE) as a novel diagnostic tool to assess the informativeness or effectiveness of different network constructions before they are embedded into forecasting models. The concept builds on graph signal processing (GSP), a rapidly growing field that generalizes classical signal processing to graph-structured data [9]. In GSP, the value of a variable observed on each node is treated as a graph signal, while graph edges encode relationships or dependencies among nodes. The GSE then quantifies how “smooth” or “irregular” this signal is with respect to the graph topology, based on the graph Laplacian quadratic form.

Applied to RV spillover networks in a rolling-window setting, GSE provides an intuitive and theoretically grounded way to assess whether a network topology is informative for capturing volatility interdependencies. Especially, we expect RV GSE to behave differently in turbulent versus stable periods because the strength of cross-market volatility linkage varies across different market conditions [10], [11]. For illustration, we compare RV GSE patterns for networks constructed from the Pearson correlation matrix and the Diebold-Yilmaz (DY) framework [10] combined with the magnetic Laplacian matrix [12]. In Figure 1, the RV GSE of the Pearson correlation–based network remains unexpectedly low during major volatility spikes such as the European Debt Crisis and the COVID-19 pandemic, but reaches a relatively high level in the comparatively stable year of 2005. In contrast, Figure 2 shows that the RV GSE of the DY-based network with the magnetic Laplacian exhibits the expected dynamics. It rises during turbulent periods, including the GFC, European Debt Crisis, and COVID-19, and remains subdued in tranquil periods. These findings demonstrate that GSE provides a powerful, model-agnostic lens for evaluating graph construction methods, and suggest that the DY-magnetic-Laplacian construction method offers a more informative and theoretically consistent representation of volatility spillover networks than the Pearson correlation approach. Further details of the GSE computation and its integration into volatility spillover analysis with theoretical justifications are provided in Sections 2.1 and 3.1.

Having identified the DY-magnetic-Laplacian method yields a more informative representation of the RV spillover network, this study next employs GSP techniques to fully exploit the topological advantages of the chosen graph. GSP provides a principled framework for analyzing and forecasting signals defined over networks, making it particularly well-suited to capture the cross-market volatility dynamics embedded in the spillover structure. The integration of the magnetic Laplacian with GSP methods is not merely a technical choice, but a necessity: it allows directional spillovers to be incorporated while ensuring that the spectral properties of the network are leveraged in the modeling process, which is ideal for the chosen DY-magnetic-Laplacian construction method. The methodological rationale is discussed in detail in Section 3.

Building on this foundation, we propose a new GSP-based forecasting method that extends the well-established Heterogeneous Autoregressive (HAR) framework for RV forecasting. The HAR model, introduced by [13] and subsequently extended by [14], [7] and others, has remained popular in financial econometrics due to its balance of simplicity and empirical effectiveness. By embedding HAR dynamics into a GSP framework (GSP-HAR), this study aims to enhance the multi-asset RV forecasting accuracy by explicitly accounting for volatility interdependencies across markets, which is not captured by the traditional HAR models.

Therefore, this paper makes two main contributions to the literature on RV forecasting. Firstly, this study employs the GSE as a diagnostic tool for evaluating the effectiveness of different volatility spillover network constructions. By examining how well the time-varying behavior of RV GSE aligns with financial market conditions, the analysis provides an intuitive and model-agnostic way to verify whether a network topology meaningfully reflects volatility interdependencies. This approach highlights a straightforward yet often overlooked principle: more accurate volatility forecasts can be achieved by first identifying the most informative network construction and then tailoring graph-based models to that structure.

Secondly, the research introduces the GSP-HAR framework, integrating GSP into the well-established HAR model. While existing HAR extensions incorporate volatility spillover effects by adding intermarket variables and thereby treating volatility linkages as exogenous covariates [8], [15], [16], the GSP-HAR framework instead embeds the volatility spillover effects within the graph spectral domain. By combining the volatility interrelationships and historical RV patterns, the proposed approach explicitly models both temporal and cross-market volatility dynamics. Empirical results demonstrate that this integration not only captures volatility spillovers more effectively, but also delivers substantial improvements in RV forecasting accuracy.

The paper is organized as follows. Section 2 provides preliminaries of pertinent concepts and models. The relevant volatility GSE calculation, the application of the magnetic Laplacian and GSP techniques and the proposed framework are developed in Section 3. Section 4 presents the conducted experiments and their results. The implementation code is available at https://github.com/MikeZChi/GSPHAR.git. Section 5 concludes the paper.

2 sec:Background↩︎

Before presenting the construction of the volatility spillover network and the proposed GSP-HAR model, we first review the fundamental components of graph theory along with the commonly used measures and representations of graph topological structure. These include the weight matrix, the degree matrix and the Laplacian matrix. Particularly, the Laplacian matrix serves as a key tool for calculating the GSE. Additionally, we briefly review the HAR model and its volatility spillover extensions, which are commonly employed in RV forecasting.

2.1 sec:Graph32Theory32Fundamentals↩︎

A graph can be denoted as \(\mathcal{G} = (\mathcal{V}, \mathcal{E})\), where \(\mathcal{V}\) is the node set and \(\mathcal{E}\) is the edge set of the graph. An edge \(e_{ij} \in \mathcal{E}\) can be formulated as \(e_{ij} = (v_i, v_j)\), which indicates that the edge \(e_{ij}\) starts from node \(v_i\), and ends at node \(v_j\). In general, we consider the graph as a directed graph, thus \(e_{ij} \neq e_{ji}\). However, if \(e_{ij} = e_{ji}\) for all \(i,j\), the graph is undirected. Examples of undirected and directed graphs can be found in Appendix 6. In a directed graph, two types of neighbors can be formulated for a node \(v_i\). The outgoing neighbors are all nodes \(v_j\) that \(v_i\) points to: \(N_i^{out} = \{v_j \;| \;e_{ij} \in \mathcal{E}\}\). Besides, the incoming neighbors are all nodes \(v_j\) that point to \(v_i\) points: \(N_i^{in} = \{v_j \;| \;e_{ji} \in \mathcal{E}\}\).

Edges can have weights to represent the importance or strength of the connection. The weight matrix \(\mathbf{W}\) is used to capture the structure of graphs. Suppose the weight attached to edge \(e_{ij} \in \mathcal{E}\) is a positive real number \(w_{ij}\). If the graph \(\mathcal{G}\) has \(N\) nodes, the formulation of the weight matrix \(\mathbf{W} \in \mathbb{R}^{N \times N}\) is presented below: \[\label{directed32weighted32adj95mx} \mathbf{W}_{ij} = \begin{cases} w_{ij}, & e_{ij} \in \mathcal{E};\\ 0, &e_{ij} \notin \mathcal{E}. \end{cases}\tag{1}\] For an undirected graph, its weight matrix is always symmetric because if \(e_{ij} \in \mathcal{E}\), \(e_{ji} = e_{ij}\) and, thus, \(w_{ji} = w_{ij}\). On the other hand, the weight matrix of a directed graph is asymmetric. The degree matrix \(\mathbf{D}\) is formulated based on the weight matrix as \(\mathbf{D} = \text{diag}(d_1, \ldots, d_N)\), where \(d_i = \sum_{j = 1}^{N} \mathbf{W}_{ij}\) measures the total connection strength of node \(v_i\) to all others.

Besides, the Laplacian matrix of the graph is formulated as \(\mathbf{L} = \mathbf{D} - \mathbf{W}\). Furthermore, the Laplacian matrix can be normalized as \(\mathbf{L}_{norm} = \mathbf{I} - \mathbf{D}^{-\frac{1}{2}}\mathbf{W}\mathbf{D}^{-\frac{1}{2}}\) to reduce the influence of node degree variations, and become more numerically stable. To ensure that the degree matrix \(\mathbf{D}\) is invertible, all its diagonal elements \(d_i (i = 1, \ldots, N)\) must be positive. This can be easily achieved by assuming that the graph contains no isolated nodes. This assumption is reasonable in the real financial world because volatility is seldom limited to individual assets or markets. Since ensuring the invertibility of the degree matrix is straightforward, the following discussions adopt the normalized Laplacian matrix, and the subscript \(_{norm}\) is omitted for notational simplicity.

First, we consider an undirected graph \(\mathcal{G}\), consisting of \(N\) nodes with no isolated nodes, a symmetric weight matrix \(\mathbf{W}\) with non-negative entries, and an invertible degree matrix \(\mathbf{D}\). A graph signal \(x\) is defined as a mapping from \(\mathcal{V}\) to \(\mathbb{R}\) for each node in the graph. Given the finite cardinality \(N\) of \(\mathcal{V}\), the graph signal can be regarded as the collection of signal values observed at the nodes of the graph, and it is represented as a vector \(\mathbf{x} \in \mathbb{R}^{N}\), where each component corresponds to the node value under specific fixed node indexing. The GSE \(E(\mathbf{x}) \in \mathbb{R}\) quantifies how much the signal changes across connected nodes. It is calculated as: \[\label{graph95signal95energy} E(\mathbf{x}) = \mathbf{x}^\top \mathbf{L} \mathbf{x}.\tag{2}\]

Thus, the GSE is capable of measuring the lack of smoothness of the graph signal \(\mathbf{x}\) over the graph structure. The proof of Theorem 1 is provided in Appendix 7.1.

2.2 sec:HAR32Model32and32Its32Extensions32on32Volatility32Spillovers↩︎

The HAR framework is designed to account for the fact that different market participants, such as quantitative traders, institutional investors and long-term investors, operate on different time horizons. The HAR model incorporates pooled RV data from these various time scales, typically daily, weekly and monthly, to represent the volatility behaviors. This design helps to capture the heterogeneous nature of volatility patterns in financial markets.

Suppose \(\mathbf{v}_t\) contains RV observations of \(N\) stock market indices on common trading day \(t\): \(\mathbf{v}_t = (\text{RV}_{1,t}, \text{RV}_{2,t}, \ldots, \text{RV}_{N,t})^\top\). For integers \(n > 1\), \(\mathbf{v}_{t-n:t-1}\) is the average-pooled data: \(\mathbf{v}_{t-n:t-1} = \frac{1}{n-1}\sum_{j = t-n} ^ {t-1} \mathbf{v}_j\). For notational simplicity, set the past daily, weekly and monthly RV pattern as \(\mathbf{v}_d = \mathbf{v}_{t-1}\), \(\mathbf{v}_w = \mathbf{v}_{t-5:t-1}\) and \(\mathbf{v}_m = \mathbf{v}_{t-22:t-1}\). \(\widehat{\mathbf{v}}_t\) denotes the forecast RV value on day \(t\). The conventional HAR model is formulated as: \[\label{HAR} \widehat{\mathbf{v}}_t = \boldsymbol{\alpha} + \boldsymbol{\beta}_{d} \mathbf{v}_d + \boldsymbol{\beta}_{w} \mathbf{v}_w + \boldsymbol{\beta}_{m} \mathbf{v}_m,\tag{3}\] where coefficients \(\boldsymbol{\beta}_{d}\), \(\boldsymbol{\beta}_{w}\) and \(\boldsymbol{\beta}_{m}\) are diagonal matrices of size \(N\) and they measure how the past daily, weekly and monthly volatility patterns influence future RVs. The intercept term \(\boldsymbol{\alpha}\) can be a vector. The HAR model cannot capture the comovement or correlations of the volatility in different stock markets, which limits its ability to produce more accurate forecasts.

The GNN-HAR model adds a \(L\)-layer spatial Graph Neural Network (GNN) component to the HAR model as an additional input variable to capture the volatility spillover effect in a nonlinear way. The GNN component relies on a precomputed volatility relational network \(\mathcal{G}\) where the nodes are markets (or assets) and the edges are the volatility linkages. The weight matrix \(\mathbf{W}\) is provided and remains fixed during the GNN-HAR training process. The GNN component aggregates past daily, weekly, and monthly volatility patterns by degree-weighted averaging over the neighborhood of each node (market), which is defined by the volatility spillover network. This aggregation is followed by learnable linear and non-parametric non-linear transformations. The model is formulated as: \[\label{GNN-HAR} \begin{align} \mathbf{H}^{(0)} =&\; (\mathbf{v}_d, \mathbf{v}_w, \mathbf{v}_m),\\ \mathbf{H}^{(1)} =&\; \text{ReLU}(\mathbf{D}^{-\frac{1}{2}} \mathbf{W} \mathbf{D}^{-\frac{1}{2}}\mathbf{H}^{(0)}\mathbf{W}^{(0)}),\\ & \ldots,\\ \mathbf{H}^{(L+1)} =&\; \text{ReLU}(\mathbf{D}^{-\frac{1}{2}} \mathbf{W} \mathbf{D}^{-\frac{1}{2}}\mathbf{H}^{(l)}\mathbf{W}^{(L)}),\\ \widehat{\mathbf{v}}_t =&\; \boldsymbol{\alpha} + \beta_d \mathbf{v}_d + \beta_w \mathbf{v}_w + \beta_m \mathbf{v}_m + \boldsymbol{\gamma} \mathbf{H}^{(L+1)} +\boldsymbol{\epsilon}_t. \end{align}\tag{4}\] Here, \(\{\mathbf{W}^{(l)}\}_{l=0}^L\) are trainable GNN parameters for the linear transformation, and ReLU\((\cdot)\) function transforms the aggregated volatility information nonlinearly. \(\boldsymbol{\gamma}\) also contains trainable parameters to measure the overall influence of related stock market indices on the forecast. The GNN-HAR model uses spatial GNNs to capture the nonlinear relationships between stock market indices and generates more accurate RV forecasts compared to the HAR model [8].

There are several other HAR model extensions, such as the VHAR model [15] and the HAR-KS model [16], which expand the linear combination of the HAR framework to incorporate historical volatility patterns of all entities in the system when generating RV forecasts for a target entity. However, these approaches are less powerful because they provide limited information about the structure of the volatility spillover network, and oversimplify the volatility spillover effect by assuming it is linear.

3 sec:Methodology32and32Implementation↩︎

In this section, the GSE is defined and analyzed for volatility spillover networks constructed using different methods. Building on the insights from the GSE analysis, we then discuss the motivation, necessity and detailed implementation of the proposed GSP-HAR framework.

Denote the RV time series as \(\{\mathbf{v}_t \in \mathbb{R}^N \}_{t = 1}^T\), which consists of the RV time series of \(N\) different stock market indices on \(T\) common trading days. Specifically, each \(\mathbf{v}_t = (v_{1,t}, v_{2,t}, \dots, v_{N,t})^\top \in \mathbb{R}^N\) represents the RV observations of \(N\) markets at trading day \(t\).

3.1 sec:Volatility32Spillover32Network32and32Its32GSE↩︎

Unlike traffic networks [17] where linkages correspond to tangible infrastructure, such as roads between cities, financial networks of stock market indices lack observable physical connections. Although the geographical distance between markets can be used to define volatility linkages, inferring the volatility spillover network of global stock market indices from data offers a more informative representation of inter-market volatility relationships. Here, three common data-driven methods for constructing volatility networks are discussed. These methods take the multi-asset RV time series \(\{\mathbf{v}_t \in \mathbb{R}^N \}_{t = 1}^T\) as input, and output the weight matrix \(\mathbf{W}\) of the volatility spillover network. The use of RV time series data, as opposed to return data, to construct the volatility spillover network is motivated by the argument that the volatility shows a higher degree of serial dependence than returns, particularly when measured at relatively high frequencies [18]. A brief summary of the four weight matrices produced by the three different volatility network construction methods is provided below, whereas their detailed calculation processes are described in Appendix 8.

• The Pearson correlation matrix method: the symmetric weight matrix is produced based on the Pearson correlation matrix according to Appendix 8.1: \(\mathbf{W} = \mathbf{W}^{P}\).

• The GLASSO precision matrix method: the symmetric weight matrix is produced based on the GLASSO precision matrix according to Appendix 8.2: \(\mathbf{W} = \mathbf{W}^{GL}\).

• The DY-symmetrization method: the symmetric weight matrix is produced based on the DY framework according to Appendix 8.3: \(\mathbf{W} = \frac{1}{2} (\mathbf{W}^{DY} + (\mathbf{W}^{DY})^\top)\).

• The DY-magnetic-Laplacian method: the asymmetric weight matrix is produced based on the DY framework according to Appendix 8.3: \(\mathbf{W} = \mathbf{W}^{DY}\). The magnetic Laplacian refers to the corresponding Laplacian matrix construction method for the asymmetric weight matrix.

As discussed in Section 2.1, in order to calculate the non-negative GSE based on the normalized Laplacian matrix \(\mathbf{L}\) (Equation 2 ), the weight matrix \(\mathbf{W}\) of the volatility spillover network must be symmetric and can only contains non-negative elements, and the degree matrix \(\mathbf{D}\) needs to be invertible. Suppose these conditions are always satisfied for the applied volatility spillover network construction methods. A rolling-window procedure is employed where the half window length, denoted as \(\tau\), is set to approximately match the number of half-year common trading days across the selected stock market indices. For a given center date \(t_0 \;(t_0 > \tau)\), the window extracts a one-year segment of RV data \(\{\mathbf{v}_t \in \mathbb{R}^N\}_{t = t_0 - \tau}^{t_0 + \tau}\). It is the input to each graph construction method to generate the weight matrix \(\mathbf{W}_{t_0}\) and, subsequently, the normalized graph Laplacian matrix \(\mathbf{L}_{t_0}\). In addition, the yearly average vector RV, \(\overline{\mathbf{v}}_{t_0} = (\overline{v}_{t_0,1}, \overline{v}_{t_0,2}, \dots, \overline{v}_{t_0,N})^\top\), is computed for each market within the window. Thus, the average GSE over the window centered at \(t_0\) of the volatility spillover network can be calculated as: \[\label{VSP95signal95energy} E(\overline{\mathbf{v}}_{t_0}) = \overline{\mathbf{v}}_{t_0}^\top \mathbf{L}_{t_0} \overline{\mathbf{v}}_{t_0},\tag{5}\] which can be regarded as a function of \(t_0\). For comparison purposes, the GSE \(E(\overline{\mathbf{v}}_{t_0})\) of each graph construction method is normalized to the interval \([0,1]\) by dividing its all-time maximum value.

As ensuring the invertibility of the degree matrix is straightforward, the following discussion focuses on the conditions that the weight matrix must satisfy to produce a non-negative GSE. The Pearson correlation matrix method, the GLASSO precision matrix method and the DY-symmetrization method yield symmetric weight matrices that contain non-negative values. However, the DY-magnetic-Laplacian method produces an asymmetric weight matrix. Thus, additional operations are needed in this method to ensure a valid normalized Laplacian matrix for the corresponding GSE calculation. The normalized magnetic Laplacian \(\mathbf{L}_m^{(q)} \in \mathbb{C}^{N \times N}\) is calculated based on the asymmetric \(\mathbf{W}^{DY}\) to capture the directional information [12]. The formulation process is detailed in Section 3.2. The GSE value can also be measured for the magnetic Laplacian \(\mathbf{L}_m^{(q)}\) through Equation 5 with \(\mathbf{L}_{t_0} = \mathbf{L}_{m,t_0}^{(q)}\).

Thus, the GSE of the DY–magnetic Laplacian method can be measured consistently with the other construction methods, which ensures comparability, and it remains a positive real number. Relevant proofs are given in Appendix 7.2 and Appendix 7.3.

The depiction of the rolling-window RV GSE \(E(\overline{\mathbf{v}}_{t_0})\), calculated as Equation 5 for each graph construction method, is displayed in Appendix 9. During turbulent market periods (e.g., financial crises), the cross-market volatility linkages become stronger [10], [11]. The volatility spikes tend to be synchronous across markets, but with different magnitudes due to the distinct nature across markets. This is expected to lead to larger differences within RV graph signals \(\overline{\mathbf{v}}_{t_0}\), which causes higher GSE levels. On the other hand, during calm market conditions, cross-market volatility interrelations are weaker [10], [11], and volatility levels are generally low and change slowly for various markets. This results in a smoother RV signal and, thus, lower RV GSE levels. It is intuitive that the RV GSE of the DY-magnetic-Laplacian volatility network best reflects the idea described above, and can be regarded as the most informative graph about the volatility interrelationships because it successfully identifies the significant financial crises and stable periods.

In addition, by producing an asymmetric weight matrix of the volatility spillover network, the DY-magnetic-Laplacian method can capture the directional transmission of volatility. Influential stock markets tend to transmit more volatility to impressionable stock markets, whereas impressionable stock markets have less influence on influential stock markets. Thus, it is chosen to construct the volatility spillover network for the proposed RV forecasting model. In Section 4.4, we provide further empirical support for the choice.

3.2 sec:GSP32Based32on32the32Magnetic32Laplacian↩︎

Since the DY-magnetic-Laplacian construction method produces the most informative volatility spillover graph, the next step is to harness its structure to enhance the accuracy of RV forecasts. According to [12], while spatial GNNs can be readily adapted to the directed volatility spillover graph due to their message-passing mechanism, they also have certain limitations. Relying on the neighborhood information aggregation and transformation, spatial methods typically aggregate information from outgoing neighbors of a node. However, they place insufficient emphasis on the information from incoming neighbors. Thus, they may lose important volatility spillover information, which limits their ability to capture the volatility interrelationships. In contrast, the magnetic Laplacian is specifically designed to account for both incoming and outgoing connections, which allows it to preserve directional information and better capture the full structure of volatility spillovers [12].

Suppose a directed graph \(\mathcal{G}\) has \(N\) nodes, which means \(\mathbf{W} \in \mathbb{R}^{N \times N}\) and \(\mathbf{W} \neq \mathbf{W}^\top\). In order to define the normalized magnetic Laplacian \(\mathbf{L}_m^{(q)}\), the symmetric variant of weight matrix \(\mathbf{W}^s\) and the corresponding degree matrix \(\mathbf{D}^s\) are firstly defined in Equation 6 and Equation 7 , respectively, \[\label{A94s} \mathbf{W}^s = \frac{1}{2} (\mathbf{W} + \mathbf{W}^\top);\tag{6}\] \[\label{D94s} \mathbf{D}^s = \text{diag}(\sum_{j = 1}^{N} \mathbf{W}^s_{1j}, \ldots, \sum_{j = 1}^{N} \mathbf{W}^s_{Nj}).\tag{7}\] Besides, to capture the directional information, the phase matrix \(\boldsymbol{\Theta}^{(q)}\) is constructed as follows, \[\label{Theta} \boldsymbol{\Theta}^{(q)} = 2 \pi q (\mathbf{W} - \mathbf{W}^\top),\tag{8}\] where \(q\) is a non-negative hyperparameter that determines how to handle the directionality of the graph. Use \(\mathrm{i}\) to represent the imaginary unit. Based on \(\mathbf{W}^s\) and \(\boldsymbol{\Theta}^{(q)}\), the normalized magnetic Laplacians is formulated as, \[\label{normalized32magnetic32Laplacian} \mathbf{L}_m^{(q)} = \mathbf{I} - ((\mathbf{D}^s)^{-\frac{1}{2}} \mathbf{W}^s (\mathbf{D}^s)^{-\frac{1}{2}}) \odot \exp(\mathrm{i} \boldsymbol{\Theta}^{(q)}),\tag{9}\] where \(\exp(\cdot)\) is the exponential transformation applied to each element of the input matrix. \(q\) determines how the directional information is perceived and processed. For example, if \(q = 0\), \(\mathbf{L}_m^{(0)} = \mathbf{I} - ((\mathbf{D}^s)^{-\frac{1}{2}} \mathbf{W}^s (\mathbf{D}^s)^{-\frac{1}{2}})\). The directionality of the graph is eliminated through symmetrizing the asymmetric weight matrix. For more discussions on how different \(q\) values fit different graph motifs and reveal corresponding topological features, see [19]–[22]. In this study, the DY-magnetic-Laplacian construction method requires \(q>0\) to utilize the direction of the volatility spillover effect.

According to Theorem 2, the normalized magnetic Laplacian matrix \(\mathbf{L}_m^{(q)}\) is Hermitian. Besides, it is positive semi-definite [12]. The spectral theorem indicates that its eigendecomposition yields an orthonormal basis of complex eigenvectors with real and non-negative eigenvalues, presented as follows, \[\label{normalized32magnetic32Laplacian32eigendecomposition} \mathbf{L}_m^{(q)} = \mathbf{U}_m \boldsymbol{\Lambda}_m \mathbf{U}_m^{\dagger},\tag{10}\] where the two complex matrices \(\mathbf{U}_m\) and \(\mathbf{U}_m^{\dagger}\) are unitary inverses of each other, i.e., \(\mathbf{U}_m \mathbf{U}_m^{\dagger} = \mathbf{U}_m^{\dagger} \mathbf{U}_m = \mathbf{I}\), and the eigenvalues in \(\boldsymbol{\Lambda}_m\) are real, non-negative numbers. The columns of \(\mathbf{U}_m\) are eigenvectors of \(\mathbf{L}_m^{(q)}\), and they form the orthonormal graph Fourier basis. Meanwhile, \(\mathbf{U}_m^{\dagger}\) is the conjugate transpose of \(\mathbf{U}_m\), and the rows of \(\mathbf{U}_m^{\dagger}\) constitute another orthonormal basis dual to the graph Fourier basis in \(\mathbb{C}^N\). Elements in \(\boldsymbol{\Lambda}_m\) represent the frequency of the corresponding graph Fourier modes.

Since the DY-magnetic-Laplacian method is most informative about the inter-market volatility relationships and market volatility conditions, the orthonormal graph Fourier basis of the corresponding normalized magnetic Laplacian matrix is expected to construct a space where the representation of RV data becomes highly effective in capturing the volatility dynamics. To obtain such RV representation in the graph spectral domain, the Graph Fourier Transformation (GFT) and its inverse (IGFT) for the RV signal \(\mathbf{v}_t \in \mathbb{R}^N\) \((t = 1, \ldots, T)\) based on the normalized magnetic Laplacian are, respectively, defined as follows, \[\begin{align} \widetilde{\mathbf{v}}_t &= \mathbf{U}_m^{\dagger} \mathbf{v}_t; \tag{11}\\ \mathbf{v}_t &= \mathbf{U}_m \widetilde{\mathbf{v}}_t. \tag{12} \end{align}\]

3.3 Model Formulation↩︎

In this section, we introduce the proposed GSP-HAR model, which combines graph signal processing techniques and the HAR framework to capture both temporal and cross-market volatility dynamics. The overall workflow and key components of the model are illustrated in Figure 3, providing a structured overview of how volatility spillover information is incorporated into the forecasting process.

Now we present the detailed step-by-step development of the proposed framework. In the first step, the weight matrix of the volatility spillover graph is produced under the DY framework and is denoted as \(\mathbf{W}^{DY}\). The formulation process is illustrated in Appendix 8.3. Based on \(\mathbf{W}^{DY}\), the corresponding normalized magnetic Laplacian matrix is calculated based on the process described from Equation 6 to Equation 9 with \(\mathbf{W} = \mathbf{W}^{DY}\). To effectively capture the directional volatility spillover effect between global stock markets, the hyperparameter \(q\) of the magnetic Laplacian is positive, and its value is determined through the hyperparameter tuning process. Since \(q > 0\), the eigendecomposition of the magnetic Laplacian \(\mathbf{L}_m^{(q)}\) returns two complex matrices, \(\mathbf{U}_m\) and \(\mathbf{U}_m^{\dagger}\), and one eigenvalue matrix, \(\boldsymbol{\Lambda}_m\), as shown in Equation 10 . Thus, the historical RV features \((\mathbf{v}_d, \mathbf{v}_w, \mathbf{v}_m) \in \mathbb{R}^{N \times 3}\) can be projected to the spectral domain through the GFT: \[\label{GSP-HAR-GFT} (\widetilde{\mathbf{v}}_d, \widetilde{\mathbf{v}}_w, \widetilde{\mathbf{v}}_m) = \mathbf{U}_m^{\dagger} (\mathbf{v}_d, \mathbf{v}_w, \mathbf{v}_m).\tag{13}\] Here, \((\widetilde{\mathbf{v}}_d, \widetilde{\mathbf{v}}_w, \widetilde{\mathbf{v}}_m)\) is a complex matrix. Its real part is denoted as \((\text{Re}(\widetilde{\mathbf{v}}_d), \text{Re}(\widetilde{\mathbf{v}}_w), \text{Re}(\widetilde{\mathbf{v}}_m))\), and its imaginary part is represented as \((\text{Im}(\widetilde{\mathbf{v}}_d), \text{Im}(\widetilde{\mathbf{v}}_w), \text{Im}(\widetilde{\mathbf{v}}_m))\). The spectral RV representation embeds the volatility interrelationships, which are encoded in the informative volatility spillover network based on the DY-magnetic-Laplacian method, into the historical RV pattern via the GFT process.

In the spectral domain, the GSP-HAR model leverages the HAR model as a linear data-driven filter within different parts of the basis. In the spectral domain, the HAR model filter takes the spectral RV representations as input signals and outputs new signals. The filtering processes for the real part and the imaginary part are independent, and their coefficients are drawn from two HAR filters so that the model can be more expressive. The filtering process of the two parts is formulated as: \[\label{real32and32imag32HAR32fixed} \begin{align} \text{Re}(\widetilde{\mathbf{v}}_t) =&\; \alpha^{Re} + \beta_{d}^{Re} \text{Re}(\widetilde{\mathbf{v}}_d) + \beta_{w}^{Re} \text{Re}(\widetilde{\mathbf{v}}_w) + \beta_{m}^{Re} \text{Re}(\widetilde{\mathbf{v}}_m);\\ \text{Im}(\widetilde{\mathbf{v}}_t) =&\; \alpha^{Im} + \beta_{d}^{Im} \text{Im}(\widetilde{\mathbf{v}}_d) + \beta_{w}^{Im} \text{Im}(\widetilde{\mathbf{v}}_w) + \beta_{m}^{Im} \text{Im}(\widetilde{\mathbf{v}}_m). \end{align}\tag{14}\] The new signal \(\widetilde{\mathbf{v}}_t \in \mathbb{C}^N\) and \(\widetilde{\mathbf{v}}_t = \text{Re}(\widetilde{\mathbf{v}}_t) + \mathrm{i} \text{Im}(\widetilde{\mathbf{v}}_t)\). The two sets, \(\{\alpha^{Re}, \beta_{d}^{Re}, \beta_{w}^{Re}, \beta_{m}^{Re}\}\) and \(\{\alpha^{Im}, \beta_{d}^{Im}, \beta_{w}^{Im}, \beta_{m}^{Im}\}\), represent independent collections of HAR filter coefficients, all of which are real-valued. The new signal is projected back to the spatial domain through the IGFT: \[\label{fixeC-GSP-HAR-IGFT} \breve{\mathbf{v}}_t = \mathbf{U}_m \widetilde{\mathbf{v}}_t,\tag{15}\] where \(\breve{\mathbf{v}}_t \in \mathbb{C}^N\) and \(\breve{\mathbf{v}}_t = \text{Re}(\breve{\mathbf{v}}_t) + \mathrm{i} \text{Im}(\breve{\mathbf{v}}_t)\). In the spatial domain, the information of the real and imaginary parts of \(\breve{\mathbf{v}}_t\) is aggregated by a \(3\) layer shallow neural network \(NN(\cdot)\) to output the final forecast \(\widehat{\mathbf{v}}_t \in \mathbb{R}^N\): \[\label{fixeC-GSP-HAR-NN} \widehat{\mathbf{v}}_t = NN(\text{Re}(\breve{\mathbf{v}}_t), \text{Im}(\breve{\mathbf{v}}_t)).\tag{16}\] The neural network can merge the volatility information in the real and imaginary parts in a nonlinear way, which effectively enhances the expressive power of the proposed GSP-HAR model. As noted by [8], such nonlinearity is crucial for accurately capturing volatility spillover effects and for improving the precision of RV forecasts.

4 Experiments↩︎

Through the intuitive GSE analysis discussed in previous sections, the DY-magnetic-Laplacian construction method is selected to construct the volatility spillover network due to its informativeness about market conditions. GSP techniques are subsequently leveraged to build a forecasting model specifically tailored to the selected graph construction method. To support this design, experiments are conducted using RV time series data from \(24\) stock market indices. The main experiments assess the performance of the proposed GSP-HAR model by comparing it with baseline models to demonstrate its validity and effectiveness. Additional experiments further confirm that the model achieves superior forecasting performance when paired with the DY-magnetic-Laplacian construction method, compared to alternative approaches.

Moreover, as shown in Figure 2, the GSE profile of the DY-magnetic-Laplacian method exhibits the potential to distinguish different financial market conditions. Building on this observation, an intuitive analysis is conducted to assess whether the GSP-HAR model, which leverages the informative DY-magnetic-Laplacian construction method, can provide early warnings of market turbulence.

4.1 Model Evaluation Tools↩︎

Various evaluation criteria are applied to evaluate the effectiveness of the proposed model. The evaluation process is conducted on an index-wise basis, as different markets operate under different economic conditions and investor behaviors. Especially, the accuracy of out-of-sample RV forecasts is measured through the commonly used Mean Squared Error (MSE) and Mean Absolute Error (MAE) criteria. Their calculation detail are shown in Equation 17 in Section 4.3. Both evaluation criteria are employed to ensure that the models perform consistently, thereby providing a more robust assessment of their forecast accuracy. In addition, the Model Confidence Set (MCS, [23]) test is applied to evaluate all tested models to determine a subset of models that have statistically superior forecasting accuracy (based on MSE or MAE) for each stock market index. Similar to the input settings of the HAR models, the forecasting windows are set as \(h = 1, 5, 22\), corresponding to the short-term (daily), mid-term (weekly) and long-term (monthly) forecasting. The software and hardware configurations are:

• Software: Python 3.10.14; NumPy 1.26.4; Statsmodels 0.14.2; PyTorch 2.2.2+cu121.

• Hardware: Operating system: Windows 11; CPU: 13th Gen Intel(R) Core(TM) i9-13900HX; GPU: NVIDIA GeForce RTX 4080 Laptop GPU.

4.2 Benchmark Datasets and Baseline models↩︎

The benchmark dataset, which is downloaded from Oxford-Man Institute’s realized library, consists of the RV data of \(24\) stock market indices with around \(3500\) common trading days from May 2002 to June 2022. Among them, some indices are commonly used indices in previous RV forecast research (e.g., [16]). The daily RV data is produced based on the \(5\)-minute intraday high frequency returns. Given that RV forecast models are trained on RV observations on common trading days of the included stock market indices, the selection of these \(24\) indices is intentionally focused on markets with more overlapping trading schedules. This approach ensures that the RV time series can be more continuous with more synchronized trading days across different markets, thus improving the data usage efficiency and enhancing the reliability and comparability of the RV forecasts produced by different models. The first \(70\%\) of data (from May 2002 to Feb 2016) are the in-sample data, whereas the last \(30\%\) of data (from March 2016 to June 2022) are the out-of-sample data. Besides, the RV is scaled up by multiplying \(100\). A descriptive analysis of the scaled RV data for each stock market index is reported in Table 7 in Appendix 10.

Figure 4 shows the volatility spillover network for \(7\) selected stock market indices for demonstration and clarity. The nodes are the abbreviation of the selected indices, and the width of edges is proportional to the scale of the volatility spillover effect calculated under the DY framework. The arrows of the edges indicate the direction of the volatility spillover effect.

The baseline models include the HAR, VHAR, HAR-KS and GNN-HAR models. Although the exact construction of some baseline models is not available, they are formulated according to the corresponding research papers. For the GNN-HAR model, different values are tried for the hyperparameter \(L\), and only the best performance is reported in each RV forecasting task.

4.3 Experimental Results↩︎

The general index-wise performance of each model is measured by the MSE and MAE calculated as follows: \[\label{MSE95MAE95loss95cal} \begin{align} \text{MSE}_i =&\; \frac{1}{T_{out}} \sum_{t=1}^{T_{out}}(\widehat{\mathbf{v}}_{i,t} - \mathbf{v}_{i,t})^2;\\ \text{MAE}_i =&\; \frac{1}{T_{out}} \sum_{t=1}^{T_{out}}|\widehat{\mathbf{v}}_{i,t} - \mathbf{v}_{i,t}|, \end{align}\tag{17}\] where \(\text{MSE}_i\) and \(\text{MAE}_i\) is the MSE and MAE of the \(i_{\text{th}}\) stock market index, respectively. \(T_{out}\) denotes the duration of the out-of-sample dataset. \(\widehat{\mathbf{v}}_{i,t}\) refers to the RV forecasting value of the \(i_{\text{th}}\) index at time \(t\), and \(\mathbf{v}_{i,t}\) is the corresponding true RV value. The index-wise MSE results, MAE results and the MCS test results for different models across various forecasting horizons (\(h = 1, 5, 22\)) are presented from Table 1 to Table 3.

3pt

3pt

3pt

According to the MSE and MAE results in the tables above, the GSP-HAR model improves the RV forecasting performance across all three forecast horizons. It consistently achieves the lowest MSE and MAE scores for more than half of the stock market indices in the short-term (\(h = 1\)), mid-term (\(h = 5\)) and long-term (\(h = 22\)) RV forecasting tasks. In addition, the MCS test is conducted to assess whether these improvements are statistically significant. Across different evaluation criteria and forecasting horizons, the GSP-HAR model remains in the MCS for more than \(20\) stock market indices. On the other hand, other models appear in the MCS less frequently. These imply the consistent and significant superiority of the GSP-HAR model in different forecasting horizons. Compared to linear HAR models, such as HAR-KS and VHAR, the GSP-HAR is informative of the volatility interrelationships and is able to capture nonlinear volatility spillover effects, which is beneficial for accurate RV forecasting [8].

In contrast to the spatial GNN-HAR model, which faces challenges in capturing directional spillover information [12], the spectral GSP-HAR model leverages GFT techniques based on the magnetic Laplacian to effectively model directional volatility spillovers. Besides, while adding additional GNN layers to the spatial GNN-HAR model yields little improvement in RV forecasting accuracy, consistent with the findings of [8], the GSP-HAR model attains optimal performance by incorporating a global perspective on the volatility spillover effect. This contrast arises because multi-layer spatial GNNs aggregate information sequentially from local neighborhoods via multi-hop message passing, whereas GSP methods transform signals into the spectral domain, where the entire network structure is encoded and global interdependencies are effectively captured.

4.4 GSP-HAR Forecasting Results with Different Laplacian Construction↩︎

Although the RV GSE analysis in previous sections indicates that the DY-magnetic-Laplacian method most effectively captures volatility interrelationships by reflecting market conditions, its superior performance in RV forecasting relative to other construction methods has yet to be empirically validated. In this section, the impact of different volatility network construction and the associated Laplacian formulation on GSP-HAR RV forecasting accuracy is analyzed. As shown in Appendix 8, there are three widely-applied, data-driven volatility spillover network formulation methods of interest: the Pearson correlation matrix, the GLASSO precision matrix and the DY framework.

According to [12], the magnetic Laplacian is a unified GSP framework which can be applied to both undirected and directed graphs by controlling the hyperparameter \(q\). When \(q = 0\), the volatility spillover is regarded as undirected, and the magnetic Laplacian matrix is identical to the common Laplacian matrix based on the symmetric weight matrix. On the other hand, the directional transmission of volatility can be processed when \(q > 0\). Thus, the GSP-HAR model can handle all different volatility spillover network construction methods by choosing the appropriate \(q\) value without changing the model architecture. The notations of GSP-HAR models with different graph inputs are summarized below:

• ‘GSP-HAR-P’: the weight matrix is symmetric: \(\mathbf{W} = \mathbf{W}^{P}\). The hyperparameter \(q\) in the GSP-HAR model is set to \(0\): \(q = 0\).

• ‘GSP-HAR-GL’: the weight matrix is symmetric: \(\mathbf{W} = \mathbf{W}^{GL}\). The hyperparameter \(q\) in the GSP-HAR model is set to \(0\): \(q = 0\).

• ‘GSP-HAR-DY-sym’: the weight matrix is symmetric: \(\mathbf{W} = \frac{1}{2} (\mathbf{W}^{DY} + (\mathbf{W}^{DY})^\top)\). The hyperparameter \(q\) in the GSP-HAR model is set to \(0\): \(q = 0\).

• ‘GSP-HAR-DY-asym’: the weight matrix is asymmetric: \(\mathbf{W} = \mathbf{W}^{DY}\). The hyperparameter \(q\) in the GSP-HAR model is positive: \(q > 0\). This is the same GSP-HAR model tested in Section 4.3.

The calculation process of \(\mathbf{W}^{P}\), \(\mathbf{W}^{GL}\) and \(\mathbf{W}^{DY}\) are detailed in Appendix 8. Similar to Section 4.3, the index-wise MSE results, MAE results and the MCS test results for different volatility network construction methods and different forecasting horizons (\(h = 1, 5, 22\)) are presented from Table 4 to Table 6.

3pt

3pt

3pt

According to the tables above, the GSP-HAR model delivers consistent RV forecasting performance across different Laplacian constructions of the volatility spillover network, which underscores the robustness of the model design. Nevertheless, the GSP-HAR-DY-asym model, which is built on an asymmetric weight matrix derived from the DY framework and a magnetic Laplacian with a positive \(q\) to capture directional spillovers, emerges as the preferred variant. It achieves the lowest forecast errors and remains in the MCS for more than half of the stock market indices across different evaluation criteria and forecasting horizons. While the MCS test often includes most GSP-HAR variants at the 25% significance level due to their similar levels of loss, the GSP-HAR-DY-asym model consistently appears most frequently, highlighting its superior performance.

This advantage likely stems from two key factors. First, the DY-magnetic-Laplacian method produces RV GSE patterns that most accurately reflect underlying market conditions. Thus, it provides graph Fourier bases that yield highly effective spectral representations of historical RV patterns for forecasting. Second, by explicitly accounting for the directional volatility transmission, the DY-magnetic-Laplacian captures the inherent imbalance between influential and impressionable stock markets, which is crucial to accurately model volatility spillovers. Therefore, the DY-magnetic-Laplacian method can be regarded as the most informative for capturing volatility interrelationships and market conditions, offering the strongest forecasting capability within the GSP-HAR framework for accurate RV prediction.

4.5 GSP-HAR-Forecasted RV GSE as Indicators of Market Turbulence↩︎

This section examines the practical utility of the RV GSE forecasts generated by the proposed GSP-HAR model, which is identified as the most effective RV forecasting approach. Specifically, the model considered here is based on the asymmetric DY volatility network combined with the magnetic Laplacian matrix with \(q>0\). The analysis proceeds by intuitively comparing the GSP-HAR-forecasted RV GSE with the ground-truth RV GSE shown in Figure 2. The aim is not to evaluate whether the forecasted GSE perfectly replicates the ground-truth series. Forecasting accuracy has already been validated in Sections 4.3 and 4.4. Rather, the focus is on assessing whether the GSP-HAR-forecasted RV GSE can function as an early-warning indicator of global market conditions.

As discussed in Section 3.1, when the DY-magnetic-Laplacian method is applied, the RV GSE function is highly informative about market dynamics, distinguishing both turbulent and stable periods. Extending this property to forecasts, the GSE derived from GSP-HAR-forecasted RV is expected to provide early signals of potential market stress. To this end, the analysis employs long-term forecasts (\(h=22\)), as longer horizons are particularly valuable: unlike short- or mid-term forecasts, they allow decision-makers greater lead time to anticipate and respond to shifts in market conditions.

To visualize both the GSP-HAR-forecasted RV GSE and the ground-truth RV GSE functions, the procedure is similar to the one described in Section 3.1. A yearly rolling-window with center date point \(t_0\) and window size \(2 \tau\) can be denoted as \(\{\mathbf{v}_t \in \mathbb{R}^N\}_{t = t_0 - \tau}^{t_0 + \tau}\). The normalized magnetic Laplacian matrix \(\mathbf{L}_{m,t_0}^{(q)}\) used for calculating the GSE is constructed from the rolling window data \(\{\mathbf{v}_t\}_{t = t_0 - \tau}^{t_0 + \tau}\) under the DY-magnetic-Laplacian method. The GSP-HAR model is trained with the forecast horizon as \(22\), the Laplacian matrix \(\mathbf{L}_{m,t_0}^{(q)}\) and the training data \(\{\mathbf{v}_t \in \mathbb{R}^N\}_{t = t_0 - \tau}^{t_0 + \tau}\). The ground-truth RV is \(\overline{\mathbf{v}}_{t_0} = (\overline{v}_{t_0,1}, \overline{v}_{t_0,2}, \dots, \overline{v}_{t_0,N})^\top\), whereas the GSP-HAR-forecasted RV \(\hat{\mathbf{v}}_{t_0}\) is generated by the trained GSP-HAR model with input data \(\{\mathbf{v}_t \in \mathbb{R}^N\}_{t = t_0 + \tau - 22}^{t_0 + \tau}\). Use \(E(\overline{\mathbf{v}}_{t_0})\) and \(E(\hat{\mathbf{v}}_{t_0})\) to denote the ground-truth RV GSE and the GSP-HAR-forecasted RV GSE, respectively. Both \(E(\overline{\mathbf{v}}_{t_0})\) and \(E(\hat{\mathbf{v}}_{t_0})\) are calculated according to Equation 5 with \(\mathbf{L}_{t_0} = \mathbf{L}_{m,t_0}^{(q)}\), and are plotted at \(t = t_0\). Since \(\hat{\mathbf{v}}_{t_0}\) is the forecast RV, it is expected that the \(E(\hat{\mathbf{v}}_{t_0})\) should act as a leading indicator of the \(E(\overline{\mathbf{v}}_{t_0})\). For instance, when \(E(\hat{\mathbf{v}}_{t_0})\) rises sharply, \(E(\overline{\mathbf{v}}_{t_0})\) should follow with an increase after a certain lag. The GSE of the GSP-HAR-forecasted RV and the ground-truth RV are presented in Figure 5 below. They are scaled by the same factor so that the values of both functions are bounded within the interval \([0, 1]\).

As shown in Figure 5, the GSP-HAR-forecasted RV GSE emerges as a leading indicator of the ground-truth RV GSE, providing early-warning signals of heightened market stress. Notably, the forecasted energy function rises sharply prior to the GFC, with similar but less pronounced upward movements observed ahead of the Federal Reserve liftoff, the Brexit referendum, and the COVID-19 pandemic. Moreover, compared to the relative stability observed in tranquil periods, the forecasted energy function exhibits greater fluctuations both before and during episodes of market turbulence. These empirical patterns provide preliminary evidence that the GSP-HAR-forecasted RV GSE is both predictive and sensitive to potential financial risks. An increase in its level and volatility may signal forthcoming market instability. Thus, beyond its strength in producing accurate RV forecasts, the GSP-HAR model also demonstrates practical value as a tool for financial risk monitoring and early-warning analysis.

5 Conclusion↩︎

This study examines the interconnections of global financial markets by analyzing volatility spillovers. It creatively discovers the relationship between RV GSE and the financial market conditions, which can be leveraged to measure the effectiveness of different volatility spillover network construction methods. In addition, this study proposes a novel RV forecasting framework that combines the HAR model and GSP techniques based on the magnetic Laplacian to leverage the most informative volatility spillover network for more accurate RV forecasting. Results from the empirical tests on data from \(24\) global stock market indices show that, compared to baseline models, the proposed GSP-HAR model significantly improves the RV forecasting accuracy. Unlike existing RV forecast approaches, where the volatility spillover effect is regarded as additional variables representing spatial information in the HAR framework, the proposed GSP-HAR model uses the GFT technique to obtain effective RV representations in the spectral domain based on the volatility spillover network structure, and leverages the HAR model as a linear, data-driven filter to obtain the spectral volatility signal. Subsequently, the IGFT converts the filtered spectral RV signal back to the spatial domain, where a shallow neural network is applied to derive the final RV forecast in a nonlinear way. Thus, the volatility spillover effect is successfully integrated into the model’s design.

In this study, the exploration of market regime identification and turbulence forecasting is limited to an intuitive, preliminary analysis based on the GSE function. While the findings suggest that the RV GSE derived from the DY-magnetic-Laplacian method can distinguish between turbulent and stable periods, and that the GSE of GSP-HAR-forecasted RV is predictive of potential market stress, these insights remain exploratory. Extending the GSP-HAR framework into more systematic applications, such as financial early-warning systems and formal market regime analysis, offers potential direction for future research.

6 sec:Undirected32vs32Directed32Graphs↩︎

7 sec:Theorems32and32Corresponding32Proofs↩︎

In general, suppose the graph signal \(\mathbf{x} \in \mathbb{R}^N\) is defined on a \(N\)-node graph \(\mathcal{G} = (\mathcal{V}, \mathcal{E})\), whose weight matrix and degree matrix are \(\mathbf{W} \in \mathbb{R}^{N \times N}\) and \(\mathbf{D} \in \mathbb{R}^{N \times N}\), respectively. The weight matrix \(\mathbf{W}\) is assumed to contain non-negative values. And the degree matrix \(\mathbf{D}\) is assumed to be invertible. The normalized Laplacian matrix is formulated as: \(\mathbf{L} = \mathbf{I} - \mathbf{D}^{-\frac{1}{2}}\mathbf{W}\mathbf{D}^{-\frac{1}{2}}\). The normalized magnetic Laplacian \(\mathbf{L}_m^{(q)}\) formulation process is described from Equation 6 to Equation 9 .

7.1 sec:Non-negative32GSE↩︎

Considering the GSE defined in Equation 2 , we have the following theorem.

For any given undirected and non-negative weighted graph \(\mathcal{G}\) with no isolated nodes, i.e. \(\mathbf{W} = \mathbf{W}^\top\), \(\mathbf{W}_{ij}\geq 0\) \((i,j=1, ..., N)\), and \(d_i>0\) \((i=1, 2, ..., N)\), any graph signal \(\mathbf{x} \in \mathbb{R}^N\) on \(\mathcal{G}\) has non-negative energy: \(E(\mathbf{x}) \geq 0\).

Proof. According to Equation 2 , \(E(\mathbf{x}) = \mathbf{x}^\top \mathbf{L} \mathbf{x}\). The normalized Laplacian matrix \(\mathbf{L} = \mathbf{I} - \mathbf{D}^{-\frac{1}{2}}\mathbf{W}\mathbf{D}^{-\frac{1}{2}}\). Set \(\widetilde{\mathbf{W}} = \mathbf{D}^{-\frac{1}{2}}\mathbf{W}\mathbf{D}^{-\frac{1}{2}}\) and, correspondingly, \(\mathbf{L} = \mathbf{I} - \widetilde{\mathbf{W}}\). According to the assumptions, \(\widetilde{\mathbf{W}}\) is symmetric (\(\widetilde{\mathbf{W}} = \widetilde{\mathbf{W}}^\top\)) and only contains non-negative values (\(\widetilde{\mathbf{W}}_{ij} \geq 0\)). And it has a row sum of \(1\): \(\sum_{j=1}^N\widetilde{\mathbf{W}}_{ij} = 1\). Its column sum is also \(1\) due to the symmetry. The GSE calculation \(E(\mathbf{x}) = \mathbf{x}^\top \mathbf{L} \mathbf{x}\) can be expanded as: \[\begin{align} E(\mathbf{x}) =&\; \mathbf{x}^\top \mathbf{I} \mathbf{x} - \mathbf{x}^\top \widetilde{\mathbf{W}} \mathbf{x}\\ =&\; \sum_{i,j} \mathbf{I}_{ij}\mathbf{x}_i\mathbf{x}_j - \sum_{i,j} \widetilde{\mathbf{W}}_{ij}\mathbf{x}_i\mathbf{x}_j\\ =&\; \frac{1}{2}(\sum_i \mathbf{x}_i^2 + \sum_j\mathbf{x}_j^2 - 2\sum_{i,j}\widetilde{\mathbf{W}}_{ij}\mathbf{x}_i\mathbf{x}_j)\\ =&\; \frac{1}{2}(\sum_i \sum_j \widetilde{\mathbf{W}}_{ij} \mathbf{x}_i^2 + \sum_j \sum_i \widetilde{\mathbf{W}}_{ij} \mathbf{x}_j^2 - 2\sum_{i,j}\widetilde{\mathbf{W}}_{ij}\mathbf{x}_i\mathbf{x}_j)\\ =&\; \frac{1}{2}\sum_{i,j}\widetilde{\mathbf{W}}_{ij}(\mathbf{x}_i^2+\mathbf{x}_j^2-2\mathbf{x}_i\mathbf{x}_j)\\ =&\; \frac{1}{2}\sum_{i,j}\widetilde{\mathbf{W}}_{ij}(\mathbf{x}_i-\mathbf{x}_j)^2 \end{align}\] Here, \(\forall\, i,j \in \{1,\ldots,N\},\;\widetilde{\mathbf{W}}_{ij} \geq 0\), and \(\forall\, \mathbf{x} \in \mathbb{R}^N,\;(\mathbf{x}_i-\mathbf{x}_j)^2 \geq 0\). Hence, \(\forall\, \mathbf{x} \in \mathbb{R}^N,\;E(\mathbf{x}) = \mathbf{x}^\top \mathbf{L} \mathbf{x} \geq 0\). ◻

Remark 1: This also indicates that, when the conditions stated at the beginning of this section are satisfied, the normalized Laplacian matrix is positive semi-definite: \(\mathbf{L} \succeq 0\). Besides, when \(\mathbf{W} = \mathbf{W}^\top\), both the weight matrix \(\mathbf{W}\) and the degree matrix \(\mathbf{D}\) are symmetric. The normalized Laplacian matrix \(\mathbf{L}\) is also symmetric. In addition, the equality \(E(\mathbf{x}) = 0\) is achieved when \(\mathbf{x}_i = \mathbf{x}_j\), which indicates that the node feature \(\mathbf{x}_i\) is constant over each graph component. A graph component is a maximal connected subgraph in which each node is reachable from others and no additional nodes are included.

7.2 sec:Hermitian32Magnetic32Laplacian↩︎

For any directed graph \(\mathcal{G}\), its weight matrix \(\mathbf{W} \in \mathbb{R}^{N \times N}\) may be asymmetric: \(\mathbf{W} \neq \mathbf{W}^\top\). In this case, the normalized magnetic Laplacian can be defined as \(\mathbf{L}_m^{(q)} \in \mathbb{C}^{N \times N}\) by the process from Equation 6 to Equation 9 .

For any graph \(\mathcal{G}\), the normalized magnetic Laplacian \(\mathbf{L}_m^{(q)} \in \mathbb{C}^{N \times N}\) is Hermitian, that is, \(\mathbf{L}_m^{(q)} = (\mathbf{L}_m^{(q)})^*\).

Proof. To prove that \(\mathbf{L}_m^{(q)} \in \mathbb{C}^{N \times N}\) is Hermitian: \(\mathbf{L}_m^{(q)} = (\mathbf{L}_m^{(q)})^*\), it suffices to show the following conjugate symmetry condition: \(\forall\, i,j, \,(\mathbf{L}_m^{(q)})_{ij} = \overline{(\mathbf{L}_m^{(q)})_{ji}}\).

According to Equation 8 , \((\boldsymbol{\Theta}^{(q)})^\top = 2 \pi q (\mathbf{W} - \mathbf{W}^\top)^\top = 2 \pi q (\mathbf{W}^\top - \mathbf{W}) = -\boldsymbol{\Theta}^{(q)}\), which indicates \(\boldsymbol{\Theta}^{(q)}_{ji} = -\boldsymbol{\Theta}^{(q)}_{ij}\). For each element in the element-wise exponential operation, \(\exp(\mathrm{i} \boldsymbol{\Theta}^{(q)}_{ij}) = \cos(\boldsymbol{\Theta}^{(q)}_{ij}) + \mathrm{i} \sin(\boldsymbol{\Theta}^{(q)}_{ij})\). \(\exp(\mathrm{i} \boldsymbol{\Theta}^{(q)})\) is a Hermitian matrix because: \[\label{Hermitian32exp40iota32Theta94q41} \begin{align} \exp(\mathrm{i} \boldsymbol{\Theta}^{(q)})_{ji} =&\; \exp(\mathrm{i} \boldsymbol{\Theta}^{(q)}_{ji})\\ =&\; \cos(\boldsymbol{\Theta}^{(q)}_{ji}) + \mathrm{i} \sin(\boldsymbol{\Theta}^{(q)}_{ji})\\ =&\; \cos(-\boldsymbol{\Theta}^{(q)}_{ij}) + \mathrm{i} \sin(-\boldsymbol{\Theta}^{(q)}_{ij})\\ =&\; \cos(\boldsymbol{\Theta}^{(q)}_{ij}) - \mathrm{i} \sin(\boldsymbol{\Theta}^{(q)}_{ij})\\ =&\;\overline{\exp(\mathrm{i} \boldsymbol{\Theta}^{(q)})_{ij}}. \end{align}\tag{18}\]

Set \(\widetilde{\mathbf{W}}^s = (\mathbf{D}^s)^{-\frac{1}{2}} \mathbf{W}^s (\mathbf{D}^s)^{-\frac{1}{2}}\), where \(\mathbf{W}^s\) and \(\mathbf{D}^s\) are calculated in Equation 6 and Equation 7 , respectively. \(\widetilde{\mathbf{W}}^s\) is symmetric: \(\widetilde{\mathbf{W}}^s = (\widetilde{\mathbf{W}}^s)^\top\). Set \(\mathbf{M} = \exp(\mathrm{i} \boldsymbol{\Theta}^{(q)})\), where \(\mathbf{M} = \text{Re}(\mathbf{M}) + \mathrm{i} \text{Im}(\mathbf{M})\), \((\text{Re}(\mathbf{M}))^\top = \text{Re}(\mathbf{M})\), and \((\text{Im}(\mathbf{M}))^\top = -\text{Im}(\mathbf{M})\) because \(\exp(\mathrm{i} \boldsymbol{\Theta}^{(q)})\) is Hermitian.

According to Equation 9 , each element of the normalized magnetic Laplacian is calculated as: \[\begin{align} (\mathbf{L}_m^{(q)})_{ij} =&\; \mathbf{I}_{ij} - (\widetilde{\mathbf{W}}^s_{ij}\odot\mathbf{M}_{ij})\\ =&\; \mathbf{I}_{ij} - (\widetilde{\mathbf{W}}^s_{ij} \text{Re}(\mathbf{M})_{ij} + \mathrm{i} \widetilde{\mathbf{W}}^s_{ij} \text{Im}(\mathbf{M})_{ij})\\ =&\; \mathbf{I}_{ij} - \widetilde{\mathbf{W}}^s_{ij} \text{Re}(\mathbf{M})_{ij} - \mathrm{i} \widetilde{\mathbf{W}}^s_{ij} \text{Im}(\mathbf{M})_{ij}\\ =&\; \mathbf{I}_{ji} - \widetilde{\mathbf{W}}^s_{ji} \text{Re}(\mathbf{M})_{ji} + \mathrm{i} \widetilde{\mathbf{W}}^s_{ji} \text{Im}(\mathbf{M})_{ji}\\ =&\;\overline{(\mathbf{L}_m^{(q)})_{ji}} \end{align}\]

Thus, the normalized magnetic Laplacian \(\mathbf{L}_m^{(q)}\) is proved to be Hermitian. ◻

7.3 sec:Real4432Non-negative32GSE32Based32on32Hermitian32Laplacian↩︎

The GSE based on the normalized magnetic Laplacian, calculated as \(E(\mathbf{x}) = \mathbf{x}^\top \mathbf{L}_m^{(q)} \mathbf{x}\), is a real and non-negative scalar for graph signal \(\mathbf{x} \in \mathbb{R}^N\).

Proof. As shown in Appendix 7.2, the normalized magnetic Laplacian \(\mathbf{L}_m^{(q)} \in \mathbb{C}^{N \times N}\) is Hermitian. Set \(\mathbf{L}_m^{(q)} = \text{Re}(\mathbf{L}_m^{(q)}) + \mathrm{i} \text{Im}(\mathbf{L}_m^{(q)})\). Correspondingly, \(E(\mathbf{x}) = \mathbf{x}^\top \text{Re}(\mathbf{L}_m^{(q)}) \mathbf{x} + \mathrm{i} \mathbf{x}^\top \text{Im}(\mathbf{L}_m^{(q)}) \mathbf{x}\). Since \(\mathbf{L}_m^{(q)}\) is Hermitian, \((\text{Re}(\mathbf{L}_m^{(q)}))^\top = \text{Re}(\mathbf{L}_m^{(q)})\) and \((\text{Im}(\mathbf{L}_m^{(q)}))^\top = -\text{Im}(\mathbf{L}_m^{(q)})\). Set \(s = \mathbf{x}^\top \text{Im}(\mathbf{L}_m^{(q)}) \mathbf{x}\) and \(s \in \mathbb{R}\). \[\begin{align} s =&\;s^\top \\ =&\; (\mathbf{x}^\top \text{Im}(\mathbf{L}_m^{(q)}) \mathbf{x})^\top \\ =&\; \mathbf{x}^\top (\text{Im}(\mathbf{L}_m^{(q)}))^\top \mathbf{x} \\ =&\; -\mathbf{x}^\top \text{Im}(\mathbf{L}_m^{(q)}) \mathbf{x} \\ =&\; -s \\ \end{align}\]

Given \(s = -s\) and \(s \in \mathbb{R}\), \(s = 0\). Hence, \(E(\mathbf{x}) = \mathbf{x}^\top \text{Re}(\mathbf{L}_m^{(q)}) \mathbf{x} + \mathrm{i} s = \mathbf{x}^\top \text{Re}(\mathbf{L}_m^{(q)}) \mathbf{x}\), which is a real scalar.

Similar to Appendix 7.2, set \(\widetilde{\mathbf{W}}^s = (\mathbf{D}^s)^{-\frac{1}{2}} \mathbf{W}^s (\mathbf{D}^s)^{-\frac{1}{2}}\). \(\widetilde{\mathbf{W}}^s\) is symmetric (\(\widetilde{\mathbf{W}}^s = (\widetilde{\mathbf{W}}^s)^\top\)) and only contains non-negative values (\(\widetilde{\mathbf{W}}^s_{ij} \geq 0\)). And it has a row sum of \(1\): \(\sum_{j=1}^N\widetilde{\mathbf{W}}^s_{ij} = 1\). Its column sum is also \(1\) due to the symmetry.

According to Equation 9 , each element of the normalized magnetic Laplacian \((\mathbf{L}_m^{(q)})_{ij} = \mathbf{I}_{ij}-\widetilde{\mathbf{W}}^s_{ij}\cos{(\boldsymbol{\Theta}_{ij})} - \mathrm{i} \widetilde{\mathbf{W}}^s_{ij} \sin{(\boldsymbol{\Theta}_{ij})}\). Hence, the element-wise representation of the real part of the normalized magnetic Laplacian is \(\text{Re}(\mathbf{L}_m^{(q)})_{ij} = \mathbf{I}_{ij}-\widetilde{\mathbf{W}}^s_{ij}\cos{(\boldsymbol{\Theta}_{ij})}\). Thus, the matrix form GSE calculation, \(E(\mathbf{x}) = \mathbf{x}^\top \text{Re}(\mathbf{L}_m^{(q)}) \mathbf{x}\), can be expanded as: \[\begin{align} E(\mathbf{x}) =&\; \sum_{i,j} \text{Re}(\mathbf{L}_m^{(q)})_{ij} \mathbf{x}_i \mathbf{x}_j \\ =&\; \sum_{i,j} \mathbf{I}_{ij}\mathbf{x}_i \mathbf{x}_j-\sum_{i,j}\widetilde{\mathbf{W}}^s_{ij}\mathbf{x}_i \mathbf{x}_j\cos{(\boldsymbol{\Theta}_{ij})} \end{align}\] Since \(\cos{(\boldsymbol{\Theta}_{ij})} \leq 1\), the above equations can be further simplified as: \[\begin{align} E(\mathbf{x}) &\geq \sum_{i,j} \mathbf{I}_{ij}\mathbf{x}_i \mathbf{x}_j-\sum_{i,j}\widetilde{\mathbf{W}}^s_{ij}\mathbf{x}_i \mathbf{x}_j\\ =&\; \frac{1}{2}(\sum_{i}\mathbf{x}_i^2 + \sum_{j}\mathbf{x}_j^2 - 2 \sum_{i,j}\widetilde{\mathbf{W}}^s_{ij}\mathbf{x}_i \mathbf{x}_j)\\ =&\; \frac{1}{2}(\sum_{i}\sum_{j}\widetilde{\mathbf{W}}^s_{ij}\mathbf{x}_i^2 + \sum_{j}\sum_{i}\widetilde{\mathbf{W}}^s_{ij}\mathbf{x}_j^2 - 2 \sum_{i,j}\widetilde{\mathbf{W}}^s_{ij}\mathbf{x}_i \mathbf{x}_j) \\ =&\; \frac{1}{2}\sum_{i,j}\widetilde{\mathbf{W}}^s_{ij}(\mathbf{x}_i^2 + \mathbf{x}_i^2 - 2 \mathbf{x}_i \mathbf{x}_j)\\ =&\; \frac{1}{2}\sum_{i,j}\widetilde{\mathbf{W}}^s_{ij}(\mathbf{x}_i - \mathbf{x}_i)^2 \end{align}\] Since \(\forall\, i,j \in \{1,\ldots,N\},\;\widetilde{\mathbf{W}}^s_{ij} \geq 0\), and \(\forall\, \mathbf{x} \in \mathbb{R}^N,\;(\mathbf{x}_i-\mathbf{x}_j)^2 \geq 0\), it follows that, \(\forall\, \mathbf{x} \in \mathbb{R}^N,\;E(\mathbf{x}) = \mathbf{x}^\top \text{Re}(\mathbf{L}_m^{(q)}) \mathbf{x} \geq 0\). ◻

Remark 2: It is important to prove that the GSE \(E(\mathbf{x})\) remains a real, non-negative scalar for the graph signal \(\mathbf{x} \in \mathbb{R}^N\) when the applied Laplacian matrix is Hermitian under the DY-magnetic Laplacian method. This means that the RV GSE analysis can be consistently applied to different volatility network construction methods. Thus, the resulting RV GSE functions are comparable.

In addition, it can be further proved that the normalized magnetic Laplacian matrix is positive semi-definite: \(\mathbf{L}_m^{(q)} \succeq 0\). For further details, readers are referred to [12].

8 sec:Volatility32Spillover32Network32Formulation↩︎

In this section, different volatility spillover network construction methods are discussed. The dataset settings are consistent with Section 3. The given RV time series data is denoted as \(\{\mathbf{v}_t \in \mathbb{R}^N \}_{t = 1}^T\), where each \(\mathbf{v}_t = (v_{1,t}, v_{2,t}, \dots, v_{N,t})^\top\) represents the RV observations of \(N\) assets (or markets) at time \(t\). Since the volatility spillover effect refers to cross-market or cross-asset transmission, the diagonal elements of the output weight matrices of different volatility spillover construction methods are set to \(0\) to exclude the own volatility component (i.e., the self-loop) from the spillover measure.

8.1 The Pearson Correlation Matrix↩︎

The Pearson correlation matrix is a symmetric matrix that measures the pairwise linear dependence between variables. Each element in the Pearson correlation matrix takes a value within \([-1, 1]\). A value of \(1\) indicates a perfect positive linear relationship, \(-1\) indicates a perfect negative linear relationship, and \(0\) indicates no linear correlation. The Pearson correlation matrix can be regarded as the weight matrix of the volatility spillover network \(\mathbf{W}^{P} \in \mathbb{R}^{N \times N}\) of \(N\) stock market indices and is computed as follows.

First, for each variable \(i \in \{1, \dots, N\}\), compute its sample mean and sample standard deviation in Equation 19 and Equation 20 , respectively. \[\label{series32sample32mean} \overline{v}_i = \frac{1}{T} \sum_{t=1}^T v_{i,t}.\tag{19}\] \[\label{series32sample32std} s_i = \sqrt{\frac{1}{T-1} \sum_{t=1}^T (v_{i,t} - \overline{v}_i)^2}.\tag{20}\] The Pearson correlation coefficient between variables \(i\) and \(j\) is then defined as \[\label{Pearson32correlation32coef} \rho_{ij} = \frac{\sum_{t=1}^T (v_{i,t} - \overline{v}_i)(v_{j,t} - \overline{v}_j)}{(T-1) \, s_i \, s_j}.\tag{21}\]

Repeating this computation in Equation 21 for all \(i,j \in \{1, \dots, N\}\) yields the symmetric correlation matrix \(\mathbf{W}^{P} = [\rho_{ij}]\), whose diagonal entries satisfy \(\mathbf{W}^P_{ii} = 0\).

8.2 The GLASSO Precision Matrix↩︎

The GLASSO precision matrix is an estimate of the inverse covariance (precision) matrix that encodes the conditional dependence structure among variables, which are assumed to follow the multivariate Gaussian distribution. Here, the covariance matrix is denoted as \(\boldsymbol{\Sigma}\) and the corresponding precision matrix is \(\boldsymbol{\Omega}\). In the precision matrix, a zero entry \(\boldsymbol{\Omega}_{ij}\) indicates conditional independence between the corresponding variables given all other variables.

Let \(\mathbf{S} \in \mathbb{R}^{N \times N}\) denote the empirical covariance matrix of the given multi-asset volatility time series \(\{\mathbf{v}_t \in \mathbb{R}^N \}_{t = 1}^T\). The GLASSO estimator of the precision matrix \(\widehat{\boldsymbol{\Omega}}\) solves the following optimization problem [24]: \[\label{GLASSO95precision95estimation} \widehat{\boldsymbol{\Omega}} = \arg\max_{\boldsymbol{\Omega} \succ 0} \{ \log\det(\boldsymbol{\Omega}) - \text{trace}(\mathbf{S}\boldsymbol{\Omega}) - \lambda \sum_{i \neq j} |\boldsymbol{\Omega}|_{ij} \},\tag{22}\] where \(\lambda > 0\) is the regularization parameter controlling the level of sparsity and \(\boldsymbol{\Omega} \succ 0\) enforces positive definiteness. According to [8], it is more beneficial to build the volatility spillover network weight matrix \(\mathbf{W}^{GL} \in \mathbb{R}^{N \times N}\) based on the estimated precision matrix in the following way: \[\label{GLASSO95adj95mx} \mathbf{W}^{GL}_{ij} = \begin{cases} 1, &\widehat{\boldsymbol{\Omega}}_{ij} \neq 0;\\ 0, &\widehat{\boldsymbol{\Omega}}_{ij} = 0 \text{or} i=j. \end{cases}\tag{23}\]

8.3 The DY Framework↩︎

The DY framework can quantify the magnitude and direction of spillovers among multiple volatility time series. It measures the spillover between time series based on the forecast error variance decomposition of a vector autoregressive (VAR) model. Suppose the given volatility time series \(\{\mathbf{v}_t \in \mathbb{R}^N \}_{t = 1}^T\) is covariance stationary. A \(p\)-order VAR model (VAR(\(p\))) for these time series can be represented as: \[\label{VAR-p} \mathbf{v}_t = \sum_{i=1}^p \phi_i \mathbf{v}_{t-i} + \boldsymbol{\epsilon}_t, \;\;\boldsymbol{\epsilon}_t \sim (0, \boldsymbol{\Sigma}_{\boldsymbol{\epsilon}}),\tag{24}\] where each \(\phi_i\) is a scalar coefficient. This model can be rewritten in a moving average (MA) form: \[\label{VAR-p-MA} \mathbf{v}_t = \sum_{i=1}^{\infty} \mathbf{A}_i \boldsymbol{\epsilon}_{t-i},\tag{25}\] where \({\mathbf{A}_i}\) defined recursively by \(\mathbf{A}_i = \sum_{j=1}^{p} \phi_j \mathbf{A}_{i-j}\) with the assumption that \(\mathbf{A}_0 = \mathbf{I}_{N \times N}\) and \(\mathbf{A}_i = 0\) for \(i < 0\). By performing a variance decomposition of the \(h\)-step-ahead forecast error, each element in the resulting matrix \(\boldsymbol{\theta}^g(h)\) given by: \[\label{var-decomp} \boldsymbol{\theta}_{ij}^g(h) = \frac{\sigma_{jj}^{-1}\sum_{t=0}^{h-1}(\mathbf{e}_i^\top \mathbf{A}_h \boldsymbol{\Sigma}_{\boldsymbol{\epsilon}} \mathbf{e}_j)^2}{\sum_{t=0}^{h-1}(\mathbf{e}_i^\top \mathbf{A}_h \boldsymbol{\Sigma}_{\boldsymbol{\epsilon}} \mathbf{A}_h^\top \mathbf{e}_i)},\tag{26}\] where \(h\) is the forecast horizon and \(\mathbf{e}_i \in \mathbb{R}^N\) is a selection vector with \(1\) in the \(i\)-th position and \(0\) elsewhere. \(\boldsymbol{\theta}_{ij}^g(h)\) is regarded as the pairwise directional volatility spillover from \(j\) to \(i\). To standardize matrix \(\boldsymbol{\theta}^g(h)\) such that each row sums to \(100\) to derive the spillover index in percentage, the elements in the normalized matrix \(\widetilde{\boldsymbol{\theta}}^g(h)\) can be defined as: \[\label{std-var-decomp} \widetilde{\boldsymbol{\theta}}_{ij}^g(h) = \frac{100 \;\boldsymbol{\theta}_{ij}^g(h)}{\sum_{j=1}^N \boldsymbol{\theta}_{ij}^g(h)}.\tag{27}\]

Here, in order to construct a more meaningful volatility spillover network weight matrix \(\mathbf{W}^{DY}\) under the DY framework, the standardized volatility spillover index matrix \(\widetilde{\boldsymbol{\theta}}^g(h)\) is transposed so that \(\mathbf{W}^{DY} = \widetilde{\boldsymbol{\theta}}^g(h)^\top\) and \(\mathbf{W}_{ij}^{DY}\) is the spillover index which measures the relative strength of volatility spillover from stock market index \(i\) to stock market index \(j\). Particularly, the diagonal elements are set to \(0\): \(\mathbf{W}^{DY}_{ii} = 0\). For notational simplicity, the dependence on \(h\) in \(\mathbf{W}^{DY}\) is omitted, although strictly it should be written as \(\mathbf{W}^{DY}(h)\). In the empirical experiment section, \(h\) is also the forecast horizon for forecasting RV so that the volatility spillover network construction is consistent with the RV forecast model.

Furthermore, the net pairwise volatility spillover from stock market index \(i\) to stock market index \(j\) is denoted as \(\overline{\mathbf{W}}^{DY}_{ij}\) and can be calculated based on \(\mathbf{W}^{DY}\): \[\label{net32pairwise32DY} \overline{\mathbf{W}}^{DY}_{ij} = \begin{cases} \mathbf{W}^{DY}_{ij} - \mathbf{W}^{DY}_{ji}, &\mathbf{W}^{DY}_{ij} \geq \mathbf{W}^{DY}_{ji};\\ 0, &\mathbf{W}^{DY}_{ij} < \mathbf{W}^{DY}_{ji}. \end{cases}\tag{28}\]

9 sec:RV32GSE32Plots↩︎

10 RV Time Series Data Summary Statistics↩︎

3pt

Table 7: Descriptive Statistics for the RV Dataset of Various Stock Market Indices

Index	Area	Mean	Std Dev	Skewness	Kurtosis
AEX	Netherlands	0.880	0.529	2.686	12.486
AORD	Australia	0.604	0.379	4.218	35.781
BFX	Belgium	0.816	0.446	2.604	12.192
BSESN	India	0.939	0.551	4.814	54.345
BVSP	Brazil	1.063	0.495	2.834	15.739
DJI	USA	0.808	0.547	3.201	19.406
FCHI	France	0.959	0.541	2.349	9.241
FTSE	UK	0.884	0.548	3.479	25.224
GDAXI	Germany	1.001	0.598	2.461	9.298
GSPTSE	Canada	0.669	0.459	3.518	21.735
HSI	China (Hong Kong)	0.816	0.383	3.214	20.980
IBEX	Spain	1.025	0.557	2.775	16.836
IXIC	USA	0.824	0.477	2.562	12.042
KS11	South Korea	0.815	0.417	2.299	10.896
KSE	Pakistan	0.832	0.513	2.542	10.339
MXX	Mexico	0.758	0.409	3.547	29.078
N225	Japan	0.824	0.432	3.053	20.018
NSEI	India	0.946	0.591	5.272	70.530
OSEAX	Norway	0.945	0.588	5.532	82.454
RUT	USA	0.762	0.451	2.953	15.755
SSEC	China	1.085	0.682	2.445	9.384
SSMI	Switzerland	0.766	0.471	4.253	33.192
SPX	USA	0.799	0.546	2.858	14.293
STOXX50E	Eurozone	1.033	0.611	2.646	12.626

References↩︎