11.1.2.1 Main notation

will denote the integer part of any nonnegative real such that .

Let be the probability space on which we will be working.

Let and denote, respectively, the cumulative distribution function (cdf) and the probability density function (pdf) of the standard normal distribution and and the cdf and pdf of the normal distribution with mean and variance .

Let be a random variable (r.v.) Pareto (type I) distributed with shape parameter , pdf denoted by and cdf defined by

11.1

and probability density function (pdf) denoted by .

Note that the inverse function of is given by

11.2

Recall that for , and for , .

We will consider i.i.d. Pareto r.v.'s in this study and denote by the Pareto sum being an -sample with parent r.v. and associated order statistics .

When dealing with financial assets (market risk data), we define the returns as

being the daily price and representing the aggregation factor.

Note that we can also write

In what follows, we will denote by .

Further comments or questions

• Is it still worth considering i.i.d. r.v.'s, whereas most recent research focus on dependent ones?

  Concerning the i.i.d. condition, note that this study fills up a gap in the literature on the sum of i.i.d. moderate heavy r.v.'s (see, e.g., (Feller, 1966); (Hahn et al., 1991), and (Petrov, 1995)). Moreover, in our practical example of log returns (the motivation of this work), the independence condition is satisfied (see, e.g., (Taylor, 1986); (Dacorogna et al., 2001)) and hence is not a restriction in this case of time aggregation.

  Another theoretical reason comes from the EVT; indeed we know that the tail index of the aggregated distribution corresponds to the one of the marginal with the heaviest tail and hence does not depend on considering the issue of dependence.

  Finally, there was still mathematically a missing “brick” when studying the behavior of the sum of i.i.d. r.v.'s with a moderately heavy tail, for which the CLT applies (for the center of distribution!) but with a slow convergence for the mean behavior and certainly does not provide satisfactory approximation for the tail. With this work, we aim at filling up the gap by looking at an appropriate limit distribution.

• Why considering Pareto distribution?

  It is justified by the EVT (see, e.g., (Leadbetter et al., 1983); (Resnick, 1987), and (Embrechts et al., 1997)). Indeed recall the Pickands theorem (see (Pickands, 1975) for the seminal work) proving that for sufficiently high threshold , the generalized Pareto distribution (GPD) (with tail index and scale parameter ) is a very good approximation to the excess cdf of a r.v. defined by :

  

  if and only if the distribution of is in the domain of attraction of one of the three limit laws. When considering risks under the presence of heavy tail, it implies that the extreme risks follow a GPD with a positive tail index (called also extreme value index) , which corresponds to say that the risks belong to the Fréchet maximum domain of attraction (see, e.g., (Galambos, 1978); (Leadbetter et al., 1983); (Resnick, 1987), or (Embrechts et al., 1997)). In particular, for ,

  11.3

  for some constant . It is then natural and quite general to consider a Pareto distribution (with shape parameter ) for heavy-tailed risks.

  A natural extension would then be considering r.v.'s with other distributions belonging to the Fréchet maximum domain of attraction.

• A last remark concerns the parameter that we consider as given in our study. A prerequisite, when working on real data, would be to estimate . Recall that there are various ways to test the presence of a heavy tail and to estimate the tail index, for example, the Hill estimator (see (Hill, 1975)) or the -estimator (see (Kratz and Resnick, 1996)) (see also (Huston McCulloch, 1997); (Beirlant et al., 2004); (Resnick, 2006), and references therein). We will not provide an inventory of the methods, except a brief recall in the next section of an important empirical EVT method used for estimating the heaviness of the tail. Let us also mention a test, easy to use in practice, for the existence of fat tails, namely, the scaling law (see (Dacorogna et al., 2001), Section 5.5). It consists of comparing the two plots, for and 2, respectively, of , ; if the scaling exponent for is larger than for , then it is a sign of the existence of a fat tail. For financial data, there are numerous empirical studies that show the existence of fat tail and that the shape parameter is between 2 and 4 for developed markets (see, e.g., (Jansen and De Vries, 1991); (Longin, 1996); (Müller et al., 1998); (Dacorogna2 et al., 2001); (Dacorogna et al., 2001), and references therein).

