On the Efficient Simulation of the Left-Tail of the Sum of Correlated Log-normal Variates

TL;DR

This paper introduces an efficient importance sampling estimator for accurately simulating the left-tail probabilities of the sum of correlated Log-normal variates, a problem relevant in wireless and financial applications.

Contribution

It presents a novel estimator combining importance sampling and control variates that achieves asymptotically vanishing relative error for the left-tail of correlated Log-normal sums.

Findings

The estimator outperforms existing methods in simulation accuracy.

It achieves asymptotically vanishing relative error in tail probability estimation.

Simulation results confirm improved efficiency and accuracy.

Abstract

The sum of Log-normal variates is encountered in many challenging applications such as in performance analysis of wireless communication systems and in financial engineering. Several approximation methods have been developed in the literature, the accuracy of which is not ensured in the tail regions. These regions are of primordial interest wherein small probability values have to be evaluated with high precision. Variance reduction techniques are known to yield accurate, yet efficient, estimates of small probability values. Most of the existing approaches, however, have considered the problem of estimating the right-tail of the sum of Log-normal random variables (RVS). In the present work, we consider instead the estimation of the left-tail of the sum of correlated Log-normal variates with Gaussian copula under a mild assumption on the covariance matrix. We propose an estimator combining an existing mean-shifting importance sampling approach with a control variate technique. The main result is that the proposed estimator has an asymptotically vanishing relative error which represents a major finding in the context of the left-tail simulation of the sum of Log-normal RVs. Finally, we assess by various simulation results the performances of the proposed estimator compared to existing estimators.

Equations

$$f(\mathbf{y})=\frac{\exp \left(-\frac{1}{2}{\left(\mathbf{y}-\mathbf{\mu }\right)}^t{\Sigma }^{-1}\left(\mathbf{y}-\mathbf{\mu }\right)\right)}{\sqrt{(2\pi {)}^N|\Sigma |}},$$

$$\underset{{\Delta }_N}{\min }{\mathbf{w}}^t\Sigma \mathbf{w}$$

$$\alpha ({\gamma }_{th})=P_f\left(\sum _{i=1}^N{X}_i\le {\gamma }_{th}\right).$$

$$\alpha ({\gamma }_{th})$$

$$={\mathbb{E}}_g\left[1_{\left({\sum }_{i=1}^N\exp \left(Y_i\right)\le {\gamma }_{th}\right)}L\left(Y_1,Y_2,\cdots ,Y_N\right)\right],$$

$$L\left(Y_1,Y_2,\cdots ,Y_N\right)=\frac{f\left(Y_1,Y_2,\cdots ,Y_N\right)}{g\left(Y_1,Y_2,\cdots ,Y_N\right)}.$$

$${\hat{\alpha }}_{IS}({\gamma }_{th})=\frac{1}{M}\sum _{k=1}^M{1}_{\left({\sum }_{i=1}^N\exp \left(Y_i^{(k)}\right)\le {\gamma }_{th}\right)}L\left(Y_1^{(k)},Y_2^{(k)},\cdots ,Y_N^{(k)}\right),$$

$$T_{{\gamma }_{th}}^{{}^{\prime }}\left(\mathbb{Y}\right)=T_{{\gamma }_{th}}\left(\mathbb{Y}\right)+\beta \left(Z_{{\gamma }_{th}}\left(\mathbb{Y}\right)-P({\gamma }_{th})\right),$$

$${\mathrm{var}}_g\left[T_{{\gamma }_{th}}^{{}^{\prime }}\left(\mathbb{Y}\right)\right]$$

$$+{\beta }^2{\mathrm{var}}_g\left[Z_{{\gamma }_{th}}\left(\mathbb{Y}\right)\right].$$

$${\beta }^{\ast }=-\frac{{\mathrm{cov}}_g\left[T_{{\gamma }_{th}}\left(\mathbb{Y}\right),Z_{{\gamma }_{th}}\left(\mathbb{Y}\right)\right]}{{\mathrm{var}}_g\left[Z_{{\gamma }_{th}}\left(\mathbb{Y}\right)\right]}.$$

$${\mathrm{var}}_g\left[T_{{\gamma }_{th}}^{{}^{\prime }}\left(\mathbb{Y}\right)\right]=\left(1-{\rho }_{T_{{\gamma }_{th}}\left(\mathbb{Y}\right),Z_{{\gamma }_{th}}\left(\mathbb{Y}\right)}^2\right){\mathrm{var}}_g\left[T_{{\gamma }_{th}}\left(\mathbb{Y}\right)\right],$$

$${\hat{\alpha }}_{IS-CV}({\gamma }_{th})=\frac{1}{M}\sum _{k=1}^M\left(T_{{\gamma }_{th}}\left({\mathbb{Y}}^{(\mathbb{k})}\right)+\beta \left(Z_{{\gamma }_{th}}\left({\mathbb{Y}}^{(\mathbb{k})}\right)-P({\gamma }_{th})\right)\right).$$

$${\mathbb{E}}_g\left[T_{{\gamma }_{th}}^{{}^{\prime }2}\left(\mathbb{Y}\right)\right]\ge {\alpha }^2({\gamma }_{th}).$$

$$\frac{\log \left({\mathbb{E}}_g\left[T_{{\gamma }_{th}}^{{}^{\prime }2}\left(\mathbb{Y}\right)\right]\right)}{\log \left(\alpha ({\gamma }_{th})\right)}\le 2.$$

$$\underset{{\gamma }_{th}\to 0}{\lim }\frac{\log \left({\mathbb{E}}_g\left[T_{{\gamma }_{th}}^{{}^{\prime }2}\left(\mathbb{Y}\right)\right]\right)}{\log \left(\alpha ({\gamma }_{th})\right)}=2.$$

