Multivariate Joint Probability Function of Earthquake Ground Motion Prediction Equations Based on Vine Copula Approach

Abstract

In the structural earthquake engineering, a single parameter is often not sufficient enough to depict the severity of ground motions, and it is thus necessary to use multiple ones. In this sense, the correlation among multiple parameters is generally considered as an importance issue. The conventional approach for developing the correlation is based on regression analysis, along with simple pair copula approaches proposed in recent years. In this study, an innovative mathematical technique—vine copula—is firstly introduced to develop the empirical model for the multivariate dependence of pseudospectral accelerations (PSAs), which are the most commonly used earthquake ground motion parameters. This advancement not only offers a more flexible way of describing nonlinear dependence among multivariate PSAs from the marginal distribution functions but also highlights the extreme dependence. The results can be conventionally acquired in the ground motion selection and seismic risk and loss assessment based on multivariate parameters.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

In the structural earthquake engineering, a single ground motion parameter (GMP) is often not sufficient enough to characterize the severity of earthquake ground motions, and it is necessary to use multiple ones. Consequently, it is critical to evaluate the correlation among multiple GMPs when they are used to select ground motions [1, 2] or to calculate aggregated seismic losses of distributed infrastructures and portfolios [3, 4]. Attempts have been paid to study the correlation during the past few years [5–7]. As pseudospectral accelerations (PSAs) play an important role in antiearthquake design of structures, their correlation at different vibration periods is widely investigated [8, 9].

The correlation is computed based on the residuals of ground motion prediction equations (GMPEs) (e.g., [10–13]) derived from large number of ground motion records. Conventionally, the correlation model of earthquake parameters (such as PSAs at various vibration periods) is developed using regression analysis (such as [6–8]) on Pearson product-moment correlation coefficients which are derived from the residuals of the GMPEs. In recent years, copula techniques have been more and more widely applied in engineering [14–18] due to the advantages in the probabilistic analysis [19, 20]. A simple bivariate copula technique is introduced by Goda and Atkinson [21] to model the PSA interperiod dependence. It demonstrates that PSAs are marginally lognormal distributed. It also validates that the conventional two-step approach is appropriate to develop the correlation model. However, these techniques can only describe a linear correlation and take a bivariate interperiod dependence of parameters into account. In addition, the study of Weatherill et al. [3] shows that the interperiod dependence of PSA results in a larger difference of losses at lower annual probabilities of exceedance. Given the detail dependence has significant impacts on the low probability of risk, the detail dependence of PSAs should be highlighted in the risk assessment.

Multivariate elliptical copulas, such as normal copula and t copula, have been used for modeling multivariate dependence [22–24], including for earthquake GMPs [21, 25]. However, multivariate normal copula is unable to capture extreme dependence because of its independent property in the tail region. In the meanwhile, the multivariate t copula has been criticized for using only a single parameter, i.e., degree of freedom, to determine tail dependence. Thus, it has limited abilities to describe the complicated tail dependence in the multivariate context. Furthermore, the vine copula approach, which is based on the decomposition of different bivariate copulas, can perform better in high-dimensional cases since it can embrace heterogeneous dependence structures among variables and can use a series of pair copulas to capture a complex relationship [26–29].

Therefore, this study aims to develop a multivariate joint probability function of PSAs using vine copula technique. Instead of using pairwise interperiod dependence, we develop a multivariate dependence structure among PSAs at different periods. We also investigate the tail dependence. In this paper, firstly, we introduce the vine copula. Then, multivariate dependence structure of PSAs is derived from a ground motion record database. Since there are two orthogonal horizontal components for each ground motion record, we use the geometric mean of two components at different periods to calculate the interperiod dependence.

2. Vine Copula Dependence

In this section, we will provide an in-depth discussion about vine copula. We will review some basic concepts of vine copula methodology, namely, the definitions, the properties (i.e., tail dependence), and some widely used bivariate copula families which are nested in trees of vine copula.

2.1. Vine Copula Construction

A $K$-dimensional copula $C:{\left[\mathrm{0,1}\right]}^K⟶\left[\mathrm{0,1}\right]$ is a cumulative distribution function on ${\left[\mathrm{0,1}\right]}^K$ with uniform margins on [0, 1]. Let $F$ be a $K$-dimensional distribution function with marginal distributions $F_1$, $F_2$,…, $F_K$. According to Sklar’s theorem [30], there exists a $K$-dimensional copula C, such that $\forall \left(x_1,x_2,\dots ,x_K\right)\in {ℝ}^K$:$$\begin{array}{}\text{(1)}& F\left(x_1,x_2,\mathrm{...},x_K\right)=C\left(F_1\left(x_1\right),F_2\left(x_2\right),\dots ,F_K\left(x_K\right)\right).\end{array}$$

If $F_i,i=1,\dots ,K$, are continuous, copula $C$ is unique. In a continuous case, after taking derivatives on both sides of equation (1), we obtain the joint density representation of $X$ as follows:$$\begin{array}{c}\text{(2)}& f\left(x_1,x_2,\dots ,x_K\right)=\frac{{\partial }^{K}F\left(x_1,\mathrm{...},x_K\right)}{\partial x_1\cdots \partial x_N}=\frac{{\partial }^{K}C\left(F_1\left(x_1\right),\dots ,F_K\left(x_K\right)\right)}{\partial F_1\left(x_1\right)\cdots \partial F_N\left(x_N\right)}\times f_1\left(x_1\right)\times \cdots \times f_K\left(x_K\right)=c\left(F_1\left(x_1\right),\dots ,F_K\left(x_K\right)\right)\times \prod _{i=1}^K{f}_i\left(x_i\right).\end{array}$$

Joe [26] noticed that a $K$-dimensional copula can be represented by $K\left(K-1\right)/2$ different bivariate copulas through the decomposition method. Studies [28, 29] introduced the vine copula method on the basis of decomposition approaches. This method builds multivariate copulas based on the product of a series of simple pair copulas. Aas et al. [31] applied this PCC (pair copula construction) method to model the multiple dependence structure of returns and produced an entire algorithm relating to model estimation and simulation.

Vine copula approach allows us to combine different families of bivariate copulas for different pairs of margins and higher order dependencies. Fischer et al. [27] compared this approach with other methods and found that the vine copula method performed better in high-dimensional cases. Other studies [32–34] pointed out the same conclusion.

Through vine copula decomposition, a multivariate density $f\left(x_1,x_2,\dots ,x_K\right)$ can be represented as a product of pair copula densities and marginal densities. For this reason, appropriate bivariate copulas can capture the variant dependence relationship of different pairs of margins instead of using only one multidimensional copula to describe the whole dependence structure of the multivariate density.

In general, an $K$-dimensional vine structure is represented by $\left(K-1\right)$ trees. The $j$-th tree $T_j$ has $\left(K+1-j\right)$ nodes and $\left(K-j\right)$ edges, with each edge corresponding to a pair copula density. The edges in $T_j$ are the nodes in $T_{j+1}$. For a $K$-dimensional multivariate distribution, the decomposition is not unique. Bedford and Cooke [28] presented the multivariate density in terms of a regular vine, termed an R-vine, as follows:$$\begin{array}{}\text{(3)}& f\left(x_1,\dots ,x_K\right)=\prod _{k=1}^{K}f\left(x_k\right)\cdot \prod _{j=1}^{K-1}\prod _{e\in E_j}{c}_{\left.e\right|D_e}\left(F_{\left.e_a\right|D_e},F_{\left.e_b\right|D_e}\right),\end{array}$$where $E_j$ is a set of edges associated with tree $T_j$; $e=\left\{a,b\right\}$; and the edge connecting variables denoted with $e_a$ and $e_b$ in $T_j$ gave the variable labeled by $D_e$, and $D_e$ denotes conditioning set for the edge $e$; see Corollary 1 in the study of Bedford and Cooke [28] for full representation and more details. This is a general representation, and it does not restrict the vine structure to a particular pattern. The selection of the vine structure depends on the observed data.

