Modiﬁed Maximum Pseudo Likelihood Method of Copula Parameter Estimation for Skewed Hydrometeorological Data

: For multivariate frequency analysis of hydrometeorological data, the copula model is commonly used to construct joint probability distribution due to its ﬂexibility and simplicity. The Maximum Pseudo-Likelihood (MPL) method is one of the most widely used methods for ﬁtting a copula model. The MPL method was derived from the Weibull plotting position formula assuming a uniform distribution. Because extreme hydrometeorological data are often positively skewed, capacity of the MPL method may not be fully utilized. This study proposes the modiﬁed MPL (MMPL) method to improve the MPL method by taking into consideration the skewness of the data. In the MMPL method, the Weibull plotting position formula in the original MPL method is replaced with the formulas which can consider the skewness of the data. The Monte-Carlo simulation has been performed under various conditions in order to assess the performance of the proposed method with the Gumbel copula model. The proposed MMPL method provides more precise parameter estimates than does the MPL method for positively skewed hydrometeorological data based on the simulation results. The MMPL method would be a better alternative for ﬁtting the copula model to the skewed data sets. Additionally, applications of the MMPL methods were performed on the two weather stations (Seosan and Yeongwol) in South Korea.


Introduction
Hydrological phenomena have multidimensional characteristics; therefore, more than one variable is required to be considered in simultaneous analyses of the hydrological phenomena [1]. The multivariate probability model accounts for the multidimensional characteristics. However, traditional multivariate probability models that are widely used in the analysis of hydrological data were derived on the basis of univariate probability distribution [2][3][4][5][6]. They assume that random variables follow the distribution model that is utilized in the derivation of multivariate probability distribution model [7]. Hence, they provide limited performances for analysis of the random variables when each random variable follows a different distribution model.
The copula model has been widely adopted as an alternative to mitigate the limitation of the traditional multivariate probability model in the frequency analysis of multidimensional data. In field of the hydrology, the copula model has been used in the analyses of various hydrometeorological phenomena, such as rainfall, flood, drought, and groundwater. For modeling rainfall events, Water 2020, 12, 1182 3 of 23 was then performed for a case study of a real data set. Performances of the three tested methods were evaluated in the case study. Results of this study may provide an alternative to fit the copula model in multivariate frequency analysis of the hydrometeorological variables.
This study is organized as follows: Section 2 describes theoretical backgrounds about copula and conventional parameter estimation methods. Section 3 describes the MMPL method. Section 4 presents the simulation procedure and its results. Section 5 includes application contents for weather stations in South Korea. Section 6 includes a discussion on the results of this paper, and Section 7 presents the conclusions.