$$\underset{{\gamma }_{th}\to 0}{\lim }\frac{{\mathrm{var}}_g\left[T_{{\gamma }_{th}}^{{}^{\prime }}\left(\mathbb{Y}\right)\right]}{{\alpha }^{2-ϵ}({\gamma }_{th})}=0.$$

$$\underset{{\gamma }_{th}\to 0}{\mathrm{limsup}}\frac{{\mathrm{var}}_g\left[T_{{\gamma }_{th}}^{{}^{\prime }}\left(\mathbb{Y}\right)\right]}{{\alpha }^2({\gamma }_{th})}<+\infty .$$

$$\underset{{\gamma }_{th}\to 0}{\mathrm{limsup}}\frac{{\mathrm{var}}_g\left[T_{{\gamma }_{th}}^{{}^{\prime }}\left(\mathbb{Y}\right)\right]}{{\alpha }^2({\gamma }_{th})}=0.$$

$$g(\mathbf{y})=\frac{\exp \left(-\frac{1}{2}{\left(\mathbf{y}-\mathbf{\mu }-\Lambda \right)}^t{\Sigma }^{-1}\left(\mathbf{y}-\mathbf{\mu }-\Lambda \right)\right)}{\sqrt{(2\pi {)}^N|\Sigma |}}.$$

$$L\left(Y_1,Y_2,\cdots ,Y_N\right)=\exp \left(-{\Lambda }^t{\Sigma }^{-1}(\mathbb{Y}-\mathbf{\mu })+\frac{1}{2}{\Lambda }^t{\Sigma }^{-1}\Lambda \right).$$

$$\alpha ({\gamma }_{th})$$

$$={\mathbb{E}}_f\left[\exp \left(-{\Lambda }^t{\Sigma }^{-1}(\mathbb{Y}-\mathbf{\mu })-\frac{1}{2}{\Lambda }^t{\Sigma }^{-1}\Lambda \right)1_{\left({\sum }_{i=1}^N\exp \left(Y_i+{\Lambda }_i\right)\le {\gamma }_{th}\right)}\right].$$

$${\mathbb{E}}_g\left[T_{{\gamma }_{th}}^2\left(\mathbb{Y}\right)\right]=\exp \left({\Lambda }^t{\Sigma }^{-1}\Lambda \right)P_f\left(\sum _{i=1}^N\exp \left(Y_i-{\Lambda }_i\right)\le {\gamma }_{th}\right).$$

$${\Lambda }_k^{\ast }=\sum _{i,j=1}^{\bar{n}}{\Sigma }_{ki}{\bar{\Sigma }}_{ij}^{-1}\left(\log \left(x\right)-\log \left(\frac{{\bar{A}}_1+\cdots +{\bar{A}}_{\bar{n}}}{{\bar{A}}_j}\right)-{\bar{\mu }}_j\right).$$

$${\mathbb{E}}_g\left[T_{{\gamma }_{th}}^2\left(\mathbb{Y}\right)\right]\le C{\alpha }^2({\gamma }_{th}){\left(\log \left(\frac{1}{{\gamma }_{th}}\right)\right)}^{\bar{n}},$$

$$\alpha ({\gamma }_{th})$$

$$\underset{{\gamma }_{th}\to 0}{\sim }\frac{\sqrt{{\Sigma }_{11}}}{\sqrt{2\pi }\log \left(\frac{1}{{\gamma }_{th}}\right)}\exp \left(-\frac{{\left(\log \left({\gamma }_{th}\right)-{\mu }_1\right)}^2}{2{\Sigma }_{11}}\right).$$

$${\bar{w}}_1$$

$${\bar{w}}_j$$

$${\Lambda }_1^{\ast }$$

$${\Lambda }_k^{\ast }$$

$$Z_{{\gamma }_{th}}\left(\mathbb{Y}\right)=1_{\left(\exp (Y_1)\le {\gamma }_{th}\right)}L\left(Y_1,Y_2,\cdots ,Y_N\right),$$

Full text

**Keywords: **Sum of correlated Log-normal, small probability values, variance reduction techniques, left-tail of the sum of correlated Log-normal variates, importance sampling, control variate, asymptotically vanishing relative error.

1 Introduction

The Log-normal distribution is encountered in many applications such as in financial engineering [13]. In performance analysis of wireless communications systems, it has shown a good fit to realistic propagation channels [24, 15, 20]. Therefore, investigating the distribution of sums of Log-normal random variables (RVs) is of primordial practical interest. In fact, the pricing of Asian or basket options is closely related to the distribution of the sum of Log-normal variates [13]. A further application where the problem of finding the distribution of sums of Log-normal variates is encountered is in the evaluation of the value at risk, defined as (1-$\alpha$) quantile of the loss distribution, for a sufficiently small value of $\alpha$ [4]. In wireless communication systems, the distribution of sums of Log-normal RVs is of major practical interest in the problem of evaluating the outage probability at the output of receivers with diversity techniques [9].

The distribution of the sum of Log-normal variates is not known in a closed form. This has lead researchers to propose various approximation techniques such as in [14, 23, 11, 12, 7]. However, the accuracy of these approximations is not ensured for all scenarios especially in the tail regions, i.e. the region we are interested in, as illustrated in [7].

Variance reduction Monte Carlo (MC) methods constitute good alternatives to efficiently estimate tail probabilities of the distribution of sums of Log-normal RVs. The literature of estimating the distribution of sums of Log-normal variates is very rich. However, most of them have concentrated in the right-tail region instead of the left-tail one, which is the region of interest in this work. For instance, a non exhaustive list of efficient variance reduction techniques have been developed in [5, 8, 18] to estimate the right-tail of the distribution of sums of independent Log-normal RVs. In the correlated setting as well, efficient simulation approaches of the right-tail of sums of correlated Log-normal variates have been proposed in [1, 10, 21].

