Copula-based composite likelihood approach for frequency analysis of short annual precipitation records

Hydrologic frequency analysis is used to estimate the magnitude of a hydrologic event with a given frequency of occurrence or to estimate the frequency of occurrence of a hydrologic event with a given magnitude (Haan 1977; Rao & Hamed 1999). Sufficiently long data are essential for accurate parameter estimation. However, it is often difficult to obtain such long-term data, and hydrological design estimates based on data of insufficient length involve large uncertainties. Therefore, analyzing these uncertainties is of great significance for improving the calculation accuracy and reducing the uncertainties of hydrological design values.

In the past, most hydrologic frequency analyses were analyzed through the use of univariate distributions. Using an assumed probability distribution to estimate the frequency analysis of a hydrologic event, one can compute the statistics for various magnitudes, which can then be used for the planning, design, and management of water resource systems.

Research on the univariate flood frequency analysis has been reviewed by Bobée & Rasmussen (1995), Rao & Hamed (1999), Singh & Strupczewski (2002), and Vittal et al. (2015). Singh & Strupczewski (2002) classified frequency analysis methods into four groups: (1) empirical, (2) phenomenological, (3) dynamic, and (4) stochastic watershed in conjunction with Monte-Carlo simulation. In addition, the suggested distributions for fitting hydrological data and a number of parameter estimation methods were summarized by Rao & Hamed (1999). Two popular methods for regional frequency analysis are the index FFA (Flood Frequency Analysis) and regression analysis (Singh & Strupczewski 2002).

Some studies on annual maximum flow frequency analysis show that the empirical probability plot often displays two or more distinct segments. Using a single probability distribution smooth function to fit all the segments may produce a compromised curve and may lead to serious flood estimation errors for large return periods. In most cases, the observed data of the smaller censored flows have no influence on the fitting of the flood frequency curve. If the smaller flows produce some (very marginal) influence, the final fitting still depends very heavily on the larger flows. Wang (1990a, 1990b, 1996a, 1996b, 1997a, 1997b) proposed using partial probability weighted moments and higher probability weighted moments to estimate flood distribution and presents the potential merits of this method.

Recently, a vast amount of research has been dedicated to investigating frequency analysis under conditions of climate change, land use change and river regulation, acting individually or together (Strupczewski & Kaczmarek 2001; Strupczewski et al. 2001; Cunderlik & Burn 2003; Khaliq et al. 2006; Leclerc & Ouarda 2007; Villarini et al. 2009, 2010; Gilroy & McCuen 2012; Salas & Obeysekera 2014; Zeng et al. 2014; Li & Tan 2015; Vasiliades et al. 2015; Xiong et al. 2015).

The above-mentioned studies still require data over a sufficiently long period to obtain sound hydrological designs based on statistical analysis. However, measuring and monitoring hydrologic events involves high costs, and the global hydrological gauge network is still unevenly distributed, with most measuring stations located in a few developed countries and regions and particularly poor gauge coverage in developing countries. The station density is relatively low in developing countries, and the length of the hydrological data is limited, resulting in critical data gaps or representativeness problems. As a result of various natural conditions and human activities, such as destruction and migration of stations, construction of water resource projects, and river regulation, these deficiencies have been exacerbated over recent decades. Inadequate data lengths can produce large uncertainties in hydrological design estimates (Chowdhary & Singh 2010). Even when employing new parameter estimation techniques, such as partial probability weighted moments and higher probability weighted moments, it is difficult to improve the estimation accuracy for shorter data sequences.

