A quantile regression model based on a linear regression equation estimates parameters for τ quantiles under the variable Y t by introducing an indicative function in the loss function [32], defined as Equation 1. min β ∑ t = 1 n ρ τ Y t − X t β , (1) where ρ τ z = z ⋅ τ − I z < 0 is the loss function at the quantile level of τ (known as the pinball loss function). Where I z < 0 is the indicative function, the value is 1 when the independent variable z < 0 ; otherwise, the value is 0 . However, this model only considers linear relationships between the variables, which cannot state the effect of non-linearity. For this reason, Taylor (2000) involved a neural network model called the quantile regression neural network approach ( Q R N N ) [33]. Cannon (2018) extended the Q R N N to the monotone composite neural network of quantile regression ( M C Q R N N ), which can mitigate the “quantile crossing” problem [34]. The comprehensive estimations of multi-quantile C o V a R can be obtained by adjusted M C Q R N N . Assuming the number of market indexes is N and their price returns are R t i , the conditional value at risk of market i in level of quantile q s can be obtained by h i R t − i , q s defined as ² Equation 2. h i R t − i , q s ≡ ∑ m = 1 M n ω m i ⋅ ψ e ω m i i ⋅ q s + ∑ j ≠ i ω m j i R t j + b m i + b i (2)

The difference between Q R N N and M C Q R N N is introducing quantile q s as an input variable with a positive weight of e ω m i i . In addition, the nonlinear activation function ψ ⋅ is assumed to be invariant and known. The parameter of each node m in the hidden layer consists of weights ω m j i j ≠ i and the intercept b m i , while the output layer parameters are ω m i and b i . These parameters are verified to be consistent and asymptotically normal under large sample and regularity conditions. Moreover, they converge to the true function at a certain rate [35, 36].

However, the loss function in the quantile regression is not differentiable everywhere. It limits the use of artificial neural networks’ regular algorithms (ANNs) in M C Q R N N . Therefore, it is necessary to adjust the form of the loss function. One approach is to add a Huber norm h u [37], which allows a smoothing approximation to the error term near the origin. Furthermore, h u is a hybrid L 1 / L 2 − n o r m , which makes it possible to use standard gradient-based optimization algorithms.

In addition, when the capacity of the neural network is large, it is prone to over-fitting problems. Choosing a modest neural network structure and hyperparameters is an effective approach often used in machine learning. In a single hidden layer network, the most important hyperparameter is the number of hidden nodes M n . Therefore, choosing the appropriate number of nodes can reduce the capacity of the neural network and avoid the over-fitting phenomenon. In addition, Bishop (1995) proposed alleviating this problem by introducing weight decay regularization [38]. Such regularization requires adding an additional penalty term to the weight parameter ω k , m h . Referring to model of Cannon [39, 40], the final estimator is set as Equation 3 included a quadratic penalty term. min h i 1 TS ∑ t = 1 T ∑ s = 1 S ρ ∑ t = 1 T ∑ s = 1 S ρ R t i ‐ h i R t ‐ i , q s + λ M n ∑ m = 1 M n ω m i 2 + λ M n N ∑ m = 1 M n e ω m ii 2 + λ M n N ∑ m = 1 M n ∑ j ≠ i ω m ji 2 , (3) where λ is the parameter that regulates the weight of the quadratic penalty term in the loss function. When λ is 0 , this regression is transformed to an ordinary M C Q R N N . In this paper, simple sampling is used to train neural networks. According to the experience of [40] and [16], we selected 50% of the samples from the sample period as the training set to train the loss function, which is Equation 3, in the MCQRNN model. Except for λ and M n hyperparameters, the intercept and weights of the neural network are trained. Otherwise, too large λ will lead to the loss of non-linear characteristics of the model, when the transfer function ψ is the sigmoidal hidden layer transfer function, such as hyperbolic tangent tanh . To balance the degree of the over-fitting and the prediction accuracy, λ is set to equal to 2 here. Furthermore, the optimization of the loss function as shown in Equation 3 can be achieved by using a quasi-Newton optimization algorithm, which is less complex and more appropriate for the computational complexity in this paper. The quantile regression neural network process is visualized in Figure 1.

The calibration details of M C Q R N N − Q Q − C o V a R are explained in this section. There are four steps involved in the systemic risk system calibration. The first step is the estimation of V a R based on C A V i a R . Next, the results are used to estimate C o V a R with M C Q R N N and Q Q for each country. In the third step, the composite risk spillover effects are calculated by resulting in an extreme risk spillover measure. Finally, the systemic risk measures are proposed based on the systemic risk network. The process of the risk modeling system is demonstrated in Figure 2 as follows:

