High level quantile approximations of sums of risks

Abstract The approximation of a high level quantile or of the expectation over a high quantile (Value at Risk (VaR) or Tail Value at Risk (TVaR) in risk management) is crucial for the insurance industry.We propose a new method to estimate high level quantiles of sums of risks. It is based on the estimation of the ratio between the VaR (or TVaR) of the sum and the VaR (or TVaR) of the maximum of the risks. We show that using the distribution of the maximum to approximate the VaR is much better than using the marginal. Our method seems to work well in high dimension (100 and higher) and gives good results when approximating the VaR or TVaR in high levels on strongly dependent risks where at least one of the risks is heavy tailed.


Introduction
Because of regulatory rules (such as Solvency 2 in Europe) or for internal risk management purposes, the estimation of high level quantiles of a sum of risks is of major interest both in nance and insurance industry. Consider an insurance company that has a portfolio of d ≥ (possibly) dependent risks which is represented as a random vector X = (X , . . . , X d ) with cumulative distribution function (c.d.f.) F(x , . . . , x d ). We assume that all the risks are almost surely positive but we do not assume that they are identically distributed. Let S denote the aggregated risk S = X + · · · + X d .
We are interested here in the computation of the Value-at-Risk (VaR) and the Tail Value-at-Risk (TVaR) of the sum, VaR p(S) = F S (p) and TVaR p(S) = − p p VaR u(S) du, for con dence levels p ∈ ( , ) near , where F S is the c.d.f. of S and F is its generalized inverse. Problems like this arise for insurance companies, for example, which are required to maintain a minimum capital requirement which is typically calculated as the VaR for the distribution of the sum at some high level of probability. Even when the distribution function F is known, good estimations for VaR p(S) are not trivial since they require a precise calculation of F S , which is given by the following integral This integral is more di cult to approximate when d is large and it is usually more e cient to apply Monte Carlo methods to estimate it (for a comprehensive introduction to Monte Carlo methods see [25]). Nevertheless, when p is near , the number of replications required to give precise estimations is also large, so new methods are always well received. Classical Extreme Value Theory (EVT) allows one to get some estimation of the VaR ( [15,26]), but EVT based methods requires an estimation of the EVT parameters, which is known to be not an easy task. Recently, in [8,9,17], some approximations on the VaR are obtained for some speci c models; see also [18] where theoretical results on the asymptotic behavior of the ratio VaR p (Xi) are given. Results for the tail distribution of the sum of dependent subexponential risks are obtained in [19] and also in [20] when risks are non-identically distributed and not necessarily positive. In [5], an algorithm to compute the distribution function of S is proposed and in [12], bounds are obtained. Nevertheless, these results may be used to estimate VaR p(S) for small dimensions (d < ) and give ranges in dimension or .
We present a method which seems to be quite accurate even for a large number of summands, in the order of several hundreds for instance (see Sections 6.2 and 6.3 for simulations in dimension and dimension ). Our method will be compared to the EVT driven ones as well as to the Monte Carlo method, especially for very high level quantiles and in dimension greater than . are not necessary. This also means that VaR p(M) can be easily calculated for any given level of con dence p, at most a simple numerical inversion is needed.
In this paper we give some conditions on X under which the Value-at-Risk and the Tail Value-at-Risk of the sum and maximum are asymptotically equivalent in the sense that there exists some ∆ ≥ such that and where a(t) ∼ b(t) when t → l, for l ∈ [−∞, ∞] means throughout this paper that lim t→l a(t) b(t) = . This result is interesting because it allows to estimate the VaR (or TVaR ) of the sum by using the VaR (or TVaR ) of the maximum, which is easier to calculate, and the estimation of ∆.
For random vectors with common marginals (Fréchet, Gumbel, Weibull) and an Archimedean copula dependence structure [3] and [2] get an asymptotic approximation of the tail of S. These results are generalized in [4] to other dependence structures. In [6], the same results are obtained in the multivariate regularly varying framework. Examples in which the limiting constant ∆ can be computed explicitly are also given in [16]. Finally, we would like to mention [22] which is related to our work, in an independent framework and for Pareto marginals.
In this paper, we consider the more general framework with non common marginals (and regularly varying tails). We emphasize that our method applies when there are dependences between risks as well as the presence of heavy tailed marginal distributions (see Section 4 for more details). This may be a typical context for risk management applications in insurance and nance. Moreover, the proposed method is tractable, even in high dimension (dimension tested).
The paper is organized as follows. In Section 2, we recall the de nition of regularly varying function and then present conditions under which the VaR and TVaR of the sum and the maximum are asymptotically equivalent. In Section 3, we give classes of random vectors satisfying our hypothesis. Section 4 is devoted to a methodology for the estimation of ∆. In Section 5, we give explicit expressions of the VaR on some speci c models (introduced in [23,27] and also considered in [13] where the expression of the VaR is derived). In Section 6, we compare our method with classical ones on several models. Conclusions are given in Section 7.

