Avoiding zero probability events when computing Value at Risk contributions: a Malliavin calculus approach
AAvoiding zero probability events when computing Value at Riskcontributions: a Malliavin calculus approach
Yuri F. Saporito · Rodrigo S. Targino
April 29, 2020
Abstract
This paper is concerned with the process of risk allocation for a generic multivariate modelwhen the risk measure is chosen as the Value-at-Risk (VaR). Making use of Malliavin calculus, we recastthe traditional Euler contributions from an expectation conditional to an event of zero probability to aratio of conditional expectations, where both the numerator and the denominator’s conditioning eventshave positive probability. For several different models we show empirically that the estimator using thisnovel representation has no perceivable bias and variance smaller than a standard estimator used inpractice.
Keywords
Risk Management · Capital Allocation · Malliavin Calculus
Let us consider X = ( X , ..., X d ) the losses (negative of the returns) of d different assets in a portfolio. Fora linear portfolio with unitary exposure to each asset, the portfolio-wide loss is defined as X = (cid:80) di =1 X i .After a risk measure of the portfolio is computed, one is usually interested in understanding how mucheach asset contributes to the overall portfolio risk, in a process known as risk allocation.Here we follow the Euler allocation principle, proposed in [21]. This principle arises in different contextsin the literature. For example, in [4] and [11] it is motivated by two (different) sets of axioms, leadingto coherent allocation principles (for a relationship between coherent risk measures and coherent capitalallocations see [2]). It has also found a wide range of applications, ranging from credit ([21], [8]) andsystemic risks ([1], [14]) to banking and insurance capital [22].In this paper we assume the risk measure of interest is the Value-at-Risk (VaR) with confidence level α . For this measure, the Euler principle states the VaR allocation to the i -th asset is given by C αi = E [ X i | X = VaR α ( X )] . (1.1)The expectations in (1.1) rarely have close form solutions; with a notable exception being whenthe losses are jointly Gaussian or, more generally, elliptically distributed [16]. Therefore, approximationstrategies need to be used if one is willing to use a generic multivariate loss model.As the VaR allocations are expectations conditional to an event of probability zero, direct MonteCarlo simulation cannot be used. One strategy, followed by [8], to avoid the zero probability conditionin the expectations is to generate samples from p ( x | x ∈ [VaR α − δ ( X ) , VaR α + δ ( X )]), which leads toapproximate allocations, henceforth named δ -allocations, C α,δi = E [ X i | X ∈ [VaR α − δ ( X ) , VaR α + δ ( X )]] . (1.2)To the best of our knowledge, the issue of computing the exact VaR allocations has only been addressedin [13] and its extension in [12]. In these papers the authors use a Metropolis Hastings algorithm to samplefrom p (( x , . . . , x d − ) | x = VaR α ( X )) allied to the fact that C αd = VaR α ( X ) − (cid:80) d − i =1 C αi , in a strategythat resembles the slice sampler employed for the computation of the ES allocations of [20].Here we take a different route. Instead of directly targeting the distribution of X | X = VaR α ( X ) webuild upon the ideas of [7] and develop a novel expression for VaR allocations, namely, C αi = E [ X i π i | X ≥ VaR α ( X )] E [ π i | X ≥ VaR α ( X )] , (1.3) School of Applied Mathematics, Getulio Vargas Foundation, Brazil Bold letters always denote vectors and we use the notation = (0 , . . . , a r X i v : . [ q -f i n . C P ] A p r Yuri F. Saporito, Rodrigo S. Targino with weights π i depending on the distribution of X and precisely described in Theorem 2.2.The reader should note that, besides the fact the expectations in (1.3), differently from (1.1), donot depend on a zero probability event, the conditioning events are precisely the same as the one inthe Expected Shortfall (ES) allocations derived from the Euler principle. Therefore, the computationallyefficient algorithms derived in [20], [19], [13] and [12] can be easily adapted to estimate the expectationsin (1.3).The remainder of the paper is organized as follows. In Section 2 we state the main result of the paper,i.e. the new expression of the VaR contributions using only expectations conditional to events of positiveprobability. In Section 3 we provide examples of VaR allocations for several multivariate models and inSection 4 we consider some simulation studies to assess the efficiency of the proposed expression for theVaR allocations. Conclusions are provided in Section 5. In this section we present and discuss the main result of this paper, which is proved in Appendix A.2.For a very short overview of Malliavin calculus with emphasis on the necessary results for our work, weforward the reader to Appendix A.1. We start with a basic definition.
Definition 2.1
For a given g ∈ C ( B ; R d ), where B ⊂ R k is closed, we define the operator A g,i : B −→ R × k as A g,i ( u ) = ∇ g i ( u ) (cid:88) j (cid:54) = i ∇ g j ( u ) , (2.1)where g ( u ) = ( g ( u ) , . . . , g d ( u )). Theorem 2.2
Assume X = g ( U ) , where U = ( U , . . . , U k ) ∼ U [0 , k and g ∈ C ([0 , k ; R d ) . Moreover,assume there exists, for each i = 1 , . . . , d , f i ∈ C ([0 , k ; R k ) such that, for any u ∈ [0 , k , A g,i ( u ) f i ( u ) = (cid:20) (cid:21) . (2.2) Then, the marginal risk allocation for the VaR risk measure is given by C αi = E [ X i π i | X ≥ VaR α ( X )] E [ π i | X ≥ VaR α ( X )] , (2.3) where the weight π i is π i = Tr ( ∇ f i ( U )) , (2.4) T r ( A ) is the trace operator of a matrix A and ∇ f i is the Jacobian matrix of f i ∈ C ([0 , k ; R k ) . We refer to (2.3) as the Malliavin formulation of the allocations. Several examples of models (i.e. g )and their correspondent weight π are studied in the following section and summarized in Table 2.1.It is possible to avoid Condition (2.2). However, this would lead to a more complex formula in detrimentof the simplicity of formula (2.3). For this alternative formulation of our result, see Appendix A.3.Nonetheless, we use the formulation of Theorem 2.2 due to the simplicity of the final formula for C αi allowing a more direct numerical implementation. Additionally, finding the function f i that solves (2.2)was feasible in realistic examples. Remark 2.3
Notice that if π = C ˜ π , for some known, constant C , then Equation (2.3) is valid with ˜ π . Remark 2.4
There are several approaches available to reduce variance of Malliavin estimators of the form(2.3). See, for instance, [7].
In this section we consider different choices of g , which imply different distributions for X . In the examplesbelow, given the model g , we find a function f that solves Equation (2.2) and then derive the weight π . voiding zero probability for VaR contributions: a Malliavin calculus approach 3 Distribution Marginal Weight π i Independent Log-Normal LN (0 , σ j ) (cid:88) j (cid:54) = i Z j + σ j σ j X j , where Z j = Φ − ( U j )Independent Exponential Exp ( λ j ) (cid:88) j (cid:54) = i (1 − X j )Independent Gamma Gamma ( α j , β j ) (cid:88) j (cid:54) = i (cid:18) α j − U j − β j (cid:19) p j ( U j )Gaussian N ( µ, LL T ) f i · Z , where (cid:20) (cid:96) i L − (cid:96) i (cid:21) f i = (cid:20) (cid:21) and L = (cid:80) dj =1 (cid:96) j Archimedean Copula ψ V ∼ LS − ( ψ ) (cid:88) j (cid:54) = i,k p j ( X j ) (cid:32) V ψ (cid:48) ( φ j ( U )) + ψ (cid:48)(cid:48) ( φ j ( U )) ψ (cid:48) ( φ j ( U )) (cid:33) − p (cid:48) j ( X j ) p j ( X j )Clayton Copula ψ ( t ) = (1 + t ) − /ϑ V ∼ Γ (1 /ϑ, (cid:88) j (cid:54) = i,k p j ( X j ) ψ ( φ j ( U )) ( − ϑ ( V ) − log U j ) + ϑ + 1) − p (cid:48) j ( X j ) p j ( X j )Gumbel Copula ψ ( t ) = e − t /ϑ V ∼ S ( ϑ , , c,
0; 1) c = (cos(0 . π/ϑ )) ϑ (cid:88) j (cid:54) = i,k p j ( X j ) ψ ( φ j ( U )) (cid:18) − ϑ V φ j ( U ) − ϑ + ( ϑ − φ j ( U ) − ϑ + 1 (cid:19) − p (cid:48) j ( X j ) p j ( X j ) Table 2.1
Table with the weights π for several distributions. Remember p j denotes the marginal density. Let us consider the case of independent random variable, X j = ϕ j ( U j ), j = 1 , . . . , d , where U , . . . , U d are iid U [0 ,
1] and ϕ j is twice differentiable. This implies that g i ( u ) = ϕ i ( u i ) and ∇ g i ( u ) = (0 , . . . , , ϕ (cid:48) i ( u i ) , , . . . , , where the non-zero entry is in the i -th position. Hence, we consider f i ( u ) = 1 d − (cid:18) ϕ (cid:48) ( u ) , . . . , ϕ (cid:48) i − ( u i − ) , , ϕ (cid:48) i +1 ( u i +1 ) , . . . , ϕ (cid:48) d ( u d ) (cid:19) to satisfy (2.2), which implies Tr( ∇ f i ( u )) = − d − (cid:88) j (cid:54) = i ϕ (cid:48)(cid:48) j ( u j )( ϕ (cid:48) j ( u j )) and, by Remark 2.3, we take π i = (cid:88) j (cid:54) = i ϕ (cid:48)(cid:48) j ( U j )( ϕ (cid:48) j ( U j )) . (3.1)Frequently, ϕ j is the cumulative distribution function (c.d.f.) of a given distribution with density p j . Ifthis density is differentiable, then π i = (cid:88) j (cid:54) = i p (cid:48) j ( U j )( p j ( U j )) . (3.2)For instance, let us consider some particular distributions:1. Log-normal: ϕ j ( u ) = e µ j + σ j Φ − ( u ) . This implies π i = (cid:88) j (cid:54) = i Z j + σ j σ j X j , where Z j = Φ − ( U j ). Yuri F. Saporito, Rodrigo S. Targino
2. Exponential: ϕ j ( u ) = 1 − e − λ j u , which gives π i = (cid:88) j (cid:54) = i (1 − X j ) .
3. Gamma: p j ( u ) = β αjj Γ ( α j ) u α j − e − β j u , which gives π i = (cid:88) j (cid:54) = i (cid:18) α j − U j − β j (cid:19) p j ( U j ) . Assume X is distributed according a multivariate Gaussian distribution, i.e. X = µ + L Z , where Z is a k -dimensional standard Gaussian random vector with k ≥ d and L is a d × k full-rank, lower triangularmatrix. In this model, it is advantageous to use the formulation of Theorem A.3. Hence, g i ( z ) = µ i + (cid:96) i · z and then ∇ g i ( z ) = (cid:96) i , where (cid:96) i is the i -th row of the matrix L . This implies that A g,i ( z ) is independent of z and A g,i ( z ) = A g,i = (cid:96) i (cid:88) j (cid:54) = i (cid:96) j = (cid:96) i L − (cid:96) i , where L = (cid:80) dj =1 (cid:96) j . Thus, we might consider f i ∈ R k constant such that A g,i f i = (cid:20) (cid:21) . (3.3)to satisfy (A.6). Since the matrix L is full rank, the rank of A g,i is 2. Therefore, there exist infinitelymany solution to A g,i x = (cid:20) (cid:21) and we could take any of them to be f i . Since f i is constant, ∇ f i = 0,which implies π i = f i · Z , (3.4)see Equation (A.7). An interesting generalization of the previous example is the class of elliptical distributions: X = µ + RL S , where S is a d -dimensional random vector uniformly distributed in the the sphere in R d , L is a d × d full-rank, lower triangular matrix and R is an one-dimensional radial random variable independent of S .If Z is a d -dimensional standard Gaussian random vector, then we may write S = Z (cid:107) Z (cid:107) . Moreover, without loss of generality, we may consider R = φ ( Z d +1 ) for some positive function φ , where Z d +1 ∼ N (0 ,
1) is independent of Z . Therefore, as in the previous example, we use the formulation ofTheorem A.3. Notice g i ( z , z d +1 ) = µ i + φ ( z d +1 ) ψ ( z )( (cid:96) i · z ) , where z = ( z , . . . , z d ) and ψ ( z ) = (cid:107) z (cid:107) − . Hence ∇ g i ( z , z d +1 ) = (cid:2) ( (cid:96) i · z ) φ ( z d +1 ) ∇ ψ ( z ) + ψ ( z ) φ ( z d +1 ) (cid:96) i , ψ ( z ) φ (cid:48) ( z d +1 )( (cid:96) i · z ) (cid:3) , voiding zero probability for VaR contributions: a Malliavin calculus approach 5 where we are specifying the ( d + 1)-entry of the vector. If we define Ψ i ( z ) = ( (cid:96) i · z ) ∇ ψ ( z ) + ψ ( z ) (cid:96) i = − ( (cid:96) i · z ) ψ ( z ) z + ψ ( z ) (cid:96) i , we find ∇ g i ( z , z d +1 ) = (cid:2) φ ( z d +1 ) Ψ i ( z ) , ψ ( z ) φ (cid:48) ( z d +1 )( (cid:96) i · z ) (cid:3) . Then, slightly abusing the notation, we write f i ( z , z d +1 ) = ( f ( z , z d +1 ) , f d +1 ( z , z d +1 )).Hence, Equation (2.2) becomes φ ( z d +1 )( Ψ i ( z ) · f ( z , z d +1 )) + ψ ( z ) φ (cid:48) ( z d +1 )( (cid:96) i · z ) f d +1 ( z , z d +1 ) = 0 , (cid:88) j (cid:54) = i ( φ ( z d +1 ) Ψ j ( z ) · f ( z , z d +1 ) + ψ ( z ) φ (cid:48) ( z d +1 )( (cid:96) j · z ) f d +1 ( z , z d +1 )) = 1 . Assuming that f ( z , z d +1 ) was already computed, we may write f d +1 using the first equation above: f d +1 ( z , z d +1 ) = − φ ( z d +1 )( Ψ i ( z ) · f ( z , z d +1 )) ψ ( z ) φ (cid:48) ( z d +1 )( (cid:96) i · z ) . (3.5)Using this formula in the second equation, we find1 = (cid:88) j (cid:54) = i ( φ ( z d +1 ) Ψ j ( z ) · f ( z , z d +1 ) + ψ ( z ) φ (cid:48) ( z d +1 )( (cid:96) j · z ) f d +1 ( z , z d +1 ))= (cid:88) j (cid:54) = i (cid:18) φ ( z d +1 ) Ψ j ( z ) · f ( z , z d +1 ) − ψ ( z ) φ (cid:48) ( z d +1 )( (cid:96) j · z ) φ ( z d +1 )( Ψ i ( z ) · f ( z , z d +1 )) ψ ( z ) φ (cid:48) ( z d +1 )( (cid:96) i · z ) (cid:19) = (cid:88) j (cid:54) = i (cid:18) φ ( z d +1 ) Ψ j ( z ) · f ( z , z d +1 ) − (cid:96) j · z (cid:96) i · z φ ( z d +1 )( Ψ i ( z ) · f ( z , z d +1 )) (cid:19) = φ ( z d +1 ) (cid:88) j (cid:54) = i (cid:18) Ψ j ( z ) · f ( z , z d +1 ) − (cid:96) j · z (cid:96) i · z ( Ψ i ( z ) · f ( z , z d +1 )) (cid:19) = φ ( z d +1 ) (cid:88) j (cid:54) = i Ψ j ( z ) − (cid:96) j · z (cid:96) i · z Ψ i ( z ) · f ( z , z d +1 )= φ ( z d +1 ) ψ ( z ) (cid:18) L − L · z (cid:96) i · z (cid:96) i (cid:19) · f ( z , z d +1 ) = φ ( z d +1 ) (cid:107) z (cid:107) (cid:18) L − L · z (cid:96) i · z (cid:96) i (cid:19) · f ( z , z d +1 ) , where, as before, L = (cid:80) j (cid:96) j . Therefore, if we define E ( z ) = L − L · z (cid:96) i · z (cid:96) i , we may choose f ( z , z d +1 ) = (cid:107) z (cid:107) φ ( z d +1 ) E ( z ) (cid:107) E ( z ) (cid:107) . Plugging this back into the formula for f d +1 in Equation (3.5), we find f d +1 ( z , z d +1 ) = −(cid:107) z (cid:107) Ψ i ( z ) · E ( z ) φ (cid:48) ( z d +1 )( (cid:96) i · z ) (cid:107) E ( z ) (cid:107) = − (cid:107) z (cid:107) φ (cid:48) ( z d +1 )( (cid:96) i · z ) (cid:107) E ( z ) (cid:107) (cid:18) ( L · (cid:96) i ) − L · z (cid:96) i · z (cid:107) (cid:96) i (cid:107) (cid:19) . The weight could then be computed by performing the calculation of Equation (A.7).
Yuri F. Saporito, Rodrigo S. Targino
Let us now assume we have a random vector X = ( X , . . . , X d ) such that its joint c.d.f. given by F ( x , . . . , x d ) = C ψ ( F ( x ) , . . . , F d ( x d )) , (3.6)where C ψ is an Archimedean copula with generator ψ . Denoting by F = LS − ( ψ ) the inverse Laplace-Stieltjes transform of ψ , we assume that both F and the marginals F i are absolutely continuous and thatthe density of F i , denoted by p i , is differentiable. Note that these regularity conditions may be relaxedusing standard mollification techniques.Following the Marshall-Olkin algorithm [15], in order to generate a sample ( U , . . . , U d ) from thecopula C ψ , one should(a) Sample V ∼ F = LS − ( ψ );(b) Sample U i iid ∼ U [0 , i = 1 , . . . , d ;(c) Define U i = ψ ( − log( U i ) / V ), i = 1 , . . . , d .Then, a sample X = ( X , . . . , X d ) from the joint distribution in (3.6) is such that X i = F − i ( U i ).In order to write each X i as a function of i.i.d. uniform random variables, we first note that theMarshall-Olkin algorithm needs k = d + 1 i.i.d. uniforms: U , . . . , U k − , from item (b) and U k , which issuch that V = H ( U k ) with H = F − . Therefore, we write the random variable X i as X i = F − i (cid:0) U i ( U ) (cid:1) , where U = ( U , . . . , U k ) and the function U i computes the pseudo-observations: U i ( u ) = ψ (cid:18) − log( u i ) H ( u k ) (cid:19) . The reader should note there is an abuse of notation here, as we use U i for both the function definedabove and the random variable from the Marshall-Olkin sampling procedure.If we further define H i = F − i , then we have that g i ( u ) = H i ( U i ( u )) and also that A g,i ( u ) = H (cid:48) i ( U i ( u )) ∇U i ( u ) (cid:88) j (cid:54) = i H (cid:48) j ( U j ( u )) ∇U j ( u ) . (3.7)The gradient of U i can easily be computed as ∂ U i ∂u j ( u ) = , if j (cid:54) = i and j (cid:54) = k, − ψ (cid:48) (cid:18) − log( u i ) H ( u k ) (cid:19) u i H ( u k ) , if j = i,ψ (cid:48) (cid:18) − log( u i ) H ( u k ) (cid:19) log( u i ) H (cid:48) ( u k ) H ( u k ) , if j = k. In a more succinct notation, for i = 1 , . . . , k −
1, define φ i ( u ) = − log( u i ) H ( u k )and ∇U i ( u ) = ψ (cid:48) ( φ i ( u )) 1 H ( u k ) (cid:18) − u i e i − φ i ( u ) H (cid:48) ( u k ) e k (cid:19) , with e i = (0 , . . . , , , , . . . , ∈ R k is the vector with zeros everywhere but at the i -th coordinate.In order to find a a function f i ( u ) satisfying Equation (2.2), we write f i ( u ) = (cid:80) k(cid:96) =1 α (cid:96) ( u ) e (cid:96) . We pick α i ( u ) = 0 ,α k ( u ) = 0 , voiding zero probability for VaR contributions: a Malliavin calculus approach 7 but other choices could be considered. We then choose α j such that (cid:88) j (cid:54) = i,k H (cid:48) j ( U j ( u )) ψ (cid:48) ( φ j ( u )) 1 H ( u k ) (cid:18) − α j ( u ) u j (cid:19) = 1 . One solution is α j ( u ) = − u j k − H ( u k ) ψ (cid:48) ( φ j ( u )) 1 H (cid:48) j ( U j ( u )) . Then, the weight π could be written as π i = (cid:88) j (cid:54) = i,k ∂α j ∂u j ( U ) , where ( k − ∂α j ∂u j ( u ) = − H ( u k ) ψ (cid:48) ( φ j ( u )) 1 H (cid:48) j ( U j ( u )) − ψ (cid:48)(cid:48) ( φ j ( u )) ψ (cid:48) ( φ j ( u )) H (cid:48) j ( U j ( u )) − H (cid:48)(cid:48) j ( U j ( u ))( H (cid:48) j ( U j ( u ))) . Moreover, since H i = F − i and p i is the density of F i , we find H (cid:48) i ( x ) = 1 p i ( H i ( x )) and H (cid:48)(cid:48) i ( x ) = − p (cid:48) i ( H i ( x )) p i ( H i ( x )) . Therefore, using the fact H i ( U i ( U )) = X i , we conclude( k − ∂α j ∂u j ( U ) = − H ( U k ) ψ (cid:48) ( φ j ( U )) f j ( X j ) − ψ (cid:48)(cid:48) ( φ j ( U )) ψ (cid:48) ( φ j ( U )) p j ( X j ) + p (cid:48) j ( X j ) p j ( X j ) . Finally, we find π i = (cid:88) j (cid:54) = i,k p j ( X j ) (cid:32) H ( U k ) ψ (cid:48) ( φ j ( U )) + ψ (cid:48)(cid:48) ( φ j ( U )) ψ (cid:48) ( φ j ( U )) (cid:33) − p (cid:48) j ( X j ) p j ( X j ) . If we define γ j ( u ) = H ( u k ) ψ (cid:48) ( φ j ( u )) + ψ (cid:48)(cid:48) ( φ j ( u )) ψ (cid:48) ( φ j ( u )) , (3.8)which depends solely on the copula C ψ , we may rewrite π as π i = (cid:88) j (cid:54) = i,k p j ( X j ) γ j ( U ) − p (cid:48) j ( X j ) p j ( X j ) . (3.9)Moreover, if we denote c ( t ) = (log ψ ( t )) (cid:48) (i.e. ψ (cid:48) ( t ) = c ( t ) ψ ( t )), we can write γ j as γ j ( u ) = 1 ψ ( φ j ( u )) (cid:18) H ( u k ) c ( φ j ( u )) + c (cid:48) ( φ j ( u )) c ( φ j ( u )) + 1 (cid:19) . For ψ ( t ) = (1 + t ) − /ϑ , with ϑ ∈ (0 , ∞ ), one gets the Clayton family of copulas. Here, c ( t ) = − ϑ t and γ j ( u ) = 1 ψ ( φ j ( u )) ( − ϑ ( H ( u k ) − log u j ) + ϑ + 1) , where H ( U k ) ∼ Gamma (1 /ϑ, In this case, ψ ( t ) = e − t /ϑ , for ϑ ∈ [1 , + ∞ ). Hence, c ( t ) = − ϑ t /ϑ − and γ j ( u ) = 1 ψ ( φ j ( u )) (cid:16) − ϑ H ( u k ) φ j ( u ) − /ϑ + ( ϑ − φ j ( u ) − /ϑ + 1 (cid:17) , where H ( U k ) ∼ S (1 /ϑ, , (cid:0) cos( π ϑ ) (cid:1) ϑ ,
0; 1), Nolan’s 1-parametrization of stable distribution [17].
Yuri F. Saporito, Rodrigo S. Targino
Depending on the chosen marginals, we can simplify the term p (cid:48) j /p j . Indeed,1. Gaussian N (0 , σ j ): p (cid:48) j ( x ) p j ( x ) = − xσ j ;2. Generalized Pareto GP D ( ξ j , β j ): p (cid:48) j ( x ) p j ( x ) = − ξ j β j (cid:18) ξ j xβ j (cid:19) − , for x ≥ Exp ( λ j ): p (cid:48) j ( x ) p j ( x ) ≡ − λ j ;4. Log-normal LN ( µ j , σ j ): p (cid:48) j ( x ) p j ( x ) = − x (cid:32) log x − µ j σ j + 1 (cid:33) . When C : [0 , d → [0 ,
1] is a copula and
U ∼ C , its survival copula is defined as (cid:101) C ( u ) = C (1 − u , . . . , − u d ), where C ( u ) = P ( U > u , . . . , U d > u d ). In other words, the survival copula is the copula associatedwith the survival function of U . Therefore, to sample (cid:101) U ∼ (cid:101) C one can simply sample U ∼ C and define (cid:101) U i = 1 − U i , for i = 1 , . . . , d .If the original copula is Archimedean, a simple calculation shows that the weight for the correspondingsurvival copula is π i = (cid:88) j (cid:54) = i,k p j ( X j ) γ j ( U ) + p (cid:48) j ( X j ) p j ( X j ) , (3.10)where γ j is given by Equation (3.8). In order to assess the efficiency of the Malliavin representation (see Theorem 2.2) when compared tothe δ -allocations of equation (1.2), we perform simulation exercises for some of the models presented inSection 3.To compare simple Monte Carlo estimators (to be described soon), we first pre-compute the VaRand use it for both methods. The error in this step is irrelevant, as we can think we are computing E [ X i | X = B ] where B is a constant, so for the remainder of this section we assume all VaRs have beenperfectly computed.Let us denote by A α,δ the conditioning event in (1.2), i.e. A α,δ = [ X ∈ [VaR α − δ ( X ) , VaR α + δ ( X )]].We also define A α = [ X ≥ VaR α ( X )], the event in the Malliavin formulation (2.3).For a sample X (1) , . . . , X ( N ) , we denote by Y (1) α , . . . , Y ( N ∗ α ) α the N ∗ α samples which belong to A α and γ (1) α , . . . , γ ( N ∗ α ) α the corresponding π ’s from Theorem 2.2. Similarly, Y (1) α,δ , . . . , Y ( N ∗ α ) α,δ denotes the N ∗ α,δ samples in A α,δ . Note that 0 ≤ N ∗ α , N ∗ α,δ ≤ N , i.e., both N ∗ α and N ∗ α,δ can be zero and in this casethe estimates would not be defined. For a fixed sample size N , we have that E [ N ∗ α ] = (1 − α ) N and E [ N ∗ α,δ ] = 2 δN , and as long as 1 − α ≥ δ , we have E [ N ∗ α ] ≥ E [ N ∗ α,δ ].We wish to compare the following estimators of C αi : (cid:98) C α,δi = 1 N ∗ δ N ∗ α,δ (cid:88) n =1 Y ( n ) α,δ,i . (4.1)and (cid:98) C αi = (cid:80) N ∗ α n =1 Y ( n ) α,i γ ( n ) α,i (cid:80) N ∗ α n =1 γ ( n ) α,i . (4.2) voiding zero probability for VaR contributions: a Malliavin calculus approach 9 . . . . d = log ( N ) E x pe c t a t i on l l l l ll l l l ll l l l ll l l l ll l l l ll l l l l . . . . d = log ( N ) E x pe c t a t i on l l l l ll l l l ll l l l ll l l l ll l l l ll l l l l . . . . d = log ( N ) E x pe c t a t i on l l l l ll l l l ll l l l ll l l l ll l l l ll l l l l . . . . d = log ( N ) E x pe c t a t i on l l l l ll l l l ll l l l ll l l l ll l l l ll l l l l − − − − − d = log ( N ) l og ( V a r i an c e ) l l l l ll l l l ll l l l ll l l l ll l l l ll l l l l − − − − − d = log ( N ) l og ( V a r i an c e ) l l l l ll l l l ll l l l ll l l l ll l l l ll l l l l − − − − − d = log ( N ) l og ( V a r i an c e ) l l l l ll l l l ll l l l ll l l l ll l l l ll l l l l − − − − − d = log ( N ) l og ( V a r i an c e ) l l l l ll l l l ll l l l ll l l l ll l l l ll l l l l Fig. 4.1
Multivariate Gaussian model of Section 4.1. Mean (top row) and variance (bottom row) of the δ -estimator (black)and the Malliavin (red) for the first allocation C α . Each line type (solid, dashed and dotted) represent a different value of α . The columns represent different values for δ . In the sequel we study the impact of δ and the sample size N in the mean and variance of bothestimators (computed using 10 000 replicates). For each combination of δ ∈ { − , − , − , − } and N ∈ { , (cid:98) . (cid:99) , , (cid:98) . (cid:99) , } we calculate the allocations for α ∈ { . , . , . } . Althoughpractitioners would not be interested in α = 0 .
5, the method proposed should produce lower variance forany α and for higher values of α the estimates are numerically more stable.It should also be noticed that regardless of the value of α , the conditioning event in (1.2) has probability2 δ , while the events in (2.3), the Malliavin formulation for the VaR allocations, have probability 1 − α .Moreover, all comparisons are naturally adjusted by their computational cost, as both estimators use thesame number of samples. In fact, we use exactly the same samples in both methods.As the results are qualitatively similar across the components of X , for each model we only presentthe results for its first component, i.e. C α . If we assume the vector of risks follows a multivariate normal distribution, i.e., X = ( X , . . . , X d ) ∼ N ( , Σ ) then the VaR contributions can be computed in closed form, as C αi = Φ − ( α ) ( Σλ ) Ti √ λ T Σλ . We alsohave that VaR α ( X ) = Φ − ( α ) √ λ T Σλ . Using the notation from Section 3.2, we have µ = and Σ = LL T .In our example we take d = 3 and L = . . . . , which leads to variances of 1 . , .
74 and 2 .
85 for each one of the components of X and correlations rangingfrom 0 .
58 to 0 . α . When the value of δ is too small, we see the δ -estimator needs more samples to produce the same error, specially for higher quantiles.The striking result is related to the variance. As expected, when δ gets smaller the variance of theestimator of the δ -allocations (for the same sample size) depletes very quickly, since there are less samplesthat satisfies the condition in the expectation, if any. On the other hand, for the Malliavin method theproportion of samples satisfying the condition is on average 1 − α , whilst for the δ method is 2 δ . It can alsobe seen from Figure 4.1 that, although all variances converge exponentially fast to zero as N increases,the Malliavin variances are consistently smaller than its counterparts in the δ formulation. . . . . d = log ( N ) E x pe c t a t i on l l l l ll l l l ll l l l ll l l l ll l l l ll l l l l . . . . d = log ( N ) E x pe c t a t i on l l l l ll l l l ll l l l ll l l l ll l l l ll l l l l . . . . d = log ( N ) E x pe c t a t i on l l l l ll l l l ll l l l ll l l l ll l l l ll l l l l . . . . d = log ( N ) E x pe c t a t i on l l l l ll l l l ll l l l ll l l l ll l l l ll l l l l − − − − − d = log ( N ) l og ( V a r i an c e ) l l l l ll l l l ll l l l ll l l l ll l l l ll l l l l − − − − − d = log ( N ) l og ( V a r i an c e ) l l l l ll l l l ll l l l ll l l l ll l l l ll l l l l − − − − − d = log ( N ) l og ( V a r i an c e ) l l l l ll l l l ll l l l ll l l l ll l l l ll l l l l − − − − − d = log ( N ) l og ( V a r i an c e ) l l l l ll l l l ll l l l ll l l l ll l l l ll l l l l Fig. 4.2
Independent Log-Normals model of Section 4.2. Mean (top row) and variance (bottom row) of the δ -estimator(black) and the Malliavin (red) for the first allocation C α . Each line type (solid, dashed and dotted) represent a differentvalue of α . The columns represent different values for δ . For this example we follow the structure of Section 3.1 and assume X i ind ∼ LN (0 , σ i ) with σ = 0 . σ = 1and σ = 2. The results are presented in Figure 4.2.Although the results are qualitatively equivalent to the Gaussian case of the previous section, somepoints are worth stressing. First, there is no close form expression for the allocations in this example.Therefore, we can only believe the estimates have converged for large enough sample sizes. In this regard,the δ -estimates for δ = 10 − can be deemed unsatisfactory for N ≤ , as the (empirical) expectedvalues are consistently off the estimate for larger δ – which is consistent with the Malliavin formulation.The estimators’ variances behave exactly like in the Gaussian case, with the Malliavin formulationperforming uniformly better in all scenarios. Apart from that, for δ ≤ − we start to see that the rateof convergence of the variance decreases substantially when compared with the Malliavin formulation oreven with larger values of δ . This simulation exercise is related to Section 3.4.1, where the dependence structure of the multivariatevector X is assumed to be given by a Clayton copula. Here we assume the copula parameter is θ = 2 andthe marginals are X i ∼ N (0 , σ i ), with σ = 1, σ = 0 . σ = 1.The mean and variance of the estimators are presented in Figure 4.3. Once again the results arequalitatively equivalent to the others, with the additional feature that for δ ≤ − there are clearlyconvergence issues for the δ -estimator when the sample size N is smaller than 10 . . For the variance,the same comments made in Section 4.2 apply. For the last computational experiment, we use the model M1 from [12], where it is assumed X i ∼ GP D ( ξ i , β j ) with ξ i = 0 . β i = 1, for i = 1 , ,
3. The results are found in Figure 4.4.As with all other examples, we see the Malliavin allocations do not show any bias for large enoughsample size N and α ∈ { . , . } , even though for N = 10 some bias appear to be present. It is worthnoticing the apparent bias in the estimation for α = 0 .
5, which we cannot guarantee it is coming fromthe δ or Malliavin method. voiding zero probability for VaR contributions: a Malliavin calculus approach 11 . . . . d = log ( N ) E x pe c t a t i on l l l l ll l l l ll l l l ll l l l ll l l l ll l l l l . . . . d = log ( N ) E x pe c t a t i on l l l l ll l l l ll l l l ll l l l ll l l l ll l l l l . . . . d = log ( N ) E x pe c t a t i on l l l l ll l l l ll l l l ll l l l ll l l l ll l l l l . . . . d = log ( N ) E x pe c t a t i on l l l l ll l l l ll l l l ll l l l ll l l l ll l l l l − − − − d = log ( N ) l og ( V a r i an c e ) l l l l ll l l l ll l l l ll l l l ll l l l ll l l l l − − − − d = log ( N ) l og ( V a r i an c e ) l l l l ll l l l ll l l l ll l l l ll l l l ll l l l l − − − − d = log ( N ) l og ( V a r i an c e ) l l l l ll l l l ll l l l ll l l l ll l l l ll l l l l − − − − d = log ( N ) l og ( V a r i an c e ) l l l l ll l l l ll l l l ll l l l ll l l l ll l l l l Fig. 4.3
Clayton copula with Gaussian marginals of model in Section 4.3. Mean (top row) and variance (bottom row)of the δ -estimator (black) and the Malliavin (red) for the first allocation C α . Each line type (solid, dashed and dotted)represent a different value of α . The columns represent different values for δ . d = log ( N ) E x pe c t a t i on l l l l ll l l l ll l l l ll l l l ll l l l ll l l l l d = log ( N ) E x pe c t a t i on l l l l ll l l l ll l l l ll l l l ll l l l ll l l l l d = log ( N ) E x pe c t a t i on l l l l ll l l l ll l l l ll l l l ll l l l ll l l l l d = log ( N ) E x pe c t a t i on l l l l ll l l l ll l l l ll l l l ll l l l ll l l l l − − − − d = log ( N ) l og ( V a r i an c e ) l l l l ll l l l ll l l l ll l l l ll l l l ll l l l l − − − − d = log ( N ) l og ( V a r i an c e ) l l l l ll l l l ll l l l ll l l l ll l l l ll l l l l − − − − d = log ( N ) l og ( V a r i an c e ) l l l l ll l l l ll l l l ll l l l ll l l l ll l l l l − − − − d = log ( N ) l og ( V a r i an c e ) l l l l ll l l l ll l l l ll l l l ll l l l ll l l l l Fig. 4.4
Survival Clayton copula with GPD marginals. Mean (top row) and variance (bottom row) of the δ -estimator(black) and the Malliavin (red) for the first allocation C α . Each line type (solid, dashed and dotted) represent a differentvalue of α . The columns represent different values for δ . For sufficiently large values of δ ( δ = 10 − and, to some extent, δ = 10 ), the δ -estimator show lessvariance than the Malliavin. For δ = 10 the δ -estimator presents smaller variance for all quantiles α ,while for δ = 10 this is only valid for α = 0 .
99. For δ ≤ − and N ≥ . the estimator based on theMalliavin representation delivers a less variable estimator. As in the previous examples, for δ sufficientlysmall ( δ = 10 − ) the variance’s convergence of the δ -estimator when N increases is vastly compromised. We use results from the Malliavin calculus literature and are able to derive a novel expression for theValue-at-Risk contributions, where we go from an expectation conditional to a zero probability event in the usual representation, to a ratio of expectations conditional to events of positive probability. The newformulation is amenable to Monte Carlo simulation with mild hypothesis on the multivariate models andthe precise formulas are provided for a wide range of models.As the expectations in the proposed formulation resemble the Expected Shortfall allocations, onecould easily adapt the sampling schemes from, for example, [20] and [12] for further computational gains.Another possible avenue for further improvements of the proposed representation is explore betterestimators for the ratio of expectations. The methods in [10] have been tested but, although simple, didnot result in substantial bias reduction.
A Appendices
A.1 Malliavin Calculus
The main result, Theorem 2.2, was written under the general setting: X = g ( U ), where U = ( U , . . . , U k ) ∼ U [0 , k .In order to simplify the application of Malliavin calculus we deploy, we consider the following setting: X = g ( Z ), where Z = ( Z , . . . , Z k ) ∼ N ( , I k ), where I k is the identity matrix in R k . These formulations are obviously equivalent. InAppendix A.2, we first prove the result under the Gaussian formulation using techniques from Malliavin calculus and thenuse this result to prove Theorem 2.2.In this section, we set the main definitions and results on Malliavin calculus required for the results in this paper. Inshort, Malliavin calculus is a differential calculus for functionals of the Brownian motion. The goal of this appendix is tomake this paper as self-contained as possible. However, the theory of Malliavin calculus is very vast and therefore we referthe reader, for instance, to [18] for more details. Applications in Mathematical Finance and Insurance can be found in [6],[7] and [9].Let ( W t ) t ∈ [0 ,T ] be a k -dimensional Brownian motion, with W t = ( W t , . . . , W kt ), and ( F t ) t ∈ [0 ,T ] the filtration generatedby it. The space of random variables in L ( Ω, F T , P ) that are differentiable in the Malliavin sense is denoted by D , . Avery important subspace of D , is the space of smooth random variables, which is defined as F = g (cid:18)(cid:90) T h ( s ) dW s , . . . , (cid:90) T h n ( s ) dW s (cid:19) , with g ∈ C ( R n ) and h ∈ L ([0 , T ]; R k ). In this case, the Malliavin derivative at time t ≤ T , which is denoted by D t , isgiven by D t F = n (cid:88) k =1 ∂ k g (cid:18)(cid:90) T h ( s ) dW s , . . . , (cid:90) T h n ( s ) dW s (cid:19) h k ( t ) , where ∂ k g is the derivative of g with respect to the k th variable. An important case for our application is F = g ( W T , . . . , W kT ),where g ∈ C ( R k ). In this case, D t X = ∇ g ( W T , . . . , W kT ) . In the multivariate case where F = ( F , . . . , F m ), the Malliavin derivative D t F is a m × k matrix where the j th row isgiven by D t F j .Moreover, the adjoint operator of D , denoted by δ and called Skorokhod integral, is defined by the integration-by-partsformula: E [ F δ ( v )] = E (cid:20)(cid:90) T D t F · v t dt (cid:21) , ∀ F ∈ D , , where x · y is the canonical inner product in R k . The domain of δ is characterized by the R k -valued stochastic processes v = ( v t ) t ∈ [0 ,T ] (not necessarily adapted to the filtration ( F t ) t ∈ [0 ,T ] ) such that (cid:12)(cid:12)(cid:12)(cid:12) E (cid:20)(cid:90) T D t F · v t dt (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) ≤ C (cid:107) F (cid:107) , ∀ F ∈ D , , where C > v and (cid:107) F (cid:107) = E [ | F | ] / . The following result regarding the Skorokhod integrals will beuseful: for F j a smooth random variable and h j ∈ L ([0 , T ]; R k ), j = 1 , . . . , m , δ m (cid:88) j =1 F j h j = m (cid:88) j =1 (cid:18) F j (cid:90) T h j ( t ) dW t − (cid:90) T D t F j · h j ( t ) dt (cid:19) . (A.1)The following are the fundamental results for the theorem presented in this paper. The proof might be found in [5] and[7]. Moreover, applications of these results for pricing and hedging of financial products in Lvy and jump-diffusion settingscan be found in [3] Theorem A.1
Let
F, G ∈ D , such that F is R m -valued, G is R -valued with D t G non-degenerate. Assume there existsa process v in the domain of δ and E (cid:20) (cid:90) T D t G · v t dt (cid:12)(cid:12)(cid:12)(cid:12) F, G (cid:21) = 1 . (A.2)voiding zero probability for VaR contributions: a Malliavin calculus approach 13 Assume further that φ ∈ C ( R ) . Then E [ φ ( F ) | G = 0] = E (cid:20) φ ( F ) δ ( v ) H ( G ) − φ (cid:48) ( F ) H ( G ) (cid:90) T D t F v t dt (cid:21) E [ δ ( v ) H ( G )] , (A.3) where H ( x ) = 1 x ≥ is the Heaviside function. Theorem A.2
Additionally to the assumptions of the theorem above, assume E (cid:20) (cid:90) T D t F v t dt (cid:12)(cid:12)(cid:12)(cid:12) F, G (cid:21) = m , (A.4) where m is the m -dimensional zero vector. Then, for any Borel measurable function φ with at most linear growth atinfinity, E [ φ ( F ) | G = 0] = E [ φ ( F ) δ ( v ) H ( G )] E [ δ ( v ) H ( G )] . (A.5) A.2 Proof of Theorem 2.2
Because of the nature of the Malliavin calculus, we consider the Gaussian noise first and then prove the Uniform case basedon this result.
Theorem A.3
Assume X = g ( Z ) , where Z = ( Z , . . . , Z k ) ∼ N ( , I k ) and g ∈ C ( R k ; R d ) . Moreover, assume thereexists, for each i = 1 , . . . , d , f i ∈ C ( R k ; R k ) such that, for any z ∈ R k , A g,i ( z ) f i ( z ) = (cid:20) (cid:21) . (A.6) Then, the weight π i takes the form π i = π i ( Z ) = f i ( Z ) · Z − Tr ( ∇ f i ( Z )) (A.7) and the marginal risk allocation for the VaR risk measure can be computed as in (2.3).Proof Set T = 1 and define F = X i and G = X − VaR α ( X ), where X = (cid:80) dj =1 X j in Theorem A.2. Notice that both r.v’ssatisfy the conditions of Theorem A.2, as X i satisfy Condition A.6. Then, choosing v to satisfy Conditions (A.2) and (A.4),we find E [ X i | X = VaR α ( X )] = E (cid:104) X i δ ( v )1 X> VaR α ( X ) (cid:105) E (cid:104) δ ( v )1 X> VaR α ( X ) (cid:105) . (A.8)Since the VaR only depends on the distribution of the r.v. and not on its specific values, we have D t VaR α ( X ) = 0 and thenCondition (A.2) becomes d (cid:88) j =1 E (cid:20) (cid:90) T D t X j · u t dt (cid:12)(cid:12)(cid:12)(cid:12) F, G (cid:21) = 1 . (A.9)By assumption of Theorem A.3, for j = 1 , . . . , d , X j = g j ( W , . . . , W k ), where ( W , . . . , W k ) is a d -dimensional standardBrownian motion. In this case, it is straightforward to conclude, for 0 ≤ t ≤ D t X j = ∇ g j ( W , . . . , W k ) . (A.10)For each i , choosing v i = ( v i , . . . , v ki ) ∈ D , and setting v t ≡ v i , we find that Conditions (A.2) and (A.4) with T = 1 aresatisfied if ∇ g i ( W , . . . , W k ) · v i = 0 and (cid:88) j (cid:54) = i ∇ g j ( W , . . . , W k ) · v i = 1 . (A.11)Weaker conditions could be pursed, but are outside the scope of this paper. Moreover, by Equation (A.1), we must have δ ( v ) = δ (cid:32) k (cid:88) m =1 v mi e m (cid:33) = k (cid:88) m =1 (cid:18) v mi (cid:90) e m dW t − (cid:90) D t v mi · e m dt (cid:19) (A.12)= v i · W − (cid:90) k (cid:88) m =1 D mt v mi dt = v i · W − (cid:90) Tr( D t v i ) dt, (A.13)where e j is the vector of zeros with 1 in the j -th position and D jt is the Malliavin derivative with respect to the j -thcomponent of the Brownian motion. This implies E [ X i | X = VaR α ( X )] = E (cid:20) X i (cid:18) v i · W − (cid:90) Tr( D t v i ) dt (cid:19) X> VaR α ( X ) (cid:21) E (cid:20)(cid:18) v i · W − (cid:90) Tr( D t v i ) dt (cid:19) X> VaR α ( X ) (cid:21) . (A.14)Moreover, since v i could take the form f i ( W , . . . , W k ), with f i ∈ C ( R k ; R k ), we have D t v i = ∇ f i ( W , . . . , W k ) and thetheorem is proved.4 Yuri F. Saporito, Rodrigo S. Targino Proof (of Theorem 2.2)
In other to deal with the uniform sample case, notice that U j = Φ ( Z j ), where Φ is the cdfof the standard Gaussian distribution. So, X j = h j ( Z ), where h j ( z ) = g j ( Φ ( z ) , . . . , Φ ( z k )) = g j ( Φ ( z )), where Φ ( z ) =( Φ ( z ) , . . . , Φ ( z k )). Notice that ∇ h j ( z ) = ∇ g j ( Φ ( z )) ∗ φ ( z ), where φ = Φ (cid:48) is the pdf of the standard Gaussian distributionand ∗ is the element-wise product. Therefore, it is easy to see that if we take f defined through the equation ˜ f i,m ( z ) = f i,m ( Φ ( z )) /φ ( z m ), where f i = ( f i, , . . . , f i,k ) (similarly for ˜ f ), then Equation (A.6) is satisfied: A h,i ( z ) ˜ f i ( z ) = [0; 1].Moreover, for m = 1 , . . . , k , ∂ m ˜ f i,m ( z ) = ∂ m f i,m ( Φ ( z )) + f i,m ( Φ ( Z )) Z m φ ( Z m ) = ∂ m f i,m ( Φ ( z )) + ˜ f i,m ( z ) z m . Hence π i = − Tr( ∇ f i ( Φ ( Z ))) . By Remark 2.3, the theorem is proved.
A.3 Another Formulation
It is possible to avoid Condition (2.2) if one uses Theorem A.1. Under the Gaussian formulation presented in AppendixA.2, not requiring Equation (A.4) would lead to the formula E [ X i | X = VaR α ( X )] = E (cid:104) ( X i δ ( v i ) − D t X i · v i ) 1 X ≥ VaR α ( X ) (cid:105) E (cid:104) δ ( v i )1 X ≥ VaR α ( X ) (cid:105) , (A.15)where v i = (cid:32) (cid:80) dj =1 ∂ g j ( W , . . . , W k ) , . . . , (cid:80) dj =1 ∂ k g j ( W , . . . , W k ) (cid:33) and we are assuming that none of the terms in the denominator are zero. One should notice that δ ( v i ) could be explicitlycomputed following the same argument as in the proof of Theorem A.3 (see Equation (A.12)). These formulas require thecomputation of the partial derivatives up to second order of the functions g i ’s. Numerically, this could be achieved byautomatic differentiation in cases where analytical calculations would be overly complicated. References
1. Brownlees, C.T., Engle, R.: Volatility, correlation and tails for systemic risk measurement. Available at SSRN (2012)2. Buch, A., Dorfleitner, G.: Coherent risk measures, coherent capital allocations and the gradient allocation principle.Insurance: Mathematics and Economics (1), 235–242 (2008)3. Daveloose, C., Khedher, A., Vanmaele, M.: Representations for conditional expectations and applications to pricing andhedging of financial products in Lvy and jump-diffusion setting. Stochastic Analysis and Applications (2), 281–319(2019)4. Denault, M.: Coherent allocation of risk capital. Journal of risk (1), 1–34 (2001)5. Ewald, C.O.: Local volatility in the heston model: a Malliavin calculus approach. Journal of Applied Mathematics andStochastic Analysis (3), 307–322 (2005)6. Fourni´e, E., Lasry, J.M., Lebuchoux, J., Lions, P.L.: Applications of Malliavin calculus to Monte-Carlo methods infinance. I. Finance and Stochastics , 391–412 (1999)7. Fourni´e, E., Lasry, J.M., Lebuchoux, J., Lions, P.L.: Applications of Malliavin calculus to Monte-Carlo methods infinance. II. Finance and Stochastics (2) (2001)8. Glasserman, P.: Measuring marginal risk contributions in credit portfolios. Journal of Computational Finance (2), 1(2005)9. Hillairet, C., Jiao, Y., Rveillac, A.: Pricing formulae for derivatives in insurance using Malliavin calculus. Probab.Uncertain. Quant. Risk (7) (2018)10. Iglehart, D.L.: Simulating stable stochastic systems, v: Comparison of ratio estimators. Naval Research LogisticsQuarterly (3), 553–565 (1975)11. Kalkbrener, M.: An axiomatic approach to capital allocation. Mathematical Finance (3), 425–437 (2005)12. Koike, T., Hofert, M.: Markov chain monte carlo methods for estimating systemic risk allocations. Risks (1), 6 (2020)13. Koike, T., Minami, M.: Estimation of risk contributions with mcmc. Quantitative Finance (9), 1579–1597 (2019)14. Mainik, G., Schaanning, E.: On dependence consistency of covar and some other systemic risk measures. Statistics &Risk Modeling (1), 49–77 (2014)15. Marshall, A.W., Olkin, I.: Families of multivariate distributions. Journal of the American statistical association (4), 53 (2017)20. Targino, R.S., Peters, G.W., Shevchenko, P.V.: Sequential Monte Carlo samplers for capital allocation under copula-dependent risk models. Insurance: Mathematics and Economics61