Spatial Multivariate GARCH Models and Financial Spillovers

: We estimate the risk spillover among European banks from equity log-return data via Conditional Value at Risk (CoVaR). The joint dynamic of returns is modeled with a spatial DCC-GARCH which allows the conditional variance of log-returns of each bank to depend on past volatility shocks to other banks and their past squared returns in a parsimonious way. The backtesting of the resulting risk measures provides evidence that (i) the multivariate GARCH model with Student’s t distribution is more accurate than both the standard multivariate Gaussian model and the Filtered Historical Simulation (FHS), and (ii) the introduction of a spatial component improves the assessment of risk proﬁles and the market risk spillovers.


Introduction
The interconnectedness of risk between banks is an increasingly hot topic.In the last decades, several countries have simultaneously faced severe economic conditions with spillover effects of risk across the EU.Due to the direct and indirect links among the banks, the stand-alone measurement of a Value-at-Risk (VaR) of each bank cannot provide a comprehensive representation of the risk (Adrian and Brunnermeier 2014; Billio et al. 2012;Rahman 2014).
Recently, multivariate GARCH models have been playing a crucial role in estimating risk interconnectedness.The constant conditional correlation GARCH model (CCC-GARCH) proposed by Bollerslev (1990) is computationally less complex than other multivariate models (see, among others, Bollerslev et al. 1988;Diebold and Nerlove 1989;Engle et al. 1990).However, it does not capture the dynamic interactions between the volatilities.The BEKK model proposed by Engle and Kroner (1995) (the acronym stands for Baba, Engle, Kraft, and Kroner) allows for the dependence of conditional variances and covariance of one variable on the lagged values of another variable, so that spillovers in variances can be modeled.However, it is highly computationally intensive due to the large number of parameters.A more parsimonious model is the Dynamic Conditional Correlation GARCH (DCC-GARCH), introduced by Engle (2002), that introduces an autoregressive process for the conditional correlation matrix, allowing us to model its dynamics with only two parameters in addition to the ones of the CCC-GARCH model.
An approach for modeling explicitly the volatility interactions is the introduction of a spatial component that accounts for the effect of direct bilateral exposures, closeness, or similarities between different financial institutions.Borovkova and Lopuhaa (2012) adopted a spatial GARCH approach to handle the spillover effects where the spatial weights are obtained from the GDP data and alternatively from the market capitalization of the US and European countries' stock market and embedded in the Extended CCC-GARCH model (E-CCC), see Jeantheau (1998).As a result, they better capture the high kurtosis of squared returns.Keiler and Eder (2013) studied the systematic risk that integrates the interaction between the micro and macro stress situations as spatial econometrics parameters.Analogously, Chen (2017) showed that when the spatial weights are derived from credit rating downgrades, the multivariate spatial BEKK model can capture the spillover effects among the southern European stock index: Portugal, Italy, Ireland, Greece, and Spain (PIIGS).Zhang et al. (2018) applied the multivariate GARCH with a dynamic panel of spatial weight matrices based on the GDP.The work studies the countries' interconnectedness of returns and uses the estimated parameters to forecast the portfolio risk of six stock indices.
Our contribution is to introduce a dynamic conditional correlation GARCH (DCC-GARCH) model with a spatial component based on the credit exposure similarity among banks derived from the EU-wide stress test data.We add the spatial components into a DCC-GARCH model to investigate whether we can better capture the spillover effects thanks to this additional information and alternative distributional assumptions of the DCC-GARCH model.We discuss and implement the estimation of the model using both Gaussian and Student's t distributional assumptions.In particular, we estimate the individual risk via Value at Risk (VaR) and spillover risk via CoVaR.The results of the different models are compared, both amongst themselves and with the one obtained thanks to Filtered Historical Simulations (FHS).Finally, the results are backtested (Abad et al. 2014;Caporin 2008;Christoffersen and Pelletier 2004;Kupiec 1995) to evaluate the accuracy of the different VaR and CoVaR estimates, showing the superiority of the spatial DCC-GARCH model compared to the other models.
Finally, we point out that the applications of spatial DCC-GARCH models not only allow us to estimate CoVaR, but also permit investors to estimate more accurately the distribution of future returns of a set of assets to develop optimal portfolio strategies.The proposed framework has relevant applications for financial regulators interested in accurately measuring risk spillovers and systemic risk in financial systems, but it also caters to risk managers who want to measure the risk related to the interconnectedness among institutions.
The remainder of this paper is organized as follows.In Section 2, the spatial DCC-GARCH model and the estimation procedure are discussed, and the methodology for the financial application is presented.Section 3 presents the data and the empirical results, and in Section 4 we present our conclusions and discuss the results in relation to the literature.

