A New Family of Bivariate Exponential Distributions with Negative Dependence Based on Counter-Monotonic Shock Method

We introduce a new family of bivariate exponential distributions based on the counter-monotonic shock model. This family of distribution is easy to simulate and includes the Fréchet lower bound, which allows to span all degrees of negative dependence. The construction and distributional properties of the proposed bivariate distribution are presented along with an estimation of the parameters involved in our model based on the method of moments. A simulation study is carried out to evaluate the performance of the suggested estimators. An extension to the general model describing both negative and positive dependence is sketched in the last section of the paper.


Introduction
Exponential distributions are undoubtedly among the most popular and used distributions in many areas of application. They play a prominent role in a variety of fields, including reliability, hydrology, engineering, telecommunication, biological and environmental sciences, among others.
However, the exponential distribution cannot be naturally extended to the bivariate or the multivariate case in a unique way. As a result, the literature on bivariate exponential distributions is vast, including many different classes and models that have been developed in the past decades, for example, refs. [1][2][3][4][5][6][7][8][9][10], among others.
It is worth mentioning that most of the bivariate exponential models proposed in the literature are restricted to the case of non-negative dependence. Very few models have negative or both positive and negative correlation but do not fully complete the range of correlation [1 − π 2 6 , 1] (see Moran [11]) and necessitate a complex structure in their construction.
The main aim of this paper is to present a new bivariate exponential model that fully covers the negative dependence. To that end, we will adopt a technique based on the counter-monotonic shock model, which was introduced by Genet et al. [12]. This procedure is quite different from the common shock method used by Marshal and Olkin [8] to define a family of bivariate exponential distributions. Indeed, the latter is limited to model the positive dependence and imposes restrictions on the correlation structure, especially when the marginal are not identically distributed. In contrast, the counter-monotonic shock technique provides a flexible framework for building negatively correlated bivariate exponential distributions. Thanks to an appropriate parametrization, the resulting model can be viewed as a family of bivariate exponential distributions with given marginals and can be provided with a dependence parameter inducing the dependence in the model.
In addition, this family of distribution is easy to simulate and includes the Fréchet lower bound, which allows describing the full negative range of correlation, namely [1 − π 2 6 , 0]. We proceed as follows: We introduce our novel class of bivariate exponential distributions based on the counter-monotonic shock model in the next section. The derivation of the probability density function of this distribution is given in Section 3. We also present in this section the joint moment generating function, monotonicity, singularity, and scaling properties. Estimation of the model parameters through the moment method is discussed in Section 4. The proposed framework will be illustrated by simulations in Section 5. Concluding comments and directions for further research are presented in Section 6.

The Model
In order to define the suggested bivariate distribution, let us first recall the notion of a counter-monotonic random pair. Let K be a joint distribution with given marginal distributions F and G. The next double inequalities are due to [13]: As pointed out by Fréchet, these bounds are themselves bivariate distributions with the same marginals F and G. The counter-monotonic concept is related to the lower Fréchet bound, and it is defined as follows. Definition 1. The random pair (X, Y) with marginal distributions F and G, respectively, is counter-monotonic if its joint distribution function is the lower Fréchet bound. Equivalently, there exists a unit uniform random variable U such that Note that the counter-monotonic notion describes the perfect negative dependence.

The New Bivariate Exponential Distribution
In the following, we introduce a new family of bivariate exponential distributions with given marginals describing the negative dependence. The idea is based on the countermonotonic shock method introduced in [12]. The principle of this approach is to link independent exponential random variables through counter-monotonic ones in order to produce negative dependence. To this end, let λ i > 0, i = 1, 2 be the marginal parameters and let θ ∈ (0, 1) denote the dependence parameter. Definition 2. Let (X 1 , X 2 ) and (Y 1 , Y 2 ) be independent random pairs such that Y i ∼Exp(θλ i ) and X i ∼Exp(λ i (1 − θ)), i = 1, 2. Denote by G i the distribution functions of Y i , i = 1, 2, respectively. Suppose further that 1.
Y 1 and Y 2 are counter-monotonic, that is, X 1 , X 2 and U are independent.
The distribution of the random pair (X, Y), defined by is called the counter-monotonic shock bivariate exponential distribution.
Note that the set of all random pairs defined by (1) will be denoted BED − (θ, Λ), where Λ = (λ 1 , λ 2 ). Clearly, the latter is a family of bivariate exponential distributions with given marginals, since by construction, X∼Exp(λ 1 ) and Y∼Exp(λ 2 ). The parameter θ ∈ (0, 1) does not affect the marginal distributions; it can be interpreted as a dependence parameter. In fact, one observes that this family reaches the independence case when θ goes to 0 and it approaches the perfect negative dependence described by the Fréchet lower bound when θ goes to 1, respectively.
Thanks to the relations Y 1 = G −1 1 (U) and Y 2 = G −1 2 (1 − U), one deduces an interesting alternative representation of (1) given by where Z is an exponential random variable, with parameter 1, which is independent of X 1 and X 2 . This presentation provides an easy way to simulate data from this model through the next steps: 1.
As illustrated in Figure 1, the proposed distribution seems to have both absolutely continuous and singular components. This interesting property will be examined in Proposition 2.

