Fitting Copula Functions to Bivariate Drought Variables

Using the identified drought events, Figure 13.5 plots show the scatter plots of the paired random variables. Figure 13.5 clearly shows the positive dependence between drought severity and drought duration; between drought severity and drought interarrival time; as well as between drought duration and drought interarrival time.

Figure 13.5 Scatter plots of paired drought variables.

As introduced in Chapter 3, we will investigate the dependence separately from the marginals. The empirical marginals (i.e., using the Weibull plotting position formula) are applied to study the dependence among drought variables. The parametric marginal distributions, fitted in the previous section, will be applied for risk analysis through joint (and/or conditional) return periods. In addition, we will choose the Archimedean and meta-elliptical copulas for the bivariate drought analysis, while only meta-elliptical copulas will be chosen for trivariate drought analysis. Given the positive dependence structure between drought variables, the selected Archimedean copulas are the Gumbel–Hougaard (which belongs to the extreme value family), Clayton, and Frank copulas. One may also choose other Archimedean copulas as candidates. The Gaussian and Student t copulas are selected as candidates from the meta-elliptical family.

Applying the pseudo-MLE to the bivariate drought variables, Table 13.6 lists the estimated copula parameters as well as the GoF test results through the improved Rosenblatt transform following the procedures discussed in Section 3.8.3 (i.e., the SnB test; Genest, et al., 2007) with the following steps briefly outlined:

1. 
   

   i. For d-dimensional random variables u = [u₁, u₂, …, u_d]u=u1u2…ud; setting Z = [Z₁, Z₂, …, Z_d],Z=Z1Z2…Zd, where Z1=u1=F1x1;Z2=∂Cu1u2∂u1;…,Zd=∂Cd−1u1…ud∂u1∂u2…∂ud−1∂Cd−1u1…ud−1∂u1∂u2…∂ud−1The null hypothesis is that u = [u₁, u₂, …, u_d]~C_θ ⇒ Z = [Z₁, Z₂, …, Z_d]u=u1u2…ud~Cθ⇒Z=Z1Z2…Zd are independent.

   
   
2. 
   

   ii. Compute the test statistic using

SnB=n3d−12d−1∑i=1n∑k=1d1−Eik2+1n∑i=1n∑j=1n∏k=1d1−maxZikZjk(13.3).

1. 
   

   iii. Generate random variables from the fitted copula function with the same sample size as the sample data.

   
   
2. 
   

   iv. Reestimate the copula parameters from the tested copula function and recompute the test statistics.

   
   
3. 
   

   v. Repeat steps ii–iv for a large number of times (N) and approximate the P-value using

   
   

   Pvalue=1N∑k=1N1Sn,kB∗>SnB, where Sn,kB∗ represents the test statistic from step iv.

From Table 13.6 we obtain that (i) all the copula candidates from Archimedean and meta-elliptical families may be applied to model drought severity (S) and drought duration (D) as well as drought severity (S) and drought interarrival time (INT); (ii) the Gumbel–Hougaard and Guassian copulas are the only two copula functions that may be applied to model drought duration and drought interarrival time.

Joint and Conditional Return Period for Bivariate Drought Analysis

In this section, we will apply the Gumbel–Hougaard copula (which belongs to the extreme value family) and the Gaussian copula to investigate risk through joint and conditional return periods. Here, we will only focus on drought severity and drought duration.

Joint Return Period of Bivariate Drought Analysis

As discussed in Section 3.10.2, the bivariate joint return period may be represented with either the “AND” case or “OR” case. Here we will focus on the “AND” case only. Equation (3.139) in Chapter 3 can be revised as follows:

TANDsd=EINT1−FSs−FDd+C12FSsFDd(13.4)

where E(INT) represents the expected drought interarrival time in years, E(INT) = 0.723 year.

Using Equation (3.139), Table 13.7 lists the “AND” case joint return periods with the fitted parametric marginal distributions. To further illustrate the computation, we will show how to compute the return period for T(S>8000 cfs. day, D>30 days)TS>8000cfs.dayD>30days:

• 
  

  F_S(S < 8000) = 0.2035FSS<8000=0.2035 from the fitted log-normal distribution S~LN2(10.1992, 1.4614)S~LN210.19921.4614.

  
  
• 
  

  F_D(D < 30) = 0.1661FDD<30=0.1661 from the fitted Weibull distribution D~Weibull(203.80, 0.8903)D~Weibull203.800.8903.

  
  
• 
  

  F_{S, D}(S ≤ 8000, D ≤ 30) = C(0.2035, 0.1661; θ = 6.2015) = 0.1479FS,DS≤8000D≤30=C0.20350.1661θ=6.2015=0.1479 from the fitted Gumbel–Hougaard copula for drought severity and drought duration.

  
  
• 
  

  The exceedance probability:

  
  
  FS>8000D>30=1−FS8000−FD30+CFSFD=1−0.2035−0.1661+0.1479=0.7783
  
  
  
  
• 
  

  The “AND” case joint return period TS>8000D>30=0.7230.7783≈0.93yr.

Figure 13. 6 shows the Joint return period of the “AND” case for drought severity and drought duration.

Table 13.7. Joint return period of drought severity and duration (“AND” case).

		D (days) Gumbel–Hougaard copula						D (days) Gaussian copula
		30	120	365	520	700	1120	30	120	365	520	700	1120
	8000	0.93	1.35	3.88	7.23	14.52	69.01	0.93	1.35	3.88	7.23	14.52	69.01
	28000	1.48	1.55	3.88	7.23	14.52	69.01	1.48	1.57	3.88	7.23	14.52	69.01
