Aggregation of Dependent Risks with Heavy-Tail Distributions∗

Straightforward methods to evaluate risks arising from several sources are specially difficult when risk components are dependent and, even more if that dependence is strong in the tails. We give an explicit analytical expression for the probability distribution of the sum of non-negative losses that are tail-dependent. Our model allows dependence in the extremes of the marginal beta distributions. The proposed model is flexible in the choice of the parameters in the marginal distribution. The estimation using the method of moments is possible and the calculation of risk measures is easily done with a Monte Carlo approach. An illustration on data for insurance losses is presented.


Introduction and Background
One of the most difficult problems of risk analysis is modelling dependence of extreme cases. This is the evaluation of multivariate catastrophic scenarios "when everything goes wrong". This question appears in all domains of risk management. For example, when companies analyse their day-to-day losses, which can be interconnected, they envision the possibility of a very large loss that would threaten the firm's economic viability. So, when deciding managerial strategies for risk mitigation and control, methods to evaluate the aggregate risk, i.e. the sum of risk from several dimensions, become an essential tool.
The assumption that losses are generated by a Gaussian statistical distribution is too restrictive when analysing risks, because normality assumes that extreme observations have too a small probability of occurrence. Moreover, it is well-known that for a multivariate Gaussian distribution, and for elliptical statistical distributions in general, large losses in each dimension would tend to occur independently. In practice, this may not be true and this phenomenon is known 1, 2 as "the devil is in the tails", referring to the fact that joint extremes mean that the distribution is heavier in the joint tails. 3 In this article, we propose a statistical multidimensional model that leads to analytical expressions for the sum of non-negative risks. This model can handle general tail-dependence between the marginal distributions that are heavy-tailed. Moreover, our model has also an analytical expression for the multivariate probability distribution and we show that it can be fitted easily. As a consequence of its simplicity, our results have many practical implications. For instance, in the example mentioned above, when direct methods to address the evaluation of risks arising from several sources with dependencies between the different dimensions are needed. Our model allows evaluating worst-case scenarios of a "perfect storm" and so, assessing the risk that several losses could occur simultaneously.
For decades, the difficulties of handling the analytical expressions of general multivariate distributions and of the distribution of sums of dependent random variables have challenged progress in the area of risk analysis and, as a consequence, the assumption of independence and normality has prevailed. Moreover, since no analytical expressions for risk measures such as the value-at-risk (VaR) and the tail value-at-risk (TVaR) are available for most models, estimation for risk analysis has remained difficult. 4,5 As opposed to a modelling approach, empirical analysis in this context refers to calculations based on sampling information, but this has the major disadvantage that large samples are required, especially when addressing confidence levels of 95% or beyond, which means analysing only 5% or even less observations.
Dependence of losses can be introduced via copulas. 6-10 Popular multivariate distributions for losses are: the multivariate Pareto distribution, 11 multivariate mixtures of exponential distributions 12 and the multivariate second kind beta distribution. 13 The widespread use of copula functions can only be understood because it is a simple way to overcome the connection between the marginal distributions and the multivariate distribution, while keeping mathematical expressions tractable. However, a suitable copula should be able to model dependence in the tails in a flexible way. Most copula families either have a fixed structure or cover only a small range of possibilities for the strength of dependence. Tail dependence has been studied asymptotically. 14 Elliptical distributions have been analysed, 15 while some authors 16 concentrated on the tail pattern for multivariate Gaussian distributions. Many contributions have focussed on the tail behaviours of copulas, 17,18 in particular Archemedian copulas. 19,20 Another concept that we address in this article is the sum of risks. The notion of aggregation in probability theory refers to any function that combines a multivariate random vector into a single random variable. By far, the sum of the losses is the most common type of aggregation, but unfortunately, analytical formulas for guillen-sarabia-prietojorda-rev1 Aggregation of Dependent Risks 3 the sum of random variables are difficult, especially when non-elliptical distributions are considered. 21,22 Although the aggregated risk could be easily obtained for independent losses, this assumption is in most cases too restrictive.
In the context of a general multivariate model, if the sum of risks is fully characterized, then it is easy to compute its expectation and additional risk measures. This result, in practice, allows managers to decide whether to retain the risk or to transfer all, or part of it, to an insurer. All over the world, this decision is the core of what chief risk officers do. The practical interest of what we propose in this article is that a simple new method is provided to approximate the aggregated risk. 23 In the following sections, a general model for dependent risks is presented, where the individual risks are modelled according to a second kind beta distribution. The beta distribution has a wide variety of shapes and it is characterized by a Pareto-like heavy-tail behaviour, which means that extremes are observed with a higher probability than in other standard distributions. In our modelling approach, dependence is introduced in the form of a common factor. We discuss the aggregated random variable for which we obtain explicit analytical formulas for the probability density function. An estimation procedure based on maximum likelihood with initial estimates for the parameters based on the method of the moments is presented and, finally, an illustration with real data is discussed.