Modelling and Inference
A common feature of financial time series is the presence of volatility clustering (see, e.g., Cont 2001). 1 Common tools used to address such features are Generalized Auto-Regressive Conditional Heteroscedasticity (GARCH) models (Bollerslev 1986), which generalize the ARCH models introduced by Engle (1982).Let r t be the return discretetime process with zero mean.The standardized disturbances ε t are independent and identically distributed (i.i.d.) with zero mean, E(ε t | ε t−1 , . ..) = 0, and unit variance, Var(ε t | ε t−1 , . ..) = 1.Then, the GARCH(p, q) process for return r t is defined as and where h t is the conditional variance, ω > 0, α k ≥ 0 and β k ≥ 0 ∀k.When studying spillover risk, it is natural to look for multivariate extensions of the GARCH model to characterize the joint evolution of stock returns.Before presenting the multivariate model it is useful to define the following quantities of interest.Assuming a market with N assets, then at time t = 1, . . ., T we have:  (Bollerslev et al. 1988;Ling and McAleer 2003) or the BEKK model, (Engle and Kroner 1995) have been extensively discussed in the literature.Using matrix notation, it is possible to characterize a multivariate GARCH as follows: and where H t , A k , B k are N × N matrices and ε t is an R N valued i.i.d.sequence of random variables with zero-mean and unit-variances (see Engle and Kroner (1995) for the restrictions required to ensure stationarity and positive semi-definiteness of the conditional covariance matrix).
Multivariate GARCH models have the drawback of having a large number of parameters, making the estimation complex and computationally challenging, hence these models are suitable only if the dimensionality N is very small.A solution to the dimensionality problem is to pose further restrictions on the multivariate process.A common restricted specification is the Constant Conditional Correlation model (CCC) proposed by Bollerslev (1990) that assumes that the conditional covariance matrix is constant over time, requiring focusing solely on the estimation of conditional variances.According to the CCC-GARCH model, Equation ( 1) is given by (3) and Equation (2) can be written as follows: where ω is the N × 1 dimensional vector of unconditional variances with ω ∈ R + , α k and β k are the N × 1 dimensional vector of ARCH and GARCH parameters of order q and p with α k,i ∈ R + 0 , β k,i ∈ R + 0 , and is the Hadamard product.The CCC-GARCH model assumes that the conditional covariance matrix, H t , can be factorized as where the correlation matrix is assumed to be constant throughout time (R t = R, ∀t) and the conditional standard deviation matrix D t is a diagonal matrix given by The generic element of conditional covariance matrix H t is constructed as where ρ ij = [R] ij is the constant conditional correlation coefficient between the ith and jth variables.
The multivariate GARCH model with a dynamic conditional correlation structure (DCC), introduced by Engle (2002), improves the dynamic relationship, assuming a timevarying correlation matrix as follows The dynamic correlation model allows R t to be time-varying, and its dynamics are modeled assuming a GARCH(1,1) process for the covariance of the standardized residuals.Hence R t is decomposed into where where γ and δ are ARCH parameters and GARCH parameters of the DCC model, respectively.By following the GARCH model from Equation (2), the generic element of the time-varying conditional covariance matrix of the standardized residuals [Q t ] ij = q ij,t can be expressed as The process is mean-reverting as long as 0 < δ < 1 and γ + δ < 1.In the particular case of γ + δ = 1, the process will follow the exponential smoother matrix of the standard residuals, as described in Engle (2002).Finally, the generic conditional correlation can be written into matrix form as in Equation ( 10).Substituting the conditional correlation matrix into Equation ( 9), the DCC is given by Restricted GARCH models beyond CCC and DCC-GARCH have been discussed by Caporin (2008) and Billio et al. (2021), who introduce spatial matrices within BEKK models for measuring risk spillover.In these approaches, the interaction components of the model are based on spatial weight matrices provided exogenously (for instance on the basis of geographical distances among assets, or some similarity metrics).These models allow easier and more accurate estimation by effectively imposing restrictions on the parameter space.An alternative approach to improve the estimation of multivariate GARCH models is to introduce sparsity in the parameter estimates by using an L 1 penalization, as suggested by Dhaene et al. (2022).

Spatial DCC-GARCH
In this work, we introduce a spatial extension of the DCC-GARCH model.The model is based on the approach of Borovkova and Lopuhaa (2012).In particular, to enrich the DCC-GARCH model with a spatial component we introduce a spatial matrix W into the vector of the conditional variances h t .The resulting conditional variance is expressed as where A 0 = (a 0,1 , . . ., a 0,N ) , A 1,k , A 2,k , B 1,k , and B 2,k are diagonal matrices.The term W = W ij is the weight matrix (i, j = 1, . . ., N) with N ∑ j=1 w ij = 1 and w ii = 0 ∀i, given by where The introduction of the spatial component results in two exogenous spatial variables in the conditional variance equation and two additional parameters a 2,i and b 2,i , which measure the influence of the aggregated lagged variances and squared returns of all the other assets.These two new variables measure the aggregated spillover effects.To complete the Spatial DCC-GARCH model, we then estimate the conditional correlation matrix following the two-step procedure described in Engle and Sheppard (2001), see Section 2.1.2.
The condition for the weak stationarity of the spatial GARCH model follows from the corresponding stationarity condition for E-CCC models, derived by Jeantheau (1998) and Conrad and Karanasos (2010) for E-CCC models.The positivity conditions on all GARCH coefficients are not necessary for the positivity of variance and in many empirical cases, these may be too restrictive, ruling out possible negative volatility feedback.One author Conrad and Karanasos (2010) studied the E-CCC models and stated necessary and sufficient conditions (in terms of the process parameters) for the positivity of variance; these conditions are summarized in Theorem 1 of their paper.It can be seen easily that our spatial DCC-GARCH(1,1) model is equivalent to the E-DCC model of order one, with the particular form of the parameter matrices A = A 1 + A 2 W and B = B 1 + B 2 W.So for the conditional variances to be positive, the conditions (C1)-(C3) of Theorem 1 of Conrad and Karanasos (2010) must apply.The proposed spatial DCC-GARCH(1,1) model is weakly stationary if the modulus of the largest eigenvalue of the matrix A 1 + B 1 + (A 2 + B 2 )W is less than 1.In that case, the unconditional variances are given by A 0 (I − (A 1 + B 1 ) − (A 2 + B 2 )W) −1 .More specifically the unconditional variance of the ith bank is given by and it is positive for a 0,i > 0, a 2,i + b 2,i > 0, (a 1,i + b 1,i ) < 1.

