On a Bivariate Copula for Modeling Negative Dependence: Application to New York Air Quality Data

In many practical scenarios, including finance, environmental sciences, system reliability, etc., it is often of interest to study the various notion of negative dependence among the observed variables. A new bivariate copula is proposed for modeling negative dependence between two random variables that complies with most of the popular notions of negative dependence reported in the literature. Specifically, the Spearman's rho and the Kendall's tau for the proposed copula have a simple one-parameter form with negative values in the full range. Some important ordering properties comparing the strength of negative dependence with respect to the parameter involved are considered. Simple examples of the corresponding bivariate distributions with popular marginals are presented. Application of the proposed copula is illustrated using a real data set on air quality in the New York City, USA.


Introduction
Copulas provide an effective tool for modeling dependence in various multivariate phenomena in the fields of reliability engineering, life sciences, environmental science, economics and finance, etc (Fontaine et al., 2020;Cooray, 2019;Joe, 2015, Ch-7). Specifically, in recent decades, bivariate copulas were used to generate bivariate distributions with suitable dependence properties (Lai and Xie, 2000;Bairamov and Kotz, 2003;Finkelstein, 2003;Durante et al., 2012;Mohtashami-Borzadaran et al., 2019). The detailed discussion of historical developments, obtained results and perspectives along with the up to date theory can be found in Durante and Sempi (2015) and Hofert et al. (2018). It should be noted that most copulas available in the literature possess some limitations in modeling negatively dependent data, which is a certain disadvantage, as negative dependence between vital variables is often encountered in real life. Lehmann (1966) introduced several concepts of negative dependence for bivariate distributions.
Later, Esary and Lehmann (1972) and Yanagimoto (1972) extended the corresponding definitions and developed stronger notions of bivariate negative dependence. See Balakrishnan and Lai (2009) for detailed discussion on popular dependence notions and their applications in the context of continuous bivariate distributions. Scarsini and Shaked (1996) provided a detailed overview of the corresponding ordering properties for the multivariate distributions. These results provide useful tools for describing the dependence properties of copulas with respect to a dependence parameter.
However, only a few bivariate copulas that allow for a simple and meaningful analysis of this kind have been developed and studied in the literature so far. The Farlie- Gumbel-Morgenstern (FGM) family of distributions exhibits negative dependence, but the Spearman's rho for this family lies within the interval [−1/3, 1/3] (Schucany et al., 1978). Bairamov and Kotz (2000) and Bekrizadeh et al. (2012) have considered the four-parameter and the three-parameter extensions of the FGM family proposed by Sarmanov (1996), with Spearman's rho lying within the interval [−0.48, 0.50], and [−0.5, 0.43], respectively. To address this issue Amblard and Girard (2009) proposed another extension, but its application is limited because of a singular component that is concentrated on the corresponding diagonal. Some other extensions of the FGM copula are discussed in Ahn (2015) and Bekrizadeh and Jamshidi (2017). Hürlimann (2015) have proposed a comprehensive extension of the FGM copula with the Spearman's rho and Kendall's tau attaining any value in (−1, 1).
However, the dependence properties and ordering properties of these copulas are not well studied in the literature. Recently, Cooray (2019) proposed a new extension of FGM family which exhibits negative dependence among the variables in a very strong sense. However, its Spearman's rho and Kendll's tau are restricted to [−0.70, 0], and [−0.52, 0], respectively. Thus, it is quite a challenging problem to construct a flexible bivariate copula with the correlation coefficient that takes any value in the interval (−1, 0). Moreover, it is not sufficient just to suggest this type of copula, but it is essential to describe its properties (including relevant stochastic comparisons) especially in the case of strong notions of dependence. In many real life scenarios, paired observations of non-negative variables possess strong negative dependence. For example, rainfall intensity and duration are jointly modeled incorporating their negative dependence for the study of derived flood frequency distribution (Kurothe et al., 1997). This paper is motivated by a real case study on air quality for New York Metropolitan area where the joint distribution of the wind speed and ozone level exhibits strong negative dependence (See Section 7). We believe that the current study meets to some extent this challenge, as we propose an absolutely continuous negatively dependent copula that satisfies most of the popular notions of negative dependence available in the literature with correlation coefficients in the interval (−1, 0). The paper is organized as follows. In Section 2, we describe the baseline (for the proposed copula) distribution and discuss some basic properties including conditional distributions and correlation coefficients. Various notions of negative dependence in the context of the proposed copula and ordering properties are considered in Section 3, and Section 4, respectively. Section 5, provides some examples of negatively dependent standard bivariate distributions. The estimation methodologies are discussed in Section 6. In Section 7, as an illustration, we provide a real case study. Finally, some concluding remarks are given in Section 8. Bhuyan et al. (2020) proposed a negatively dependent bivariate life distribution that possesses nice closed-form expressions for the joint distributions and exhibits various strong notions of negative dependence reported in the literature. Most importantly, the correlation coefficient may take any value in the interval (−1, 0). One of the marginal distribution is Exponential and the other belongs to skew log Laplace family (Dixit and Khandeparkar, 2017). We utilize the negative dependence structure inherent in this model and formulate a copula with strong negative dependence. The joint distribution function and the marginal distributions are given by
The survival copula, is the functionC which couples the joint survival function to its marginal survival functions. It is easy to show thatC is a copula, and is related to the copula C via the Nelsen (2006, p-32) for details. The survival copula and the density function of the proposed copula C θ (u, v) are given bȳ   and respectively.

