Estimating the parameters of a dependent model and applying it to environmental data set

In this paper, a new dependent model is introduced. The model is motivated using the structure of series-parallel systems consisting of two series-parallel systems with a random number of parallel sub-systems that have fixed components connected in series. The dependence properties of the proposed model are studied. Two estimation methods, namely the moment method, and the maximum likelihood method are applied to estimate the parameters of the distributions of the components based on observing the system's lifetime data. A Monte Carlo simulation study is used to evaluate the performance of the estimators. Two real data sets are used to illustrate the proposed method. The results are useful for researchers and practitioners interested in analyzing bivariate data related to extreme events.


Introduction and some preliminaries
In these decades, modeling dependence has been discussed in many articles by copulaapproaches. Copulas are useful tools for characterizing dependence between dependent random variables. In statistical concepts, a copula is a function that couple's joint distribution function to its marginal distributions. Let us recall some preliminaries of the copula theory that will be used in the following sections. Let (X, Y) be a continuous random vector with bivariate distribution function F(., .) and marginal distributions F 1 (.) and F 2 (.), respectively. Applying Sklar's Theorem (Sklar,[33]), the joint distribution function and the joint survival function of (X, Y) are given by where D(., .) andD(., .) are the unique copula and survival copula of (X, Y), respectively. In real life, extreme events happen and most of them are dependent phenomenon which their corresponding dependent models are related to extreme-value copulas. We remind that the copula D(., .) is called extreme value if D(u, v) = D r (u 1/r , v 1/r ) for all r > 0 and u, v ∈ (0, 1). Gülpinar and Katata [13] used the extreme-value copulas for modeling of oil and gas supply disruption risks. We refer to Pikhands [30], Kotz and Nadarajah [18] and Goudendorf and Segers [12] for more details.
One of the most important concepts in dependence modeling is knowing its dependency structure which illustrates the functional relationship between dependent random variables. According to Nelsen [27], we recall some of them here that will be used in this paper.
(i) The random variable X is right tail decreasing in Y (denoted by RTI(Y | X)) if and only if v − D(u, v) 1 − u is non-increasing in u for all v ∈ (0, 1).
(ii) The random variable X is said to be stochastically increasing (denoted by SI(X | Y)) in Y if and only if ∂ ∂v D(u, v) is non-increasing in v.
(iii) The random vector (X, Y) is said to be right corner set increasing (denoted by RCSI(X, Y)) if and only if P(X > x, Y > y | X > x , Y > y ) is non-decreasing in x and y for all x and y or equivalently when ∂ 2 ∂u∂v log(D(u, v)) ≥ 0.
Also, nowadays, researchers are interested in bivariate symmetry that has various applications in different fields of science. Klement and Mesiar [17], Nelsen [28], Durante [5], Dehghani et al. [3] and Alikhani-Vafa and Dolati [1] studied bivariate symmetry and its applications in statistical models. In the bivariate case, some authors defined bivariate symmetry as ex-changeability that is D(u, v) = D (v, u), for all u, v ∈ (0, 1). Nelsen [26] argued that radial symmetry is an appropriate definition for bivariate symmetry which is defined as D(u, v) =D(u, v), for all u, v ∈ (0, 1).
Dehghani et al. [3] constructed measures of radial asymmetry based on L 2 distance of D(., .) andD(., .) and the L 2 distance of the diagonal sections of D(., .) andD(., .) that are given by In the literature, several methods have been proposed to construct flexible copulas for dependent random variables. Marshall and Olkin [22] considered the component-wise maximum (minimum) of N independent and identical bivariate random vectors when N has a geometric distribution and constructed their corresponding bivariate copula. Many authors used this method to construct new flexible dependent families of bivariate distributions or bivariate copulas; see, for examples, Mirhosseini et al. [23], Kundu and Gupta [21], Zhang et al. [36] and Roozegar and Nadarajah [31]. Another method was proposed by Durrleman et al. [7] to construct a new bivariate flexible family of copulas that is applying a distortion function on a bivariate copulas. Some more recent developments of this method can be found in Genest and Rivest [11], Morillas [25], Durante et al. [6] and Dolati et al. [4].
In this paper, a new dependent model motivated by the structure of two dependent series-parallel systems with random number of series sub-systems that have fixed components connected in parallel is introduced. It's dependence properties are studied. It should be noted that this model is useful for analyzing environmental data which contains extreme values. Furthermore, we focus on estimating the parameters of the lifetime distribution of components based on lifetime data of systems. Many authors such as Balakrishnan et al. [2], Ng et al. [29], Zhang et al. [37], Kulkarni and Rajarshi [19] and Yang et al. [35] studied estimation problems with this point of view in the case of independence and Eryilmaz [8] exploited this kind of estimation in the case of dependence. Note that there are a lot of papers that have studied estimation methods for dependent models, see for example, Sajid et al. [32] (for a bivariate discrete model) and Kulkarni et al. [38] (for positively quadrant dependent).
The remainder of the paper is organized as follows. The proposed model is introduced in Section 2. In Section 3, two estimation methods, the method of moments and the maximum likelihood method are used for estimating parameters of the proposed model. Based on simulation results, the performances of the estimators are evaluated in Section 4. In Section 5, real data applications are revisited.