S	92000	3.62	3.62	4.20	7.25	14.52	69.01	3.62	3.62	4.54	7.41	14.54	69.01
(cfs.day)	180000	7.48	7.48	7.51	8.29	14.58	69.01	7.48	7.48	7.71	9.38	15.34	69.03
	300000	14.64	14.64	14.64	14.68	16.47	69.01	14.64	14.64	14.67	15.36	19.48	69.65
	800000	71.44	71.44	71.44	71.44	71.44	79.63	71.44	71.44	71.44	71.45	72.02	103.25

Figure 13.6 Joint return period of the “AND” case for drought severity and drought duration.

Conditional Return Period of Bivariate Drought Variables

There are two commonly applied approaches to study the conditional return period:

T(X₁>x₁| X₂>x₂)TX1>x1X2>x2 and T(X₁>x₁| X₂ = x₂)TX1>x1X2=x2. Here, we will investigate both conditional return periods for drought severity and drought duration, with the use of drought duration as the conditioning variable.

T(S>s| D>d)TS>sD>d In this case, the exceedance conditional probability of S given D exceeding a given duration (d) can be written through the copula as follows:

PS>sD>d=PS>sD>dPD>d=1−FSs−FDd+C12FSsFDd1−FDd(13.5)

The conditional return period can then be written as follows:

TS>sD>d=EINT1−FDd1−FSs−FDd+C12FSsFDd(13.6)

Equation (13.5) may also tell whether there exists the right tail increasing (RTI) property. The RTI property exists if the conditional exceedance probability is a nondecreasing function of drought duration for all drought severity values.

Using Equations (13.5) and (13.6) and the Gumbel–Hougaard copula as an illustrative example, Table 13.8 lists the conditional exceedance probability and conditional return period. Figure 13.7 plots the conditional exceedance probability and conditional return period. From Table 13.8 and Figure 13.7, it is seen that the exceedance probability is a nondecreasing function of duration, i.e., with the increase of drought duration, the exceedance probability of S>s|D>d is nondecreasing. The RTI property indicates that it is more likely for the drought severity exceeding a given threshold conditioned on a higher drought duration than that conditioned on a lower drought duration. Using S > 8000 cfs.day in Table 13.8 as an example, we have the following:

PS>8000D>30<PS>8000D>120=PS>8000D>365=PS>8000D>520=P(S>8000∣D>700=PS>8000D>1120

From Table 13.8 and Figure 13.7, it is also seen that for a given drought duration, the exceedance probability decreases with the increase of drought severity. To illustrate the computation, we will show the procedure to compute P(S>8000| D>30)PS>8000D>30 and T(S>8000| D>30)TS>8000D>30:

• 
  

  Previously we have computed P(S>8000, D>30) = 0.7783PS>8000D>30=0.7783 for the “AND” case.

  
  
• 
  

  The exceedance conditional probability is as follows: PS>8000D>30=PS>8000D>30PD>30=0.77831−0.1661=0.933

  
  
• 
  

  The conditional return period is as follows:

  
  
  TS>8000D>30=EINT1−FD30PS>8000.D>30=0.7231−0.16610.7783≈1.11yr
  
  

Table 13.8. Conditional exceedance probability and conditional return period using drought duration as the conditioning variable.

		D (days) P(S>s| D>d)PS>sD>d						D (days) T(S>s| D>d)TS>sD>d
		30	120	365	520	700	1,120	30	120	365	520	700	1,120
	8,000	0.93	1.00	1.00	1.00	1.00	1.00	1.11	2.52	20.82	72.27	291.60	6586.07
	28,000	0.59	0.87	1.00	1.00	1.00	1.00	1.77	2.88	20.82	72.27	291.60	6586.07
S	92,000	0.24	0.37	0.92	1.00	1.00	1.00	4.34	6.75	22.54	72.46	291.61	6586.07
(cfs.day)	180,000	0.12	0.18	0.52	0.87	1.00	1.00	8.98	13.97	40.30	82.84	292.81	6586.08
	300,000	0.06	0.09	0.27	0.49	0.88	1.00	17.55	27.32	78.54	146.80	330.78	6586.36
	800,000	0.01	0.02	0.05	0.10	0.20	0.87	85.66	133.33	383.31	714.17	1434.63	7600.17

Figure 13.7 Conditional exceedance probability and conditional return period of S>s ∣ D>dS>s∣D>d.

T(S>s| D = d)TS>sD=d. In this case, the drought duration is the fixed conditioning variable. The exceedance conditional probability may be written as follows:

PS>sD=d=1−PS≤sD=d=1−CFSFD=FDd=1−∂C12FSFD∂FDdFDd(13.7)

The conditional return period can then be written as follows:

TS>sD=d=EINTPS>sD=d=EINT1−∂C12FSFD∂FDdFDd(13.8)

According to Nelson (2006), the stochastic increasing (SI) property exists if P(S>s| D = d)PS>sD=d is a nondecreasing function of drought duration for all drought severity values.

Using Equations (13.7) and (13.8) and the Gumbel–Hougaard copula as an illustrative example, Table 13.9 lists the exceedance conditional probability and conditional return period. Figure 13.8 plots the exceedance probability and conditional return period. Figure 13.8 clearly shows that for D = d, the exceedance probability P(S>s| D = d)PS>sD=d is a nondecreasing function of drought duration, i.e., P(S>s| D = d₁) ≤ P(S>s| D = d₂), d₁ < d₂PS>sD=d1≤PS>sD=d2,d1<d2.

Table 13.9. Exceedance conditional probability and conditional return period of S>s ∣ D = dS>s∣D=d.

		D (days) P(S>s| D = d)PS>sD=d					D (days) T(S>s| D = d)TS>sD=d
		30	120	365	520	700	30	120	365	520	700
	8000	0.36	0.99	1.00	1.00	1.00	2.02	0.73	0.72	0.72	0.72
	28000	0.00	0.29	1.00	1.00	1.00	286.86	2.45	0.72	0.72	0.72
S	92000	2.73E-06	0.00	0.57	0.98	1.00	2.65E+05	1.60E+03	1.27	0.74	0.72
(cfs.day)	180000	2.08E-08	3.45E-06	0.01	0.39	0.97	3.48E+07	2.10E+05	67.70	1.83	0.74
	300000	2.78E-10	4.61E-08	1.44E-04	0.01	0.43	2.60E+09	1.57E+07	5.01E+03	80.12	1.67
	800000	1.33E-14	2.19E-12	6.85E-09	4.32E-07	3.82E-05	5.43E+13	3.31E+11	1.06E+08	1.67E+06	1.89E+04

Figure 13.8 Exceedance conditional probability and conditional return period plot.

Dynamic Return Period for a Given Drought Episode

In the previous sections, we have investigated the joint and conditional return periods of bivariate drought variables, namely drought severity and drought duration. Following De Michele et al. (2013), we will investigate the evolution of drought within a given drought event (or simply called drought episode). It is worth mentioning that the copula function fitted to the drought severity and drought duration will not be applicable here. The empirical copula will be applied to study the dynamic return period for the given drought episode.

As discussed in De Michele et al. (2013), the dynamic return period is estimated through the Survival Kendall Distribution (also called DSKRP). As introduced in Section 4.5.1, the Kendall distribution may be considered as univariate realization of the copula function. In the case of bivariate analysis, the Kendall distribution may be simply written as follows:

and the survival Kendall distribution K¯Ct may be written as follows:

K¯Ct=PC¯FX1x1FX2x2≥t(13.10)

In Equation (13.10), C¯ represents the survival copula, and FXi¯xi=1−FXixi,i=1,2.

The DSKRP can then be written as follows:

TDSKRP=μ1−K¯Ct(13.11)

In Equation (13.11), μμ represents the average interarrival time of the drought event (μ = 0.723 yr)μ=0.723yr).