The estimation of the distribution of the left-tail of sums of Log-normal RVs has received less interest than that of the right-tail region. Estimating the left-tail of the distribution of sums of Log-normal variates is motivated by the problem of evaluating outage probabilities in the performance analysis of wireless communication systems operating over a Log-normal fading environment. It is only recently that researchers have accorded some interest to the estimation of the left-tail region [9, 4, 17]. In fact, the authors in [4] have used the well-known importance sampling (IS) technique; namely the exponential twisting approach [22], and have proved that the proposed estimator achieves the asymptotic optimality property under the independent and identically distributed (i.i.d) assumption. However, given that the application of the exponential twisting technique requires the knowledge of the moment generating function (MGF) which is out of reach for the Log-normal distribution, the authors in [4] have considered instead an estimator of the MGF [3]. In [9], two unified IS approaches have been proposed using the well-known hazard rate twisting technique [8, 18]. In particular, for the Log-normal setting, the first estimator was shown to achieve the asymptotic optimality criterion under the assumption of independent and not necessarily identically distributed sums of Log-normal RVs. The asymptotic optimality property holds again using the second IS scheme under the i.i.d assumption. Comparisons between the estimator of [4] and the two latter estimators have been performed in [9]. The sole work the authors are aware of it that deals with the efficient simulation of the left-tail of the sum of correlated Log-normal variates is in [17]. In fact, an IS estimator has been proposed by shifting the mean of the corresponding multivariate normal distribution. This estimator was shown to achieve the asymptotic optimality property under a mild assumption.

We aim in the present paper to further improve the mean-shifting IS estimator of [17]. More precisely, we consider the left-tail simulation of sums of correlated Log-normal variates, that is we aim to estimate the probability that a sum of correlated Log-normal variates is less than a sufficiently small threshold, and we provide an improved estimator of the IS scheme proposed in [17]. Our methodology is based on the asymptotic behavior results provided in [17] and considers the case where the left-tail of the distribution of the sum is determined by only one dominant component. In this setting, we improve the IS estimator of [17] by combining it with a control variate type of variance reduction technique. The introduced control variate is function of the dominant component that characterizes the tail behavior of the left-tail of the sum. The main result of the present work is that the improved estimator using the proposed control variate technique has the asymptotically vanishing relative error property which is the most desired property in the context of rare event simulations. Such a result represents a relevant contribution for the problem of estimating the left-tail of sums of Log-normal RVs since as it was mentioned above the existing estimators in the corresponding literature were only proved to achieve a weaker property; namely the asymptotic optimality criterion [9, 4, 17]. Simulation results show that our proposed approach yields a substantial amount of variance reduction compared to the mean-shifting IS approach of [17], a reduction that increases as we decrease the probability of interest, i.e. as we decrease the threshold value to zero.

The rest of the paper is organized as follows. In section II, we describe the problem setting and we review the concepts of IS and control variate. Section III is devoted to presenting the IS scheme of [17]. The main idea of the paper is provided in Section IV where we show how to improve the IS approach of [17] through the use of control variate. In the same section, the main result proving the asymptotically vanishing relative error property is provided. Finally, some selected simulation results are provided in Section V to assess the performance of the proposed approach which achieves a substantial amount of computational gain over the mean shifting IS of [17].

2 Problem Setting

We consider a random vector $\mathbb{Y}=(Y_{1},Y_{2},\cdots,Y_{N})^{t}$ with $N$-dimensional multivariate normal distribution with mean vector $\boldsymbol{\mu}$ and covariance matrix $\mathbb{\Sigma}$. We assume that $\mathbb{\Sigma}$ is positive definite. The probability density function (PDF) of the random vector $\mathbb{Y}$ is given as follows

[TABLE]

where $|\mathbb{\Sigma}|$ denotes the determinant of the matrix $\mathbb{\Sigma}$ and $\boldsymbol{y}^{t}$ is the transpose of the vector $\boldsymbol{y}$. Similarly to [17], we define the following quantities. The elements of $\mathbb{\Sigma}$ and its inverse $\mathbb{\Sigma}^{-1}$ will be denoted by $\Sigma_{ij}$ and $\Sigma_{ij}^{-1}$, respectively. Moreover, we define $A_{k}=\sum_{j=1}^{N}{\Sigma_{kj}^{-1}}$. Let $\bar{\boldsymbol{w}}$ be the unique solution of

[TABLE]

with $\Delta_{N}=\{\boldsymbol{w}\in\mathbb{R}^{N}\text{ such that }w_{i}\geq 0,i=1,2,\cdots,N,\text{ and}\sum_{i=1}^{N}{w_{i}}=1\}$. Let us denote by $\bar{n}$ the cardinal of the set $\{1,2,\cdots,N,\text{ such that }\bar{w}_{i}\neq 0\}$ and $\bar{I}=\{\bar{k}(1),\cdots,\bar{k}(\bar{n})\}$ the corresponding indexes. Then, we denote by $\bar{\boldsymbol{\mu}}\in\mathbb{R}^{\bar{n}}$ the vector with entries $\bar{\mu}_{i}=\mu_{\bar{k}(i)}$ and $\mathbb{\bar{\Sigma}}$ the $\bar{n}\times\bar{n}$ matrix with elements $\bar{\Sigma}_{ij}=\Sigma_{\bar{k}(i)\bar{k}(j)}$. Moreover, we denote by $\bar{\mathbb{\Sigma}}^{-1}$ the inverse of $\bar{\mathbb{\Sigma}}$ with elements $\bar{\Sigma}_{ij}^{-1}$. Finally, the row sums of $\bar{\mathbb{\Sigma}}^{-1}$ is denoted by $\bar{A_{k}}=\sum_{j=1}^{\bar{n}}{\bar{\Sigma}_{kj}^{-1}}$. As it was mentioned in [17], we assume without loss of generality that $\bar{I}=\{1,2,\cdots,\bar{n}\}$.

