Valentine Wanjiku Mwai, Cyprian Ondieki Omari, Simon Maina Mundia
Department of Statistics and Actuarial Science, Dedan Kimathi University of Technology, Nyeri, Kenya.
DOI: 10.4236/jmf.2025.154031   PDF    HTML   XML   128 Downloads   1,667 Views   Citations

Abstract

The global financial landscape is increasingly becoming interconnected, with financial markets exhibiting complex interdependencies. This increases the possibility of market risk spreading from one market to another, as market shocks often propagate across asset classes especially during periods of economic uncertainty. Failure to adequately capture the characteristics of univariate return series and the dependence structure between them, may lead to significant underestimation of the market risk forecasts. The standard multivariate Generalized Autoregressive Conditional Heteroskedasticity models assume that financial data follow a normal distribution, an assumption that fails to capture the heavy tails, skewness, and non-linear dependencies commonly observed in asset returns. Thus, this study models the dependency structures among a portfolio of financial asset classes and forecasts the Value-at-Risk and Expected Shortfall using multivariate Generalized Autoregressive Conditional Heteroskedasticity Vine Copula approach. The multivariate GARCH model captures the dynamic volatilities and conditional correlations among assets, then vine copulas are used to model the remaining non-linear and tail dependence relationships between the standardized residuals. The empirical results indicated that the financial return series exhibit complex dependence patterns that vary across asset classes and evolve over time, reflecting the diverse behaviors of financial markets under varying economic conditions. Among the models considered, the constant conditional correlation stationary vine copula demonstrates superior performance in dependence modelling. Backtesting results for one-day-ahead Value-at-Risk and Expected Shortfall indicate that constant conditional correlation vine copula models significantly outperformed the constant and dynamic conditional correlation models with normal and Student-t innovations over a one-day horizon. In contrast, dynamic conditional correlation vine copula models generally exhibit poor predictive accuracy and fail to meet key backtesting tests. Overall, the empirical findings of this study indicated that the constant conditional correlation regular vine copula offers the most reliable and precise framework for modelling and forecasting portfolio market risk.

Keywords

Portfolio Market Risk, Dependency Modelling, Backtesting, Multivariate GARCH, Constant Conditional Correlation, Dynamic Conditional Correlation, Regular Vine, D-Vine, C-Vine, S-Vine, Value at Risk, Expected Shortfall

Share and Cite:

1. Introduction

The financial markets are pivotal to the global economy, serving as platforms where various assets are traded, including stocks, bonds, currencies, commodities and recently cryptocurrencies. These markets facilitate the efficient allocation of capital, provide liquidity, play an important role in the processes of determining prices, and allow the transfer and redistribution of risk, thereby influencing economic growth and the stability of global financial systems. Financial markets operate within a framework of economic, political, and social influences, such as recessions, inflation, natural disasters, market sentiments, geopolitical events and changes in interest rate, which can introduce substantial uncertainty and risk, thus significantly influencing price fluctuations in the financial markets.

Investing in these markets means accepting a certain amount of exposure to risk. Portfolio market risk is the possibility of experiencing a financial loss in a portfolio due to adverse price movements. Accurately measuring market risk is important in effective risk management and some common measures of market risk are Value at Risk and Expected Shortfall. VaR assesses the highest possible loss a portfolio may incur over a set period and at a determined confidence level, whereas ES indicates the anticipated loss that exceeds the VaR threshold, thus offering thorough understanding of tail risk. These measures form the foundation of risk management practices, ensuring portfolios are prepared for potential adverse scenarios.

Forecasting portfolio market risk accurately relies on analyzing financial time series data, characterized by several unique features that influence risk modelling and management strategies. The data often exhibits volatility clustering, in which minor price changes typically follow small changes and large variations in price are more likely to follow large price changes. The data also has fat tails than that expected from the Normal (Gaussian) distribution. If normal distribution is assumed in modelling financial returns, the number and magnitude of crashes and booms maybe underestimated. The leverage effect is also observed in financial markets, where losses have a greater influence on future volatilities than gains. Another characteristic is asymmetry, where extreme negatives returns are more frequent than extreme positive returns. Additionally, interdependency is a fundamental characteristic of financial market data, reflecting the intricate connections between different markets, asset classes, and economic sectors.

The GARCH model, introduced by [1], is widely used to model volatility clustering in financial time series. Extensions of the GARCH model such as the Threshold GARCH (TGARCH), Glosten-Jagannathan-Runkle GARCH (GJR-GARCH), and Exponential GARCH (EGARCH), further improve volatility estimation by capturing asymmetry and the leverage effect. However, these models are primarily univariate, focusing on a single asset or return series at a time. Modelling multivariate financial time series data extends beyond univariate models to capture the interdependencies and analyze time series jointly. Unlike univariate models that focus on individual asset returns, multivariate GARCH consider the dynamic interactions and correlations among several assets simultaneously, leading to improved forecasts. This allows the joint volatility and covariance structures to be represented more accurately, which is important in risk and portfolio management. There is a wide range of MGARCH models including CCC, BEKK, and DCC models. Recent studies have compared MGARCH variants in portfolio risk and pricing applications, highlighting the benefits of modelling dynamic correlations. [2] and [3] found that the CCC model outperformed DCC in Bitcoin option pricing and portfolio selection, respectively. [4] showed that DCC and GO-GARCH provided more reliable VaR estimates than CCC by capturing time-varying volatility. Similarly, [5] demonstrated that multivariate GARCH models consistently outperformed univariate models in forecasting portfolio risk with greater accuracy and efficiency. However, MGARCH models still assume a normal distribution, and linear correlations and fail to capture dependency structures that exist in financial assets such as tail dependencies. These linear relationships can be inadequate for accurately modelling the relationships in financial markets.

Copulas, introduced by [6], link multivariate distributions to their individual distribution functions. This allows for a more flexible representation of dependencies beyond linear correlation, by modelling non-linear and tail dependencies between random variables. The dependency structure between financial assets also fluctuates over time due to market conditions, economic events, and investor sentiment. Since the seminal work of [7] on using copulas for credit risk and default correlations, copula-based models have gained popularity in finance for dependency modeling, portfolio optimization, and risk management, with recent findings by [8] that bivariate copulas are superior to MGARCH models in modelling asymmetric dependencies and enhancing hedging effectiveness.

Numerous parametric bivariate copulas were developed to address different types of dependencies. Among these, the elliptical and Archimedean copulas are the most commonly used for modelling bivariate dependence. While bivariate copulas are effective in modelling dependencies between two assets, financial portfolios are often composed of multiple assets, necessitating the application of multivariate copulas in modelling the dependency structures across these assets simultaneously. The multivariate copulas such as Gaussian (Normal), Student-t and hierarchical Archimedean copulas are commonly used in this context. However, some limitations of the multivariate elliptical copulas are that has only one parameter, which controls both the shape and tail dependence, reducing its flexibility in capturing diverse dependency structures. Moreover, these copulas assume that all marginal distributions are identical, which can be restrictive. Additionally, the flexibility of these extensions is reduced as they require additional parameter restrictions.

The pair-copula constructions (PCC), introduced by [9], decomposes complex multivariate dependencies into simpler bivariate copula components addressing the limitations of multivariate copula models. This probabilistic approach to constructing multivariate distributions using pair-copulas was further expanded and organized systematically by [10] and [11] through the development of a graphical framework known as regular vine. This R-vine copula model is hierarchical and breaks down the joint density into a series of marginal densities and pair-copulas. This approach offers greater flexibility in modelling multivariate distributions since it permits the integration of pair-copulas from various families within a vine copula, accommodating any possible dependency structure. [12] introduced the drawable and canonical vine copulas into the finance and insurance literature. [13], further suggested that the choice between D-vine and C-vine depends on the context. For instance, C-vines are ideal when one variable strongly impacts all others, while D-vines focus on modelling variables with temporal order or some form of sequential relationship.

Recent studies have demonstrated the effectiveness of vine copulas in modelling complex financial dependencies and enhancing risk forecasts. [14] applied R-vine copulas to six currency exchange rates, showing their ability to capture interconnectedness during periods of market stress and identifying the C-vine specification as the best fit. [15] employed the R-vine framework to Turkish equities to model sectoral interdependencies surrounding the 2008 global financial crisis, underscoring the importance of adaptive portfolio decisions under varying market conditions. [16] combined R-vine copulas with GARCH to model banking stocks in the ISE100, reporting significant improvements in VaR and Expected Shortfall forecasts compared to GARCH models. [17] also demonstrated that regular vine copulas provided superior risk measurement when compared with DCC-GARCH and elliptical copulas on Istanbul Stock Exchange equities. [18] identified D-vine decompositions as providing the best fit for modelling dependency structures among major stock indices, whereas [19] found that R-vine outperformed both C-vine and D-vine in capturing dependencies in international equity markets. While standard vine copula models are limited to capturing cross-sectional dependence, [20] proposed the stationary S-vine copula, which extends the framework to jointly model both cross-sectional and serial (temporal) dependencies, thereby enhancing its applicability to time series financial data. [21] further demonstrated the usefulness of S-vines for ESG-based hedging strategies, showing that assets with extreme ESG scores exhibit stronger non-Gaussian dependencies.

The integration of copula functions with MGARCH models has also emerged as a powerful approach for jointly modelling the time-varying volatilities and complex dependency structures of financial asset returns. [22] proposed the Copula-based MGARCH (CMGARCH) framework combining several MGARCH specifications with Archimedean copulas to improve joint density forecasts. Recent applications have demonstrated the versatility of copula-MGARCH models in various financial contexts. [23] examined risk spillovers and dynamic dependency structures between energy and BRICS equities markets by integrating DCC-GARCH with bivariate copulas, finding significant interconnections and effective modelling of positive dependencies. [24] analyzed the volatility and dependence between Bitcoin and the five largest stablecoins, showing that copula-based models outperformed BEKK in capturing the stable interconnection of these markets. Similarly, [25] employed a Markov Switching vine copula-MGARCH model to forecast Value-at-Risk and Expected Shortfall for a portfolio of equities, bonds, and gold, achieving more accurate risk forecasts compared to standard MGARCH and copula-MGARCH models. Despite the documented advantages of copula-MGARCH models, much of the existing literature has focused on bivariate or low-dimensional settings, often overlooking higher-dimensional portfolios and time-varying dependence during periods of market stress.

In this study the Multivariate GARCH Vine copula approach is used to model the dependence structure of financial assets considered and forecast the one-day ahead VaR and ES of an equally weighted portfolio. This approach integrates the ARMA-EGARCH for modelling heteroscedasitity in return series, DCC and CCC MGARCH for modelling time-varying and constant conditional correlation, and then fitting the vine copula to model the dependence structure in between the twelve financial assets in the portfolio. The one-day ahead VaR and ES are forecasted via the Monte Carlo based simulation using the DCC-Vine and CCC-Vine copula approaches. Finally, back-testing techniques are used to assess the accuracy of DCC-Vine and CCC-Vine models in forecasting VaR and ES. The back-testing results are compared to DCC and CCC with student-t and normal innovations.

The rest of the paper is structured as follows: Section 2 presents the methodology employed including the Multivariate GARCH models, copula dependence measures, vine copulas, and VaR and ES forecasting and their back-testing procedures. Section 4 gives the empirical results. Finally, Section 5 concludes the paper and presents recommendations for future research.

2. Methodology

This section presents the methodology of the study in detail. Section 2.1 presents the Multivariate GARCH model specification for volatility and correlation estimation. Section 2.2 provides details of copula functions and Section 2.3 introduces the copula-based dependence measures. Section 2.4 introduces the vine copula function for dependency modelling. Section 2.5 presents the MGARCH-Vine copula model structure. Measures of portfolio market risk are presented in section 2.6. Finally, section 2.7 is about back-testing techniques.

2.1. Multivariate GARCH model specification

The Multivariate Generalized Autoregressive Conditional Heteroskedasticity model extends the univariate GARCH to multiple time series data. In a univariate setting, a GARCH(p,q) model assumes that the variance of the time series is influenced by past values of the squared returns (or innovations) and past variances. For a series of returns $\left\{r_t\right\}$ the GARCH model is defined as:

$\begin{array}{l}{r}_t=\mu +e_t\\ e_t=h_t^{1/2}{z}_t\\ h_t={\alpha }_0+{\displaystyle \sum }_{i=1}^q{\alpha }_i{e}_{t-i}^2+{\displaystyle \sum }_{j=1}^p{\beta }_j{h}_{t-j}\end{array}$ (1)

where $\mu$ represents the average of the asset’s log return, $e_t$ is the innovation at time $t$ , $h_t$ stands for the variance, $z_t$ is a white noise process with a mean of zero and variance of one, ${\alpha }_0>0$ , ${\alpha }_i\ge 0$ , ${\beta }_j\ge 0$ and ${\alpha }_i+{\beta }_j\le 1$ .

In order to generalize the univariate GARCH model to a multivariate framework, we consider a multivariate time series $r_t$ , which represents the vector of log-returns of $d$ assets. This vector $r_t$ has a dimension of $d\times 1$ and is assumed to be conditionally heteroskedastic, meaning its variance changes over time depending on past information. Let ${ℰ}_{t-1}$ represent the information set created by the observed series $r_t$ up until time $t-1$ . The standard MGARCH framework is given by:

$r_t={\mu }_t+{ϵ}_t$ (2)

where $r_t$ represents a $d\times 1$ vector of returns, ${\mu }_t$ denotes the mean vector, and ${ϵ}_t$ represents the error terms vector associated with time $t$ ;

${\epsilon }_t=H_t^{1/2}{\eta }_t$ (3)

where $H_t$ is an $d\times d$ conditional covariance matrix of $r_t$ , $H_t^{1/2}$ is any $d\times d$ matrix at time $t$ that scales the innovations of ${\epsilon }_t$ , obtainable via a Cholesky factorization of $H_t$ and ${\eta }_t$ represents a vector of iid errors with dimension $d\times 1$ , such that $\mathbb{E}\left[{\eta }_t\right]=0$ and $\mathbb{E}\left[{\eta }_t{{\eta }^{\prime }}_t\right]=I$ . There are two primary categories of MGARCH models: those that model $H_t$ directly and those that focus on conditional variances and correlations, decomposing the covariance matrix into correlations and standard deviations.

2.1.1. CCC-GARCH Model

[26], proposed CCC model, one of the simplest multivariate correlation models. Since the matrix of conditional covariances in this model doesn’t change over time, it can be written as:

$H_t=D_t{R}_t{D}_t$ (4)

where $R_t=\left[{\rho }_{ij}\right]$ is the $d\times d$ positive definite conditional correlation matrix, with ${\rho }_{ij}=1$ for $i=1,\cdots ,d$ and is time invariant such that $R_t=R$ and $D_t=\text{diag}\left(h_{1t}^{1/2},\cdots ,h_{dt}^{1/2}\right)$ is $d\times d$ matrix consisting of the standard deviations of the ${\epsilon }_t$ elements, with each $h_{it}$ being a univariate GARCH process. This suggests that the conditional covariance matrix’s off-diagonal components are provided by:

${\left[H_t\right]}_{ij}=h_i^{1/2}{h}_j^{1/2}{\rho }_{ij},\forall i\ne j$ (5)

where $h_{it}^{1/2}$ and $h_{jt}^{1/2}$ are the standard deviations of asset $i$ and $j$ at time $t$ . Consequently, the vector form of $h_t$ can be given as:

$h_t=a+{\displaystyle \sum }_{j=1}^q{B}_j{ϵ}_{t-j}^{\left(2\right)}+{\displaystyle \sum }_{j=1}^p{C}_j{h}_{t-j}$ (6)

where $a$ represents a $d\times 1$ vector comprising constant values, while the matrices $B_j$ and $C_j$ are diagonal $d\times d$ matrices, and ${ϵ}_{t-j}^{\left(2\right)}={ϵ}_{t-j}\odot {ϵ}_{t-j}$ . The formulation in Equation (6) enables each $h_{it}$ to be influenced by its historical variances and prior squared returns. For $H_t$ to maintain positive definiteness, $R$ needs to be positive definite, and the vectors $a$ , along with the matrices $B_j$ and $C_j$ , must also be positive.

The Extended CCC (ECCC) model, which was introduced by [27], relaxed the requirement that the matrices $B_j$ and $C_j$ are diagonal. The relaxation allows each conditional variance equation to incorporate the past variances and squared innovations from all series. For example, $h_{it}$ in the first order ECCC model is equivalent to

$h_{it}=c_i+a_{11}{\epsilon }_{1,t-1}^2+\cdots +a_{1d}{\epsilon }_{d,t-1}^2+b_{11}{\epsilon }_{1,t-1}^2+\cdots +b_{1d}{\epsilon }_{d,t-1}^2$ (7)

where $i=1,\cdots ,d$ and $a_{ij}$ and $b_{ij}$ are parameters that determine how past squared returns and past volatilities of asset $i$ , respectively, influence the current volatility of asset $i$ . An advantage of the ECCC over CCC model is that the autocorrelation structure for the squared returns is more flexible. Although the CCC and ECCC models’ parameter structures are relatively simple and computationally effective, the assumption that $R_t$ remains constant can be overly limiting. Therefore, it may be beneficial to generalize the CCC model by allowing $R_t$ to vary over time while preserving the previous decomposition.

2.1.2. DCC-GARCH Model

The Dynamic Conditional Correlation model, introduced by [28], extends the CCC model by permitting correlations to varying over time, which helps in capturing the dynamic co-volatility among assets. Using a proxy matrix $P_t$ , the model breaks down $H_t$ into time-varying conditional standard deviations $D_t$ and the correlation matrix $R_t$ , such that

$R_t=P_t^{\ast -1/2}{P}_t{P}_t^{\ast -1/2}$ (8)