Univariate Distributions and basic Risk Measures
We assume a heavy-tailed distribution for the univariate losses. Our choice for the univariate model is a second kind beta distribution that is also known as the beta prime or the inverted beta distribution. The choice is motivated by the fact that this is a statistical distribution for positive values that has three parameters. It has a flexible form with one single mode for small values and a long tail. Definition 1. A random variable X is known as a second kind beta distribution when its probability density function (pdf) is of the form, and f (x; p, q, λ) = 0 if x < 0, where p, q, λ > 0 and B(p, q) denotes the beta function. We represent a second kind beta distribution as X ∼ B2(p, q, λ).
A second kind beta distribution random variable corresponds to the Pearson VI distribution 23,24 in the classical Pearson systems of distributions. If we set p = 1 in (1), we obtain a Pareto II distribution Pa(q, λ) and in this case, many risk measures can be obtained in a closed form.
The cumulative distribution function (cdf) of the second kind beta distribution is given by, where IB(x; p, q) represents the incomplete beta function.
Our objective is to calculate basic risk measures in the univariate case. We first define the basic measures and then we show their corresponding expression for a random variable that follows a second kind beta distribution.
Definition 2. Let X be a random variable with cdf F X (x), the value at risk at the α level is defined by, 4, 13, 25 For the second kind beta distribution the V aR can be computed using a simple expression. 22 Proposition 1. Let X ∼ B2(p, q, λ) be a second kind beta distribution, the value at risk at the α level, denoted as V aR[X; α], is given by, where IB −1 (α; p, q) is the inverse of the incomplete beta function and we also include a reference to the parameters of the second kind beta distribution.
The V aR and the T V aR are two of the most popular risk measures.
Definition 3. Let X be a random variable with cdf F X (x), the tail value at risk at the α level is defined by, 13 Proposition 2. If X is distributed according to a second kind beta distribution, X ∼ B2(p, q, λ), the T V aR can be computed as follows, where E[X] = λp q−1 and q > 1 and we also include a reference to the parameters of the second kind beta distribution.
Proofs follow immediately from the definitions given above.

Multivariate extension
In this section, we consider a vector of dependent risks, for which observations in a fixed period of time are available. The corresponding random variables are denoted by {X i , i ∈ N}, where X i is a positive amount. We usually say in insurance that X i refers to the loss corresponding to the ith claim or to the ith risk.
We assume that {X i } is a sequence of dependent random variables with second kind beta distributions.
Usually, the assumption of the independence in the classical risk models leads to simpler formulas than the dependence assumption. However, the formulas for the aggregated risks obtained in the case of dependence between risks are not simple in general.
A stochastic multivariate representation in terms of ratios of independent gamma random variables is our starting point. Lemma 1. Suppose that G p and G q denote independent Gamma random variables, then, by definition, the new random variable resulting from the ratio between the previous two random variables follows a second kind beta distribution: Based on the previous result, the multivariate extension is obtained in the following lemma.
In model (8), each marginal component is distributed according to a second kind beta distribution. Indeed, X i is a second kind beta distribution with shape parameters p i and q 0 , scale parameter λ i and location parameter τ i , i = 1, 2, . . . , d.
In the rest of the paper we assume τ i = 0, i = 1, 2, . . . , d. Moreover, it should be noted that the random variable G 0 is present in all components, and so marginal distributions are dependent. Theorem 1. The pdf of the multivariate variable X = (X 1 , . . . , X d ) where the components are defined in (8), is given by .
Proof. The proof is direct taking into account that the joint probability density function can be written as where f Gi (z) represents the pdf of a gamma random variable G i , i = 0, 1, . . . , d.

Parameter Estimation
Let X 1 , X 2 , . . . , X n be a random sample of d−dimensional vectors from (8). The maximum likelihood estimators of the parameters are given by, where the likelihood function is regular. In order to find the maximum likelihood parameter estimates in (10), a numerical optimization requires initial values. Using the mean, the mode and the variance, the following expressions can be obtained where µ i , σ 2 i , m i , are, respectively, the mean, the variance and the mode of the i-th random variable. Obtaining parameter estimates with these expression is generally known as the method of moments procedure.
With the previous starting values that result from the method of moments, loglikelihood maximization in (10) is done via a standard optimization routine, which is implemented by the authors in R.

Distribution of the Aggregated Risks
Our main result is that, for the particular case where λ i = λ j , i = j, the distribution of the aggregated risks in (8), its VaR and its TVaR can be obtained in a closed form. So, risk analysis calculations become straightforward.
Theorem 2. The pdf of the S n = n i=1 X i , where the components are defined in (8) and λ 1 = . . . = λ n = λ is given by and The proof of this result can be found in Sarabia et al. 12 Lemma 3. If λ 1 = λ 2 = . . . = λ d = λ the distribution of the convolution S d = X 1 + . . . + X d is a second kind beta distribution S d ∼ B2(p 1 + p 2 + . . . + p d , q, λ). Then, the V aR and the T V aR are available in a simple closed form: V aR[S X ; α, p, q 0 , λ] = λ IB −1 (α; p 1 + p 2 + . . . + p d , q 0 ) 1 − IB −1 (α; p 1 + p 2 + . . . + p d , q 0 ) , jorda-rev1 Aggregation of Dependent Risks 7 and T V aR[X; α; p, q 0 , λ] = E[X] , Theorem 3. Let X be a random vector defined by the stochastic representation (8), and consider the distribution of the aggregated risks S d = X 1 + . . . + X d . Then, the probability density function of S d is given by, Proof. The distribution of the sum can be written as S d = U/V , where U = λ 1 G 1 +. . .+λ d G g and V = G 0 . Then, U is the distribution of the sum of independent gamma random variables and its probability density function is, 26 and 0 elsewhere, where C is given in (20) and ρ and δ k were defined previously. Finally, making the change of variable S d = U/V , we obtain the expression for f S d (x; p, q 0 , λ).