Let us consider the Log-normal random vector $\mathbb{X}=(X_{1},X_{2},\cdots,X_{N})^{t}$ such that $Y_{i}=\log(X_{i})$, $i=1,2,\cdots,N$. Our aim is to efficiently estimate the left-tail of the sum of $X_{i}$’s, that is the probability that the sum of correlated Log-normal RVs with Gaussian copula falls below a sufficiently small threshold:

[TABLE]

As it was mentioned previously, the above quantity is of paramount practical interest in the performance evaluation of wireless communication systems as it corresponds to the probability that the communication system is in outage, i.e. the probability that the system fails to operate correctly. The simplest method to estimate $\alpha(\gamma_{th})$ is to use the naive MC estimator. However, it is well-known that the naive estimator is computationally expensive when estimating rare events, i.e. region of sufficiently small values of $\alpha(\gamma_{th})$ which is the region of interest of the present work.

To construct computationally efficient estimators, variance reduction techniques represent good alternatives [19]. When appropriately employed, variance reduction techniques are known to yield accurate estimate with a smaller computational effort than the naive MC sampler. In the present work, we will consider two instances of variance reduction techniques which are IS and control variate.

2.1 Review of IS

IS has been extensively employed for the efficient simulation of rare events [8, 18, 4, 1]. The main idea of IS is to introduce an IS distribution $g(\cdot)$ under which sampling is performed instead of the original PDF $f(\cdot)$. In fact, the probability of interest $\alpha(\gamma_{th})$ could be re-written as follows

[TABLE]

where $\mathbb{1}_{\left(\cdot\right)}$ denotes the indicator function, whereas $\mathbb{E}_{f}\left[\cdot\right]$ and $\mathbb{E}_{g}\left[\cdot\right]$ are the expectation operators under the PDFs $f(\cdot)$ and $g(\cdot)$, respectively. $L(\cdot)$ is the likelihood ratio which is defined as

[TABLE]

Following the above probability change of measure, the IS estimator is defined as:

[TABLE]

where $\{(Y_{1}^{(k)},Y_{2}^{(k)},\cdots,Y_{N}^{(k)})^{t}\}_{k=1}^{M}$ are independent realizations of the random vector $\mathbb{Y}$ under the PDF $g(\cdot)$. The crucial step when using IS is the choice of the IS distribution $g(\cdot)$ that results in a variance reduction. In fact, a good IS distribution encourages the sampling of important realizations, that is samples that belong to the rare set $\{\sum_{i=1}^{N}{\exp\left(Y_{i}\right)}\leq\gamma_{th}\}$, and also tries to maintain the likelihood ratio $L(\cdot)$ constant in the rare region. In the following section, we will review the IS technique developed in [17] for the estimation of $\alpha(\gamma_{th})$.

2.2 Review of Control Variate

Control variate is also a variance reduction technique that can, if adequately used, yield a substantial amount of variance reduction. This method could be combined with IS in order to further reduce the variance. Suppose we want to estimate $\alpha=\mathbb{E}_{g}\left[T_{\gamma_{th}}\left(\mathbb{Y}\right)\right]$ where $T_{\gamma_{th}}\left(\mathbb{Y}\right)=\mathbb{1}_{\left(\sum_{i=1}^{N}{\exp\left(Y_{i}\right)}\leq\gamma_{th}\right)}L\left(Y_{1},\cdots,Y_{N}\right)$ is an IS estimator. Let $Z_{\gamma_{th}}\left(\mathbb{Y}\right)$ be a control variate with a known expected value $P(\gamma_{th})$. The idea of control variate technique is to consider the following estimator of $\alpha(\gamma_{th})$

[TABLE]

where $\beta\in\mathbb{R}$. Obviously, the above estimator is an unbiased estimator of $\alpha(\gamma_{th})$. The variance of $T^{{}^{\prime}}_{\gamma_{th}}\left(\mathbb{Y}\right)$ is given by

[TABLE]

Hence, the optimal value of $\beta$ is the value that minimizes the variance of $T^{{}^{\prime}}_{\gamma_{th}}\left(\mathbb{Y}\right)$ and given, through a simple computation, by

[TABLE]

By plugging the optimal value $\beta^{*}$ into (2.2), we easily get that the variance of $T^{{}^{\prime}}_{\gamma_{th}}\left(\mathbb{Y}\right)$ is given by

[TABLE]

where $\rho_{T_{\gamma_{th}}\left(\mathbb{Y}\right),Z_{\gamma_{th}}\left(\mathbb{Y}\right)}$ is the correlation coefficient between the two RVs $T_{\gamma_{th}}\left(\mathbb{Y}\right)$ and $Z_{\gamma_{th}}\left(\mathbb{Y}\right)$. From the above result, we conclude that, in order to further reduce the variance of $T_{\gamma_{th}}\left(\mathbb{Y}\right)$, the RV $Z_{\gamma_{th}}\left(\mathbb{Y}\right)$ has to be selected such that it is highly correlated with the RV $T_{\gamma_{th}}\left(\mathbb{Y}\right)$. In the next section, we will describe how the RV $Z_{\gamma_{th}}\left(\mathbb{Y}\right)$ is selected in order to achieve a substantial amount of variance reduction. The estimator of $\alpha(\gamma_{th})$ following the control variate technique with any value of $\beta$ is as follows

[TABLE]

Note that the optimal value $\beta^{*}$ is generally unknown and has to be estimated via sample covariance and sample variance estimators. However, it is worth mentioning that estimating $\beta^{*}$ using the same simulated data used in getting the estimate introduces some dependence in the above estimate. However, it was shown in [16] that working with the estimated value of $\beta^{*}$ yield the same estimator’s performances as working with $\beta^{*}$ for large number of samples.