11.3.1.1 How to fit for the best mean behavior of aggregated heavy-tailed distributed risks?

Let us start by studying the behavior of the trimmed sum when writing down the sum of the i.i.d. -Pareto r.v.'s (with ), , as

11.10

Much literature, since the 1980s, has been concerned with the behavior of trimmed sums by removing extremes from the sample; see, for example, Hall (1984), Mori (1984), and Hahn et al. (1991).

The main issue is the choice of the threshold , in order to use the CLT but also to improve its fit since we want to approximate the behavior of by a normal one.

We know that a necessary and sufficient condition for the CLT to apply on is to require the summands , , to be -r.v.'s. But we also know that requiring only the finitude of the second moment may lead to a poor normal approximation, if higher moments do not exist, as occurs, for instance, with financial market data. In particular, including the finitude of the third moment provides a better rate of convergence to the normal distribution in the CLT (Berry–Esséen inequality). Another information that might be quite useful to improve the approximation of the distribution of with its limit distribution is the Fisher index, defined by the ratio , which is a kurtosis index. The skewness of and measures the closeness of the cdf to . Hence we will choose based on the condition of existence of the fourth moment of the summands of (i.e., the first order statistics).

The following Edgeworth expansion involving the Hermite polynomials points out that requiring the finitude of the fourth moments appears as what we call the “optimal” solution (of course, the higher order moments exist, the finer the normal approximation becomes, but it would imply too strong conditions and difficult to handle). If denotes the cdf of the standardized defined by , then

11.11

uniformly in , with

The rate of convergence appears clearly as whenever , .

Note that in our Pareto case, the skewness and the excess kurtosis are, respectively,

Therefore we set (but prefer to keep the notation so that it remains general) to obtain what we call an “optimal” approximation. Then we select the threshold such that

which when applied to our case of -Pareto i.i.d. r.v.'s (using (11.17)) gives

11.12

This condition allows then to determine a fixed number as a function of the shape parameter of the underlying heavy-tailed distribution of the 's but not of the size of the sample. We can take it as small as possible in order to fit for the best both the mean and tail behaviors of . Note that we look for the smallest possible to be able to compute explicitly the distribution of the last upper order statistics appearing as the summands of the second sum . For this reason, based on condition (11.12), we choose

11.13

Let us summarize in Table 11.1 the necessary and sufficient condition on (for ) to have the existence of the th moments for the upper order statistics for , respectively, using (11.12) written as .

	0	1	2	3	4	5	6	7

We deduce the value of the threshold satisfying (11.13) for which the fourth moment is finite according to the set of definition of : We notice from Table 11.2 that we would use Zaliapin et al.'s decomposition only when . When considering, as they do, , we would rather introduce the decomposition , with varying from 2 to 5 depending on the value of , to improve the approximation of the distribution of , if we omit the discussion on their conditions.

with							]2,4]
=	7	6	5	4	3	2	1

11.3.1.2 Some properties of order statistics for Pareto random variables

First we apply known results on distribution of order statistics (see, e.g., (David and Nadaraja, 2003)) when considering Pareto distributions. Next we compute conditional distributions of order statistics, as well as conditional moments, to apply them to the Pareto case.

• Distribution of Pareto order statistics

  For -Pareto r.v.'s, the pdf of () and the pdf of the order statistics , (), with , are expressed, respectively, as

  11.14

  and, for ,

  11.15

  When considering successive order statistics, for , for , , with , we obtain

  11.16

  Moments of -Pareto order statistics satisfy (see also, e.g., (Zaliapin et al., 2005); Theorem 1)

  11.17

  and, for ,

  

