Portfolio Tail Risk: A Multivariate Extreme Value Theory Approach

This paper develops a method for assessing portfolio tail risk based on extreme value theory. The technique applies separate estimations of univariate series and allows for closed-form expressions for Value at Risk and Expected Shortfall. Its forecasting ability is tested on a portfolio of U.S. stocks. The in-sample goodness-of-fit tests indicate that the proposed approach is better suited for portfolio risk modeling under extreme market movements than comparable multivariate parametric methods. Backtesting across multiple quantiles demonstrates that the model cannot be rejected at any reasonable level of significance, even when periods of stress are included. Numerical simulations corroborate the empirical results.


Introduction
Events of the Global Financial Crisis, the Eurozone sovereign debt crisis, and the ongoing COVID-19 pandemics increase the investors' awareness of the need to assess the risk of extreme losses and the necessity of protection against future market meltdowns. To quantify financial risks, investors and regulators commonly resort to risk measures based on quantiles of distributions of losses or returns, such as Value at Risk (VaR) or Expected Shortfall (ES). However, the most significant problem with conventional approaches to assessing the quantile-based risk measures is their poor performance in characterizing the "tail" behavior of the distribution of returns. Any model aiming to mitigate severe losses should specifically focus on the tail risk.
The EVT approach seems like a natural choice for risk modeling: its implementation is relatively easy, and the method is based on a few assumptions required for the asymptotics to work. Although straightforward in the univariate case, EVT becomes of limited practical use when applied to a portfolio of financial instruments or several risk factors at once. Investors usually deal with many risk factors simultaneously, and hence an appropriate method requires a multivariate approach rather than a univariate one based on portfolio-level returns. As shown by Longin and Solnik [22], Campbell et al. [23], and many authors since, multivariate modeling is particularly vital during bear markets, when risk factors exhibit extreme co-movements. It also allows for risk attribution and additive decomposition of VaR, which is essential for risk management and economic capital allocation in financial institutions.
Defining portfolio-level dynamics of risk factors under extreme market conditions has so far been a daunting task. A seemingly obvious technique involves a multivariate version of the EVT based on the multidimensional limiting relations (see, for example, [14,[24][25][26][27]). However, the model complexity increases faster than linearly with the number of risk factors in that case. Alternatively, the joint distribution of returns can be seen as a product of marginal distributions and a copula. McNeil and Frey [28] and Nyström and Skoglund [29], for example, describe the copula approach to the assessment of the extreme co-dependence structure of risk factors. The copula-based methods introduce an additional model risk, inherent in the assumption of a specific analytical form of the co-dependence function. In addition, they become quite intractable with an increase in dimensionality. Moreover, a typical copula method for multivariate EVT, such as the one described in Nyström and Skoglund [29], requires an additional simulation step to retrieve the residuals from the joint distribution, given the fitted marginals and parameters of the copula. The dimensionality problem in the multivariate setting can be tackled by vine pair copulas, as in Aas et al. [30], Chollete et al. [31], Yu et al. [32], or Trucíos et al. [33].
This paper introduces a multivariate EVT method for risk assessment based on separate estimations of the univariate model. Instead of estimating joint n-dimensional distributions (using multidimensional limiting relations, copulas, or otherwise), the method proposed here works with n orthogonal series that are approximately independent and identically distributed. These series are obtained from the eigenvectors of the joint return series, filtered to remove autocorrelation and heteroskedasticity. The filtering can be achieved, for instance, by assuming that the conditional covariance matrix follows a stationary n-dimensional model from the orthogonal GARCH family, as in Alexander [34] or van der Weide [35]. The tails of n independent univariate series of residuals are then fitted to a generalized Pareto distribution, allowing for closed-form expressions for VaR and ES. As an illustration, the technique is applied to a sequence of daily returns on stocks comprising the Dow Jones Industrial Average. For an out-of-the-sample analysis, daily VaR and ES estimates are compared to the actual portfolio losses over a more extended period, covering the extreme market co-movements during the last quarter of 2008. The results indicate that the method performs well in capturing extreme events jointly across the entire portfolio. The in-sample tests show that the proposed approach is better suited for portfolio risk modeling under extreme market movements than comparable multivariate parametric methods, such as orthogonal GARCH with normal or t-distribution. In addition, backtesting across multiple quantiles demonstrates that the method cannot be rejected at any reasonable level of significance, even for out-of-sample windows that contain stress periods.
This paper contributes to the literature in several important ways. To the best of our knowledge, it is the first method that extends the EVT to the multivariate case by combining the relative simplicity of univariate EVT and orthogonalization of return residuals. By doing so, we can capture tail correlations and extreme co-movements. The proposed method is universal, as it can be applied to any asset class, and it can incorporate any volatility model from the GARCH family.
The remainder of the paper is organized as follows: Section 2 presents the theoretical background behind the proposed multivariate EVT approach and the estimation methodology used in this paper. Section 3 describes the data, provides an example of VaR and ES estimation, and compares the method's in-sample forecasting ability with other similar approaches. Section 4 illustrates the out-of-sample performance of the method through backtesting. Section 5 demonstrates the model performance on a set of simulated data. Concluding remarks are given in Section 6.

Theoretical Framework
Two important theorems are related to the univariate EVT. The first one shows that order statistics of independent and identically distributed (i.i.d.) random variables converge in distribution to the generalized extreme value (GEV) distribution [36,37]. The usual approach then is to fit the GEV distribution to the sequence of maxima of portfolio returns or their residuals, as described for example in Smith [14]. An alternative method is based on exceedances over threshold. The following theorem, first stated by Pickands [38], gives the asymptotic form of conditional distribution beyond a threshold: to be the distribution of exceedances of Z t over the threshold u. Let z F be the right endpoint of the distribution F, possibly a positive infinity. Then, if F is such that GEV distribution is well-defined, there are constants ξ ∈ R and β := β(u) ∈ R + such that is the generalized Pareto (GP) distribution, while the subscript + denotes the positive part of a function, i.e., f + (x) ≡ max ( f (x), 0) for any function f .
The challenges associated with implementation of EVT-based models are well known in the literature. Since the Theorem 1 defines a limiting distribution for the extreme returns, we need a sufficient number of observations and a large enough threshold u. There are different methods of making this choice, and some of them are examined in Bensalah [13] and Gabaix et al. [39]. In addition, the limit theorems hold only if the extreme observations Z t are i.i.d. Therefore, we cannot apply the results of univariate EVT to returns on financial assets directly without taking into account their autocorrelation and heteroskedasticity.

Estimation Methodology
The method works in several steps. We start from a multivariate series of returns in a portfolio consisting of n assets. Next, we apply orthogonal GARCH to obtain a set of random variables that are approximately i.i.d. Then, we fit the tails of n independent univariate residual series to a generalized Pareto distribution to obtain the quantiles for each component. Finally, we substitute these quantiles into a formula for VaR and ES. Sections 2.2.1-2.2.6 describe the methodology in detail, while Section 2.2.7 provides its summary.

Orthogonalization
The ultimate goal is to adapt the EVT approach to a portfolio consisting of n assets in order to capture their tail risk and their extreme co-movements. Before applying any of the results of EVT outlined in Section 2.1, we have to construct a set of cross-sectionally uncorrelated random variables. A natural choice is to work with the eigenvectors of the unconditional covariance matrix of returns. Definition 1. Define ε t to be an n-dimensional random vector whose components ε t,i have zero mean for each i = 1, 2, . . . , n. Let V ∞ = E(ε t ε t ) be the n-by-n unconditional covariance matrix of ε t and P be the corresponding orthogonal matrix of normalized eigenvectors, where the prime symbol ( ) denotes the transpose of a matrix. The eigenvalue decomposition of V ∞ is given by where Λ is a diagonal matrix of the eigenvalues of V ∞ , ordered by descending values, Then, is called the vector of principal components of ε t , for any t. The i-th element of the vector z t is called the i-th principal component of ε t .
Note that E (z t ) = 0 and var (z t ) = 1 n , which follows from E (ε t ) = 0 and V ∞ = LL , respectively. Hence, z t are cross-sectionally uncorrelated and each component has a unit variance. Since ε t = Lz t , each coordinate of ε t can be written as a linear combination of the principal components, where L ij are the elements of L. The method can be easily adapted to fewer than n risk factors (see Appendix A).
An alternative linear transformation that can be used to separate multivariate returns into additive subcomponents is the independent component analysis. This technique has found its successful applications in finance: Back and Weigend [40] use it to extract the structure from stocks returns, while Moody and Wu [41] and Moody and Yang [42] apply it to foreign exchange rates. The main goal of independent component analysis is to minimize statistical dependence between the transformed vectors. Unlike principal component analysis, these vectors are neither orthogonal nor ranked by the variance they explain, and the linear transformation matrix analogous to L cannot be obtained in a closed form. Therefore, independent component analysis is not entirely suitable for the problem treated in this paper. A good comparison between the two methods can be found in Hyvärinen and Oja [43].

Filtering
Orthogonalization transforms a cross-sectionally correlated series into a set of uncorrelated ones. We also have to filter out any autocorrelation and heteroskedasticity from the series. As the end-result, we will obtain sequences of conditional residuals that are orthogonal, and (approximately) serially uncorrelated and identically distributed.
Specifically, we can assume that, for each asset i = 1, 2, . . . , n, the series of log returns y t,i := ln(S t,i /S t−1,i ), where S t,i denotes the price at time t, follows a process: where µ t,i is the conditional mean of y t,i . Conditionally on the information available at t − 1, the vector of residuals ε t := [ε t,1 ε t,2 . . . ε t,n ] has a zero mean and a conditional covariance matrix V t : where, for any t, the matrix V t is positive definite and measurable with respect to the information set F t−1 , a σ-algebra generated by the past residuals {ε 1 , ε 2 , . . . , ε t−1 }. To capture the volatility clustering effect, we can assume that the conditional covariance matrix follows a model from the GARCH family. The standard GARCH(p,q) is sufficient to capture most of the clustering, and-to some extent-excess kurtosis. The leverage effect can be taken into account, for example, by assuming that the conditional residuals follow an asymmetric distribution, such as skewed Student's t. Alternatively, we can model the asymmetry explicitly in the equation followed by the conditional covariance matrix. Without losing generality of the proposed method, we can assume that the conditional covariance V t follows a multivariate asymmetric GARCH(p,q) of Glosten et al. [44], also known as the (multivariate) GJR-GARCH(p,q), or simply GJR(p,q): where Ω, A 1 , . . . , A p , Θ 1 , . . . , Θ p , B 1 , . . . , B q are constant, positive semidefinite n-by-n matrices of parameters, and for any t. As usual, the matrices A s in (6) measure the extent to which volatility shocks in previous periods affect the current volatility, while A s + B s measure the rate at which this effect fades away. The terms proportional to matrices Θ s capture the impact of asymmetric return shocks to volatility. For any t, the unconditional covariance matrix of ε t is given by where 1 n denotes an n-by-n identity matrix. Covariance stationarity of the GJR(p,q) process (6) is assured by setting the matrix to be positive definite. Note that the multivariate GJR process (6) contains the ordinary multivariate GARCH as its special case. There are also many plausible alternatives. For example, the APARCH model of Ding et al. [45] includes ARCH, GARCH, and GJR-GARCH models as special cases, as well as four other ARCH extensions. These four extensions are the TS-GARCH model of Taylor [46] and Schwert [47], Log-ARCH model of Geweke [48] and Pentula [49], N-ARCH model of Higgins and Bera [50], and T-ARCH model of Zakoïan [51]. Similarly, the EGARCH model of Nelson [52] can be used to capture the leverage effect in a slightly different manner. The method developed in this paper can be applied straightforwardly with any of the multivariate models from the GARCH family. Cross-sectional correlations are reflected in the off-diagonal terms of matrices V t and E t . This, in turn, makes the parameter matrices non-diagonal. In total, one would have to estimate (1 + 2p + q)(n + 1)n/2 different parameters. Clearly, this number explodes as we increase the number of assets in the portfolio. However, this is only one aspect of the problem. The other is that we cannot apply the results of univariate EVT to conditional residuals ε t,i directly. However, we can switch to the orthonormal basis of principal components by applying the linear transformation (2) to the conditional residuals ε t . In this basis, Equation (6) reads: where M := L −1 M L −1 for any square matrix M. In particular, be the conditional covariance matrix of principal components. As in Alexander [34], it is reasonable to assume that the matrix V t is diagonal, since the eigenvectors z t are orthogonal. Then, the process given by Equation (7) can be estimated separately for each principal component. This gives a set of n independent scalar equations of the form where, in general, M i := M ii is the i-th diagonal element of the matrix M, i being 1, 2, . . . , n for the first, second, . . ., n-th principal component, respectively.
Once we estimate the parameters Ω, for t ≥ max{p, q}, to retrieve the series of conditional covariance matrices in the original basis of log returns. This allows us to forecast VaR and ES in a multivariate framework, and for an arbitrary portfolio.

Estimating Independent Univariate Excess Distributions
Theorem 1 states that, for a large class of underlying excess distributions, the exceedances over threshold converge in distribution to GP as the sample size increases and the threshold is raised. Thus, the GP is the natural model for the unknown excess distribution above the threshold u, and we may conjecture that for any x satisfying 0 ≤ x < z F − u. Assuming that we have a set of realizations for the variable Z t , we can choose a sensible threshold u and estimate parameters ξ and β. If there are T u out of a total of T observations that exceed the threshold, the GP distribution will be fitted to the T u exceedances.
In the literature, several estimators have been used to fit the GP parameters. The most popular ones are the maximum likelihood (ML) and the Hill estimator. The ML estimator assumes that, if the tail under consideration strictly follows a GP distribution, the likelihood function can be written in a closed-form. The estimators of the parameters ξ and β are then obtained using the standard ML approach. Provided that ξ > −1/2, the ML estimator of the parameters is consistent and asymptotically normal as the length of the series tends to infinity. On the other hand, the Hill estimator is based on a combination of the ML method and a semi-parametric result, which describes the scaling of F when ξ > 0 (see, for instance, [53]). Here, we will use the ML estimator due to its performance and universality, as discussed in Nyström and Skoglund [15].

Tails of Univariate Distributions
By combining Theorem 1 and Equation (10), we can write for z > u. The only additional element we require to construct a tail estimator is F(u).
The simplest choice would be to use the non-parametric method and take the obvious empirical estimator, F(u) = 1 − T u /T. Alternatively, we can find T u that is closest to a predetermined F(u). Thus for example, in a sample of T = 1000 observations, F(u) = 0.90 will correspond to T u = 100. The threshold is then set to u = Z 900 , if {Z t } T t=1 are ordered in an increasing way. Combining the empirical estimate F(u) with the ML estimates of the GP parameters, we obtain the tail estimator: Note that, when the scale parameter β tends to infinity, G ξ,β (·) vanishes and the tail estimator converges to the empirical one for any z. Thus, the tail estimator in (11) can be viewed as the non-parametric estimator augmented by the tail behavior captured by the GP distribution.

Estimating Univariate Var and Es
We formally define VaR in the following way.
Value at Risk is the α quantile of the distribution F: where α ∈ (0, 1) and F −1 is the inverse of F.
The usual critique of VaR as a risk measure is that it is not coherent and ignores the structure of losses beyond a specific quantile. (See, for example, Artzner et al. [54,55], Acerbi and Tasche [56,57] and Szego [58].) In order to overcome these drawbacks, Expected Shortfall (ES) is often used as an obvious alternative. It is defined as the conditional expectation of loss that surpasses a fixed level of VaR. Formally, under the same assumptions as in Definition 2, we can define ES as: As such, ES takes into account tail risk and satisfies the sub-additivity property, which assures its coherence as a risk measure.
Let u + ≡ u be the upper-tail threshold, and let the lower-tail threshold u − be defined symmetrically, that is, by F(u − ) = 1 − F(u + ). Then, for a given upper-tail probability α + > F(u + ) or a given lower-tail probability α − < F(u − ), the general form of the VaR estimate is where the subscript + (−) refers to parameters in the upper (lower) tail. Similarly, the general form of the ES estimate is (see, for example, [12]). As the tail becomes heavier, the tail index increases (equivalently, ξ ± → 0), and ES becomes progressively greater than VaR.

Portfolio-Level VaR and ES
Given the upper-and lower-tail quantiles α ± , the quantile of i-th principal component h time steps ahead is given by where τ = T + h and F −1 i (·) is the inverse of the univariate probability function for the set of principal components {z t,i } T t=1 , given by Equations (12) and (13). Our final goal is to estimate VaR and ES for a portfolio of n assets. Denote by a the vector of portfolio positions, in monetary units. This, among other things, facilitates the treatment of portfolios with short positions. Then, h-steps-ahead portfolio VaR is given by where Q ± τ is a real symmetric matrix that can be decomposed as The intuition behind Formula (16) is fairly simple. The first term, a µ τ , represents the forecast of the expected portfolio return. The second term is the upper-or lower-tail bound of the confidence interval for total portfolio return. It is obtained in three steps: first, we use (15) to calculate the quantiles z ± τ,i , i.e., individual contribution to VaR for each principal component, and stack them into a diagonal matrix D ± τ ; second, we use L to transform this matrix back to the original basis of log returns; finally, we take a square root of the matrix "sandwich" consisting of the vector of portfolio holdings a and transformed quantiles Q ± τ to obtain the bounds of confidence intervals around the expected return. If all principal components happen to be identically distributed (which is unlikely, even if all come from the same family of distribution functions), then the square-root term is simply the α-quantile of this distribution times the h-steps-ahead forecast of the portfolio variance.
A similar rationale can be applied to derive the forecast of portfolio ES: where, by analogy with Equations (14)- (16),

Summary of the Methodology
We summarize the methodology through the following procedure: Step 1. Start by performing the orthogonalization on portfolio returns. In other words, use the eigenvalue decomposition of the unconditional covariance matrix V ∞ of return residuals ε t to obtain the matrix L following Definition 1. Then, use Equation (2) to obtain the set of orthogonal vectors z t . This step assures cross-sectional independence and allows working with univariate series.
Step 2. Filter out the heteroskedasticity on each component from step 1 to obtain serially i.i.d. random variables. In other words, use the set of orthogonal vectors z t from step 1, calculate E t = z t z t and run a GARCH-like process such as one given by Equation (8) on E t . This results in a set of univariate estimates of the conditional variance V t,i . This step is a prerequisite for using the univariate EVT.
Step 3. Use the estimator from Equation (11) to fit the tails of each of the independent series from step 2. This gives the parameters of GP distribution.
Step 4. Substitute the parameters from step 3 into Equation (12) or (13). This step allows us to calculate the quantile for each component, Equation (15).
Substitute the quantiles from step 4 into the closed-form expressions from Equations (16) and (17) to obtain portfolio VaR and ES.
Steps 1 and 2 are analogous to the orthogonal GARCH of Alexander [34]. Steps 3 and 4 are typical for a univariate EVT approach. A combination of steps 1-4 with step 5 is essential for the multivariate EVT approach.
As discussed, the method is based only on a few assumptions that allow the exploitation of asymptotical properties. From Theorem 1, we can see that it is sufficient to have a set of independent univariate random variables {Z t } T t=1 drawn from the same distribution function F. We can obtain the set of univariate random variables by merely applying a simple eigenvalue decomposition of the unconditional covariance matrix V ∞ in step 1. This step is well defined as long as the original multivariate series has a finite unconditional covariance matrix, a reasonable assumption for the financial time series. The set of i.i.d. random variables having the same distribution function F can be obtained through a GARCH-like filter in step 2. The filter can be applied provided that the covariance stationarity condition has been met. This condition is merely a restriction on the model parameters. Other assumptions are not required for the method to work since it relies on the transformed random variables' asymptotic property.

Illustration of the Method
The empirical results that follow are based on daily closing prices of the thirty stocks that constituted the Dow Jones Industrial Average at the end of 2011, gathered from Thomson Reuters Datastream. The prices are adjusted for stock splits and dividend payments. The observation period covers all trading days between and, including 2 January 2001, and 30 December 2011, there is a total of 2766 returns per stock. Figure 1 displays the returns of an equally-weighted portfolio, which will be used in the later examples.
Following the procedure described in Section 2.2, we first estimate the upper-and the lower-tail parameters of the univariate GP distribution, ξ ± and β ± , for each of the principal components. The upper and lower thresholds, u + and u − , are determined by F(u + ) = 0.90 and F(u − ) = 0.10, respectively. (g choice of threshold u ± is of minor importance for this illustration. For instance, using F(u − ) = 0.05 instead of 0.10 at α = 0.001 changes the level of VaR forecast by less than a percent in relative terms.) The lower-tail VaR and ES forecasts can be found easily from Equations (16) and (17). The conditional covariance matrix V τ,i in (15) is obtained using the steps described in Equations (6)- (8).
Here, we illustrate this procedure, assuming that the conditional covariance matrix in Equation (6) follows a multivariate GJR(1,1) process, i.e., p = q = 1, or its multivariate GARCH(1,1) restriction, in which Θ s ≡ 0. Table 1 shows ten-day VaR and ES forecasts for several confidence levels reported as a percentage of the total portfolio value. (We choose the convention that the negative values obtained in Equation (15)   We can compare different multivariate volatility models in terms of their ability to fit the conditional distributions of principal components. Here, we use three different distributional assumptions: normal, Student's t, and GP. The results for normal and Student's t-distribution are calculated from the usual procedure of Alexander [34]. The GP results are based on our multivariate EVT method, as described in Section 2. The results are summarized in Table 2. The second column of Table 2 shows the cumulative R 2 , i.e., the fraction of variance explained by each of the thirty eigenvectors. The remaining columns contain Kolmogorov-Smirnov test statistics for the univariate probability distribution of standardized eigenvectors z t,i . The methods based on multivariate EVT appear to be the only ones that cannot be rejected at any reasonable significance level for all thirty eigenvectors-everything else being equal. Due to the fat tails of returns, conditional normality can be rejected for almost all eigenvectors. The assumption of Student's t conditional distribution can be safely rejected for the first principal component, which is the one that explains almost half of the variations of returns in this sample. It is also rejected at more restrictive levels of significance for the remaining eigenvectors.
To visualize the model's ability to fit the lower tail, we can plot the distribution of exceedances F u (x). Figure 2 shows log 10 F u (x) for the first four principal components, which jointly contribute to more than two-thirds of variations in portfolio returns. Each graph compares the empirical distribution function tails and three different distribution functions obtained from the orthogonal-GJR(1,1) residuals. The fitted GP distribution G ξ,β (·) is obtained from the model proposed in this paper. Consistently with Table 2, it provides the best fit to the actual data.

Out-of-Sample Performance
To analyze the out-of-sample performance of the proposed method, we use a 1000-day rolling-window sequence of one-day VaR forecasts and compare them to the corresponding actual returns of the equally-weighted portfolio for each day. In other words, we use the sequence of returns between 3 January 2001, and 14 January 2008, to obtain the VaR forecast for 15 January 2008. Then, we shift the window by one day and use the sequence of returns between 4 January 2001, and 15 January 2008, to obtain the VaR forecast for 16 January 2008, and so on, until the last sequence becomes the one corresponding to the period between 28 December 2004, and 29 December 2011, which is used to obtain the VaR forecast for 30 December 2011. We use five different VaR quantiles (α = 0.001, 0.005, 0.010, 0.050, 0.100) and perform a backtesting procedure on the same models as in Section 3: orthogonal multivariate GARCH(1,1) and GJR(1,1) processes, under three distributional assumptions (normal, Student's t, and GP). Table 3 compares the actual and expected number of VaR violations for each quantile. We may notice that the number of violations is above the expected value, indicating that the models, on average, underestimate the tail risk. The assumption of conditional normality displays far too many exceptions in the extremes. The models based on the multivariate EVT method with GP distribution have violations whose number is closest to the expected. Table 3. Actual and expected number of lower-tail VaR violations for an equally-weighted portfolio of 30 stocks.

mv GARCH
where N obs k and N exp k are, respectively, the observed and the expected number of violations in the k-th bin. The Q statistic converges in distribution to a χ 2 with K − 1 degrees of freedom. The results of the test are summarized in Table 4 Table 4 indicates that the proposed multivariate EVT method cannot be rejected whenever the comparable methods can. For example, if the orthogonal multivariate GJR process is used, the model based on conditionally normal residuals can be rejected at 5% significance, but the models based on conditionally t-distributed residuals and the multivariate EVT method cannot. However, if the orthogonal multivariate GARCH process is used, then, at 5% significance, the models based on conditionally normal and t-distributed residuals can both be rejected, while the multivariate EVT method cannot. Similar results are obtained with Christoffersen [60] conditional coverage test across separate quantiles (available upon request). In Appendix B, we illustrate the method using the exchange-rate data. We obtain results similar to those shown in Table 4.

Numerical Simulations
To demonstrate the proposed method's ability to incorporate various data patterns, we will apply it to simulated data in this section. To this end, we will simulate co-dependent multivariate series using different copulas: a Gaussian copula with the ''ordinary" (i.e., Pearson's) linear correlation, a Gaussian copula with Kendall's τ correlation, a t copula with Kendall's τ correlation, and three Archimedean copulas-Clayton, Frank, and Gumbel. The copula parameters will be calibrated using the stock returns from Section 3. We will also use the stock data to obtain the marginal distributions, applying kernel smoothing on the empirical distribution functions. For the Gaussian copulas and the t copula, we apply the full 30-dimensional multivariate structure. Since the applied Archimedean copulas are bivariate, we use two equally-weighted sub-portfolios containing 15 return series each. For every copula, we simulate 10,000 data points in each series.
We again follow the procedure outlined in Section 2.2, with the same choice of thresholds. We then apply Equations (16) and (17) to estimate lower-tail VaR and ES. For the sake of tractability, we focus on the multivariate GJR(1,1) process for the conditional covariance matrix. Table 5 summarizes the ten-day VaR and ES forecasts for the same confidence levels as in Table 1, again reported as a percentage of the total portfolio value. We can notice that the numbers do not change drastically with respect to those reported in Table 1. They are also similar across different joint distributional assumptions.
The method's flexibility in fitting the lower tail across various joint distributions of the data can be visualized in Figure 3. We show the log 10 F u (x) of the first principal component for each of the six copula types obtained from the orthogonal-GJR(1,1) residuals. The GP distribution G ξ,β (·) represents an excellent fit to the simulated data.  The out-of-sample performance of the proposed method on the simulated data can be assessed following the same idea as in Section 4. Again, we apply a 1000-day rolling-window sequence of one-day VaR forecasts. We use the same five VaR quantiles as in Section 4 and perform a backtesting procedure using the orthogonal multivariate GJR(1,1) processes, under six joint distributional assumptions (i.e., simulated copulas). Table 6 shows the comparison between the actual and the expected number of VaR violations across quantiles. In all cases, the number of violations is very close to the expected value. We can see this more formally through Pearson's Q test results, shown in Table 7. The proposed multivariate EVT method cannot be rejected for any of the six joint distributional assumptions used to obtain the simulated datasets.

Conclusions
This paper proposes a procedure for assessing portfolio tail risk based on Extreme Value Theory, without the necessity to use multivariate limiting relations. The method applies separate estimations of the univariate EVT model. It works with a set of orthogonal conditional residuals obtained from the principal components of the covariance matrix of returns. Autocorrelation, heteroskedasticity, and asymmetry inherent in the original return series can be removed by assuming a GARCH-type process for conditional variances of principal components. In this way, we can obtain a set of independent and identically distributed random variables, which is a prerequisite for any univariate EVT approach. The tails of the univariate distributions have been modeled by a generalized Pareto distribution of peeks over the threshold, while the interiors are fitted with an empirical distribution function. The estimation of parameters can be obtained using maximum likelihood, keeping the method computationally as fast as its orthogonal-GARCH counterparts. The method also provides closed-form expressions for Value at Risk and Expected Shortfall.
As an illustration, the method has been applied to a sequence of daily stock returns of Dow 30 companies. In-the-sample tests show that the proposed approach is more reliable for modeling extreme market movements than comparable multivariate parametric methods, such as orthogonal GARCH-type models with normal or t-distribution. Backtesting on 1000 out-of-the-sample observations indicates that VaR forecasts obtained by the proposed method perform well across multiple quantiles: the values of test statistics suggest that the null hypothesis of an equal number of observed and predicted VaR violations could not be rejected for all relevant VaR confidence levels.
The importance of these results is further highlighted because the backtesting period covers extreme market movements that had followed the events of the Global Financial Crisis.
Numerical simulations have corroborated the empirical results. Using six different joint distributional assumptions, where various copula families drive the co-dependence in the data, we show that the proposed multivariate EVT method represents an excellent fit for the simulated data tails. The out-of-the-sample performance has also been confirmed via simulations: the method cannot be rejected for any of the six joint distributional assumptions.
Even though EVT is a natural choice for tail risk modeling, its major shortcoming is the complexity of extension to the multivariate case. This paper proposes a way to overcome this problem by combining the simplicity of univariate EVT and orthogonal GARCH while capturing tail correlations and extreme co-movements. The proposed method can be applied to any asset class, and it is straightforward to modify its filtering process to include any volatility model from the GARCH family. Hence, it has many potential applications in volatility forecasting, asset allocation under VaR or ES constraints, risk management of complex portfolios, or assessing individual assets' contributions to systemic risk.

Conflicts of Interest:
The author declares no conflict of interest.

Appendix A. Working with Fewer Than n Principal Components
Definition 1 implies that the fraction of total variation in ε t explained by the j-th principal component is given by This property leads to another convenient feature of the principal component approach. Namely, if low-ranked components do not add much to the overall explained variance, which is often the case when the number of risk factors is relatively low, we can work with a reduced number of m principal components, where m < n. The first m components will then explain ∑ m j=1 λ j ∑ n k=1 λ l 1 of the variation in ε t . In that case, L is replaced by an n-by-m matrix L m , where P m is an n-by-m matrix of the first m normalized eigenvectors, and is a diagonal matrix of the first m eigenvalues. The m-dimensional vector of the first m principal components of ε t is then given by for any t, where L −1 m is the pseudo-inverse of L m . It is straightforward to generalize the multivariate EVT approach presented in the paper to the case of m < n principal components. Using definition (A1),  (7), (8), (15), and (16) maintain the same form.

Appendix B. Illustration Using Exchange-Rate Data
As an additional robustness check of our method, we present here the results for the exchange-rate data. On average, currency portfolios have lower rates of return than stocks but are significantly more volatile and notoriously fat-tailed. Here, we use daily data on 20 major currencies against the US Dollar (USD), retrieved from Refinitiv-Eikon. We used the following currencies: Australian Dollar, Brazilian Real, Canadian Dollar, Swiss Franc, Chinese Yuan, Euro, Pound Sterling, Hong Kong Dollar, Indian Rupee, Japanese Yen, South Korean Won, Mexican Peso, Norwegian Krone, New Zealand Dollar, Russian Ruble, Swedish Krona, Singapore Dollar, Turkish Lira, New Taiwan Dollar, and South African Rand. The data consists of middle exchange rates between 1 December 2000, and 1 December 2020, a total of 5217 returns per currency. Figure A1 shows the returns of an ''equally-weighted" portfolio, which we will apply in our analysis, analogous to Section 3. We use a 2000-day rolling-window sequence of one-day VaR forecasts to analyze the out-of-sample performance on the exchange-rate data. As in Section 4, we compare these results to the corresponding actual returns of the equally-weighted portfolio for each day using the same five quantiles (α = 0.001, 0.005, 0.010, 0.050, 0.100) and the same multivariate methods. The results are summarized in Tables A1 and A2. Table A1 compares the actual and expected number of VaR violations for each quantile. The multivariate EVT method with GP distribution has the number of violations that is again closest to the expected. We can see this more formally from the backtesting results in Table A2. Again, the proposed multivariate EVT method cannot be rejected whenever the comparable methods can. This result indicates that the method is robust across different asset classes, observation periods, and backtesting horizons.