2.3 Performance Metrics

Many criteria have been used in practice in order to measure the efficiency of an unbiased estimator such as the asymptotically vanishing relative error, the bounded relative error, and the asymptotic optimality [2]. For any estimator $T^{{}^{\prime}}_{\gamma_{th}}\left(\mathbb{Y}\right)$ of $\alpha(\gamma_{th})$ with $\mathbb{Y}$ is distributed according to the PDF $g(\cdot)$, we have from the non-negativity of the variance of $T^{{}^{\prime}}_{\gamma_{th}}\left(\mathbb{Y}\right)$

[TABLE]

Using the fact that $\log\left(\alpha(\gamma_{th})\right)<0$, we get

[TABLE]

We say that the estimator $T^{{}^{\prime}}_{\gamma_{th}}\left(\mathbb{Y}\right)$ achieves the asymptotic optimality property if

[TABLE]

An equivalent definition of the asymptotic optimality criterion is: $\forall\epsilon>0$, we have

[TABLE]

Two interesting interpretations could be deduced from the asymptotic optimality property. First of all, when $\alpha(\gamma_{th})^{2}\rightarrow 0$ with an exponential rate, the second moment of $T^{{}^{\prime}}_{\gamma_{th}}\left(\mathbb{Y}\right)$ converges to zero with the same exponential rate. This is the best exponential rate that the second moment may converge with. Second, when the asymptotic optimality property holds, the number of simulation runs $M$ required to meet a fixed accuracy requirement satisfies $M=o\left(\alpha(\gamma_{th})^{-\epsilon}\right)$ for all $\epsilon>0$. Such a result ensures that an estimator with the asymptotic optimality criterion will certainly yield a substantial amount of variance reduction compared to naive MC simulation which requires a number of runs of the order of $\alpha(\gamma_{th})^{-1}$ to meet the same accuracy requirement.

A stronger criterion than the asymptotic optimality is the bounded relative error. In fact, this property holds when

[TABLE]

When the bounded relative error property holds, the number of simulation runs needed to meet a fixed accuracy requirement remains bounded regardless of how small $\alpha(\gamma_{th})$ is.

A further stronger criterion is the asymptotically vanishing relative error which holds when

[TABLE]

This means that the number of samples is getting smaller as we decrease the probability of interest while ensuring a fixed accuracy requirement.

3 Mean-Shifting Approach

In this section, we review the mean-shifting IS scheme proposed in [17]. The main idea is to consider an IS distribution resulting from shifting the mean of the multivariate normal distribution. More precisely, the IS PDF is chosen to be a multivariate normal with mean vector $\boldsymbol{\mu}+\mathbb{\Lambda}$ and covariance matrix $\mathbb{\Sigma}$ with $\mathbb{\Lambda}\in\mathbb{R}^{N}$:

[TABLE]

The likelihood ratio following this IS scheme is given as follows

[TABLE]

Hence, the probability of interest $\alpha(\gamma_{th})$ is re-written as

[TABLE]

By simple computation, it was shown in [17] that the second moment of $T_{\gamma_{th}}\left(\mathbb{Y}\right)$ is given by

[TABLE]

The remaining step is to determine the value of $\mathbb{\Lambda}$ that guarantees a variance reduction compared to naive MC estimator. Obviously, the optimal value is to minimize the second moment of the RV $T_{\gamma_{th}}\left(\mathbb{Y}\right)$ with respect to $\mathbb{\Lambda}$. However, the second moment in (21) is not known explicitly. The authors in [17] have then proposed to find $\mathbb{\Lambda}$ that minimizes an asymptotic equivalent given by replacing the probability in (21) by an equivalent expression given in [17]. The value of $\mathbb{\Lambda}$ proposed in [17] is as follows

[TABLE]

With this value of $\mathbb{\Lambda}$, it was proven in [17] that the second moment of the IS estimator satisfies:

[TABLE]

where $C$ is a constant that does not depend on $\gamma_{th}$. Thus, the asymptotic optimality property (14) holds. This result is based on the asymptotic behavior of $\alpha(\gamma_{th})$ as $\gamma_{th}\rightarrow 0$, see Theorem 1 of [17].

4 Improved Algorithm Using Control Variate

The results that will be shown in the present section represent the main contribution of our work. In fact, we aim here to combine the IS scheme presented in the previous section with a control variate technique to further achieve a variance reduction. We assume that the following assumption holds:

Assumption A: there exist an index $i\in\{1,2,\cdots,N\}$ such that $\sqrt{\Sigma_{ii}}<\rho_{ij}\sqrt{\Sigma_{jj}}$ for all $j\neq i$, where $\rho_{ij}=\frac{\Sigma_{ij}}{\sqrt{\Sigma_{ii}\Sigma_{jj}}}$ denotes the correlation coefficient between $Y_{i}$ and $Y_{j}$.

With no loss of generality, we may suppose that $i=1$. Under the above assumption, it was proven in [17] that the asymptotic behavior of the left-tail of the sum is dominated by only one component that corresponds to the index $i=1$ in assumption A, that is

[TABLE]

Moreover, under assumption A we have that $\bar{n}=1$ and the solution of the minimization problem (2) is given by

[TABLE]

Under assumption A, the value of $\mathbb{\Lambda}^{*}$ in (22) simplifies to

[TABLE]

Now, we present how the control variate technique could be combined with the previous IS scheme in order to achieve a further variance reduction. The control variable $Z_{\gamma_{th}}\left(\mathbb{Y}\right)$ is selected as follows:

[TABLE]

where the random vector $\mathbb{Y}$ is distributed according to $g(\cdot)$. The expected value $P(\gamma_{th})$ of the control variable $Z_{\gamma_{th}}\left(\mathbb{Y}\right)$ is given through straightforward computation by

[TABLE]