$P_t=\left(1-a-b\right)S+a{\eta }_t{{\eta }^{\prime }}_t+b{P}_{t-1}$ (9)

where $P_t^{\ast }=\text{diag}\left(p_{1t},\cdots ,p_{dt}\right)$ , while $S$ denotes standardized errors’ unconditional covariance matrix, such that

$S=\text{Cov}\left({\eta }_t{{\eta }^{\prime }}_t\right)=\frac{1}{T}{\displaystyle {\sum }_{t=1}^T{\eta }_t{{\eta }^{\prime }}_t}$

and parameters $a$ and $b$ represent scalars, where $a+b\le 1$ , $a\ge 0$ and $b\ge 0$ . While the aspect of DCC having two parameters can be seen as an advantage, it might also pose a challenge when many assets are considered, since all $\frac{d\left(d-1\right)}{2}$ processes are correlation are confined to the same structure. [29], introduced a Quadratic Flexible DCC model which applies the BEKK structure to $R_t$ . Thus, $P_t$ is given as:

$P_t=A^{\prime }SA+B^{\prime }{\eta }_{t-1}{{\eta }^{\prime }}_{t-1}B+C{P}_{t-1}{C}^{\prime }$ (10)

However, the Quadratic Flexible DCC has $\frac{3d\left(d+1\right)}{2}$ parameters, which

becomes impractical for large systems due to the curse of dimensionality. To reduce the parameters, assets can be grouped by characteristics such as industry or sector and a block diagonal structure applied to the coefficient matrices according to the grouping.

The Asymmetric DCC GARCH model by [30] incorporates asymmetry in the correlation, since the standard DCC model only permits asymmetries in the variance and not the correlations. Thus, $P_t$ is defined as:

$P_t=\left(S-BS{B}^{\prime }-CS{C}^{\prime }-G{S}^{-}{G}^{\prime }\right)+B^{\prime }{\eta }_{t-1}{{\eta }^{\prime }}_{t-1}B+C{P}_{t-1}{C}^{\prime }+G^{\prime }{\eta }_{t-1}^{-}{\eta }_{t-1}^{-}\text{'}G$ (11)

where the parameter matrices $B$ , $C$ , and $G$ are of dimension $d\times d$ ; ${\eta }_t^{-}=1_{\left\{{\eta }_t<0\right\}}\odot {\eta }_t$ denotes the element-wise product of ${\eta }_t$ and $S^{-}=\text{Cov}\left({\eta }_t^{-}{\eta }_t^{-}\text{'}\right)$ .

2.1.3. Estimating Parameters of Multivariate GARCH Models

The parameters of the Multivariate DCC and CCC models are estimated in two step, where in step one the univariate conditional variance parameters (${\beta }_v$ ) are estimated while the correlation parameters (${\beta }_c$ ) are estimated in the step two, given ${\beta }_v$ . When ${ϵ}_t$ follows a multivariate normal distribution with mean vector zero and covariance matrix $H_t$ , the d-dimensional density function is:

$f\left({\epsilon }_t|H_t\right)=\frac{1}{{\left(2\pi \right)}^{d/2}{\left|H_t\right|}^{1/2}}\text{exp}\left(-\frac{1}{2}{{\epsilon }^{\prime }}_t{H}_t^{-1}{\epsilon }_t\right)$ (12)

The likelihood equation for $T$ observations is:

$L\left(\beta \right)={\displaystyle \prod }_{t=1}^T\frac{1}{{\left(2\pi \right)}^{d/2}{\left|H_t\right|}^{1/2}}\text{exp}\left(-\frac{1}{2}{{ϵ}^{\prime }}_t{H}_t^{-1}{ϵ}_t\right)$ (13)

Thus, the log-likelihood equation is given as

$l\left(\beta \right)=-\frac{dT}{2}\log 2\pi -\frac{1}{2}{\displaystyle \sum }_{t=1}^T\log \left|H_t\right|-\frac{1}{2}{\displaystyle \sum }_{t=1}^T{{\epsilon }^{\prime }}_t{H}_t^{-1}{\epsilon }_t$ (14)

where $\beta$ represents all parameters in the models such as elements of $A$ , $B$ and $C$ , $\left|H_t\right|$ and $H_t^{-1}$ represent the determinant and inverse of $H_t$ , respectively, $d$ denotes total assets, and $T$ represents the total observations. If $H_t=D_t{R}_t{D}_t$ , then the function of log-likelihood in Equation (14) can be expressed as

$\begin{array}{c}l\left(\beta \right)=-\frac{dT}{2}\log 2\pi -\frac{1}{2}{\displaystyle \sum }_{t=1}^T\log \left|D_t{R}_t{D}_t\right|-\frac{1}{2}{\displaystyle \sum }_{t=1}^T{{\epsilon }^{\prime }}_t{D}_t^{-1}{R}_t^{-1}{D}_t^{-1}{\epsilon }_t\\ =-\frac{1}{2}{\displaystyle \sum }_{t=1}^T\left(d\log 2\pi +2\log \left|D_t\right|+\log \left|R_t\right|+{{\epsilon }^{\prime }}_t{D}_t^{-1}{R}_t^{-1}{D}_t^{-1}{\epsilon }_t\right)\end{array}$ (15)

In step one, $R_t$ in Equation (15) is replaced by an identity matrix such that the $l\left(\beta \right)$ corresponds to a sum of $d$ univariate log-likelihood functions:

$\begin{array}{c}{l}_v\left(\beta \right)=-\frac{1}{2}{\displaystyle \sum }_{t=1}^T\left(d\log \left(2\pi \right)+2\log \left|D_t\right|+\text{log}\left|I_n\right|+{{\epsilon }^{\prime }}_t{D}_t^{-1}{I}_n^{-1}{D}_t^{-1}{\epsilon }_t\right)\\ =-\frac{1}{2}{\displaystyle \sum }_{t=1}^T\left(d\log \left(2\pi \right)+2\log \left|D_t\right|+{{\epsilon }^{\prime }}_t{D}_t^{-1}{I}_n{D}_t^{-1}{\epsilon }_t\right)\\ =-\frac{1}{2}{\displaystyle \sum }_{t=1}^T{\displaystyle \sum }_{i=1}^d\left(d\log \left(2\pi \right)+\log \left(h_{i,t}\right)+\frac{{\epsilon }_{i,t}^2}{h_{i,t}}\right)\end{array}$ (16)

In step two, the log-likelihood function given ${\theta }_v$ in step one is

$l_c\left(\beta \right)=-\frac{1}{2}{\displaystyle \sum }_{t=1}^T\left(\log \left|P_t\right|+{{\epsilon }^{\prime }}_t{P}_t^{-1}{\epsilon }_t-{{\epsilon }^{\prime }}_t{\epsilon }_t\right)$ (17)

Thus, ${\widehat{\beta }}_v=\arg \max \left\{l_v\left(\beta \right)\right\}$ and ${\widehat{\beta }}_c=\arg \max \left\{l_c\left(\beta \right)\right\}$ .

2.2. Copula Functions

A copula is a function that creates a joint distribution by coupling univariate marginal distributions. The concept of copula introduced by [6], separates the effect of the dependency between variables from the impact of each marginal variable.

Definition 3.1. The function $C\left(u_1,\cdots ,u_d\right)$ , which maps ${\left[0,1\right]}^d\to \left[0,1\right]$ , is a copula and satisfied the following properties:

1) $C\left(u_1,\cdots ,u_d\right)$ is grounded, meaning that for any $C\left(u_1,\cdots ,u_d\right)=0$ if any $u_i=0$ .

2) $C\left(u_1,\cdots ,u_d\right)=C\left(u_1,\cdots ,u_{i-1},u_{i+1},\cdots ,u_d\right)$ if any $u_i=1$ .

3) $C\left(u_1,\cdots ,u_d\right)$ is d-increasing for any two points $l$ and $m$ in ${\left[0,1\right]}^n$ , where $l\le m$ and the volume ($V_C\left(\left[l,m\right]\right)$ ) of the box formed by the vertices $l$ and $m$ is non-negative. For instance, when $d=2$ ,

$V_C\left(\left[l,m\right]\right)=C\left(l_2,m_2\right)-C\left(l_1,m_2\right)-C\left(l_2,m_1\right)+C\left(l_1,m_1\right)\ge 0$

Theorem 3.1. (Sklar’s Theorem) Let $F_X\left(x_1,\cdots ,x_n\right)$ be the joint CDF of the random variables $X_1,\cdots ,X_d$ with marginal CDFs $F_1,\cdots ,F_d$ . A d-dimensional copula $C$ exists such that the following relationship holds for every $x\in {ℝ}_d$ :

$F_X\left(x_1,\cdots ,x_d\right)=C\left(F_1\left(x_1\right),\cdots ,F_d\left(x_d\right)\right)$ (18)

The copula $C$ is said to be unique if all the individual CDFs ($F_1,\cdots ,F_d$ ) are continuous.

Theorem 3.2. Let $C$ be a multivariate Copula. For every $u_i$ in ${\left[0,1\right]}^d$ , where $i=1,\cdots ,d$ , the partial derivative $\frac{\partial C\left(u_1,\cdots ,u_d\right)}{\partial u_i}$ exist

$\frac{\partial C\left(u_1,\cdots ,u_d\right)}{\partial u_i}\ge 0$

Theorems 3.1 and 3.2 are proved in [31].

Definition 3.2. (Copula density) For each $u_1,\cdots ,u_d$ in ${\left[0,1\right]}^d$ , the copula density $c\left(u_1,\cdots ,u_d\right)$ is given by:

$c\left(u_1,\cdots ,u_d\right)=\frac{{\partial }^{n}C\left(u_1,\cdots ,u_d\right)}{\partial u_1\cdots u_d}$ (19)

Differentiating Equation (18) yields:

$f\left(x_1,\cdots ,x_d\right)={\displaystyle \prod }_{i=1}^d{f}_i\left(x_i\right)\cdot c\left(F_1\left(x_1\right),\cdots ,F_d\left(x_d\right)\right)$ (20)

where $f_i$ is the probability density function of $X_i$ for $i=1,\cdots ,d$ . Therefore, the copula density can now be given as the ratio of the joint density $f$ to the product of the individual densities:

$c\left(u_1,\cdots ,u_d\right)=\frac{f\left(u_1,\cdots ,u_d\right)}{{{\displaystyle \prod }}_{i=1}^d{f}_i\left(u_i\right)}$ (21)

Given two random variables $X_1$ and $X_2$ , a bivariate copula function $C:{\left[0,1\right]}^2\to \left[0,1\right]$ couples the marginal CDFs of the two variables to form their joint distribution. The two main bivariate copula families are the Archimedean and elliptical Copulas.

2.2.1. Elliptical Copulas

The Gaussian (normal) and Student-t Copulas are mostly used for elliptical Copulas and they model symmetric dependence structures effectively. Based on the Sklar’s Theorem, elliptical copula models’ general equation is as follows:

$F_{X_1,X_2}\left(x_1,x_2\right)=C\left(F_{X_1}\left(x_1\right),F_{X_2}\left(x_2\right)\right)$ (22)

The Gaussian copula belongs to the normal distribution and its general form is:

$\begin{array}{l}{C}_{\rho }^{Gauss}\left(F_{X_1}\left(x_1\right),F_{X_2}\left(x_2\right)\right)={\Phi }_{\rho }\left(x_1,x_2\right)\\ ={\displaystyle {\int }_{-\infty }^{x_1}{\displaystyle {\int }_{-\infty }^{x_2}\frac{1}{2\pi {\left(1-{\rho }^2\right)}^{1/2}}\exp \left[-\frac{x_1^2-2\rho x_1{x}_2+x_2^2}{2\left(1-{\rho }^2\right)}\right]\text{d}{x}_1\text{d}{x}_2}}\end{array}$ (23)

where ${\Phi }_{\rho }\left(.,.\right)$ represents the CDF of bivariate standard normal and $\rho$ is the correlation coefficient between $x_1$ and $x_2$ .

The student-t copula (or t-copula) belongs to the Student-t distribution expressed is:

$\begin{array}{l}{C}_{r,\rho }^t\left(F_{X_1}\left(x_1\right),F_{X_2}\left(x_2\right)\right)=t_{r,\rho }\left(x_1,x_2\right)\\ ={\displaystyle {\int }_{-\infty }^{x_1}{\displaystyle {\int }_{-\infty }^{x_2}\frac{1}{2\pi {\left(1-{\rho }^2\right)}^{1/2}}{\left[1+\frac{x_1^2-2\rho x_1{x}_2+x_2^2}{r\left(1-{\rho }^2\right)}\right]}^{-\frac{r+2}{2}}\text{d}{x}_1\text{d}{x}_2}}\end{array}$ (24)

where $t_{r,\rho }$ represents the CDF of the bivariate t-distribution with $r$ degrees of freedom and $\rho$ as the correlation coefficient.

Although elliptical copulas are simple to use and the simulations based on these models are simple to perform, they require numerous parameters to be estimated, generally lack closed-form expressions, and they assume symmetry, but real-world loss distributions are often skewed. These shortcomings highlight the need for more flexible alternatives such as Archimedean or vine copula constructions, which can better accommodate the complexities of real-world dependence patterns.

2.2.2. Archimedean Copula

Archimedean copulas address the symmetric tail dependence limitation elliptical copula by modelling asymmetric dependence in either the left or right tail. Archimedean copulas also have closed-form expressions and are not constructed using Equation (18), but are instead associated with the Laplace transforms of bivariate distribution functions. These copulas, which include Clayton, Frank, Joe, and Gumbel Copulas, can also capture a wide range of dependence. The construction of the Archimedean copulas relies primarily on the generator function $\phi :\left[0,1\right]\to \left[0,\infty \right]$ and the pseudo-inverse of the generating function define by:

${\phi }^{\left[-1\right]}=\{\begin{array}{ll}{\phi }^{\left[-1\right]},\hfill & 0\le t\le \phi \left(0\right)\hfill \\ 0,\hfill & \phi \left(0\right)<t<\infty \hfill \end{array}$

The general Archimedean copula is given by:

$C_{\phi }\left(F_1\left(x_1\right),F_2\left(x_2\right)\right)={\phi }^{\left[-1\right]}(\phi \left(F_1\left(x_1\right)\right)+\phi \left(F_2\left(x_2\right)\right)$ (25)

The Gumbel copula utilizes a generator function defined as $\phi ={\left(-\text{ln}\left(t\right)\right)}^{\theta }$ and is expressed as:

$C_{GU}\left(F_1\left(x_1\right),F_2\left(x_2\right)\right)=\exp \left[-{\left({\left(-\ln F_1\left(x_1\right)\right)}^{\theta }+{\left(-\ln F_1\left(x_1\right)\right)}^{\theta }\right)}^{\frac{1}{\theta }}\right]$ (26)

with parameter $\theta$ limited to the interval $\left[1,\infty \right)$ . It is asymmetric, does not account for negative dependence, exhibits upper tail-dependence (${\lambda }^{+}=2-2^{\frac{1}{\theta }}$ ) but contains no lower tail dependence.

The Clayton copula has a generating function given by $\phi =\frac{1}{\theta }\left(t^{-\theta -1}\right)$ and can be represented as:

$C_c\left(F_1\left(x_1\right),F_2\left(x_2\right)\right)={\left(F_1{\left(x_1\right)}^{-\theta +F_1{\left(x_1\right)}^{-\theta }}-1\right)}^{-1/\theta }$ (27)

with $\theta$ the dependence parameter limited in the interval $\left(0,\infty \right)$ . The Clayton Copula is asymmetric, ${\lambda }^{-}=2^{-\frac{1}{\theta }}$ , ${\lambda }^{+}=0$ and cannot account for negative dependence.

The Frank copula generator function corresponds to $\phi =-\text{ln}\frac{\text{exp}\left(-\theta t\right)-1}{\text{exp}\left(-\theta \right)-1}$ and its equation is given as:

$C_F\left(F_1\left(x_1\right),F_2\left(x_2\right)\right)=-{\theta }^{-1}\ln \left(1+\frac{\exp \left(-\theta F_1\left(x_1\right)\right)\cdot \exp \left(-\theta F_2\left(x_2\right)\right)}{\exp \left(-\theta \right)-1}\right)$ (28)

with single dependency parameter $\theta \in \left(0,\infty \right)$ . Similar to elliptical copulas, the Frank copula is symmetrical and ${\lambda }^{-}=0$ .

The Joe copula is characterized by the generating function $\phi =-\text{ln}\left[1-{\left(1-t\right)}^{\theta }\right]$ and is define as follows:

$\begin{array}{l}{C}_J\left(F_1\left(x_1\right),F_2\left(x_2\right)\right)\\ =1-\left({\left(1-F_1\left(x_1\right)\right)}^{\theta }+{\left(1-F_2\left(x_2\right)\right)}^{\theta }-{\left(1-F_1\left(x_1\right)\right)}^{\theta }{\left(1-F_2\left(x_2\right)\right)}^{\theta }\right)\end{array}$ (29)

with parameter $\theta \ge 1$ . The Joe copula is asymmetric, ${\lambda }^{+}=2-2^{\frac{1}{\theta }}$ and ${\lambda }^{-}=0$ .

The drawback of the Archimedean and Student-t copulas in tail dependence modelling is their reliance on a single parameter. This single parameter constrains them to symmetric dependence or to specific types of tail dependence that cannot be independently adjusted for upper and lower tails. As a result, they may fail to capture the asymmetric tail dependence often observed in financial markets, where extreme losses and gains can exhibit different dependency structures. This limitation reduces their effectiveness in accurately modelling risk during market stress and in constructing robust portfolios.

2.3. Copula-Based Dependence Measures

The degree to which variables $X$ and $Y$ are related is indicated by a measure of dependence. The dependence coefficients lie from −1 and 1, where the value of −1 means a perfect negative dependency, 0 means no dependency and 1 indicates a perfect positive dependency. The Pearson’s correlation coefficient is commonly used measure of dependency; however, it only focuses on linear dependencies between variables, it is sensitive to outliers, and assumes the data is homoscedastic, meaning it may be inaccurate if the variability of one variable changes with respect to the other. Therefore, copula-based dependency measures, particularly Kendall’s tau and tail dependency, offer more flexibility and accuracy by capturing non-linear dependencies and the probability of extreme co-movements. The dependence coefficients lie from −1 and 1, where the value of −1 means a perfect negative dependency, 0 means no dependency and 1 indicates a perfect positive dependency.

2.3.1. Kendall’s Tau

Kendall’s Tau is a non-parametric measure of concordance between two pairs. If there is a tendency for high values of a variable to be linked to high values of the other, then the two variables are concordant; the same is true for low values. Kendall’s tau always lies in the interval [−1, 1] and is equal to 0 if the variables are independent. The likelihood that the variables will move in the same direction is greater (lower) than the likelihood that they will move in opposite ways if Kendall’s tau is positive (negative).

Definition 3.3. (Concordance, Discordance) Consider two observations from a continuous vector $\left(X,Y\right)$ : $\left(x_i,y_i\right)$ and $\left(x_j,y_j\right)$ . Both pairs are considered concordant and discordant if $\left(x_i-x_j\right)\left(y_i-y_j\right)>0$ and $\left(x_i-x_j\right)\left(y_i-y_j\right)<0$ , respectively.

Definition 3.4. (Kendall’s $\tau$ ) Let $\left(x_1,y_1\right),\left(x_2,y_2\right),\cdots ,\left(x_n,y_n\right)$ be a sample of n observations from a vector $\left(X,Y\right)$ of continuous random variables, with at least two distinct pairs of observations, each pair can either be concordant or discordant. Assume that the concordant and discordant pairs are denoted by $c$ and $d$ , respectively, $\tau$ is calculated as:

$\tau =\frac{c-d}{c+d}$ (30)

Alternatively, $\tau$ is given as the probability of concordance less discordance as:

$\begin{array}{c}\tau =P\left(\left(X_1-X_2\right)\left(Y_1-Y_2\right)>0\right)-P\left(\left(X_1-X_2\right)\left(Y_1-Y_2\right)<0\right)\\ =2P\left(\left(X_1-X_2\right)\left(Y_1-Y_2\right)>0\right)-1\end{array}$ (31)

where $\left(X_1,Y_1\right)$ and $\left(X_2,Y_2\right)$ are iid random vectors with a joint distribution $F$ . In terms of the copula $C$ , $\tau$ is represented as follows:

$\tau =4{\displaystyle {\int }_0^1{\displaystyle {\int }_0^{1}C}}\left(F_1\left(x_1\right),F_2\left(x_2\right)\right)\text{d}C\left(F_1\left(x_1\right),F_2\left(x_2\right)\right)-1$ (32)

Table 1 presents formula and the range of $\tau$ for different bivariate copulas, where $\rho =\frac{\mathrm{cov}\left(XY\right)}{{\sigma }_x{\sigma }_y}$ is the correlation coefficient between $x$ and $y$ , $\theta$ represents the parameter in each copula, $D_1\left(\theta \right)=\frac{1}{\theta }{\displaystyle {\int }_0^{\theta }\frac{t}{\text{exp}\left(t\right)-1}\text{d}t}$ is the order one Debye function and $DG$ is the digamma function such that $DG\left(z\right)=\frac{\text{d}}{\text{d}z}\left(\ln \Gamma \left(z\right)\right)$ .

Table 1. Kendall’s tau for various bivariate copulas.

Copula	Kendall’s τ Formula	Range of τ
Gaussian	$\frac{2}{\pi }\text{arcsin}\left(\rho \right)$	[−1, 1]
Student-t	$\frac{2}{\pi }\text{arcsin}\left(\rho \right)$	[−1, 1]
Gumbel	$1-\frac{1}{\theta }$	[0, 1]
Clayton	$\frac{\theta }{\theta +2}$	[0, 1]
Frank	$1-\frac{4}{\theta }\left(1-D_1\left(\theta \right)\right)$	[−1, 1]
Joe	$1+\frac{2}{2-\theta }\left[DG\left(2\right)-DG\left(\frac{2}{\theta }+1\right)\right]$	[0, 1]

Kendall’s tau correlation is an effective metric for assessing the direction and strength of the relationship between two variables. However, it may not be the best measure of dependence, particularly when the variables’ relationship is not monotonic, there are tied ranks, or if outliers are present in the data. Additionally, while Kendall’s Tau captures overall dependence, it doesn’t specifically measure the likelihood of joint extreme values.

2.3.2. Tail Dependence Measures

Tail dependence measures, such as upper and lower tail dependence coefficients, directly assess the probability that one variable will experience an extreme value given that another does. The tail dependence coefficients measure the strength of dependency in the joint upper and lower tails of bivariate or multivariate distributions and are derived using conditional probability, hence their values are between 0 and 1.

Suppose $X_1$ and $X_2$ are random variables with $F_1$ and $F_2$ as their marginal distributions. The dependence coefficient in the lower tail is defined as

${\lambda }^{-}=\underset{q\to 0}{\text{lim}}P\left(X_2\le F_2^{-1}\left(q\right)|X_1\le F_1^{-1}\left(q\right)\right)$ (33)

provided the limit ${\lambda }^{-}\in \left[0,1\right]$ exist and $q$ is the quantile. Thus, dependency in the lower tail is the likelihood that $X_2$ does not exceed its $q^{th}$ -quantile given that $X_1$ does not exceed its $q^{th}$ -quantile. The lower tail dependency is present if ${\lambda }^{-}\in \left(0,1\right]$ , and $X_1$ and $X_2$ are asymptotically independent if ${\lambda }^{-}=0$ . Similarly, the upper tail dependence coefficient is given as

${\lambda }^{+}=\underset{q\to 1}{\text{lim}}P\left(X_2\ge F_2^{-1}\left(q\right)|X_1\ge F_1^{-1}\left(q\right)\right)$ (34)

provided ${\lambda }^{+}\in \left[0,1\right]$ exist. Thus, given that $X_1$ exceeds its $q^{th}$ -quantile, the likelihood that $X_2$ exceeds its $q^{th}$ -quantile is the upper tail dependency, also known as the extremal dependence. When $F_1$ and $F_2$ are continuous, ${\lambda }^{-}$ and ${\lambda }^{+}$ can be represented in terms of a bivariate copula $C$ as follows:

$\begin{array}{c}{\lambda }^{-}=\underset{q\to 0}{\text{lim}}\frac{P\left[X_1\le F_1^{-1}\left(q\right),X_2\le F_2^{-1}\left(q\right)\right]}{P\left[X_2\le F_2^{-1}\left(q\right)\right]}\\ =\underset{q\to 0}{\text{lim}}\frac{C\left(q,q\right)}{q}\end{array}$ (35)

The bivariate copula shows dependency in the lower tail, if ${\lambda }^{-}\in \left(0,1\right]$ and independence if ${\lambda }^{-}=0$ .

$\begin{array}{c}{\lambda }^{+}=\underset{q\to 1}{\text{lim}}\frac{P\left[X_1\ge F_1^{-1}\left(q\right),X_2\ge F_2^{-1}\left(q\right)\right]}{P\left[X_2>F_2^{-1}\left(q\right)\right]}\\ =\underset{q\to 1}{\text{lim}}\frac{S_C\left(q,q\right)}{1-q}=\underset{q\to 1}{\text{lim}}\frac{1-2q-C\left(q,q\right)}{1-q}\end{array}$ (36)

where $S_C$ represents the survival copula. The copula $C$ has dependency in the upper tail if ${\lambda }^{+}\in \left(0,1\right]$ and independent ${\lambda }^{+}=0$ . Table 2 summarizes the tail dependence characteristics of various bivariate copulas, providing an overview of each copula’s upper and lower tail dependencies.

Table 2. Upper and lower tail dependence coefficients for different bivariate copulas.

Copula	Lower tail (${\lambda }^{-}$ )	Lower tail (${\lambda }^{+}$ )
Gaussian	0	0
Student-t	$2t_{r+1}\left(-\sqrt{\frac{\left(r+1\right)\left(1-\rho \right)}{1+\rho }}\right)$	$2t_{r+1}\left(-\sqrt{\frac{\left(r+1\right)\left(1-\rho \right)}{1+\rho }}\right)$
Frank	0	0
Clayton	$2^{-\frac{1}{\theta }}$	0
Gumbel	0	$2-2^{\frac{1}{\theta }}$
Joe	0	$2-2^{\frac{1}{\theta }}$

2.4. Vine Copulas

[11], presented the pair-copula decompositions graphically using a series of nested trees with undirected edges, referred to as vine trees. Each vine copula is structured according to a Regular Vine framework, and the two specific varieties of R-vine are C-Vine (Canonical Vine) and D-Vine (Drawable Vine).

Definition 3.5. (Tree) $T=\left(N,E\right)$ has N nodes and E is the edge connecting these nodes. A node’s degree, represented as $d\left(v\right)$ , is the number of neighbouring nodes it has, $v\in N$ .

Definition 3.6. (Regular Vine (R-Vine)) A regular vine tree sequence on $d$ elements consists of a set of trees $V=\left(T_1,\cdots ,T_{d-1}\right)$ , satisfying the following conditions:

1) For each $i=1,\cdots ,d-1$ , the trees $T_i=\left(N_i,E_i\right)$ must be connected.

2) The tree $T_1$ has a node set $N_1=1,\cdots ,d$ and an edge set $E_1$ .

3) For $i\ge 2$ , the tree $T_i$ has a node set $N_i=E_{i-1}$ and an edge set $E_j$ .

4) For $i=2,\cdots ,d-1$ and $\left\{a,c\right\}\in E_i$ , it must hold that $\left|a\cup c\right|=1$ (proximity condition).

As long as the relevant edges in the preceding tree, $T_{i-1}$ , have at least one shared node, the proximity condition guarantees that two nodes in tree $T_i$ , where $i\ge 2$ , can be connected.

Figure 1 illustrates an example of 5-dimensional R-vine with 5 variable, 10 edges and 4 trees. The unconditional bivariate copula models all pairs of the variables in $T_1$ , while at the lower trees ($T_2,\cdots ,T_4$ ) dependency is modelled using the conditional bivariate copulas. For example, given the third variable, the conditional bivariate copula $c_{25|3}$ represents the dependency between the second and fifth variables. This 5-dimensional R-vine’s joint density is as follows:

$\begin{array}{c}{f}_{12345}\left(x\right)=c_{41}\left(F_4\left(x_4\right),F_1\left(x_1\right)\right)\cdot c_{13}\left(F_1\left(x_1\right),F_3\left(x_3\right)\right)\cdot c_{32}\left(F_3\left(x_3\right),F_2\left(x_2\right)\right)\\ \cdot c_{35}\left(F_3\left(x_3\right),F_5\left(x_5\right)\right)\cdot c_{43|1}\left(F\left(x_4|x_1\right),F\left(x_3|x_1\right)\right)\\ \cdot c_{12|3}\left(F\left(x_1|x_3\right),F\left(x_2|x_3\right)\right)\cdot c_{25|3}\left(F\left(x_2|x_3\right),F\left(x_5|x_3\right)\right)\\ \cdot c_{42|13}\left(F\left(x_4|x_1,x_3\right),F\left(x_2|x_1,x_3\right)\right)\cdot c_{15|23}\left(F\left(x_1|x_2,x_3\right),F\left(x_5|x_2,x_3\right)\right)\\ \cdot c_{45|132}\left(F\left(x_4|x_1,x_3,x_2\right),F\left(x_5|x_{21},x_3,x_2\right)\right)\cdot {\displaystyle \prod }_{i=1}^5{f}_i\left(x_i\right)\end{array}$ (37)

Definition 3.7. An R-vine can be categorized into two specific types based on its structure:

1) C-Vine, if every tree $T_i$ contains a distinct node connected by $d-i$ edges. The node in tree $T_1$ with $d-1$ edges is referred to as the root node.

2) D-Vine, if $T_1$ nodes’ have a maximum of two edges.

Canonical vines (C-vines)

A d-dimensional C-vine copula has $\frac{d!}{2}$ possible tree structures. For example,

for a 5-dimensional, there are 60 decompositions to order the variables.

Figure 1. An R-vine tree representation for $d=5$ .

Figure 2 shows a C-vine copula with each node representing a variable and each edge corresponds to unconditional or conditional bivariate copula family. It has 4 trees, $T_i,i=1,2,3,4$ , each tree has $5-i$ edges and $6-i$ nodes and the root node is variable 4. This 5-dimensional C-vine has the joint density expressed as:

$\begin{array}{c}{f}_{12345}\left(x\right)=c_{41}\left(F_4\left(x_4\right),F_1\left(x_1\right)\right)\cdot c_{42}\left(F_4\left(x_4\right),F_2\left(x_2\right)\right)\cdot c_{43}\left(F_4\left(x_4\right),F_3\left(x_3\right)\right)\\ \cdot c_{45}\left(F_4\left(x_4\right),F_5\left(x_5\right)\right)\cdot c_{12|4}\left(F\left(x_1|x_4\right),F\left(x_2|x_4\right)\right)\\ \cdot c_{13|4}\left(F\left(x_1|x_4\right),F\left(x_3|x_4\right)\right)\cdot c_{15|4}\left(F\left(x_1|x_4\right),F\left(x_5|x_4\right)\right)\\ \cdot c_{23|41}\left(F\left(x_2|x_4,x_1\right),F\left(x_3|x_4,x_4\right)\right)\cdot c_{25|41}\left(F\left(x_2|x_4,x_1\right),F\left(x_5|x_4,x_1\right)\right)\\ \cdot c_{35|412}\left(F\left(x_3|x_4,x_1,x_2\right),F\left(x_5|x_4,x_1,x_2\right)\right)\cdot {\displaystyle \prod }_{i=1}^5{f}_i\left(x_i\right)\end{array}$ (38)

Drawable vines (D-vines)

Similarly to C-vines, there are $\frac{d!}{2}$ possible D-vine decompositions. However,

unlike C-vines, the D-vine copulas structure does not require a root variable. The variables are linked sequentially, one following the other, with unconditional bivariate copulas used to model dependency between variables in the tree one and conditional pair-copulas employed in the subsequent trees where the conditioning sets in the D-vine structure are different.

Figure 2. Tree structure of a C-vine for $d=5$ .

Figure 3 represents a five-dimensional drawable vine copula with four trees, fourteen nodes, and ten edges. This five-dimensional D-vine copula’s joint PDF is as follows:

$\begin{array}{c}{f}_{12345}\left(x\right)=c_{42}\left(F_4\left(x_4\right),F_2\left(x_2\right)\right)\cdot c_{23}\left(F_2\left(x_2\right),F_3\left(x_3\right)\right)\cdot c_{35}\left(F_3\left(x_3\right),F_5\left(x_5\right)\right)\\ \cdot c_{51}\left(F_5\left(x_5\right),F_1\left(x_1\right)\right)\cdot c_{43|2}\left(F\left(x_4|x_2\right),F\left(x_3|x_2\right)\right)\\ \cdot c_{25|3}\left(F\left(x_2|x_3\right),F\left(x_5|x_3\right)\right)\cdot c_{31|5}\left(F\left(x_3|x_5\right),F\left(x_1|x_5\right)\right)\\ \cdot c_{45|23}\left(F\left(x_4|x_2,x_3\right),F\left(x_5|x_2,x_3\right)\right)\cdot c_{21|35}\left(F\left(x_2|x_3,x_5\right),F\left(x_1|x_3,x_5\right)\right)\\ \cdot c_{41|235}\left(F\left(x_4|x_2,x_3,x_5\right),F\left(x_1|x_2,x_3,x_5\right)\right)\cdot {\displaystyle \prod }_{i=1}^5{f}_i\left(x_i\right)\end{array}$ (39)

Figure 3. Tree structure of a D-vine for $d=5$ .

2.4.1. Stationary Vine Copulas

The S-vine is a extension of the D-vine that assumes stationarity in the dependence structure. To simplify the expression of the S-vine, strict stationarity conditions are imposed on the D-vine [20]. These conditions ensure that dependency structure between variables does not change with time shifting. If $C\left(V\right)$ represents a vine copula, then it is considered to be translation invariant if the partial copulas associated with any edge $E_i$ , $i=1,\cdots ,d$ are identical to those of the edge $E_i+s$ , where $s$ represents a time shift by $s$ steps. According to the translation invariance condition, a D-vine copula in Figure 3 transitions into an S-vine if the copulas illustrating the relationships between the pairs (4,3), (2,3), (3,5), (5,1) and partial copulas $c_{43|2},c_{25|3},c_{31|5},c_{45|23},c_{21|35}$ are identical. Therefore, the S-vine has of four-pair copulas, offering a streamlined alternative to the D-vine configuration, which requires ten-pair copulas.

2.4.2. Parameter Estimation of Vine Copulas

The parameters of the vine copulas are estimated using the Full Maximum Likelihood (FML) and sequential estimation method. Consider a three dimension vine copula with the density function given as:

$\begin{array}{c}{f}_{123}\left(x|\theta \right)=c_{31}\left(F_3\left(x_3\right),F_1\left(x_1\right)\right)\cdot c_{12}\left(F_1\left(x_1\right),F_2\left(x_2\right)\right)\\ \cdot c_{32|1}\left(F\left(x_3|x_1\right),F\left(x_2|x_1\right)\right)\cdot {\displaystyle \prod }_{i=1}^3{f}_i\left(x_i\right)\end{array}$ (40)