Aas et al. [31] identified two types of vines with particular structures, namely, D-vine and C-vine. As shown in equation (4), each tree of D-vine is a path, and each node is connected to no more than two other nodes in each tree. Meanwhile, C-vine shows a star structure. In each tree, there exists one unique node connecting with all other nodes (equation (5)):$$\begin{array}{}\text{(4)}& f\left(x_1,\dots ,x_K\right)=\prod _{k=1}^{K}f\left(x_k\right)\prod _{j=1}^{K-1}\prod _{i=1}^{K-j}{c}_{i,\left.j+i\right|i+1,\dots ,i+j-1}\left(F\left(\left.x_i\right|x_{i+1},\dots ,x_{i+j-1}\right),F\left(\left.x_{j+i}\right|x_{i+1},\dots ,x_{i+j-1}\right)\right),\text{(5)}& f\left(x_1,\dots ,x_K\right)=\prod _{k=1}^{K}f\left(x_k\right)\prod _{j=1}^{K-1}\prod _{i=1}^{K-j}{c}_{j,\left.j+i\right|1,\dots ,j-1}\left(F\left(\left.x_j\right|x_1,\dots ,x_{j-1}\right),F\left(\left.x_{j+i}\right|x_1,\dots ,x_{j-1}\right)\right).\end{array}$$

It has been shown that vine copula decomposition offers a great deal of flexibility in modeling a complex dependence structure, especially in relation to the tail dependence, compared with the traditional multivariate copula (see [27, 35, 36] for details). In high-dimensional cases ($K\ge 3$), various pairs of variables may exhibit heterogeneous dependence patterns that traditional multivariate copulas are unable to capture. Vine copulas can solve this problem and provide good fit for high-dimensional data. We can fit the standardized residuals with an appropriate vine structure and select a suitable bivariate copula to describe the dependence patterns for each pair of margins.

In this paper, we fit the margins of the PSA residuals ${\epsilon }_k$ with normal distributions and select the vine structure and the pair copula families using the Akaike information criterion. Finally, the parameters of the overall model are estimated using the maximum likelihood method:$$\begin{array}{}\text{(6)}& \left({\epsilon }_1,{\epsilon }_2,\mathrm{...},{\epsilon }_K\right)\sim \text{VC}\left({\mathcal{V}}_{\mathcal{S}},B\left({\mathcal{V}}_{\mathcal{S}}\right),\theta \left(B\right)\right),\end{array}$$where ${\mathcal{V}}_{\mathcal{S}}$ denotes the vine structure, $B\left({\mathcal{V}}_{\mathcal{S}}\right)$ is a collection of pair copula families, and $\theta \left(B\right)$ denotes the parameters of the copulas, respectively.

2.2. Tail-Dependence Coefficients

Venter [37] investigated the tail concentration functions for different copula families and suggested to select copulas for a given dataset using the tail concentration characteristic. The tail dependence is expressed in terms of a conditional probability that one variable $X_1$ will incur a large loss (or gain), given that another variable $X_2$ also experiences a large loss (or gain). Considering two random variables $X_1$ and $X_2$, with joint continuous cumulative distribution function $F$, copula $C$, and margins $F_{X_1}\text{\hspace{0.17em}and\hspace{0.17em}}{F}_{X_2}$, the lower tail-dependence coefficient is defined as$$\begin{array}{}\text{(7)}& {\lambda }_L=\underset{u⟶0^{+}}{\lim }\mathbb{P}\left(F_{X_1}\left(x_1\right)\le \left.u\right|F_{X_2}\left(x_2\right)\le u\right).\end{array}$$

The upper tail-dependence coefficient is defined as$$\begin{array}{}\text{(8)}& {\lambda }_U=\underset{u⟶1^{-}}{\lim }\mathbb{P}\left(F_{X_1}\left(x_1\right)>\left.u\right|F_{X_1}\left(x_1\right)>u\right).\end{array}$$

If $F_{X_1}$ and $F_{X_2}$ are continuous, ${\lambda }_L$ and ${\lambda }_U$ can be expressed in terms of a copula representation. For a lower tail-dependence coefficient, formula (7) can be rewritten as$$\begin{array}{}\text{(9)}& {\lambda }_L=\underset{u⟶0^{+}}{\lim }\frac{\mathbb{P}\left(F_{X_1}\left(x_1\right)\le u,F_{X_2}\left(x_2\right)\le u\right)}{\mathbb{P}\left(F_{X_2}\left(x_2\right)\le u\right)}=\underset{u⟶0^{+}}{\lim }\frac{C\left(u,u\right)}{u}.\end{array}$$

Analogously, for an upper tail-dependence coefficient, we have$$\begin{array}{}\text{(10)}& {\lambda }_U=\underset{u⟶1^{-}}{\lim }\frac{\mathbb{P}\left(F_{X_1}\left(x_1\right)>u,F_{X_2}\left(x_2\right)>u\right)}{\mathbb{P}\left(F_{X_2}\left(x_2\right)>u\right)}=\underset{u⟶1^{-}}{\lim }\frac{1-2u+C\left(u,u\right)}{1-u}.\end{array}$$

If $C$ is radially symmetric, ${\lambda }_L={\lambda }_U$ (see [20] for the proof). Intuitively, if ${\lambda }_L$ and ${\lambda }_U$ exist and fall in $\left(\mathrm{0,1}\right]$, $X_1$ and $X_2$ show lower or upper tail-dependence. If ${\lambda }_L$ and ${\lambda }_U$ are equal to 0, one can say that the two variables are independent in the tails; hence, extreme events seem to occur independently. We can describe different tail-dependent behavior by choosing an appropriate copula model.

In the aggregated seismic losses assessment of distributed infrastructures and portfolios, the differences between the losses considering and without considering the dependence of pseudospectral accelerations (PSAs) tend to be increasing with decreasing annual probability of exceedance [3]. Hence, the tail dependence among pseudospectral accelerations (PSAs) should be paid more attention to the loss estimation at extreme events. Herein, we can use the tail-dependence coefficient to measure the concordance between the extreme events of different random parameters.

2.3. Bivariate Copula Families

In this paper, we use the vine copula approach to measure the interperiod dependence structure of PSAs. Vine copula is a “pair copula construction” method; hence, we focus on the selectable two-dimensional copula families. The elliptical copulas related to an elliptical distribution are the most widely used in many research fields, (see [21]).

2.3.1. Gaussian Copula