Step 1: Estimation of value at risk with C A V i a R

Because it is challenging to select common macro-state variables for all indexes, the linear quantile regression of value at risk is no longer suitable for measuring the tail risk of stock markets. Hence, C A V i a R is adopted in the first step [41]. Given the asymmetric effects of the rise and fall of each index, A S − C A V i a R at the quantile level of p l can be estimated in the following Equation 4 [42] ³ : C A V i a R l , t j β = β 1 + β 2 C A V i a R l , t − 1 j β + β 3 R t − 1 j + + β 4 R t − 1 j − , (4) where R t − 1 j + is the absolute value of market index j ’s return when the lagged log-return is larger than zero, and the rest is 0 . R t − 1 j − is the absolute value when the lagged log-return is minus zero. β = β 1 , β 2 , β 3 , β 4 is the vector of the parameters. By utilizing the differential evolution algorithm [43] and the loss function similar to Equation 1, the C A V i a R l , t j is an appropriate risk indicator for individual stock market j at the quantile level of p l .

Step 2: Estimation of C o V a R with M C Q R N N − Q Q

First, the M C Q R N N method is adopted to estimate h ^ i R t ‐ i , q s based on Equation 3. Additionally, given the quantiles of stand-alone q 1 , q 2 , . . . , q S and quantiles of exposed risk p 1 , p 2 , . . . , p L , the quantile-on-quantile conditional value at risk ( Q Q − C o V a R ) can be obtained by the M C Q R N N . Q Q − C o V a R [16] is estimated similar to C o V a R [12] and T E N E T [44]. To embed the dependency among financial markets, the estimation of C o V a R by M C Q R N N is introduced. Following [16], the definition of C o V a R is adjusted as Equation 5 to adapt to a multivariate model. P r o b R t i < C o V a R l , t s , i | R t j = C A V i a R l , t j , ∀ j ≠ i = q s . (5)

Assume the C o V a R of market index i is predictable via function h i R t − i , q s , and other indexes’ current return R t j j ≠ i , h ^ i is the estimator of the function and can be trained by the M C Q R N N algorithm [34]. Therefore, the C o V a R in the condition of each risk state ( p l , q s ) can be obtained by Equation 6. C o V a R l , t s , i = h i C A V i a R l , t − i , q s . (6)

It should be noted that, in distinction to Keilbar and Wang, this paper estimates C o V a R l , t s , i for all quantiles of the whole range ⁴ . Thus, the C o V a R of each market index i can be expressed as a three-dimensional surface at each time point by setting the X and Y axes to denote different quantile levels p and q , respectively. The z-axis values are C o V a R l , t s , i . There is an advantage of reflecting both the non-linear relationship between C o V a R l , t s , i varying with the risk condition C A V i a R l , t j and the C o V a R at each of its own risk levels q .

Step 3: Calculation of composite risk spillover effects

To examine the margin impact of index j on the C o V a R l , t s , i , the partial derivative is taken, which is named d l s , t j i . According to the format of the function h ^ i estimated by M C Q R N N , the risk spillover from j to i is expressed as Equation 7, d ls , t j i = ∂ h ^ i ∂ R t j CAVia R l , t ‐ i , q s = ∑ m = 1 M n ω ^ m i ω ^ m ji ⋅ ψ ′ e ω ^ m ii ⋅ q s + ∑ j ≠ i ω ^ m ji CAVia R l , t j + b ^ m i , (7)

where the derivative of the transfer function is set as ψ ′ = 1 2 1 − tanh 2 z 2 , and all parameters are trained by M C Q R N N . This measure should be based on the accuracy of the C o V a R estimation, but the Q R N N , as a non-linear neural network model, is prone to over-fitting at a single quantile. Consequently, the M C Q R N N adopted in this paper is able to reduce the impact of potential fitting error at a single quantile on the overall risk spillover. Considering the risk spillovers at different quantiles, the composite risk spillover can be defined as δ t j i as Equation 8. δ t j i = 1 S L ∑ s = 1 S ∑ l = 1 L 1 + C A V i a R l , t j ) ⋅ 1 + C o V a R l , t s , i ) ⋅ d l s , t j i , i f   j ≠ i   0 ,                                                 i f   j = i , (8)

where the d l s , t j i is the absolute value of risk spillover, and the risk-weighted composite risk spillover from index j to i is defined as δ t j i . Because the absolute value of the risk indicators ( C A V i a R , C o V a R ) in the more extreme state is larger, the risk spillover in the extreme risk state couple p l , q s accounts for more weight in the aggregate indicator.