The FML method estimates all vine structure parameters simultaneously. The joint log-likelihood for the three-dimensional vine copula is given as:

$\begin{array}{c}l\left(\theta \right)={\displaystyle \sum }_{t=1}^T\log c_{31}\left(F_3\left(x_3\right),F_1\left(x_1\right);{\theta }_{31}\right)+{\displaystyle \sum }_{t=1}^T\log c_{12}\left(F_1\left(x_1\right),F_2\left(x_2\right);{\theta }_{12}\right)\\ +{\displaystyle \sum }_{t=1}^T\log c_{32|1}\left(F\left(x_3|x_1\right),F\left(x_2|x_1\right);{\theta }_{32|1}\right)\end{array}$ (41)

where $\theta =\left({\theta }_{31},{\theta }_{12},{\theta }_{32|1}\right)$ represents all the parameters that need to be estimated. By maximising $l\left(\theta \right)$ and solving

$\frac{\partial l\left(\theta \right)}{\partial \theta }=0$ (42)

$\widehat{\theta }$ is obtained. However, in practice for some high dimensional cases, it is computationally expensive.

The sequential estimation method is a stepwise approach that estimates parameters for all pair copulas in the vine structure sequentially. This makes it computationally efficient and suitable for high-dimensional datasets. In step one, ${\theta }_{31}$ and ${\theta }_{12}$ are separately estimated, by maximizing $l\left(\theta \right)$ in Equation (41) by using the bivariate sub-sample $\left(F_1\left(x_{1t}\right),F_2\left(x_{2t}\right)\right)$ and $\left(F_2\left(x_{2t}\right),F_3\left(x_{3t}\right)\right)$ for $t=1,\cdots ,T$ , respectively. The corresponding estimates are given by:

$\begin{array}{l}{\widehat{\theta }}_{31}=\arg \underset{{\theta }_{31}}{\max }{\displaystyle \sum }_{t=1}^T\log c_{31}\left(F_3\left(x_3\right),F_1\left(x_1\right);{\theta }_{31}\right)\\ {\widehat{\theta }}_{12}=\arg \underset{{\theta }_{12}}{\max }{\displaystyle \sum }_{t=1}^T\log c_{12}\left(F_1\left(x_1\right),F_2\left(x_2\right);{\theta }_{12}\right)\end{array}$ (43)

In the second step, the pseudo conditional copula observations are formed as $F\left(x_{3,t|1,{\widehat{\theta }}_{31}}\right)=C_{3|1}\left(F\left(x_{3,t}|x_{1,t},{\widehat{\theta }}_{31}\right)\right)$ and $F\left(x_{1,t|2,{\widehat{\theta }}_{12}}\right)=C_{1|2}\left(F\left(x_{1,t}|x_{2,t},{\widehat{\theta }}_{12}\right)\right)$ for

$t=1,\cdots ,T$ to be used to estimate ${\theta }_{32|1}$ by maximizing

$l\left({\theta }_{13|2}\right)={\displaystyle \sum }_{t=1}^T\log c_{32|1}\left(F\left(x_3|x_1\right),F\left(x_2|x_1\right);{\theta }_{32|1}\right)$ (44)

over ${\theta }_{32|1}$ . Thus, ${\widehat{\theta }}_{32|1}$ is given as

${\widehat{\theta }}_{32|1}=\text{arg}\underset{{\theta }_{32|1}}{\text{max}}l\left({\theta }_{32|1}\right)$ (45)

2.5. Vine Copula Multivariate GARCH Model

This section introduces hybrid models that combine MGARCH models with vine copulas. The MGARCH will model the time-varying volatilities and conditional correlation and vine copula will model the dependency structure separately but concurrently for non-normal distributions. The VC-MGARCH model eliminates linear correlation from the variables and forms dependent but uncorrelated innovations controlled by a vine copula, while MGARCH controls the correlation. The VC-MGARCH models are defined in a identical way to convectional MGARCH models, but the multivariate distribution function of the residuals or innovations are modelled using a vine copula function.

2.5.1. Structure of VC-MGARCH Model

For $d$ assets, let $r_t={\left(r_{1,t},\cdots ,r_{d,t}\right)}^{\prime }$ be the log-returns vector, ${\eta }_t={\left({\eta }_{d,t},\cdots ,{\eta }_{d,t}\right)}^{\prime }$ be the uncorrelated dependent errors and ${\xi }_t={\left({\xi }_{1,t},\cdots ,{\xi }_{d,t}\right)}^{\prime }$ be random variables from the copula distribution, such that:

${\xi }_t|{ℰ}_{t-1}~F_{1,\cdots ,d}\left({\xi }_t\right)=C_{1,\cdots ,d}\left(F_1\left({\xi }_{1,t}\right),\cdots ,F_d\left({\xi }_{d,t}\right);{\theta }_t\right)$

where ${\theta }_t$ is the parameter associates with the copula $C\left(.,.\right)$ . Instead of assuming a simple multivariate normal distribution for ${\eta }_t$ in Equation (3), [22], proposed the transformation of these residuals using a copula-based structure:

${\eta }_t={\Sigma }_t^{-1/2}{\xi }_t$ (46)

where ${\xi }_t$ are the transformed residuals that follow copula distribution, with covariance matrix $\mathbb{E}\left[\xi {\xi }^{\prime }|F_{t-1}\right]={\Sigma }_t={\left[{\sigma }_{ij,t}\right]}_{i,j=1}^d$ , and $\mathbb{E}\left[{\xi }_t|F_{t-1}\right]=0$ . The VC-MGARCH model assumes a dependent copula for ${\xi }_t$ while keeping ${\eta }_t$ uncorrelated, that is

$C\left(u_1,\cdots ,u_d\right)\ne {\displaystyle {\prod }_{i=1}^d{u}_i}$ and ${\sigma }_{ij}\ne 0$

where $u_i=F_i\left({\xi }_i\right)$ . In contrast, the convectional MGARCH models assume that ${\eta }_t$ is independent and normally distributed such that $C\left(u_1,\cdots ,u_d\right)={\displaystyle {\prod }_{i=1}^d{u}_i}$ and ${\sigma }_{ij}=0$ . If ${\xi }_i$ and ${\xi }_j$ are two random variables with CDFs $F_i$ and $F_j$ , joint distribution $F_{ij}$ and the first two moments are well-defined and finite, then ${\sigma }_{ij}$ can be obtained using the Hoeffding Lemma, given by

${\sigma }_{ij}\equiv \mathrm{cov}\left({\xi }_i{\xi }_j\right)={\displaystyle {\iint }_{R^2}\left[F_{ij}\left({\xi }_i,{\xi }_j\right)-F\left({\xi }_i\right)F\left({\xi }_j\right)\right]\text{d}{\xi }_i\text{d}{\xi }_j}$ (47)

Using the Sklar’s theorem, ${\sigma }_{ij}$ is given by;

${\sigma }_{ij,t}\left({\theta }_t\right)={\displaystyle {\iint }_{R^2}\left[C_{ij}\left(F_i\left({\xi }_i\right),F_j\left({\xi }_j\right);{\theta }_{ij,t}\right)-F\left({\xi }_i\right)F\left({\xi }_j\right)\right]\text{d}{\xi }_i\text{d}{\xi }_j}$ (48)

For $d=3$ , the random variables ${\xi }_1$ , ${\xi }_2$ and ${\xi }_3$ have marginal CDFs given as $F_1$ , $F_2$ , and $F_3$ , then the vine copula joint CDF is as follows;

$F_{123}\left({\xi }_{1,t},{\xi }_{2,t},{\xi }_{3,t}\right)={\displaystyle {\int }_{-\infty }^{{\xi }_{1,t}}{C}_{32|1}}\left(F_{3|1}\left({\xi }_{3,t}|z_1;{\theta }_{32}\right),F_{2|1}\left({\xi }_{2,t}|z_1;{\theta }_{12}\right);{\theta }_t\right)\text{d}{F}_1\left({\xi }_1\right)$ (49)

where $C_{32|1}(.)$ is the conditional copula, $F_{3|1}$ and $F_{2|1}$ is the conditional CDF of ${\xi }_3$ and ${\xi }_2$ , respectively, given ${\xi }_1$ and ${\theta }_{32}$ and ${\theta }_{12}$ are the parameters defining the dependency between ${\xi }_3$ and ${\xi }_2$ and between ${\xi }_1$ and ${\xi }_2$ , respectively. The 3-dimensional vine copula joint density function is

$\begin{array}{c}{f}_{123}\left({\xi }_{1,t},{\xi }_{2,t},{\xi }_{3,t}\right)={\displaystyle \prod }_{i=1}^3\left[f\left({\xi }_i\right)\right]\cdot c_{31}\left(F_3\left({\xi }_3\right),F_1\left({\xi }_1\right)\right)\cdot c_{12}\left(F_1\left({\xi }_1\right),F_2\left({\xi }_2\right)\right)\\ \cdot c_{32|1}\left(F_{3|1}\left({\xi }_3|{\xi }_1\right),F_{2|1}\left({\xi }_2|{\xi }_1\right)\right)\end{array}$ (50)

where $f_i\left({\xi }_i\right)$ are the marginal densities of the residuals ${\xi }_i$ , $c_{31}$ and $c_{12}$ are unconditional copula densities and $c_{32|1}$ is the conditional copula density.

The $d$ -dimensional vine copula joint PDF will be given by:

$\begin{array}{l}f\left({\xi }_1,\cdots ,{\xi }_d\right)={\displaystyle \prod _{i=1}^d{f}_i\left({\xi }_i\right)}{\displaystyle \prod _{j=1}^{d-1}{\displaystyle \prod _{i=1}^{d-j}{c}_{i,i+j|i+1,\cdots ,i+j-1}(F_{i|i+1,\cdots ,i+j-1}\left({\xi }_i|{\xi }_{i+1,\cdots ,i+j-1}\right),}}\\ F_{i+j|i+1,\cdots ,i+j-1}\left({\xi }_{i+1}|{\xi }_{i+1,\cdots ,i+j-1}\right))\end{array}$ (51)

where $c_{i,i+j|i+1,\cdots ,i+j-1}$ represent the conditional copula densities for the pair of residuals ${\xi }_i$ and ${\xi }_{i+j}$ , and $F_{i+j|i+1,\cdots ,i+j-1}$ is the conditional density functions. The density function of MGARCH errors ${ϵ}_t$ can be derived by rearranging ${ϵ}_t=H_t^{1/2}{\eta }_t=H_t^{1/2}{\Sigma }_t^{-1/2}{\xi }_t$ into ${\xi }_t={\Sigma }_t^{1/2}{H}_t^{-1/2}{ϵ}_t$ and utilizing the change of variable technique. Then, the density function of ${ϵ}_t$ is

$\begin{array}{c}g\left({ϵ}_t\right)=f\left({\Sigma }_t^{1/2}{H}_t^{-1/2}{ϵ}_t\right)\cdot \text{det}\left({\Sigma }_t^{1/2}{H}_t^{-1/2}\right)\\ =f_{1,\cdots ,d}\left({\xi }_1,\cdots ,{\xi }_d\right)\cdot \left|{\Sigma }_t^{1/2}{H}_t^{-1/2}\right|\end{array}$ (52)

For $d=3$ ,

$\begin{array}{c}g\left({ϵ}_t\right)=f_{123}\left({\xi }_{1,t},{\xi }_{2,t},{\xi }_{3,t}\right)\cdot \text{det}\left({\Sigma }_t^{1/2}{H}_t^{-1/2}\right)\\ ={\displaystyle \prod }_{i=1}^3\left[f\left({\xi }_i\right)\right]\cdot c_{31}\left(F_3\left({\xi }_3\right),F_1\left({\xi }_1\right)\right)\cdot c_{12}\left(F_1\left({\xi }_1\right),F_2\left({\xi }_2\right)\right)\\ \cdot c_{32|1}\left(F_{3|1}({\xi }_3|{\xi }_1),F_{2|1}({\xi }_2|{\xi }_1)\right)\cdot \left|{\Sigma }_t^{1/2}{H}_t^{-1/2}\right|\end{array}$ (53)

where $\left|{\Sigma }_t^{1/2}{H}_t^{-1/2}\right|$ denotes the Jacobian associated with the transformation from ${\xi }_t$ to ${ϵ}_t$ .

2.5.2. Parameter Estimation of VC-MGARCH Model

The VC-MGARCH model’s parameters are estimated using the Maximum Likelihood Estimation method. Let $\beta$ be the MGARCH parameters, which are used to parameterize, while $\theta$ represent the parameters of the vine copula. The log-likelihood of ${ϵ}_t$ is

$\begin{array}{c}{l}_{1,\cdots ,d}^{{ϵ}_t}\left(\beta ,\theta \right)=l_{1,\cdots ,d}^{{\xi }_t}\left(\theta \right)+\text{ln}\left|{\Sigma }_t^{1/2}\left(\theta \right)H_t^{-1/2}\left(\alpha \right)\right|\\ ={\displaystyle \sum }_{t=1}^T\ln f_{1,\cdots ,d}\left({\xi }_1,\cdots ,{\xi }_d;\theta \right)+\text{ln}\left|{\Sigma }_t^{1/2}\left(\theta \right)H_t^{-1/2}\left(\alpha \right)\right|\end{array}$ (54)

For $d=3$ , the log-likelihood associated with ${ϵ}_t$ is expressed as follows:

$\begin{array}{c}{l}_{123}^{{ϵ}_t}\left(\beta ,\theta \right)={\displaystyle \sum }_{t=1}^T{\displaystyle \sum }_{i=1}^d\ln f\left({\xi }_i\right)+{\displaystyle \sum }_{t=1}^T\ln c_{12}\left(F_1\left({\xi }_1\right),F_2\left({\xi }_2\right)\right)\\ +{\displaystyle \sum }_{t=1}^T\ln c_{32}\left(F_3\left({\xi }_3\right),F_2\left({\xi }_2\right)\right)+{\displaystyle \sum }_{t=1}^T\ln c_{32|1}\left(F_{3|1}\left({\xi }_3|{\xi }_1\right),F_{2|1}\left({\xi }_2|{\xi }_1\right)\right)\\ +\text{ln}\left|{\Sigma }_t^{1/2}\left(\theta \right)H_t^{-1/2}\left(\alpha \right)\right|\end{array}$ (55)

The log-likelihood functions are maximized to get the corresponding estimates $\widehat{\theta }$ and $\widehat{\beta }$ :

$\left(\widehat{\beta },\widehat{\theta }\right)=\text{arg}\underset{\beta ,\theta }{\text{max}}{l}_{123}^{{ϵ}_t}\left(\beta ,\theta \right)$ (56)

2.6. Forecasting Portfolio Market Risk

In this section, the risk measures for forecasting market risk, Value at Risk and Expected Shortfall, are presented.

2.6.1. Value at Risk

One commonly used metric for evaluating portfolio risks is VaR. With a probability of $1-\alpha$ , it calculates the possible loss of a portfolio $Z$ that is not exceeded, that is

$Va{R}_{\alpha }\left(Z\right)=\text{min}\left\{z:P(Z>z)\le \alpha \right\}=\text{min}\left\{z:P\left(Z\le z\right)\ge 1-\alpha \right\}=F_Z^{-1}\left(1-\alpha \right)$ (57)

where $F_Z$ is the CDF of $Z$ and $F_Z^{-1}$ is the standardized quantile of $F_Z$ . VaR is a quantile-based risk metric that has several benefits, including the usage of a single value (expressed in monetary amounts or percentages) for the specified level of risk and ease of comparison and interpretation. However, two primary criticisms of VaR have emerged. Firstly, the sub-additivity property is violated, which means that the combined portfolio’s VaR may exceed the total individual VaRs of its components. Secondly, when the underlying distribution contains heavy tails, VaR may underestimate risk since it only shows the smallest loss that is not surpassed with a probability of $1-\alpha$ . In such cases, the VaR may be much smaller than the projected loss in the $\alpha$ -quantile. ES was proposed to address this challenge by mitigate the described underestimation of risks associated with VaR.

2.6.2. Expected Shortfall

Expected Shortfall focuses on the tail portion when losses exceed the VaR threshold. ES is the average losses over $Va{R}_{\alpha }\left(Z\right)$ given the confidence level ($1-\alpha$ ) and time $t$ . It is given as follows:

$\begin{array}{c}E{S}_{\alpha }\left(Z\right)=\mathbb{E}\left(Z|Z\ge Va{R}_{\alpha }\left(Z\right)\right)=\frac{{\displaystyle {\int }_{Va{R}_{\alpha }\left(Z\right)}^{\infty }zf\left(z\right)\text{d}z}}{P\left(Z>Va{R}_{\alpha }\left(Z\right)\right)}\\ =\mathbb{E}\left(Z+Va{R}_{\alpha }\left(Z\right)-Va{R}_{\alpha }\left(Z\right)|Z\ge Va{R}_{\alpha }\left(Z\right)\right)\\ =\mathbb{E}\left(Va{R}_{\alpha }\left(Z\right)|Z\ge Va{R}_{\alpha }\left(Z\right)\right)+\mathbb{E}\left(Z-Va{R}_{\alpha }\left(Z\right)|Z\ge Va{R}_{\alpha }\left(Z\right)\right)\\ =Va{R}_{\alpha }\left(Z\right)+\mathbb{E}\left(Z-Va{R}_{\alpha }\left(Z\right)|Z\ge Va{R}_{\alpha }\left(Z\right)\right)\end{array}$ (58)

If $u=F_Z$ for $Va{R}_{\alpha }\left(Z\right)\le z\le \infty$ , then $\text{d}u=f\left(z\right)\text{d}z$ , $F\left(Va{R}_{\alpha }\left(Z\right)\right)=\alpha$ , $F\left(\infty \right)=1$ and $z=F^{-1}\left(u\right)=Va{R}_u$ . Then, Equation (58) becomes:

$E{S}_{\alpha }\left(Z\right)=\frac{1}{1-\alpha }{\displaystyle {\int }_{\alpha }^{1}Va{R}_u\text{d}u}$ (59)

2.6.3. Forecasting Value at Risk and Expected Shortfall Using MGARCH-Vine Copula Model

In forecasting the market risk measures, VaR and ES, of the portfolio of $d$ assets using the MGARCH-Vine Copula model, the following steps were used:

1) Compute the return vector $r_t={\left(r_{1t},\cdots ,r_{dt}\right)}^{\prime }$ , $r_{it}=\text{log}\left(\frac{P_{i,t}}{P_{i,t-1}}\right)$ , for $i=1,\cdots ,d$

and $t=1,\cdots ,T$ .

2) Fit the MGARCH models to the returns of the portfolio and extract the standardized residuals ${\eta }_t=\left({\eta }_{1,t},\cdots ,{\eta }_{d,t}\right)$

${\eta }_t=H_t^{-1/2}{\epsilon }_t$ (60)

3) Transform the standardized residuals to uniform distribution such that, $u_{it}=F_i\left({\eta }_t\right)~U{\left[0,1\right]}^d$ .

4) Fit a vine copula model to the marginal residuals $u_t={\left(u_{1,1},\cdots ,u_{d,t}\right)}^{\prime }$ obtained in step 3.

5) Simulate $M$ uniform random numbers, ${\widehat{u}}_{t+1}^{\left(m\right)}=\left({\widehat{u}}_{1,t}^{\left(m\right)},\cdots ,{\widehat{u}}_{d,t}^{\left(m\right)}\right)$ , using the estimate vine copula model and apply the inverse probability transformations of the marginal distributions to convert the simulated marginals into standardized residuals ${\widehat{\xi }}_{m,t+1}=\left({\widehat{\xi }}_{t+1}^{\left(m\right)},\cdots ,{\widehat{\xi }}_{d,t+1}^{\left(m\right)}\right)$ , for ${\widehat{\xi }}_{i,t+1}^{\left(m\right)}=F_i^{-1}\left({\widehat{u}}_{i,t}^{\left(m\right)}\right)$ .

6) Compute the return forecasts as

${\widehat{r}}_{t+1}^{\left(m\right)}={\widehat{\mu }}_t+{\widehat{H}}_t^{1/2}{\widehat{\Sigma }}_t^{-1/2}{\widehat{\xi }}_{t+1}^{\left(m\right)}$ (61)

7) The estimated return of the portfolio, denoted as ${\widehat{r}}_{t+1,P}$ , is expressed as

${\widehat{r}}_{P,t+1}=w^{\prime }{\widehat{r}}_{t+1}^{\left(m\right)}$ (62)

where $w={\left(w_1,\cdots ,w_d\right)}^{\prime }$ is the vector of portfolio weights and ${\displaystyle {\sum }_{i=1}^d{w}_i}=1$ .

8) Finally, compute the $Va{R}_{\alpha }$ and $E{S}_{\alpha }$ forecast

$Va{R}_{t+1|t}\left(\alpha \right)=F_{{\widehat{r}}_{P,t+1}}^{-1}\left(1-\alpha \right)$ (63)

$E{S}_{t+1|t}\left(\alpha \right)=\mathbb{E}\left({\widehat{r}}_{P,t+1}|{\widehat{r}}_{P,t+1}\ge Va{R}_{t+1|t}\left(\alpha \right)\right)$ (64)

where $1-\alpha$ is the confidence interval and $\alpha$ is the significance level.

2.7. Backtesting Techniques for Portfolio Market Risk Measures Based on Value-at-Risk and Expected Shortfall

Backtesting techniques were used to evaluate how well a risk model performs by checking its accuracy and comparing different model specifications. Several methods are available for backtesting VaR and ES, each with its own merits and demerits. For VaR backtesting, two primary methodologies are typically used: the Unconditional and Conditional Coverage Test. The VaR coverage test is built upon the idea of the violation process. This implies that if the financial loss on a particular day exceeds the VaR forecast, it is considered a violation of the VaR threshold. This process is given by