In the bivariate case, the Gaussian copula is defined by the following expression:$$\begin{array}{}\text{(11)}& C_{\rho }^G\left(u,v\right)={\Phi }_{\rho }\left({\Phi }^{-1}\left(u\right),{\Phi }^{-1}\left(v\right)\right)=\frac{1}{2\pi \sqrt{1-{\rho }^2}}{\int }_{-\infty }^{{\Phi }^{-1}\left(u\right)}{\int }_{-\infty }^{{\Phi }^{-1}\left(v\right)}\exp \left(\frac{-\left(s^2-2\rho st+t^2\right)}{2\left(1-{\rho }^2\right)}\right)\text{d}s\text{d}t,\end{array}$$where ${\Phi }_{\rho }$ is the bivariate normal cumulative distribution function with linear correlation coefficient $\rho \in \left[\mathrm{0,1}\right]$ and $\Phi$ is the standard normal cumulative distribution function and ${\Phi }^{-1}$ is its inverse function. We can see that the bivariate Gaussian copula density is symmetric, and it has weak capability to capture skewness in the dependence structure. If we go far into the tail, the extreme events tend to be independent, even though we choose a very high correlation. Actually, the Gaussian copula shows asymptotical independence in the tail regions, more details can be found in [20].

2.3.2. $t$ Copula

The $t$ copula corresponds to a Student t-distribution. It is defined as$$\begin{array}{}\text{(12)}& C_{\nu ,\rho }^t\left(u,v\right)=t_{\nu ,\rho }\left(t_{\nu }^{-1}\left(u\right),t_{\nu }^{-1}\left(v\right)\right)=\frac{\Gamma \left(\nu /2+1\right)}{\Gamma \left(\nu /2\right)\pi \nu \sqrt{1-{\rho }^2}}{\int }_{-\infty }^{t_{\nu }^{-1}\left(u\right)}{\int }_{-\infty }^{t_{\nu }^{-1}\left(v\right)}\left(+\frac{-\left(s^2-2\rho st+t^2\right)}{2\left(1-{\rho }^2\right)}\right)\text{d}s\text{d}t,\end{array}$$where $t_{\nu ,\rho }$ is the cumulative distribution function of a two-dimensional $t$ distribution, $\nu$ are the degrees of freedom, and $\rho$ is a measure of dependence. The $t$ copula also has symmetric shape, upper and lower tail dependences are identical, and they are completely determined by $\nu$ and $\rho$. When $\nu$ gets large, the $t$ copula decays to a Gaussian copula. The expression of ${\lambda }_L$ and ${\lambda }_U$ is as follows:$$\begin{array}{}\text{(13)}& {\lambda }_L={\lambda }_U=2t_{\nu +1}\left(-\frac{\sqrt{\left(\nu +1\right)\left(1-\rho \right)}}{\sqrt{\rho +1}}\right),\end{array}$$where $t_{\nu +1}$ is the cumulative distribution function of a univariate t distribution with $\nu +1$ degrees of freedom [20].

Archimedean copulas, defined by their generator functions, are also used intensively. Generally, if a function $\phi :\left[\mathrm{0,1}\right]⟶\left[0,\infty \right]$ with a continuously decreasing and convex derivative, it can be considered as a generator function of an Archimedean copula. By definition, a $n$-dimensional Archimedean copula has the following expression:$$\begin{array}{}\text{(14)}& C\left(u_1,u_2,\dots ,u_K\right)={\phi }^{-1}\left(\phi \left(u_1\right)+\phi \left(u_2\right)+\cdots +\phi \left(u_K\right)\right),\end{array}$$

Different generator functions create different Archimedean copulas. More details about the generator function can be found in the studies of Joe [38] and Nelsen [39]. In the bivariate case, the copula function is defined by$$\begin{array}{}\text{(15)}& C\left(u,v\right)={\phi }^{-1}\left(\phi \left(u\right)+\phi \left(v\right)\right),\text{\hspace{0.17em}if\hspace{0.17em}}\phi \left(u\right)+\phi \left(v\right)\le \phi \left(0\right),\end{array}$$where $\phi \left(u\right)$ is a $C^2$ function with $\phi \left(1\right)=0,{\phi }^{\prime }<0$ and ${\phi }^{″}>0$.

2.3.3. Frank Copula

The Frank copula is defined by$$\begin{array}{}\text{(16)}& C^F\left(u,v;\theta \right)=-\frac{1}{\theta }\ln \left(\frac{\left(e^{-\theta u}-1\right)\left(e^{-\theta v}-1\right)}{e^{-\theta }-1}+1\right),\hspace{1em}\theta \in \left(-\infty ,0\right)\cup \left(0,+\infty \right).\end{array}$$

The generator function is ${\phi }_{\theta }\left(t\right)=-\ln \left(\left(e^{-\theta t}-1\right)/\left(e^{-\theta }-1\right)\right)$.

Similar to the Gaussian copula, the Frank copula is symmetric in both tails and it is not sensitive to the relationship between the extreme values in both upper and lower tails. It shows asymptotic independence in the tails, whereas it has a strong dependence in the center of the distribution. This means that the Frank copula fails to capture tail-dependence behavior; hence, it is suitable to use the Frank copula when the tail dependence of a given dataset is relatively weak. If $\theta ⟶0$, the Frank copula will decay to an independent copula as a special case.

Joe [38] provided some examples of two-parameter bivariate copulas, such as Joe’s BB1, BB4, and BB8 copula. Two-parameter copulas are distinguished from other bivariate copula families mentioned above. They show a high flexibility in modeling bivariate dependence structures through two parameters, especially in modeling an asymmetric tail-dependence behavior. More details about these copulas can be found in the study of Joe et al. [40, 41].

2.3.4. Joe’s BB1 Copula

The BB1 copula is defined as$$\begin{array}{}\text{(17)}& C^{\text{BB}1}\left(u,v;\theta ,\delta \right)={\left\{1+{\left[{\left(u^{-\theta }-1\right)}^{\delta }+{\left(v^{-\theta }-1\right)}^{\delta }\right]}^{1/\delta }\right\}}^{-\left(1/\theta \right)},\end{array}$$where $\theta >0\text{\hspace{0.17em}and\hspace{0.17em}}\delta \ge 1$. The lower tail-dependence coefficient ${\lambda }_L$ is $2^{-\left(1/\delta \theta \right)}$, and the upper tail-dependence coefficient ${\lambda }_U$ is equal to $2-2^{1/\delta }$ which is independent of $\theta$. The concordance of two variables increases as $\theta$ increases. The Gumbel copula is the limiting case of BB1 copula as $\theta ⟶0$. Obviously, BB1 copula decays to the Clayton copula when $\delta$ is equal to 1 and decays to independent copula as $\theta ⟶0$ and $\delta ⟶1$.

2.3.5. Joe’s BB8 Copula

The BB8 copula is defined as$$\begin{array}{}\text{(18)}& C^{\text{BB}8}\left(u,v;\theta ,\delta \right)={\delta }^{-1}\left[1-{\left\{1-{\left[1-{\left(1-\delta \right)}^{\theta }\right]}^{-1}\left[1-{\left(1-\delta u\right)}^{\theta }\right]\left[1-{\left(1-\delta v\right)}^{\theta }\right]\right\}}^{\left(1/\theta \right)}\right],\end{array}$$where $\theta \ge 1\text{\hspace{0.17em}and\hspace{0.17em}}0\le \delta \le 1$. The independent copula is obtained as $\theta ⟶1$ or $\delta ⟶0$. BB8 copula decays to the Joe copula when $\delta$ is equal to 1, while the Frank copula is obtained as $\theta ⟶\infty$, and the single parameter can be calculated by using the formula $1-{\left(1-\delta \right)}^{\theta }$. BB8 copula dose not exhibit tail dependence except when $\delta$ is equal to 1.