The proposed model
In these decades studying the corresponding distribution of systems with dependent components is one of the most interesting topics in distribution theory and reliability. It is known that systems are designed with various structures. A system with a series structure succeeds if all of its components work correctly and a system with a parallel structure succeeds if at least one of its components works correctly. The structure of series-parallel (parallel-series) systems are composed of fixed numbers of parallel (series) sub-systems that are connected in series (parallel). Several models have been proposed for systems lifetimes. In this paper, two series-parallel systems with random number of sub-systems that have fixed number of components are considered.

Example 2.3 (Geometric distribution):
Suppose that N has a geometric distribution that where η m (x, y) is given in (4).

Dependence properties
In recent years, copulas have been used for describing the dependence between dependent random variables. In this section, the corresponding copula function of (2) is presented and some of its dependence properties are studied. In addition, measures of radial asymmetry proposed in Section 1 are studied. By exploiting Sklar's Theorem forH T 1 ,T 2 (., .) in equation (2), the corresponding survival copulaĈ g θ (., .) is given bŷ where D(., .) is the corresponding copula of F(., .). In the special case, if D(., .) is an extremevalue copula, thenĈ g θ (., .) simplifies aŝ In dependent models, it is important to know the attitude of two relative random variables, and the concepts of dependency are used for this purpose. The next proposition illustrates conditions that (7) is RTI, SI and RCSI.
Proposition 3.1: Let (T 1 , T 2 ) be a random vector distributed as (7). Then is decreasing in t 2 for all t 1 , is increasing in t 2 for all t 1 , The next proposition displays the measures of radial asymmetry as introduced in Section 1, for our proposed model in (7).

Some special cases
In this subsection, dependence structures and radial asymmetry measures are investigated for model (7) when N has some distribution as given in Examples 2.
In the special case, letD(t 1 , t 2 ) = t 1 t 2 then RTI(T 1 | T 2 ), SI(T 1 | T 2 ) and RCSI(T 1 , T 2 ). Furthermore, by applying Proposition 3.2 the corresponding radial asymmetry measures are given by It should be noted that ifD(x, y) = xy, then ψ = J = 0 which means that the model is always radially symmetric. Case II: (Shifted Bernoulli distribution) Let N be as in Example 2.2. Then, applying Proposition 3.1 we conclude that In the special case, considerD(t 1 , t 2 ) = t 1 t 2 then RTI(T 1 | T 2 ), SI(T 1 | T 2 ) and RCSI(T 1 , T 2 ). In addition, by exploiting Proposition 3.2 the corresponding radial asymmetry measures are found to be Figure 1 illustrates the measures of radial asymmetry whenD(x, y) = xy. As it is observed ψ and J first increase and then decrease as θ increases. It means that by increasing θ the model gets more asymmetric and after reaching θ = 0.6 it tends to be more symmetric.
Case III: (Geometric distribution) Let N be as in Example 2.3. Then, by Proposition 3.1 we have In the special case, letD(t 1 , t 2 ) = t 1 t 2 then RTI(T 1 | T 2 ), SI(T 1 | T 2 ) and RCSI(T 1 , T 2 ). Moreover, exploiting Proposition 3.2 the corresponding measures of radial asymmetry are given by It is seen that ψ and J increase as θ increases. It means that by increasing θ the model tends to be more asymmetric.