$1_{t+1}\left(\alpha \right)=\{\begin{array}{ll}1,\hfill & \text{for}{x}_{t,t+1}\le -Va{R}_T\left(\alpha \right)\hfill \\ 0,\hfill & \text{for}{x}_{t,t+1}\ge -Va{R}_T\left(\alpha \right)\hfill \end{array}$ (65)

where, $x_{t,t+1}$ represents the realized portfolio loss from time $t$ to $t+1$ . Using Equation (65), the forecasted VaR values are converted into a vector of 1’s and 0’s, where 1 indicates an exceedance.

2.7.1. Kupiec’s Unconditional Coverage Test

The Unconditional Coverage test, by [32] is often referred to as the proportion of failures test. It compares the actual exceedance rate with the expected rate, based on the assumption that exceedances occur independently and follow a Bernoulli process with $\left(1-\alpha \right)$ success probability. The null hypothesis $H_0:\mathbb{E}\left[I_{t+1}\right]=\left(1-\alpha \right)$ while the alternative hypothesis $H_A:\mathbb{E}\left[I_{t+1}\right]\ne \left(1-\alpha \right)$ . The test statistic for UC test is:

$L{R}_{uc}=2\ln \left[{\left(\frac{x}{n}\right)}^x{\left(1-\frac{x}{n}\right)}^{n-x}\right]-2\ln \left[{\alpha }^{n-x}{\left(1-\alpha \right)}^x\right]$ (66)

where $1-\alpha$ represents the anticipated exceedance probability, $x$ indicates the actual number of exceedances observed, and $n$ refers to the sample size. This statistic follows an asymptotic ${\chi }_1^2$ distribution and large values lead to rejection of the null hypothesis.

2.7.2. Christoffersen’s Conditional Coverage Test

The conditional coverage test, proposed by [33], extends the unconditional test by examining both the frequency and independence of exceedances. It combines a Markov independence test with the unconditional test to determine if exceedances occur randomly over time. The independence test uses a likelihood ratio test to evaluate if the chance of a violation in the VaR on a specific day is influenced by the outcome of the previous day.

Let $n_{ij}$ , $i$ and $j=0,1$ , denote the total observations where state j occurred on a single day, assuming that state i occurred on the day before. For instance, $n_{01}$ specifies the number of days for which $I_{t+1}=1$ and $I_t=0$ .

Also, let ${\pi }_{01}=\frac{n_{01}}{n_{00}+n_{01}}$ be the likelihood of a VaR breach occurring on the next day (1), assuming no breach occurred on the current day (0), ${\pi }_{11}=\frac{n_{11}}{n_{10}+n_{11}}$ the conditional probability of having a violation tomorrow given that today there is also a violation and ${\pi }_{11}+{\pi }_{01}=\pi$ . The null hypothesis $H_0:{\pi }_{01}={\pi }_{11}$ and the alternative hypothesis $H_A:{\pi }_{01}\ne {\pi }_{11}$ . The independence test’s test statistic is:

$L{R}_{ind}=2\ln \left[{\left(1-{\pi }_{01}\right)}^{n_{00}}{\pi }^{n_{01}}{\left(1-{\pi }_{11}\right)}^{n_{10}}{\pi }_{11}^{n_{11}}\right]-2\ln \left[{\left(1-\pi \right)}^{n_{00}+n_{10}}{\pi }^{n_{01}+n_{11}}\right]~{\chi }_1^2$ (67)

where $1-{\pi }_{01}={\pi }_{00}$ is the conditional probability no violation followed by non-violation and $1-{\pi }_{11}={\pi }_{10}$ is the conditional likelihood of non-violation occurring after a VaR breach.

Thus, the conditional coverage’s test statistic is:

$L{R}_{cc}=L{R}_{ind}+L{R}_{uc}~{\chi }_2^2$ (68)

The VaR estimates are considered appropriate if the test statistic is less than the critical value of the ${\chi }_2^2$ distribution, in which case the null hypothesis is not rejected.

The VaR estimates are appropriate if the statistic value is less than the ${\chi }_2^2$ distribution’s critical value, thus $H_0$ is not rejected.

2.7.3. Exceedance Residuals Tests

To test the ES estimates, [34], introduced a back testing procedure based on exceedance residuals. This method relies on the idea that the realized return $r_{t+1}$ for the next period exceeds the $E{S}_{\alpha }^{t+1}$ at time $t$ , given that $r_{t+1}$ exceeds the $Va{R}_{\alpha }^{t+1}$ :

$e_{t+1}=\frac{r_{t+1}-\widehat{E}{S}_{\alpha }^{t+1}}{{\widehat{\sigma }}_{t+1}}$ (69)

where ${\widehat{\sigma }}_{t+1}$ is the conditional volatility estimate. The null hypothesis $\mathbb{E}\left(e_{t+1}\right)=0$ against a one-sided alternative hypothesis $\mathbb{E}\left(e_t\right)<0$ . To test if $H_0$ holds a non-parametric bootstrap test proposed by [35] can be used. If $H_0$ is rejected, it indicates that the risk is underestimated.

2.7.4. Expected Shortfall Regression Tests

The ES regression backtest, introduced by [36], is an alternative ES backtesting method that only requires ES estimates as the inputs. The test uses a regression framework to determine if ES forecasts are unbiased and effective predictors of actual losses. The ES regression backtesting test is based on the regression model:

$r_t=\delta +\gamma \left(E{S}_t^{\alpha }\right)+{\epsilon }_t$ (70)

where $r_t$ is the actual return observed at time $t$ , $E{S}_t^{\alpha }$ denotes the forecasted expected shortfall, $\delta$ and $\gamma$ are the intercept term and slope coefficient, respectively, and $e_t$ is the error term. The hypotheses of this test are: $H_0:\delta =0$ and $\gamma =1$ against $H_A:\delta \ne 0$ or $\gamma \ne 1$ . To test $H_0$ , the Wald-type test statistic $W$ is used:

$W=\left(\left(\widehat{\delta },\widehat{\gamma }\right)-\left(0,1\right)\right){\widehat{\Omega }}^{-1}{\left(\left(\widehat{\delta },\widehat{\gamma }\right)-\left(0,1\right)\right)}^{\prime }$ (71)

where the estimated parameters’ covariance matrix is $\widehat{\Omega }$ and $W~{\chi }_q^2$ , in this case $q=2$ . If $W$ exceeds the critical value of the ${\chi }_q^2$ distribution at the chosen significance level, the null hypothesis is rejected, indicating that the ES forecasts may be misspecified.

3. Empirical Results

3.1. Data description

The dataset used in this study consists of a portfolio of twelve assets that include: two cryptocurrencies (Bitcoin (BTC) and Ethereum (ETH)), three currency exchange rates against the US Dollar (the Euro (EUR), the Great British Pound (GBP), and the Japanese Yen (JPY)), three commodities (gold, silver, and natural gas), and four stock indices (the Deutscher Aktienindex (DAX), the Standard and Poor’s 500 (S&P500), the National Association of Securities Dealers Automated Quotations (NASDAQ), and the Financial Times Stock Exchange 100 (FTSE 100)). The data set was downloaded from https://www.investing.com/markets/. The dataset contains daily average closing prices from January 1, 2004, to December 31, 2024, for stock indices, commodity prices, and foreign exchange rates, resulting in a total of 4341 observations, excluding weekends and public holidays. For the cryptocurrencies, the dataset includes 2557 daily average prices starting from January 1, 2018, to December 31, 2024. The selected samples encompass both the 2008 Global Financial Crisis (GFC) and the COVID-19 pandemic periods. All hypothesis tests in this study are conducted at the conventional 5% significance level, which is widely accepted in empirical financial econometrics. This significance level strikes a practical balance between minimizing false positives, which could overstate the risk linkages between assets, and ensuring that real effects are detected, such as volatility clustering and tail dependence, which are critical for accurate risk assessment.

The dataset was first analyzed by plotting the historical price series of the twelve selected financial assets over the sample period from 2004 to 2024. Figure 4 represents the price movements of major equity indices, exchange rates, cryptocurrencies and commodities over the sample period. The stock indices (S&P 500, DAX, NASDAQ, FTSE) general exhibit volatility fluctuations and clustering over the period with a sharp decline in 2008 attributed to the Global Financial Crisis and early 2020 also the COVID-19 global pandemic. Subsequently, a steady recovery post COVID-19 period is recorded as the global economy was showing signs of recovery from the effects of the pandemic. The exchange rates prices also fluctuate over time, with a notable decline in 2022 due to hikes in the Federal Reserve rates which strengthened the US dollar against the EUR and GBP, while JPY prices trends upwards from 2021, reflecting yen depreciation against the USD. The prices of gold and silver increased in 2020 as safe-haven assets, and natural gas prices went up and were more volatile in 2022, associated to the increased demand and reduction in supply of the liquefied natural gas following the Russian-Ukraine war. The cryptocurrencies market were not adverse affected by the COVID-19 pandemic compared to stock, currency and commodities markets. In late 2020 and into 2021, Bitcoin and Ethereum illustrated an increase in price and reached new all-time high, which as a result of increased adoption and retail participation. However, in November 2021 through 2022, their prices decreased due to increase in interest rates, macroeconomic uncertainties and reduced liquidity in the financial markets. The asset prices were transformed into log-returns to ensure stationarity and prepare the data for volatility modelling. Figure 5 presents the log-returns plots for the twelve assets in the portfolio. The daily log-returns fluctuates around zero and are characterized by volatility clustering, where periods of high volatility follow each other (as in 2008 and early 2020) and periods of low volatility are followed by small changes.

Figure 4. Time plots of the Stock indices, Currency exchange rates and Commodities for the period starting from 2004 to 2024 and Cryptocurrencies 2018 to 2024.

Figure 5. Log-returns time plots of the Stock indices, Currency exchange rates and Commodities for the period starting from 2004 to 2024 and Cryptocurrencies 2018 to 2024.

Table 3 presents a descriptive summary of statistics, results from selected statistical tests, and the correlation matrix for the log-return series of a portfolio that includes stock indices, currency exchange rates, commodity prices, and cryptocurrency prices. This analysis covers the period from January 2004 to December 2024. The returns of the assets in the portfolio have a mean close to zero, with most assets showing positive means, except for natural gas. Stock indices yielded higher returns than other financial assets, but they also experienced greater volatility across all considered time intervals. The skewness for all financial assets was negative, except for the euro (EUR), which was right-skewed. This indicates that the left tail of the distribution for most assets is longer than the right tail, reflecting asymmetry. All asset returns exhibit positive excess kurtosis, indicating heavy tails in the log-returns distribution. The Jarque-Bera test confirms that the log-returns do not follow a normal distribution as assumed. Further validation from the Mardia and Henze-Zirkler tests also identifies non-normality in the return series. With all p-values falling below 0.05, we reject the null hypothesis that the returns conform to the multivariate normality assumption, revealing significant skewness, fat tails, and deviations from the multivariate normal distribution. The Augmented Dickey-Fuller test established that all return series are stationary. Results from the Ljung-Box test at lags 10 and 20 indicated the presence of serial autocorrelation in the log-returns of the S&P 500, NASDAQ, FTSE, JPY, GBP, gold, and natural gas, as their p-values were below the 5% significance level. In contrast, other assets did not show significant autocorrelation, suggesting the need for time-series models that account for serial correlation. Finally, the correlation matrix reveals that the assets in the portfolio exhibit both positive and negative correlations, with the S&P 500 and NASDAQ showing the highest positive correlation.

Table 3. Descriptive statistics, test statistics, and correlation matrix for portfolio log-returns of commodities, exchange rates, and cryptocurrencies.

Statistics	SP500		NASDAQ		DAX		FTSE		EUR			JPY		GBP		Gold		Natural gas		Silver
Descriptive Statistics for the log-returns (2004-2017)
No. of observations	3020		3020		3020		3020		3020			3020		3020		3020		3020		3020
Mean	0.0003		0.0004		0.0002		0.0001		0.0001			0.0001		0.0000		0.0001		−0.0002		0.0001
Std Dev	0.0107		0.0120		0.0123		0.0104		0.0060			0.0064		0.0059		0.0115		0.0299		0.0206
Min	−0.0947		−0.0923		−0.0727		−0.0925		−0.0280			−0.0378		−0.0841		−0.0981		−0.2115		−0.1949
Max	0.1025		0.1037		0.1067		0.0847		0.0319			0.0622		0.0304		0.0859		0.2677		0.1247
Skewness	−0.3744		−0.1705		−0.2187		−0.4063		0.0592			0.0032		−1.1024		−0.3999		0.8021		−1.0160
Kurtosis	10.8539		6.0660		5.2710		8.2872		2.2881			5.2820		14.9810		6.3402		7.0664		8.1238
Statistics	SP500	NASDAQ		DAX		FTSE		EUR		JPY	GBP		Gold		Natural gas		Silver		BTC		ETH
Descriptive Statistics for the log-returns (2018-2024)
No. of observations	1321	1321		1321		1321		1321		1321	1321		1321		1321		1321		1321		1321
Mean	0.0006	0.0007		0.0003		−0.0001		−0.0001		0.0000	0.0000		0.0005		−0.0008		0.0001		0.0008		0.0006
Std Dev	0.0121	0.0152		0.0120		0.0101		0.0046		0.0056	0.0057		0.0093		0.0359		0.0182		0.0362		0.0466
Min	−0.0999	−0.0973		−0.1305		−0.1151		−0.0206		−0.0386	−0.0369		−0.0474		−0.1160		−0.0085		−0.4973		−0.5896
Max	0.0897	0.0959		0.1041		0.0867		0.0212		0.0316	0.0310		0.0578		0.1873		0.0724		0.1774		0.2308
Skewness	−0.1798	−0.1671		−0.6864		−1.0553		0.0009		−0.3267	−0.1484		−0.2729		−0.0034		−0.4442		−1.1269		−1.0479
Kurtosis	9.7180	3.9579		16.1757		18.0339		1.5254		5.7428	3.9202		3.7802		4.5144		4.6055		16.5418		13.2428
Statistics	SP500	NASDAQ		DAX		FTSE		EUR		JPY	GBP		Gold		Natural gas		Silver		BTC		ETH
Descriptive Statistics for the log-returns (2004-2024)
No. of observations	4341	4341		4341		4341		4341		4341	4341		4341		4341		4341		2557		2557
Mean	0.0004	0.0005		0.0002		0.0000		0.0000		0.0001	0.0000		0.0002		−0.0003		0.0001		0.0008		0.0006
Std Dev	0.0111	0.0130		0.0122		0.0103		0.0056		0.0062	0.0058		0.0108		0.0319		0.0199		0.0362		0.0466
Min	−0.0999	−0.0973		−0.1305		−0.1151		−0.0279		−0.0386	−0.0841		−0.0981		−0.2133		−0.1949		−0.4973		−0.5896
Max	0.1025	0.1037		0.1068		0.0867		0.0372		0.0622	0.0310		0.0859		0.2677		0.1247		0.1774		0.2308
Skewness	−0.2958	−0.1647		−0.3520		−0.5901		0.0615		−0.0704	−0.8266		−0.3875		0.4485		−0.8858		−1.1269		−1.0479
Excess Kurtosis	10.5846	5.3965		8.3203		11.0099		2.4466		5.9384	11.8422		6.2306		6.1200		7.4952		16.5418		13.2428
Jarque Bera	21198	5446.3		13468		23537		1146.2		6757.4	28027		7258.8		7219.3		11055		28646		18612
p-value	0.0000	0.0000		0.0000		0.0000		0.0000		0.0000	0.0000		0.0000		0.0000		0.0000		0.0000		0.0000
ADF-test	−16.611	−16.437		−16.202		−16.232		−15.558		−15.997	−16.332		−17.22		−15.612		−16.839		−10.072		−10.708
p-value	<0.01	<0.01		<0.01		<0.01		<0.01		<0.01	<0.01		<0.01		<0.01		<0.01		<0.01		<0.01
Ljung-Box(10)	41.472	28.362		11.741		37.246		7.7132		42.471	20.815		21.432		28.315		18.548		14.676		15.689
p-values	0.0000	0.0016		0.3028		0.0000		0.6568		0.0000	0.0224		0.0183		0.0016		0.0464		0.1443		0.3089
Ljung-Box(20)	66.457	50.354		30.478		55.547		19.364		66.874	35.144		38.404		44.712		31.108		19.425		23.56
p-values	0.0000	0.0002		0.0625		0.0000		0.4983		0.0000	0.0194		0.0079		0.0012		0.0538		0.4944		0.2622
Correlation matrix
SP500	1.0000	0.9222		0.5762		0.5310		0.2117		0.2849	0.2663		0.0387		0.0666		0.1141		0.3377		0.3505
Nasdaq				0.5059		0.4465		0.1630		0.2442	0.2198		0.0225		0.0447		0.0915		0.3346		0.3453
DAX						0.8363		0.1031		0.2273	0.2063		0.0450		0.0609		0.1141		0.2553		0.2474
FTSE								0.1120		0.2175	0.1045		0.0816		0.0529		0.2117		0.2274		0.2240
EUR										−0.3114	0.6424		0.3358		0.0565		0.3191		0.1370		0.1469
JPY											−0.1961		−0.2683		0.0535		−0.1824		−0.0408		−0.0363
GBP													0.2674		0.0676		0.2809		0.1785		0.1804
Gold															0.0565		0.8021		0.1309		0.0971
Natural Gas																	0.0571		0.0162		0.0313
Silver																			0.1309		0.1306
BTC																					0.8572
ETH																					1.0000

3.2. Parameter Estimates for the Fitted MGARCH Models

To estimate the parameters of the ARMA-MGARCH models, the dataset was divided into two parts: 70% (3038 observations) for training and 30% (1303 observations) for testing. The training set was utilized to fit the models and estimate the parameters, while the testing set was reserved for forecasting risk measures and conducting back-testing of the results. In this section, we present the empirical results for the estimated mean and volatility components of the fitted ARMA-MGARCH models, along with the chosen appropriate innovation distributions. The empirical results of the model adequacy diagnostic checks for these fitted models are also presented.

3.2.1. Results for the Fitted ARMA-GARCH Model

The ARMA(p,q) model, where the orders of p and q range from zero to two, was applied to the log-return data to characterize the mean component. The optimal ARMA lag order for each asset was determined based on the smallest AIC and BIC values. Table 4 displays the AIC values for various ARMA orders, with bold values indicating the lowest AIC, suggesting the most suitable model. For the S&P 500, the optimal mean model was ARMA(2,1). The FTSE was best represented by ARMA(2,2), while both the JPY and Silver followed an ARMA(1,1) process. The optimal model for DAX and GBP log-returns was a constant model, specifically an ARMA(0,0) process with a non-zero mean, implying that the returns were essentially white noise. The MA(1) model was identified as the best fit for NASDAQ and EUR, whereas the MA(2) model was most suitable for Bitcoin (BTC). The AR(1) model was the best-fitting option for gold and Ethereum (ETH), and the AR(2) model was found to be the best for natural gas. Overall, there was no single model applicable to all assets, as each exhibited a different ARMA structure.

Table 4. AIC values for ARMA (p,q) specifications.

	ARMA(0,0)	ARMA(0,1)	ARMA(0,2)	ARMA(1,0)	ARMA(1,1)	ARMA(1,2)	ARMA(2,0)	ARMA(2,1)	ARMA(2,2)
SP500	−21111.79	−21114.23	−21113.75	−21114.03	−21112.94	−21115.59	−21113.94	−21115.61	−21113.65
NASDAQ	−20062.90	−20064.70	−20064.13	−20064.54	−20063.82	−20062.12	−20064.10	−20062.12	−20060.12
DAX	−20516.73	−20514.95	−20516.11	−20514.93	−20512.94	−20514.12	−20515.88	−20513.88	−20516.19
FTSE	−21695.36	−21694.21	−21693.29	−21694.18	−21692.20	−21691.27	−21693.12	−21691.13	−21696.87
EUR-USD	−25752.06	−25753.17	−25752.24	−25753.06	−25752.53	−25750.31	−25752.14	−25750.17	−25748.55
USD-JPY	−25157.78	−25162.15	−25160.17	−25162.16	−25160.16	−25158.16	−25160.16	−25158.16	−25157.43
GBP-USD	−25483.57	−25483.23	−25481.39	−25483.21	−25481.22	−25479.40	−25481.43	−25479.42	−25481.94
Gold	−21404.17	−21404.64	−21402.90	−21404.68	−21402.66	−21400.90	−21402.87	−21400.87	−21399.18
Natural gas	−14004.70	−14018.11	−14017.73	−14018.73	−14017.88	−14014.60	−14017.89	−14015.86	−14013.90
Silver	−17215.84	−17215.28	−17214.72	−17215.34	−17219.91	−17212.74	−17214.86	−17212.84	−17216.27
BTC	−4813.89	−4819.77	−4820.70	−4820.43	−4819.61	−4818.72	−4820.36	−4818.48	−4817.96
ETH	−4171.09	−4178.70	−4176.90	−4178.85	−4176.86	−4174.89	−4176.86	−4174.87	−4172.96

First, a standard GARCH(1,1) model was used for all financial assets to determine the most appropriate error distribution. The innovation distributions considered in this study included Student-t distribution, normal distribution, generalized error distribution, and skewed Student-t distribution. For each return series, a GARCH(1,1) model was fitted using each of these four distributions as the assumed error term. Table 5 presents the AIC and BIC values of the fitted GARCH(1,1) model for the respective innovation distributions. The model with the lowest AIC and BIC values indicates the most suitable innovation distribution for the given series. With the exception of natural gas and cryptocurrencies, the generalized error distribution was identified as the most appropriate error distribution for all other financial assets in the portfolio. The skewed Student-t distribution was selected for natural gas, while the Student-t distribution was used for cryptocurrencies. This can be attributed to the fact that most financial assets in the markets exhibit similar trends and patterns of volatility movements. In contrast, natural gas and cryptocurrency markets are influenced more by their specific market dynamics.

Table 5. AIC and BIC values of the fitted ARMA-GARCH(1,1) model for the various error distributions.

	Model	NORM		STD		SSTD		GED
	ARMA-GARCH(1,1)	AIC	BIC	AIC	BIC	AIC	BIC	AIC	BIC
SP500	ARMA(1,2)	−6.5645	−6.5520	−6.6604	−6.6460	−6.6648	−6.6487	−6.6728	−6.6585
NASDAQ	ARMA(0,1)	−6.1122	−6.1032	−6.1869	−6.1761	−6.1915	−6.1790	−6.2014	−6.1907
DAX	ARMA(0,0)	−6.2138	−6.2067	−6.2793	−6.2704	−6.2829	−6.2722	−6.3015	−6.2925
FTSE	ARMA(2,2)	−6.3163	−6.3038	−6.7134	−6.6973	−6.7178	−6.6999	−6.7207	−6.7046
EUR-USD	ARMA(1,0)	−7.6441	−7.6352	−7.6720	−7.6613	−7.6715	−7.6590	−7.6843	−7.6736
USD-JPY	ARMA(1,0)	−7.4682	−7.4593	−7.5594	−7.5487	−7.5601	−7.5476	−7.5714	−7.5607
GBP-USD	ARMA(0,0)	−7.5604	−7.5532	−7.6195	−7.6107	−7.6190	−7.6083	−7.6237	−7.6147
Gold	ARMA(1,0)	−6.3537	−6.3446	−6.4613	−6.4505	−6.4609	−6.4484	−6.4781	−6.4674
Natural gas	ARMA(2,0)	−4.2269	−4.2161	−4.3237	−4.3111	−4.3239	−4.3096	−4.3174	−4.3049
Silver	ARMA(1,1)	−5.1716	−5.1609	−5.3123	−5.2998	−5.3131	−5.2988	−5.3356	−5.3231
BTC	ARMA(0,2)	−3.6711	−3.6475	−3.9532	−3.9258	−3.9517	−3.9203	−3.9412	−3.9137
ETH	ARMA(1,0)	−3.2125	−3.1929	−3.4326	−3.4090	−3.4311	−3.4037	−3.4283	−3.4048

The volatility component of the ARMA-GARCH model was analyzed by fitting several types of GARCH models, including the Standard GARCH(1,1), EGARCH(1,1), GJR-GARCH(1,1), FGARCH(1,1), IGARCH(1,1), and CSGARCH(1,1). Table 6 displays the AIC and BIC values for each fitted GARCH model corresponding to the return series. The model with the lowest AIC and BIC values was deemed the most suitable and is highlighted in bold. The EGARCH(1,1) model recorded the lowest AIC and BIC values for seven out of the twelve assets in the portfolio, making it the optimal choice. Following to the minimal difference between the AIC and BIC values of the various GARCH model specifications fitted to the portfolio returns data, the ARMA-EGARCH model was chosen as the optimal GARCH specification for all return series, along with their respective innovation distributions. In the following sections of the study, the ARMA-ECARCH(1,1) model was employed to model within the multivariate GARCH framework.

Table 6. AIC and BIC values for fitted ARMA-GARCH(1,1) model specifications.

SGARCH

GJR-GARCH

EGARCH

IGARCH

CSGARCH

FGARCH

AIC

BIC

AIC

BIC

AIC

BIC

AIC

BIC

AIC

BIC

AIC

BIC

SP500

−6.6728

−6.6585

−6.6807

−6.6646

−6.6814

−6.6653

−6.6714

−6.6589

−6.6746

−6.6567

−6.6847

−6.6686

NASDAQ

−6.2013

−6.1907

−6.2092

−6.1967

−6.2122

−6.1997

−6.1999

−6.1909

−6.2022

−6.1879

−6.2151

−6.2025

DAX

−6.3015

−6.2925

−6.3144

−6.3037

−6.3166

−6.3059

−6.2999

−6.2928

−6.3029

−6.2904

−6.3180

−6.3073

FTSE

−6.7207

−6.7046

6.7314

−6.7135

−6.7324

−6.7146

6.7174

−6.7031

−6.7240

−6.7044

−6.7339

−6.7160

EUR-USD

−7.6843

−7.6736

−7.6852

−7.6727

−7.6828

−7.6702

−7.6834

−7.6745

−7.6836

−7.6693

−7.6803

−7.6677

USD-JPY

−7.5714

−7.5607

−7.5721

−7.5596

−7.5732

−7.5607

−7.5694

−7.5604

−7.5704

7.5561

−7.5731

−7.5606

GBP-USD

−7.6237

−7.6147

7.6231

−7.6124

−7.6244

−7.6137

−7.6208

−7.6137

−7.6233

−7.6108

−7.6222

−7.6115

Gold

−6.4781

−6.4674

−6.4778

−6.4652

−6.4834

−6.4709

−6.4765

−6.4675

−6.4775

−6.4632

−6.4821

−6.4696

Natural gas

−4.3174

−4.3049

−4.3183

−4.3040

−4.3242

−4.3099

−4.3150

−4.3043

−4.3168

−4.3007

−4.3225

−4.3082

Silver

−5.3356

−5.3231

−5.3354

−5.3211

−5.3385

−5.3242

−5.3346

−5.3239

−5.3357

−5.3196

−5.3374

−5.3231

Bitcoin

−3.9517

−3.9203

−3.9543

−3.9150

−3.9680

−3.9287

−3.9341

−3.9066

−3.9382

−3.8989

−3.9384

−3.8952

Ethereum

−3.4306

−3.4032

−3.4360

−3.4007

−3.4507

−3.4154

−3.4251

−3.4016

−3.4324

−3.3971

−3.4346

−3.3954

Table 7 presents the estimated parameters and residual diagnostic tests results for the ARMA-EGARCH(1,1) model. Most of the parameters in the model are statistically significant at the 5% level, with the exception of the ARCH component and the mean for certain assets. The persistence of the conditional volatility parameter was found to be positive and notably high (above 0.9) for all assets, indicating that volatility shocks have a lasting impact on future volatility. Gold and Silver demonstrate the highest persistence values, suggesting that their volatility is more enduring compared to other financial assets in the portfolio. Additionally, an asymmetric effect is observed, as the ARCH effect parameter estimate is significant for most assets, except for Gold and Silver. This indicates that negative shocks tend to have a greater impact on volatility than positive shocks. To validate the adequacy of the fitted ARMA-EGARCH(1,1) model, the Ljung-Box Q (LBQ) and LM-ARCH tests were applied to the standardized residuals of the ten assets in the portfolio. These diagnostic tests evaluate whether the model has adequately addressed serial correlation and conditional heteroskedasticity in the return series. For most of the portfolio assets, the p-values from both tests were above the 5% significance level, indicating no significant autocorrelation or remaining ARCH effects in the residuals. This suggests that the ARMA-EGARCH(1,1) model has effectively captured the dynamics of the data. The only exception was the DAX index, which exhibited evidence of remaining serial correlation. As a result, the ARMA-EGARCH(1,1) model with a generalized error distribution was chosen as the most appropriate univariate model for use in the multivariate GARCH framework. This choice is supported by its superior performance, indicated by lower AIC and BIC values across most assets, satisfactory residual diagnostics for the majority of return series, and its ability to capture asymmetries in volatility, which are characteristic of financial returns.

Table 7. The parameter estimates and residual diagnostic tests results of the fitted optimal ARMA-EGARCH(1,1) with generalized error distribution.

Parameter	SP500	NASDAQ	DAX	FTSE	EUR	JPY	GBP	Gold	Natural gas	Silver
Mean (μ)	0.0005	0.0006	0.0005	0.0002	0.0000	0.0000	0.0000	0.0002	−0.0002	0.0000
	(0.0001, 0.0001)	(0.0001, 0.0000)	(0.0002, 0.0000)	(0.0001, 0.1998)	(0.0000, 0.7680)	(0.0000, 0.8099)	(0.0000, 0.7011)	(0.0000, 0.0000)	(0.0002, 0.0802)	(0.0000, 0.9998)
ARCH term (α0)	−0.3389	−0.3917	−0.5371	−0.3644	−0.0762	−0.3368	−0.1493	−0.0882	−0.1817	−0.1094
	(0.1044, 0.0159)	(0.0254, 0.0000)	(0.0224, 0.0000)	(0.0571, 0.0000)	(0.0286, 0.0053)	(0.0099, 0.0000)	(0.0031, 0.0000)	(0.0039, 0.0000)	(0.0138, 0.0000)	(0.0039, 0.0000)
ARCH effect (α1)	−0.1103	−0.1149	−0.3759	−0.1126	−0.0107	−0.0299	−0.0161	0.0084	−0.0127	0.0099
	(0.0388, 0.0231)	(0.0171, 0.0000)	(0.0158, 0.0000)	(0.0172, 0.0000)	(0.0250, 0.6605)	(0.0135, 0.0267)	(0.0120, 0.1800)	(0.0090, 0.3538)	(0.0145, 0.3803)	(0.0106, 0.3490)
GARCH effect (β)	0.9642	0.9565	0.9588	0.9619	0.9927	0.9673	0.9856	0.9904	0.9744	0.9864
	(0.0112, 0.0000)	(0.0028, 0.0000)	(0.0024, 0.0000)	(0.0059, 0.0000)	(0.0031, 0.0000)	(0.0010, 0.0000)	(0.0002, 0.0000)	(0.0004, 0.0000)	(0.0018, 0.0000)	(0.0007, 0.0000)
Leverage effect (γ)	0.2126	0.2065	0.1965	0.2002	0.0945	0.1742	0.1060	0.0908	0.1747	0.1229
	(0.0595, 0.0054)	(0.0220, 0.0000)	(0.0217, 0.0000)	(0.0225, 0.0000)	(0.1436, 0.4926)	(0.0214, 0.0000)	(0.0216, 0.0000)	(0.0112, 0.0000)	(0.0205, 0.0000)	(0.0230, 0.0000)
Shape	1.0587	1.1048	1.1165	1.1968	1.2523	1.0546	1.2148	1.0133	1.1601	0.9017
	(0.0422, 0.0000)	(0.0399, 0.0000)	(0.0435, 0.0000)	(0.0496, 0.0000)	(0.1529, 0.0000)	(0.0473, 0.0000)	(0.0488, 0.0000)	(0.0368, 0.0000)	(0.0543, 0.0000)	(0.0386, 0.0000)
Residual diagnostic tests results p-values
LBQ(10)	0.8322	0.7255	0.1329	0.7542	0.7992	0.7566	0.9254	0.9481	0.2791	0.6241
LBQ(20)	0.7937	0.6744	0.2423	0.3137	0.2313	0.8850	0.9363	0.9557	0.0930	0.8949
ARCH-LM(10)	0.1112	0.0237	0.0079	0.1340	0.4559	0.9624	0.5781	0.3787	0.3184	0.0613
ARCH-LM(20)	0.3427	0.1876	0.0057	0.0670	0.7222	0.9255	0.9638	0.3179	0.7240	0.2008

The values under the parenthesis are the standard errors and p-values, respectively.

3.2.2. Parameter Estimates Results of the Fitted CCC(1,1)-EGARCH Model

The Constant Conditional Correlation (CCC) EGARCH model was fitted using a two-step procedure. In the first step, each asset in the portfolio was modeled using the ARMA-EGARCH(1,1) specification with a generalized error distribution, which was identified as the optimal univariate model. In the second step, the residuals obtained from the ARMA-EGARCH(1,1) model were used to estimate the parameters of the CCC model. This process resulted in the computation of a constant conditional correlation matrix that reflects the average pairwise correlations across the entire study period. Table 8 displays the estimated constant conditional correlations among the assets for the period of study. The equity indices showed strong positive correlations, particularly between the S&P 500 and NASDAQ, primarily due to their shared exposure to U.S. market conditions. The correlations among currency pairs varied, with the EUR and GBP exhibiting the highest correlation within the currency markets, which can be attributed to the economic ties between the Eurozone and the United Kingdom (UK). Conversely, the EUR and JPY displayed the highest negative correlation among the assets in the portfolio, indicating that as the price of EUR-USD increased, the price of USD-JPY decreased. Natural gas demonstrated the weakest correlations with all the financial assets in the portfolio, as its prices are mainly influenced by demand, supply, and weather factors, which are unrelated to the stock and currency markets. Overall, most of the assets were positively correlated with one another, suggesting that they moved in tandem during the sample period. The correlation matrix is based on the assumption that the relationships between asset returns remain constant over time. However, since correlations are expected to vary, the next section will explore the interconnections between assets from a dynamic correlation perspective.

Table 8. Constant correlation matrix estimated from the CCC(1,1)-EGARCH model with the generalized error distribution.

	SP500	NASDAQ	DAX	FTSE	EUR	JPY	GBP	Gold	Natural gas	Silver
SP500	1.0000	0.9226	0.5871	0.5377	0.2101	0.2861	0.2709	0.0296	0.0651	0.1160
NASDAQ			0.5223	0.4477	0.1570	0.2499	0.2229	0.0130	0.0453	0.0898
DAX				0.8307	0.0948	0.2317	0.1981	0.0298	0.0563	0.1707
FTSE					0.1091	0.2106	0.0999	0.0676	0.0508	0.2215
EUR-USD						-0.3073	0.6413	0.3255	0.0505	0.3163
USD-JPY							-0.1914	-0.2748	0.0647	-0.1849
GBP-USD								0.2522	0.0689	0.2686
Gold									0.0555	0.7980
Natural Gas										0.0577
Silver										1.0000

The volatilities of all twelve assets were estimated using the CCC(1,1)-EGARCH model, and the results are illustrated in Figure 6, which shows the time-varying conditional standard deviations for each asset. Natural Gas exhibits the highest conditional volatility, reflecting high market sensitivity, while the EUR is the least volatile, reflecting relative stability among the given assets. The equity markets show volatility spikes during the 2008 GFC and COVID-19 pandemic, while the currency markets show less pronounced spikes compared to the equities.

Figure 6. Plot of conditional volatility estimated using the CCC(1,1)-EGARCH model.

3.2.3. Parameter Estimates Results of the Fitted DCC(1,1)-EGARCH Model

The Dynamic Conditional Correlation (DCC) EGARCH model was used to capture the time-varying interdependencies among the financial assets in the portfolio. Initially, the return series for each asset was modeled using an ARMA-EGARCH(1,1) specification with a generalized error distribution. The residuals obtained from this model were then utilized in the DCC model. To identify the most suitable multivariate volatility model, we estimated the Standard, Asymmetric, and Flexible DCC models with both multivariate normal and Student’s t error distributions. Table 9 presents the log-likelihood, AIC and BIC values for the different DCC models considered, along with the estimated parameters of the fitted DCC(1,1)-EGARCH model with Student-t errors. The DCC(1,1)-MVT-EGARCH model was identified as the best fit, as it produced the highest log-likelihood and the lowest AIC and BIC values. In the selected model, the DCC parameters ($a$ and $b$ ) govern the short-term response of correlations to shocks and their persistence over time. The sum of $a+b$ being close to one indicates strong persistence of correlations, consistent with slow mean reversion. Additionally, both parameters are statistically significant at the 5% level, confirming the presence of strong time-varying correlations among the asset returns.

Table 9. Model selection results for DCC(1,1)-EGARCH specifications based on AIC and BIC, along with estimated coefficients, standard errors, and p-values of the selected DCC(1,1)-MVT-EGARCH model for the ten portfolio assets.

DCC-EGARCH models	Log-Likelihood	AIC	BIC
DCC-MVN	105868	−69.619	−69.387
DCC-MVT	108016.5	−71.033	−70.799
Asymmetric DCC-MVN	105871.9	−69.621	−69.387
Asymmetric DCC-MVT	108016.5	−71.032	−70.796
Flexible DCC-MVN	106073.6	−69.758	−69.538
DCC-MVT-EGARCH model parameter estimates
Parameter	Estimate	Standard Error	p−value
a	0.0121	0.0016	0.0000
b	0.9779	0.0042	0.0000
Shape	4.0000	0.2505	0.0000

Figure 7 illustrates the time-varying conditional correlation from the Dynamic Conditional Correlation (DCC) model between the S&P 500 and various other assets. Generally, the correlation among stock indices, commodities, and currency exchange rates in the portfolio fluctuates over time, but tends to follow similar patterns. For example, during the 2008 Global Financial Crisis (GFC), the correlation between the S&P 500 and assets such as NASDAQ, DAX, FTSE, JPY, and natural gas appears to have increased compared to the period from 2005 to 2007. This indicates that during the 2008 GFC, these assets experienced a greater degree of co-movement. In contrast, the correlations between the S&P 500 and the euro, British pound, gold, and silver were lower during the 2008 GFC. Since both the Constant Conditional Correlation (CCC) and DCC models may not fully capture non-linear relationships and tail dependencies, a more robust approach, such as the vine copula, is necessary.

Figure 7. Time-varying conditional correlation between S&P 500 and NASDAQ, DAX, FTSE, EUR, JPY, GBP, Gold, Natural gas, and Silver.

3.2.4. Results of the Diagnostic Checks of the Fitted MGARCH Models

To evaluate the adequacy of the fitted Constant Conditional Correlation (CCC) and Dynamic Conditional Correlation (DCC) models, several tests were conducted. These included the Ljung-Box test for serial correlation, the ARCH-LM test to check for remaining ARCH effects on squared standardized residuals, and Q-Q plots for a visual assessment of the distributional assumptions of the residuals. Table 10 displays the p-values for the Ljung-Box tests of the standardized residuals (Q(20)), the squared residuals (Q²(20)), and the ARCH-LM(20) tests under both the CCC and DCC models. With the exception of the DAX in the DCC model, all other assets showed no signs of ARCH errors, as the p-values from the ARCH-LM tests were above the 5% significance level. The Ljung-Box (Q²) test indicated that the residuals for the portfolio assets, apart from DAX, do not exhibit serial correlation. Failing the ARCH-LM and Ljung-Box tests suggests that the DCC model does not adequately capture the dynamic behaviour of the DAX series. This inadequacy could lead to inaccurate forecasts of volatility and correlation, potentially resulting in the underestimation or overestimation of the DAX’s contribution to overall portfolio risk and return, especially during periods of high volatility or significant autocorrelation in the DAX.

Table 10. ARCH-LM and Ljung-Box test p-values for standardized residuals and squared residuals of portfolio asset log returns.

	SP500	NASDAQ	DAX	FTSE	EUR	JPY	GBP	Gold	Natural gas	Silver
Q(20)
CCC	0.7940	0.7530	0.2241	0.3282	0.2221	0.9048	0.9374	0.9715	0.0828	0.8878
DCC	0.7862	0.6744	0.2313	0.1564	0.2314	0.8850	0.9359	0.9550	0.0889	0.8953
Q2(20)
CCC	0.4556	0.4133	0.0429	0.1292	0.6664	0.8752	0.8013	0.3359	0.1780	0.7269
DCC	0.3071	0.1555	0.0035	0.0783	0.6398	0.9273	0.9676	0.2463	0.7510	0.2166
ARCH-LM
CCC	0.4801	0.4474	0.0702	0.1036	0.7516	0.8744	0.7841	0.3293	0.1656	0.7094
DCC	0.3414	0.1876	0.0057	0.0618	0.7229	0.9254	0.9633	0.3179	0.7286	0.2006

3.3. Dependence Modelling Using Vine Copulas

This section presents the empirical findings from the dependency modeling procedure that utilizes various vine copula structures, including R-Vine, C-Vine, D-Vine, and S-Vine copula specifications. Each of these structures provides a distinct method for capturing the complex dependencies among financial assets.

First, the standardized residuals obtained from the MGARCH models are transformed into a uniform distribution over the interval $\left[0,1\right]$ using the probability integral transform, which employs the cumulative distribution function of the generalized error distribution. This transformation ensures that the marginal distributions of the return series are uniform, a crucial requirement for copula modeLling. Next, the construction of the vine copula is carried out using a sequential selection method proposed by Gruber (2015). In this method, trees are built incrementally, selecting pairs with the strongest dependencies first. This approach guarantees that the initial trees, which have a significant impact on model fit, are accurately specified. The Kendall’s tau correlation coefficient is then used to measure these dependencies, which is beneficial because it operates independently of the underlying distribution, making it particularly useful for combining different copula families. The empirical Kendall’s tau correlation for each pair of financial assets is calculated, and the spanning tree that maximizes the sum of the absolute values of the empirical Kendall’s tau is selected. After determining the optimal tree structure for each vine copula, the subsequent step is to select the most appropriate bivariate copula families and estimate their parameters. Finally, goodness-of-fit tests such as the Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC) are employed to identify the best-fitting models.