Asymptotic results on the VaR and the TVaR of the sum and the maximum
In this section, we will rst recall the de nition of regularly varying functions.
De nition 1. Let f be a positive measurable function on R+. We say that f is regularly varying at in nity of for any x > . Similarly, we say that f : R+ → R+ is regularly varying at if we replace t → ∞ by t → .
Regularly variation of f at a > is de ned as regularly variation at in nity for the function f (a − /t) (see [15], page 565). If ρ = then f is said to be slowly varying.
Examples of regularly varying distributions are Pareto, Cauchy, Burr and stable with exponent α < .

De nition 2.
A random variable X with distribution function F is said to have a regularly varying upper tail if its survival function F is regularly varying at in nity.
Let δ be the real valued function de ned by δ(t) = F S (t)/F M (t). Throughout this paper we will consider the condition The following result is somewhat a folklore theorem, it links the Value-at-Risk of the sum and the maximum in case where the survival function of the maximum, F M , is regularly varying. The result still holds for the TVaR . Recall that we do not assume that the marginal distributions are either identically distributed or independent. Proposition 2.1. Let X = (X , . . . , X d ) be a vector of positive random variables (r.v.s). Suppose that assumption (2.1) holds and that F M is regularly varying with index −ρ. Then, Classes of random vectors that satisfy the assumptions of Proposition 2.1 will be given in Section 3 while in Section 4 we will provide a method to estimate ∆.

Random vectors where the limit ∆ exists
In this section we explore several situations in which the limit ∆ exists and Proposition 2.1 applies.
. Multivariate regular framework Alink et al. ([3], [2] and [4]) studied the asymptotic behavior of the tail of the sum when the marginals of the vector X = (X , . . . , X d ) are identically distributed as one of the three extreme value families: Gumbel, Fréchet or Weibull and when the dependence within the vector is given by an Archimedean copula. Then Barbe et al.
( [6]) generalized these results under the framework of the multivariate regular variation distributions. Their main contribution is the explicit calculation of the limit where F is the common distribution function of the marginal risks X , . . . , X d .
This kind of results suggest that we may approximate the VaR (and TVaR) of the sum simply by the VaR (and TVaR) of X . This point will be detailed in Section 6.4 where it will be shown empirically that maximum based estimation gives indeed better results than F based one.
Let us recall the de nition of multivariate regularly varying random vectors.
Using the L norm, X = |X | + · · · + |X d |, the L ∞ norm, where |µ| is the total mass of the measure µ. So that, when X is multivariate regularly varying Proposition 2.1 applies.
We are also interested in random vectors whose coordinates are not identically distributed. Results for identically distributed marginals will not lead to results for arbitrary marginals. This is the purpose of the next section where di erent kinds of dependence structure are also considered.

. Examples where condition (2.1) holds.
In this section we show that condition (2.1) holds for three classes of multivariate distributions, namely those for which -a regularly varying marginal clearly dominates the other marginals, -the dependence structure is the survival of a regularly varying Archimedean and the marginals are regularly varying, -the dependence structure is regularly varying Archimedean and the marginals are regularly varying.
We now state our result when one marginal is regularly varying and dominates the others.
Proof. The proof splits into two parts, one showing that F M ∼ F at in nity, one showing that F S ∼ F at in nity. The result follows in combining these two asymptotic equivalences.
(i) If X exceeds t, so is M, and if M exceeds t, at least one of the X i does. Therefore, we have This ensures that under the assumptions of Proposition 3.1, F M ∼ F at in nity.
Furthermore, for any positive ϵ, decomposing the event {S > t} according to whether max ≤i≤d X i ≤ tϵ or not, we have Applying Bonferroni's inequality, we obtain Since F is regularly varying of index ρ say, and dominates the other F i , we obtain that Since ϵ is arbitrary, we have lim sup We now consider a dependence structure between the components of the random vector given by an Archimedean copula or the survival copula of an Archimedean copula, and give a su cient condition for condition (2.1) to hold. We rst recall the de nitions of Archimedean copulas and survival copulas.