Therefore, to overcome the problem of insufficient data for frequency analysis, additional information has been used to reduce the uncertainty in hydrologic estimates. For example, Vogel & Stedinger (1985) lengthened short hydrologic records by employing the cross-correlation among nearby longer records. They also developed improved unbiased estimators of the mean and variance of the streamflow at the sites with short records, and these estimators could obtain equal or lower variance than the simple at-site sample maximum likelihood estimates. In addition, to increase the precision of flood quantile estimators for a relatively short data series, historical or paleofloods data have been incorporated into flood frequency analyses (Stedinger & Cohn 1986; Reis & Stedinger 2005; Halbert et al. 2016), and joint cumulative distribution functions (bivariate or trivariate normal distributions, extreme value distributions, Gumbel mixed model and bivariate Gamma distribution) have been adopted to simultaneously consider data from adjacent stations (Escalante-Sanboval & Raynal-Villasenor 1998; Yue 2000, 2001, 2002; Yue et al. 2001; Escalante-Sandoval 2007).

The conventional multivariate approaches are limited due to the marginal distributions and the joint distribution belonging to the same family and the mathematical formulation becoming increasingly complicated with the increasing number of parameters (Chen et al. 2010; Li et al. 2013). To overcome these shortcomings, copula functions have been proposed and widely used in multivariate hydrological frequency analyses (Zhang & Singh 2006, 2007a, 2007b; Song & Singh 2010; Salvadori et al. 2011; Ma et al. 2013; Salvadori et al. 2013; Zhang et al. 2013, 2015; Fu & Butler 2014; Huang et al. 2015; Chen et al. 2016).

Chowdhary & Singh (2009, 2010) noted that only the concurrent periods of data sets have been used in either conventional or copula-based multivariate hydrological frequency analysis, whereas a majority of non-concurrent data in the extensive amount of staggered hydrological data remains unutilized for hydrological analysis. These non-concurrent data offer untapped potential for reducing parameter uncertainty by considering a multivariate framework which can simultaneously consider partially concurrent information in an integrated manner. However, such studies are currently still limited. Sandoval & Raynal-Villasenor (2008) constructed a trivariate logistic model with generalized extreme value distribution marginals to analyze samples with equal or unequal record lengths. Their aim was to reduce the bias or uncertainty in flood estimates by combining observations from several sites with an insufficiently long record to increase the available information and to provide a regional at-site design event. Chowdhary & Singh (2009, 2010) proposed a copula-based composite likelihood approach (CBCLA) to analyze various lengths of multivariate data and to reduce uncertainty in hydrologic design estimates. The results show that the CBCLA can provide necessary flexibility for admitting arbitrary marginals and can significantly reduce uncertainty in design flood quantiles for a relatively short flood series by utilizing associated downstream flood data. This approach can offset the impact of dwindling hydrological observation networks around the world by adding information that can be derived from existing networks.

Combining the CBCLA with the P-III distribution for short annual precipitation records has not yet been thoroughly investigated. The primary aim of this paper is to use the CBCLA to investigate the uncertainty of annual precipitation design values for a relatively short series fitted with the P-III distribution. In this study, the concurrent and exclusive parts of a multivariate data set with various lengths were combined in an integrated manner using copulas and the maximum likelihood function, and the uncertainty was quantified in terms of confidence intervals. Examples using simulated and measured data are presented in this paper. The precipitation sequences of four stations, namely Xi'an, Dali, Lintong and Huayin, in Weihe River basin (WRB) in China were selected as a case study. The data length at Huayin station was 28 years and is thus short series. The data lengths at the other three stations were 77, 53 and 50 respectively, and they are long series that can provide additional information. The confidence intervals and standard errors of the design values at Huayin station were derived and compared with the results from univariate estimations. The results show that the CBCLA can reduce the uncertainty in the precipitation design values for a relatively short data series.

The Archimedean copula family is a popular model for hydrologic analyses (Nelsen 2007) because it is easy to construct and it can be applied to both positively and negatively correlated hydrologic variables. Therefore, three one-parameter Archimedean copulas, Gumbel–Hougaard (GH), Clayton and Frank, were selected to determine the joint probability distribution of the precipitation variables in this study. The analytical expressions of these copulas have been well documented in previous pioneering studies (Shiau 2006; Zhang & Singh 2007a, 2007b).