Table 11 shows the pairwise Kendall’s tau correlation coefficients of the standardized residuals of various financial assets, reflecting different magnitudes and directions of their pairwise dependencies. Similar to the correlation matrix derived from the Pearson’s and CCC-EGARCH models, the S&P 500 and NASDAQ exhibit the strongest pairwise correlations. Additionally, when summing the Kendall’s tau values for the S&P 500 row, it yields the highest overall value.

Table 11. Estimated Kendal’s tau values for the transformed standardized residuals of stock indices, exchange rates, commodities, and cryptocurrencies.

	SP500	NASDAQ	DAX	FTSE	EUR	JPY	GBP	Gold	Natural gas	Silver	BTC	ETH
SP500	1.0000	0.7566	0.3478	0.2692	0.1537	0.0753	0.2054	0.0543	0.0556	0.1214	0.1794	0.1879
NASDAQ			0.2991	0.1978	0.1310	0.0564	0.1809	0.0432	0.0419	0.1096	0.1767	0.1803
DAX				0.5316	0.0724	0.0603	0.1380	0.0263	0.0203	0.1100	0.1174	0.1141
FTSE					0.0430	0.0832	0.0134	0.0276	0.0319	0.1104	0.0808	0.0853
EUR						−0.2911	0.5088	0.2494	−0.0099	0.2563	0.1253	0.1286
JPY							−0.2389	−0.2783	0.0600	−0.1982	−0.0358	−0.0302
GBP								0.2252	0.0211	0.2365	0.1342	0.1416
Gold									−0.0145	0.5777	0.0796	0.0800
Natural Gas										−0.0090	0.0066	0.0172
Silver											0.1067	0.1052
BTC												0.6522
ETH												1.0000

3.3.1. Estimates for Regular Vine Copulas

The estimation of the R-Vine copula follows a two-step procedure to ensure an optimal fit to the data. In the first stage, the most suitable pair-copula families are selected based on the Akaike Information Criterion and the Bayesian Information Criterion. These criteria help identify the copula structures that best capture the dependence patterns while balancing model complexity and goodness of fit. The bivariate copulas considered in this study include: Gaussian (G), Student’s-t, Gumbel (G), Clayton (C), Frank (F), and Joe(J). Once the optimal pair-copula families are chosen, the second stage involves estimating their parameters using the maximum likelihood estimation method. This step ensures that the dependence structure is accurately quantified, allowing for a precise characterisation of relationships between financial assets.

The first and second trees in this procedure are illustrated in Figure 8. The tree structures demonstrate both negative and positive dependencies. Each tree consists of several nodes and edges, where each node represents an asset or a group of assets from various financial markets. The edges in the graph indicate the dependencies between pairs of assets, with labels specifying the chosen bivariate copula family and the accompanying values representing the estimated Kendal’s tau correlation between the assets. In the pair of S&P 500 and NASDAQ, the selected bivariate copula is Student-t, with a Kendall’s tau of 0.75. In the first tree, the Student-t copula is predominant, while the second tree employs a mix of different Gaussian and Archimedean copula families.

Figure 8. First and second tree of the R-vine structure.

3.3.2. Estimates for C-Vine and D-Vine Copula Specification

Using the same selection and estimation procedures as described for the R-vine in Section 3.3.1, the C-vine and D-vine copulas were also selected and estimated. In the C-vine structure, the root node for each level (or tree) is determined by calculating the sum of the absolute values of Kendall’s tau for each variable. The variable with the highest total dependence is then selected as the root node. This approach ensures that the chosen node has the strongest overall dependence on the other variables in the dataset. In this study, the S&P 500 was identified as the root node for the first tree level, as it showed the highest cumulative dependence among all asset pairs. Consequently, all other financial assets in the portfolio are directly linked to the S&P 500, which serves as the foundation for subsequent dependency modelling.

The C-vine and D-vine copula specifications for the first and second trees are illustrated in Figure 9 and Figure 10, respectively. In the D-vine structure, each node in the initial tree is connected to a maximum of two other nodes. The S&P 500 is connected to NASDAQ using a Student-t copula with a Kendall’s tau of 0.75, and to DAX with a Student-t copula and a tau of 0.35. Similar to the first tree of the R-vine, most asset pairs in the first and second trees were modelled using the Student-t copula, which reflects symmetric tail dependence and heavy-tailed joint distributions. In the C-vine model, there are notable exceptions, such as the S&P 500-Silver pair, which is modelled using a Frank copula, and the S&P 500-Natural Gas pair, which is fitted with a Gaussian copula. Additionally, both the ETH-S&P 500 and BTC-S&P 500 pairs are fitted with a Survival Gumbel copula. In the D-vine model, deviations from the t-copula specification are also present, particularly in the FTSE-ETH and BTC-ETH pairs, both modelled with a Survival Gumbel copula, and in the BTC-Natural Gas pair, which is modelled using a Joe copula. However, the R-vine, D-vine, and C-vine models capture the relationships between assets at a specific point in time without accounting for variations in dependence over time. Therefore, an S-vine model is considered as an alternative.

Figure 9. First and second tree of the C-vine copula structure.

Figure 10. First and second tree of the D-vine copula structure.

3.3.3. Estimates of S-Vine Copula Specification

The Stationary vine copula captures both cross-sectional and temporal dependencies. At time $t$ , the dependencies of the financial assets form a cross-sectional structure while at time $t+1$ , the structure is maintained but a connection between time $t$ and $t+1$ is introduced. Figure 11 presents the dependence structure of the financial assets in the portfolio at two contiguous points in time. The first time point is presented by node S&P 500, while second S&P 500-1. For each cross-sectional time point the same R-vine copula structure was modelled. Additionally, this R-vine structure is similar to the one in Figure 8. The serial dependence was modelled with a bivariate copula (Clayton with Kendall’s tau = 0.07) between S&P 500 and DAX in the first tree, connecting the time $t$ and $t+1$ . Additionally, in Tree 1, stronger unconditional dependencies are captured, with several large positive $\tau$ values observed, such as 0.75 for S&P 500 and NASDAQ (2,1) with a t-copula, and 0.65 for BTC and ETH (11,12) with a Gumbel 180 copula. Most of the assets pair in the first tree were modelled using a Student-t copula.

3.3.4. Comparison of the Vine Copula Specifications