2.3.6. Tawn Copula

The Tawn copula has been introduced by Tawn [42], regarded as an extension of the Gumbel copula with three parameters, and it can be expressed by$$\begin{array}{}\text{(19)}& C^{Ta}\left(u,v;\theta ,{\delta }_1,{\delta }_2\right)=\left(1-{\delta }_1\right)t-\left(1-{\delta }_2\right)\left(1-t\right)+{\left[{\left({\delta }_1\left(1-t\right)\right)}^{\theta }+{\left({\delta }_{2}t\right)}^{\theta }\right]}^{\left(1/\theta \right)},\end{array}$$where $t=\log \left(v\right)/\log \left(uv\right)$, $0\le {\delta }_1,{\delta }_2\le 1$, and $\theta \in \left[1,+\infty \right)$. If ${\delta }_1={\delta }_2$, the Gumbel copula is obtained from the Tawn copula. ${\delta }_1$ and ${\delta }_2$ can be interpreted as skewness parameters. Hence, the Tawn copula can be divided into two types, and each type has one of the asymmetry parameters fixed to 1 so that the corresponding copula density is either left or right skewed (i.e., ${\delta }_1=1$ or ${\delta }_2=1$). And the upper tail-dependence coefficient has following expression:$$\begin{array}{c}\text{(20)}& {\lambda }_U={\delta }_1+1-{\left({\delta }_1^{\theta }+1\right)}^{\left(1/\theta \right)}\text{or\hspace{0.17em}}{\lambda }_U={\delta }_2+1-{\left({\delta }_2^{\theta }+1\right)}^{\left(1/\theta \right)}.\end{array}$$

3. Multivariate Dependence Structure of Ground Motion Parameters

3.1. Ground Motion Data and Residuals

In this section, we investigate the multivariate dependence of pseudospectral accelerations (PSAs) using vine copula technique. The residuals of PSAs calibrated with equation (21) are used to model the multivariate dependence structure based on vine copula:$$\begin{array}{}\text{(21)}& r_{\text{PSA}}=\ln \left({\text{PSA}}_{\text{obs}}\right)-f\left(M,R,\theta \right),\end{array}$$where $r_{\text{PSA}}$ represents the residual values of pseudospectral accelerations (PSAs); ${\text{PSA}}_{\text{obs}}$ is the observed PSAs; and $f\left(M,R,\theta \right)$ is the predicted median PSA values in logarithmic space calculated through ground motion prediction equations (GMPEs) [10] as the function of magnitude M, fault-to-site distance R, and other sets of variables $\theta$. To illustrate, this study only focuses on the 5%-damped pseudospectral accelerations (PSAs) at vibration periods of 0.1, 0.2, 1.0, 2.0, and 4.0 sec. But the conclusion drawn in the paper can also be applied to the PSAs at other different vibration periods. We use the same set of ground motion record data for the analysis as that used for the development of ground motion prediction equations (GMPEs) used by Campell and Bozorgnia [10]. It consists of 1550 pairs of components of accelerograms, which are used to calculate ${\text{PSA}}_{\text{obs}}$. $f\left(M,R,\theta \right)$ is computed using the GMPEs proposed by Campell and Bozorgnia [10].

3.2. Vine Copula-Based Multivariate Dependence Structure

3.2.1. Multivariate Copula Calibration

The lognormality of parameter PSAs has been well documented in many literatures (see [21] for more details).

First, we use the empirical cumulative distribution function to transform residual values of PSA into the so-called pseudoobservations in the domain of [0, 1] by the following formula:$$\begin{array}{c}\text{(22)}& u_{ki}\left(r_{ki}\right)=\frac{1}{N+1}\sum _{i=1}^N{1}_{\left\{r_{k\cdot }\le r_{ki}\right\}},k=1,\dots ,5,\\ i=1,\dots ,N,\end{array}$$where $1_{\left\{\cdot \right\}}$ denotes the indicator function which takes value of 1 if $r_{k\cdot }\le r_{ki}$ and value of 0 otherwise and $r_1$, $r_2$, $r_3$, $r_4$, and $r_5$ denote residuals of PSAs at vibration periods of 0.1, 0.2, 1.0, 2.0, and 4.0 sec, respectively. Then, the residuals $r_k$ are transformed to uniform data $u_k$ in domain [0, 1], and the pure dependence structure among the residuals at five vibration periods can be captured by the copula model using these uniform data, eliminating the effect of margins at each period. Figure 1 displays the bivariate dependence feature for each pair of $u_1$, $u_2$, $u_3$, $u_4$, and $u_5$. The histograms on the diagonal panel indicate that transformed data are perfectly uniformly distributed. In the meanwhile, the bivariate contour plots and scatter plots reveal that different pairs of variables show heterogeneous relationship characters, i.e., asymmetry and dependence in tail regions. Notice that the traditional multivariate normal copula and t copula fail to handle these complex dependence patterns.

[figure omitted; refer to PDF]

In order to check whether vine copula shows a more flexibility for modeling multiple dependences than elliptical copula in the high-dimensional case, we fit the residuals with the normal copula, t copula, and D-vine copula, respectively. The estimation results of five-dimensional normal copula and t copula are reported in Table 1. We can observe that all parameters are highly significant, and the off-diagonal elements $\rho$ of correlation matrix are similar between these two multivariate elliptical copulas. We then fit the residuals with a five-dimensional D-vine copula. The structure of D-vine is defined by using equation (4), and the pair copula families have been restricted to the bivariate t copula for the reason of comparison. Since the five-dimensional t copula is nested in the D-vine copula structure, the likelihood ratio test can be performed between these two copulas [22]. The results of D-vine copula are reported in Table 2. In each tree, the residuals of PSAs are connected as a path from the shorter period to longer period. And for each pair, tail-dependence coefficients can also be calculated.

Table 1

Five-dimensional normal copula and t copula.

Normal	Par.	Kendall’s $\tau$	t copula	Par.	Kendall’s $\tau$
${\rho }_{\mathrm{1,2}}$	0.8755${}^{\ast \ast \ast }$	0.6789	${\rho }_{\mathrm{1,2}}$	0.8819${}^{\ast \ast \ast }$	0.6785
	(0.0040)			(0.0050)
${\rho }_{\mathrm{1,3}}$	0.3624${}^{\ast \ast \ast }$	0.2361	${\rho }_{\mathrm{1,3}}$	0.3705${}^{\ast \ast \ast }$	0.2416
	(0.0200)			(0.0210)
${\rho }_{\mathrm{1,4}}$	0.1639${}^{\ast \ast \ast }$	0.1048	${\rho }_{\mathrm{1,4}}$	0.1623${}^{\ast \ast \ast }$	0.1038
	(0.0230)			(0.0240)
${\rho }_{\mathrm{1,5}}$	0.0282	0.0180	${\rho }_{\mathrm{1,5}}$	0.0207	0.0132
	(0.0240)			(0.0250)
${\rho }_{\mathrm{2,3}}$	0.4516${}^{\ast \ast \ast }$	0.2983	${\rho }_{\mathrm{2,3}}$	0.4500${}^{\ast \ast \ast }$	0.2972
	(0.018)			(0.0190)