Properties of the New Bivariate Exponential Distribution
The following section will be consecrated to investigating the properties of the new family of bivariate exponential distribution based on the counter-monotonic shocks. We will start with the joint survival function of the distribution and then derive the corresponding joint probability density function. It will then be followed by an analysis of the product moment of the distribution, the coefficient of correlation, and the moment generating function. Monotonicity and scaling properties will also be discussed.

Survival and Density Functions
We now study the joint survival function associated with the new family BED − (θ, Λ) and then deduce the corresponding joint probability density function. Proposition 1. The survival function of (X, Y)∼BED − (θ, Λ) is given bȳ where x + = max(x, 0) for any x ∈ R.
Proof. Using the fact that the random variables X 1 , X 2 and U are independent, one has from (1), for all (x, y) ∈ R 2 + , which completes the proof of the proposition.
Proof. Using the relation between the K θ andK θ the result can be immediately deduced.
Let us now recall the beta function and the incomplete beta function defined, respectively, by These functions are linked to the beta distribution. In fact, a random variable Z follows a beta distribution with parameters a > 0 and b > 0 if its distribution is defined, for It could be readily seen that By virtue of the identity, The singular component of the survival function is then given by where the normalized constant α is obtained by tending (x, y) to (0, 0) in the previous formula. This leads to α = B 1 θ , 1 θ . Similarly, the continuous part of the survival function is given bȳ Putting (6)-(8) together, we have the desired decomposition. This ends the proof of Proposition 2.
The density function of the proposed family of distribution can be obtained immediately from the previous result. and Proof. It is easily seen that f 1 is the density function corresponding to the continuous part of the survival functionK s,θ described by (9). Next, from Proposition 2, the singular component of the survival function of (X, Y) is given by This ensures that the survival function of X in the singular part isF 0 ( This ends the proof of Corollary 2.
One can check that which guarantees that f θ is a density function. In addition, one observes that the probability that (X, Y) lies in the singular part represented by the curve

Monotonicity
In the following, we show that the family BED − (θ, Λ) is ordered in terms of θ in the negative quadrant dependence ordering. This means that the parameter θ can be considered as a dependence parameter for the family BED − (θ, Λ).
The above result shows that the strength of the dependence of the pair (X, Y) in BED − (θ, Λ) decreases with θ ∈ [0, 1]. Hence, the covariance as well as the correlation of random pairs in BED − (θ, Λ) decrease with respect to θ. In addition, one sees that, for all (x, y) ∈ R 2 ,K θ (x, y) ≤K 0 (x, y) =F(x)Ḡ(x). The latter outlines that the components of (X, Y) ∈ BED − (θ, Λ) are negatively quadrant dependent. In particular, the correlation of any (X, Y) ∈ BED − (θ, Λ) is negative.

Product Moment and Correlation Structure
Recall the partial derivatives of the Beta function for x > 0 and y > 0

Proposition 4.
For θ ∈ (0, 1), the product moment of (X, Y)∼BED − (θ, Λ) is given, for any positive integers i and j, by Define a(x) = −θ −1 ln(1 − e −θx ). Hence, where the second line follows from the fact that Similarly, one has Finally, we have This completes the proof of Proposition 4.

Corollary 3.
The correlation coefficient of (X, Y)∼BED − (θ, Λ) is given by Proof. Using Equation (10), one has for i = j = 1, To complete the proof of the corollary, we will make use of the following properties of the beta function y) and B(x, y + 1) = y x + y B(x, y).
It results that Substituting (12) in (11), we get the following expression for the covariance of (X, Y), The result follows by using the fact that corr(X, Y) = λ 1 λ 2 cov(X, Y).
Let Γ(F, G) be the space of all bivariate random variables with given exponential marginal distributions F and G with parameters λ 1 and λ 2 , respectively. The minimal correlation in the space Γ(F, G) is calculated by using the lower Fréchet bound. This means that, for all (X, Y) ∈ Γ(F, G), one has, where U is a random variable uniformly distributed over [0, 1]. The above minimal correlation is given by Hence, the full correlation range of negative dependence in the space Γ(F, G) is [1 − π 2 /6, 0]. Note the proposed family of distribution BED − (θ, Λ) describes this full negative range of correlation since it includes the Fréchet lower bound. In particular, one has

