Estimation of extreme Component-wise Excess design realization: a hydrological application

Estimation of extreme Component-wise Excess design realization: a hydrological application E. Di Bernardino∗and F. Palacios-Rodr´ıguez†‡ Abstract The classic univariate risk measure in environmental sciences is the Return Period (RP). The RP is traditionally defined as “the average time elapsing between two successive realizations of a prescribed event”. The notion of design quantile related with RP is also of great importance. The design quantile represents the “value of the variable(s) characterizing the event associated with a given RP”. Since an individual risk may strongly be affected by the degree of dependence amongst all risks, the need for the provision of multivariate design quantiles has gained ground. In contrast to the univariate case, the design quantile definition in the multivariate setting presents certain difficulties. In particular, Salvadori et al. (2011) define the design realization as the vector that maximizes a weight function given that the risk vector belongs to a given critical layer of its joint multivariate distribution function. In this paper, we provide the explicit expression of the aforementioned multivariate risk measure in the Archimedean copula setting. Furthermore, this measure is estimated by using Extreme Value Theory techniques and the asymptotic normality of the proposed estimator is studied. The performance of our estimator is evaluated on simulated data. We conclude with an application on a real hydrological data-set. Keywords: Multivariate design quantile, Extreme Value Theory, Return Period, Archimedean copula, hydrological application. 1 Introduction The Return Period (RP) is a frequently used concept in environmental sciences which allows us to identify dangerous events and find the means for rational decision-making and risk assessment. The RP is traditionally defined as “the average time elapsing between two successive realizations of a prescribed event” (Singh et al. (2007)). The notion of design quantile related with RP is also of great importance. The design quantile represents the “value of the variable(s) characterizing the event associated with a given RP” (Salvadori et al. (2013)). For instance, Ahmed et al. (2016) use univariate return periods to construct seasonal drought maps for different climatic seasons in Balochistan province (Pakistan). In addition, environmental risks could be involved by several variables that are often correlated (see e.g., Pappad`a et al. (2016)). For this reason, it is crucial to evaluate risks in a multivariate setting. In the univariate case, the design quantile is usually defined without ambiguity. However, in the multivariate setting, the sole tool provided by the RP may be insufficient to identify a design realization (see e.g., Serfling (2002), Vandenberghe et al. (2012)). ∗D´epartement IMATH, Laboratoire C´edric EA4629, CNAM, Paris, France. †Departamento de Estad´ıstica e Investigaci´on Operativa, Facultad de Matem´aticas, Universidad de Sevilla, Seville, Spain. ‡Correspondence to: F. Palacios-Rodr´ıguez, Departamento de Estad´ıstica e Investigaci´on Operativa, Facultad de Matem´aticas, Universidad de Sevilla, Calle Tarfia sin n´umero, Seville, 41012, Spain. Telephone: +34 95 4557927 - Fax: +34 95 4622800. E-mail: fp[email protected]. 1 Recently, several approaches have been proposed to deal with the problem of how to calculate the critical design event(s) in the multivariate case. Discussions concerning the meaning and use of return periods under non-stationary and multivariate conditions can be found in Serinaldi (2015b) and in Serinaldi and Kilsby (2015). Serinaldi (2015a) points out certain critical aspects that are often overlooked and that should be carefully taken into account for a correct interpretation of return periods. Di Bernardino et al. (2015) define the design quantile by using critical layers of the multivariate distribution function and of the multivariate survival function (see also the references therein). Torres et al. (2015) provide a definition of design quantile based on rotated directional distribution functions. On the other hand, Chebana and Ouarda (2011) provide a parametric estimator for critical layers and use it to study a real rainfall data-set. Multivariate frequency analysis is also used in spatial statistics in order to model the spatial variability of hydrological random variables (see e.g., De Paola and Ranucci (2012), De Paola et al. (2013)). Fawcett and Walshaw (2016) highlight the shortcomings of standard methods for estimating design quantiles and provide an estimation procedure which substantially increases the precision of design quantile estimates. Furthermore, Salvadori et al. (2011) estimate the Component-wise Excess design realization to study the behaviour of a dam by using a parametric Gumbel setting. In the present paper, we aim - to provide the explicit expression of the Component-wise Excess design realization given in Salvadori et al. (2011) by assuming a general Archimedean copula dependency between the risks considered; - to construct an extreme estimation procedure for the obtained Component-wise Excess design realization by using Extreme Value Theory (EVT); - to study and to compare the performance of our estimator on simulated data and with the estimator of Salvadori et al. (2011) on the same real dam data-set as studied in the aforementioned paper. The structure of the paper is as follows. The Component-wise Excess design realization given in Salvadori et al. (2011) is introduced in Section 2.1. In Section 2.2, we provide the analytical expression of the Component-wise Excess design realization in an Archimedean copula framework. The proposed extreme non-parametric estimator for the Component-wise Excess design realization is presented in Section 3. In addition, a Central Limit Theorem for this estimator is provided under certain regularly varying conditions. The performance of our estimator is analysed in Section 4. Finally, our estimator is compared with that of Salvadori et al. (2011) on the real dam data-set in Section 5. Proofs, notation, definitions, and auxiliary results are located in the Appendix. 2 Component-wise Excess design realization 2.1 Notation and preliminaries Let X= (X1, . . . , Xd) be a random vector with distribution function Fand survival distribution function F. In this multivariate setting, events that have equal probability of exceedance define iso-hyper-surfaces, otherwise known as critical layers. Salvadori et al. (2011) introduce the following definition: Definition 2.1 (Critical layers in Salvadori et al. (2011)).Let X= (X1, . . . , Xd)be a random risk vector with joint distribution function F. For α∈(0,1) and d≥2, the critical layer ∂L(α)at level αis defined as ∂L(α) = {x∈Rd:F(x) = α}. 2 The critical layer ∂L(α) partitions Rdinto three non-overlapping and exhaustive regions: ∂L<(α) = {x∈ Rd:F(x)< α},∂L(α) = the critical layer itself, and ∂L>(α) = {x∈Rd:F(x)> α}. The multivariate design quantile can be defined with respect to one of the three regions above. Furthermore, by taking into account the notion of critical layers in Definition 2.1, Salvadori et al. (2011) provide the following definition of design realization. Definition 2.2 (Design realization).Let ω:∂L(α)→[0,∞)be a weight function. The design realization δω∈∂L(α)is defined as δω(α) = arg max x∈∂L(α) ω(x) where ∂L(α) = {x∈Rd:F(x) = α}is the critical layer at level α, with α∈(0,1) and d≥2. Definition 2.2 is based on the idea of introducing a suitable function that “weights” the realizations lying on the critical layer of interest. Salvadori et al. (2011) point out that a realization lying on the critical layer may be of major interest when all of its marginal components are exceeded with the greatest probability. Furthermore, Salvadori et al. (2014) provide practical guidelines for coastal and off-shore engineering, where they illustrate how to calculate the design realization in Definition 2.2 by using two different weight functions. In the following, the multivariate survival distribution function Fof the considered risk vector is chosen as the weight function. Definition 2.3 (Constrained optimization problem for Component-wise Excess design realization).The Component-wise Excess design realization δCE of level αis defined as δCE(α) = arg max x∈∂L(α) P(X∈[x,∞)).(1) Definition 2.3 involves searching for the point(s) x= (x1, . . . , xd)∈∂L(α) in order to maximize the probability that a dangerous realization y= (y1, . . . , yd) verifies all the following component-wise inequalities: y1≥x1, . . . , yd≥xd. In the present paper, in order to model the dependence structure between the components of the random vector X, we assume an Archimedean copula model. Indeed, the Archimedean copulas constitutes one of the most useful classes of copulas in environmental science (see e.g., Pappad`a et al. (2016)). For instance, Saad et al. (2015) propose a multivariate flood-risk model based on nested Archimedean Frank copulas and Clayton copulas in a hydro-meteorological context in order to determine the 2011 Richelieu River flood-causing meteorological factors. Zhang et al. (2016) obtain useful information of the possible results of reservoir regulation on drought evolution under precipitation variations in two cascade reservoirs located in China by using the Gumbel copula. As will be proved below, by considering an Archimedean copula, the constrained optimization problem in Definition 2.3 becomes tractable (see Section 2.2). Let θbe a vector of dependence parameters. Recall that a d-dimensional Archimedean copula with generator ϕθand inverse generator ϕ−1 θis defined by C(v1, . . . , vd) = ϕ−1 θ(ϕθ(v1)+. . .+ϕθ(vd)), for all (v1, . . . , vd)∈[0,1]d. The generator ϕθis a continuous, convex and strictly decreasing function from [0,1] to [0,∞] such that ϕθ(1) = 0. If ϕθ(0) = +∞,ϕθis called a strict generator. For a thorough theoretical review of copulas, see Nelsen (2006), McNeil and Neˇslehov´a (2009), and Jaworski (2013). In the present paper, we will suppose the following assumption setting: Assumption 2.1 (Assumption setting). 3 a) The d−dimensional random risk vector Xfollows an Archimedean copula with a twice differentiable generator ϕθand a vector of dependence parameters θ. In addition, the d−dimensional random risk vector has partially strictly-increasing joint distribution function F. b) The marginal Xiconsidered has an absolutely continuous distribution Fi. Denote F−1 ithe leftcontinuous inverse of Fi, for i∈ {1, . . . , d}. c) The marginal Xiconsidered has a heavy tailed distribution, for i∈ {1, . . . , d}(see details in Section 3). d)x→ϕθ(1 −1/x)∈2RV−ρ,β(AY), with ρ∈[1,+∞) and β < 0. Condition a) in Assumption 2.1 is useful to obtain a more tractable constrained optimization problem for Component-wise Excess design realization in Definition 2.3 (see Proposition 2.1). Furthermore, condition c) is used in Section 3 to construct the extreme quantile estimator by using classic extrapolation techniques for heavy tailed distributions. Condition c) is also motivated by our applications in hydrology. Indeed, in these real-life applications, we can easily observe heavy tailed distributions (Pavlopoulos et al. (2008) and Papalexiou et al. (2013)). However, it is crucial to notice that our procedure can also be adapted in the case where Xifollows a light tailed distribution (de Haan and Ferreira (2006)). Notions of regularly varying (RV) and second-order regularly varying (2RV) function are recalled in Appendix A. 2.2 Expression of Component-wise Excess design realization under Archimedean copula assumption In this section, we solve the maximization problem in Equation (1) in order to obtain the closed-form expression of the Component-wise Excess design realization in our Archimedean copula framework. Proposition 2.1 (Component-wise Excess design realization in the Archimedean copula setting).Let X= (X1, . . . , Xd)the considered random risk vector. Assume that Xsatisfies conditions a) and b) in Assumption 2.1. The solution δCE(α)of the maximization problem in Equation (1) is given by δCE(α) = F−1 1ϕ−1 θϕθ(α) d, . . . , F−1 dϕ−1 θϕθ(α) d.(2) The proof of Proposition 2.1 is set out in Appendix C. Trivially, we remark that, in the degenerate univariate case (d= 1), the solution of the maximization problem in Equation (1) is given by the union of the univariate marginal quantiles. An illustration of Proposition 2.1 can now be presented. We consider the following bivariate models: •The Ali-Mikhail-Haq copula with dependence parameter θ= 0.5 and Uniform margins; •The Clayton copula with dependence parameter θ= 3 and Exponential margins with λ= 1; •The Gumbel copula with dependence parameter θ= 4 and Pareto margins with δ1= 1 and δ2= 2. In Figure 1, we show critical layers ∂L(α) (red lines), the survival functions F(x1, x2) (black surfaces), δCE(α) (green stars) for different values of αand the associated maximum points (blue dots) for the three models considered. 4 x1 x2 Survival distribution Ali−Mikhail−Haq Copula with theta= 0.5 and Uniform margins x1 x2 Survival distribution Clayton copula with theta=3 and Exponential margins x1 x2 Survival distribution Gumbel copula with theta=4 and Pareto margins Figure 1: Critical layers ∂L(α) (red lines), survival functions (black surfaces), δCE(α) (green stars) and the associated maximum points (blue dots) for the models considered. We take α= 0.1 and 0.5 (left panel); α= 0.6 and 0.8 (centre panel); α= 0.8 and 0.9 (right panel). 3 Proposed Component-wise Excess design realization estimator and consistency result Assume that the considered risk vector Xsatisfies conditions in Assumption 2.1. 3.1 Estimation procedure Using Proposition 2.1, we now propose an extreme non-parametric estimation procedure for each component of the Component-wise Excess design realization δCE(α) in Equation (2) for extreme value α:= αn→1, for n→ ∞, where nis the sample size considered and d≥2. From conditions a) and b) in Assumption 2.1, the i-th component of δCE(α) in Equation (2) can be written as the (1 −p)-quantile of the random variable Xiwith p= 1 −ϕ−1 θϕθ(α) d, for i∈ {1, . . . , d}. Let α:= αn→1, therefore p:= pn→0, as n→+∞. The goal of this section is to estimate xi pn:= UXi1 pn,for i∈ {1, . . . , d},(3) for pn:= 1 −ϕ−1 θϕθ(αn) d→0, as n→+∞, and UXi(t) := F−1 i1−1 t, for t > 1. Therefore, the final estimator is based on the following “plug-in procedure”: we first deal with the estimation of the random level pn, and secondly with the estimation of xi pnin (3). These two steps are detailed below. First step: estimation of the random risk level pnLet Vi=Fi(Xi), for i∈ {1, . . . , d}. The vector V:= (V1, . . . , Vd) has the considered Archimedean copula as its distribution function, that is, P[max{V1, . . . , Vd} ≤ v] = ϕ−1 θ(d ϕθ(v)) , v ∈(0,1). 5 Let Yd := max{V1, . . . , Vd}, therefore pn= 1 −F−1 Y(αn). The random risk level pncan be written as a function of the αn-quantile of Y, that is, pn= 1 −UY1 1−αn= 1 −xY 1−αn. Under conditions a), b) and d) in Assumption 2.1, from Proposition B.1 in Appendix B, Y∈MDA (Ψ−1). Therefore, by using Theorem 1.1.13 in de Haan and Ferreira (2006), it is verified that UY∈ERV−1(see Section B.3 in de Haan and Ferreira (2006) for the definition of ERV). Following the same approximation technique as in Equation (3.1.6) in de Haan and Ferreira (2006), we obtain xY 1−αn=UY1 1−αn∼UYn kY+aYn kYkY n(1−αn)γY−1 γY ,(4) where kYis an intermediate sequence such that kY=kY(n)→ ∞,kY(n)/n →0, as n→ ∞. Let Yn−kY,n be the (n−kY)-th order statistic of the sample (Y1, . . . , Yn). The natural estimator of UYn kY is its empirical counterpart Yn−kY,n. In order to estimate aY(n/kY) and γY, we consider the probabilityweighted moment estimators defined in Equations (3.6.9) and (3.6.10) in de Haan and Ferreira (2006). For the sake of clarity, these estimators are laid out below. We take bγY= 1 −Pn 2Qn−1−1 (5) and baYn kY:= bσPWM =PnPn 2Qn−1−1 ,(6) with Pn:= 1 kYPkY−1 i=0 Yn−i,n −Yn−kY,n and Qn:= 1 kYPkY−1 i=0 i kY(Yn−i,n −Yn−kY,n). The consistency results for the estimators Yn−kY,n,bγYand baY(n/kY) can be found in Theorems 2.4.1 and 3.6.1 in de Haan and Ferreira (2006). Hence, from Equation (4), we obtain bpn= 1 −bxY 1−αn= 1 −Yn−kY,n +baYn kYkY n(1−αn)bγY−1 bγY.(7) Second step: estimation of the extreme quantile xi pnFrom condition c) in Assumption 2.1, there exists γi>0 such that for all x > 0, UXi∈RVγi(+∞) (Appendix A). In this case, as n→ ∞, xi pn=UXi1 pn∼UXin ki ki n pnγi ,(8) where kiis an intermediate sequence such that ki=ki(n)→ ∞,ki(n)/n →0, as n→ ∞. We denote the (n−ki)-th order statistic of the sample (Xi 1, . . . , Xi n) as Xi n−ki,n. We estimate the tail index γiin 6 Equation (8) by using the Hill estimator (Hill (1975)): bγi=1 ki ki−1 X j=0 log Xi n−j,n −log Xi n−ki,n.(9) We can therefore estimate xi pnin (8) by bxi pn=Xi n−ki,n ki n pnbγi .(10) Finally, by using a plug-in technique with Equations (7) and (10), we propose the following extreme estimator for the i-th component of δCE(α) in Equation (2): bxi bpn=Xi n−ki,n ki nbpnbγi .(11) Remark 3.1. Notice that Yis an unobservable random variable in the real-life applications (see, e.g., Section 5). In this case, one has previously to construct a pseudo-observed version of Y. More precisely, consider a real-data set of observations {X(k)= (X(k) 1, . . . , X(k) d)}k∈{1,...,n}of the random risk vector X. Firstly, we define pseudo-observations {V(k)= (V(k) 1, . . . , V (k) d)}k∈{1,...,n}by setting every component i for observation number kto V(k) i=1 n+ 1 n X j=1 1nX(j) i≤X(k) io, with i= 1, . . . , d,k∈ {1, . . . , n}. One can check that for any i= 1, . . . , d,k∈ {1, . . . , n},V(k) i∈(0,1). Therefore, we get the desired univariate sample {b Yk:= max{V(k) 1, . . . , V (k) d}}k∈{1,...,n}. Consistency of this method in a non-parametric setting has been established in Einmahl et al. (2001), Einmahl and Segers (2009). We remark that, since Ycan not be observed, in our procedure we neglect the uncertainty induced by the margins. More precisely, since the uncertainty induced by the estimation of b Ykis not taken into account in our procedure, then the main result in the following Theorem 3.1 is only valid under full knowledge of the margins. Then, the confidence intervals provided in Theorem 3.1 can be only understood as optimistic approximations. However, this issue is beyond the scope of the present paper and may be addressed in future work. 3.2 Consistency result In this section, we obtain the Central Limit Theorem for the proposed extreme estimator in (11). Theorem 3.1 provides two results: the consistency for the estimated random risk level pn(see Equation (12)); the consistency of the final extreme quantile estimator bxi bpn(see Equation (14)). Theorem 3.1. Assume that X= (X1, . . . , Xd)satisfies conditions in Assumption 2.1. Let Yd = max{V1, . . . , Vd}with Vi d =Fi(Xi), for i∈ {1, . . . , d}. Let αn→1,pn= 1 −xY 1−αnand bpnremain as in Equation (7). Suppose: 1. kY=kY(n)→ ∞,kY/n →0and √kYAY(n/kY)→λ∈R, for n→ ∞; 7 2. n(1 −αn) = o(kY)and log (n(1 −αn)) = o(√kY), for n→ ∞. Let dY n:= kY/(n(1−αn)) and qγ(t) := Rt 1sγ−1log (s)ds. Let vn=√kY aYn kYqγY(dY n)and bvn=√kY baYn kYqbγY(dY n). Then, for n→ ∞,bvn(bpn−pn)d →Θ1(12) where Θ1= Γ + B+Λ+ λ β−1with Ba standard normal distribution and Γand Λnormal distributions defined as in Theorems 3.6.1 and 4.3.1 in de Haan and Ferreira (2006). Moreover, assume that: 3. UXi∈2RVγi,τi(Ai),γi>0and τi<0; 4. ki=ki(n)→ ∞,ki/n →0and √kiAi(n/ki)→λi∈R, for n→ ∞; 5. npn=o(ki)and log (n pn) = o(√ki), for n→ ∞. Let di n:= ki/(n pn). Define xi pnand bxi bpnas in Equations (8) and (11) respectively. If, for n→ ∞, √ki log (di n)vn→0,(13) then it is verified that √ki log (di n) bxi bpn xi pn−1!d →Θ2,(14) where Θ2is a normal random variable with mean λi/(1 −τi)and variance γ2 i. For the sake of brevity, the proof of Theorem 3.1, together with several necessary auxiliary results are set out in Appendices B and C (see e.g., Proposition B.2). It should be borne in mind that the quantity di nin the convergence rate of Equation (14) remains unknown. However, an adaptive version of this consistency result is provided below. Corollary 3.1. Let b di n:= ki/(nbpn). Conditions of Theorem 3.1 imply log ( b di n) log (di n) P →1, for n→ ∞. Hence an equivalent statement of Equation (14) is √ki log (b di n) bxi bpn xi pn−1!d →Θ2, where Θ2is defined in Theorem 3.1. Proof of Corollary 3.1 is given trivially by using the convergency bpn pn P →1 for n→ ∞(see proof of Theorem 3.1 in Appendix C). Note that the latter form proposed in Corollary 3.1 is more useful for constructing confidence intervals for bxi bpn(see Section 5). 4 Simulation study A simulation and comparison study is implemented to investigate the finite sample performance of our estimator in this section. 8 The estimation procedure presented in this paper involves different notation progressively introduced in Section 3. To improve clarity, a comprehensive scheme of our extreme estimation procedure is presented in Algorithm 1. The nature of the considered parameters is specified; in particular, we distinguish between tuning parameters and estimated/calculated quantities. Firstly, we simulate under the chosen Archimedean copula model. Simulations Choose the sample size n. Choose the dimension d≥2. Simulate a sample (uk 1, . . . , uk d) from a d-Archimedean copula, with k∈ {1, . . . , n}. Build the sample Yk:= max{uk 1, . . . , uk d}, for k∈ {1, . . . , n}. Get marginal distributions sample xk i=F−1 i(uk i) for i= 1, . . . , d and k= 1, . . . , n. Secondly, we follow Algorithm 1 below to obtain our estimator. Algorithm 1 Comprehensive scheme for the estimator bxi bpn Input parameters Choose the extreme risk level α=αn. Select the margin Xi, for i∈ {1, . . . , d}. Estimations Choose kY(n) and estimate - the tail index bγYby using the estimator in Equation (5); - the scale sequence baYn kYby using the estimator in Equation (6); - the intermediate order statistic Yn−kY,n. A specific value for kY(n)is chosen by following the stability strategy proposed in Remark 4.1. Obtain bpnas in Equation (7) by using quantities bγY,baYn kYand Yn−kY,n. Choose ki(n) and estimate - the tail index ˆγiby using the Hill estimator in Equation (9); - the intermediate order statistic Xi n−ki,n. A specific value for ki(n)is chosen by following the stability strategy proposed in Remark 4.1. Obtain bxi bpnas in Equation (11) by using quantities bpn, ˆγiand Xi n−ki,n. We now consider the following models in dimensions d= 2 and d= 5: i) The Gumbel copula with dependence parameter θ= 3 and Fr´echet margins with β= 3 (i.e., Fi(t) = exp{−t−β}). ii) The Joe copula with dependence parameter θ= 2 and Pareto margins with δ1= 1 and δ2= 2 (i.e., Fi(t)=1−(δ1/(t+δ1))δ2). It should be borne in mind that assumptions of Theorem 3.1 are verified in the considered cases i) and ii). 9 points that have the greatest probability of being component-wise exceeded by an extreme realization with return period larger than 1000 years. Therefore, these points could be interpreted as a “safety lowerbound”. That is, the structure under design should, at least, support realizations having multivariate size the Component-wise Excess design realization δCE(α) for the millinery return period α≈0.946537. Hence, the underestimation of these quantities can represent a major risk for dam managers and for environmental practitioners. Furthermore, the target of the dam manager is to maintain a high water level, in order to achieve the maximum benefit concerning the production of electrical energy. From the pair (Q, V ) of the calculated Component-wise Excess design realization, it is possible to obtain the associated flood hydrograph with peak Qand volume V. By using the flood hydrograph, one can calculate the maximum level of the dam associated with (L, Q, V ), and one can check whether or not the crest level of the dam is exceeded by the reservoir level (see Salvadori et al. (2011)). The maximum water level obtained in Table 2 in Salvadori et al. (2011) for the associated values of (Q, V ) in δS CE (see the first row in Table 1) is 782.08 mabove sea-level. Since the values obtained by using bxbpnare very similar to those of Salvadori et al. (2011), we can compare our Lrealization with the maximum water level of 782.08 mabove sea-level. In this work, we focus on the estimation of the multivariate quantile of the hydrological load acting on the structure, i.e., the spillway of the dam. However, in order to evaluate the safety of the dam, one has to consider the interaction between the hydrological load and the structure. Volpi and Fiori (2014) point out the importance of considering the structure in hydraulic design and/or risk assessment problems in a multivariate environment and advise against the uncritical use of design event-based approaches. Indeed, the relationship between the structure and the hydrological loads acting on it is neglected in the study of the present paper. Requena et al. (2016) and references herein highlight the importance of considering the specific structure when designing or assessing flood risks in a multivariate context. Salvadori et al. (2016) illustrate the structural approach in a real sea storms data-set. Appendix A Definitions of regular varying functions In the following, the notions of regularly varying (RV) and second-order regularly varying (2RV) function are introduced. These definitions are useful in Section 3. For more details, the reader is referred to Section B.1 in de Haan and Ferreira (2006) and Mao and Hu (2012). Definition A.1 (RV function).A measurable function, h:R+→Rthat is eventually positive, is said to be of regular variation at infinity with index γ∈R\{0}, denoted by h∈RVγ(+∞), if for any x > 0, lim t→∞ h(tx) h(t)=xγ.(A.1) If (A.1) holds with γ= 0 for any x > 0, then his said to be slowly varying at infinity and is written as h∈RV0(+∞). Definition A.2 (2RV function).A measurable function, h:R+→Rthat is eventually positive, is said to be of second-order regular variation with the first-order parameter γ∈Rand the second-order parameter 16 τ≤0, denoted by h∈2RVγ,τ (A), if there exists some ultimately positive or negative function A(t)with A(t)→0as t→ ∞ such that lim t→∞ h(tx) h(t)−xγ A(t)=H(x),with H(x) = xγZx 1 sτ−1ds, ∀x > 0. Here, Ais referred to as an auxiliary function of h. It is verified that |A| ∈ RVτ(+∞). B First- and second-order tail behaviour of Yd = max{V1, . . . , Vd} Let Yd = max{V1, . . . , Vd}be the random variable introduced in Section 3. The (first-order) tail behaviour of Yis studied in the following proposition. The interested reader is also referred to Theorem 1.1.8 in de Haan and Ferreira (2006). Proposition B.1 (The von Mises condition for Y).Assume that X= (X1, . . . , Xd)satisfies conditions a) and b) in Assumption 2.1. Assume that ϕθ∈RVρ(1), with ρ∈[1,+∞). Let Yd = max{V1, . . . , Vd} with Vi d =Fi(Xi), for i∈ {1, . . . , d}. Then, FYverifies the von Mises condition with tail index γY=−1. In particular, Y∈MDA(ΨγY=−1). Proof. From Section 3, we obtain F′ Y(t)=(ϕ−1 θ)′[d ϕθ(t)] d ϕ′ θ(t) and F′′ Y(t)=(ϕ−1 θ)′′[d ϕθ(t)] d2(ϕ′ θ(t))2+ (ϕ−1 θ)′[d ϕθ(t)] d ϕ′′ θ(t). The limit in Equation (1.1.30) in de Haan and Ferreira (2006) can now be calculated. lim t↑1 (1 −FY(t))F′′ Y(t) (F′ Y(t))2= lim t↑1 (1 −ϕ−1 θ(d ϕθ(t))) (ϕ−1 θ)′′[d ϕθ(t)] d2(ϕ′ θ(t))2+ (ϕ−1 θ)′[d ϕθ(t)] d ϕ′′ θ(t) (ϕ−1 θ)′[d ϕθ(t)] d ϕ′ θ(t)2 = lim t↑1(1 −ϕ−1 θ(d ϕθ(t))) (ϕ−1 θ)′′[d ϕθ(t)] (ϕ′ θ(ϕ−1 θ(d ϕθ(t))))2 + lim t↑1(1 −ϕ−1 θ(d ϕθ(t))) ϕ′ θ(ϕ−1 θ(d ϕθ(t))) ϕ′′ θ(t) d(ϕ′ θ(t))2. Given the assumptions ϕθ∈RVρ(1), then ϕ′ θ∈RVρ−1(1). Therefore, lim t↑1(1 −ϕ−1 θ(d ϕθ(t))) (ϕ−1 θ)′′[d ϕθ(t)] (ϕ′ θ(ϕ−1 θ(d ϕθ(t))))2=ρ−1 and lim t↑1(1 −ϕ−1 θ(d ϕθ(t))) ϕ′ θ(ϕ−1 θ(d ϕθ(t))) ϕ′′ θ(t) d(ϕ′ θ(t))2=−ρ+ 1. Finally, we conclude limt↑1(1−FY(t))F′′ Y(t) (F′ Y(t))2= 0 ⇒γY=−1.Hence the result follows. The (second-order) tail behaviour of Yis studied in the following proposition. 17 Proposition B.2. Assume that X= (X1, . . . , Xd)satisfies conditions a), b) and d) in Assumption 2.1. Let Yd = max{V1, . . . , Vd}with Vi d =Fi(Xi), for i∈ {1, . . . , d}. Then, 1−UY∈2RV−1,β(AY). Proof. From Section 3, it is known that FY(1 −1/t)=1−ϕ−1 θ(d ϕθ(1 −1/t)) . Since x→ϕθ(1 −1/x)∈2RV−ρ,β(AY), from Proposition 2.4 in Mao and Hu (2012), we can write FY(1 −1/t)=1−ϕ−1 θhd c t−ρ1 + 1 βAY(t) + o(AY(t))i. In addition, x→ϕθ(1 −1/x)∈2RV−ρ,β(AY) implies ϕθ∈RVρ(1). From Remark C in Di Bernardino and Rulli`ere (2014), it is verified that 1 −ϕ−1 θ(1/x)∈RV−1/ρ(+∞). We now obtain FY(1 −1/t)=(dc)1/ρ t−11 + 1 βAY(t) + o(AY(t))1/ρ . By using Taylor expansion, FY(1 −1/t)=(dc)1/ρ t−11 + 1 β ρ AY(t) + o(AY(t)). It can be observed that |˜ AY|:= |1 ρAY| ∈ RVβsince assumptions. From Proposition 2.4 in Mao and Hu (2012), the result is given. C Proofs Proof of Proposition 2.1 Let α∈(0,1). Let Cand Cbe, respectively, the copula and the joint survival function associated to the random vector V= (V1, . . . , Vd). Note that the copula version of the constrained optimization problem in Equation (1) can be written as δCE(α) = arg max v∈∂LC(α) C(v1, . . . , vd),(C.1) with ∂LC(α) = {v∈[0,1]d:C(v) = α}. Equivalently, one can write the constrained optimization problem (C.1) as arg max s C(ϕ−1 θ(s1), . . . , ϕ−1 θ(sd)) s.t.: d X i=1 si=ϕθ(α), si≥0,for i= 1, . . . , d. (C.2) 18 Our aim is to find the solution of the constrained optimization problem in (C.2). From Theorem 2.21 in Boche and Jorswieck (2007), if C(ϕ−1 θ(s1), . . . , ϕ−1 θ(sd)) is a schur-concave function1, the global maximum for the problem in (C.2) is achieved by s∗=ϕθ(α) d,...,ϕθ(α) d. We now study the schur-concavity of the function C(ϕ−1 θ(s1), . . . , ϕ−1 θ(sd)). We now prove that the d-dimensional survival copula Cassociated to an Archimedean copula is a schurconcave function. To this end, it is helpful to realize that by using the symmetry property, one can take d= 2 without loss of generality. The interested reader is referred to Marshall et al. (2011) (Section A.5) for further details. That is, it is sufficient to prove that bivariate survival copula Cis a schur-concave function. In addition, every Archimedean copula is schur-concave (see Lemma 10.2.2 in Durante (2006) and Proposition 4.11 in Dolati and Dehgan Nezhad (2014)). Furthermore, from Proposition 10.1.7 in Durante (2006), a bivariate copula Cis schur-concave if and only if ˆ C(u, v)2associated to Cis also a schur-concave function. Since, in our case Cis an Archimedean copula (see condition a) in Assumption 2.1) then ˆ C(u, v) is also a schur-concave function. Remark that C(u, v) = ˆ C(1 −u, 1−v), for (u, v)∈[0,1]2(see Nelsen (2006)). Then to obtain the desired result, we have to prove that ˆ C(f(u), f(v))), with f(u) = 1 −u,u∈[0,1], is a schur-concave function. This last statement holds true because ˆ C(f(u), f(v)) is a composition of a an increasing schur-concave function ( ˆ C) and a concave function (f). The interested reader is referred to Marshall et al. (2011) (Section B.2) for further details. Using similar arguments, since Cis a decreasing schur-concave function and ϕ−1 θa convex function, C(ϕ−1 θ(s1), . . . , ϕ−1 θ(sd)) is a schur-concave function. Finally, by taking si=ϕ(vi), for i= 1, . . . , d, from Theorem 2.21 in Boche and Jorswieck (2007), the global maximum in problem (C.2) is achieved by v∗=ϕ−1 θϕθ(α) d, . . . , ϕ−1 θϕθ(α) d. By using the Probability Integral Transform Theorem (see Section 1.5.8.3 in Denuit et al. (2005)) for each margin, we obtain the result. More precisely, the global optimum point for the constrained optimization problem in Equation (1) is given by F−1 1(v∗ 1), . . . , F−1 d(v∗ d).□ Proof of Theorem 3.1 Firstly, note that bvn(bxY 1−αn−xY 1−αn)d →Θ1,(C.3) where Θ1= Γ + B+Λ+ λ β−1with Ba standard normal and, Γ and Λ normal distributions as defined in Theorems 3.6.1 and 4.3.1 in de Haan and Ferreira (2006). Indeed, from Proposition B.2 in Appendix B, condition d) in Assumption 2.1 implies that 1 −UY∈ 2ERV−1,β(AY) and a(t) = (1 −UY(t)). Therefore, the asymptotic result in Equation (C.3) results from Theorems 3.6.1 and 4.3.1, and Corollary 4.3.2 in de Haan and Ferreira (2006). Consequently, we obtain bvn(bpn−pn)d → −Θ1, (Theorem on page 24 in Serfling (1980)). 1A real valued function g:A⊆Rn→R, is schur-concave (resp. schur-convex) on Aif for all a,b∈A,a≺bimplies g(a)≥g(b) (resp. g(a)≤g(b)), where ≺represents the majorization order (see Definition A.1 in Marshall et al. (2011). 2Let (X, Y ) be a bivariate vector with margin survival distributions FXand FY. Assume that (X, Y ) follows an Archimedean copula C. Then, ˆ C: [0,1]2→[0,1] is defined as ˆ C(u, v) = u+v−1 + C(1 −u, 1−v). Furthermore, P[X≥x, Y ≥y] = ˆ C(FX(x), F Y(y)). For more details see Section 2.6 in Nelsen (2006). 19 Under conditions 3, 4 and 5 in Theorem 3.1, and by applying Theorem 4.3.8 in de Haan and Ferreira (2006), we determine that √ki log (di n) bxi pn xi pn−1!d →Θ2, where Θ2is a normal random variable with mean λi/(1 −τi) and variance γ2 i(see Theorem 3.2.5 in de Haan and Ferreira (2006)). We now write bxi bpnas a function of bxi pn. That is, we can write bxi bpn=Xi n−ki,n ki nbpn pnpn!bγi =bxi pnbpn pn−bγi . Therefore, we obtain √ki log (di n) bxi bpn xi pn−1!=√ki log (di n)bpn pn−bγi bxi pn xi pn−1!+√ki log (di n)"bpn pn−bγi −1#.(C.4) On the other hand, bpn pn =−Θ1 vn +1+oP1 vn. Hence, we obtain √ki log (di n)"bpn pn−bγi −1#=√ki log di n"−Θ1 vn +1+oP1 vn−bγi −1#. By using Taylor expansion and condition (13), it is verified that √ki log (di n)"bpn pn−bγi −1#P →0.(C.5) In addition, Equation (C.5) implies that bpn pn−bγiP →1.(C.6) Finally, by using (C.5) and (C.6) in Equation (C.4), from Slutsky’s Theorem, we attain the result. □ Acknowledgements The authors thank the associated editor and the referees whose comments helped to improve a previous version of this paper. Furthermore, the authors thank Gianfausto Salvadori and Fabrizio Durante for fruitful discussions. This work was partly supported by a grant from the Junta de Andaluc´ıa (Spain) for research group (FQM- 328) and by a pre-doctoral contract (Palacios Rodr´ıguez, F.) from the “V Plan Propio de Investigaci´on” of the University of Seville. 20 References Ahmed, K., Shahid, S., bin Harun, S., and Wang, X. J. (2016). Characterization of seasonal droughts in Balochistan Province, Pakistan. Stochastic Environmental Research and Risk Assessment, 30(2):747–762. Boche, H. and Jorswieck, E. A. (2007). Majorization and Matrix-monotone Functions in Wireless Communications. Now Publishers Inc. Delft. Cai, J. J., Einmahl, J. H. J., de Haan, L., and Zhou, C. (2015). Estimation of the marginal expected shortfall: the mean when a related variable is extreme. Journals of the Royal Statistical Society. Series B. Statistical Methodology, 77(2):417—-442. Chebana, F. and Ouarda, T. B. M. J. (2011). Multivariate quantiles in hydrological frequency analysis. Environmetrics, 22:63–78. de Haan, L. and Ferreira, A. (2006). Extreme Value Theory. An Introduction. Springer Series in Operations Research and Financial Engineering. Springer: New York. De Michele, C., Salvadori, G., Canossi, M., Petaccia, A., and Rosso, R. (2005). Bivariate Statistical Approach to Check Adequacy of Dam Spillway. Journal of Hydrologic Engineering, 10(1):50–57. De Paola, F. and Ranucci, A. (2012). Analysis of spatial variability for stormwater capture tanks assessment. Irrigation and Drainage, 61(5):682–690. De Paola, F., Ranucci, A., and Feo, A. (2013). ANTECEDENT MOISTURE CONDITION (SCS) FREQUENCY ASSESSMENT: A CASE STUDY IN SOUTHERN ITALY. Irrigation and Drainage, 62:61–71. Denuit, M., Dhaene, J., Goovaerts, M., and Kaas, R. (2005). Actuarial Theory for Dependence Risks: Measures, Orders and Models. Wiley. Di Bernardino, E., Fern´andez-Ponce, J., Palacios-Rodr´ıguez, F., and Rodr´ıguez-Gri˜nolo, M. (2015). On Multivariate Extensions of the Conditional Value-at-Risk measure. Insurance: Mathematics and Economics, 61:1—-16. Di Bernardino, E. and Rulli`ere, D. (2014). On tail dependence coefficients of transformed multivariate Archimedean copulas. Working paper, https://hal.archives-ouvertes.fr/hal-00992707v1. Dolati, A. and Dehgan Nezhad, A. (2014). Some Results on Convexity and Concavity of Multivariate Copulas. Iranian Journal of Mathematical Sciences and Informatics, 9(2):87–100. Durante, F. (2006). New results on copulas and related concepts. PhD thesis. Universit`a degli Studi di Lecce. Italy. Durante, F. and Okhrin, O. (2015). Estimation procedures for exchangeable Marshall copulas with hydrological application. Stochastic Environmental Research and Risk Assessment, 29(1):205–226. Einmahl, J., De Haan, L., and Piterbarg, V. (2001). Nonparametric Estimation of the Spectral Measure of an Extreme Value Distribution. The Annals of Statistics, 29(5):1401–1423. Einmahl, J. and Segers, J. (2009). Maximum Empirical Likelihood Estimation of the Spectral Measure of an Extreme-Value Distribution. The Annals of Statistics, 37(5B):2953–2989. Fawcett, L. and Walshaw, D. (2016). Sea-surge and wind speed extremes: optimal estimation strategies for planners and engineers. Stochastic Environmental Research and Risk Assessment, 30(2):463–480. Hill, B. M. (1975). A simple general approach to inference about the tail of a distribution. The Annals of Statistics, 3(5):1163—- 1174. Jaworski, P. (2013). The Limiting Properties of Copulas Under Univariate Conditioning. In Copulae in mathematical and quantitative finance, volume 213, pages 129–163. Lecture Notes in Statistics, Springer, Heidelberg. Mao, T. and Hu, T. (2012). Second-order properties of the Haezendonck-Goovaerts risk measure for extreme risks. Insurance: Mathematics and Economics, 51:333–343. Marshall, A., Olkin, I., and Arnold, B. C. (2011). Inequalities: Theory of Majorization and Its Applications. Second Edition. Springer, New York. McNeil, A. and Neˇslehov´a, J. (2009). Multivariate Archimedean copulas, d-monotone functions and l1−norm symmetric distributions. The Annals of Statistics, 37(5B):3059–3097. Nelsen, R. B. (2006). An Introduction to Copulas. Springer Series in Statistics. Springer: New York. Papalexiou, S. M., Koutsoyiannis, D., and Makropoulos, C. (2013). How extreme is extreme? an assessment of daily rainfall distribution tails. Hydrology and Earth System Sciences, 17:851–862. 21 Pappad`a, R., Perrone, E., Durante, F., and Salvadori, G. (2016). Spin-off Extreme Value and Archimedean copulas for estimating the bivariate structural risk. Stochastic Environmental Research and Risk Assessment, 30(1):327–342. Pavlopoulos, H., Picek, J., and Jureˇckov´a, J. (2008). Heavy tailed durations of regional rainfall. Applications of Mathematics, 53(3):249–265. Requena, A. I., Chebana, F., and Mediero, L. (2016). A complete procedure for multivariate index-flood model application. Journal of Hydrology, 535:559–580. Saad, C., El Adlouni, S., St-Hilaire, A., and Gachon, P. (2015). A nested multivariate copula approach to hydrometeorological simulations of spring floods: the case of the Richelieu River (Qu´ebec, Canada) record flood. Stochastic Environmental Research and Risk Assessment, 29(1):275–294. Salvadori, G., De Michele, C., and Durante, F. (2011). On the return period and design in a multivariate framework. Hydrology and Earth System Sciences, 15:3293–3305. Salvadori, G., Durante, F., De Michele, C., Bernardi, M., and Petrella, L. (2016). A multivariate copula-based framework for dealing with hazard scenarios and failure probabilities. Water Resources Research, 52(5):3701–3721. Salvadori, G., Durante, F., and Perrone, E. (2013). Semi-parametric approximation of Kendall’s distribution function and multivariate Return Periods. Journal de la Soci´et´e Fran¸caise de Statistique, 154(1):151—-173. Salvadori, G., Tomasicchio, G. R., and D’Alessandro, F. (2014). Practical guidelines for multivariate analysis and design in coastal and off-shore engineering. Coastal Engineering, 88:1—-14. Serfling, R. (1980). Approximation Theorems of Mathematical Statistics. John Wiley & Sons, New York. Serfling, R. (2002). Quantile functions for multivariate analysis: approaches and applications. Statistica Neerlandica, 56 (2):214— -232. Serinaldi, F. (2015a). Can we tell more than we can know? The limits of bivariate drought analyses in the United States. Stochastic Environmental Research and Risk Assessment. DOI 10.1007/s00477-015-1124-3. Serinaldi, F. (2015b). Dismissing return periods!. Stochastic Environmental Research and Risk Assessment, 29(4):1179–1189. Serinaldi, F. and Kilsby, C. G. (2015). Stationarity is undead: Uncertainty dominates the distribution of extremes. Advances in Water Resources, 77:17–36. Singh, V., Jain, S., and Tyagi, A. (2007). Risk and Reliability Analysis. ASCE Press, Reston, Virginia. Torres, R., Lillo, R. E., and Laniado, H. (2015). A directional multivariate value at risk. Insurance: Mathematics and Economics, 65:111–123. Vandenberghe, S., van den Berg, M. J., Gr¨aler, B., Petroselli, A., Grimaldi, S., De Baets, B., and Verhoest, N. E. C. (2012). Joint return periods in hydrology: a critical and practical review focusing on synthetic design hydrograph estimation. Hydrology and Earth System Sciences Discussions, 9:6781—-6828. Volpi, E. and Fiori, A. (2014). Hydraulic structures subject to bivariate hydrological loads: Return period, design, and risk assessment. Water Resources Research, 50(2):885–897. Zhang, R., Chen, X., Cheng, Q., Zhang, Z., and Shi, P. (2016). Joint probability of precipitation and reservoir storage for drought estimation in the headwater basin of the Huaihe River, China. Stochastic Environmental Research and Risk Assessment. DOI 10.1007/s00477-016-1249-z. 22