Furthermore, to investigate the DSKRP for a given drought episode, the average running drought intensity (I) will be applied. The average running drought intensity is computed as the average drought deficit starting from the initiation of a drought episode until day k into the drought. With this in mind, the new bivariate drought variable is given by pair as (I_k, k), k = 1, 2, …, mIkk,k=1,2,…,m, where m represents the total number of days of the drought episode. To illustrate the DSKRP method, we will use the recent 21-day drought episode identified (i.e., September 7–27, 2016) as an example. Table 13.10 lists the daily streamflow and streamflow deficit during this dry period.

As introduced earlier, the running average I_k can be computed as follows:

Ik=∑i=1kdeficitik,k=1,2,…,21.

For example, for the drought period of day 3, we have the following: I3=259.38+582.88+664.573=502.28cfs.

Figure 13.9 plots the flow deficit running average, as well as DSKRP and return period computed from the univariate flow deficit (RPFD) for the recent dry period. The plot on the left shows that the running average fluctuates within the drought episode. This fluctuation may reflect the severity of the state of drought on a given day. The plot on the right shows the DSKRP and the RPFD. It shows that within the drought episode, DSKRP and RPFD for the state of drought share a similar pattern and reflect the fluctuation of flow deficit. Table 13.11 lists the computed survival Kendall distribution and the corresponding DSKRP. As an illustration, here we will show how to compute DSKRP for (I₁, 1) = (259.38, 1)I11=259.381:

1. 
   

   1. To compute the exceedance empirical marginal probability, use Weibull plotting position formula as follows:

   
   
   FIk¯I1=1−FIkIk≤I1=1−121+1≈0.9545;F¯k1=1−121+1≈0.9545.
   
   
   
   
2. 
   

   2. To compute the empirical survival copula, i.e., C¯, apply the following formula:

   
   
   C¯FIkI1Fk1=FIk¯I1+F¯k1−1+CFIkI1Fk1=0.9545+0.9545−1+C122122=0.9545+0.9545−1+0.0476≈0.9567
   Here, C(F_Ik(I₁), F_k(1))CFIkI1Fk1 is estimated using the empirical copula formula given as Equation (2.59). The first pair (I₁, 1) = (259.38, 1)I11=259.381 is the smallest pair among all (I_k, k), k = 1, 2, …Ikk,k=1,2,…. and we have C122122=121≈0.0476.
   
   
   
3. 
   

   3. To compute the empirical survival Kendall distribution K¯C, Equation (13.10) is applied to compute the empirical survival Kendall distribution. Again using C¯FIkI1Fk1 as an example, we have C¯FIkI1Fk1≈0.957, which is the largest one among all the pairs, then, K¯C0.957=121+1≈0.045.

   
   
4. 
   

   4. To compute DSKRP and RPFD, DSKRP = E(INT)/1−K¯Ct; RPFD=E(INT)/[1-F_IkFIk].

Figure 13.9 Flow deficit running average and DSKRP for 2016 event from September 7–27.

Table 13.11. Results for estimating DSKRP and RPFD.

