Tail approximations for sums of dependent regularly varying random variables under Archimedean copula models

A well-known problem in applied probability is the evaluation of the probability that the sum of $n$ random variables (rvs) exceeds a certain level $s$. This problem finds its way in many areas of application such as actuarial science, finance, quantitative risk management and, reliability. Different methods can be used to tackle this problem depending on the type of distributions of the random variables and their interaction as well as the values of $n$ and $s$.

In this paper, we consider the case of $n$ positive heavy-tailed random variables that are linked through dependence structures based on Archimedean copulas. An $n$-dimensional copula $C$ is a multivariate distribution on $\left[0,1\right]^{n}$ with uniform margins. Following [Ling, 1965], any Archimedean copula can be simply written as

$$C(u_{1},...,u_{n})=\Phi\left(\Phi^{\leftarrow}(u_{1})+...+\Phi^{\leftarrow}(u_{n})\right),\text{ }(u_{1},...,u_{n})\in\left[0,1\right]^{n},$$

(1)

where $\Phi$ is a non-increasing function referred to as the generator of the copula $C$ and $\Phi^{\leftarrow}$ is the generalized inverse function of $\Phi$ defined as $\Phi^{\leftarrow}(x)=\inf\{t\in\mathbb{R}^{+}:\Phi(t)\leq x\}$. The conditions under which a generator $\Phi$ defines a proper $n$-dimensional copula are given in detail in [McNeil and Nešlehová, 2009].

Let $\mathbf{U}=\left(U_{1},...,U_{n}\right)$ be a random vector distributed as an Archimedean copula $C$. For $i=1,\ldots,n$, we assume that the survival distribution function of the $i$-th rv of the sum is regularly varying, i.e. it satifies $\bar{F}_{i}(x)=1-F_{i}(x)=x^{-\alpha_{i}}l_{i}(x)$ where $\alpha_{i}>0$ and $l_{i}$ is a slowly varying function at infinity. We assume without loss of generality that $\alpha_{1}\leq\ldots\leq\alpha_{n}$. Throughout the paper we shall consider two sums that are linked through the Archimedean copula $C$ in the following way

$$S_{n}^{X}=\sum_{i=1}^{n}X_{i}\quad\text{and}\quad S_{n}^{Y}=\sum_{i=1}^{n}Y_{i}\text{,}$$

where

$$\mathbf{X}=(X_{1},...,X_{n})\overset{d}{=}(F_{1}^{\leftarrow}(U_{1}),...,F_{n}^{\leftarrow}(U_{n}))\quad\text{and}\quad\mathbf{Y}=(Y_{1},...,Y_{n})\overset{d}{=}\left(\overline{F}_{1}^{\leftarrow}(U_{1}),...,\overline{F}_{n}^{\leftarrow}(U_{n})\right),$$

with $F_{i}^{\leftarrow}(u_{i})=\inf\{t\in\mathbb{R}^{+}:F_{i}(t)\geq x\}$ since $F_{i}$ is a non-decreasing function. Note that the multivariate cumulative distribution function of $\mathbf{X}$ is given by

$$\Pr(X_{1}\leq x_{1},...,X_{n}\leq x_{n})=C(F_{1}(x_{1}),...,F_{n}(x_{n})),\text{ }(x_{1},...,x_{n})\in\mathbb{R}^{n},$$

while the multivariate survival distribution function of $\mathbf{Y}$ is given by

$$\Pr(Y_{1}>y_{1},...,Y_{n}>y_{n})=C\left(\overline{F}_{1}(y_{1}),...,\overline{F}_{n}(y_{n})\right),\text{ }(y_{1},...,y_{n})\in\mathbb{R}^{n}.$$

The Archimedean copula $C$ is referred to as the copula of $\mathbf{X}$ and as the survival copula of $\mathbf{Y}$.