• Conditional distribution of order statistics. Application to Pareto r.v.'s

  Now straightforward computations lead to new properties that will be needed to build Normex. We express them in the general case (and labeled), with the notation and , and then for Pareto r.v.'s.

  We deduce from (11.14) and (11.15) that the pdf of given , for , is, for ,

  11.18

  and that the joint pdf of given , for , is, for ,

  11.19

  Using (11.14) and (11.16) provides, for ,

  11.20

  Then we can compute the first conditional moments. We obtain, using (11.18) and the change of variables ,

  11.21
  11.22

  For , via (11.19) and the change of variables and , it comes

  11.23

  Moreover, the joint conditional distribution of given , for , denoted by , or when no ambiguity exists, is, for ,

  11.24

  from which we get back the well-known result that are independent of and when and are given and that the order statistics form a Markov chain.

11.3.2.1 A conditional decomposition

Whatever the size of the sample is, because of the small magnitude of , we are able to compute explicitly the distribution of the last upper order statistics appearing as the summands of the second sum defined in (11.10). The choice of allows also to obtain a good normal approximation for the distribution of the trimmed sum . Nevertheless, since and are not independent, we decompose the Pareto sum in a slightly different way than in (11.10) (but keeping the same notation), namely,

11.25

and use the property of conditional independence (recalled in Section 11.1.2) between the two subsums and conditional on (for ).

Then we obtain the following approximation of the distribution of , for (i.e., when the th moment of the largest order statistics does not exist).

Comments

1. The distribution can also be expressed as

   

   and, for

   

   where the convolution product can be numerically evaluated using either the recursive convolution equation , for , (it will be fast, being small) and or, if , the explicit expression (12) (replacing by ) given in Ramsay (2006).

2. Note that we considered i.i.d. Pareto r.v.'s only as an example to illustrate our method intended to be extended to unknown distributions using the CLT for the mean behavior and heavy-tailed distributions of the Pareto type for the tail. Since the exact distribution of the Pareto sum of i.i.d. Pareto r.v.'s is known, we will be able to judge about the quality of the approximation proposed in Theorem 11.2 when comparing with the exact distribution of . We will then compare the respective associated risk measures.
3. Finally recall the following result by Feller (see 1966) on the convolution closure of distributions with regularly varying tails, which applies in our Pareto example but may also be useful when extending the method to distributions belonging to the Fréchet maximum domain of attraction.

Note that this lemma implies the result given in (11.7), and as a consequence in the Pareto case, we have

11.3.2.2 On the quality of the approximation of the distribution of the Pareto sum

To estimate the quality of the approximation of the distribution of the Pareto sum , we compare analytically the exact distribution of with the distribution defined in Theorem 11.2. It could also be done numerically, as, for instance, in Furrer (2012) with the distance between two distributions and defined by , with . We will proceed numerically only when considering the tail of the distributions and estimating the distance in the tails through the VaR measure (see Section 11.4.3). When looking at the entire distributions, we will focus on the analytical comparison mainly for the case (with some hints for the case ). Note that it is not possible to compare directly the expressions of the VaR corresponding to, respectively, the exact and approximative distributions, since they can only be expressed as the inverse function of a cdf. Nevertheless, we can compare the tails of these two distributions to calibrate the accuracy of the approximative VaR since

Moreover, we will compare analytically our result with a normal approximation made on the entire sum (and not the trimmed one) since, for , the CLT applies and, as already noticed, is often used in practice.

Since Normex uses the exact distribution of the last upper order statistics, comparing the true distribution of with its approximation simply comes back to the comparison of the true distribution of i.i.d. r.v.'s with the normal distribution (when applying the CLT). Note that, when extending Normex to any distribution, an error term should be added to this latter evaluation; it comes from the approximation of the extreme distribution by a Pareto one.