Conditional Copulas
The conditional copula of U given V = v, is as follows. For 0 < v ≤ θ (1+θ) , whereas for θ (1+θ) < v < 1, The conditional mean and variance of U | V = v are given by respectively.
Remark 2.1. The regression of U on V = v is linearly decreasing in v for 0 < v ≤ θ θ+1 , and independent of v for θ θ+1 < v < 1. Also, it is interesting to note that the conditional variance of U | V = v is an increasing function of v and bounded from above by θ+1 (θ+2) 2 (θ+3) .
The conditional copula of V given U = u, is given by The conditional mean and variance of V | U = u, are given by Remark 2.2. The regression of V on U = u is strictly decreasing in u.
One can use the conditional copula of U given V = v, provided in (5) and (6), to simulate from the proposed copula C θ , given by (3), using the following steps.
Step I. Simulate v i and u * i independently from standard uniform distribution.
Step II.
Step III. Repeat Step I and Step II n times to obtain independently and identically distributed real- A similar algorithm can be elaborated to simulate from C θ based on the conditional copula of V given U , provided in (7).The associated R programme for the aforementioned algorithm are provided in the Supplementary material. The Scatter plots based on 500 simulated observations using the aforementioned algorithm for four different values of θ are given in Figure 3. As expected, the data points are getting closer to the diagonal v = −u for higher values of θ.

Basic Properties
In this Subsection, we present three important propositions related to the proposed copula. The detailed proofs are presented in Appendix A.

Measures of Dependence
Measures of dependence are commonly used to summarize the complicated dependence structure of bivariate distributions. See Joe (1997, Ch-2), Nelsen (2006, Ch-5) and Hofert et al. (2018, Ch-2) for a detailed review on measures of dependence and its associated properties. In this section, we derive the expressions of the Kendall's tau and the Spearman's rho for the proposed copula C θ . Essentially, these coefficients measure the correlation between the ranks rather than actual values of X and Y .
Therefore, these coefficients are unaffected by any monotonically increasing transformation of X and Y .
Definition 2.1. Let X and Y be the continuous random variables with the dependence structure described by the copula C. Then the population version of the Spearman's rho for X and Y is given Proposition 2.4. Let (X, Y ) be a random pair with copula C θ . The Spearman's rho is given by which is a decreasing function in θ and takes any values between -1 and 0. Proposition 2.5. Let (X, Y ) be a random pair with copula C θ . Then the Kendall's tau is given by which is a decreasing function in θ and takes any values between -1 and 0.
In Figure 4, we have plotted the Spearman's rho and the Kendall's tau against the dependence parameter θ. It is easy to see that the Spearman's rho is less than the Kendall's tau for all θ > 0.