To approximate the probabilities $z_{X}(s)=\Pr\left(S_{n}^{X}>s\right)$ and $z_{Y}(s)=\Pr\left(S_{n}^{Y}>s\right)$ when $s$ is large, one could use functions that are asymptotically equivalent to these probabilities. However there are few results concerning the asymptotic behaviours of these probabilities. Actually it strongly depends on the tails of the Archimedean copulas (see [Charpentier and Segers, 2009] for a fine analysis of the several types of tails based on characteristics of the Archimedean generator) and strong assumptions have to hold to characterize these behaviours. For example, if the upper-tails of the Archimedean copula $C$ are independent, i.e $\lim_{u\rightarrow 1}\Pr(U_{i}>u|U_{j}>u)=0$ for all $i\neq j$ and if there exist $n-1$ non-negative constants $c_{2},...,c_{n}$ such that $c_{i}=\lim_{s\rightarrow\infty}\bar{F}_{i}(s)/\bar{F}_{1}(s)$, then it can be shown that

$$\lim\limits_{s\rightarrow\infty}\frac{\Pr(S_{n}^{X}>s)}{\overline{F}_{1}(s)}=1+\sum_{i=2}^{n}c_{i}.$$

If the lower-tails of the Archimedean copula $C$ are rather independent, i.e. $\lim_{u\rightarrow 0}\Pr(U_{i}<u|U_{j}<u)=0$ for all $i\neq j$, then

$$\lim\limits_{s\rightarrow\infty}\frac{\Pr(S_{n}^{Y}>s)}{\overline{F}_{1}(s)}=1+\sum_{i=2}^{n}c_{i}$$

(see e.g. [Jessen and H., 2006] or [Yuen and Yin, 2012]). [Sun and Li, 2010] studied the asymptotic behaviours of $z_{X}(s)$ and $z_{Y}(s)$ under the assumption of identically distributed marginals and specific upper or lower-tail dependence. Let $\beta>0$ and $l_{\Phi}$ be a slowly varying function at infinity. If $1-\Phi(x)=x^{\beta}l_{\Phi}\left(x^{-1}\right)$, they proved that

$$\lim\limits_{s\rightarrow\infty}\frac{\Pr(S_{n}^{X}>s)}{\Pr(X_{1}>s)}=q_{n}^{C}(\alpha,\beta),$$

where $\alpha$ denotes the common tail index and

$$q_{n}^{C}(\alpha,\beta)=\int_{\sum_{i=1}^{n}v_{i}^{-1/\alpha}>1}\frac{\partial^{n}}{\partial v_{1}...\partial v_{n}}\sum_{1\leq i_{1},...,i_{j}\leq n}\left((-1)^{j-1}(v_{i_{1}}^{1/\beta}+...+v_{i_{j}}^{1/\beta})^{\beta}\right)\mathrm{d}v_{1}...\mathrm{d}v_{n}.$$

If the generator rather satisfies $\Phi(x)=x^{-\beta}l_{\Phi}\left(x\right)$, they proved that

$$\lim\limits_{s\rightarrow\infty}\frac{\Pr(S_{n}^{Y}>s)}{\Pr(Y_{1}>s)}=q_{n}^{D}(\alpha,\beta),$$

where

$$q_{n}^{D}(\alpha,\beta)=\int_{\sum_{i=1}^{n}v_{i}^{-1}>1}\frac{\partial^{n}}{\partial v_{1}...\partial v_{n}}\left(v_{1}^{-\alpha/\beta}+...+v_{n}^{-\alpha/\beta}\right)^{-\beta}\mathrm{d}v_{1}...\mathrm{d}v_{n},$$

(see also [Wüthrich, 2003]). Although $q_{n}^{C}(\alpha,\beta)$ and $q_{n}^{D}(\alpha,\beta)$ are known, they do not have closed-form expressions and they can not be easily computed.

In this paper, we aim to provide two numerical methods for approximating $z_{X}(s)$ and $z_{Y}(s)$ for different values of $n$ and $s$ and choices of parameters and functions: $((\alpha_{1},l_{1}),\ldots,(\alpha_{n},l_{n}))$ and $\Phi$.