The method of moments
This method may be the oldest and simplest method for finding point estimators. It is known that by this method, the population moments are estimated from sample moments. In what follows, by applying the method of moments based on the observed data, the parameters are estimated. For estimating the dependency parameter θ , the corresponding Kendalls tau moment estimator is given bỹ Suppose that τ (θ) is the corresponding Kendall's tau ofC g θ (u, v) by applying the method of moments, the moment estimatorθ is the root of the following equatioñ Furthermore, by exploiting Theorem 7.1 of Hoeffding [14] asymptotic normality ofτ is found to be and (T * 1 , T * 2 ) is an independent copy of (T 1 , T 2 ). The corresponding asymptotic variance (9) via the copula is expressed in the Supplemental Material. It should be noted that by applying the dellta method as given in Gut [15, p. 349], the moment estimatorθ follows an asymptotic normal distribution.
In the next step the marginal parameters α r for r = 1, 2 are going to be estimated. For this purpose the marginal distribution of T r for r = 1, 2 are required. By (2) the corresponding marginal distribution is given bȳ Based on (10), the sth moment of T r is found to be By plugging the corresponding moment estimator of the dependency parameterθ in equation (11) and applying the method of moments, the sth moment estimator of α r ∈ (0, +∞) for r = 1, 2 and s ≥ 1 is given bỹ (12), then the moment estimator of α r for r = 1, 2 is found to beα where Proof: Applying the definition of the probability generating function, we have Let y = w m in (14), then By exploiting the series ratio test, the series +∞

Example 4.1:
Let N be as in Example 2.1, andC g θ (u, v) be the Gumbel copula given bŷ Thenθ = n 0 and in view of (12) we havẽ

Example 4.2:
Let N be as in Example 2.2, andC g θ (u, v) be the Ali-Mikhail-Haq (AMH) copula given bŷ Then it's corresponding Kendall's tau is found to be see, Kumar [20]. By applying the method of the moments, the moment estimator of θ , saỹ θ, is the root of the following equatioñ It should be noted that for 0 < θ < 1, τ AMH (θ ) is a non-decreasing function with respect to θ and 0 < τ AMH (θ ) < 0.33, by applying the intermediate value Theorem, equation (16) has only one root ifτ < 0.001 orτ > 0.33. After estimatingθ by exploiting (12) the components estimator α r for r = 1, 2 is given bỹ andθ is the unique root of (16).