Estimation of the Multivariate Spatial GARCH(1,1) Model
We follow a two-step procedure for the DCC-GARCH estimation, as described in Engle and Sheppard (2001) and Engle (2002).The first step is devoted to the estimation of ( 16) where the exogenous variable Y t,i is not observable since it is a function of the conditional variance of the other assets.Hence, following Borovkova and Lopuhaa (2012), we start by estimating the standard univariate GARCH(1,1) models without the external regressors to obtain the initial parameters a 0 0,i , a 0 1,i , b 0 1,i and the estimated variances h 0 1,i , . . ., h 0 T,i .Then we use an iterative procedure in which we alternate the following two steps:

•
Compute the exogenous variables (Y t−1,i ) given the weights w ij and the initially estimated variances h 0 1,i , . . ., h 0 T,i ; • Estimate the complete set of parameters a 1 0,i , a 1 1,i , b 1 1,i , a 1 2,i , b 1 2,i and the new estimated variances h 1 1,i , . . ., h 1 T,i according to Equation ( 16).We iterate this procedure until the percentage variation of the estimate is less than a small threshold.For more details please refer to Borovkova and Lopuhaa (2012).
In the second step, as in Engle (2002), we maximize the quasi log-likelihood that, when the standardized error t follows a multivariate Gaussian distribution is where θ = (θ 1 , θ 2 ) is the set of parameters of the multivariate distribution, with subsets θ 1 = (A 0 , A 1 , A 2 , B 1 , B 2 ) being the spatial GARCH parameters estimated in the first step, and θ 2 = (γ, δ) the parameters of the time-varying conditional correlation that are estimated in the second step.Excluding D t and other additive and multiplicative constants, we maximize the following function: The quasi-log-likelihood function under the Student's t distribution is where ν is the degrees of freedom, θ 2 = (γ, δ, ν) is the set of parameters estimated in the second step, and Γ(•) is the Gamma function.The estimation of the model is implemented in R using the packages rugarch and rmgarch for the estimation of univariate GARCH models in the first step, and the DCC-GARCH model in the second step, respectively.
Concerning the complexity of estimation, we see that the spatial models add two parameters (a i,2 ,b i,2 ) for each asset.Hence the number of additional parameters scales linearly with the size of the dataset considered.We also see that Student's t model has one extra parameter compared to the Gaussian model, and that in the limit for ν → ∞, the former converges to the latter.Moreover, the spatial model nests the non-spatial DCC-GARCH models, where the coefficients of the spatial components are restricted to zero.The Spatial DCC-GARCH model can therefore be considered parsimonious in terms of the number of parameters, especially compared to VEC GARCH or the BEKK model.One drawback of the proposed model is that it requires the exogenous identification of a spatial matrix.

Spatial Weight Matrix
To estimate the spatial DCC-GARCH described in Section 2.1.1,we need to specify the weight matrix W which incorporates the spatial structure defined a priori.The most intuitive way to compute the weights is to consider the geographical distance between the issuers' market cities.However, according to Borovkova and Lopuhaa (2012), the obtained weights are not economically meaningful, and as an alternative, they consider a different set of information and compute distance in terms of GDP and market capitalization.In our work, we investigate whether the banks' similarity of the structure of credit exposure provides some benefit in catching risk spillover effects.Hence we propose to consider the cosine similarity between exogenous information relative to the credit exposure of each bank derived from the EU-wide stress test under the European Banking Authority (EBA).The higher the cosine similarity the stronger the closeness of banks's credit exposure.Suppose two attribute vectors of length L, U i,L = (u i,1 , u i,2 , . . . ,u i,L ) and U j,L = u j,1 , u j,2 , . . ., u j,L which describe the credit exposure information of bank i and j with i, j = 1, . . ., N. We define the cosine similarity as follows: We set C ii = 0 ∀ i and we normalize the rows of C by dividing each element by the sum of the row.Doing so, we obtain the matrix W that is the spatial weight matrix used in Equation ( 15).