The dependence parameters for these copulas were estimated using Kendall's tau τ (KTE) and maximum likelihood estimation (MLE) (Genest & Favre 2007). For the KTE, the dependence parameter can be estimated from the relationship between τ and ⁠. The MLE is obtained by maximizing the likelihood function, which involves the parameters of both the marginal distributions and the copula function.

The root mean square error (RMSE), Akaike information criterion (AIC) and Bayesian information criterion (BIC) were used to represent the bias of the probability distributions (Posada & Buckley 2004; Zhang & Singh 2007a, 2007b; Ma et al. 2013). Additionally, a bootstrap version of the K-S test and A-D test based on Rosenblatt's transformation was employed to test the goodness-of-fit for the copulas (Dobrić & Schmid 2007). The expressions and calculation steps of these statistics have been documented in previous studies (Song & Singh 2010).

The MLM is applicable to the parameter estimation of complex density functions. The arrangement of a typical composite event, comprising some overlapping and some exclusive periods of X and Y, is shown in Figure 1. and are the full lengths of the two data series, respectively. is the length of the concurrent period, and and are the lengths of the individual nonconcurrent periods. Correspondingly, is the concurrent event, and are the nonconcurrent events. These three events constitute a composite event e_C, and the likelihood function of e_C is called the composite likelihood function (Chowdhary & Singh 2010). Let represent the pdfs of marginal distributions, and be the bivariate pdf of (X, Y), which takes the same form as f(x,y) in Equation (2). Here, ⁠, and are the parameter vectors of these distributions, respectively; ⁠, and r represent the numbers of parameters; is the dependence parameter vector of the bivariate pdf; and is a parameter number. For the sake of brevity, the pdfs are denoted as f(x), f(y), and f(x, y), respectively.

The estimators contain the new parameters of the marginal distributions for both the short and long series and the dependence parameter of the copula function. These parameters are used in the calculations in the following sections.

To calculate the confidence intervals of the design values, the standard errors of the distribution parameters should be estimated first. An asymptotic variance-covariance matrix usually provides a good description of the parameter error (Rueda 1981). In this section, following Chowdhary & Singh (2010), we summarize the procedures of the variance-covariance matrix for parameter estimation.

The asymptotic variance-covariance matrix of parameter estimation is the inverse of the Fisher information matrix, which has elements that are equal to the expectations of the second partial derivatives of the log-likelihood function distribution with respect to the parameters (Kamwi 2005; Nadarajah 2006). The main diagonal elements of the variance-covariance matrix correspond to the asymptotic variances of the parameters.

After deriving the standard errors of the design values estimated using Equations (24) and (25), the confidence intervals under confidence level were calculated using ⁠, where is the quantile of the standard normal distribution for confidence levels equal to ⁠, is the design value for the return period T which is calculated using and ⁠, and is the standard error (⁠ and ⁠) of ⁠.

In this section, CBCLA was used in a simulated series to validate its applicability for improving the accuracy of parameter estimation for shorter series. Assuming that two data series, X and Y, were P-III distributed and that the joint probability of their concurrent part was constructed following a Frank copula, the Monte-Carlo method was used to generate a data set with length N = 80 and the following distribution parameters: for X; for Y; and the Frank copula dependence parameter is ⁠, which corresponds to a Kendall correlation coefficient of ⁠. Based on the arrangement of the composite event shown in Figure 1, we constructed eleven different composite event cases S1, S2, S3, S4, S5, S6, S7, S8, S9, S10 and S11, as shown in Table 1, using the simulated data. Here, variables X and Y are called the ‘longer series’ and the ‘shorter series,’ respectively. The objective here was to see whether using CBCLA decreases the uncertainty in the design values for Y.