To reflect the sensitivity of average-level composite risk spillover to the two-sided quantiles, we decompose the δ ¯ j i = ∑ t = 1 T δ t j i into two partial spillover indicators δ s j i and δ l j i , as Equation 9. δ ¯ s j i = 1 T S ∑ t = 1 T   ∑ l = 1 L 1 + C A V i a R l , t j ⋅ 1 + C o V a R l , t s , i ⋅ d l s , t j i , i f   j ≠ i 0 , i f   j = i δ ¯ l j i = 1 T S ∑ t = 1 T   ∑ s = 1 S 1 + C A V i a R l , t j ⋅ 1 + C o V a R l , t s , i ⋅ d l s , t j i , i f   j ≠ i 0 , i f   j = i (9)

The partial spillover indicators are used to reflect the average impact of the two-sided risk state on the spillover from market j to market i in the sample period. The partial spillover indicator of p l side δ l j i is the composite spillover effect at the exposed risk state of p l , which reflects the sensitivity to the risk state of the entire system. Correspondingly, the δ s j i indicates the response of spillover effect to the change in stand-alone risk status q s . In addition, it is necessary to consider the spillover of market j to i under the most extreme conditions, that is, q 1 or p 1 . Hence, we consider the conditional partial spillover indicators under the fixed q 1 quantile or fixed p 1 as Equation 10. δ ¯ 1 s j i = 1 T ∑ t = 1 T d 1 s , t j i and δ ¯ l 1 j i = 1 T ∑ t = 1 T d l 1 , t j i . (10)

The non-linear characteristics of inter-market risk spillover can be analyzed by comparing the two types of partial spillover indices under different quantiles, and the mutation quantile of spillover can be captured. In addition, to reflect the extreme risk spillover, the δ ¯ 11 j i = ∑ t T d 11 , t j i is defined as spillover at a single quantile. If this indicator is calculated by Q R N N , it will be the same as the spillover defined by [18].

Last but not least, similar to the partial spillover indicators, four partial C o V a R terms can be defined as the average C o V a R when the one-side quantile takes different values under the condition of the other fixed, that is, C o V a R ¯ s i and C o V a R ¯ l i . To reflect the tail risk in the extreme condition, the traditional C o V a R estimated by Q R N N and M C Q R N N are also calculated as C o V a R ¯ 1 s i , C o V a R ¯ l 1 i . Taking the C o V a R ¯ s i and C o V a R ¯ 1 s i as an example, Equation 11 is consistent with the partial spillover δ ¯ s j i . C o V a R ¯ s i = 1 T L ∑ l = 1 L ∑ t = 1 T C o V a R l , t s , i C o V a R ¯ 1 s i = 1 T ∑ t = 1 T C o V a R l , t s , i | l = 1 . (11)

The quantile crossing problem in the estimation of C o V a R can be described and counted as Equation 12, ∃ s 1 < s 2 , C o V a R ¯ s 1 i < C o V a R ¯ s 2 i . (12)

For each market, the total number of “quantile crossing” problems can be accumulated by a monotonicity test. Because the C o V a R t is more adaptable in financial markets, we employ C o V a R 1 , t s , i , C o V a R ¯ 1 s i , and C o V a R ¯ s i to accumulate numbers of “quantile crossings” for each market.

Step 4: Network analysis for composite systemic risk

The average measures of composite systemic risk will be gained in the final step. First, because the risk spillover of market index j on market index i involves two quantiles p l and q s , the average indicators are composited and averaged. Referring to previous works [16, 45], the composite systemic suffering indicator Γ t i should be aggregated as Equation 13. The composite index is weighted by the own risk of each spillover emitter. This is because if market j has higher own risks, market i will bear more systemic risk spillover. Γ t i = 1 S L ∑ s = 1 S ∑ l = 1 L ∑ j ≠ i 1 + C A V i a R l , t j ⋅ d l s , t j i . (13)

Second, the composite systemic overflow indicator Φ t j can be defined as Equation 14. Φ t j = 1 S L ∑ s = 1 S ∑ l = 1 L ∑ i ≠ j 1 + C o V a R l , t s , i ⋅ d l s , t j i . (14)

To reveal the trend of the spillover effect, the total composite overflow indicator Π t can be calculated as Equation 15. Π t = 1 S L ∑ s = 1 S ∑ l = 1 L ∑ j = 1 N ∑ i ≠ j 1 + C A V i a R l , t j ⋅ 1 + C o V a R l , t s , i ⋅ d l s , t j i (15)

The total overflow indicator can be regarded as the weighted sum of Φ t j , which reflects the changing trend of risk spillover time of each market in the sample.