Financial Application: CoVaR
Financial institutions use VaR to measure the standalone risk.However, the measurement of individual risk is not able to explain the linkages between other financial institutions and the financial system.Systemic risk is the possibility that an event at the institutional level could trigger severe instability or collapse of an entire industry or economy.The work Adrian and Brunnermeier (2014) introduces CoVaR to help regulators to measure risk spillovers.
The Value at Risk (VaR) at level q ∈ (0, 1) of a random variable r with cumulative distribution function F r (.) is defined as where 100(1 − q)% denotes the confidence level of the VaR. 2 Restricting our analysis to continuous probability distribution functions, VaR can be implicitly defined as the q-quantile of the probability distribution function , is implicitly defined by the q-quantile for a continuous probability distribution function of the financial system S conditional on some event related to C(r i ), where r i is the return of institution i such that Pr r S ≤ −CoVaR The CoVaR can capture the contribution of systemic risk by conditioning the VaR to a stressed situation.It captures the spillover of risk between a particular institution and the financial system, and it is commonly used to assess the systemic risk of a bank in a financial system.Inspired by this idea, we concentrate our attention on a CoVaR pairwise analysis between institutions in order to quantify the spillover between couples of banks. 3 The conditioning event C(r i ) in the original paper by Adrian and Brunnermeier ( 2014) is defined as the return of the conditioning asset i being equal to its negative VaR, that is C(r i ) := (r i |r i = −VaR(r i )).In this work, we follow the alternative approach of Girardi and Ergün (2013) that considers as a conditioning event the return r i being smaller or equal than the following quantity: C(r i ) := (r i |r i ≤ −VaR(r i )).This formulation allows us to consider more severe distress events and improves the consistency of the measure with respect to the dependence parameter, allowing for backtesting.Following Girardi and Ergün (2013), the redefined Let f r j , r i be the bivariate probability distribution function of future returns, estimated using the DCC-GARCH model with either Gaussian or Student's t innovations, CoVaR j|i q is implicitly defined as the quantity that solves We compute the integral (23) on a grid of 100 values for CoVaR j|i q to find the approximated solution under the different distributional assumptions. 4

CoVaR Based on Filtered Historical Simulations (FHS)
In order to compare our result with a model-free approach, we consider the Filtered Historical Simulations (FHS).
FHS is a well-known tool for multivariate forecasting and simulation of time series that avoids the need for distributional assumptions for the returns' joint dynamic, relying instead on past realizations.The main novelty of this approach compared to historical simulation is to rescale the innovation by the volatility that prevails on a specific day, allowing therefore to reflect the current market conditions (Barone-Adesi et al. 2002;Giannopoulos and Tunaru 2005;Gurrola-Perez and Murphy 2015).To provide a distribution-free benchmark model for the analysis, we compute the VaR and CoVaR via FHS.Consider a time window of length T and let r t be the series of historical returns with t ∈ [1, T].The volatility weighted returns series can be computed as follows: z t = r t × σT+1 / σt , where σt is the volatility estimated with an Exponentially Weighted Moving Average procedure (EWMA) with decay factor λ = 0.9 at time t and σT+1 is the one-day-ahead estimate of volatility at the end of the estimation period.In practice, implementing FHS for the estimation of VaR and CoVaR requires the following steps: • compute the residual (or devol) time series, dividing the returns by EWMA estimated volatility σt .This allows us to sample from approximately serially independent and identically distributed data; • compute the estimated empirical distribution of rT+1 (revol), multiplying the devol time series z t by the latest estimate of volatility σT+1 and assigning to each of the possible outcome a weight 1/T, • estimate VaR i q,T+1 and CoVaR j|i q,T+1 by computing the empirical quantile of ri,T+1 and rj,T+1 | ri,T+1 < −VaR i q,T+1 , respectively.The FHS approach has the advantage of being non-parametric, although it has the drawback of requiring a large number of observations to accurately estimate risk, especially for the CoVaR.For this reason, it is not suitable for small values of q.For instance, with q = 0.01 the expected number of exceedances of the CoVaR for an estimation window of 10,000 daily observations (approx 40 years) is 1, while for q = 0.05 it is 25.

Backtesting VaR and CoVaR
In order to test the goodness of our VaR and CoVaR estimates we estimate a time series of length τ of one-day-ahead estimates, each computed on an estimation window of T = 1000 daily observations.We consider tests based on the number of violations and specifically unconditional and conditional coverage tests (Christoffersen and Pelletier 2004;Kupiec 1995), as well as tests based on asymmetric loss functions for the VaR and CoVaR (Caporin 2008).The model that provides estimates of VaR and CoVaR with the correct number and distribution of exceedances and/or lower loss function values will be considered the more accurate.

Tests Based on the Number of Violations
In order to determine the accuracy of the proposed model, we consider two tests based on the number of violations.
Denote by • r i t the ex-post realized returns of institution i with t = 1, . . ., τ; • VaR i q,t the ex-ante Value-at-Risk forecasts at t − 1 for time t, where q is the expected coverage; • I i t a sequence of violation for a given interval of the Value-at-Risk forecast: The first test is the Kupiec test or unconditional coverage (UC) test (Kupiec 1995).The null hypothesis that the observed failure rate p is equal to the failure rate, suggested by the confidence level of VaR, q, is tested.Thus, the null hypothesis assumes that the observed violation rate is equal to the expected violation rate.If the null hypothesis is rejected, the model is considered inaccurate at the 95% confidence level.
The conditional coverage (CC) test proposed by Christoffersen and Pelletier (2004) indicates that the number of violations must be independently distributed along the testing period where the dependence can be described as a first-order Markov sequence with a transition probability matrix given by where π 01 is the probability that, conditional on today being a non-violation, the next period is a violation, and π 11 is the probability that, conditional on today being a violation, the next period is a violation.The hypothesis to test for the conditional coverage property is H 0 : π 01 = π 11 which assesses the independence of failures on consecutive time periods.Girardi and Ergün (2013) proposed the backtesting of CoVaR j|i q,t via a straightforward application of the standard Kupiec and Christoffersen tests considering the violations (i.e., r j t ≤ −CoVaR j|i q,t ) for those time periods in which r i t ≤ −VaR i q,t .Having that in mind we compute a second hit sequence, I j|i t , on the sub-sample in which r i t ≤ −VaR i qt as follows: where the number of observations of the second hit sequence is equal to the number of violations of the first hit sequence.Hence for the tests on CoVaR, the sequence of violation I j|i t can be used instead of I i t .