${\rho }_{\mathrm{2,4}}$	0.2560${}^{\ast \ast \ast }$	0.1648	${\rho }_{\mathrm{2,4}}$	0.2464${}^{\ast \ast \ast }$	0.1585
	(0.0220)			(0.0230)
${\rho }_{\mathrm{2,5}}$	0.1077${}^{\ast \ast \ast }$	0.0687	${\rho }_{\mathrm{2,5}}$	0.0919${}^{\ast \ast \ast }$	0.0586
	(0.0240)			(0.0250)
${\rho }_{\mathrm{3,4}}$	0.7813${}^{\ast \ast \ast }$	0.5709	${\rho }_{\mathrm{3,4}}$	0.7754${}^{\ast \ast \ast }$	0.5649
	(0.0080)			(0.0080)
${\rho }_{\mathrm{3,5}}$	0.5905${}^{\ast \ast \ast }$	0.4021	${\rho }_{\mathrm{3,5}}$	0.5776${}^{\ast \ast \ast }$	0.3920
	(0.0140)			(0.0150)
${\rho }_{\mathrm{4,5}}$	0.8111${}^{\ast \ast \ast }$	0.6023	${\rho }_{\mathrm{4,5}}$	0.8069${}^{\ast \ast \ast }$	0.5977
	(0.0070)			(0.0070)
			$\upsilon$	20.6782${}^{\ast \ast \ast }$
				(3.6160)

The subscripts 1, 2, 3, 4, and 5 denote PSAs at 0.1, 0.2, 1, 2, and 4 sec, respectively. ${}^{\ast \ast \ast }$Statistical significance at the 1% level.

Table 2

Five-dimensional D-vine copula restricted with pair t copula.

Copula	Family	Par. 1 (ρ)	Par. 2 (υ)	${\lambda }_U$	${\lambda }_L$	Kendall’s $\tau$
$C_{\mathrm{1,2}}$	T	0.8819${}^{\ast \ast \ast }$	5.5335${}^{\ast \ast \ast }$	0.5505	0.5505	0.6875
		(0.0940)	(0.8112)
$C_{\mathrm{2,3}}$	T	0.4488${}^{\ast \ast \ast }$	16.4769${}^{\ast \ast }$	0.0192	0.0192	0.2963
		(0.0197)	(6.9061)
$C_{\mathrm{3,4}}$	T	0.7776${}^{\ast \ast \ast }$	30.0000${}^{\ast \ast \ast }$	0.0579	0.0579	0.5671
		(0.0079)	(12.2744)
$C_{\mathrm{4,5}}$	T	0.8082${}^{\ast \ast \ast }$	22.9407${}^{\ast \ast \ast }$	0.1242	0.1242	0.5992
		(0.0547)	(14.2751)
$C_{\left.\mathrm{1,2}\right|3}$	T	−0.0568${}^{\ast \ast }$	16.3974${}^{\ast \ast \ast }$	0.0004	0.0004	−0.0362
		(0.0274)	(6.8677)
$C_{\left.\mathrm{2,4}\right|3}$	T	−0.1824${}^{\ast \ast \ast }$	30.0000${}^{\ast \ast \ast }$	0.0000	0.0000	−0.1167
		(0.1497)	(12.9389)
$C_{\left.\mathrm{3,5}\right|4}$	T	−0.1272${}^{\ast \ast \ast }$	30.0000${}^{\ast \ast \ast }$	0.0000	0.0000	−0.0812
		(0.2539)	(11.9427)
$C_{\left.\mathrm{1,4}\right|\mathrm{2,3}}$	T	−0.1135${}^{\ast \ast \ast }$	30.0000${}^{\ast \ast }$	0.0000	0.0000	−0.0724
		(0.6598)	(11.2320)
$C_{\left.\mathrm{2,5}\right|\mathrm{3,4}}$	T	−0.1469${}^{\ast \ast \ast }$	30.0000${}^{\ast \ast \ast }$	0.0000	0.0000	−0.0938
		(0.6181)	(13.0847)
$C_{\left.\mathrm{1,5}\right|\mathrm{2,3,4}}$	T	−0.0504${}^{\ast }$	30.0000${}^{\ast }$			−0.0362
		(0.0262)	(15.6792)

The subscripts 1, 2, 3, 4, and 5 denote PSAs at 0.1, 0.2, 1, 2, and 4 sec, respectively. T denotes the pair t copula. ${}^{\ast }$Statistical significance at the 10% level; ${}^{\ast \ast }$statistical significance at the 5% level; ${}^{\ast \ast \ast }$statistical significance at the 1% level.

Then, we calculate quantitative measure, i.e., Akaike information criterion (AIC) and Bayesian information criterion (BIC). The three fitted multivariate copula models for residuals of PSAs are compared based on two criterions: (1) the lower values of AIC and BIC mean better level of goodness of fit; (2) the greater log-likelihood value indicates the better level of goodness of fit. The results are displayed in Table 3. The five-dimensional D-vine copula shows the lowest AIC and BIC values and the biggest log-likelihood value. We also perform the likelihood ratio test between two nested models, i.e., t copula and D-vine copula. The statistic equals to 25.978 (2922.2330–2896.2550) with degree of freedom 9 (20–11), and the associated $p$ value is almost zero. It implies that five-dimensional t copula can be rejected in favor of D-vine copula for our data.

Table 3

AIC values, BIC values, and log-likelihood values of five-dimensional normal copula, five-dimensional t copula, and five-dimensional D-vine copula restricted with pair t copula.

	Normal copula	t copula	D-vine copula with pair t copula
AIC	−5736.3750	−5770.5100	−5804.4660
BIC	−5683.0450	−5711.8470	−5834.3630
Log-likelihood value	2878.1870	2896.2550	2922.2330

The AIC value, BIC value, and log-likelihood values are shown in bold if the corresponding copula is preferred.

We restrict the structure of vine copula to a D-vine structure and limit the pair copula to bivariate t copula in order to compare three different copulas in the previous section. Then we fit transformed residual data with vine copula without restrictions. The selection of vine structure and the choice of pair copula families are data oriented. We choose the best vine copula model associated with the smallest AIC value. The suitability for selecting vine copula model has been shown (see [41]). And the independent test for bivariate copula has also been performed. The independent test for bivariate copula $T=\sqrt{9n\left(n-1\right)/2\left(2n+5\right)}\left|\widehat{\tau }\right|$, where n is the number of observations and $\widehat{\tau }$ is the empirical Kendall’s $\tau$ of the data u and $v$. The test statistic T is asymptotic normal distributed under the null hypothesis that the bivariate variable is independent (see [43]). Figure 2 illustrates the intuitive graphs for four trees of the vine structure. The nodes denote the margins, and the edges denote the bivariate copula between two linked nodes. Furthermore, the labels on the edges are the families and Kendall’s $\tau$ of corresponding pair copulas. Figure 2 also indicates that the fitted vine is D-vine structure and no node plays a major role in the whole structure. However, its nested pair copulas exhibit a heterogeneous characteristic, which is not only described by bivariate t copula but also by other two-parameter copulas.

[figure omitted; refer to PDF]