Date	IkIk	FIkFIk (cfs)	FkFk	Rank (Ik, k)Ikk	C(FIk, Fk)CFIkFk	ĈFIkFk	K¯Ct	DSKRP (yrs)	RPFD (yrs)
7-Sep-16	259.38	0.045	0.045	1	0.048	0.957	0.045	0.757	0.757
8-Sep-16	421.13	0.091	0.091	2	0.095	0.913	0.091	0.795	0.795
9-Sep-16	502.28	0.136	0.136	3	0.143	0.870	0.136	0.837	0.837
10-Sep-16	561.55	0.273	0.182	4	0.190	0.736	0.273	0.994	0.994
11-Sep-16	553.60	0.227	0.227	4	0.190	0.736	0.273	0.994	0.936
12-Sep-16	544.24	0.182	0.273	4	0.190	0.736	0.273	0.994	0.884
13-Sep-16	585.78	0.318	0.318	7	0.333	0.697	0.318	1.060	1.060
14-Sep-16	635.49	0.364	0.364	8	0.381	0.654	0.364	1.136	1.136
15-Sep-16	669.12	0.409	0.409	9	0.429	0.610	0.409	1.224	1.224
16-Sep-16	702.43	0.545	0.455	10	0.476	0.476	0.455	1.326	1.591
17-Sep-16	722.83	0.682	0.500	11	0.524	0.342	0.500	1.446	2.272
18-Sep-16	726.36	0.727	0.545	12	0.571	0.299	0.727	2.651	2.651
19-Sep-16	714.73	0.591	0.591	11	0.524	0.342	0.545	1.591	1.767
20-Sep-16	694.93	0.500	0.636	10	0.476	0.340	0.636	1.988	1.446
21-Sep-16	693.78	0.455	0.682	10	0.476	0.340	0.636	1.988	1.326
22-Sep-16	721.46	0.636	0.727	14	0.667	0.303	0.682	2.272	1.988
23-Sep-16	741.86	0.773	0.773	17	0.810	0.264	0.773	3.181	3.181
24-Sep-16	780.15	0.864	0.818	18	0.857	0.175	0.818	3.977	5.302
25-Sep-16	799.46	0.955	0.864	19	0.905	0.087	0.909	7.953	15.907
26-Sep-16	794.61	0.909	0.909	19	0.905	0.087	0.864	5.302	7.953
27-Sep-16	768.38	0.818	0.955	18	0.857	0.084	0.955	15.907	3.977

Marginal Distribution of Maximum Drought Intensity

In the previous section, we studied bivariate hydrological drought frequency analysis with the use of drought severity and drought duration as an example. In this section, we will study the trivariate drought frequency analysis by applying both vine and meta-elliptical copulas. The variables considered are drought severity (S), drought duration (D), and maximum drought intensity (MDI, i.e., the maximum flow deficit of a drought episode). With S and D fitted by log-normal and Weibull distributions, here we only need to investigate MDI. The histogram of MDI in Figure 13.10(a) clearly shows that its density function is skewed to the left with a long left tail. Thus, to reduce the complexity of fitting univariate distribution, the meta-Gaussian transformation is applied, which is the same as the preparation of marginals for the meta-elliptic copula approach as follows: Variable (XX)→empirical distribution (F_nFn, e.g., Weibull plotting-positon formula)→X^T~Φ⁻¹(F_n)XT~Φ−1Fn. The empirical distribution and density after transformation are shown in Figure 13.10(b)–(c).

Figure 13.10 Plots to study the maximum drought intensity (MDI): (a) histogram of MDI; (b) empirical distribution of DMI; (c) histogram and fitted N(0,1) for transformed MDI.

Vine-Copula Approach to Model Trivariate Drought Variables

Table 13.12 lists Kendall’s correlation coefficient for S, D, and MDI. As expected, all three drought variables are positively dependent. According to the degree of dependence, we will use severity as the center variable to construct the three-dimensional vine copula using the following structure: T₁ : D − S − MDI;  T₂ : D ∣ S − MDI ∣ ST1:D−S−MDI;T2:D∣S−MDI∣S. As introduced in Chapter 5, there are three bivariate copula functions involved in this analysis: (S and D); (S and MDI); and (D|S, MDI|S). All three bivariate copula functions are estimated separately and allowed to be fitted with different copula functions. In the bivariate analysis, we have investigated the Gumbel–Hougaard and Gaussian copulas to S and D. Applying the same copula candidates as the bivariate analysis, the best-fitted copula function will be selected to model the positively dependent S and MDI. Results for the copula candidates are listed in Table 13.13.

As shown in Table 13.13, all copula candidates may be applied to model S and MDI. Among all the copula candidates, the Frank, Gaussian, and Student t copulas yield very similar performance and outperform the Gumbel–Hougaard and Clayton copulas. In addition, the degree of freedom νν is estimated as 3.07E + 06; this high degree of freedom suggests the convergence to the Gaussian copula. Based on the MLE and SnB test statistics, the Frank copula (with the largest MLE and smallest SnB statistics) is chosen to model S and MDI. Figure 13.11 compares the pseudo-observations (empirical CDF) with those simulated from the fitted Frank copula. In Figure 13.11, we also provide a plot of comparison with observed variables in the real domain. Comparisons show the appropriateness to apply the fitted copula functions visually.

Figure 13.11 Comparison of pseudo-observations and real observations with those simulated from the Gumbel–Hougaard (S and D) and Frank (S and MDI) copulas for T1.

With Gumbel–Hougaard and Frank copulas chosen for T₁, the copula candidate for T₂ is selected, based on the conditional copulas (i.e., CS,DGHFDFS=FSs;CMDI,SFrankFMDIFS=FSs). From the computed conditional copulas, we compute the Kendall correlation coefficient for T₂ as follows: τnCS,DGHFDFS=FSsCMDI,SFMDIFS=FSs≈−0.27. With the negative Kendall’s tau computed, only Frank, Gaussian, and Student t copulas are used to evaluate for T₂, with the results listed in Table 13.14. The table reveals the following:

1. 
   

   i. All three copulas may properly model the variables for T₂.

   
   
2. 
   

   ii. The Student t copula again converges to the Gaussian copula.

   
   