Backtesting Based on Loss Functions
The backtesting based on the confidence level of VaR estimates shows the accuracy of an individual model; however, the comparison between the different models can be limited.To overcome the drawback, Lopez (1999) proposed backtesting based on a loss function.The method focuses on the magnitude of the failure when the violation occurs.Thus, the VaR estimates under the loss function can provide the model's performance as a numerical score.The value of the loss function at time t can be given as where g(•) and h(•) are the loss functions applied to exceedances and values within the VaR, respectively.Finally L i = ∑ T t=1 l i t is defined as the total loss.The best model can be identified by the lowest total loss.Other works by Abad et al. (2014), Caporin (2008), and Cesarone and Colucci (2016) show several alternative specifications for the loss functions g(•) and h(•), defined from the regulator and investor's point of view.In the regulator's view, we consider the size of the loss only if the violation occurs: On the contrary, the investor is interested in both sides, as an overestimation of VaR may trigger limitations from the risk management, or lead to higher capital requirements imposed by the regulator.In particular we consider the functions g(r t,i , VaR i q,t ) = |r t + VaR i q,t | and h(r t,i , VaR i q,t ) = q 1 − q g(r t,i , VaR i q,t ).
We underline that the resulting loss function l i t is strictly related to the Koenker loss function used for the estimation of quantile regression, defined as where (•) + = min(X, 0) and (•) − = min(−X, 0).In case of independent and identically distributed returns the minimization arg min ξ l(X, ξ, α) is the value at risk.For further details we refer to Koenker and Bassett (1978), Rockafellar and Uryasev (2013), and Giacometti et al. (2021).
We extend the backtesting procedure to the case of the CoVaR as before, estimating the measure l j|i t on a sub-sample in which r i t ≤ −VaR i,t q .

Data
We consider ten years of weekly data from seven representative banks in Italy, France, Germany, the United Kingdom, the Netherlands, Spain, and Belgium.The data span from 20 September 2010 to 18 September 2020, including 2566 daily equity log-returns.The data are downloaded from Refinitiv Eikon.We perform a rolling analysis with an estimation window of T = 1000 daily observations (approximately 4 years), forecasting one-day-ahead VaR and CoVaR, for a total of τ = 1566 out-of-sample daily observations.We use the same windows for both the DCC-GARCH models and the FHS estimation.
Table 1 reports the descriptive statistics and tests the output of the log returns for the out-of-sample period.We see that all the banks in the sample with the exception of KBC Group had negative average returns.The series have typically negative skewness and excess kurtosis, as expected from equity time series.The results of the Engle ARCH test (Engle 1982) indicate that the null hypothesis of homoscedastic returns is rejected, suggesting the need for GARCH models.The autocorrelograms and partial autocorrelograms of the returns, not reported for brevity, do not highlight relevant serial correlation structure, while the autocorrelogram of squared residuals (also omitted for brevity) show significant and persistent autocorrelations, confirming the heteroscedasticity of the data.Next, we consider the correlation between the banks.Figure 1 shows that correlations are positive and high.Figure 2 studies the evolution of correlations, computed using 6-month rolling windows.The right panel represents the dynamics of the 21 pair-wise correlations over time, while the left panel shows the Frobenius norm of the correlation matrices to provide a synthetic representation.We see that the Frobenius norm changes over time, suggesting that a CCC-GARCH model is not appropriate for the dataset.On the contrary, the time-varying correlation matrix is consistent with the assumptions of a DCC-GARCH, and the high variability in the individual correlations leaves space for spatial models, which could better characterize the multivariate stochastic process.

Spatial Weight Data
For the construction of the spatial matrix, we analyze the data from the EU-wide stress test under the European Banking Authority (EBA).The EBA stress test aims to evaluate financial institutions' resilience to adverse market conditions.It also provides the overall assessment of systematic risk in the European banking system.In the EU-wide stress test analysis report, we consider the base scenarios for each bank and the relative credit exposure information: exposure values, risk exposure amounts, stock of provision, and leverage ratio under the internal ratings-based (IRB) approach or Standardized approach (STA) referred to credit exposure in specific asset classes, 5 as presented on the EBA's website (EBA 2021).We compute for each bank the vector of percentage exposure in each class with respect to the total exposure and the similarity between couples of vectors for the different banks, as in Equation ( 21).This indicator provides a broad view of the similarity between the banks' credit structures and exposures.
We then rescale the values such that each row sums to 1, as shown in Table 2.The spatial matrix is based on the EU-wide stress test of 2018 and is kept fixed for the entire analysis.To ascertain the matrix weight's consistency over time, we perform the inequality test for couples of matrices by Jennrich (1970) on matrices computed in different years.In particular, we compare the normalized cosine similarity matrix weight from the EUwide stress test of 2014 vs. 2016, 2016 vs. 2018, and 2018 vs. 2014.We do not reject the null hypotheses at a 1% significance level, suggesting that the spatial components of the EU-wide stress test do not change significantly over time.