The estimated results are reported in Table 4, and all parameters are highly significant. The pair copula families and copula parameters are estimated by the joint maximum likelihood method. Compared to the sequential estimation method, this method can provide more precise results since all the parameters are estimated simultaneously instead of only the bivariate scenario involved. Kendall’s $\tau$ and corresponding tail-dependence coefficients are displayed in Table 4. The AIC and BIC scores decrease significantly compared with other three multivariate copulas mentioned previously (i.e., multivariate normal copula, multivariate t copula, and D-vine restricted with pair t copula). Hence, it performs better at modeling PSA dependence structure at different vibration periods. It is notable that some pairs reveal strong tail dependence. In Tree 1, the dependence between PSAs at 0.1 s and 0.2 s is captured by the survival Tawn copula [38]. This two-parameter copula can describe a lower tail dependence. It implies that PSAs at 0.1 s and 0.2 s show very strong co-movement probability in left tail region, i.e., an extreme small PSA at 0.1 s tends to accompany with an extreme small PSA at 0.2 s, while PSAs at 0.2 s and 1 s vibration periods show a slight symmetric tail dependence which is captured by t copula. PSAs at 1 s and 2 s do not exhibit any unusual character in tail regions since normal copula has been chosen for this pair. The dependence between PSAs at 2 s and 4 s vibration periods is described by BB1 copula (Clayton–Gumbel) which can capture the asymmetric tail dependence [38]. The co-movement of the extreme small values is greater than that of the extreme large values. The remaining pair copulas do not show any unusual character in tail regions. The density of corresponding pair copulas for PSA residuals at different vibration periods is illustrated in Figures 3–5. The traditional multivariate normal copula and t copula fail to describe multivariate distribution as comprehensive as this one.

Table 4

Estimated copula parameters of vine copula for PSA residuals at different vibration periods.

Copula	Family	Par. 1	Par. 2	${\lambda }_U$	${\lambda }_L$	Kendall’s $\tau$
$C_{\mathrm{1,2}}$	Tawn1 (180°)	3.4608${}^{\ast \ast \ast }$	0.9346${}^{\ast \ast \ast }$	—	0.7512	0.6745
		(0.0940)	(0.0139)
$C_{\mathrm{2,3}}$	T	0.4522${}^{\ast \ast \ast }$	17.1688${}^{\ast \ast }$	0.0173	0.0173	0.2987
		(0.0197)	(8.0943)
$C_{\mathrm{3,4}}$	N	0.7799${}^{\ast \ast \ast }$	—	—	—	0.5695
		(0.0079)	—
$C_{\mathrm{4,5}}$	BB1z (180°)	0.3726${}^{\ast \ast \ast }$	2.0394${}^{\ast \ast \ast }$	0.4016	0.5952	0.5867
		(0.0547)	(0.0627)
$C_{\left.\mathrm{1,3}\right|2}$	I	—	—	—	—	—
		—	—
$C_{\left.\mathrm{2,4}\right|3}$	F	−1.1098${}^{\ast \ast \ast }$	—	—	—	−0.1218
		(0.1497)	—
$C_{\left.\mathrm{3,5}\right|4}$	BB8 (270°)	−1.4737${}^{\ast \ast \ast }$	−0.7718${}^{\ast \ast \ast }$	—	—	−0.0922
		(0.2539)	(0.1529)
$C_{\left.\mathrm{1,4}\right|\mathrm{2,3}}$	BB8 (90°)	−1.7090${}^{\ast \ast \ast }$	−0.5839${}^{\ast \ast }$	—	—	−0.0791
		(0.6598)	(0.2916)
$C_{\left.\mathrm{2,5}\right|\mathrm{3,4}}$	BB8 (90°)	−1.8660${}^{\ast \ast \ast }$	−0.6173${}^{\ast \ast \ast }$	—	—	−0.1055
		(0.6181)	(0.2283)
$C_{\left.\mathrm{1,5}\right|\mathrm{2,3,4}}$	I	—	—		—	—
AIC		−5909.0300
BIC		−5834.3600
Log-likelihood		2968.5130

The subscripts 1, 2, 3, 4, and 5 denote PSAs at 0.1, 0.2, 1, 2, and 4 sec, respectively. N, T, F, I, Tawn, BB1, and BB8 the denote normal copula, t copula, Frank copula, independent copula, Tawn copula, BB1 copula, and BB8 copula, respectively. The values in brackets indicate that the fitted pair copula is rotated by original copula, i.e., 90°, 180°, and 270°. The standard errors are in parentheses. Independent copula, normal copula, and Frank copula do not exhibit an upper and lower tail dependence (i.e., ${\lambda }_U$ and ${\lambda }_L$). ${}^{\ast \ast }$Statistical significance at the 5% level; ${}^{\ast \ast \ast }$statistical significance at the 1% level.

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

3.2.2. Joint Distribution Modeled by Vine Copula

The joint distribution of PSAs at five different vibration periods are obtained by combining the marginal distribution of PSAs at each period and its best fitted vine copula distribution. In particular, the joint density is the product of marginal density and corresponding vine copula density defined by equation (2). The margins represent the information at each period, and the copula contains the information about pure dependence structure of PSAs at five vibration periods. The joint density is defined by the following equation:$$\begin{array}{c}\text{(23)}& f\left({\epsilon }_1,\dots ,{\epsilon }_5\right)=f_1\left({\epsilon }_1\right)\cdot f_2\left({\epsilon }_2\right)\cdot f_3\left({\epsilon }_3\right)\cdot f_4\left({\epsilon }_4\right)\cdot f_5\left({\epsilon }_5\right)\hspace{1em}⟶\text{\hspace{0.17em}marginals}\cdot c_{\mathrm{1,2}}\left(u_1,u_2\right)\cdot c_{\mathrm{2,3}}\left(u_2,u_3\right)\cdot c_{\mathrm{3,4}}\left(u_3,u_4\right)\cdot c_{\mathrm{4,5}}\left(u_4,u_5\right)\hspace{1em}⟶\text{\hspace{0.17em}edges\hspace{0.17em}of\hspace{0.17em}tree\hspace{0.17em}}1\\ \cdot c_{\left.\mathrm{1,3}\right|2}\left(u_{\left.1\right|2},u_{\left.3\right|2}\right)\cdot c_{\left.\mathrm{2,4}\right|3}\left(u_{\left.2\right|3},u_{\left.4\right|3}\right)\cdot c_{\left.\mathrm{3,5}\right|4}\left(u_{\left.3\right|4},u_{\left.5\right|4}\right)\hspace{1em}⟶\text{\hspace{0.17em}edges\hspace{0.17em}of\hspace{0.17em}tree\hspace{0.17em}}2\\ \cdot c_{\left.\mathrm{1,4}\right|\mathrm{2,3}}\left(u_{\left.1\right|\mathrm{2,3}},u_{\left.4\right|\mathrm{2,3}}\right)\cdot c_{\left.\mathrm{2,5}\right|\mathrm{3,4}}\left(u_{\left.2\right|\mathrm{3,4}},u_{\left.5\right|\mathrm{3,4}}\right)\hspace{1em}⟶\text{\hspace{0.17em}edges\hspace{0.17em}of\hspace{0.17em}tree\hspace{0.17em}}3\\ \cdot c_{\left.\mathrm{1,5}\right|\mathrm{2,3,4}}\left(u_{\left.1\right|\mathrm{2,3,4}},u_{\left.5\right|\mathrm{2,3,4}}\right),\hspace{1em}⟶\text{\hspace{0.17em}edges\hspace{0.17em}of\hspace{0.17em}tree\hspace{0.17em}}4\end{array}$$where $f_i\left(\cdot \right)$ denotes the marginal distribution of residuals ${\epsilon }_i$ of PSAs at i-th period, i = 1, …, 5, i.e., a normal distribution. The vine copula model parameters (i.e., corresponding series of pair copula densities) are reported in Table 4.