Example 4.3:
Let N be as in Example 2.3 andC g θ (u, v) be the AMH copula. Then, similar to Example 4.2, the moment estimatorθ is the unique root of (16) and the components estimators α r for r = 1, 2 is given byα It is interesting to note thatα r for r = 1, 2 is unbiased, that is E(α r ) = α r . Moreover, the corresponding variance ofα r for r = 1, 2 is found to be that is less than any other linear unbiased estimator. Thus, the moment estimatorα r , is the best linear unbiased estimator (BLUE), too. Furthermore, the covariance betweenα 1 and α 2 is given by this shows that there is a straight relation between the covariance of (T 1 , T 2 ) and (α 1 ,α 2 ). Furthermore, we readily find that Cov(α 1 , Furthermore, exploiting the multiple central limit theorem implies the asymptotic normality of (α 1 ,α 2 ) as k → ∞, that is Note that the asymptotic distribution ofα i for i = 1, 2 is given by , and an asymptotic 100 where z q is the qth upper quantile of the standard normal distribution.

Maximum likelihood method
In dependence modeling once the copula function is chosen, three maximum likelihood methods can be applied to calibrate the model and estimate the parameters of the copulas.
In the first method the marginal and dependency parameters are estimated simultaneously. The second method is inference function for margins (IFM) which the dependency and marginal parameters are estimated in two stages, see for example, Joe [15]. The third method is the pseudo maximum likelihood method (PML) that the empirical distribution functions are used to estimate the marginals and after that the dependency parameters are estimated, see, for example, Genest et al. [10]. It is known that the maximum likelihood method finds the parameter value that maximizes the log-likelihood function when the observations are given. In view of the copula-based parametric model of the random vector (T 1 , T 2 ) in (2), we have Then, the corresponding bivariate density function of (17) is given by where H r (t) = 1 − g θ (1 − (1 − e −t/α r ) m ) for r = 1, 2 and the corresponding pdf of T r for r = 1, 2 is found to be Moreover, the corresponding log-likelihood function of (18) is given by The MLE of α 1 , α 2 and θ are estimated by maximizing (19) with respect to α 1 , α 2 and θ call themα 1 ,α 2 andθ , respectively. The observed Fisher information matrix ofα 1 ,α 2 and θ is given by . (20) We recall that I −1 (α 1 ,α 2 ,θ) can be used as the asymptotic variance-covariance matrix of the MLE of the parameters. Thus, by applying the asymptotic properties of MLE as k → ∞, we conclude that It should be noted that the asymptotic distribution ofα i for i = 1, 2 is found to be and an asymptotic 100(1 − γ )% confidence interval for α i , i = 1, 2 is given by where z q is the qth upper quantile of the standard normal distribution. Furthermore, the asymptotic distribution ofθ is given by where Var(θ) = − ∂ 2 l(α 1 , α 2 , θ) ∂θ 2 −1 and the asymptotic 100(1 − γ )% confidence interval for θ is found to be ⎛ For evaluating the efficiency and the attitude of the moment estimatorsα 1 ,α 2 andθ and the maximum likelihood estimatorsα 1 ,α 2 andθ , we use numerical methods that are investigated in the next section by applying Monte Carlo simulation.

Simulation results
To evaluate the performance of moment estimators and maximum likelihood estimators of the preceding section a Monte Carlo simulation study is carried out based on 10,000 replications for different sample sizes. The mean squared error (MSE) and bias of the estimators have been calculated to evaluate the efficiency of the estimators. Furthermore, the corresponding bootstrap confidence interval of the estimators is proposed. We use the following algorithm to generate data from (T 1 , T 2 ).
are the random vector generated from the proposed model.

MSE and bias of estimators
In this subsection, the corresponding MSE and bias of moment estimators and maximum likelihood estimators of Examples 2.1, 2.2 and 2.3 are investigated for different sample sizes when m = 10, 20. Figures 3 and 4 display the MSE and bias of the moment and the maximum likelihood estimators of α 1 , α 2 and θ when N has a pmf as given in Example 2.1 when n 0 = 2, α 1 = 2 and α 2 = 5 for m = 10 and m = 20, respectively. Note that the corresponding loglikelihood is presented in the Supplemental Material. It is observed that the MSE of the moment estimators and the maximum likelihood estimators decrease when the sample size increases as we expected. In addition, the MSE of the moment estimators are less than  the corresponding MSE of maximum likelihood estimators. Thus, the moment estimators are better than maximum likelihood estimators. Moreover, the maximum likelihood estimators are more biased than moment estimators. Figures 5 and 6 show the MSE and bias of the moment and the maximum likelihood estimators when N has a pmf as given in Example 2.2 and α 1 = 0.5, α 2 = 1.5, θ = 0.75 for m = 10 and m = 20, respectively. The corresponding log-likelihood function is presented in the Supplemental Material. The MSE of the moment estimators and the maximum likelihood estimators decrease when the sample size increases in both cases that m = 10, 20. In addition, for m = 10, 20 the corresponding MSE of the moment estimatorsα 1 and α 2 are less than the corresponding MSE of maximum likelihood estimatorsα 1 andα 2 , respectively. So, we conclude that they are better estimators. Furthermore, the maximum likelihood estimatorθ is better than the moment estimatorθ for m = 10, 20. Moreover, it should be mentioned that the maximum likelihood estimators are more biased than the moment estimators.    Furthermore, for m = 10, 20 the moment estimatorsα 1 andα 2 are better than the maximum likelihood estimatorsα 1 andα 2 , respectively. In addition, for sample sizes smaller than 30 the maximum likelihood estimatorθ is better than the moment estimatorθ but for sample sizes more than 30 the moment estimatorθ is a better choice. Note that the maximum likelihood estimators are more biased than the moment estimators.

The confidence interval estimators
The confidence interval estimators of α 1 , α 2 and θ and their corresponding length are studied in this subsection for different sample sizes.  Tables 1 and 2 illustrate the 95% and 99% confidence interval estimators of α 1 , α 2 and θ when N has a pmf as given in Example 2.1 when n 0 = 2, α 1 = 2, α 2 = 5 for m = 10 and m = 20, respectively. In addition, the corresponding length of 95% and 99% confidence interval estimators are displayed in Figures 9 and 10 for different sample sizes for m = 10 and m = 20, respectively. As we can see the length of the 95% and 99% confidence intervals decrease by increasing sample sizes. Moreover, both the 95% and 99% confidence intervals of the moment estimators are usually better than the maximum likelihood estimators since their corresponding length are smaller. But, it should be noted that for m = 10 and sample sizes less than 10 the corresponding confidence intervals of the maximum likelihood estimatorsα 1 andα 2 are better than their corresponding confidence intervals of the moment estimatorsα 1 andα 2 , respectively. In addition, for m = 20 and sample sizes less than 20 the corresponding confidence interval of the maximum likelihood estimatorsα 1 andα 2 are better choices than the corresponding confidence intervals ofα 1 andα 2 , respectively.
The 95% and 99% confidence interval estimators of α 1 , α 2 and θ when N has a pmf as given in Example 2.2 when α 1 = 0.5, α 2 = 1.5 and θ = 0.75 are displayed in Tables 3  and 4 for m = 10 and m = 20, respectively. Furthermore, the corresponding length of the   Figures 11 and 12 for m = 10 and m = 20, respectively. As we expected by increasing sample sizes the corresponding length of the 95% and 99% confidence interval decrease. Moreover, the corresponding length of the 95% and 99% confidence intervals of the moment estimatorsα 1 andα 2 are better than the maximum likelihood estimatorsα 1 andα 2 . But, the corresponding length of the 95% and 99% confidence intervals of the maximum likelihood estimatorθ is smaller than the moment estimatorθ so it is a better estimator. Tables 5 and 6 displays the 95% and 99% confidence interval estimators of α 1 , α 2 and θ when N has a pmf as given in Example 2.3 when α 1 = 0.5, α 2 = 1.5 and θ = 0.75 for m = 10 and m = 20, respectively. Furthermore, the corresponding length of the 95% and 99% confidence intervals are illustrated in Figures 13 and 14 for different sample sizes when  m = 10 and m = 20, respectively. Note that the corresponding 95% and 99% confidence interval length of the moment estimatorsα 1 andα 2 are smaller than the corresponding 95% and 99% confidence interval length of the maximum likelihood estimatorsα 1 and α 2 , respectively. In addition, the corresponding 95% and 99% confidence interval of the maximum likelihood estimatorθ are better estimators than the corresponding 95% and 99% confidence interval of moment estimators since their lengths are smaller.

Real data analysis
In this section, we illustrate the proposed estimation procedure explained in the previous sections for two environmental data sets.

Temperature data
In order to illustrate the proposed procedure, we selected a real data set from the US national centers for environmental information publicly. These data are available from https://www.ncdc.noaa.gov. To establish the assumptions of the model proposed in Example 2.1, we considered the maximum temperature of the first day of each month from 1940 to 2018 recorded at two nearby US weather stations. Station I is La Guardia Airport in New York (Z (1) i,j ) and Station II is Newwark Liberty International Airport in New Jersey (Z (1) i,j ) for i = 1, . . . , 79 and j = 1, . . . , 12. The easyfit software was applied for data of both stations and it was concluded that they both follow a Johnson S B (bounded system) distribution. For more details about Johnson S B distribution we refer to Johnson [16]. Furthermore, the i,j and Z (2) i,j for i = 1, . . . , 79 and j = 1, . . . , 12, respectively. Where (.) is the cdf of standard normal distribution. By exploiting R software for the transformed  data we conclude that the data of both stations are converted to exponential data with parameters 1.010981 and 1.007963, respectively. In addition, the P-values of the corresponding Kolmogorov-Smirnov goodness of fit test for the converted data are 0.3789 and 0.1232, respectively. Furthermore, the Q-Q plot, P-P plot, histogram and theoretical densities, the theoretical and empirical CDFs are drawn for both converted data sets in Figures 15 and 16, respectively, which provides further evidence that both converted data follow an exponential distribution.
Consider the converted maximum temperature of the first day of each month in a year from 1940 to 2018 in both stations and take the maximum value between them for each year and call them X i and Y i for i = 1, . . . , 79 in each station, respectively. On the other hand   Figure 17 which illustrates positive dependence and upper tail dependence between T 1 and T 2 . Now we are going to fit an appropriate copula to T 1 and T 2 . Since N is generated from a discrete uniform distribution the dependency parameter is more than 1. Moreover, there is tail dependence between T 1 and T 2 . Hence Joe, Gumbel, Galambos and HuslerReiss copulas can be appropriate suggestions. The P-value of the corresponding bootstrap goodness of fit test and the AIC criteria of the suggested copulas are displayed in Table 7. Since the P-values of the corresponding goodness of fit tests for Gumbel, Galambos and HuslerReiss copulas are more than 0.05 so they can be suitable for these data but the AIC criteria of Gumbel copula is less than the others so it is the best choice.
After fitting Gumbel copula to T 1 and T 2 we applied the proposed estimation process for estimating the moment estimators and the maximum likelihood estimators of α 1 , α 2 and θ. The moment estimatorsα 1 ,α 2 ,θ and the maximum likelihood estimatorsα 1 ,α 2 ,θ and their corresponding confidence intervals at the level of 95% and 99% are illustrated in Table 8.

Wind data
The annual maximum wind speed data from Whitmore and Gentleman [34]    .0028 are fitted to X i and Y i for i = 1, . . . , 76, respectively which confirms that the base random variables marginals follow exponential distributions. For more details about generalized exponential distribution we refer to Gupta and Kundu [13]. Furthermore, the P-values of the corresponding Kolmogorov-Smirnov goodness of fit tests are 0.667 and 0.2677, respectively. For more verification that X i and Y i for i = 1, . . . , 76 follow a generalized exponential distribution the Q-Q plot, P-P plot, histogram and theoretical densities, the theoretical and empirical CDFs are illustrated in Figures 18 and 19, respectively. In order to fit the proposed model in Example 2.2, the random variable N = 2 is generated from a shifted Bernoulli distribution with parameter 0.25. The minimum value between Vancouver and Regina is picked up as T 1 and the minimum value between Toronto and Montreal is picked up as T 2 for the years 1947 to 1984. The corresponding Kendall's plot and scatter plot of (T 1 , T 2 ) are displayed in Figure 20 as we can see the data are positively dependent but there is no tail dependence between them. Furthermore, N is generated from a shifted Bernoulli distribution thus the dependency parameter is always less than 1 and more than 0 so we have to choose a copula which it's dependency parameter is less than 1 and more than 0.   The AMH and FGM are recommended for these data and the corresponding P-value of the bootstrap goodness of fit tests and the AIC criteria are displayed in Table 9. The P-value's are more than 0.05 hence both copulas are suitable since the corresponding AIC of AMH copula is less than FGM so it is a better choice. Applying the expressed estimation process the moment estimators and the maximum likelihood estimators of α 1 , α 2 , θ and their corresponding confidence intervals at the level of 95% and 99% are illustrated in Table 10.

Conclusion
In this paper, a new distribution for bivariate data has been constructed which motivated using a component-wise extreme of independent bivariate vectors. Some of its dependency properties have been investigated. In addition, the method of moments and the maximum likelihood method have been used to estimate the unknown parameters. Furthermore, the performance of the estimators has been evaluated by simulation results. Moreover, the estimators have been applied for two environmental data sets. It should be noted that the proposed model can be extended for constructing a multivariate distribution. Let Z (r) i,j be the lifetimes of the jth component respected to the ith sub-system corresponding to system r for j = 1, . . . , m, i = 1, . . . , N and r = 1, . . . , s. Furthermore, assume that Z i,j ) for j = 1, . . . , m and i = 1, . . . , N follow a multivariate cdf F r 1 ,...,r q (., . . . , .) for r 1 , r 2 , . . . , r q = 1, . . . , s such that r 1 = r 2 = · · · = r q . Moreover, assume that (Z (1) i,j , . . . , Z  Investigating the properties of the proposed multivariate model and analyzing multivariate data requires an independent study, and the authors hope to report the findings in future work.