Connections with notions of Negative Dependence
As discussed in Subsection 2.3, the Spearman's rho and the Kendall's tau measure the correlation between two random variables. However, it is possible that these random variables may have the strong correlation, but possess the weak association with respect to different notions of dependence or vice versa. In this section, we discuss several relevant notions of negative dependence, namely Quadrant Dependence, Regression Dependence and Likelihood Ratio Dependence, etc., and explore, whether the corresponding properties are satisfied by the proposed copula or not. First, we provide the definitions of the aforementioned dependence notions as discussed in Nelsen (2006) and Balakrishnan and Lai (2009).
Definition 3.1. Let X and Y be continuous random variables with copula C. Then for all (x, y) ∈ R 2 , where R 2 is the domain of joint distribution of X and Y , or equivalently a copula C is said to be NQD if for all (u, v) ∈ I 2 , C(u, v) ≤ uv.

Y is left tail increasing in
is a nonincreasing function of y for all x.
6. Y is stochastically decreasing in X denoted as SD(Y | X), (also known as negatively regression is a nonincreasing function of x for all y.
7. X is stochastically decreasing in Y denoted as SD(X | Y ), (also known as negatively regression is a nonincreasing function of y for all x. 8. Let X and Y be continuous random variables with joint density function h(x, y). Then X and Y are negatively likelihood ratio dependent, denote by NLR(X,Y), if h(x 1 , y 1 )h(x 2 , y 2 ) ≤ h(x 1 , y 2 )h(x 2 , y 1 ) for all x 1 , x 2 , y 1 , y 2 ∈ I such that x 1 ≤ x 2 and y 1 ≤ y 2 . Now, in the following theorems, we establish that the proposed copula C θ satisfies all the aforementioned dependence properties. The detailed proofs are provided in Appendix B.
Theorem 3.1. Let X and Y be two random variables with copula C θ . Then (i) X and Y are LTI(Y | X), (ii) X and Y are LTI(X | Y ), (iii) X and Y are RTD(Y | X), and (iv) X and Y are RTD(X | Y ).
Theorem 3.2. Let X and Y be two random variables with copula C θ . Then (i) X and Y are SD(Y | X), and (ii) X and Y are SD(X | Y ).
Theorem 3.3. Let X and Y be two random variables with copula C θ . Then X and Y are NLR.
Remark 3.1. Two random variables X and Y with copula C θ are NQD. This directly follows from Theorem 3.3. See the interrelationships between different concepts of negative dependence summarised in (Balakrishnan and Lai, 2009, p-130) for details.

Ordering Properties
In Section 3, several negative dependence properties of the proposed copula C θ has been investigated for the fixed θ > 0. In this section, we discuss the ordering properties of the proposed copula C θ , which provides a precise (and also intuitively expected) notion for one bivariate distribution being more positively or negatively associated than another. For this purpose, we first recall the definitions of the dependence orderings for bivariate distributions. These definitions describe the strength of dependence of a copula with respect to its dependence parameter θ. Lehmann (1966) was first to introduce the NQD and NRD notions. Following this notions, Yanagimoto and Okamoto (1969) introduced the ordering properties as defined below.
Definition 4.1. Let F and G be two bivariate distributions with the same marginals. Then F is said to be smaller than G in the NQD sense denoted as F ≺ N QD G if F (x, y) ≥ G(x, y) ∀x and y.
Definition 4.2. Let F and G be two bivariate distributions with the same marginals, and let (U, V ) and (X, Y ) be two random vectors having the distributions F and G, respectively. Then F is said to be smaller than G in the NRD sense, denoted by F ≺ N RD G or (U, V ) ≺ N RD (X, Y ) if, for any for any u, v ∈ I, where F V |U denote the conditional distribution of V given U = u and F −1 V |U denote its right-continuous inverse. Equivalently, F ≺ N RD G if and only if G −1 Y |X F V |U (y | x) | x is decreasing in x for all y (Fang and Joe, 1992).
Later, Kimeldorf and Sampson (1987) have introduced and studied in detail the notion of the Negatively Likelihood Ratio dependence ordering that is described in the following definition. Let the random variables X and Y have the joint distribution G(x, y). For any two intervals I 1 and I 2 of the real line, let us denote I 1 ≤ I 2 if x 1 ∈ I 1 and x 2 ∈ I 2 imply that x 1 ≤ x 2 . For any two intervals I and J of the real line let G(I, J) represent the probability assigned by G to the rectangle I × J. and (X, Y ) be two random vectors having the distributions F and G, respectively. Then F is said to be smaller than G in the NLR dependence sense, denoted by if F (I 1 , J 1 )F (I 2 , J 2 )G(I 1 , J 2 )G(I 2 , J 1 ) ≥ F (I 1 , J 2 )F (I 2 , J 1 )G(I 1 , J 1 )G(I 2 , J 1 ) whenever I 1 ≤ I 2 and J 1 ≤ J 2 . When the densities F and G exist and denoted by f and g, respectively, then the aforementioned condition equivalently is written as f (x 1 , y 1 )f (x 2 , y 2 )g(x 1 , y 2 )g(x 2 , y 1 ) ≥ f (x 1 , y 2 )f (x 2 , y 1 )g(x 1 , y 1 )g(x 2 , y 1 ) whenever x 1 ≤ x 2 and y 1 ≤ y 2 .
In the following theorems, we derive the sufficient conditions under which one bivariate distribution will be more negatively associated than another. The detailed proofs of the following theorems are presented in Appendix C.

Examples
Traditionally, bivariate life distributions available in the literature are positively correlated (Balakrishnan and Lai, 2009). However, in many real life scenarios, paired observations of non-negative variables are negatively correlated (Bhuyan et al., 2020). For example, the rainfall intensity and the duration are jointly modeled incorporating their negative dependence for the study of the corresponding flood frequency distribution (Kurothe et al., 1997). Gumbel (1960) and Freund (1961) have proposed the bivariate Exponential distributions with lower bound of the correlation coefficient as −0.4. In this section, several specific families of bivariate distributions are generated using the proposed copula (3) with different choices for marginal distribution. For modelling purposes, the Lognormal, Weibull, and Gamma distributions are popular among practitioners in the fields of engineering, medical science, and environmental science (Sharma et al., 2016;Pobočíková et al., 2017;Ramos et al., 2019). We consider these choices as baseline distribution. We first define a bivariate Weibull and bivariate Gamma distribution. Then we consider a case when the marginals are different, one from the Lognormal and another from the Weibull family. It should be noted that the resulting bivariate distributions can be described implementing all notions of negative dependence discussed in Section 3 and 4.

Estimation Methodology
In a classical parametric setting, a straightforward approach is to estimate the dependence parameter and the parameters associated with the marginals using maximum likelihood method. This method is theoretically valid but there are some practical limitations. Firstly, the estimation of the dependence parameter θ depends on the parametric assumptions made on the marginals and the estimate of θ will be biased if the marginals are misspecified. The second drawback is computational as the log-likelihood function involves potentially large number of parameters and high-dimensional optimization is known to be challenging. See Hofert et al. (2018, Ch-4) for details. To avoid aforementioned computational burden Joe (1997) proposed a two stage method known as inference function for margins (IFM). This estimation method is based on two separate maximum likelihood estimations of the univariate marginal distributions, followed by an optimization of the bivariate likelihood as a function of the dependence parameter. Similar to maximum likelihood estimate, the estimate of θ based on IFM may be biased if the margins are partially misspecified (Hofert et al., 2018, p-136). Although the IFM has computational edge, it is less efficient compared to the maximum likelihood estimate (Joe, 2015, Ch-5).
We propose to use a method that close in spirit to the method of inference function for margins (IFM) but avoids the issue with misspecified marginals for the estimation of θ. In contrast to IFM, we do not maximize the bivariate likelihood. Instead, we determine the dependence parameter using method of moments (Hofert et al., 2018, p-141). The method of fitting a bivariate distribution with marginals F η i (·), indexed by parameter η i for i = 1, 2, involves the following steps: (i) Obtain the estimatesη i for i = 1, 2 using maximum likelihood method.
These steps are easy to execute and familiar to the practitioners of different fields of science. This method allows the copula to adequately approximate the dependence structure of the bivariate data, which is of prime concern from a practical point of view.

Exploratory Data Analysis
For an illustrative data analysis based on the proposed copula, we consider a data set on daily air quality measurements for 153 days in the New York Metropolitan Area from May 1, 1973, to September 30, 1973. Information on average wind speed (in miles per hour) and mean ozone level (in parts per billion), were obtained from the New York State Department of Conservation and the National Weather Service, USA. This data set is openly available in R software. See Chambers et al. (1983, Ch 2-5) for the detailed description of the data. Ozone in the upper atmosphere protects the earth from the sun's harmful rays. On the contrary, exposure to ozone also can be hazardous to both humans and some plants in the lower atmosphere. Variations in weather conditions play an important role in determining ozone levels (Khiem et al., 2010;Topcu et al., 2003). In general, concentration of the ozone level is affected by a wind speed. High winds tend to disperse pollutants, which in turn, dilute the concentration of the ozone level. However, stagnant conditions or light winds allow pollution levels to build up and thereby, the ozone level too becomes larger. Environmental scientists and meteorologists are interested in the study of the effect of a wind speed on the distribution patterns of ozone (Gorai et al., 2015) levels. For our analysis, we consider 116 observations discarding the missing values and presented the scatter plot of average wind speed versus ozone levels in Figure 5

Modeling Wind Speed and Ozone Level
In the field of engineering and environmental science, Lognormal, Weibull, and Gamma distributions are widely used for modeling wind speed recorded in the same location (Shepherd, 1978;Monjean and Robyns, 2015;Pobočíková et al., 2017;Dhiman et al., 2020;Ramadan et al., 2020). These distributions are also used for modeling the level of various pollutants and ozone level (Sharma et al., 2016;Souza et al., 2018;Mishra et al., 2021). Therefore, we consider these three models for estimation of the parameters associated with the marginal distributions of the wind speed and the mean ozone level. Based on the Akaike information criterion, the Gamma distribution fits both marginals better as compared with other choices. The maximum likelihood estimates of the shape and the scale parameters are obtained as 7.171 and 1.375, respectively, for the wind speed, and the same for the mean ozone levels are 1.7 and 24.775, respectively. The estimate of the dependence parameter is obtained asθ = 0.765. Therefore, the joint distribution of wind speed and the mean ozone level is represented by the bivariate Gamma distribution provided in Example 5.2, and presented graphically in Figure 5(b). Following Balakrishnan and Ristic (2016), we then use bootstrap based Kolmogorov-Smirnov test to check whether the Gamma distribution is a good fit for the marginals. Also, we evaluate the goodness of fit of the proposed copula based on Kolmogorov-Smirnov statistic utilising the bootstrap algorithm proposed by Genest et al. (2006). We find the proposed model fits the data reasonably well. The R programme related to the proposed estimation methodology are provided in the Supplementary material. In Figure 5(c) we present the contour plot of the the distribution of wind speed and mean ozone level. It indicates that the concentration of mean ozone level varies from 6-30 ppb when the wind speed is within 7-16 mph. The estimated conditional distributions of the mean ozone level keeping the wind speed fixed at the empirical first decile (5.7 mph), median (9.7 mph), and ninth decile (14.9 mph) are presented in Figure 5(d). It is easy to see that the distribution of the mean ozone level decreases stochastically (in the sense of the usual stochastic order) as the wind speed increases. This visual representation of the regression dependence property indicates that the ozone level distributions below the level of 60 ppb differ significantly with wind speed. This can assist in formulating policies and guidelines to choose between locations to avoid health hazards related to high ozone levels.

Concluding Remarks
We construct the new flexible bivariate copula for modeling negative dependence between two random variables. Its correlation coefficient takes any value in the interval (−1, 0), which was not the case for other copulas reported in the literature. It is important to note that the Spearman's rho and the Kendall's tau have a simple one-parameter form with negative values in the full range. The properties of the proposed copula is an agreement with most of the popular notions of negative dependence available in the literature, namely quadrant Dependence, regression dependence and likelihood ratio dependence, etc. It is an interesting problem to consider a semi-parametric generalisation of the proposed copula and investigate its associated properties. Another possible direction of future research could be a multivariate extension of the proposed copula using the approaches considered by Fischer and Köck (2012) and (Mazo et al., 2015).
For an illustrative data analysis based on the proposed copula, we consider a data set on daily air quality measurements for New York Metropolitan Area. Based on the observed data, we find that wind speed and ozone levels strongly dependent in the NQD sense. We consider three different models (Lognormal, Weibull, and Gamma distributions) for estimation of parameters associated with the marginal distributions of the wind speed and the mean ozone level. It is shown that the Gamma distributions fits better for both marginals and that the distribution of the mean ozone level decreases stochastically (in the sense of the usual stochastic order) as the wind speed increases. The  scope of the proposed copula goes far beyond this particular application. For example, biomedical researchers can utilize the proposed copula in studying the negative association between BMI and glycated proteins (Espasandín-Domínguez et al., 2019). One can also extend the proposed copula to asses and model the nonlinear and asymmetric negative dependence over time in security and commodity markets (Liu et al., 2017).
Appendix A A1. Proof of Proposition 2.1.

A3. Proof of Proposition 2.3.
To establish the absolute continuity of the proposed copula C θ , it is required to show for every (u, v) ∈ I 2 .
Case I.
Case II. For 0 < u < 1, and θ 1+θ < v < 1, we have Therefore, the results follows by combining Case I and II.
(ii) In view of Theorem 5.2.5 in (Nelsen, 2006, p-192), the necessary and sufficient condition for LTI(X | Y ) is that, C(u,v) v is nondecreasing in v, for any u in I.
(iii) To establish RTD(Y | X), it is sufficient to show that v−C(u,v) (1−u) is a nondecreasing function in u for any v ∈ I (Nelsen, 2006, Theorem 5.2.5, p-192).
To establish SD(Y | X) property of the proposed copula C θ , we utilise the geometric interpretation of the stochastic monotonicity given in Corollary 5.2.11 of (Nelsen, 2006, p-197). Therefore, it is sufficient to show that C θ (u, v) is a convex function of u. Similarly, SD(X | Y ) can be established by showing C θ (u, v) is a convex function of v.
Hence C θ (u, v) is a convex function of u.

B3. Proof of Theorem 3.3.
To established the NLR between X and Y with copula C θ , we need to show c θ (u 1 , v 1 )c θ (u 2 , v 2 ) ≤ c θ (u 1 , v 2 )c θ (u 2 , v 1 ) holds for all u 1 ≤ u 2 , and v 1 ≤ v 2 , where c θ (u, v) is the copula density given in (4). Note that for the proposed copula C θ , the aforementioned condition holds with equality for all u 1 ≤ u 2 and v 1 ≤ v 2 in I.
The results directly follow from Proposition 2.1.

C3. Proof of Theorem 4.3.
Let θ 1 ≤ θ 2 . Now, it is easy to verify that the condition provided in Definition 4.3 holds for any choice of u 1 , u 2 , v 1 , v 2 , where u 1 ≤ u 2 , v 1 ≤ v 2 .