3.2.3. Vine Copula-Based Joint Distribution in Earthquake Engineering

In the seismic hazard and risk assessment of calculating the aggregated losses of portfolios or infrastructures, for different types of structures, it is necessary to use different fragility functions characterized by different ground motion parameters. In this sense, multiple ground motion parameters are embedded in the seismic hazard and risk assessment. Herein, we present an example to illustrate the application and performance of the proposed vine copula-based multivariate joint distribution function by adopting joint exceedance probability of multiple ground motion parameters. The joint exceedance probability of multiple parameters is defined as a probability that a set of parameters (X₁ to X_n) simultaneously and, respectively, exceeds a set of certain values (x₁ to x_n), shown as follows:$$\begin{array}{}\text{(24)}& P_{X_{1,}{X}_2,\dots ,X_n}\left(x_1,x_2,\dots ,x_n\right)=P\left(X_1>x_1,X_2>x_2,\dots ,X_n>x_n\right).\end{array}$$

In the example, an earthquake scenario is assumed: (1) a hypothetical site is located at a distance of 30 km to a point strike-slip earthquake source; (2) an earthquake with moment magnitude of 6 occurs at this source; and (3) soil condition at the site is characterized with V_S30 = 720 m/s. We investigate the PSAs at 0.1 s, 0.2 s, 1 s, 2 s, and 4 s at this site. The median values of the PSAs in logarithmic space, namely, $f\left(M,R,\theta \right)$, are calculated using the GMPEs proposed by Campell and Bozorgnia [10]. Based on the proposed vine copula-based multivariate joint probability function, we use the Monte Carlo method to generate 20,000 realizations of the jointly distributed multivariate residuals of PSAs ($r_{\text{PSA}}$) at the given site. We obtain the final realizations of the PSAs in logarithmic space by summing up the mean PSAs values in logarithmic space and the residual realizations.

Figure 6 describes the joint exceedance probability of PSAs based on the above realizations and demonstrates the effects of the multivariate joint distribution of the PSAs. The joint exceedance probability is calculated through equation (24), where n = 5, X₁ to X₅ is PSAs at 0.1 s, 0.2 s, 1 s, 2 s, and 4 s, respectively, and x₁ to x₅ herein indicates a certain level of PSAs at 0.1 s, 0.2 s, 1 s, 2 s, and 4 s, respectively. In this case, the level value is the mean value minus or plus a number times of standard deviation for each PSA investigated. Other two cases are also considered in Figure 6 for comparison purpose by assuming that the residuals of the PSAs are independent without correlation or perfectly dependent. The results imply that the vine copula-based multivariate distribution function proposed can properly characterize the joint distribution of multiple ground motion parameters. The joint exceedance probability of the ground motion parameters are underestimated or overestimated, respectively, if their correlations are ignored or they are assumed to be perfectly correlated. Especially, the difference among three cases become larger in the tail region, that is at a level of large values of ground motions, suggesting that the proposed copula-based multivariate distribution model is necessary to apply in the analysis, especially important in the extreme region.

[figure omitted; refer to PDF]

4. Conclusion

In this study a multivariate joint probability function of PSAs at different vibration periods is calibrated using the vine copula technique. The dependence structure is developed based on a large set of ground motion data consisting of 1550 ground motion records. We show that vine copula can not only better capture the multivariate dependence of PSAs at different vibration periods but also capture their tail dependence which is critical to the losses estimation at low-probability high-impact risks.

In particular, (1) in this study, the vine copula performs better than normal and t copula according to the results of AIC, BIC, and likelihood ration tests; (2) among all the investigated vine copula structures, the best fitted one is a D-vine structure; (3) no residuals of PSAs play a major role but are connected as a path from the shorter period to longer period; and (4) it is observed that the bivariate copulas may show asymmetric tail-dependence property which the normal and t copula could not capture.

The proposed vine copula-based correlation model in this study can be conventionally used in the probabilistic aggregated seismic loss assessment of portfolios or infrastructures.

Acknowledgments

This work was financially supported by the National Natural Science Foundation of China (Grants nos. 51708460 and 71703123) and the Fundamental Research Funds for the Central Universities (2682017CX004 and 2452019117).

References

[1] J. W. Baker, "Vector-valued ground motion intensity measures for probabilistic demand analysis," 2005. Ph.D. dissertation

[2] W. Du, S. Long, C. L. Ning, "An algorithm for selecting spatially correlated ground motions at multiple sites under scenario earthquakes," Journal of Earthquake Engineering,DOI: 10.1080/13632469.2019.1688736, 2019.

[3] G. A. Weatherill, V. Silva, H. Crowley, P. Bazzurro, "Exploring the impact of spatial correlations and uncertainties for portfolio analysis in probabilistic seismic loss estimation," Bulletin of Earthquake Engineering, vol. 13 no. 4, pp. 957-981, DOI: 10.1007/s10518-015-9730-5, 2015.

[4] S. Akkar, Y. Cheng, "Application of a Monte-Carlo simulation approach for the probabilistic assessment of seismic hazard for geographically distributed portfolio," Earthquake Engineering & Structural Dynamics, vol. 45 no. 4, pp. 525-541, DOI: 10.1002/eqe.2667, 2016.

[5] W. Du, "Empirical correlations of frequency-content parameters of ground motions with other intensity measures," Journal of Earthquake Engineering, vol. 23 no. 7, pp. 1073-1091, DOI: 10.1080/13632469.2017.1342303, 2019.

[6] W. Du, G. Wang, "Empirical correlations between the effective number of cycles and other intensity measures of ground motions," Soil Dynamics and Earthquake Engineering, vol. 102, pp. 65-74, DOI: 10.1016/j.soildyn.2017.08.014, 2017.

[7] Y. Cheng, A. Lucchini, F. Mollaioli, "Correlation of elastic input energy equivalent velocity spectral values," Earthquakes and Structures, vol. 8 no. 5, pp. 957-976, DOI: 10.12989/eas.2015.8.5.957, 2015.

[8] G. P. Cimellaro, "Correlation in spectral accelerations for earthquakes in Europe," Earthquake Engineering & Structural Dynamics, vol. 42 no. 4, pp. 623-633, DOI: 10.1002/eqe.2248, 2013.

[9] J. W. Baker, N. Jayaram, "Correlation of spectral acceleration values from NGA ground motion models," Earthquake Spectra, vol. 24 no. 1, pp. 299-317, DOI: 10.1193/1.2857544, 2008.

[10] K. W. Campbell, Y. Bozorgnia, "NGA ground motion model for the geometric mean horizontal component of PGA, PGV, PGD and 5% damped linear elastic response spectra for periods ranging from 0.01 to 10 s," Earthquake Spectra, vol. 24 no. 1, pp. 139-171, DOI: 10.1193/1.2857546, 2008.

[11] W. Du, "An empirical model for the mean period (Tm) of ground motions using the NGA-West2 database," Bulletin of Earthquake Engineering, vol. 15 no. 7, pp. 2673-2693, DOI: 10.1007/s10518-017-0088-8, 2017.