The distribution parameters of Y in different composite cases were estimated using CBCLA, and the design values and 95% confidence widths corresponding to return periods of 50, 100 and 200 years were calculated based on these parameters. Table 2 compares the calculation results with the estimations from the univariate case. For both the univariate and composite case, the standard errors and confidence widths increased with increasing return periods. For the same size of Y, the standard errors and confidence widths obtained based on CBCLA were less than that of the estimations using the univariate likelihood method. For example, in the case of N_Y = 30, the standard errors and confidence widths in cases S1 and S2 are smaller than those in the univariate case. Therefore, integrating additional observations with a series of insufficient length can increase the available information, can obtain more accurate parameter estimations from CBCLA, and can improve the accuracy of the design values.

Table 2 shows that for a constant length of Y, the standard errors and confidence widths decrease as the length of the additional series X increases. For example, in the case of N_Y = 40, the standard errors and confidence widths in cases S3, S4 and S5 gradually decrease, which indicates that the reduction in uncertainty increases as the length of the additional series increases. Therefore, CBCLA can reduce the uncertainty in the design values for inadequate series and can provide a more precise design basis for engineering projects.

To illustrate the application of the CBCLA as outlined above, we analyzed the annual precipitation data collected in the WRB, which is the longest tributary of the Yellow River in China. The river originates from Gansu province and flows 818 km eastward into the Yellow River, and the catchment covers an area of approximately 134,800 km². The WRB is located in the transitional zone between arid and humid areas with an average annual precipitation of 572 mm. Due to the typical temperate continental monsoon climate, the total annual precipitation exhibits significant spatial and seasonal variations.

In this study, the observed annual precipitation series at four rainfall gauges (Figure 2), namely, Xi'an, Dali, Lintong and Huayin, were used as a case study. All the data were obtained from the National Climate Center of China Meteorological Administration. The observation periods at these four stations were 1932–2008, 1956–2008, 1959–2008 and 1981–2008, respectively. All these data have low first-order series correlation coefficients (−0.1140, −0.0785, 0.0480 and 0.0503 respectively) and no missing records. Anderson's test of independence showed that these data have an independent structure at a 90% confidence level. The characteristics of the annual precipitation data are listed in Table 3.

Among the four stations, Huayin (the shortest data series) served as the base station, and the other three stations were regarded as reference stations to provide additional information. Three composite events were constructed based on the four series and were labelled case I, case II, and case III. The relevant information on the annual precipitation data series and the arrangement of composite events are summarized in Table 4.

The P-III distribution was applied to fit each data series on a univariate basis. The MLM was used to estimate the distribution parameters. Using the estimated parameters, we compared the theoretical and empirical probabilities of the observed data (Figure 3). The empirical nonexceedance probabilities were calculated using the Weibull formula recommended by the Ministry of Water Resources of China (MWR 2006). The P-III distribution in Figure 3 shows a satisfactory fitting of the observed data series.

The K-S test and A-D test were used to evaluate the goodness-of-fit of the precipitation observations. The computed sample values, Dn and ⁠, are listed in Table 5. The critical values of Dn and are also shown in Table 5 for simulation times N = 5,000 and significance levels of α = 0.20, 0.15, 0.10, 0.05, 0.02, 0.01. Because all the sample statistics are smaller than the corresponding critical values at various significance levels, the P-III distribution is acceptable as a marginal distribution.

As shown in Table 4, the length of the concurrent periods in all composite events is 28 years. Common measures, such as the Pearson correlation coefficient γ, Kendall's τ, and Spearman's ρ, were used to investigate the dependence of the concurrent variables. The results of the correlation coefficients and the corresponding p-values shown in Table 6 reveal the presence of significant dependent relationships between the concurrent periods in cases I, II, and III.