Scaling Property
Analogously to the univariate case, the proposed family of distributions enjoys the scaling property.
Proof. The proof is straightforward, and therefore, it is omitted.

Moment Generating Function
Here, we derive an explicit expression for the moment generating function of the family (X, Y)∼BED − (θ, Λ) in terms of the beta function. Proposition 6. The moment generating function of (X, Proof. For the sake of easy reference, we restate b 1 (x) and b 2 (x) introduced earlier in (5) where , , .

Proposition 7.
For all θ ∈ (0, 1), one has Proof. In fact, it is well known (see, e.g., Theorem 8 on p. 52 of [14]) that S 12 is a consistent and asymptotic Gaussian estimator of population covariance h λ 1 ,λ 2 (θ), namely The result is then derived from the Delta method and Slutsky's lemma applied toθ = h −1 λ 1 ,λ 2 (S 12 ). The previous proposition will be useful to establish a confidence interval for θ. To this end, let us first compute The mixed moments E(X i Y j ), i, j = 1, 2 involved in the previous formula can be obtained from Proposition 4 as follows Using the properties of the partial derivatives of the beta function (see Appendix A) and putting (11), (13), (19)-(22) and together, it follows that where with ψ (x) is the derivative of the Digamma function ψ(x). From Proposition 7, one can build an asymptotic confidence interval for the dependence parameter θ. Thus, the (1 − α) × 100% confidence interval for θ is given by the following formula:θ ± z α/2 h (θ) −1 σ(θ,λ 1 ,λ 2 ) √ n .

Simulation Study
In this section, we will illustrate the performance ofθ, the estimator of the dependence parameter θ. We will be examining the finite-sample accuracy of our estimates for different sample sizes. An asymptotic confidence interval for θ will also be provided.
To assess the performance of the moment-based estimators, the marginal parameters were held fixed at Λ = (2, 4), and various values of θ were chosen to cover a broader range of negative dependence. Different sample sizes, n, are considered, and each scenario was replicated 500 times. As estimation of the parameters λ 1 and λ 2 is standard, we focus on the results of estimating θ.
Recall that the method of moments, used to estimate the model parameters λ 1 and λ 2 , requires estimation of the dependence parameter θ ∈ (0, 1) by solving for the unique root of Table 1 exhibits the estimateθ of the dependence parameter θ, bias, mean squared error (MSE), and confidence interval estimations for θ. It illustrates that simulation results obtained by the method of moments are consistent. In fact,θ provides a good estimator for the dependence parameter θ, bias, and MSE ofθ decrease as the sample size increases. As expected, the confidence intervals get narrower as the sample size grows. Moreover, we observed from many more simulations, not presented here, that the estimatorθ performs very well regardless of the choice of λ 1 and λ 2 .

Conclusions
The main purpose of this paper was to introduce a new class of bivariate exponential distribution that fully covers the negative dependence. The concept of counter-monotonic shock was used to create the negative dependence among the exponential components of the model. The basic features of this class were studied, and moment-based estimators of model parameters were derived. This family of distributions could be easily interpreted and simulated. Moreover, the proposed method can be adapted to derive a general model describing both negative and positive dependence. The Fréchet family of distribution can be a good candidate to construct such a model. In this direction, provided exponential marginal distributions F and G, this family of distributions is defined, for any θ 0 ∈ [0, 1], by Let us consider independent random variables U, Y 1 , Y 2 , and Z such that U∼U [0,1] , Y i ∼Exp(λ i (1 − θ)), i = 1, 2 and Z∼Bernoulli(θ 0 ) with (θ 0 , θ) ∈ [0, 1] 2 . The key idea allowing to build such a model is based on the fact that Indeed, one can construct a bivariate exponential random variable by It can be verified that the marginal distributions of (X, Y) are fixed, namely, X∼Exp(λ 1 ) and Y∼Exp(λ 2 ). This general model will be explored in a forthcoming paper.