Suppose . Applying the CLT gives the normal approximation , with and , where in the case of a Pareto sum, , and . We know that applying the CLT directly to leads to nonsatisfactory results even for the mean behavior, since, for any , the quantity , involving the third moment of and appearing in the error (11.11) made when approximating the exact distribution of by a normal one, is infinite for any . The rate of convergence in is reduced to . When , even if the rate of convergence improves because , we still have (because the fourth moment of does not exist), which means that we cannot get a rate of order .

Now let us look at the rate of convergence when approximating with .

Considering the exact distribution of the Pareto sum means taking, at given and for any , with i.i.d. r.v.'s with parent r.v. with finite th moment and pdf defined, for , by

Let us look at the three first moments of . The direct dependence is on (and ) and indirectly on since . We have

11.31

(note that , for any that we consider, and any ) and

11.32

using the expressions of and given in Theorem 11.2. A straightforward computation of the third centered moment of provides

11.33

where denotes the antiderivative of the function , that is, if ,

whereas, if ,

and, if ,

For simplicity, let us look at the case and consider the Berry–Esséen inequality. For , we would use the Edgeworth expansion, with similar arguments as developed later. Various authors have worked on this type of Berry–Esséen inequality, in particular to sharpen the accuracy of the constant appearing in it. In the case of Berry–Esséen bounds, the value of the constant factor has decreased from 7.59 by Esséen (1942) to 0.4785 by Tyurin (2010), to 0.4690 by Shevtsova (2013) in the i.i.d. case, and to 0.5600 in the general case. Note also that these past decades, much literature ((Stein, 1972, 1986); (Chen and Shao, 2004); (Cai, 2012); (Pinelis, 2013), etc.) has been dedicated to the generalization of this type of inequality, such as the remarkable contribution by Stein.

We can propose the following bound.

Note that depends on the two parameters and . We represent this function on a same plot for a given value of but for various values of , namely, , and 1000, respectively, to compare its behavior according to the parameter . Then we repeat the operation for different , namely, for , respectively (Figure 11.3).

We observe that the bound is an increasing then decreasing function of , with a maximum less than , which is decreasing with and . The -coordinate of the maximum is proportional to , with the proportion decreasing with . The interval on the -axis for which the error is larger than has a small amplitude, which is decreasing with .

We show in Table 11.3 the values of the coordinates of the maximum of computed on R for and (corresponding to aggregating weekly returns to obtain yearly returns), 100, 250 (corresponding to aggregating daily returns to obtain yearly returns), 500, 1000, respectively.

52	101	4.9	86	4.9	78	4.9
100	196	4.6	166	4.6	150	4.6
250	494	4.2	417	4.1	376	4.0
500	990	3.9	834	3.7	751	3.5
1000	1984	3.6	1667	3.3	1501	3.0

Hence the result of Proposition 11.2.

Indeed, we have

Note that the Berry–Esséen inequality has been proved by Petrov to hold also for probability density functions (see (Petrov, 1956) or (Petrov, 1995)). It has been refined by Shevtsova (2007), and we will use her result to evaluate . We need to go back to the pdf of the standardized sum of i.i.d. r.v.'s with pdf , which can be expressed as

It is straightforward to show by induction that

Then, since , we can write

11.36

Since we consider a sum of i.i.d. r.v.'s () with parent r.v. having a finite th moment, we obtain via (Petrov, 1956) and (Shevtsova, 2007) that there exists a constant such that

11.37

where is defined in (11.34).

Hence, combining (11.36) and (11.37) gives

from which we deduce that

As in the case (Proposition 11.2), this bound could be computed numerically.

11.4.3.1 Presentation of the study

We simulate with parent r.v. -Pareto distributed, with different sample sizes, varying from (corresponding to aggregating weekly returns to obtain yearly returns) through (corresponding to aggregating daily returns to obtain yearly returns) to representing a large size portfolio.