where $\Phi(\cdot)$ is the cumulative distribution function of the standard normal distribution. The main idea of the above choice of the control variable is that, with the values of $\mathbb{\Lambda^{*}}$ in (4) and under the PDF $g(\cdot)$, $X_{1}$ still represents the dominant component in the sense that it determines the asymptotic behavior of the left-tail of the sum. Hence, it is likely that each realizations that belongs to the set $\{X_{1}\leq\gamma_{th}\}$ will belong to the set $\{\sum_{i=1}^{N}{X_{i}}\leq\gamma_{th}\}$ for a sufficiently small threshold value $\gamma_{th}$.

The covariance between $T_{\gamma_{th}}\left(\mathbb{Y}\right)$ and $Z_{\gamma_{th}}\left(\mathbb{Y}\right)$ which is useful in the computation of the optimal value $\beta^{*}$ in (9) is given by

[TABLE]

The variance of $Z_{\gamma_{th}}\left(\mathbb{Y}\right)$ is given in a closed-form expression. In fact, via a simple computation, we have

[TABLE]

Now, using the value of $\Lambda_{1}^{*}$ in (4), we get

[TABLE]

The following lemma is very useful to study the performance of the estimator $T^{{}^{\prime}}_{\gamma_{th}}\left(\mathbb{Y}\right)$

Lemma 1**.**

There exists a constant $C_{1}$ such that

[TABLE]

with $P_{1}=\mathbb{E}_{g}\left[\mathbb{1}_{\left(\exp(Y_{1})\gamma_{th}\right)}\right]=\frac{1}{2}$ and $P_{2}=\mathbb{E}_{g}\left[\mathbb{1}_{\left(\sum_{i=1}^{N}{\exp\left(Y_{i}\right)}\leq\gamma_{th}\right)}\right]$.

Proof.

Via Cauchy Schwartz inequality, it follows that

[TABLE]

Now, let us compute $\mathbb{E}\left[L^{4}\mathbb{1}_{\left(\exp\left(Y_{1}\right)\leq\gamma_{th}\right)}\right]$. Using a similar computation as (4), we get

[TABLE]

Now, we use the asymptotic behavior of $\Phi(\cdot)$ [1]

[TABLE]

By combining (31), (4), and (4), we get

[TABLE]

From the expression of $\mathbb{\Lambda^{*}}$ in (4), we observe that $\mathbb{\Lambda^{*}}^{t}\mathbb{\Sigma}^{-1}\mathbb{\Lambda^{*}}=\frac{(\log(\gamma_{th})-\mu_{1})^{2}}{\Sigma_{11}}$. Hence, it follows that

[TABLE]

Thus, the proof is concluded. ∎

In the next lemma, we study the asymptotic behavior of $P_{1}-P_{2}$.

Lemma 2**.**

For $i\in\{2,\cdots,N\}$, let us define $a_{i}=\frac{\Sigma_{i1}}{\Sigma_{11}}-1$ and $c_{i}=\exp\left(\mu_{i}-\frac{\Sigma_{i1}}{\Sigma_{11}}\mu_{1}\right)$ with $a_{i}>0$ from Assumption A. Let $c_{i_{0}}\gamma_{th}^{a_{i_{0}}}=\max_{i}{c_{i}\gamma_{th}^{a_{i}}}$, then, for a sufficiently small $\gamma_{th}$, there a constant $C_{2}$ such that

[TABLE]

Proof.

Let us first re-write $P_{1}-P_{2}$ as follows:

[TABLE]

Using the value of $\mathbb{\Lambda^{*}}$ in (4), the above expression could be expressed as

[TABLE]

with $P(\cdot)$ is the probability measure under which $\mathbb{Y}$ in this case is a multivariate Gaussian vector with zero mean and covariance matrix $\mathbb{\Sigma}$. Let us denote by $\lambda$ the minimum eigenvalue of $\mathbb{\Sigma}^{-1}$, then the above expression could be upper-bounded by

[TABLE]

with $d_{1}=\frac{1}{\sqrt{|\mathbb{\Sigma}|\lambda^{N}}}$ and $\tilde{P}(\cdot)$ is the probability measure under which $\mathbb{Y}$ is now an independent Gaussian random vector with zero mean and covariance matrix $I_{N}/\lambda$ where $I_{N}$ denotes the identity matrix of order $N$. Note that for a sufficiently small threshold, the index $i_{0}$ is independent of $\gamma_{th}$. In the other hand, the probability in the right-hand side could be written as

[TABLE]

Let $\tilde{f}(\cdot)$ be the univariate normal PDF with mean zero and variance $1/\lambda$. Then, the quantity $I_{1}(\gamma_{th})$ is expressed as follows

[TABLE]

Let $Z=\sum_{i=2}^{N}{\exp(Y_{i})}$ and $\tilde{f}_{Z}(\cdot)$ its corresponding PDF, then we have

[TABLE]

Let $g_{\gamma_{th}}(z)=\Phi\left(\sqrt{\lambda}\log\left(1-c_{i_{0}}\gamma_{th}^{a_{i_{0}}}z\right)\right)$. Via a simple computation, we prove that for all $z\leq\frac{1}{2c_{i_{0}}\gamma_{th}^{a_{i_{0}}}}$

[TABLE]

Therefore, it follows that

[TABLE]

In the other hand, we have

[TABLE]

Therefore, by combining (46) and (4), we get

[TABLE]

In the other hand, the asymptotic behavior of the right-tail of the sum of Log-normal variates is given by [6]

[TABLE]

Hence, the proof is concluded. ∎

The results in Lemma 1 and Lemma 2 serve to study the goodness of the proposed estimator $T_{\gamma_{th}}^{{}^{\prime}}\left(\mathbb{Y}\right)$. The next proposition provides interesting results on the correlation coefficient as well as the optimal value $\beta^{*}$.