Case study
Let us consider losses X i from three sources, i = 1, 2, 3. Without loss of generality, we assume they are all greater than zero. Note that losses can be defined as excesses above a certain threshold level that is pre-specified. Here all losses are above zero. We are interested in the stochastic behaviour of the random vector (X 1 , X 2 , X 3 ) and the distribution of the sum X 1 + X 2 + X 3 . Our data set was obtained from 482 individual customers of an insurance company, that have been insured for a maximum of 10 years. All customers have two types of policies: a motor insurance coverage and a homeowners coverage. In this sample, all customers had had losses in the two insurance policies during the last decade. Losses were classified in three major sources: bodily injuries (BI), property damage (PD), and some home insurance (HO) losses, corresponding to X 1 , X 2 and X 3 , respectively.
We summarize the aggregated overall loss, but we do not take into account that exposure to risk may have been lower than ten years. Moreover, we do not correct for inflation, hence these losses are expressed in current monetary units (thousand Euros). Descriptive statistical measures are presented in Table 1 and histograms are presented in Figures 1 and 2 . A multivariate second kind beta distribution is fitted and maximum likelihood parameter estimates are obtained with the proposed set of initial estimates. There was no single mode because no exact values were repeated in the sample. So to compute the mode for each random variable, samples had to be grouped in intervals and the middle point of the most frequent interval was selected. A sensitivity analysis was carried out and no substantial change was obtained when increasing the number of intervals defined. a To evaluate the risk of the sum of the three types of losses via the V aR and the T V aR, two possible methods are possible. First, empirical analysis of the sample and second, fitting a multivariate model and finding the risk measures of the sum with expressions given in Section 5.
We start by calculating the empirical sample distribution and we obtained an empirical V aR equal to 24.239, 47.769 and 306.126 thousand Euros for confidence levels equal to 95%, 99% and 99.9%, respectively. These empirical results can be interpreted as a rough measure. For instance in the case of a 99.9% confidence, we estimate that one in every thousand clients will have accumulated losses above 306 thousand Euros over a decade. Conversely, the figure corresponding to the 95% jorda-rev1 Aggregation of Dependent Risks 9 confidence means that one customer in every twenty customers has an accumulated loss in the three coverages of motor and home insurance above 24 thousand Euros.
The corresponding T V aR estimates are equal to 63.736, 187.657 and 319.974 thousand Euros for tolerances equal to 95%, 99% and 99.9%, respectively. Our criticism about the empirical risk measures is that they may not be correctly extrapolating the large values due to the lack of data. Given that we only have 482 observations, it is difficult to find risk measures with confidence levels beyond 95%, as they would be based on less than twenty data. In this case, there are plenty of available softwares that extrapolate and interpolate values.
As an alternative, the fitted multivariate beta of the second kind distribution was used to determine the corresponding risk measures. So 10, 000 risk estimates from Monte Carlo simulations with repeated independent samples of size equal to one million observations were generated assuming that losses were represented according to the multivariate pdf. For each simulated sample, risk measures at the 95%, 99% and 99.9% confidence level were obtained. The average and the standard deviation of the 10, 000 risk estimates and the 95% confidence level to evaluate the estimation errors of the risk measures. The fitted risk measures are shown in Table  2  This analysis shows that the empirical T V aR in general overestimates the fit- ted T V aR, which means that the empirical procedure may be producing distorted results due to the extrapolation in the classical empirical analysis when the sample is not too large. The empirical V aR only underestimates the fitted V aR for confidence levels under 95% when the extrapolation of the empirical procedure is not necessary because there are still 25 observations in the tail.

Conclusions and further work
We have studied some classes of multivariate dependent second kind beta distributions with applications in risk aggregation. Individual risks are represented according to second kind beta distribution, which is able to accommodate a variety of shapes and is characterized by a Pareto tail behavior. The second kind beta distribution is considered a suitable model in risk management because it has heavy tails and many of the risk measures are available in a closed form. The dependence arises in the model in the form of a common factor.
Our case study indicates that when sample size is moderate, or confidence levels are high, this aggregation model can be useful to estimate risk measures of dependent heavy-tailed risks, because the empirical classical measures may underestimate or overestimate the tail associations between components of the multivariate random vector.