Lastly, the adjusted adjacency matrix is defined as Equation 16. A t = 0 δ ¯ t 12 ⋯ δ ¯ t 1 N δ ¯ t 21 0 ⋯ δ ¯ t 2 N ⋮ ⋯ ⋱ ⋮ δ ¯ t N 1 δ ¯ t N 2 ⋯ 0 . (16)

The adjusted adjacency matrix accounts are the risk spillover indicators. Systemic spillover effects are thus determined by the marginal effects of the M C Q R N N procedure, as well as by the V a R and C o V a R of the considered countries.

We select the stock market indices of N = 18 representative countries or regions from the Asia–Pacific Economic Cooperation (APEC) organization varying from January 2012 to December 2021 as sample data. Due to the weekday effect and excessive short-term volatility noise in daily frequency data, the daily returns of stock indexes taken from the WIND database are transformed into weekly data. The logarithmic return is calculated using the closing price on the last trading day of each week. Adopting weekly frequency data can effectively avoid the problem of time differences between the markets, and the data for each week can simply be treated as contemporaneous. The indexes and their abbreviations of regions for the 18 markets are shown in Table 1.

As shown in Table 2, the total number of observations is 9,396 because some markets have missing samples due to the holidays or other factors. The average weekly return of each market is all positive, which indicates that the indices prices of Asia–Pacific stock markets at the end of 2021 are higher than they were in 2011. It is worth noting that the minimum value of all samples is −20.13%, which is the weekly return of the Chile index in the fourth week of March 2020. The maximum value occurred in the next week, which is 15.81% in Japan. The World Health Organization recognized the COVID-19 outbreak as a global pandemic on 11 March 2020. After 2 weeks of declines, most global stock markets rebounded sharply in the last week of March. Except for the stock markets of Mexico, mainland China, and Hong Kong, all the financial markets have a kurtosis of 3 or more in their return distributions. The return curves show a sharp peak pattern, which indicates that the outliers are more dispersed. In addition, the skewness of all the markets in the sample is negative, indicating that the return series of each market is left skewness; that is, the probability of negative extreme value is higher than positive. The augmented Dickey–Fuller (ADF) value of each market return is negative and less than the test critical value at the 1% significant level, rejecting the null hypothesis of a unit root. All Jarque–Bera (JB) statistics are significant at the 1% level, which rejects the null hypothesis of Gaussian distribution for the market returns.

In this part, we analyze the average overflow of each market to the systemic Φ ¯ j , and their suffering Γ ¯ i , relying on the methodology in ⁸ Section 2.4.

The suffering indicators Γ ¯ i and the overflow indicators Φ ¯ j , which are respectively calculated by single- Q R N N , multi- Q R N N , single- M C Q R N N , and multi- M C Q R N N , are drawn as bar charts in Figure 11. As the legend shows, the suffering indicators Γ ¯ i are in light color and on the left side of each market, while the overflow indicators Φ ¯ j are darker and on the right side. In addition, the C o V a R calculated by Q R N N and M C Q R N N are illustrated as light and dark orange in the secondary axis. It is obvious that the light red bars stand out, indicating that China suffers the highest risk overflow when the Γ ¯ i calculated by Q R N N at a single quantile. However, when the algorithm is substituted by M C Q R N N with multiple quantiles, China is no longer the highest suffering market. In addition, the C o V a R obtained by Q R N N and M C Q R N N are relatively close, and the C o V a R of each market states no significant correlation with both Φ ¯ j and Γ ¯ i indicators.

Although Figure 11 shows differences in index calculations under the four algorithms, the comprehensive index is smaller than that calculated by a single sub-site. Whether under single or multi-component sites, the exponents obtained by the M C Q R N N algorithm are all smaller than those obtained by the Q R N N algorithm. In order to investigate whether the results of algorithms affect the ranking of systemic risk in each market, we also list the top ten markets of systemic risk index under various algorithms. As shown in Table 4, markets with higher risk overflow Φ ¯ j are HK, CA, US, and SG. This result is consistent with the global financial markets’ practical experiences. In contrast, JP, CN, HK, and PH suffer more systemic risk because of the higher Γ ¯ i . On the one hand, this may be related to the fact that these countries are more dependent on trade exports and have poor economic resilience, resulting in suffering more risk overflow. On the other hand, higher Γ ¯ i may also indicate that investors in these stock markets react strongly to the shocks. Compared with M C Q R N N , the Γ ¯ i of China is driven higher than it calculated by Q R N N , which means that the M C Q R N N method weakens the impact of the extreme condition.