3. 
   

   iii. SnB test statistics increase significantly, compared to those for T₁.

   
   
4. 
   

   iv. With the increase of SnB statistics, the corresponding P-value decreases.

   
   
5. 
   

   v. Based on the P-value as well as the computed MLE, the Gaussian copula is the fitted copula for T₂.

To this end, we have completed the construction of the vine copula for the drought variables as follows:

The joint distribution may then be computed using Equation (5.60). Figure 13.12 compares pseudo-observations (through parametric conditional copula) with simulations from the Gaussian copula of T₂. Figure 13.12 also compares the empirical copula and the parametric joint distribution from the vine copula. It is shown that the vine copula fits the trivariate drought variable reasonably well.

Figure 13.12 Comparison plots for T₂ and joint CDF.

Simulation from the Fitted Vine Copula

Following the simulation algorithms (Aas et al., 2009), Figure 13.13 shows the comparison of observations with drought variables simulated from the fitted vine copula. Here we will again illustrate how to simulate the random variable from the vine copula with the fitted GH–Frank–Gaussian copulas.

1. 
   

   1. Generate three independent, uniformly distributed random variables: w = [0.7372, 0.7869, 0.6537]w=0.73720.78690.6537, where

   
   
   
   

   w(1) = U(1); w(2) = C₁₂(U(2)| U(1)); w(3) = C_3 ∣ 12(U(3)| U(1), U(2))

   w1=U1;w2=C12U2U1;w3=C3∣12U3U1U2.
   In this example, U(1) = F_D(d); U(2) = F_S(s); U(3) = F_MDI(mdi)U1=FDd;U2=FSs;U3=FMDImdi.
   
   
   
2. 
   

   2. Compute U(2)U2, i.e., F_D(d)FDd:

   
   

   From step 1, we have U(1) = F_D(d) = w(1) = 0.7372U1=FDd=w1=0.7372 and C12U2U1=CSDGHFSFDd=0.7372=w2=0.7869.

   
   

   With the fitted Gumbel–Hougaard copula (θ = 6.2θ=6.2) for drought severity and duration, we will need to compute U(2) (i.e., F_S(s)FSs) using the following:

   
   
   CSDGHFSFDd=0.7372=∂CSDGHFDFS6.2∂FDFD=0.7372
   
   
   
   

   To be consistent with the discussion in Chapter 5, we assign the conditional copula as the h function:CSDGHFSFDd=0.7372=hGHFSFD6.2=0.7869. The conditional copula for the Gumbel–Hougaard copula is listed as No. 4 in Table 4.2. As seen in Table 4.2, U(2) needs to be solved for numerically using the root-finding (e.g., bisection method) technique. Using the bisection method, we have the following:

   
   
   
   

   U(2) = F_S(s) = h⁻¹(0.7869, 0.7372; 6.2) = 0.777.

   U2=FSs=h−10.78690.73726.2=0.777.
   
   
   
   