We consider different shape parameters, namely, , respectively. Recall that simulated Pareto r.v.'s 's ( can be obtained simulating a uniform r.v. on and then applying the transformation .

For each and each , we aggregate the realizations 's (). We repeat the operation times, thus obtaining realizations of the Pareto sum , from which we can estimate its quantiles.

Let denote the empirical quantile of order of the Pareto sum (associated with the empirical cdf and pdf ) defined by

Recall, for completeness, that the empirical quantile of converges to the true quantile as and has an asymptotic normal behavior, from which we deduce the following confidence interval at probability a for the true quantile: , where can be empirically estimated for such a large . We do not compute them numerically: being very large, bounds are close.

We compute the values of the quantiles of order , ( indicating the chosen method), obtained by the main methods, the GCLT method, the Max one, Normex, and the weighted normal method, respectively. We do it for various values of and . We compare them with the (empirical) quantile obtained via Pareto simulations (estimating the true quantile). For that, we introduce the approximative relative error:

We consider three possible order : , (threshold for Basel II) and (threshold for Solvency 2).

We use the software R to perform this numerical study, with different available packages. Let us particularly mention the use of the procedure Vegas in the package R2Cuba for the computation of the double integrals. This procedure turns out not to be always very stable for the most extreme quantiles, mainly for low values of . In practice, for the computation of integrals, we would advise to test various procedures in R2Cuba (Suave, Divonne, and Cuhre, besides Vegas) or to look for other packages. Another possibility would be implementing the algorithm using altogether a different software, as, for example, Python.

11.4.3.2 Estimation of the VaR with the various methods

All codes and results are obtained for various and are given in Kratz (2013) (available upon request) and will draw conclusions based on all the results.

We start with a first example when to illustrate our main focus, when looking at data under the presence of moderate heavy tail. We present here the case in Table 11.4.

		(%)
95%	103.23	104.35	102.60	103.17	109.25
		1.08	0.61	0.06	5.83
99%	119.08	111.67	117.25	119.11	118.57
		6.22	1.54	0.03	0.43
99.5%	128.66	114.35	127.07	131.5	121.98
		11.12	1.24	2.21	5.19
		(%)
95%	189.98	191.19	187.37	189.84	197.25
		0.63	1.38	0.07	3.83
99%	210.54	201.35	206.40	209.98	209.74
		4.36	1.96	0.27	0.38
99.5%	222.73	205.06	219.14	223.77	214.31
		7.93	1.61	0.47	3.78
		(%)	(%)	(%)
95%	454.76	455.44	446.53	453.92	464.28
		0.17	1.81	0.18	2.09
99%	484.48	471.5	473.99	483.27	483.83
		2.68	2.17	0.25	0.13
99.5%	501.02	477.38	492.38	501.31	490.98
		4.72	1.73	0.06	2.00
		(%)	(%)	(%)
95%	888.00	888.16	872.74	886.07	900.26
		0.02	1.72	0.22	1.38
99%	928.80	910.88	908.97	925.19	927.80
		1.93	2.14	0.39	0.11
99.5%	950.90	919.19	933.23	948.31	937.89
		3.33	1.86	0.27	1.37

Let us also illustrate in table 11.5 the heavy tail case, choosing, for instance, , which means that . Take, for example, respectively, to illustrate the fit of Normex even for small samples. Note that the weighted normal does not apply here since .

		(%)
95%	246.21	280.02	256.92	245.86
		13.73	4.35	0.14
99%	450.74	481.30	455.15	453.92
		6.78	0.97	0.71
99.5%	629.67	657.91	631.66	645.60
		4.48	0.31	2.53
		(%)
95%	442.41	491.79	456.06	443.08
		11.16	3.09	0.15
99%	757.82	803.05	762.61	761.66
		5.97	0.63	0.51
99.5%	1031.56	1076.18	1035.58	1032.15
		4.33	0.39	0.06

11.4.3.3 Discussion of the results

• Those numerical results are subject to numerical errors due to the finite sample of simulation of the theoretical value, as well as the choice of random generators, but the most important reason for numerical error of our methods resides in the convergence of the integration methods. Thus, one should read the results, even if reported with many significant digits, to a confidence we estimate to be around 0.1% (also for the empirical quantiles).
• The Max method overestimates for and underestimates for ; it improves a bit for higher quantiles and . It is a method that practitioners should think about when wanting to have a first idea on the range of the VaR, because it is very simple to use and costless in terms of computations. Then they should turn to Normex for an accurate estimation.
• The GCLT method () overestimates the quantiles but improves with higher quantiles and when increases.
• Concerning the CLT method (), we find out that:
  • The higher the quantile, the higher the underestimation; it improves slightly when increases, as expected.
  • The smaller , the larger the underestimation.
  • The estimation improves for smaller-order of the quantile and large , as expected, since, with smaller order, we are less in the upper tail.
  • The VaR estimated with the normal approximation is almost always lower than the VaR estimated via Normex or the weighted normal method. The lower and , the higher the difference with Normex.
  • The difference between the VaR estimated by the CLT and the one estimated with Normex appears large for relatively small , with a relative error going up to 13%, and decreases when becomes larger.
• With the Weighted Normal method (), it appears that:
  • The method overestimates the 95% quantile but is quite good for the 99% one. In general, the estimation of the two upper quantiles improve considerably when compared with the ones obtained via a straightforward application of the CLT.
  • The estimation of the quantiles improve with increasing n and increasing .
  • The results are generally not as sharp as the ones obtained via Normex, but better than the ones obtained with the max method, whenever .
• Concerning Normex, we find out that:
  • The accuracy of the results appears more or less independent of the sample size , which is the major advantage of our method when dealing with the issue of aggregation.
  • For , it always gives sharp results (error less than 0.5% and often extremely close); for most of them, the estimation is indiscernible from the empirical quantile, obviously better than the ones obtained with the other methods.
  • For , the results for the most extreme quantile are slightly less satisfactory than expected. We attribute this to a numerical instability in the integration procedure used in R. Indeed, for very large quantiles (%), the convergence of the integral seems a bit more unstable (due to the use of the package Vegas in R), which may explain why the accuracy decreases a bit and may sometimes be less than with the Max method. This issue should be settled using other R-package or software.
• We have concentrated our study on the VaR risk measure because it is the one used in solvency regulations both for banks and insurances. However, the expected shortfall, which is the only coherent measure in presence of fat tails, would be more appropriate for measuring the risk of the companies. The difference between the risk measure estimated by the CLT and the one estimated with Normex would certainly be much larger than what we obtain with the VaR when the risk is measured with the expected shortfall, pleading for using this measure in presence of fat tails.

11.4.3.4 Normex in practice

From the results we obtained, Normex appears as the best method among the ones we studied, applicable for any and . This comparison was done on simulated data. A next step will be to apply it on real data.

Let us sketch up a step-by-step procedure on how Normex might be used and interpreted in practice on real data when considering aggregated heavy-tailed risks.

We dispose of a sample , with unknown heavy-tailed cdf having positive tail index . We order the sample as and consider the aggregated risks that can be rewritten as .

1. Preliminary step: Estimation of , with standard EVT methods (e.g., Hill estimator (Hill, 1975)), -estimator (Kratz and Resnick, 1996); (Beirlant et al., 1996)), etc.); let be denoted as an estimate of .
2. Define with (see (11.13)); the largest order statistics share the property of having a th infinite moment, contrary to the first ones. Note that is independent of the aggregation size.
3. The first order statistics and the last ones being, conditionally on , independent, we apply the CLT to the sum of the first order statistics conditionally on and compute the distribution of the sum of the last ones conditionally on assuming a Pareto distribution for the r.v.'s because of (11.3).
4. We can then approximate the cdf of by defined in Theorem 11.2, which provides a sharp approximation, easily computable whatever the size of the sample is.
5. We deduce any quantile of order of as , which allows an accurate evaluation of risk measures of aggregated heavy-tailed risks.