Copula Model
The copula model was first suggested by Reference [29]. A copula function is a joint multivariate distribution that binds each univariate distribution on [0,1]. The d-dimensional copula function can be expressed simply as C : [0, 1] d → [0, 1] and is defined in Equation (1): where H is the joint probability, C is the copula function, x i is a random variable, F i is the marginal distribution of the i-th random variable that can be considered as the cumulative distribution function (CDF) of i-th random variable, and u i is a non-exceedance probability of marginal distribution variable of the copula function. Archimedean copulas have been broadly and frequently used for modeling hydrometeorological data because it is easy to derive their probability functions and estimate their parameters [12,[30][31][32][33][34]. A bivariate Archimedean copula C is a solution of the functional equation where ϕ is the generator which should be continuous strictly decreasing in [0, 1] with ϕ(0) = ∞ and ϕ(1) = 0, and the inverse ϕ −1 should be strictly monotonic [8,35]. Among the copula models, the multivariate extreme value copulas were developed for modeling the dependence structure between rare events, such as extreme rainfall, flood, and drought [36]. The Gumbel copula model, one of the extreme value copula, is the most common choice to model the dependence due to its simplicity [37][38][39][40]. The selection of Gumbel copula would be a good choice to evaluate the performances of the proposed parameter estimation for the copula model in the current study because data of extreme events often has high skewness. Hence, the Gumbel copula model is used for the skewed hydrometeorological data in the current study. In addition, it is worth noting that the Gumbel copula model is the only extreme value copula in Archimedean family [35,41].
The bivariate Gumbel copula function and its generator function are given in Equations (3) and (4).
where θ is a copula parameter.

Conventional Parameter Estimation Methods
Depending on the assumptions made for the copula models considered, the copula estimation method can be classified into parametric (e.g., the maximum likelihood estimation), semiparametric (e.g., the MPL method), and nonparametric methods. The maximum likelihood estimations usually include the exact maximum likelihood (EML), IFM, and MPL methods. The EML method is considered a more accurate method than the IFM method, but the EML method requires high dimensional derivations and the solutions from optimization approaches sometimes tend to be stuck in local optima [42]. The NP method is available for a few copula models because there is no existence of an explicit form of a relation between the copula parameter and dependence measure, τ, for all the copulas [10]. Thus, the IFM and MPL methods are considered conventional estimation methods for the parameter of copula in this study.

Inference Function for Margin (IFM)
The IFM method is a parametric estimation method and is also known as a two-step maximum likelihood method. In the first step, the marginal distributions are fitted by any fitting method of univariate probability distribution for each random variable. In the second step, the parameter of the copula model is estimated by using the maximum likelihood method. In the bivariate case, the log-likelihood function of the copula model is given by Equation (5): where n is the number of data, andF(x) andĜ(y) are marginal probability distributions which are estimated in the first step, respectively. u i and v i are the non-exceedance probabilities from each marginal distribution, and c is the density copula function of the copula function C.
In the IFM method, application of an inappropriate distribution model for the marginal distribution can lead to inaccurate estimation of non-exceedance probabilities and incorrect estimation of the copula parameter.

Maximum Pseudo-Likelihood Method (MPL)
The MPL method is a semi-parametric method because it includes both parametric and nonparametric procedures. The MPL method also consists of two steps like the IFM method. The copula parameter θ is estimated by the maximum likelihood method, which is a parametric approach. However, u i and v i are estimated by order statistics of each sample data which is a nonparametric approach. Because parametric and nonparametric approaches are used in the MPL method, it is called a semiparametric method. In the bivariate case, the log-likelihood function of the MPL method is given by Equation (6) and the copula parameter θ is obtained by maximizing the log-likelihood function.
where the R i and S i are the ranks of each variable X and Y, respectively. The MPL method uses the empirical distribution function (Weibull plotting position formula) for each marginal distribution. Thus, the method can attenuate impact from the wrong assumption of marginal distribution type in parameter estimation of the copula model [43]. Because marginal distributions, obtained by the Weibull plotting position, are unable to represent all distribution types, application of the MPL method can cause loss of information of marginal distributions. Performance of this method is worse than the IFM and EML methods when marginal distribution information is available.

Modified Maximum Pseudo-Likelihood Method (MMPL)
In the original MPL method, the probability obtained by the Weibull plotting position formula is used for random variables of the copula model. As mentioned above, the Weibull plotting position formula cannot reproduce the proper characteristics of some type of distributions.
Many hydrometeorological variables follow skewed distribution types, such as Gamma, Weibull, Gumbel, and generalized extreme value (GEV). In the current study, the MMPL method is proposed in which the Weibull plotting position formula is replaced by several plotting position formulas which consider skewness of data set. The MMPL method proposed in this study is designated for modeling extreme events.
Many studies have examined plotting position formulas for extreme value analysis [44][45][46]. In this study, six plotting position formulas were derived based on the extreme value distribution which were utilized for the MMPL method on behalf of the existing (Weibull) plotting position formula of the original MPL method: (1) three formulas which employ additional parameters to consider the embedded skewness [47][48][49]; and (2) three formulas which employ sample coefficient of skewness in the plotting position formula [50][51][52]. The six plotting position formulas utilized in this study are shown in Table 1.

MMPL-GD (Goel and De) [51]
where P i : exceedance probability, i: order, n: sample size, and C s : coefficient of skewness.
Reference [48] suggested the general plotting position formula for various probability distributions and defined as Equation (7): where i is order, n is sample size, and α is the parameter of the plotting position formula. The Weibull formula takes α = 0 and reproduces for the uniform distribution. The GEV takes α = 0.44 [48]. The unbiased plotting position formula based on reduced variates was developed by Reference [47], which takes α = 0.4. Reference [46] proposed the plotting position formula for extreme value (EV) nested distributions, which takes α = 0.25. Reference [50] introduced an unbiased plotting position formula for the GEV distribution by using the probability weighted moment method to estimate the exact plotting positions. Similarly, Reference [51] suggested that the plotting position formula including a coefficient of skewness for the GEV distribution, and Reference [52] derived a plotting position formula using the adaptation of the theoretical reduced variates of the GEV distribution and suggested the plotting position formula by using a genetic algorithm.

Simulation Experiment
In order to assess the performance of the MMPL methods, the Monte-Carlo simulation experiment was performed for given conditions in a variety of correlation coefficient (dependence) between two variables, sample size, marginal distribution, and coefficient of skewness. Particularly, coefficient of skewness of marginal distribution was set to have a positive value. The population of marginal distribution is explained in detail in Section 4.1.1, Section 4.1.2.

Simulation Design
First, since the copula function is a joint distribution from the univariate distributions based on the dependence structure, the dependence between variables x and y is determined by Kendall's tau (τ). Kendall's tau is a measure of the correlation coefficient based on rank and can depict nonlinear dependence. Spearman's rho (ρ) is also a rank-based correlation coefficient; however, Kendall's tau is more robust in the case of small samples and less sensitive to ties in data [1]. Kendall's tau is defined as Equation (8): In the simulation, 0.10, 0.15, . . . , 0.25, . . . , 0.50, . . . , and 0.75 values of Kendall's tau, which cover from weak to very strong positive correlations, are utilized for the dependence of the random variables because hydrometeorological variables often have a wide range of positive dependences. The corresponding Gumbel copula parameters (θ) are 1.11, 1.18, . . . , 1.33, . . . , 2.0, . . . , and 4.0, respectively. These values can be calculated by the equation (τ = θ−1 θ ) between the Gumbel copula parameter and τ [35]. Four sample sizes, i.e., n = 30, 50, 100, and 200, are employed in the simulation. These numbers represent the variety of sample sizes for hydrometeorological variables, from a small to a large number of data. By using random sampling of the Gumbel copula, 5000 sets of (u i , v i ) are generated for each of the 56 cases (14×4; the number of dependence cases × the number of sample size cases).
A main advantage of the MPL is that it can decrease the impact of selecting the wrong distribution type for margin on parameter estimation of the copula model. The simulation cases should be designed to evaluate the performances of parameter estimation methods taking into consideration the cases with unknown and known marginal distributions. In this study, two distribution models are selected for sampling of the marginal distribution: (1) Wakeby and (2) GEV distribution models.
The parameter estimation of the copula model is performed by eight different methods: (1) IFM method, (2) original MPL method (Weibull plotting position formula), (3)-(5) MMPL methods with plotting position formulas without coefficient of skewness, and (6)-(8) MMPL methods with plotting position formulas including coefficient of skewness. The copula parameter θ is estimated by the maximum likelihood method for every data set. The flowchart of the simulation design is shown in Figure 1.

Case of Unknown Marginal Distribution
The simulation using the Wakeby distribution represents the case where the marginal distribution is unknown, while the simulation using the GEV distribution represents the case where the marginal distribution is known. The Wakeby distribution for marginal distribution can represent a broad range of statistical characteristics of data because of a large number of parameters in the Wakeby distribution. Owing to its flexibility, the samples from the Wakeby distribution can cover a number of distribution models and potentially be applied to many other fields of science and environment where a multi-parameter distribution is required to model the data [53]. Thus, the Wakeby distribution can be a good choice in the case of unknown marginal distribution. The Wakeby distribution is given by quantile function form as Equation (9) [54].
where F is a standard uniform random variable, and m, a, b, c, and d are the parameters. The five parameter sets of the Wakeby distribution are used to represent various coefficients of skewness: 0.0, 1.0, 2.0, 3.0, and 4.0. These parameter sets were suggested by Reference [55]. The parameter sets and the basic statistics (mean: µ, standard deviation: σ, coefficient of variation: CV, coefficient of skewness: C s , coefficient of kurtosis: C k ) of each case are shown in Table 2.  parameter and [35]. Four sample sizes, i.e., = 30, 50, 100, and 200, are employed in the simulation. These numbers represent the variety of sample sizes for hydrometeorological variables, from a small to a large number of data. By using random sampling of the Gumbel copula, 5000 sets of ( , ) are generated for each of the 56 cases (14 × 4; the number of dependence cases × the number of sample size cases).

No. Parameters Statistics
A main advantage of the MPL is that it can decrease the impact of selecting the wrong distribution type for margin on parameter estimation of the copula model. The simulation cases should be designed to evaluate the performances of parameter estimation methods taking into consideration the cases with unknown and known marginal distributions. In this study, two distribution models are selected for sampling of the marginal distribution: (1) Wakeby and (2) GEV distribution models.
The parameter estimation of the copula model is performed by eight different methods: (1) IFM method, (2) original MPL method (Weibull plotting position formula), (3)-(5) MMPL methods with plotting position formulas without coefficient of skewness, and (6)-(8) MMPL methods with plotting position formulas including coefficient of skewness. The copula parameter is estimated by the maximum likelihood method for every data set. The flowchart of the simulation design is shown in Figure 1.

Case of Known Marginal Distribution
Because the GEV distribution has been widely applied to frequency analysis of extreme value data, it is employed for the marginal distribution in the case of known marginal distribution. The probability density and cumulative distribution functions of the GEV distribution are defined by Equations (10) and (11) [56]. where x 0 , α, and β are the location parameter, scale parameter, and shape parameter, respectively. The five parameter sets of the GEV distribution are given by Table 3. In this case, the shape parameters are −0.2, −0.1, 0.0, 0.1, and 0.2, which covers various coefficients of skewness, and their corresponding coefficients of skewness are 0.25, 0.64, 1.14, 1.91, and 3.54, respectively, by using Equation (12).
where Γ is gamma function. Table 3. The parameter sets of the GEV distribution for the case assuming known marginal distribution.

No.
Parameters Statistics

Finding Appropriate Marginal Distribution
The marginal univariate probability distribution should be determined to induce the copula variable (u, v) for the IFM method. Thus, in this study, five univariate distribution candidates were examined for marginal distribution, which are the gamma (GAM), GEV, Gumbel (GUM), generalized logistic (GLO), and Weibull (WBU) for the case of unknown marginal distribution (Wakeby). The Kolmogorov-Smirnov (K-S) distance statistic D is calculated for every sample data set to choose the most appropriate probability distribution from the five probability distribution candidates. The K-S distance statistic D for a given cumulative distribution function F(x) is given as Equation (13).
where F n is the empirical distribution with sample size of n.
For the case of known marginal distribution, the GEV distribution is applied to the marginal distribution. For the MPL and MMPL methods, different plotting position formulas are applied for (X, Y) to infer (u, v), which are explained in Section 3.

Case of Unknown Marginal Distribution (Wakeby)
The selection ratios of all employed distributions for each parameter set are shown in Table 4. The selection ratios of the GEV and GLO distributions are 53.6% and 41.6% for the first parameter set, respectively. For the second parameter set, the GLO distribution is the most frequently selected. The WBU and GEV distributions led to the highest selection ratios for the third and fourth parameter sets, respectively. The selection ratio of GEV distribution is the highest for the fifth parameter set. Figures 2-4 present the averaged values of parameter estimates from 5000 sample sets for three dependence coefficients values (Figure 2: τ =0.25; Figure 3: τ =0.50; Figure 4: τ =0.75). These three dependence coefficients are selected to describe the performance of the parameter estimation methods from weak to very strong correlation for illustration purposes. The results of these three correlation coefficients hold for the other cases with similar correlation coefficients. In Figures 2-4, the x-axis Water 2020, 12, 1182 9 of 23 indicates sample size and the gray line indicates the true value of the copula parameter. The black line is the MPL method, the dotted-line is the IFM method, the blue lines are the three MMPL methods without coefficient of skewness, and the red lines are the three MMPL methods with coefficient of skewness.
Water 2020, 12, x FOR PEER REVIEW 9 of 23 The black line is the MPL method, the dotted-line is the IFM method, the blue lines are the three MMPL methods without coefficient of skewness, and the red lines are the three MMPL methods with coefficient of skewness.   In general, for all combinations of τ = 0.25 and τ = 0.50 in Figures 2 and 3, the parameter estimates by the MMPL methods are overall closest to the true value, while the MPL and IFM methods provide overestimation and underestimation for the copula parameter, respectively. The IFM method also fails to become closer to the true value of the copula parameter as sample size increases. In addition, for the combinations of τ = 0.75 in Figure 4, the MMPL methods underestimate the copula parameter, while the MPL method is closest to the true value as sample size increases. The MPL method shows a similar pattern for all combinations of different coefficients of skewness, while the MMPL methods do not. The IFM method generally provides the worst parameter estimates among all employed parameter estimation methods. In other words, the IFM method shows the severe underestimation of the copula parameter for many combinations of coefficient of skewness showing the parameter estimates below the limits of figures for eight combinations (1_0, 2_1, 3_0, 3_1, 4_0, 4_1, 4_2, and 4_3). In general, for all combinations of = 0.25 and = 0.50 in Figures 2 and 3, the parameter estimates by the MMPL methods are overall closest to the true value, while the MPL and IFM methods provide overestimation and underestimation for the copula parameter, respectively. The IFM method also fails to become closer to the true value of the copula parameter as sample size increases. In addition, for the combinations of = 0.75 in Figure 4, the MMPL methods  In the case of τ = 0.25 (Figure 2), the MMPL-K method generally leads to the best performance which provides the closest parameter estimates to the true value for 10 combinations in which the sum of the coefficients of skewness is larger than or equal to three, for example, 2_1, 2_2, 3_1, 3_2, 3_3, 4_0, 4_1, 4_2, 4_3, and 4_4 combinations. On the other hand, the MMPL-G method leads to the best performance for four combinations in which the sum of the coefficients of skewness is less than or equal to two, for example, 0_0, 1_0, 1_1, and 2_0 combinations. The parameter estimates by the IFM method are the closest to the true value when sample size is smaller than or equal to 50, but the IFM method converges to the true value of θ only when combinations of coefficient of skewness are 0_0, 2_0, 2_2, 3_0, and 3_2, while the other MPL and MMPL methods converge to the true value of θ for all combinations of coefficient of skewness. Note that the MMPL methods always show a better performance than that of the original MPL method for all combinations and sample sizes.
In the case of τ = 0.50 (Figure 3), the performance of parameter estimates from all methods show similar patterns to that obtained for the case of τ = 0.25. The MMPL-G method shows the best performance for four combinations in which the sum of the coefficients of skewness is smaller than or equal to 2 (0_0, 1_0, 1_1, and 2_0). Alternately, the MMPL-K method generally shows the best performance when the sum of the coefficients of skewness is larger than 4 even though it slightly underestimates the parameters when the sample size is larger than 50. The IFM method shows a worse performance than the previous case. The performance of the IFM method is only the best for a sample size of 30 in most combinations. As in the case of τ = 0.25, the MMPL methods show a better performance than that of the original MPL method for all combinations and sample sizes.
In the case of τ = 0.75 (Figure 4), in which random variables are unrealistically highly dependent, all the MMPL methods underestimate the copula parameter. The MMPL-A method leads to the best performance when sample size is less than or equal to 50, while the MPL method shows the best performance when sample size is larger than or equal to 100. The MMPL-A is generally the best MMPL method among the tested MMPL methods. As in the previous cases, the IFM method resulted in an underestimation for all coefficient of skewness combinations and shows very poor performance, especially for eight combinations (1_0, 2_1, 3_0, 3_1, 4_0, 4_1, 4_2, and 4_3) in which the parameter estimates lie below the lower limits of figures. The MPL method behaves well as sample size increases and shows the best performance when sample size is larger than or equal to 100.

Case of Known Marginal Distribution (GEV)
Figures 5-7 present the averaged value of parameter estimates from 5000 sample sets for three values of dependence coefficients ( Figure 5: τ =0.25; Figure 6: τ =0.50; Figure 7: τ =0.75) as shown in the case of unknown marginal distribution. The performance of the IFM method is generally better than the MPL and MMPL methods. This is the expected result due to the assumption of the IFM method that the IFM method fits the copula model with a known marginal distribution model. The MMPL methods always provide better performance than the MPL method for the shape parameter combinations of τ = 0.25 and τ = 0.50, while the MPL method shows better performance than the MMPL method when τ = 0.75 and larger sample size (n ≥ 100).  In the case of τ = 0.25 ( Figure 5), the IFM method shows the best performances for all shape parameter combinations except four combinations: 0.1_0.1, 0.2_0.0, 0.2_0.1, and 0.2_0.2. In these four combinations, the MMPL-K method shows to the best performance for large sample size, while the IFM method shows the best performance only for relatively small sample size (n ≤ 50). The IFM method generally shows the best performance when the sum of the shape parameters is smaller than or equal to 0.2, while the MMPL-K method shows the best performance when the sum of the shape parameters is larger than or equal to 0.2 and sample size is larger than 50. Among the MMPL methods, the MMPL-K method shows the best performance, except when the sum of the shape parameters is smaller than or equal to −0.1. For these combinations, the MMPL-G method shows slightly better performance than the MMPL-K method. Note that the original MPL method shows the worst performance for all combinations and sample sizes.  In the case of τ = 0.50 (Figure 6), the IFM method generally gives the closest parameter estimate to the true value for all shape parameter combinations, especially at small sample sizes. Similar to the case of unknown marginal distribution, overall performance of the MMPL-G and MMPL-K are better than the other tested MMPL methods. Even though the performance of the IFM method is greatly improved compared to the case of unknown marginal distribution, the MMPL methods provide closer parameter estimates to the true value than the IFM method for larger sample size and shape parameters. As in the case of τ = 0.25, the original MPL method shows the worst performance for all combinations. In the case of = 0.25 (Figure 5), the IFM method shows the best performances for all shape parameter combinations except four combinations: 0.1_0.1, 0.2_0.0, 0.2_0.1, and 0.2_0.2. In these four combinations, the MMPL-K method shows to the best performance for large sample size, while the IFM method shows the best performance only for relatively small sample size ( ≤ 50). The IFM method generally shows the best performance when the sum of the shape parameters is smaller than or equal to 0.2, while the MMPL-K method shows the best performance when the sum of the shape parameters is larger than or equal to 0.2 and sample size is larger than 50. Among the MMPL methods, the MMPL-K method shows the best performance, except when the sum of the shape parameters is In the case of τ = 0.75 (Figure 7), the IFM method generally provides better performance than the MPL and MMPL methods, which, respectively, overestimate and underestimate the copula parameter for all shape parameter combinations and sample sizes. The exception is that the MMPL-A method shows the best and second-best performance when sample sizes are 30 and 50, respectively. Among the MMPL methods, the MMPL-A method provides the best performance for all combinations of shape parameters and sample sizes even though it generally underestimates the copula parameter. As in case of unknown marginal distribution, the MPL method behaves well as sample size increases, and the performance is quite comparable to the IFM method, especially for large sample sizes.

Data and Application Methodology
To investigate the difference between performances of the conventional estimation methods and the proposed estimation methods of the Gumbel copula for real hydrometeorological data, hourly rainfall data of 64 weather stations in South Korea were utilized. The locations of the weather stations are shown in Figure 8. The data were collected from the Korea Meteorological Association (https://data.kma.go.kr/cmmn/main.do), and the record lengths are from 32 to 56 years. Annual maximum volume (AMV) events were extracted from hourly recorded data, and then the rainfall volume and rainfall duration were used as random variables of bivariate Gumbel copula for each weather station. The scatter plot of all AMV events and the box-plot of Kendall's tau (τ) from 64 stations are illustrated in Figure 9.     To determine the type of marginal distribution for each variable, the chi-square goodness-of-fit test at 5% significance level (α = 0.05) and K-S distance statistic D were employed. Five univariate probability distributions, namely GAM, GEV, GUM, GLO, and WBU, were used as marginal distribution candidates. For each variable of employed two stations, Table 5 presents the p-values of the chi-square test and K-S distance statistics for each distribution candidate. For the Seosan station, the GEV distribution was selected for both marginal distributions. For the Yeongwol station, the GEV and WBU distributions were selected for rainfall volume and rainfall duration, respectively. The selected marginal distributions cannot be rejected by chi-square test and provide minimum K-S statistics among all distribution candidates.  Table 6 presents the estimates of the copula parameters from all employed parameter estimation methods. For the Seosan station, the parameter estimates by the MPL and IFM methods are 1.352 and 1.290, respectively. The parameter estimates of the six employed MMPL methods range from 1.296 to 1.326. For the Yeongwol station, the parameter estimates by the MPL and IFM methods are 1.785 and 1.651, respectively. The parameter estimates of the six employed MMPL methods range from 1.696 to 1.737. Joint probability of the estimated copula from three methods (MPL, IFM, and MMPL-K) are illustrated in Figure 10, which represents the joint occurrence probability of rainfall volume and rainfall duration. Since the MMPL-K method is most frequently selected to give the best performance from the results of the simulation experiment, it is illustrated in Figure 10 among the employed MMPL methods. As the probability increases, the differences among quantiles by the three methods increase. For example, a rainfall event with 523.6 mm of rainfall volume and 174.7 h rainfall duration corresponds to a 100 year return period, which corresponds to a 0.99 non-exceedance probability, for each marginal distribution of the Seosan station. The joint occurrence probabilities of the rainfall event (523.

Discussion
The simulation experiment is designed by a variety of conditions, such as correlation, sample size, the unknown and known marginal distributions, and coefficients of skewness of data. The results of the simulation experiments showed, the proposed parameter estimation method led to improvements in the performances of parameter estimation for the copula model when margins follow skewed distribution models as compared to the conventional MPL method.
In the case of the unknown marginal distribution, the parameter estimate of the MMPL method led to the best performance, and the MPL and IFM methods follow. When marginal distribution information is unavailable, the MMPL and MPL methods have the advantage compared to the IFM method, because they fit the copula model regardless of the marginal distribution. Even though the IFM method shows better performance than the MPL and the MMPL methods at small sample size, the IFM method generally shows the worst performance for most combinations of coefficient of skewness. Moreover, it often shows severe underestimation of the parameter estimate for large sample sizes. In the case of known marginal distribution, the performance of the IFM method is

Discussion
The simulation experiment is designed by a variety of conditions, such as correlation, sample size, the unknown and known marginal distributions, and coefficients of skewness of data. The results of the simulation experiments showed, the proposed parameter estimation method led to improvements in the performances of parameter estimation for the copula model when margins follow skewed distribution models as compared to the conventional MPL method.
In the case of the unknown marginal distribution, the parameter estimate of the MMPL method led to the best performance, and the MPL and IFM methods follow. When marginal distribution information is unavailable, the MMPL and MPL methods have the advantage compared to the IFM method, because they fit the copula model regardless of the marginal distribution. Even though the IFM method shows better performance than the MPL and the MMPL methods at small sample size, the IFM method generally shows the worst performance for most combinations of coefficient of skewness. Moreover, it often shows severe underestimation of the parameter estimate for large sample sizes. In the case of known marginal distribution, the performance of the IFM method is expected to be better than the MPL and MMPL methods. The MMPL method provides the best performance in fitting the copula model in some cases with the large values of shape parameters and sample size based on the simulation experiments. The large value of the shape parameter can lead to large uncertainty in the parameter estimation of GEV distribution [57]. The performance of the MMPL method may be better than the IFM method for large shape parameters (large coefficient of skewness) because some employed plotting positions take coefficient of skewness into consideration.
The MMPL methods employed in this study are classified into two groups: adopting plotting position formula without coefficient of skewness and with coefficient of skewness. According to the results of simulation experiments, the MMPL-G method, which belongs to the former group, is considered the best parameter estimation method for data sets with small coefficients of skewness. The MMPL-K method, which belongs to the latter group, is considered the best MMPL method for data sets with large coefficients of skewness. The MMPL methods with plotting position formulas that consider coefficient of skewness might be a good choice when the sum of coefficients of skewness is larger than or equal to three based on the results of simulation experiments. Thus, the MMPL-K method would be the best alternative among the tested MMPL and MPL methods for fitting the copula model to the skewed data of hydrometeorological variables.
The MMPL method underestimates the copula parameter when the dependence of two variables is too strong. For example, when Kendall's τ is 0.75 (see Figures 4 and 7), the MPL and IFM methods provide good performances, except for the MMPL-A method, which, in this case, is considered the best method among the employed MMPL methods. Seventy-five hundredths is an unrealistically large coefficient at the applied 64 stations as shown in Figure 9. All Kendall's τ are smaller than 0.5 at the applied 64 stations. For the skewed data, such as extreme precipitation and flood, the data sets with a very strong correlation may be infrequent. Thus, the MMPL method would be an alternative to the copula fitting methods for skewed data of hydrometeorological variables.
In this study, the bias is employed as a measure of accuracy compared to given true values of the copula parameter in the simulation experiment. For application of the real data sets, the parameter of population copula is unknown. Thus the bias may not be a good evaluation measure for real data sets. To overcome this drawback, goodness-of-fit tests are often used, and many of them examine goodness-of-fit by comparing empirical and estimated models. Because the empirical copula of goodness-of-fit tests, such as Kolmogorov-or Cramer-von Mises-type statistics, are based on the Weibull plotting position, the conventional goodness-of-fit tests are inappropriate for the MMPL methods. The goodness-of-fit tests should be modified or proposed for measuring goodness-of-fit of the copula model fitted by the MMPL method in the future.

Conclusions
In this study, a new parameter estimation method for the copula model was proposed. The proposed method is a modified version of the MPL method (MMPL). Six MMPL methods were proposed, and their performances and applicability were evaluated for the bivariate Gumbel copula model via the simulation experiments and case study to the extreme precipitation events. According to the results of the current study, the conclusions are as follows: 1.
The MMPL methods suggested in this paper prefers to estimate the parameters in a multivariate frequency analysis using the copula model for hydrometeorological data. The MMPL methods provide better performance than the original MPL method when the values of Kendall's tau (τ) are moderate (τ = 0.25 and 0.50) in the simulation experiments regardless of unknown and known marginal distribution. However, for the case of τ = 0.75, which is very rare as shown in the applications, the original MPL method performs better than the MPL method, especially for large sample sizes showing convergence to the true value as the sample size increases. Among the MMPL methods, the MMPL-K and MMPL-G methods can be the best methods depending on the statistical characteristics of applied data.

2.
The original MPL method generally overestimates the parameters and shows the worst performance among the applied methods, but the parameter estimates converge to the true values as the sample size increases regardless of unknown and known marginal distributions and combinations of coefficient of skewness. However, this method shows comparably competitive to the best method, particularly when the sample size is large (n ≥ 100) and τ = 0.75. 3.
The IFM method, as expected, shows the worst performances except for small sample size of 30 and underestimates the parameters severely as the Kendall's tau increases in the case of unknown marginal distribution. However, in the case of known marginal distribution (GEV), the IFM method shows the best performance, while the MMPL-K method provides better performance than the IFM method when the sum of the shape parameters is larger than or equal to 0.2 in large sample sizes (n ≥ 50) for the cases of moderate Kendall's tau values. In addition, the IFM method generally performs the best even though Kendall's tau is unrealistically large (τ = 0.75).
Finally, the suggested MMPL methods generally perform better than the original MPL method, except with very strong and unrealistic Kendall's tau value (τ = 0.75), regardless of unknown and known marginal distributions. In addition, the MMPL methods perform marginally better than the IFM method when the sum of the shape parameters is larger than or equal to 0.2 in large sample sizes (n ≥ 50), except τ = 0.75, even in the case of known marginal distribution. In conclusion, the MMPL-K and MMPL-G methods are appropriate methods of fitting the Gumbel copula model for the case of moderate Kendall's tau values. Particularly, the MMPL-K method generally provides more accurate parameter estimates in the case of unknown marginal distribution when the sum of the coefficients of skewness is larger than or equal to three with n > 50, and in the case of known marginal distribution when the sum of the shape parameters is larger than 0.2 with n > 50.