Our first method is based on conditional Monte Carlo and is ideally suited when $s$ and/or $n$ are large. The classical Monte Carlo method is easy to implement and can be applied in complex situations such as high dimensional calculations. However, it is well known that it is inadequate for small probability simulation since the relative errors (variation coefficients) are too large. [Asmussen and Glynn, 2007] introduced relative error as a measure of efficiency of an estimator and several definitions of efficient estimators. An unbiased estimator $Z(s)$ of the probability $z(s)$, with relative error $e(Z(s))=\sqrt{\mathbb{E[}Z^{2}(s)]}/z(s)$, is called (i) a logarithmically efficient estimator if $\lim\sup_{s\rightarrow\infty}e(Z(s))\ [z(s)]^{\epsilon}=0$ for all $\epsilon>0$; (ii) an estimator with bounded relative error if $\lim\sup_{s\rightarrow\infty}e(Z(s))<\infty$; (iii) an estimator with vanishing relative error $\lim\sup_{s\rightarrow\infty}e(Z(s))=0$.

For sums of independent random variables, the most widely used alternatives to crude Monte Carlo computation of rare-event probabilities are conditional Monte Carlo and importance sampling. [Asmussen and Binswanger, 1997] propose a logarithmically efficient algorithm based on conditional Monte Carlo simulation using order statistics.

[Boots and Shahabuddin, 2001] use importance sampling to simulate ruin probabilities for subexponential claims and [Juneja and Shahabuddin, 2002] use importance sampling based on hazard rate twisting to simulate heavy-tailed processes. [Asmussen and Kroese, 2006] propose two algorithms which use importance sampling and conditional Monte carlo and study their efficiency in the Pareto and Weibull case.

Estimating tail distribution of the sums of dependent random variables via simulation requires a specific expression for the dependence structure or a closed form expression for the conditional distribution functions, the case of elliptic distributions is an example. For an elliptic dependence structure, [Blanchet and Rojas-Nandayapa, 2011] proposed a conditional Monte Carlo estimator for the tail distribution of the sum of log-elliptic random variables and proved that it has a logarithmically efficient relative error. The sum of the log-elliptic random variables was also estimated by [Kortschak and Hashorva, 2013] using the simulation method introduced by [Asmussen and Kroese, 2006] and favorable results are presented especially in the multivariate lognormal case. [Asmussen et al., 2011] and [Blanchet et al., 2008] focus on the efficient estimation of sums of correlated lognormals using importance sampling and conditional Monte Carlo strategies. [Chan and Kroese, 2010], [Chan and Kroese, 2011] use conditional Monte Carlo notably in a credit risk setting under the t-copula model to estimate rare-event probabilities.

In this paper, we introduce four different estimators of the probabilities $z_{X}(s)=\Pr\left(S_{n}^{X}>s\right)$ and $z_{Y}(s)=\Pr\left(S_{n}^{Y}>s\right)$ using techniques of conditional Monte Carlo simulation. The main idea to build our estimators is to first isolate the known probabilities $\Pr(M_{n}^{X}>s)$ or $\Pr(M_{n}^{Y}>s)$ where $M_{n}^{A}$ correspond to the maximum element of a given vector $\mathbf{A}$ (because $\Pr(M_{n}^{X}>s)$ or $\Pr(M_{n}^{Y}>s)$ have closed-form expressions in our framework), and then simulate conditionally on the values taken by these maxima. Two effective simulation techniques of vectors of Archimedean copula proposed in [Brechmann et al., 2013] and [McNeil and Nešlehová, 2009] will be used. We show that most of our estimators have bounded relative errors.

Our second method is based on analytical expressions of the survival multivariate distribution function and provides sharp, deterministic and numerical bounds of the probabilities using the same ideas as developed in [Cossette et al., 2014]. This method performs very well for cases when $n$ is relatively small and effectively completes the conditional Monte Carlo method.