The joint distribution was constructed using the GH, Clayton and Frank copulas. The dependence parameters estimated based on the MLE and KTE method are listed in Table 7, and the RMSE, AIC and BIC for three copulas are also compared in this table. According to the smallest values of these three statistics, the best-fit copula and parameter estimation methods for the three composite cases are GH (MLE), Clayton (KTE) and GH (KTE). Table 8 lists the goodness-of-fit results of the best-fit copulas for the three cases. Given a significance level of α = 0.05, all the test statistics based on the observed data are smaller than the corresponding critical values, indicating that GH, Clayton and GH copulas are acceptable for modeling the bivariate joint distributions of the concurrent data in cases I, II, III, respectively.

Considering P-III as a marginal distribution and the best-fit copulas selected above, the estimates of the distribution parameters and the variance-covariance matrix based on the CBCLA were obtained from Equations (8) and (25), respectively. The precipitation design values along with the associated standard errors and widths of the 95% confidence intervals for 10, 20, 50, 100, 200 and 500 year return periods based on these parameters and the dispersion characteristics are shown in Table 9. The precipitation design values, standard errors and 95% confidence interval widths derived for only a univariate basis are also given in Table 9.

In addition, the standard errors and confidence widths obtained from the composite cases were smaller than those derived from the univariate case (see Table 9), indicating that the CBCLA yields a more precise parameter estimation. Table 10 compares the differences in the uncertainty reductions for the standard errors and confidence widths under different situations. For example, for composite cases I, II, and III, the standard error and confidence interval for a return period of T = 50 years decreased 10.71, 2.48, and 2.73%, respectively, compared with the univariate case when using CBCLA.

Furthermore, the reductions in the standard errors and confidence widths in case I (10–20%) are greater than those in cases II and III (1.5–3%) due to the length of the exclusive period in the additional data. For the same concurrent period length, the exclusive period in case I is longer than those in cases II and III (Table 4). The increased information for a shorter series is more significant when the length of the associated data is longer. As a result, integrating the multivariate data set in the CBCLA is more advantageous than the univariate analysis.

Overall, the uncertainty in the design values decreases when using the CBCLA, which validates the merits of using this approach to improve the accuracy and precision of design values with insufficient data.

To summarize, this study presented a CBCLA for the parameter estimation of shorter series with a P-III distribution. Using this method improved the accuracy of the parameter estimation, as reflected in the confidence intervals of the design precipitation values. Four annual precipitation series in the WRB of China were selected as a case study. The following conclusions were drawn:

1. Three bivariate composite events were constructed based on four precipitation series of different lengths. The bivariate data set in each composite event consisted of two parts, i.e., the concurrent period and the exclusive part in the longer series. The marginal distributions of the concurrent and exclusive part were fitted using a P-III distribution, and the corresponding parameters were estimated using the MLM. The dependence structure of the overlapping records was modeled using Archimedean copulas, and a copula-based bivariate composite likelihood function for parameter estimation was established based on the likelihood function. Finally, the precision of the estimates calculated with this approach was compared with the precision of the univariate maximum likelihood estimates. The results show that the uncertainty in the precipitation design values due to inadequate data decreased and that this reduction in uncertainty became more significant as the length of the exclusive data series increased.

2. Monte-Carlo simulation results showed that using CBCLA with P-III marginals achieved improved design values for short sequences in terms of both standard error and confidence intervals.

3. The main advantage of this approach is that it considers the unused data in a longer data series to build a multivariate data set with different lengths, which cannot be implemented using conventional methods. This integration enhances the information for short series and yields improved parameter estimation and estimation precision. Furthermore, copulas were used to build joint distribution functions in this approach to avoid the assumption that the marginal distribution must belong to the same type, as in the traditional multivariate distribution. The marginal and bivariate distributions used the P-III distribution and Archimedean copulas, respectively. Further studies on this approach with other marginals and copula models related to other hydrologic designs should be conducted.

The present study is financially supported by the National Natural Science Foundation of China (Grant Nos 51479171, 51179160 and 51579059). The authors also wish to express their cordial gratitude to the editor and anonymous reviewers for their helpful comments which have greatly helped to improve the quality of this paper.