Risk management is an old subject for banking industry, regulators and academia. After the “Black Monday”, financial professionals with quantitative background worried about the firm-wide risk management. Thus, big banks started evaluating risks with quantitative tools. In the early 90s, most banks’ methods were already close to the formal concept of Value at Risk (VaR). In 1994, JP Morgan published its extensively developed quantitative methodology on evaluating the VaR and made a free access, the RiskMetrics, to the public. Thereafter, US Securities and Exchange Commission (SEC) required all big banks on Wall Street to disclose their quantitative risk evaluation of trading activities associated with derivatives. VaR evaluation became a more formal and common requirement after the release of Basel Accord II, which included detailed requirements on both VaR estimation and backtesting. The release of Basel III, together with the release of regional regulation frameworks such as Comprehensive Capital Analysis and Review (CCAR, the regulation framework released and used in USA), has made VaR a popular risk measure in the practice of risk evaluation.

This paper presents a hybrid method on the one-day-ahead prediction of high-dimensional portfolios VaR. In time series analysis, a popular way to validate the predictive power of a methodology is to use historical data sets sequentially with moving windows of a certain size. An analysis conducted with moving window is often referred as a rolling analysis. This paper also presents real data rolling analyses with sample portfolios and the corresponding backtesting results, utilizing publicly available historical data sets and different moving windows.

2.1. ARMA and multivariate GARCH models

The Autoregressive-Moving-Average (ARMA) model is a combination of the Autoregressive model and the Moving-Average model. The general form of the ARMA(p, q) model follows:

where X_t, a current value from a time series at time point t, can be expressed as the result of some general effect θ, the historical values of the series X₁, . . . , X_t–1 and the error terms (also referred to innovation terms or residual terms) ?₁, . . . , ?_t.

The Generalized Autoregressive Conditional Heteroskedasticity (GARCH) model developed by Bollerslev (1986), is the extension of the Autoregressive Conditional Heteroskedasticity (ARCH) model (Engle, 1982). The multivariate version of GARCH has many different specifications mainly to resolve the dependency problem. In this paper, Constant Conditional Correlation (CCC) specification (Bollerslev, 1990), its enhanced version, Dynamic Conditional Correlation (DCC) specification (Engle, 2002), and Generalized Orthogonal GARCH (GO-GARCH) model (Van derWeide, 2002) are used in modeling procedures.

The standard univariate GARCH (p, q) model has the form as:

where α’s and β’s are coefficients, and ?’s are the error terms from the ARCH model. GARCH models explicitly account for the effects from error patterns and the historical variances. The formula of univariate GARCH (1, 1) is given as:

The multivariate version of GARCH pushes the edge even further to account for the dependency relationship, or covariance, of the multiple volatility that is included in the multivariate time series model. Assume that we have a time-varying N by 1 vector x_t and that we also have a mean vector μ_t. With given ?_{t – 1} which includes all available information on the time horizon up to the time point t – 1, we can model the vector x_t:

Note that

where H_t is the conditional covariance matrix of the vector x_t and Ω_t is an i.i.d. random vector with E(Ω_t) = 0, and VaR(Ω_t) = I_N. The conditional variance matrix H_t can also be defined as (Bauwens et al., 2006):

For the bivariate case, the correlation between variables, or two time series, is presumed to be static over the time. The stochastic relationship does not change between two variables over time. Thus, a simple CCC-GARCH(1, 1) is expressed as:

where h_ii,t indicates the variance of the i^th entry at time t, h_{i j,t} indicates the covariance between i^th entry and j^th entry at time t, ρ_{i j} is the static correlation coefficient for i^th entry and j^th entry, ω’s, α’s, β’s are coefficients, and ?’s stand for error terms.

According to CCC, the covariance entirely depends on the prediction of the variance terms, once a fixed correlation estimation is provided with a certain historical time horizon. Therefore, CCC accounts for only a partial effect of the interaction between variables in the system. But a major advantage of this specification is the unrestricted applicability for large systems of time series (Franke et al., 2008).

As a successor of CCC, DCC specification is proposed, which bridges the gap of time varying interaction between variables (Engle, 2002). The idea is similar to the CCC specification. The DCC specification estimates the covariance and variance separately, but allows the existence of the time-varying correlations instead of assuming that the correlation remains the same. The correlation matrix R_t can be expressed as:

with

where $v_t={({?}_{1,t}/\sqrt{h_{11,t}},\dots ,{?}_{d,t}/\sqrt{h_{dd,t}})}^{?}$, Ω is a symmetric positive definite parameter matrix, γ and δ are scalar parameters, and Q_t is positive definite. Once the changing correlation matrix is found at a certain time point of the multivariate time series, the corresponding variance-covariance matrix at the time point can be defined by the same logic of CCC specification.