De nition 4. (Archimedean Copulas) A generator is a function ψ from [ , ] to [ , ∞] such that
(i) ψ is decreasing with ψ( ) = , (ii) the rst d derivatives of ψ ← exists, (iii) for any k = , , . . . , d and any t positive, where ψ ← denotes the pseudo-inverse of ψ de ned by The Archimedean copula C with generator ψ is the distribution function on [ , ] d de ned by De nition 5. (Survival copula) Given a copula C, we de ne: with (U , . . . , U d ) having C as distribution function. C * is a copula known as the survival copula of C.
We can now extend the result of Alink et al. [2,3,4] and Barbe et al. [6], to a situation where the marginal distributions are not identical.
Let the dependence structure of X be given by one of the following: (i) a survival copula of an Archimedean copula with generator ψ which is regularly varying at with negative index, (ii) an Archimedean copula with generator ψ which is regularly varying at with negative index.
Then X is multivariate regularly varying, and condition (2.1) holds.
Note that since a regularly varying function is ultimately positive, the a i are nonnegative. If a i is positive, then (3.3) implies that F i is regularly varying with the same index of regular variation as h. In particular, if one marginal tail is regularly varying and dominates the others, then h could be the corresponding survival function.
Proof. (i) First we assume that the dependence of X is the survival copula of an Archimedean copula with generator ψ which is regularly varying at with negative index. Using Bonferroni's identity and agreeing that a sum over an empty set is 0, we have, for any x , . . . , x d positive, If I is such that a i = for some i in I, we de ne We then have as t tends to in nity.
If I is such that all a i , i ∈ I, are positive, we have Since h is regularly varying with index −ρ say, we have Since ψ is regularly varying at with index −θ say, the uniform convergence theorem (Theorem 1.2.1 in [10]) and (3.7) ensure that as t tends to in nity. Since ψ ← is regularly varying with index − /θ, we then have, using (3.6) and the uniform convergence theorem, as t tends to in nity. Note that if we take the limit of the right hand side of (3.8) as one of the a i tends to , we obtain . Therefore, as long as we agree that / = ∞ and /∞ = , we may capture (3.5) in (3.8). Then, considering (3.4), and using that at least one of the a i does not vanish, we obtain as t tends to in nity. It then follows from Theorem 6.1 in [24] that the distribution of X is multivariate regularly varying. Condition (2.1) then follows from (3.1).
(ii) We now assume that the dependence of X is an Archimedean copula with generator ψ which is regularly varying at with negative index. By de nition, for any x , . . . , x d positive, Since ψ is regularly varying at with index −θ say, and h is regularly varying with index −ρ say, then − ψ ← is regularly varying at with index θ − and ψ • h is regularly varying at in nity with index −ρθ. Then, by using the same arguments as above we can conclude that as t tends to in nity and X is multivariate regularly varying.
Notice that Proposition 3.2 implies that a random vector X with regularly varying marginals, for example Pareto distributed marginals not necessarily with same scale or shape parameters, and dependence structure given by one of the copulas listed below, satis es the assumption of Proposition 2.1. Possible dependence structures are: -independence (recall that the independent copula is an Archimedean copula with generator ψ(t) = − ln(t), and thus regularly varying at 1 with index − ), -Gumbel copula with parameter θ ≥ (which is an Archimedean copula with generator ψ(t) = − ln(t) θ and thus regularly varying at with index −θ), -survival copula of a Clayton copula with parameter θ > (which is an Archimedean copula with generator ψ(t) = (t −θ − )/θ and thus regularly varying at with index −θ).