Proposition 1**.**

The correlation coefficient $\rho_{T_{\gamma_{th}}\left(\mathbb{Y}\right),Z_{\gamma_{th}}\left(\mathbb{Y}\right)}$ between $T_{\gamma_{th}}\left(\mathbb{Y}\right)$ and $Z_{\gamma_{th}}\left(\mathbb{Y}\right)$ goes to $1$ whereas the optimal value $\beta^{*}$ approaches $-1$ as $\gamma_{th}\rightarrow 0$.

Proof.

Let us first prove that $\beta^{*}$ approaches $-1$ as $\gamma_{th}$ goes to zero. In fact, from the expression of $\beta^{*}$ in (9), we have

[TABLE]

From the results of Lemma 1 and Lemma 2, we have that $\mathbb{E}_{g}\left[Z_{\gamma_{th}}^{2}\left(\mathbb{Y}\right)\right]\sim\mathbb{E}_{g}\left[T_{\gamma_{th}}^{2}\left(\mathbb{Y}\right)\right]$ as $\gamma_{th}\rightarrow 0$. Moreover, the results in (4) and (31) implies that $\alpha^{2}/\mathbb{E}_{g}\left[Z_{\gamma_{th}}^{2}\left(\mathbb{Y}\right)\right]\rightarrow 0$ as $\gamma_{th}$ goes to zero. In particular, this shows that $\mathrm{var}_{g}\left[Z_{\gamma_{th}}\left(\mathbb{Y}\right)\right]\sim\mathbb{E}_{g}\left[T_{\gamma_{th}}^{2}\left(\mathbb{Y}\right)\right]$ as $\gamma_{th}\rightarrow 0$. By combining these results and the fact that $\alpha(\gamma_{th})\sim P(\gamma_{th})$, we conclude that $\beta^{*}$ approaches $-1$ as $\gamma_{th}\rightarrow 0$.

Now, let us prove that $\rho_{T_{\gamma_{th}}\left(\mathbb{Y}\right),Z_{\gamma_{th}}\left(\mathbb{Y}\right)}$ approaches $1$ as $\gamma_{th}\rightarrow 0$. The expression of $\rho_{T_{\gamma_{th}}\left(\mathbb{Y}\right),Z_{\gamma_{th}}\left(\mathbb{Y}\right)}$ is given by

[TABLE]

Given that $\mathrm{var}_{g}\left[Z_{\gamma_{th}}\left(\mathbb{Y}\right)\right]\sim\mathrm{var}_{g}\left[T_{\gamma_{th}}\left(\mathbb{Y}\right)\right]$ as $\gamma_{th}\rightarrow 0$ (this follows from Lemma 1 and 2 and the fact that $\alpha^{2}/\mathbb{E}_{g}\left[Z_{\gamma_{th}}^{2}\left(\mathbb{Y}\right)\right]\rightarrow 0$ and $\alpha(\gamma_{th})\sim P(\gamma_{th})$ as $\gamma_{th}$ goes to [math]), it follows that

[TABLE]

Hence, following the same arguments as in the first part of the proof, we conclude that $\rho_{T_{\gamma_{th}}\left(\mathbb{Y}\right),Z_{\gamma_{th}}\left(\mathbb{Y}\right)}$ approaches $1$ as $\gamma_{th}$ goes to zero. ∎

The main conclusions that can be drawn from the above results are the following. Since $\rho$ approaches $1$, the improved estimator is guaranteed to achieve a variance reduction, with respect to the mean shifting IS scheme, that increases as we decrease the threshold value. The second interesting point is that we may select $\beta$ to be equal to $-1$ instead of working with the optimal unknown value $\beta^{*}$. In fact, in addition to retrieving the same estimator’s performances as with $\beta=\beta^{*}$, we avoid, when working with $\beta=-1$, the approximation of $\beta^{*}$ using simulated data which might causes some statistical issues especially when the number of samples is not large enough [16].

The next theorem exhibits the main result of the present work. In fact, we prove that the estimator $T_{\gamma_{th}}^{{}^{\prime}}\left(\mathbb{Y}\right)$ has an asymptotically vanishing relative error. Such a result represents a relevant contribution to the field of left-tail simulation of the sum of Log-normal variates since previous works have just achieved the weaker property; namely the asymptotic optimality criterion.

Theorem 1**.**

The estimator $T_{\gamma_{th}}^{{}^{\prime}}\left(\mathbb{Y}\right)$ has the asymptotically vanishing relative error property with $\beta=-1$, that is

[TABLE]

Proof.

By replacing $\beta=-1$ into (2.2), it follows that

[TABLE]

Hence, by dividing the above expression by $\alpha^{2}(\gamma_{th})$ and via the use of Lemma 1 and Lemma 2, we get

[TABLE]

Hence, the proof is concluded. ∎

5 Simulation Results

In this section, we provide some selected simulation results to validate some of the results proven in our work and also to investigate the amount of variance reduction achieved by employing the control variate technique compared to the IS scheme of [17]. The problem parameters are given as follows: the covariance matrix of the $4$-dimensional Gaussian vector is given by

[TABLE]

whereas the mean vector is $\mu=(4,4,4,4)^{t}$. Note that with the above choice of the covariance matrix, assumption A is satisfied and hence our intuition regarding the choice of the control variable in (27) is applicable. We define now some performance metrics that serve to compare the proposed estimator with that of [17]. We define the squared coefficient of variation of the estimator $T^{{}^{\prime}}_{\gamma_{th}}\left(\mathbb{Y}\right)$ as the ratio of its variance to its squared mean

[TABLE]

Similarly, we define the squared coefficient of variation of the estimator $T_{\gamma_{th}}\left(\mathbb{Y}\right)$ as

[TABLE]