The generalized orthogonal GARCH or GO-GARCH (Van der Weide, 2002) is a special case of the early BEKK specification (named after Baba, Engle, Kraft, and Kroner), which is one of several ground-breaking specifications of multivariate GARCH models. With the assumption that the errors in the standard GARCH model are dependent on some uncorrelated component and can be modeled by that component through a linear mapping. The logic of GO-GARCH can be expressed as:

where Z is a linear mapping factor, and y_t is a unobservable factor. Notice that the Z is assumed to be invertible. If y_t follows a GARCH(1, 1), then

where $\mathrm{\Omega }={\sum }_{i=1}^d(1-{\alpha }_i-{\beta }_i)e_i{e}_i^{?}$, λ_i = Ze_i, ω_i = (Z⁻¹)^?e_i. Here, e_i is a vector that has only the i^th entry equal to 1, and the rest are all zero.

With the above assumption and definition, Z can have the following decomposition:

where R is a rotation matrix. The unconditional covariance matrix is then

This specification balances the computational cost and the complexity of covariances over time (Broda and Paolella, 2009).

To help understand the modeling process with the combination of the GARCH and ARMA models, let us take DCC-GARCH(1, 1)-ARMA(1, 1) as an example. With DCC-GARCH(1, 1), the variance and covariance of the multivariate time series can be modeled as:

where h_ii,t is the variance of i^th time series at time t, h_{i j,t} is the covariance of the i^th and j^th time series at time t, and ρ_{i j,t} is the changing correlation for the covariance. At the mean time, the multivariate mean of the time series can be modeled with the ARMA(1, 1) model. More details are included in Section 3.

2.2. Introduction to copulas

Modeling returns’ dependency of different assets is important, especially for portfolios related to a volatile capital market such as future market and derivatives market. The conventional techniques for modeling the dependence among different variables are multivariate distributions, usually the multivariate normal distribution. However, when modeling the dependency among returns from different assets in a portfolio, the dependency is restricted by a specific theoretical distribution assumption, and the multivariate distribution does not provide a flexible fitting. Furthermore, for higher dimensional portfolios, the joint distribution is difficult to estimate. In addition, if the multivariate normal distribution is used to fit the returns, the tail part is automatically assumed to be not fat and the returns are assumed to be symmetrically distributed. Unfortunately, these assumptions are not verified by the market. A popular alternative for dependency quantification is to utilize copulas, which provides different structures to simplify the estimation of a multivariate joint distribution and also allows different marginal distributions to appear in the structures. For random variables x_i’s, an informal explanation of copula technique can be symbolized as:

where C_θ(F(x₁), . . . , F(x_n)) is a copula expression of the multivariate joint distribution F(x₁, . . . , x_n), θ is a parameter of the copula structure, and F(x_i)’s are marginal distributions. For time-varying scenarios, copula allow models to capture the information of time-varying dependency.

Copula can be defined in different ways for different situations. Commonly used copulas include Gaussian copula, t copula, Gumbel copula, Clayton copula, Ali-Mikhail-Haq copula, and those extensions of Archimedean copulas in hierarchical structures and nested structures. For completeness, definitions of some popular copulas are listed below. The general definition of copula is:

3.1. Modeling methodology

Statistically, VaR_1–α, or sometimes VaR_α, which indicates the VaR at 1 – α confidence level, can be calculated using the lowest 100 × (1 – α) percentile of the returns or log returns. For example, the simplest nonparametric method to calculate VaR of an asset is the historical simulation method. Basically, it takes the value at the lowest 100 × α percentile of historical returns, or historical log returns, directly as VaR_α. Moreover, the log returns X_t, which is the return time series, can be modeled with the formula:

where μ_t is the mean of log returns by the time point t – 1 or the estimated log returns at time point t, σ_t is the square root of the associated variances predicted from the GARCH model at time point t, and ε_t is the error terms from the ARMA-GARCH model with distribution assumption for the included noises at time point t. Note that the distribution assumptions can be flexible, since we only care about the fit of the tails.

For example, if the multivariate GARCH(1, 1) and the ARMA(1, 1) are applied to model the returns, the variance-covariance matrix and the error vector on the right hand side of the above equation can be produced from the multivariate GARCH(1, 1) model with distribution assumption for noises included in the error terms. At the mean time, ARMA(1, 1) can produce the modeling results for mean return vector on the right hand side of the equation.

Inspired by the properties of ARMA models, GARCH models and the concept of copulas, the error terms from ARMA-GARCH models for multivariate cases can be modeled and simulated with copulas and marginal distributions, and a sequential modeling process can be implemented to model and predict VaR:

• Step 1: Model the log return time series with multivariate ARMA-GARCH.

• Step 2: Fit marginal distributions and a copula using the errors from the multivariate ARMA-GARCH.