Approximation of the limit ∆
In this section, we assume that the limit ∆ exists and we show how to estimate it using samples of X. We will use this estimation to approximate VaR −p (S), for di erent values of p close to using Proposition 2.1.
Recall that δ is the real valued function de ned by δ(t) = F S (t)/F M (t) and continue to denote by ∆ its limit at in nity if it exists.
If a sample of X is available, the function δ can be estimated using the empirical cumulative distribution function (e.c.d.f.) of S and M. As we assume that F M can be easily calculated by the c.d.f. F of the portfolio, at least two versions of the empirical delta should interest us: where F S and F M are the e.c.d.f.s of S and M respectively, based on the sample of X. The rst version δ may be more tractable statistically, while the secondδ has the nice property thatδ ≥ . In order to obtain some insight on the convergence of δ to its limit ∆, we plot, in Figure 1, functions δ andδ for four di erent models which are multivariate regularly varying. Pareto distributions with di erent tail indexes and di erent dependence structures. From top-left to bottom-right, we nd: 1) independent Pareto distributions with tail index one; 2) the tail index is still one but dependence is given by a Gumbel copula of parameter . ; 3) independent Pareto distributions: ve with tail index one and the other ve with tail index ; 4) the same as 3) but dependence is given by a Gumbel copula of parameter . . In the rst model (sum of 10 independent Pareto distributions with tail index 1) we notice that the limit δ(t) seems to be but the convergence is not fast enough to consider using this limit to approximate VaR p(S) even for higher con dence levels p. For the second model (sum of 10 Pareto distributions with tail index 1 and dependence structure given by a Gumbel copula of parameter 1.5) the convergence is a lot faster, δ(t) seems to be close to its limit for t greater than the VaR at the % con dence level. The two models in the lower side (sum of 10 Pareto distributions: ve with tail index 1 and ve with tail index 3 both in the independent and Gumbel copula dependent case) behave the same as the ones in the upper side. The models on the right side of Figure 1 correspond to cases where our method will be applicable: the limit ∆ is reached by δ(t) for t near the VaR . . These models exhibit a strong dependence combined with at least one of the marginal risks with a very heavy tail. Even if this is a limitation of our method we should remark that this kind of models are also those where Monte Carlo methods are less e cient to approximate the VaR or the TVaR , so that it may be interesting to have an alternative method of approximation. converges in distribution to a Gaussian process with zero mean and covariance given by for t ≤ t . Thus, given any sequence < t < · · · < t k , the vector converges in law to a centred Gaussian vector with covariances given by for any i ≤ j. As a consequence converges to a normal distribution with zero mean and variance If we assume that the values t i are large enough, the approximation δ(t i ) ≈ ∆ holds for each i = , . . . , k and the variance (4.1) can be approximated by In practice we should plot points S (i) , δ(S (i) ) where S ( ) < · · · < S (n) is the ordered sample of S and then choose a threshold in such a way that the approximation δ(S (n−i) ) ≈ ∆ holds for any ≤ i ≤ k. The choice of the threshold is a recurrent and di cult problem in EVT, for which few theoretical results exist and are generally hardly applicable in practice. We propose then to estimate ∆ by ∆ = k k i= δ(S (n−i) ). (4.2) As an example, the behavior of δ(x) for the Pareto-Clayton model, which will be described in Section 5, may be seen on Figure 2. The estimation ∆ is represented by the solid line while dashed lines are for the estimated % con dence interval. See also Figure 3 for the shape of the δ function and the limit ∆.

