When the well known univariate risk measure analysis is generalized into the multivariate setting, a number of complex theoretical and applied problems appear. This emerging theory just extends the univariate case for VaR estimate by using copulas, theory of extreme value, Monte Carlo method, historical simulation, variance and covariance analysis, historical simulation, among many others. However, some of the univariate translations became unrealistic and are based on inappropriate assumptions. For example, under the restriction of a perfect dependence, the simple summation method computes the total risk by summation of the stand alone risks; a preservative method with statical benefit (see for example, [Embrechts et al. (2013) ] and [Li et al. (2015) ]). In a similar way, for a large number of assets, the variance covariance method fails, because the estimation of the corresponding matrix is extremely cumbersome due to the high amount of correlations, see [McNeil and Embrechts (2015) ]. Now, in the context of risk management, the multivariate theory of the extreme value (EVT) (see [McNeil (1999) ]) and the popular technique of multivariate copulas (see [Embrechts et al. (2002) ]) are useful in some scenarios for VaR estimation in portfolios. In particular, copulas method attains a robust structure for dependence in financial time series by producing joint distributions with known non gaussian marginal distributions. Modelling the marginal distributions via copulas allows VaR computations with a better performance than the classical methods; but it involves some untractable assumptions in the context of risk measures which are difficult to elucidate; a similar quotation for the multivariate extreme value theory are also addressed by [Jin and Lehnert (2018) ] and [Barone et al.(2015) ].
Inspired by the above discussions and some interesting insights collected by [Jianping et al. (2012) ] for risk models in banking industry, this work proposes risk measures based on the Wasserstein barycenter. The new method considers a reliable risk measure based on distances among probabilistic models. The underlying suitable probability laws obey for example, opinions, beliefs and estimates of data sources, in the context of the financial risk. Explicitly, a far away concept in probability theory is brought into the financial models by proposing the named Fréchet risk measures; which are calibrated by certain metrization of the probability measure space. In this case, the well studied metric of Wasserstein supports the method and provides fundamental connections for the rising concept of barycenter in the sense of Agueh and Carlier in [1] . A seminal work for a number of generalizations and applications, see for example, [Bigot et al. (2018) ], [Alvarez et al. (2018) ], [Le Gouic and Loubes (2017) ] and the references therein. A crucial aspect underlies here for the new method: the proposed barycenter remains invariant under a class of location and scatter set of (finite or infinite) set of probabilities.
The Wasserstein metric has enriched notably the risk management literature, see for example [Kiesel et al. (2016) ] and [Feng Yu and Erik Schlogl (2018) ]. In particular, as a canonical metric under well defined assumptions, the named robust risk management has been studied under the Wasserstein metric. As a consequence, this work represents risk measures via statistical functionals by hybridizing robustness and continuity under the Wasserstein metric. Then, several financial applications of the Wasserstein metric can be obtained in real time series, where the classical methods provide an excellent scenario for the correctness of the predictions.
The above discussion is organized in this paper as follows: preliminaries about Wasserstein barycenter are given in Section 2, in order to establish the existence and uniqueness theorem for barycenter of distributions under a class of location and scatter distributions. Then Section 3 defines the Wasserstein Barycenter in risk measures and results for VaR and CVaR estimation are given under the addressed family of scale and location distributions. Finally, Section 4 applies the the results in a portfolio consisting of two assets (the Nasdaq and the S&P500 stock indices).
This apart provides the necessary background about the barycenter in a Wasserstein space. First the quadratic transportation cost is presented, in order to introduce the 2-Wasserstein distance. Then, the basic multivariate Value At Risk models are revisited in terms of the Wasserstein barycenter. Finally, following [Alvarez et al. (2018) ] and [Cédric Villani. (2008) ], the fundamental results on Wasserstein barycenters for measures under a class of scatter and location distributions are provided. Start with P 2 (R d ) as the set of all probability measures defined on R d with a finite second order moment. Denote P 2,ac (R d ) as the subset of absolutely continuous measures and consider (Ω, σ, P ) as a generic probability space. If µ, ν in P (R d ), are two measures, then P(µ, ν) will denote the set of all probability measures π in the product set R d × R d . Here, µ and ν are the corresponding first and second marginals. Now, for two measures µ, ν in P(R d ), the quadratic transportation cost between µ and ν (also referred as the transportation cost with a quadratics cost function) is defined as follows
The transportation cost with quadratics cost function endows the set P 2 (R d ) with the metric called 2-Wasserstein distance or Monge-Kantorovich distance metric, which is given by W 2 (µ, υ) = T 2 (µ, ν) 1 2 . When d = 1, the 2-Wasserstein distance in the real line is just given by the quantile-like expression:
here F −1 ν and F −1 µ ) are the quantile function of ν and µ, respectively. Now, in the Euclideam space, the barycenter of points
In fact, is the unique minimizer
Motivated by the Euclidean version, the Wasserstein barycenter can be defined as follows.
Definition 2.1. Let µ 1 , ..., µ N be random probability measures over R d , endowed with positive weights λ 1 , ...λ N , with N j=1 λ j = 1. The measure µ ∈ P 2 (R d ) is a Wasserstein barycenter, if µ is a minimizer of the functional
This fact will be denoted by µ B (λ) ∈ Bar((µ j , λ j ) 1≤j≤N )
Definition 2.2. We say that the measure µ ∈ P 2 (R d ) is a Wasserstein barycenter for the random probability measures µ 1 , ..., µ N over R d , endowed with positive weights λ 1 , .
We will write µ B (λ) ∈ Bar((µ j , λ j ) 1≤j≤N )
Empirical consistency of the Wasserstein barycenter has been studied in [Agueh and Carlier (2011) ], [Boissard et al. (2015) ] and [Le Gouic and Loubes (2017) ].
For introducing a fundamental result, we recall the definition of a location-scatter class.
Definition 2.3. If M + d×d denotes the set of d × d positive definite matrices and X 0 is a random vector with measure µ 0 ∈ P 2,ac (R d ), then the set F(µ 0 ) of probability laws defined by F(µ 0 ) := {l(AX 0 + m) : A ∈ M + d×d , m ∈ R d } is a location-scatter family induced by positive definite affine transformations from µ 0 .
The Wasserstein barycenters of measures on a location-scatter family satisfies the following remarkable property, see [Alvarez et al. (2018) ].
Lemma 2.1 (Theorem 3.10., [Alvarez et al. (2018)] ). Let µ 0 ∈ P 2,ac (R d ), and µ ∈ W 2 (P 2 (R d )), assume that for every ω ∈ Ω, the measure µ ω ∈ F(µ 0 ). Then the unique barycenter, µ of µ also belongs to F(µ 0 ). The mean ofμ is m := m ω P (dω), and the covariance matrix, Σ, is the only positive definite root of the equation
This result means that the Wasserstein barycenters are closed respect the locationscatter class.
An interesting case follows for N Gaussian measures on R d :
Lemma 2.2 (Theorem 2.5., [Alvarez et al. (2018)] ). Consider N Gaussian measures µ 1 , ..., µ N on R d with corresponding means m 1 , ..., m N and positive definite covariances Σ 1 , ..., Σ N , and let λ 1 , ..., λ N be positive weights with N j=1 λ j = 1. Then the unique barycenter of the normal measures µ 1 , ..., µ N is the Gaussian distribution with mean m λ = N j=1 λ j m j and covariance matrix Σ, which is the only positive definite root of the equation
According to [Alvarez et al. (2018) ], Wasserstein barycenters inherits the strong computational problems of the classical optimal transportation. However, in the real line some explicit distributions can be obtained.
.., F −1 N be the quantile functions corresponding to µ 1 , ..., µ N in the real line. Thus the barycenter of µ 1 , ..., µ N is the probability with quantile function
This notable aspect will be used in the context of risk measures.
This section proposes the Wasserstein Barycenter Risk Measures at a confidence level α. The research considers risk measures such as Value at Risk and Conditional Value at Risk for a loss random variable defined by
Here X 1 , ..., X N are loss random variables attributed to risk types endowed with positive weights ω 1 , ...ω N (such that N j=1 ω j = 1) and over a fixed time period T. Now, for computation of V aR α (X + ), a joint law for the random vector (X 1 , ..., X N ) ′ is required. The Wasserstein barycenter can be regarded as the aggregate model for certain set of probability measures. It is also suitable for reaching an "average" distribution. The procedure also considers and optimal selection for the positive weights. They are connected with the source credibility for every prior. Moreover, the weights must be chosen equal when all priors remains acceptable. The equality also holds under unknowing performance reliability of the competing laws.
We are in position for definition of the Wasserstein Barycenter Value-at-Risk.
Definition 3.1. Given the aggregate position X + , a set of measures M = (µ 1 , ..., µ N ), a set of weights λ = (λ 1 , ..., λ N ) ∈ R N −1 and a set of quantiles F = (F −1 µ 1 , ..., F −1 µ N ) with α ∈ (0, 1). The Wasserstein Barycenter Value-at-Risk is defined as:
Next we use the notable property that the barycenter of distributions of location and scale families belongs to the same class. This allows to derive a closed-form formulas for the Wasserstein Barycenter risk measures for location and scale distributions.
Theorem 3.1. Let µ 1 , ..., µ N be location and scale measures with corresponding means m 1 , ..., m N and standard deviations σ 1 , ..., σ N ; then the Wasserstein Barycenter Value-at-Risk (V aR α (X + , λ)) is given by
Proof. It follows straightforwardly from 3.1 and 2.3.
The Wasserstein Barycenter Conditional Value-at-Risk is established next:
Theorem 3.2. Let µ 1 , ..., µ N be location and scale measures with corresponding means m 1 , ..., m N and standard deviations σ 1 , ..., σ N , then the Wasserstein Barycenter Conditional Value-at-Risk (T CV aR α (X + , λ)) is given by
where Z = xq−m λ σ λ , g Z (.) and G Z (.) are the density and cumulative distribution functions of the standard random variable,
Proof. Note that
We now illustrate Theorems 3.1 and 3.2 under the Gaussian distribution.
Normal Distribution: Let µ 1 , ..., µ N be Normal measures with corresponding means m 1 , ..., m N and standard deviations σ 1 , ..., σ N ; then the V aR α (X + , λ) and the T CV aR α (X + , λ) are given by
here φ stands for the standard Gaussian distribution and Φ −1 is the inverse of the standard Gaussian distribution,
Randomness in financial markets has promoted important research about robust measures of market risk. This problem motivates a profuse study about market risk. An issue involving the risk of loss for an investment under multifactor movements in a market. Some dynamical risk factors consider the interest and exchange rates, commodity risks and capital, among others. Thus, this section focus on estimation of the aggregation VaR for a risk portfolio ruled by Nasdaq and S&P500 stock indices. The Nasdaq log-returns and the S&P500 log-returns will be denoted as X 1 and X 2 , respectively. In this case, the portfolio log-return, X + , has the form X + = λ 1 X 1 + λ 2 X 2 . Here λ = (λ 1 , λ 2 ) and λ 1 and λ 2 are the portfolio weighs of the assets 1 and 2, X 1 and X 2 . Without loss of generality, a portfolio under equal weights, in both indices, is considered. However, it is not a strict restriction and they can change freely. Finally, for the marginal returns a normal distribution is proposed and a one-day period VaR will be considered.
Next, the Wasserstein Barycenter VaR is computed by using 6. In this case, each stock is ruled by Normal distribution. We follow the 2972 daily closing prices given by [Landsman et al.(2003) According to Table 1 , the return series distributions of Nasdaq and S&P500 have small asymmetry, but strong kurtosis, in particular the first one. Note also that both series present positive means (annualized).
Thus, the results of Section 3 can be used for estimation of the the aggregation Valueat-Risk by using the equation
where Φ −1 holds for the inverse of the standard Gaussian distribution, and the mean and the standard deviation are computed via m λ = N j=1 λ j m j , and σ λ = N j=1 λ j σ j , respectively.
The model computation includes both "unfiltered" and "filtered" forms. The filtered model case has also considered the volatility changes of the instrument. Such model will be referred as Normal*. In the unfiltered Normal VaR (Normal), all the σ j 's, j = 1, , , N receive the same value of the sample standard deviation. But in the filtered Normal VaR (Normal*) the σ j 's are estimated with a Exponentially Weighted Moving Average (Ewma) model, where σ t = (1 − ζ)x 2 i + ζσ 2 t−1 . Kupiec test evaluates the performance by computing the exceptions number in the corresponding test period. In this case, H 0 : represents the null hypothesis and 1 − α is the probability of an exception occurrence. If m is the number of observations for the test period and x denotes the expected frequency of exceptions, then h = x m is the difference between the observed frequency of losses and VaR. The test statistics,
The corresponding test statistics is given by,
.
It rejects the null hypothesis, with a 95% confidence level, for LR > χ 2 (1). In that case, the VaR estimations are not statistical meaningful generated by the particular VaR model, see [McNeil (1999) ].
The dataset under consideration is now divided into sample and test periods, with a selected window of 750 observations; and since there are 2971 observations available, then 2220 VaR tests can be performed at each level can be performed. The corresponding results are presented in Table 2 . For all the α levels, the Wasserstein Barycenter-N* model showed the best performance in the VaR estimation. Moreover, for α = 0.1, 0.05, the Wasserstein Barycenter-N model also provided a satisfactory behavior. In terms of the Kupiec test, applied to the number of exceptions for the Wasserstein Barycenter-N * model, the null hypothesis were not rejected for all the α levels under consideration. In particular, high p values of 0.2837, 0.9223, 0.8653 and 0.1668 were obtained for Wasserstein Barycenter-N * model, and 0.8323 and 0.07130 for the Wasserstein Barycenter-N model at α = 0.1, 0.05 levels.
According to [Jianping et al. (2012) ], most of the banks consider simple methods for assessing risk. Complex methods involving copulas, simulation, hybrids models and advanced probability theory are implemented by few large banks.
Next apart will provide a summary of the variance-covariance and simple summation methods. Then a performance comparison with the proposed model is given. At the end, the new method is also confronted with the sophisticated hybrid copulas methodology given by [Landsman et al.(2003) ].
In the referred context of [Jianping et al. (2012) ], a research by the IFRI and CRO Forum's24 showed that 60% or more of the studied banks consider simple approaches as the var-cov method in order to aggregate risk. And at least a little more advanced approaches (such as the supported methods by simulation) were used by only 20% or less of the financial institutions in the survey. In fact, the summation setting prevails in most of the banks, a general model including variance-covariance method and simple summation. In the latest case the risks can be seen as explicitly weighted, meanwhile for the simple integration, the risks appear implicitly weighted.
This section computes the total VaR aggregated under several methods and different confidence levels. For correctness of our results different performances with published results were studied. The methodology starts with a series of univariate portfolio returns; then by the use of simple sum and variance-covariance methods, the aggregate VaR is estimated. These results are briefly described in the next lines.
Simple Summation method: The integration of N risks is intuitively reached by aggregating the risks under summation of the particular V aR α (X i ) of each risk X i , i = 1, . . . , N. Then the total aggregated VaR, V aR α (X + ), is expressed as:
see [Embrechts et al. (2013) ] and [Li et al. (2015) ].
Now, as usual, the Gaussian model supports several methods in probability and statistics studies. In particular, the VaR of multivariate Gaussian laws are a common parametric method for multivariate VaR models. The technique supposes a multivariate Gaussian distribution (with mean µ and covariance Σ ) for the returns of the components in the portfolio. The method is characterized as follows:
Variance-Covariance method: If σ + = √ λΣλ ′ . and µ + = λµ are respectively the deviation and expected portfolio return, then, the estimation of the Value at Risk for the corresponding multivariate Gaussian distribution returns is given by
as before, Φ −1 represents the inverse of the standard Gaussian distribution. The covariance matrix and the mean vector in the Var-Cov approach are frequently unknown, then the model requires extra estimates taken from the current observations, see [Jianping et al. (2012) ] and [Li et al. (2015) ]. (38) A summary of the results are given next:
• According Kupiec test, the basic multivariate VaR approaches does not predict future losses properly. Exception number in the test period is small under the Simple Summation; a fact explaining future loss overestimation. In contrast, the exception number under the variance-covariance method is large; promoting an underestimation of the future losses.
• The Wasserstein Barycenter-N* and Wasserstein Barycenter-N model exhibits a remarkable performance for future losses predictions. Moreover, Wasserstein Barycenter methods are stronger respect the other VaR models, because in the same reference time they provide a small exception probability and then a high-level capital reserve is not required.
In the set of the analyzed methods, the VaR Forecasting at all confidence levels were achieved in high performance by the proposed Wasserstein Barycenter-N* model. In fact, the Wasserstein Barycenter-N exhibited a better VaR forecasting than Var-Covar and Simple Summation methods in all the confidence levels. Empirical results also demonstrated the known fact that the Simple Approach provides an upper bound for the true VaR. In particular, for a confidence level of 0.1%, a VaR of 0.0387 derived by Wasserstein Barycenter-N* was the third part of the value 0.0978 based on the Simple Summation method. In such comparison context, our approach provides a number of possibilities for a wide bank class. Thus, under a conservative Wasserstein Barycenter VaR compared with the general average, the proposed method, indexed by different types of weights, allow several criteria in order to improve the results. In contrast, the Var-covar method is preferably optimistic. Finally, the notable closed property for the barycenter in the class of location-scale distribution, just require a profuse knowledge of the prior barycenter under the selected distribution. Then, in the Gaussian case, an exact formula for the value at risk Wasserstein Barycenter can be derived and applied into the complete reference class of distributions. This opens an interesting perspective for manage risk studies of series based on non Gaussian models, because a complete mathematical description of the Wasserstein Barycenter VaR can be found under the selected prior distribution.
Robust multivariate methods including Copulas and ARCH models also involve VaR estimation. However, under a big number of assets, such models produces biased parameter estimations and demand a high computational cost. For completeness we contrast our methods with those derived by [Jianping et al. (2012) ]. (32) Note that the Wasserstein Barycenter-N* model gives a satisfactory VaR estimation for α = 0.1, α = 0.05 and α = 0.01. When α = 0.01 the estimation equals the result of the hibrid SJC Copula + GARCH-E method; for α = 0.005 it is near to the hibrid SJC Copula + GARCH-E model, which was proposed by [Jianping et al. (2012) ] as the best approach. In fact, note that SJC Copula + GARCH-E method succeeds in the VaR forecasting at a 99.5% confidence level, but it fails for 90% and 95%. Now, kupiec test given in Table 2 highlights that our methods show a high performance in all the confidence levels. Moreover, inclusion of GARCH models should improved our methods, but for space reasons, we will leave that study for a future work. Given the robustness of GARCH models, the VaR estimation will require a few violentios and then a minor capital reserve should be demand.
We end this section by showing the behavior of our method in a crucial modern problem. Explicitly, we will research the impact of COVID-19 in the 2020 Stock Market Crashs measured in the Wall Street indexes of NASDAQ Composite, S&P 500 and Dow Jones Industrial Average. These indexes reported historical loss levels, only registered in the financial crisis of 2008.
As usual, the complete data set is split in sample and test periods. Sample data ranges from January 4th, 2010 to December 31st, 2018. The 2264 daily return for every each stock index also refers as the required historical data for a plausible VaR estimation. The test data, for detecting the performance of the value at risk models, ranges from January 2nd, 2019 to July 31st, 2020. The VaR is estimated for each day in the test period (399 days), with the information offered by the 2663 observations ahead of it. Finally, the performance of the VaR model is measured by a comparison of the current loss and the estimated value at risk. The division for the sample and test periods are summarized in Table 5 .
Basic descriptive statistics analysis for the interval 2010-2019 and the 2020 year are given in Tables 6 and 7, respectively.
Note that 2020 exhibits extremal high volatilities, in fact all the indexes highlight a notorious difference between that year and the interval 2010-2019. In addition, the minimum returns of the three indicators were reached on March 16th, 2020.
Next we compare the forecasting performance of the Wasserstein Barycenter-N* method with the basic multivariate VaR models. The best challenge for both methods resides in the financial market behavior of the 2020 year and its high volatility. First, Figure 1 shows the forecast of the next trading day VaR for 2019 and 2020 by using both methods. Under a one day holding period, the models were computed for a 99% confidence level. As it can be seen the basic multivariate VaR models provide poor estimates for future losses. The VaR estimates rises improperly VaR estimates and unacceptable test period exceptions. However, the VaR estimation given by the Wasserstein Barycenter-N* VaR model is highly accurate. The new approach fits the volatile movements of the returns and predicts future losses notably, in comparison with the basic multivariate VaR approach. Moreover, according to Figure 1 the variance-covariance approach cannot follow the strong volatility since February 2020; and always presents underestimation. In terms of Simple Summation model, the test period shows lower violations but a strong conservative characteristic is noted. Finally, our proposed methods reach satisfactory VaR forecast in normal periods and extremal periods for high volatility of the periods perform The results indicate that the proposed approach provides satisfactory forecasts of VaR not only for the "normal periods" but also for the periods of high volatility due to the COVID-19 stock market influence.
This work has proposed a new multivariate risk measure model based on the Wasserstein barycenters of probability measures under a location and scatter class. The method was compared with basic and sophisticated models under a portfolio characterized by S&P500 and Nasdaq stock indices. The new model was also compared in United States market indices of high volatility during the current COVID-19. Kupiec test was used for assessing the performance of the existing and new methods.
The new approach is based on a notable property: Wasserstein barycenters of measures supported on location and scatter family belong to the same class. Then the paper proposed exact formulae for the Wasserstein Barycenter Value at Risk and the Wasserstein Barycenter Conditional Value at Risk for the addressed location and scatter family. This promotes a new setting for building robust risk measures and aggregation of different VaR in multiple financial markets. The closed form formulae are easily programmed for applications.
The technique can be used for non Gaussian distributions, opening an interesting future research for more robust risk measures.