Empirical Results
We estimate the DCC-GARCH(1,1) models using the equity log returns, considering four specifications that differ in terms of distribution and inclusion of the spatial component: Gaussian DCC (GaussDCC), spatial Gaussian DCC (GaussSpDCC), Student's t DCC (tDCC), and spatial Student's t (tSpDCC). 6The procedure is numerically stable, and only for a small percentage of the estimation window, do the DCC-GARCH models fail to converge.In such cases, we carry over the result from the previous estimation window.
Table 3 reports the Akaike, Bayesian, Shibata, and Hannan-Quinn information criteria, averaged across the 1566 estimation windows.Information criteria allow us to assess the quality of the model in relation to the data, controlling for both the quality of the fit and the number of parameters.Therefore, we use them to assess whether the inclusion of a distribution with more parameters (Student's t, compared to the Gaussian) and the use of the spatial component actually improve the quality of the fit or, on the contrary, the added complexity of the model affects negatively the estimation.We see that the models based on the Student's t distribution have lower information measures according to all the measures considered (denoting a better model for the data), and that the introduction of the spatial component improves the performance of the model.To further illustrate the characteristics of the model, Table 4 shows the average value of the coefficients computed on the 1566 rolling windows (for brevity we report the averages of the parameters related to each bank).We notice that a relatively large part of the spatial parameters (a 2,i ,b 2,i ) are statistically significant at the 95% significance level, suggesting that the spatial components have some explanatory power, motivating their introduction in the model.We also see that the DCC parameters γ and δ are almost always statistically significant (confirming the presence of time-varying correlations).The ν parameter is also significant and has values close to 6, suggesting that the innovations have fat-tails and that the Student's t model is more suitable than the Gaussian one (that is the limit of the Student's t model with ν → ∞) Table A1 in the Appendix A reports the coefficients of the tSpDCC model for the first estimation window as an example.

VaR and CoVaR
In this section, we use the spatial DCC-GARCH model to compute pair-wise CoVaRs and study risk spillover in the European banking system.We report in Table 5 the estimates of the average VaR i 5% (diagonal elements) and CoVaR j|i 5% (off-diagonal elements) for the four DCC-GARCH models and for the Filtered Historical Simulations (FHS) estimates.The one-day-ahead forecast of VaR i 5% is computed using the conditional variance estimate and the parametric distribution of the model.The corresponding CoVaR 5% is computed numerically, according to (23) using the time-varying covariance matrices.We observe that, as expected, the CoVaR is always greater in absolute values than the VaR figures on the diagonal.Furthermore, we see that the estimates of the Value at Risk are similar across the five considered models.On the contrary, the estimates of CoVaR are significantly larger for the DCC-GARCH models based on the Student's t distribution, suggesting that the Gaussian model may potentially underestimate the risk of joint distress and risk spillover, as expected by the stylized fact of fat tails and high tail correlations in financial time series (see, e.g., Cont 2001).Finally, we see that the estimates of CoVaR of the tSpDCC model are slightly smaller than the tDCC model.The FHS estimates yield similar results to the other models in terms of VaR, and lie in the middle between the Gaussian and Student's t models in terms of CoVaR. Figure 3 shows the out-of-sample equity log returns and the estimate of the Value at Risk.We see that the dynamics are similar for all the models, and that the two models with a spatial component share some similar dynamics in specific time periods.Figure 4 shows the VaR j 5% of each bank j together with the average CoVaR j|i 5% for i = j (for brevity we report only the results estimated using the tSpDCC model).We see that for all the banks, the CoVaR is always higher than the VaR.Focusing on the dynamics, we see two main spikes in the series that affected all banks: one in mid-2016 (corresponding to the Brexit Referendum), and one when the COVID-19 crisis started in March 2020.The former shock was short-lived, and risk measures returned to normal levels quickly, while the Covid crisis had more long-lasting effects, with CoVaR decreasing slowly in the following months (although with differences across banks, for instance, Intesa San Paolo risk levels returned to a normal level quicker than ING Group).In other periods, the dynamics of CoVaR are diversified across banks, with some institutions (in particular Intesa San Paolo) characterized by several spikes (likely related to idiosyncratic or regional shocks) while other institutions such as Credit Agricole or Barclays characterized by a more stable risk profile.The last panel of Figure 4 shows the average ratio between CoVaR j|i 5% and VaR j 5% for each couple of institutions, and the rescaled average log-returns of the seven banks included as reference.We see that the shocks associated with Brexit and Covid had the effect of increasing the ratio, thus increasing the risk spillover in the system.The effect is persistent, as the ratio remained higher for the periods after the shocks, despite the levels of risks reduced more quickly after the event, suggesting the presence of long-lasting spillover effects.