• Step 3: Get the prediction of mean and variance from the multivariate ARMA-GARCH.

• Step 4: Simulate the error terms from the results of Step 2.

• Step 5: Acquire a set of predicted log return values with results from Step 3 and Step 4, and determine the value at 100 × α percentile of the log return values, which is the VaR at the corresponding time point.

To avoid any confusion, the error terms from GARCH, or ARMA-GARCH refer to the error terms of the mean process. GARCH models inherit the error terms from the setup of the ARCH model, and the error terms applied to model the variance in the ARCH model ultimately refer back to the error terms in the mean process, which is the ARMA model in our modeling process.

3.2. Copula fit

The copula fit of the error terms’ dependency structure is the key to the modeling process. We assume the error distribution to be multivariate normal; however, the choice of copula is flexible depending on the emphasis on the tail distribution of the shocks. The estimation of a copula based multivariate distribution involves both the estimation of the copula parameters θ and the estimation of margins F(x_i). Suppose δ_i is the parameter of F(x_i), the parameter vector is Ω = (δ₁, δ₂, . . . , δ_d, θ)^T . Then the log-likelihood function is:

where c(u₁, . . . , u_d) = ∂^dC_θ(u₁, . . . , u_d)/(∂u₁ . . . ∂u_d), and x_t = (x_1,t, x_2,t, . . . , x_d,t). In order to obtain the MLE, we firstly need to estimate the parameter vector Ω. By Full Maximum Likelihood (FML) Estimation, estimating Ω is equal to solve the equation:

However, it costs too much time in practice if the dimension is not moderate. Another popular estimating method is called the Inference for Margins (IFM), unlike the former one step FML method, this method estimates the parameter from the marginal distributions δ_i in the first step; second, it estimates the dependence parameter θ through the pseudo log likelihood function. The main advantage of IFM method is that it highly reduces the computational and notational complexity (Franke et al., 2008).

This work is not limited with the gaussian copula or t copula. As a visual demonstration of the goodness of fit, we show an example of a 3-dimensional 250-day dependency structure with returns of MSFT, AMZN, and AAPL in Figure 1. Note that the simulated multivariate distribution of the error terms reasonably covers the entire area of the realized errors.

3.3. Model validation

To examine the adequacy of the VaR measures produced by our copula-based models, we turn to a conventional “backtesting” validation method, which is used through our entire modeling process. The backtesting procedure evaluates the quality of the forecast of a risk model by comparing the actual results with those generated by the VaR models. For example, consider an event that the loss of a portfolio exceeds its reported VaR, VaR(α) at quantile α. Denoting the profit or loss of the portfolio between time t to t + 1 as r_t,t+1, the “hit” function I_t+1(α) can be defined as (Christoffersen, 1998):

Thereupon, the hit function sequence, e.g., (1, 0, 1, 0, 0, 0, 0, 1), tells the times that the loss of the portfolio has exceeded the labeled VaR(α) during the observed period.

Backtests can often be classified by if they examine an unconditional coverage property or independence property of a VaR measure. However, the test is called an Unconditional Coverage (UC) test if we are only interested in the times of the reported VaRs being violated. According to Christoffersen and Pelletier (2004), the interim between two VaR violations should exhibit no duration dependency, which is specified as the independence property. Therefore, the Conditional Coverage (CC) test is a joint test that checks both the violated times and the independence property.

To examine the stability of the modeling methodology and its accuracy, different example portfolios were established with assets from global financial markets including stocks, index futures, and commodities. A key assumption is that the parameter estimates should be constant over time when analyzing financial time series data using the Multivariate GARCH and Copula models. However, in a real financial market, considerable changes make it unreasonable to make this assumption. Therefore, we compute the parameter estimates over a rolling window of fixed sizes, e.g., 250, 500, through the samples. If the parameters are truly constant over the entire time horizon, then the estimates over the rolling windows should not be too different. If the parameters oscillate too much at some point during the period, then the rolling estimates should capture this volatility. Details of the example portfolios can be found in Table 1. We implemented mainly DCC-GARCH for 3-dimensional cases and GO-GARCH for 5-dimensional and 9-dimensional cases. Different copulas, simulation sizes and moving windows were tested. Details of the modeling specification and corresponding backtesting results can be found in Tables 2–4. Note that the key observation should be on the backtesting results, which include the p-values for both unconditional coverage and conditional coverage tests as well as the exceedance.

With our 3- and 5-dimensional experimental data sets, the modeling process provides a good and stable performance as indicated with two backtesting results, as well as the exceedance. In addition, it is important to pay attention to the balance of the results in the evaluation of such process. The exceedance should not be too far away from 95%, because capital reserve requirements for financial institutions will be evaluated according to the VaR models and VaR Stress Testing performance. If a model is claimed to be at a 95% confidence level but always delivers results close to 99%, the efficiency of capital allocation will be reduced due to the 4% difference. From this perspective, the 3- and 5-dimensional modeling examples also present a good performance. The exceedance of the modeling results are reasonably close to the 95% target level.