[12] Y. Cheng, A. Lucchini, F. Mollaioli, "Proposal of new ground-motion prediction equations for elastic input energy spectra," Earthquakes and Structures, vol. 7 no. 4, pp. 485-510, DOI: 10.12989/eas.2014.7.4.485, 2014.

[13] Y. Cheng, A. Lucchini, F. Mollaioli, "Ground-motion prediction equations for constant-strength and constant-ductility input energy spectra," Bulletin of Earthquake Engineering, vol. 18 no. 1, pp. 37-55, DOI: 10.1007/s10518-019-00725-x, 2019.

[14] J. Zhang, H. W. Huang, C. H. Juang, W. W. Su, "Geotechnical reliability analysis with limited data: consideration of model selection uncertainty," Engineering Geology, vol. 181, pp. 27-37, DOI: 10.1016/j.enggeo.2014.08.002, 2014.

[15] J. C. Huffman, A. W. Stuedlein, "Reliability-based serviceability limit state design of spread footings on aggregate pier reinforced clay," Journal of Geotechnical and Geoenvironmental Engineering, vol. 140 no. 10,DOI: 10.1061/(asce)gt.1943-5606.0001156, 2014.

[16] X.-S. Tang, D.-Q. Li, Z.-J. Cao, K.-K. Phoon, "Impact of sample size on geotechnical probabilistic model identification," Computers and Geotechnics, vol. 87, pp. 229-240, DOI: 10.1016/j.compgeo.2017.02.019, 2017.

[17] X. S. Tan, D. Q. Li, C. B. Zhou, K. K. Phoon, "Copula-based approaches for evaluating slope reliability under incomplete probability information," Structural Safety, vol. 52, pp. 90-99, DOI: 10.1016/j.strusafe.2014.09.007, 2015.

[18] X.-S. Tang, D.-Q. Li, G. Rong, K.-K. Phoon, C.-B. Zhou, "Impact of copula selection on geotechnical reliability under incomplete probability information," Computers and Geotechnics, vol. 49, pp. 264-278, DOI: 10.1016/j.compgeo.2012.12.002, 2013.

[19] C. L. Ning, Y. Cheng, X. H. Yu, "A simplified approach to investigate the seismic ductility demand of shear-critical reinforced concrete columns based on experimental calibration," Journal of Earthquake Engineering,DOI: 10.1080/13632469.2019.1605949, 2019.

[20] A. J. McNeil, R. Frey, P. Embrechts, Quantitative Risk Management: Concepts, Techniques and Tools, 2005.

[21] K. Goda, G. M. Atkinson, "Interperiod dependence of ground-motion prediction equations: a copula perspective," Bulletin of the Seismological Society of America, vol. 99 no. 2A, pp. 922-927, DOI: 10.1785/0120080286, 2009.

[22] A. Barban, L. D. Persio, "Multivariate option pricing with pair-copulas," Journal of Probability, vol. 2014,DOI: 10.1155/2014/839204, 2014.

[23] D. Huang, C. Yang, B. Zeng, G. Y. Fu, "A copula-based method for estimating shear strength parameters of rock mass," Mathematical Problems in Engineering, vol. 2014,DOI: 10.1155/2014/693062, 2014.

[24] C. D. Michele, G. Salvadori, G. Passoni, R. Vezzoli, "A multivariate model of sea storms using copulas," Coastal Engineering, vol. 54 no. 10, pp. 734-751, DOI: 10.1016/j.coastaleng.2007.05.007, 2007.

[25] Y. Xu, X.-S. Tang, J. P. Wang, H. Kuo-Chen, "Copula-based joint probability function for PGA and CAV: a case study from Taiwan," Earthquake Engineering & Structural Dynamics, vol. 45 no. 13, pp. 2123-2136, DOI: 10.1002/eqe.2748, 2016.

[26] H. Joe, "Families of m-variate distributions with given margins and m (m-1)/2 bivariate dependence parameters," Lecture Notes-Monograph Series, vol. 28, pp. 120-141, 1996.

[27] M. Fischer, C. Köck, S. Schlüter, F. Weigert, "An empirical analysis of multivariate copula models," Quantitative Finance, vol. 9 no. 7, pp. 839-854, DOI: 10.1080/14697680802595650, 2009.

[28] T. Bedford, R. M. Cooke, "Probability density decomposition for conditionally dependent random variables modeled by vines," Annals of Mathematics and Artificial Intelligence, vol. 32 no. 1-4, pp. 245-268, 2001.

[29] T. Bedford, R. M. Cooke, "Vines--a new graphical model for dependent random variables," The Annals of Statistics, vol. 30 no. 4, pp. 1031-1068, DOI: 10.1214/aos/1031689016, 2002.

[30] M. Sklar, Fonctions de Repartition an Dimensions et leurs Marges, 1959.

[31] K. Aas, C. Czado, A. Frigessi, H. Bakken, "Pair-copula constructions of multiple dependence," Insurance: Mathematics and Economics, vol. 44 no. 2, pp. 182-198, DOI: 10.1016/j.insmatheco.2007.02.001, 2009.

[32] D. Allen, M. McAleer, A. Singh, "Risk measurement and risk modelling using applications of vine copulas," Sustainability, vol. 9 no. 10,DOI: 10.3390/su9101762, 2017.

[33] H. Ji, H. Wang, B. Liseo, "Portfolio diversification strategy via tail-dependence clustering and Arma-Garch vine copula approach," Australian Economic Papers, vol. 57 no. 3, pp. 265-283, DOI: 10.1111/1467-8454.12126, 2018.

[34] H. Ji, H. Wang, J. Xu, B. Liseo, "Dependence structure between China’s stock market and other major stock markets before and after the 2008 financial crisis," Emerging Markets Finance and Trade,DOI: 10.1080/1540496X.2019.1615434", 2019.

[35] D. Berg, K. Aas, "Models for construction of multivariate dependence: a comparison study," The European Journal of Finance, vol. 15 no. 7-8, pp. 639-659, DOI: 10.1080/13518470802697428, 2009.

[36] C. Czado, Pair-copula Constructions of Multivariate Copulas. Copula Theory and its Applications, pp. 93-109, 2010.

[37] G. G. Venter, "Tails of copulas," Proceedings of the Casualty Actuarial Society, vol. 89 no. 171, pp. 68-113, 2002.

[38] H. Joe, Multivariate Models and Multivariate Dependence Concepts, 1997.

[39] R. B. Nelsen, An Introduction to Copulas, 2007.

[40] H. Joe, T. Hu, "Multivariate distributions from mixtures of max-infinitely divisible distributions," Journal of Multivariate Analysis, vol. 57 no. 2, pp. 240-265, DOI: 10.1006/jmva.1996.0032, 1996.

[41] H. Joe, H. Li, A. K. Nikoloulopoulos, "Tail dependence functions and vine copulas," Journal of Multivariate Analysis, vol. 101 no. 1, pp. 252-270, DOI: 10.1016/j.jmva.2009.08.002, 2010.

[42] J. A. Tawn, "Bivariate extreme value theory: models and estimation," Biometrika, vol. 75 no. 3, pp. 397-415, DOI: 10.2307/2336591, 1988.

[43] C. Genest, A.-C. Favre, "Everything you always wanted to know about copula modeling but were afraid to ask," Journal of Hydrologic Engineering, vol. 12 no. 4, pp. 347-368, DOI: 10.1061/(asce)1084-0699(2007)12:4(347), 2007.