3. 
   

   3. Compute U(3)U3, i.e., F_MDIFMDI:

   
   

   From the vine structure, we have the following:

   
   
   w3=C3∣12U3U1U2=∂C13∣2C23U3U2C12U1U2∂C12U1U2=∂D,MDI∣SGaussianCMDI,SFrankFMDIFSs11.47CS,DGHFDFSs6.2−0.418∂CS,DGHFDFSs6.2
   and U(3)U3 can then be computed with the following three steps:
   
   
   
   
   
   1. 
      

      a. Compute the conditional copula CS,DGHFDFS. From the first two steps, we have F_D = 0.7372, F_S = 0.777FD=0.7372,FS=0.777; CS,DGHFSFD6.2 can then be computed by substituting F_D = 0.7372, F_S = 0.777, θ = 6.2FD=0.7372,FS=0.777,θ=6.2 into No. 4 conditional copula in Table 4.2. We obtain the following:

      
      
      CS,DGHFDFS6.2=C0.73720.7776.2=0.2791
      
      
      
      
   2. 
      

      b. Compute the conditional copula of CMDI,SFrankFMDIFSs11.47 with the use of the meta-Gaussian copula fitted to T₂ by setting the h function as follows:

      
      
      w3=hgaussianCMDI,SFrankFMDIFS11.47CS,DGHFDFS6.2−0.418=0.6537
      and we have the following: CMDI,SFrankFMDIFS11.47=hgaussian−10.65370.2791−0.418
      
      
      

      For the Gaussian copula, its conditional copula is the univariate normal distribution. The derivation of the conditional copula is given as Equation (7.42). In this particular problem, Equation (7.42a) can be rewritten as follows:

      
      
      hgaussianCMDI,SFrankFMDIFS11.47CS,DGHFDFS6.2−0.418
      
      
      ~ΦΦ−1CMDI,SFrankFMDIFS11.47−ρΦ−1CS,DGHFDFS6.21−ρ20.5
      LethMDI,S=CMDI,SFrankFMDIFS11.47 we have:
      
      Φ−1hMDI,S=Φ−1w31−ρ20.5+ρΦ−1(CS,DGHFDFSs6.2
      
      
      
      

      =(1 − (−0.418)²)^0.5Φ⁻¹(0.6537) − 0.418Φ⁻¹(0.2791) = 0.571

      =1−−0.41820.5Φ−10.6537−0.418Φ−10.2791=0.571
      and hMDI,S=CMDI,SFrankFMDIFS11.47=Φ0.571=0.716.
      
      
      
   3. 
      

      c. Compute F_MDIFMDI from CMDI,SFrankFMDIFS11.47.

      
      

      In steps b, we have computed CMDI,SFrankFMDIFS11.47=0.716. Using the conditional Frank copula in Table 4.2, we have the following:

      
      
      
      

      U(3) = F_MDI(mdi) = h⁻¹(0.716, 0.777; 11.47) = 0.8421

      U3=FMDImdi=h−10.7160.77711.47=0.8421.
      
      

To this end, we have finished one simulation:

[F_D, F_S, F_MDI] = [0.7372, 0.777, 0.8421]

FDFSFMDI=0.73720.7770.8421.

To repeat the preceding procedure, we can simulate the random variables of size N.

By applying the fitted marginal distributions we obtain the simulation in the real domain as follows:

S_simu = lognormal⁻¹(0.777, 10.1992, 1.4614²) = 8.19 × 10⁴cfs. day

Ssimu=lognormal−10.77710.19921.46142=8.19×104cfs.day

D_simu = weibull⁻¹(0.7372; 203.8, 0.89) = 282.26 day

Dsimu=weibull−10.7372203.80.89=282.26day

MDIsimutranformed=Φ−10.8421=1.0031

Applying the linear interpolation of F_MDIFMDI to [MDI, empirical distribution of MDI], we have the following:

MDI_simu = 1534.9 cfsMDIsimu=1534.9cfs. Finally, we have the corresponding simulation in the real domain as follows:

[D_simu, S_simu, MDI_simu] = [282.26 day, 8.19 × 10⁴cfs. day, 1534.9cfs]

DsimuSsimuMDIsimu=282.26day8.19×104cfs.day1534.9cfs.

Figure 13.13 Comparison of observed drought variables with simulations from the fitted vine copula.

Figure 13.14 compares the sample Kendall’s tau of the observed trivariate drought variables with those computed using the simulated trivariate variables from the vine copula. Comparison shows that the dependence structure of the drought variables is well preserved. The fitted vine copula may be applied further for risk analysis.

Figure 13.14 Comparison of sample Kendall’s tau and those simulated from the fitted vine copula.

Joint and Conditional Return Period through Vine Copula

In this section, we will proceed with risk analysis through the joint and conditional return period. In the case of joint return period, we will only investigate the “AND” case.

Joint Return Period “AND” Case

Similar to the bivariate case, the trivariate return period of the “AND” case can be given as follows:

TS>s∩D>d∩MDI>mdi=EINTPS>s∩D>d∩MDI>mdi=EINTCFS>FSs∩FD>FDd∩FMDI>FMDImdi(13.12)

where

CFS>FSs∩FD>FDd∩FMDI>FMDImdi=1−FSs−FDd−FMDImdi+CS,DGHFSsFDd+CS,MDIFrank(FSs,FMDImdi+CD,MDIFDdFMDImdi−CS,D,MDIGH−Frank−GaussianFSsFDdFMDImdi(13.12a)

To compute using Equation (13.12a), we will need to first evaluate C_{D, MDI}(F_D, F_MDI)CD,MDIFDFMDI. C_{D, MDI}(F_D, F_MDI)CD,MDIFDFMDI is the 2-margin of the trivariate copula as follows:

CD,MDIFDFMDI=CS,D,MDI1FDFMDI=∫01CS,D,MDIFDFMDItftdt(13.13)

where f(t) = 1ft=1 for the uniformly distributed random variable on [0, 1].

Figure 13.15 compares the empirical copula and that computed using Equation (13.13). In Figure 13.15, we also plot the joint CDF of D and MDI estimated separately using the Gumbel–Hougaard copula (θ = 2.22)θ=2.22). The goodness-of-fit study ensures the appropriateness with SnB = 0.05; P_value = 0.91SnB=0.05;Pvalue=0.91. From Figure 13.15, we see that there is minimal difference between the computed joint CDF (i.e., for D and DMI) using Equation (13.13) and that computed directly from the Gumbel–Hougaard copula (the maximum absolute difference is 0.028). To reduce the complexity of computation, we will directly apply the Gumbel–Hougaard copula to D and MDI. Setting MDI ≥ 1526.1 cfsMDI≥1526.1cfs (i.e., equivalent to F_MDI(mdi ≤ 1526.1) = 0.8FMDImdi≤1526.1=0.8), Figures 13.16 and 13.17 plot the joint exceedance probability (i.e., equivalent to survival copula) and the corresponding joint return period. Table 13.15 lists the sample results for computing the exceedance joint probability and its joint return period. As shown in Figures 13.16 and 13.17, the exceedance probability shows the concave shape and the joint return period shows the convex shape with MDI ≥ 1.53E + 03 cfsMDI≥1.53E+03cfs as the example graphed. It shows the low exceedance probability and high return period for a long duration with low drought severity, as well as for the short duration with high drought severity. The moderate drought duration and drought severity yield a shorter return period than does high drought severity with a short duration. We should note that with the change of plotted examples for MDI (or D or S), the shape may change accordingly.

Figure 13.15 Comparison of joint CDF computed for D and MDI.

Figure 13.16 Exceedance probability for C(D ≥ d, S ≥ s, MDI ≥ 1.53E + 03cfs)CD≥dS≥sMDI≥1.53E+03cfs.

Figure 13.17 Joint return period “AND” case for D ≥ d, S ≥ s, MDI ≥ 1.53E + 03D≥d,S≥s,MDI≥1.53E+03.

Table 13.15. Exceedance joint CDF and joint return period (“AND” case) with MDI ≥ 1.53E + 03 cfsMDI≥1.53E+03cfs.

Frequency domain			FDFD
			0.2	0.5	0.9	0.95	0.99
Exceed. joint CDF	FSFS	0.2	0.199	0.192	0.078	0.042	0.009
		0.5	0.199	0.192	0.078	0.042	0.009
		0.9	0.136	0.170	0.078	0.042	0.009
		0.95	0.109	0.156	0.077	0.041	0.009
		0.99	0.084	0.142	0.074	0.040	0.009

Real domain			D (days)a
			38	135	520	699	1133
Joint return period (yrs)	S	7860	3.626	3.774	9.281	17.402	83.468
	(cfs.day)	26880	3.633	3.774	9.281	17.402	83.468
		174920	5.298	4.265	9.292	17.402	83.468
		297460	6.620	4.634	9.424	17.424	83.468
		805320	8.595	5.078	9.727	17.875	83.546

^a The duration is rounded to the nearest integer number.

Conditional Return Period with the Constructed Vine Copula

Here we will investigate the following cases for the conditional return period as examples:

1. 
   

   (i) D>d ∩ MDI>mdi ∣ S ≤ sD>d∩MDI>mdi∣S≤s

   
   
2. 
   

   (ii) D>d ∩ MDI>mdi ∣ S = sD>d∩MDI>mdi∣S=s

   
   
3. 
   

   (iii) D>d ∪ MDI>mdi ∣ S ≤ sD>d∪MDI>mdi∣S≤s

   
   
4. 
   

   (iv) D>d ∪ MDI>mdi ∣ S = sD>d∪MDI>mdi∣S=s

   
   
5. 
   

   (iv) D>d ∣ MDI ≤ mdi ∩ S ≤ sD>d∣MDI≤mdi∩S≤s

   
   
6. 
   

   (vi) D>d ∣ MDI = mdi ∩ S = sD>d∣MDI=mdi∩S=s

Cases (i) and (ii): D>d ∩ MDI>mdi ∣ S ≤ SD>d∩MDI>mdi∣S≤S and D>d ∩ MDI>mdi ∣ S = sD>d∩MDI>mdi∣S=s

Cases (i) and (ii) investigate the impact of drought severity on drought duration and MDI during the drought episode under the condition of both D and MDI exceeding the corresponding critical level.

The conditional exceedance probability P(D>d ∩ MDI>mdi| S ≤ s)PD>d∩MDI>mdiS≤s for case (i) may be written as follows:

PD>d∩MDI>mdiS≤s=1−PD≤dS≤s−PMDI≤mdiS≤s+PD≤dMDI≤mdiS≤s=1−CDSFDFSFS−CMDI,SFMDIFSFS+CD,MDI,SFDFMDIFSFS(13.14)

and the corresponding conditional return period can simply be given as follows:

TD>d∩MDI>mdiS≤s=EINTPD>d∩MDI>mdiS≤s(13.14a)

The conditional exceedance probability P(D>d ∩ MDI>mdi| S = s)PD>d∩MDI>mdiS=s may be written as follows:

PD>d∩MDI>mdiS=s=1−PD≤dS=s−PMDI≤mdiS=s+PD≤dMDI≤mdiS=s(13.15)

where

PD≤dS=s=∂CDSFDFS∂FSFS=FSs;PMDI≤mdiS=s=∂CMDI,SFMDIFS∂FSFS=FSs(13.15a)

PD≤dMDI≤mdiS=s=CFDFMDIFS=FSs=∂CFDFMDIFS∂FSFS=FSs(13.15b)

Applying F_S = [0.2, 0.5, 0.9, 0.95, 0.99]FS=0.20.50.90.950.99 for the conditioning severity, we have the drought severity estimated as S≈[7860, 26880, 174920, 297460, 805320]S≈786026880174920297460805320cfs.day.

Table 13.16 lists and Figures 13.18 and 13.19 plot the conditional return period for cases (i) and (ii) using S = 174920 cfs.day as an example.

Table 13.16. Conditional exceedance probability and conditional return period: cases (i) and (ii).

		Case (i)				Case (ii)
		MDI>mdiMDI>mdi (cfs)				MDI>mdiMDI>mdi (cfs)
			646	1106	1537	646	1106	1537
Conditional exceedance prob.	D > d (day)	38	0.573	0.295	0.005	1.00	1.00	0.45
		135	0.295	0.131	0.001	0.99	0.99	0.44
		520	0.023	0.007	1.84E−05	0.41	0.41	0.12
		699	0.008	0.002	4.27E−06	0.20	0.20	0.043
		1133	0.001	2.38E−04	2.46E−07	0.037	0.037	0.005
Conditional return period	D > d (day)	38	1.26	2.45	160.04	0.72	0.72	1.62
		135	2.45	5.52	636.34	0.73	0.73	1.65
		520	31.24	103.24	3.94E+04	1.78	1.78	6.22
		699	85.11	316.69	1.69E+05	3.66	3.66	16.67
		1133	658.89	3038.20	2.93E+06	19.45	19.46	151.40

Figure 13.18 Conditional exceedance probability and conditional return period: case (i).

Figure 13.19 Conditional exceedance probability and conditional return period: case (ii).

Cases (iii) and (iv): D>d ∪ MDI>mdi ∣ S ≤ s; D>d ∪ MDI>mdi ∣ S = SD>d∪MDI>mdi∣S≤s;D>d∪MDI>mdi∣S=S

Cases (iii) and (iv) again investigate the impact of drought severity on drought duration and MDI but under different conditions, i.e., at least one drought variable (D or MDI) exceeding the corresponding critical level.

The conditional exceedance probability P(D>d ∪ MDI>mdi| S>s)PD>d∪MDI>mdiS>s of case (iii) may be rewritten with the following set of equations as follows:

PD>d∪MDI>mdiS≤s=1−PD≤d∩MDI≤mdiS≤s=1−CFDdFMDImdiFSsFSs(13.16)

The corresponding joint return period is then written as follows:

TD>d∪MDI>mdiS≤s=EINT1−PD>d∪MDI>mdiS≤s(13.16a)

The conditional exceedance probability P(D>d ∪ MDI>mdi| S = s)PD>d∪MDI>mdiS=s of case (iv) can be written as follows:

PD>d∪MDI>mdiS=s=1−PD≤d∩MDI≤mdiS=s=1−CFDdFMDImdiFSs(13.17)

Here it is obvious that C_{D, MDI ∣ S}CD,MDI∣S in Equation (13.17) is nothing but the copula function at T₂ of the vine structure. The corresponding conditional return period can be written as follows:

TD>d∪MDI>mdiS=s=EINT1−CFDdFMDImdiFSs(13.17a)

Figures 13.20 and 13.21 plot the conditional exceedance probability and the conditional return period for cases (iii) and (iv) using S = 174920 cfs.day as an illustrative sample. Table 13.17 lists the sample results.

Figure 13.20 Conditional exceedance probability and the corresponding conditional joint return period for case (iii).

Figure 13.21 Conditional exceedance probability and corresponding conditional joint return period for case (iv).

Table 13.17. Conditional exceedance probability and conditional return period: cases (iii) and (iv).

		Case (iii)				Case (iv)
		MDI>mdiMDI>mdi (cfs)				MDI>mdiMDI>mdi (cfs)
			646	1106	1537	646	1106	1537
Conditional exceedance prob.	D > d (day)	38	0.98	0.93	0.79	1.00	1.00	1.00
		135	0.93	0.76	0.46	1.00	1.00	1.00
		520	0.81	0.49	0.06	1.00	1.00	0.74
		699	0.79	0.46	0.03	1.00	1.00	0.60
		1,133	0.78	0.45	0.02	1.00	1.00	0.48
Conditional return period	D > d (day)	38	0.74	0.78	0.92	0.72	0.72	0.72
		135	0.78	0.95	1.58	0.72	0.72	0.72
		520	0.90	1.48	11.52	0.72	0.72	0.98
		699	0.91	1.56	21.46	0.72	0.72	1.20
		1,133	0.93	1.61	45.04	0.72	0.72	1.51

Cases (v) and (vi): D>d ∣ MDI ≤ mdi ∩ S ≤ sD>d∣MDI≤mdi∩S≤s and D>d ∣ MDI = mdi ∩ S = sD>d∣MDI=mdi∩S=s

Cases (v) and (vi) investigate the combined impact of maximum drought intensity and severity on drought duration.

For case (v), its conditional exceedance probability may be written as follows:

PD>dMDI≤mdi∩S≤s=PMDI≤mdiS≤s−PD≤dMDI≤mdiS≤sPMDI≤mdiS≤s=CMDI,SFRANKFMDIFS−CFDFSFMDICMDI,SFRANKFMDIFS(13.18)

and the corresponding conditional joint return period can be written as follows:

TD>dMDI≤mdi∩S≤s=EINTPD>dMDI≤mdi∩S≤s(13.18a)

For case (VI), its conditional exceedance probability can be written as follows:

PD>dMDI=mdi∩S=s=1−PD≤dMDI=mdi∩S=s=1−∂CMDI,D∣SCMDI,SFMDIFSsCD,SFDdFSs∂(CMDI,SFMDImdiFSs(13.19)

The corresponding conditional return period may be estimated through the vine copula as follows:

TD>dMDI=mdi∩S=s=EINTPD>dMDI=mdi∩S=s(13.19a)

The copula functions in Equation (13.19) directly reflect the constructed vine copula, i.e., C_{MDI, D ∣ S}CMDI,D∣S is Gaussian copula of T₂, and C_{D, S}CD,S and C_{MDI, S}CMDI,S are the Gumbel–Hougaard and Frank copula of T₁.

Table 13.18 lists the sample results using [S = 174920 cfs.day, MDI = 1526 cfs], [S = 25000 cfs.day, MDI = 800 cfs], and [S = 10000 cfs.day, MDI = 1000 cfs] as the illustration samples. Figures 13.22 and 13.23 plot the conditional exceedance probability and the corresponding conditional return period for cases (v) and (vi).

Table 13.18. Sample results for cases (v) and (vi).

		Case (v)				Case (vi)
		S ≤ s (cfs. day); MDI ≤ mdiS≤scfs.day;MDI≤mdi (cfs)				S = s; MDI = mdiS=s;MDI=mdi
			S = 174920; MDI = 1526	S = 25000; MDI = 800	S = 10000; MDI = 1000	S = 174920; MDI = 1526	S = 25000; MDI = 800	S = 10000; MDI = 1000
Conditional exceedance prob.	D > d (day)	38	0.81	0.68	0.28	1.00	1.00	0.60
		135	0.48	0.09	3.16E-03	1.00	0.58	0.002
		520	0.01	1.07E-06	2.68E-08	0.56	6.99E-06	3.80E-10
		699	2.24E-04	1.24E-08	3.09E-10	0.01	5.14E-08	1.18E-12
		1133	9.16E-09	5.04E-13	1.26E-14	1.40E-07	6.09E-13	0
Conditional return period	D > d (day)	38	0.89	1.06	2.59	0.72	0.72	1.21
		135	1.50	8.19	229.10	0.72	1.24	369.24
		520	51.20	6.75E+05	2.69E+07	1.29	1.03E+05	1.90E+09
		699	3234.94	5.85E+07	2.34E+09	63.30	1.41E+07	6.12E+11
		1133	7.89E+07	1.43E+12	5.74E+13	5.16E+06	1.19E+12	Inf

Figure 13.22 Conditional exceedance probability and conditional return period: case (v).

Figure 13.23 Conditional exceedance probability and conditional return period: case (vi).