For 9-dimensional examples, especially those examples with skewed-t and the provided NAC copula structure, the modeling process demonstrates a very conservative projection on VaR. The modeling process highly depends on copula structures; therefore, it is clear that copula’s capability of explaining high-dimensional dependency decreases along with the increase of dimensionality. The modeling process tends to overestimate the risk in 9-dimensional cases. One potential way to improve this part is to compare the performance among different copulas in order to select the best copula for the corresponding dependency structure. Another simpler way would be to increase the simulation size.

5.1. Marginal distributions in multi-dimensional copulas

Copula methods are introduced to decompose the continuous joint distribution into individual marginal distributions so that the dependence structure can be easily studied. This raises concerns on the selection of marginal distribution and fitting of the distributions using the copula structures. Commonly used dispersion structures include: autoregressive of order 1, exchangeable, Toeplitz, and unstructured. Some researchers use one of the correlation matrices (usually autoregressive of order one and simple exchangeable Archimedean (Yan, 2007)) as a direct proxy for multivariate dependency. However, it is widely acknowledged that volatility parameters and variables are not always normally distributed. Many financial related variables often have fat tails and exhibit “tail dependence”. In our study, the skewed-t distribution fits the real log-return density function well most of the time. However, the skewed-t tends to underestimate peak areas.

According to the definition of Normal Mean-Variance Mixture Representation of Skewed-t Distribution, if we fix the degree of freedom, it is possible to fit financial data well with elliptical distributions. In practice, the normal mean-variance mixture formatted skewed-t distribution presents good flexibility to fit different types of portfolios. In addition, skewed-t distribution has non-zero tail dependence since the covariance matrix ∑ is positive semi-definite. The skewed t distribution also gives heavier tails at either upper or lower part of the margins. Therefore, we prefer skewed-t distribution to skewed normal as a margin to fit copulas.

5.2. GO-GARCH Model

A DCC-GARCH model is designed to overcome the drawbacks on the covariance estimation that early specifications such as BEKK and CCC models present (Van der Weide, 2002); however, the DCC specification is still not perfect due to the limitation on the possible failure of finding an invertible matrix. DCC assumes the existence of at least one orthogonal matrix that links the observed variables linearly to a set of components similar to the concepts of latent variables. Since the corresponding latent components are assumed to be independent, it is now possible to get dynamically time-varying conditional covariance for each of the multivariate time series systems presenting heteroskedasticity. Nevertheless, it is not guaranteed that a linkage or an invertible matrix will always exist. For example, if the diagonal elements of H, the covariance matrix, are not distinct, the transformed orthogonal matrix will no longer coincide with the linkage Z. Orthogonal GARCH (O-GARCH) and Generalized Orthogonal GARCH (GO-GARCH) were proposed successively in-between 1994 to 2002. GOGARCH allows Z to be given by any possible invertible matrix. Instead of using all of the historical information, the estimation of Z depends on the recent time point. Thus, GO-GARCH introduces a much easier and balanced computation method on dynamic conditional multivariate time series systems, especially high dimensional time series systems. For our exhibitions, we prefer GO-GARCH in high dimensional cases. It is important to mention that Z is chosen to be constant over time for the prediction of covariance. Since the direction of Z does not change during the process, this constant assumption holds.

Through the investigation on applicability of the proposed modeling methodology, we conclude that the proposed modeling process is a good alternative approach to study and predict VaR. With multivariate time series and copula, the proposed modeling methodology keeps a good balance between the complexity and flexibility. At the meantime, it helps reduce the computational cost. Furthermore, a methodology combining different techniques provides an opportunity to conduct the model selection over a system of models which can be easily constructed. The copula technique plays an important role in the description of the dependency relationship. More importantly, it enables us to directly model the VaR with a regard to the interactions among assets at a high dimensional level. During the study and experimental modeling investigation, we also noticed that the effect of the mean time series on VaR was weak when comparing to the variance effect. The fit of the marginal distributions is also not as important as the fit of GARCH and copula. Additionally, there was no universal solution. Therefore, when modeling with the methodology presented in this paper, we should use a greater amount of caution and patience to adjust different components such as moving window size and copula selection. With the careful adjustment, a model selection process based on backtesting results should be conducted. For future studies on the high-dimension dependency structures, a careful evaluation is highly recommended for the difference between the simulated error terms and the realized error terms. As the dimensionality increases, the difference for fitted multivariate dependency structures can be very vulnerable in terms of precision and accuracy, compared to the realized dependency structure.