Some Explicit Calculations
In this section we will consider a simpli ed model in order to obtain explicit formulas for F S and F M and to better understand the scope and the limitations of our ∆ estimation. The model is described by the following compound process: let Λ be a positive random variable and let X = (X , . . . , X d ) be a random vector such that for each x , . . . , x d ≥ .
That means that conditionally on the value of Λ the marginals of X are independent and exponentially distributed. In general, the nal distribution of X does not have independent marginals and they are not exponential either. Actually the dependence structure of X and its marginal distributions will depend on the distribution of Λ.
Some particular Λ distributions de ne some well-known models in which the explicit calculation of F S and F M is possible. For example when Λ is Gamma distributed, then the marginals of X are of Pareto type with dependence given by a survival Clayton copula. When Λ is Levy distributed the marginals will be Weibull distributed with a Gumbel survival copula. These models have been studied in [23,27] and used in [1] where explicit formulas for ruin probabilities have been derived. In [11,21], explicit results for the minimum of some risk indicators are obtained for this kind of models. We also would like to mention that the computation of the VaR for this model is given in [13].
Let us consider the case where Λ is Gamma(α, β) distributed with density In this case, the X i 's are Pareto(α, β) distributed with tail given by and the dependence structure is described by a survival Clayton copula with parameter /α. Through this paper we will refer to this model as a Pareto-Clayton vector with parameters (α, β). This model is a particular Multivariate Pareto of type II with location parameters µ i = and scale parameters σ i = β for i = , . . . , d (see [27]). As already noticed in Section 3.2, this model satis es the hypothesis of Proposition 3.2 so that the limit ∆ exists.
In the Pareto-Clayton model, the exact distribution function of S = d i= X i can be calculated. Conditionally on Λ = λ, the sum S is Gamma distributed with parameters ( /λ, d), distribution also known as the Erlang distribution. Then, as here we are assuming that Λ is Gamma(α, β) distributed, the total distribution of S is the result of compounding two Gamma distributions, more precisely It is well known that the result of this compound distribution is the so-called Beta prime distribution (see [14]). The c.d.f. of S can be expressed in terms of F β , the c.d.f. of the Beta(dβ, α) distribution, as Naturally, the inverse of F S can also be expressed in function of the inverse of the Beta distribution In this example, the δ function is explicitly calculated (see Figure 3). Moreover, computer algebra softwares allow us to calculate explicitly the limit ∆ for speci ed parameters. Vertical lines are displayed at the VaR of the sum at con dence levels %, %, . %, . %. The limit ∆ ≈ .
is represented by the horizontal line.
In order to see how fast the function δ converges to its limit ∆, we plot the function p → δ(VaR p(S)) for di erent values of the parameter α and di erent dimensions d (see Figure 4). We remark that δ(x) is already very close to ∆ when x = VaR . (S), for α ≤ . . The lower the value α, the atter the tail of δ and thus the limit ∆ is attained rapidly. Remark that the lower the levels of α, the heavier the tails of the Pareto marginals. Finally, this plot con rms the intuition that for heavier marginals the tail of the sum is better approximated by the tail of the maximum. A similar phenomenon in the i.i.d. case has been noted in [7] when approximating F S by F .

Some numerical examples
In this section we show how the ideas presented in the above section can help to estimate in practice the VaR and the TVaR of a sum at con dence levels close to 1. We compare the estimation done via the ∆-limit estimation (New in the tables below) as described in Section 4 with other common quantile estimation methods, with the same sample size: 1. The direct Monte Carlo quantile estimation (MC). 2. The quantile estimation from a GPD tted distribution where parameters are estimated using maximum likelihood method (GPD 1). 3. The quantile estimation from a GPD tted distribution where parameters are estimated using the moment method (GPD 2). 4. The high quantile estimate based on a method by Weissman [26] (Weiss.).
We rst consider the Pareto-Clayton model presented in Section 5 (dimension and ), where exact values for the Value-at-Risk are computable (see Section 6.1 and Section 6.2). Then, we test our method with a di erent model where exact values are not known. In order to study the performance of our estimator and to compare it with the main competitors, we consider the root-mean-squared error (RMSE) loss function. When n estimations have been performed, it is de ned by where VaR p(S i ) represents the estimate of VaR p(S) for any of the di erent methods presented above, on the ith sample. In the case where the exact value is not known, in Section 6.3, we compare our results on a sample of size with several methods (1-4 above) to a Monte Carlo quantile estimation based on a very large sample of size × . This last estimation is considered as the exact VaR value in the RMSE computation.

. Pareto-Clayton model dimension 2
Here we consider the model presented in Section 5. We rst consider d = and α = which corresponds to a model with Pareto marginals with α = and dependence given by a survival Clayton copula with parameter θ = .
In Table 1, the exact VaR at di erent con dence levels (from % to . %) is presented. In Table 2 and Table 3, we present the RMSE criterion in percentage of the real value based on simulations at di erent con dence levels. At each simulation a sample of size in Table 2 and size in Table 3 is used to estimate the VaR. On each method (New, GPD 1, GPD 2 and Weiss) the threshold used on each estimation corresponds to the empirical 95% quantile. Clearly, in term of RMSE, our method performs better than classical methods at each con dence level, even for very high levels. When increasing the size of the sample ( instead of ) classical methods improve but our method still produces the best results.