The outline of the paper is as follows. In Section 2, we present the two simulation techniques related to Archimedean copulas which are later used to develop the proposed estimators. These estimators are introduced, described and discussed in Section 3. Some results on their asymptotic efficiency are also given. Section 4 explains how to derive and compute the numerical bounds following the approach proposed in [Cossette et al., 2014]. Section 5 illustrates the accuracy of both methods through a numerical study and discusses implementation issues.

The classical simulation method for a dependent vector relies on conditional distributions. Consider a random vector $\mathbf{U}$ with density function $c(u_{1},...,u_{n})$ which can be decomposed as the product of conditional densities

$$c(u_{1},...,u_{n})=c_{n|n-1,\ldots,1}(u_{n}|u_{n-1},...,u_{1})...c_{2|1}(u_{2}|u_{1})c_{1}(u_{1}).$$

The classical procedure of simulating vector such a vector $\mathbf{U}$ is then: simulate $u_{1}$ based on $c_{1}(u_{1})$, simulate $U_{2}$ based on $c_{2|1}(u_{2}|u_{1})$, …, simulate $U_{n}$ based on $c_{n|n-1,...,1}(u_{n}|u_{n-1},...,u_{1})$. Hence, the realization of vector $\mathbf{U}$ is created by calculating $(n-1)$ times the inverses of conditional distribution functions. However it can be difficult and take quite an amount of time when the distribution function of $\mathbf{U}$ is an Archimedean copula.

Another method for an Archimedean copula could be to consider the mixed exponential or frailty representation often used to model dependent lifetimes and discussed in, notably, [Marshall and Olkin, 1988], [McNeil, 2008] and [Hofert, 2008]. In this case, the Archimedean generator $\Phi$ is the Laplace-Stieltjes transform of a non-negative random variable. Such a method thus requires to invert the Laplace-Stieltjes transform $\Phi$ which can not always be evaluated explicitly.

To circumvent these problems, one can resort to two effective simulation techniques proposed in [Brechmann et al., 2013] and [McNeil and Nešlehová, 2009]. More precisely, [Brechmann et al., 2013] use the Kendall distribution function while [McNeil and Nešlehová, 2009] suggest a simulation method which relies on $\ell_{1}$-norm symmetric distributions.

In this section, we propose four different estimators of $z_{X}(s)$ and $z_{Y}(s)$. All estimators rely on a similar idea which is to decompose the probability of interest into different components. The known components are exactly evaluated and the other ones, which are our main concern, are estimated by simulation.

To begin, we need to establish basic notations. Throughout, $\mathbf{X}_{-i}$ denotes the vector $\mathbf{X}=(X_{1},...,X_{n})$ with the $i$-th component $X_{i}$ removed and $M_{i}^{X}$ (or $M_{i}\{\mathbf{X}\}$) corresponds to the $i^{th}$ element of the vector $\mathbf{X}$ after rearranging the elements of $\mathbf{X}$ in a non-decreasing order. Obviously, $M_{1}^{X}$ and $M_{n}^{X}$ correspond to the minimum and maximum element of vector $\mathbf{X}$, respectively. The same convention holds for $\mathbf{Y}=(Y_{1},...,Y_{n})$.

In what follows, we encounter frequently the evaluation of the probabilities $\Pr(M_{n}^{X}>s)$ and $\Pr\left(M_{n}^{Y}>s\right)$ which rely on the marginal distributions $F_{1}$,…,$F_{n}$ and the Archimedean generator $\Phi$. They are obtained as follows

$$\Pr(M_{n}^{X}>s)=1-C\left(F_{1}(s),...,F_{n}(s)\right)$$

and

$$\Pr(M_{n}^{Y}>s)=1-\overline{C}\left(\overline{F}_{1}(s),...,\overline{F}_{n}(s)\right),$$

with $\overline{C}\left(u_{1},...,u_{n}\right)=\Pr\left(U_{1}>u_{1},...,U_{n}>u_{n}\right)$.