The squared coefficient of variation, also referred as squared relative error in some references, is an interesting performance metric that serves to indicate the number of samples needed to achieve a fixed accuracy requirement. In fact, from the central limit theorem, the number of samples should be proportional to the squared coefficient of variation in order to maintain the width of the confidence interval constant.

We define also the amount of variance reduction between the mean-shifting IS approach of [17] and our proposed approach as

[TABLE]

In a first experiment, we aim to validate that the correlation between $T_{\gamma_{th}}\left(\mathbb{Y}\right)$ and $Z_{\gamma_{th}}\left(\mathbb{Y}\right)$ approaches $1$ as $\gamma_{th}$ goes to [math]. To this end, we plot in Fig. 1 the correlation coefficient $\rho_{T_{\gamma_{th}}\left(\mathbb{Y}\right),Z_{\gamma_{th}}\left(\mathbb{Y}\right)}$ as a function of $\gamma_{th}$. This figure clearly shows that the correlation coefficient $\rho_{T_{\gamma_{th}}\left(\mathbb{Y}\right),Z_{\gamma_{th}}\left(\mathbb{Y}\right)}$ approaches $1$ as the event of interest becomes rare and rare, i.e. as we decrease the threshold value $\gamma_{th}$. Such a result supports the performance of the proposed control variate technique in ensuring a variance reduction with respect to the mean shifting approach of [17] and interestingly this reduction is increasing as we decrease the value of $\gamma_{th}$.

In a second experiment, we want to investigate the value of $\beta^{*}$ as a function of the $\gamma_{th}$. We plot in Fig. 2 the approximate value of $\beta^{*}$ as a function of $\gamma_{th}$. From this figure, we easily observe that $\beta^{*}$ approaches $-1$ as $\gamma_{th}$ decreases which validate the result in Proposition 1. Thus, instead of working with the optimal unknown value $\beta^{*}$, this result suggests to work with a fixed value of $\beta$ equal to $-1$.

In Fig. 3, we plot the estimated value of $\alpha(\gamma_{th})$ as a function of $\gamma_{th}$ using the mean-shifting IS estimator of [17] as well as our proposed estimator with $\beta=\beta^{*}$ and $\beta=-1$.

Note that we employ the version of [16] where $\beta^{*}$ is estimated using the same simulated data as for the estimation. While this could result in a dependence across the replicants, it was noted in [16] that when $M$ is sufficiently large, the dependence could be ignored and we retrieve the same performance when using the approximate value of $\beta^{*}$ instead of the optimal unknown value. From Fig. 3, we observe that the three curves coincide perfectly and thus the three approaches yield accurate estimates of $\alpha(\gamma_{th})$.

Now, we focus on studying the efficiency of the three approaches. To this end, we plot in Fig. 4 the squared coefficient of variation given by the mean-shifting approach and the two versions of the control variate techniques; namely when $\beta$ is equal to the approximate value of $\beta^{*}$ and when $\beta$ is equal to $-1$. We point out from this figure that the squared coefficient of variation corresponding to the mean-shifting IS scheme is increasing as we decrease the threshold values. This is expected from our analysis since the mean-shifting approach is only asymptotically optimal. More precisely, from (4) and (4), the quantity $CV\left(T_{\gamma_{th}}\left(\mathbb{Y}\right)\right)$ is equivalent to a constant times $\log\left(\frac{1}{\gamma_{th}}\right)$. Hence, in order to meet a fixed accuracy requirement, i.e. maintain the width of the confidence interval fixed, the number of simulation runs required by the mean-shifting approach should be of the order of $\log\left(\frac{1}{\gamma_{th}}\right)$.

On the other hand, the squared coefficient of variation of the control variate estimator, using the two values of $\beta$, is decreasing as we decrease the threshold values. Such an observation is expected from the fact that our proposed estimator $T^{{}^{\prime}}_{\gamma_{th}}\left(\mathbb{Y}\right)$ with $\beta=-1$ has the asymptotically vanishing relative error property as it was proven in Theorem 1.

Hence, the number of samples needed by the proposed control variate technique with both values of $\beta$ is getting smaller in order to meet a certain accuracy measured by the width of the confidence interval. Moreover, we point out that asymptotically $CV(T^{{}^{\prime}}_{\gamma_{th}}\left(\mathbb{Y}\right))$ with $\beta=-1$ approaches the optimal value of $CV(T^{{}^{\prime}}_{\gamma_{th}}\left(\mathbb{Y}\right))$ which in accordance with the result proven in Proposition 1.

Finally, we aim to quantify the amount of variance reduction $\xi$ achieved by both versions of the proposed control variate estimators compared to the mean-shifting one. We plot then in Fig. 5 the value of $\xi$ as a function of $\gamma_{th}$. Interestingly, we observe that, for both versions of the control variate approaches, the amount of variance reduction is increasing as decrease the probability of interest. Moreover, the value of $\xi$ with $\beta=-1$ approaches, as we decrease $\gamma_{th}$, the optimal value corresponding to $\beta=\beta^{*}$ which is $\frac{1}{1-\rho_{T_{\gamma_{th}}\left(\mathbb{Y}\right),Z_{\gamma_{th}}\left(\mathbb{Y}\right)}^{2}}$.

6 Conclusion

In this paper, we considered the problem of estimating the probability that a sum of correlated Log-normal variates with Gaussian copula falls below a certain threshold. We developed a variance reduction technique based on a combination of a control variate approach with a previously developed mean-shifting importance sampling approach. Under a mild assumption on the covariance matrix of the corresponding multivariate Gaussian random vector, we proved that the proposed estimator has the asymptotically vanishing relative error property. This result is a major contribution in the field of left-tail estimation of sum of Log-normal random variables since previous estimators were only proven to satisfy weaker properties. We showed some selected simulations in order to validate our theoretical results and to quantify the performances of the proposed estimator with and without the employment of the control variate technique.