. Pareto-Clayton model dimension 10
We consider again the Pareto-Clayton model but here d = and α = which corresponds to a model with Pareto marginals with α = and dependence given by a survival Clayton copula with parameter θ = . Results are presented in Tables 4, 5 and 6. As above, on each method (New, GPD 1, GPD 2 and Weiss) the threshold used on each estimation corresponds to the empirical 95% quantile. We mention that even in dimension , the estimation remains e cient for high level quantiles.   We also remark that our method is more e cient than classical ones from level . .

. A model with di erent Pareto marginals and Gumbel copula
We apply now our method to a model where the exact value of VaR p(S) is not known. The model is composed of 150 marginals Pareto(α i , β i ) distributed with parameters α i = − i mod ( ) / and β i = −i mod ( ) for i = , . . . , , where i mod (j) denotes the reminder of i divided by j. The model is then composed of fty Pareto marginals of tail index 0.5, fty of tail index 1 and fty with tail index 1.5, and di erent scale parameters within , , . . . , . The dependence structure is given by a Gumbel copula of parameter . . Recall that for this model, by the comments that follow Proposition 3.2, the limit ∆ exists. Table 7 presents the VaR estimation based on a classical Monte Carlo quantile estimation with a sample of size × . We assume this estimation is the "real VaR" in the computation of the RMSE presented in Table 8. On each method (New, GPD 1, GPD 2 and Weiss) the threshold used on each estimation corresponds to the empirical % quantile. It is notable that our method is very stable with respect to others and is more e cient to approximate the VaR p from p = . .  . Comparison of the method using max(X) vs X The method of estimation of the Value-at-Risk of the sum proposed in this paper relies on the convergence of the function δ(t) = F S (t)/F M (t). When the convergence is assured and it is fast enough, it has been shown that the proposed method gives accurate and stable estimations of the VaR at high levels. In theory, similar results could be obtained if the maximum M is replaced by X where X is assumed to have the heaviest tail in the vector X. In this section we compare numerically the estimation of the VaR using, on one side, δ(t) = F S (t)/F M (t) and, on the other side, δ ′ (t) = F S (t)/F X (t), i.e we compare the approximation of VaR −p (S) by VaR −p/∆ (M) and VaR −p/∆ ′ (X ) where ∆ and ∆ ′ are the approximated limits of δ(t) and δ ′ (t) respectively estimated using (4.2).
We rst consider the model (X , . . . , X ) where X is Pareto distributed with α = . and X , . . . , X are Pareto distributed with α = . The dependence structure is given by a Gumbel copula with parameter 2. Empirical δ and δ ′ functions are displayed in Figure 5.
The δ function becomes almost horizontal before the VaR of the sum at the % con dence level whereas δ ′ does not seem to be close to the limit on the displayed range. Then, the estimation of the VaR using δ ′ seems to be not accurate. This is con rmed by Table 9 where some VaR estimations are presented. From now on, the threshold used for the ∆ and the ∆ ′ approximations using formula (4.2) corresponds to the 95% empirical quantile and for each estimation a sample of size is generated. Even in the case where all the marginal risks are equal the use of the max seems to give better results. We consider the model (X , . . . , X ) where all the X i 's are Pareto distributed with the same index α = .  The dependence structure is given by a Gumbel copula with parameter 2. Empirical δ and δ ′ functions are displayed in Figure 6. As above the δ function seems to converge faster than δ ′ but in this case the di erence is not as important as in Figure 5. In Table 10 some VaR estimations are presented. Again, estimations provided by using the estimation of ∆ are of better quality than the ones provided by using the estimation of ∆ ′ .
Mathematically speaking, some work remains to be done to understand why the approximation of F S by F M is so much better than that by F . This will be the object of further investigations.

Conclusion
In this paper, we give some conditions under which the tail distribution of the sum can be approximated by using the tail of the maximum of a vector. We show how the VaR or the TVaR on high levels for the sum can be approximated, by rst estimating a limiting constant ∆. The models in which our results can be applied include those where marginals are regularly varying and such that dependence is given by an Archimedean copula or survival copula. We do not require the marginals to be identically distributed and the method works for very high dimensions d (d = for exemple). Our method gives a good approximation for the VaR and the TVaR when the convergence of δ(x) to ∆ is fast enough. This generally happens when at least one of the  marginal risks is strongly heavy tailed and when the dependence is strong. In particular, the method is not suitable e.g. for the case of two independent Pareto distributions. We also remark that the models for which our method applies correspond generally to those where Monte Carlo approximations are less e cient and there so is a real need for alternative methods.