Backtesting Results
We now present two backtesting analyses in order to assess the quality of the estimation of VaR and CoVaR based on different models.First, we report the results for the backtesting of VaR i 0.05 with the Kupiec (unconditional) and Christoffersen (conditional) coverage tests.Table 6 provides the p-value of the unconditional coverage (UC) and conditional coverage (CC) tests.We see that, regarding the UC test, in most cases, the null hypothesis of correct exceedances is not rejected at the 95% confidence level, meaning that the models identify correctly the expected number of exceedances.The spatial component, according to this statistic does not provide benefits, leading typically to lower p-values and higher rejection rates.Concerning the CC test, the null hypothesis of correct and independent exceedances is rejected in a higher number of cases, indicating the presence of some residual clustering of the exceedances.Summing up, according to the CC tests, the introduction of the spatial component improves the estimation of Value at Risk, leading to higher p-values and lower rejection rates.The FHS estimation has good performance in terms of the UC test, while it is the worst one considered in terms of the CC test (the null hypothesis is rejected for all banks except for Banco de Sabadell).Table 7 reports a summary of the UC and CC tests for the CoVaR.Panels A and B show the average number of exceedances, the average p-value across all the combinations of assets, and the percentage of p-values that are smaller than 5% (i.e., the null hypothesis is rejected at the 95% confidence level) for the UC test and CC test, respectively. 7We see that for the GaussDCC and GaussSpDCC based on the Gaussian distribution, the null hypotheses for both the UC and CC tests are rejected for all the bilateral CoVaRs.On the contrary, for the Student's t models the null hypotheses are never rejected for both the UC and CC tests.Finally, we observe that the spatial model tSpDCC performances are aligned with those of tDCC.The results of the FHS are slightly better than the Gaussian DCC-GARCH models but are worse than the Student's t DCC-GARCH models.The results, although they do not highlight relevant differences between the spatial and non-spatial models, do confirm that the models based on the Gaussian distribution are not suitable to properly measure spillover risk, while Student's t models have much better results when compared to FHS estimation of CoVaR. .The null hypotheses are "Correct Exceedances" (unconditional test) and "correct & independent exceedances" (conditional test).Averages and percentages are computed across all the couples of assets.For q = 5% the expected number of exceedances is 3.9.Finally, we assess the quality of the estimation of VaR and CoVaR using the loss function methodology described in Section 2.2.4.Table 8 reports the value of the loss function for the investor and the regulator (lower values are better).We see that, consistently with the results for CC and UC tests, the introduction of the spatial component has a positive effect for both the models based on Gaussian and Student's t distribution.Contrary to expectations, the Student's t models do not outperform the Gaussian models according to this metric, showing similar or slightly worse performance.The similarity between Gaussian and Student's t distribution may be due to the fact that we are considering a low confidence level (95% for q = 0.05): for a higher confidence level (e.g., 99%, q = 0.01) the shape of the tails may matter more, and the Student's t model may perform better.In this analysis, we do not consider higher confidence levels as the backtesting of the CoVaR would become impossible.Looking at Table 9, we see the average value of the loss functions for the CoVaR estimations.The results highlight once again a beneficial effect of the spatial component, that leads to improvements regardless of the distribution and the loss function used.We also see that the average loss is clearly smaller for the Student's t models, highlighting their ability to better estimate the spillover risk compared to the Gaussian models.The results for the FHS show that such a model is not as accurate as tSpDCC and tDCC, although it performs better than Gaussian DCC-GARCH models.

Panel
Overall, the results suggest that the introduction of a spatial component improves the performance of the DCC-GARCH model, both in terms of the information criteria, and in terms of the out-of-sample estimation of CoVaR, as confirmed by the backtesting.The results also confirm the better fit of the Student's t models compared to the Gaussian models.