After fitting various vine copulas to the data, we compared their performance using log-likelihood and information criteria. Table 12 displays the log-likelihood, AIC, and BIC for the different vine copula specifications examined in this study for the periods 2018-2024 and 2004-2018. For the period 2018-2024, the CCC-SVine model exhibited the highest log-likelihood and the lowest AIC value among the vine copulas considered, indicating it provided the best fit for modelling the dependence between the assets in the portfolio. In contrast, for the period 2004-2018, the DCC-SVine model achieved the lowest AIC value and the highest log-likelihood value, making it the preferred model for that time-frame.

Figure 11. First and second trees of the S-vine copula structure.

Table 12. Log-likelihood, AIC, and BIC values for fitted R-vine, C-vine, D-vine, and S-vine models over two periods: 2018-2024 (including BTC and ETH) and 2004-2018 (excluding BTC and ETH).

	2018-2024			2004-2018
Model	Log-likelihood	AIC	BIC	Log-likelihood	AIC	BIC
CCC-RVine	4790.42	−9414.83	−8984.32	9189.72	−18213.45	−17713.88
DCC-RVine	4793.17	−9422.34	−8997.02	9203.16	−18240.32	−17740.74
CCC-CVine	4779.01	−9382.03	−8925.58	9154.91	−18139.81	−17628.20
DCC-CVine	4782.92	−9389.84	−8933.39	9169.16	−18168.33	−17656.72
CCC-DVine	4748.35	−9320.69	−8864.25	9172.36	−18170.72	−17647.07
DCC-DVine	4752.64	−9331.27	−8880.01	9184.65	−18193.30	−17663.63
CCC-SVine	5038.13	−9590.25	−8329.83	9447.19	−18520.38	−17394.83
DCC-SVine	5035.07	−9590.15	−8345.29	9466.17	−18556.34	−17424.78

3.4. Results of Portfolio Market Risk Forecasts Based on Value-at-Risk and Expected Shortfall

In this section, we forecast the one-day-ahead out-of-sample Value at Risk (VaR) and Expected Shortfall (ES) using the methodology outlined in Section 3.7. The out-of-sample period spans from February 7, 2018, to December 20, 2024, encompassing a total of 1303 observations and 5,000 Monte Carlo simulations. Figure 12 and Figure 13 illustrate the VaR and ES forecasts at both 95% and 99% confidence levels for an equally weighted portfolio. These forecasts are derived from various Constant Conditional Correlation (CCC) and Dynamic Conditional Correlation (DCC) Vine copula specifications and are compared against the actual portfolio losses.

In Figure 12, the actual portfolio losses fluctuate below the Value at Risk (VaR) and Expected Shortfall (ES) thresholds, indicating that the Constant Correlation Coefficients (CCC) Vine models provide a relatively conservative estimate of downside risk. During the COVID-19 period, however, there was a significant spike in volatility, resulting in higher losses that clearly exceeded both the VaR and ES forecasts. This violation suggests that the model underestimated the severity of the tail event during that extreme period, which is a common limitation of constant correlation models; they often fail to capture time-varying dependencies and volatility clustering during crises. Outside of this period of crisis, the model performs reasonably well, with actual losses remaining within the VaR and ES bounds. Notably, the VaR and ES forecasts from the Dynamic Conditional Correlation (DCC)-based Vine copula models are more extreme compared to those generated by the CCC-Vine copulas. This implies that DCC-based Vine copula models, as shown in Figure 13, tend to be more responsive to changes in market conditions and may significantly overestimate portfolio market risk at both the 5% and 1% quantile levels—especially during periods of large return fluctuations, such as early 2020 (the COVID-19 period).

Figure 12. VaR and ES forecasts plots for an equally-weighted portfolio at 95% and 99% confidence levels from the CCC-based Vine copula model.

Figure 13. VaR and ES forecasts plots for an equally-weighted portfolio at 95% and 99% confidence levels from the DCC-based Vine copula model.

3.5. Backtesting Results of Portfolio Market Risk Forecasts Based on Value-at-Risk and Expected Shortfall

The backtesting procedure was carried out to assess the accuracy and reliability of the forecasted Value at Risk (VaR) and Expected Shortfall (ES) generated by eight MGARCH-Vine copula models. This evaluation also involved a comparison with the benchmark Constant Conditional Correlation (CCC) and Dynamic Conditional Correlation (DCC) models using both Normal and Student-t innovations. The focus was on comparing the predicted portfolio risk measures against the actual realized portfolio returns during the out-of-sample period. Specifically, we examined whether the observed exceedances (i.e., returns falling below the forecasted VaR) occurred at the expected frequency, and whether the ES forecasts adequately represented the average loss.

The 95% and 99% confidence levels were considered and the expected exceedances for the out-of-sample period, with $T=1303$ observations, were observed to be 65 for VaR at the 95% level and 13 at the 99% level. To determine whether the models accurately estimated VaR at these confidence levels, we calculated the actual exceedances—the number of times actual losses exceeded the estimated VaR at both the 95% and 99% levels. We then evaluated these exceedances using Kupiec’s unconditional coverage test and Christoffersen’s conditional coverage test. Table 13 presents the expected and actual exceedances of the Value-at-Risk (VaR) forecasts at the 95% and 99% confidence levels for the CCC, DCC, and various Vine copula models examined in this study.

Table 13. Number of exceedances of the VaR_95% and VaR_99% estimates across CCC, DCC and Vine copula models.

Confidence level	95%
	CCC-GARCH	CCC-GARCH	DCC-GARCH	DCC-GARCH	CCC-RVine	DCC-RVine
Expected exceedance	(Normal)	(Student-t)	(Normal)	(Student-t)
65	81	49	71	42	75	25
	CCC-CVine	DCC-CVine	CCC-DVine	DCC-DVine	CCC-SVine	DCC-SVine
	75	27	76	26	77	26
Confidence level	99%
	CCC-GARCH	CCC-GARCH	DCC-GARCH	DCC-GARCH	CCC-RVine	DCC-RVine
Expected exceedance	(Normal)	(Student-t)	(Normal)	(Student-t)
13	32	9	19	6	9	1
	CCC-CVine	DCC-CVine	CCC-DVine	DCC-DVine	CCC-SVine	DCC-SVine
	8	1	8	1	9	1

At the 95% confidence level, the number of exceedances for the VaR at this level is slightly higher than the expected value for the CCC-Vine copula and the CCC and DCC models with a Normal error distribution. Conversely, it is lower than the expected value for the DCC-Vine copulas and the CCC and DCC models with Student-t innovations. The CCC-RVine, CCC-CVine, CCC-DVine, CCC-SVine, and DCC (Normal innovations) models demonstrated good performance, with exceedance counts ranging from 71 to 77. However, the DCC-Vine models continued to overestimate risk, reporting only 25 to 27 violations, which is less than half of the expected number. This substantial deviation suggests poor estimation of tail risk. Among all models, the DCC (Normal), CCC-RVine, and CCC-Cvine models were considered the best performers at the 95% level, based on their exceedance behaviour.

At the 99% confidence level, the actual exceedances were lower than expected for most models, with the exception of the CCC and DCC models that utilized Normal innovations. The CCC-RVine, CCC-SVine, and CCC model with Student-t innovations recorded 9 exceedances compared to an expected 13, making them the closest to the expected number of violations among all other models. This finding suggests a reasonable estimation of extreme losses, avoiding a conservative or overly aggressive approach. In contrast, the DCC-based Vine copula models (DCC-RVine, DCC-CVine, DCC-DVine, and DCC-SVine) significantly overestimated tail risk, with only one actual exceedance. This indicates that the DCC-Vine copula specifications failed to effectively capture the magnitude and frequency of extreme losses in the portfolio. Therefore, based on the number of exceedances at the 99% confidence level, the CCC-RVine, CCC-SVine, and CCC (Student-t) models emerge as the best-performing models, as they align most closely with the expected violation count.

Table 14 presents the Kupiec unconditional ($L{R}_{uc}$ ) and Christoffersen conditional coverage ($L{R}_{cc}$ ) test statistics, along with their corresponding p-values in parentheses, at 95% and 99% confidence levels for the various CCC, DCC, and Vine copula models under consideration. The models retained in Table 14 are those whose number of VaR exceedances were reasonably close to the expected values at both the 95% and 99% confidence levels, indicating that they satisfied the preliminary exceedance criterion. The Kupiec unconditional coverage test is used to evaluate the null hypothesis that the actual observed exceedances are equal to the expected number predicted by the VaR confidence level. The null hypothesis, denoted as $H_0$ , is rejected if the actual exceedances are significantly different from the expected number. Based on the p-values ($L{R}_{uc}$ ), the CCC (Student-t), DCC (Normal), CCC-RVine, CCC-CVine, CCC-DVine, and CCC-SVine models yield exceedances that were not significantly different from the expected number, indicating that the null hypothesis of correct unconditional coverage was not rejected at the 5% level of significance. The DCC with Normal innovations, CCC-RVine, and CCC-CVine were identified as the best fit models under the 95% confidence level, while the CCC with Student-t innovations, CCC-RVine, and CCC-SVine were viewed as the best under the 99% confidence level due to having the highest p-values.

Value-at-Risk (VaR) violations often cluster during periods of heightened market volatility, which can compromise the reliability of risk forecasts. To address this potential time-dependency in exceedances, the Christoffersen Conditional Coverage test was employed. This test assesses not only whether the number of observed violations aligns with the expected frequency but also whether these violations occur independently over time. The conditional coverage p-values were evaluated at the 5% significance level. The CCC-RVine, CCC-CVine, CCC-DVine, and CCC-SVine are the only models for which the null hypothesis was not rejected, indicating that these models exhibit independence in exceedances and adequate coverage. For both the unconditional and conditional coverage tests, the CCC-Based Vine copulas were models for which the null hypothesis could not be rejected. This implies that the CCC-Vine copula models are sufficient for forecasting one-day-ahead VaR. Notably, the CCC-RVine copula model shows the smallest deviation in the number of exceedances and the highest unconditional and conditional coverage $p$ -values, which suggest that during the specified backtesting period, CCC-RVine performed best in forecasting the portfolio VaR.

Table 14. Backtesting results of Value at Risk Conditional and unconditional Coverage tests p-values at the 95% and 99% confidence levels for CCC, DCC, and Vine copula models.

Confidence level	95%		99%
Model	$L{R}_{uc}$	$L{R}_{cc}$	$L{R}_{uc}$	$L{R}_{cc}$
CCC-GARCH (Student-t)	4.5923 (0.0321)	6.6791 (0.0355)	1.41204 (0.2347)	1.53733 (0.4636)
DCC-GARCH (Normal)	0.5379 (0.4633)	5.9915 (0.4321)	2.4207 (0.1197)	3.5977 (0.1654)
CCC-RVine	1.49796 (0.2210)	1.6132 (0.4464)	1.41204 (0.2347)	1.53733 (0.4636)
CCC-CVine	1.49796 (0.2210)	1.61321 (0.4464)	2.27458 (0.1315)	2.37350 (0.3052)
CCC-DVine	1.80962 (0.1786)	1.88758 (0.3891)	2.27458 (0.1315)	2.37350 (0.3052)
CCC-SVine	2.14923 (0.1426)	3.44842 (0.1783)	1.41204 (0.2347)	1.53733 (0.4636)

The backtesting of Expected Shortfall (ES) was conducted using two methods: exceedance residuals and the ES Regression method, to assess the accuracy of the forecasted ES values. Table 15 presents the p-values for both methods at 95% and 99% confidence levels. In the exceedance residual test, the null hypothesis that the mean of excess violations of VaR is equal to zero was not rejected for the CCC-RVine, CCC-CVine, CCC-DVine, CCC-SVine, and CCC and DCC with Student-t innovations at the 99% confidence level, as the p-values were greater than the 5% significance level. This suggests that these models successfully forecast the portfolio ES. The ES regression tests provide a stricter evaluation of the joint accuracy of VaR and ES forecasts. The asymptotic version of the test at the 99% level indicated that CCC-RVine, CCC-DVine, CCC-SVine, CCC (Student-t), and DCC (Normal) performed well, with p-values significantly above 0.05, while CCC-CVine models did not pass the test. The bootstrap version of the regression test, which is more robust in finite samples as it does not rely on asymptotic distributional assumptions, confirmed these findings, showing that CCC-RVine, CCC-DVine, CCC-SVine, CCC (Student-t), and DCC (Normal) maintained acceptable performance. In contrast, CCC-CVine produced p-values below 0.05, indicating marginal or failed performance. At the 95% confidence level, all models except for DCC with Normal innovations passed the Exceedance Residual tests (both two-sided and one-sided), with p-values well above 0.05. However, the ES Regression test (both asymptotic and bootstrap) again highlighted that only DCC with Normal error distribution, CCC-RVine, CCC-CVine, CCC-DVine, and CCC-SVine met the criteria with p-values above the 5% significance level.

Table 15. Backtesting results of Expected Shortfall p-values at the 95% and 99% confidence levels for CCC, DCC, and Vine copula models.

Confidence level	95%
Model	CCC-GARCH	DCC-GARCH	CCC-RVine	CCC-CVine	CCC-DVine	CCC-SVine
	(Student-t)	(Normal)
Excedance Residual - twosided simple	0.796	0.005	0.213	0.148	0.187	0.203
onesided simple	0.687	0.022	0.872	0.886	0.681	0.873
ES Regression - twosided asymptotic	0.0175	0.1621	0.5770	0.2787	0.1304	0.2110
twosided bootstrap	0.0610	0.2660	0.559	0.378	0.164	0.282
Confidence level	99%
Model	CCC-GARCH	DCC-GARCH	CCC-RVine	CCC-CVine	CCC-DVine	CCC-SVine
	(Student-t)	(Normal)
Excedance Residual - twosided simple	0.975	0.018	0.974	0.101	0.986	0.969
onesided simple	0.7222	0.01	0.736	0.101	0.769	0.734
ES Regression - twosided asymptotic	0.8734	0.7316	0.4184	0.0000	0.5020	0.2671
twosided bootstrap	0.77	0.312	0.277	0.0020	0.453	0.1310

Overall, the CCC specifications (CCC-RVine, CCC-DVine, and CCC-SVine) successfully passed VaR and ES backtests procedures at both the 99% and 95% confidence levels. In contrast, DCC-Vine copula models generally underperform, with consistently low p-values and frequent failure in the backtesting tests. This suggests that while time-varying dependence structures may offer theoretical flexibility, they introduce estimation challenges when combined with high-dimensional and flexible vine copulas, which undermines the accuracy of VaR and ES forecasts. These results are consistent with the findings of [37], in that vine copula with constant correlation outperforms the models with dynamic, time-varying correlations in forecasting VaR and ES. Therefore, among the models tested, the CCC-RVine copula model recorded the highest p-values across all VaR and ES backtesting procedures, underscoring its superior accuracy and robustness in out-of-sample risk forecasting. This suggests that the CCC-RVine model provides the most accurate and reliable forecasts of potential extreme losses in this study.

4. Conclusion

This study combines multivariate GARCH models with vine copulas to improve portfolio risk modelling for twelve financial assets selected from cryptocurrencies, stock indices, exchange rates, and commodities. The analysis uses daily data, with an in-sample period from January 2004 to January 2018 and an out-of-sample period from February 2018 to December 2024. The ARMA-EGARCH(1,1) model with generalized error distributions was identified as the best-fitting univariate specification for integration into the MGARCH framework. The Constant Conditional Correlation (CCC) and Dynamic Conditional Correlation (DCC) MGARCH models were employed to capture volatility clustering and conditional correlation, while remaining non-linear and tail dependencies were modelled using vine copula constructions. The multivariate CCC-GARCH model estimated a constant correlation over time, whereas the multivariate DCC-GARCH models captured time-varying correlation patterns, particularly during the 2008 Global Financial Crisis. During this period, the correlation between assets such as the S&P 500 and NASDAQ increased, while it decreased between Gold and the S&P 500. The residuals derived from both the multivariate CCC-GARCH and DCC-GARCH models were transformed into uniform variables through the probability integral transform and used to estimate the dependence structure via copula modelling. The R-vine, C-vine, D-vine, and S-vine copulas were utilized alongside bivariate copulas such as Gaussian, Student’s t, Gumbel, Clayton, Frank, and Joe. The CCC-SVine structure exhibited the lowest Akaike Information Criterion (AIC) value, indicating superior performance during the 2018-2024 period when the portfolio included Bitcoin and Ethereum, thereby capturing stronger and more flexible tail dependencies. The empirical results demonstrated that the CCC-Vine copula models outperformed both standard CCC and DCC models with normal and Student’s t innovations, as well as the DCC-Vine copula model at both 95% and 99% confidence levels. The CCC-RVine model yielded more accurate Value at Risk and Expected Shortfall forecasts, successfully passing all out-of-sample backtesting tests. In contrast, the DCC-Vine model consistently failed the backtesting tests, highlighting its inadequacy in this study. This suggests that the greater flexibility of vine copulas in capturing complex, non-linear dependencies may diminish the marginal benefit of modeling dynamic linear correlations through DCC, making the simpler CCC models more effective for out-of-sample forecasting when combined with vine copulas. This paper contributes to the literature by providing empirical evidence that integrating multivariate GARCH and vine copulas can significantly enhance risk assessment in portfolios comprising diverse assets. For practitioners, the findings underscore the practical value of utilizing flexible copula-based methods to account for non-linear and tail dependencies. For future researchers could consider extending the multivariate GARCH-copula framework by incorporating regime-switching structures to capture dynamic changes in tail dependencies experienced in financial time series. Additionally, applying this MGARCH-Vine copula approach to higher-frequency intra-day portfolio data, and calculating risk measures over longer horizons than one day could further explore its scalability and robustness in various market conditions. Future studies could also investigate applying MGARCH-Vine copula models for portfolio optimization, leveraging improved risk forecasts to design more robust investment strategies.