Discussion and Conclusions
As the spillover effects of risk become a problem in interconnected banking systems, this study introduces a spatial DCC-GARCH(1,1) to provide a more accurate measurement of joint tail risk in a parsimonious way.The model aims to improve the standard DCC-GARCH model (Engle 2002), without introducing the estimation complexity typical of VEC GARCH and BEKK model (Bollerslev et al. 1988;Engle and Kroner 1995;Ling and McAleer 2003).After discussing the estimation of the model, we perform an empirical analysis on the equity log returns of seven large European banks, using a matrix that reflects the similarity in credit structure and exposures.We compare four multivariate GARCH specifications (with and without spatial components, and with two alternative distributions for the innovations).The comparison of information criteria suggests that the proposed spatial model with Student's t innovation provides a better fit compared to the alternatives.
A common approach to measuring the spillover risk is the usage of pairwise CoVaR Adrian and Brunnermeier (2014), which measures the tail risk of an institution conditional to the distress of another institution.We estimate pairwise CoVaRs using our spatial DCC-GARCH(1,1) model.Compared to other GARCH-based estimation procedures for CoVaR (see, e.g., Girardi and Ergün 2013), our framework allows us to consider a time-varying correlation matrix thanks to the DCC component, and a network dimension thanks to the spatial component.Indeed, as highlighted in the literature, the network component of risk is more and more relevant (see, e.g., Billio et al. 2012;Diebold and Yılmaz 2014), and spatial GARCH models allow us to include it while maintaining the estimation feasible.
To test the reliability of the VaR and pairwise CoVaR, we first examine their accuracy via the UC and CC tests considering different GARCH specifications and a non-parametric FHS approach.The results show that the Student's t spatial DCC-GARCH(1,1) model (tSpDCC) provides the lowest rejection rate for CoVaR 5% compared to other models.Second, we investigate the models via backtesting based on loss functions.The analysis confirms that the Student's t spatial DCC GARCH(1,1) model outperforms the other DCC-GARCH specifications, as well as the filtered historical simulations model in terms of the estimation of CoVaR 5% .Overall, from a methodological perspective, we conclude that the multivariate GARCH model with the Student's t and the spatial component obtained thanks to the proposed similarity matrix can improve the assessment of credit risk profiles and the credit risk market's spillover.
Concerning the economic analysis and interpretation of the empirical results, the spillover analysis shows that the dynamics of CoVaR were diversified across European banks, and that in the out-of-sample period (2014-2020) there were two main shocks common to all the institutions: the Brexit Referendum (mid-2016) and the COVID-19 Crisis (started in the first half of 2020).Both periods were associated with spikes in VaR and CoVaR, and persistent increases in the ratio of CoVaR over VaR (denoting therefore a long-lasting increase in the interconnectedness and risk spillovers).The persistent increase in CoVaRs after Brexit is consistent with Li (2020), which studies the behavior of European stock markets in a multivariate time-varying setting, finding that market co-volatility continues to be substantial and persists after Brexit despite the fact that the market adjusted quickly to the shock.Similarly, the increase of CoVaR during the Covid Crisis confirms results from the literature that show how spillover and interconnectedness increased in the first part of the Covid period.In particular, Aslam et al. (2021) studied twelve European markets using the methodology from Diebold and Yılmaz (2014) on high-frequency data and found more stable spillovers in the Covid period compared to the previous, and Foglia et al. (2022) show an increase of volatility connectedness during the Covid period across 30 major Eurozone banks.
Finally, we point out that the proposed framework not only has relevant applications for financial regulators, but it is also relevant for the asset management industry.Indeed, the proposed model allows us to estimate the joint distribution of future returns, providing asset managers with reliable inputs for optimal portfolio strategies (Meucci 2005), and risk managers with data useful to measure the risk of investment funds and to conduct stress-testing analyses based on hypothetical scenarios (see, e.g., Alexander and Sheedy 2008;Koliai 2016).
Future works may further extend the current work by testing alternative spatial matrices and applying the model to other markets and larger datasets.

Figure 1 .
Figure 1.Correlations, marginal distributions, and bivariate scatterplots with LOESS local regression lines of equity daily log-returns.

Figure 2 .
Figure 2. Frobenius norm of rolling correlation matrix (left panel), and 21 pair-wise rolling correlations (right panel).Rolling correlations are computed using a 6-month window of daily observations.

Figure 3 .
Figure 3. Out-of-sample equity log-returns of the seven banks considered in the study and VaR estimates.

Figure 4 .
Figure 4. Out-of-sample equity log-returns, VaR estimated using the tSpDCC model (blue lines), and the average of the CoVaRs for each bank (red lines).The last panel reports the average ratio between the CoVaR and the VaR (moving average over 10 days), and the (rescaled) average log-return of the seven banks for reference.

Table 3 .
Akaike, Bayesian, Shibata, and Hannan-Quinn information criteria.The table reports the average across estimation windows.The smallest value for each measure is highlighted in bold.

Table 4 .
The table reports the average values of the coefficients of the different specifications across the estimation windows of the multivariate GARCH.The percentages of estimation windows in which the coefficients are significantly different from 0 at the 95% confidence level are reported in italics.For brevity we report the values of a 0,i , a 1,i , a 2,i , b 1,i , b 2,i averages across the seven banks considered in the analysis.

Table 5 .
Estimated VaR i 5% (diagonal elements) and CoVaR j|i 5% (off-diagonal elements) for the spatial and non-spatial DCC GARCH models with Student's t and Gaussian innovations, and for FHS.The reported values are the average of the out-of-sample estimates computed across 1566 daily rolling windows.

Table 6 .
p-Values of the unconditional (Panel A) and conditional (Panel B) for VaR i 5% .The null hypotheses are "Correct Exceedances" (unconditional test) and "correct & independent exceedances" (conditional test).p-values lower than 5% are highlighted in bold.

Table 7 .
Average exceedances, p-values, and percentage of p-values smaller than 0.05 of the unconditional (Panel A) and conditional (Panel B) for CoVaR j|i 5%

Table 8 .
Backtesting based on loss function for VaR i 5% .Panel A reports the investor's point of view and Panel B reports the regulator's point of view (see Section 2.2.4).For each line the best value is highlighted in bold.

Table 9 .
Backtesting based on loss function for CoVaR j|i 5% .The null hypotheses are "Correct Exceedances" (unconditional test) and "correct & independent exceedances" (conditional test).Panel A reports the investor's point of view and Panel B reports the regulator's point of view (see Section 2.2.4).For each line, the best value is highlighted in bold.
Parameters of the Student-t Spatial DCC-GARCH(1,1) model estimated on the first rolling window (September 2010-September 2014).p-value in brackets, values significant at 95% confidence level in bold.