A Bivariate Extension of Type-II Generalized Crack Distribution for Modeling Heavy-Tailed Losses

Abstract

As an extension of the (univariate) Birnbaum–Saunders distribution, the Type-II generalized crack (GCR2) distribution, built on an appropriate base density, provides a sufficient level of flexibility to fit various distributional shapes, including heavy-tailed ones. In this paper, we develop a bivariate extension of the Type-II generalized crack distribution and study its dependency structure. For practical applications, three specific distributions, GCR2-Generalized Gaussian, GCR2-Student’s t, and GCR2-Logistic, are considered for marginals. The expectation-maximization algorithm is implemented to estimate the parameters in the bivariate GCR2 models. The model fitting results on a catastrophic loss dataset show that the bivariate GCR2 distribution based on the generalized Gaussian density fits the data significantly better than other alternative models, such as the bivariate lognormal distribution and some Archimedean copula models with lognormal or Pareto marginals.

1. Introduction

Researchers have examined diverse classes of distributions to study various facets of problems. Most loss datasets in the context of actuarial loss modeling share some common characteristics, such as being skewed to the right, unimodal (multimodal in certain situations), and having a thin left tail and a moderate to extremely thick right tail. In recent years, various classes of heavy-tailed distributions, including the subexponential distribution class, have been studied for modeling heavy-tailed (or extreme) data [1,2,3]. On the other hand, a different stream of problems, such as the periodic vibrations in commercial aircraft, have motivated the introduction of the Birnbaum–Saunders (BS) distribution [4]. The BS distribution models the total time elapsed until a critical threshold is exceeded by fatigue accumulated on a subject (material) of interest, causing the failure event (or a crack) of the material to occur. Due to its ability to fit right-skewed data, the BS distribution is highly effective for modeling numerous scenarios, e.g., situations where there is an accumulation of a certain factor that drives a quantifiable characteristic to surpass a critical threshold. See [5] for details on the theoretical properties and applications of the BS distribution.

Various extensions of the BS distribution have been discussed in the literature. In [6], an extended version of Birnbaum–Saunders distribution family is introduced using the density of elliptical distribution in place of the standard normal density that quantifies the amount of the stress per cycle of material use in the BS setting. In [7], the extreme value version of the generalized Birnbaum–Saunders (GBS) distribution, whose tail thickness is determined by that of the auxiliary distribution (i.e., an elliptical distribution), has been discussed. Some applications of the extreme value BS models can be found in [5,8].

The (three-parameter) Gaussian crack (lifetime) distribution introduced in [9] is another important extension of the BS distribution. The Gaussian crack (CR) distribution is a two-component mixture of inverse Gaussian and length-biased inverse Gaussian distributions with a weight parameter p, and it features increased flexibility to fit various datasets due to the additional mixture weight parameter. However, the Gaussian crack distribution relies on the standard normal base density, and thus it lacks heavy-tailedness. The limited applicability of the Gaussian crack distribution for modeling heavy-tailed data such as insurance losses motivates the construction of a large class of generalized crack (GCR) distributions [10]. The GCR distribution class contains the Gaussian crack distribution as a specific member, and each member of the class is built on a specific choice of a base-density function that determines the tail characteristics of the resulting GCR distribution. In [10,11], the GCR distributions with the Student’s t and the generalized Gaussian base-density functions are applied to catastrophic losses and heavy-tailed precipitation time series, respectively.

In [12], the GCR distribution class has recently been further extended to the class of Type-II generalized crack (GCR2) distributions in which an additional shape parameter $\tau$ is included to increase flexibility over the GCR class. The key distributional properties, such as the tail characteristics of GCR2 distribution, depend on the specification of the base-density function and the shape parameters involved in each model.

In the literature, several important modeling frameworks allow for the creation of new distribution families from given ones, including the Azzalini method [13], Lehmann-type distributions [14], and Topp-Leone families [15]. While each framework possesses its own distributional characteristics, these general frameworks are common in the sense that the specification of the baseline distribution function (or a density function) plays a crucial role in determining key distributional properties of the constructed model. For instance, the Lehmann-type I setting with Stoppa (baseline) distribution function renders a model that can be effectively used for actuarial data with extreme observations [16]. In this sense, the modeling approach used in the construction of the GCR/GCR2 distribution class is in line with these general frameworks. Typical parametric distributions built under the aforementioned general frameworks have simple (closed) forms, and thus, they may not perform well when data features complex (i.e., multimodal) shapes. Due to its mixture structure, the GCR/GCR2 distribution with an appropriate baseline density can be advantageous over the simple models in such cases.

Regarding applications for bivariate data with heavy-tailed marginals, a bivariate GCR distribution in the form of a four-component mixture of independent models has been constructed in [17]. The authors demonstrate that bivariate GCR models can exhibit useful dependence structures and serve as valuable models for diverse real-world situations.

This paper aims to extend the univariate GCR2 models to bivariate cases by employing the mixture model structure used in [17]. Three specific examples of GCR2 distributions, i.e., a newly constructed GCR2 model based on the logistic density in addition to the GCR2-t and the GCR2-GG models introduced in [12], are used as marginals to effectively model heavy-tailed insurance/catastrophic loss data. We investigate some theoretical properties of the proposed bivariate models, such as the conditional distribution and the dependence structure, and discuss the expectation-maximization (EM) algorithm for model estimation. The applicability of the proposed bivariate GCR2 models and the estimation method is illustrated through a model fitting to a real disaster loss dataset.

The rest of this paper is organized as follows. Section 2 gives a brief review of the origin, definition, examples, and key properties of the univariate Type-II generalized crack distributions for the reader’s convenience. In Section 3, the bivariate Type-II generalized crack distribution is introduced with some detailed discussions of its theoretical properties. A method of model estimation and its application on a real catastrophic loss dataset are presented in Section 4 and Section 5, respectively. Finally, Section 6 provides some concluding remarks.

2. Type-II Generalized Crack Distribution

2.1. Birnbaum–Saunders Distribution

Birnbaum–Saunders distribution (also known as the fatigue life distribution) introduced in [4] is one of the most popular distributions for modeling the time of a crack (or a failure) occurrence when a material specimen is used repeatedly and experiences material fatigue due to the gradual accumulation of stress/damages. The failure event occurs when the accumulated stress on the material specimen hits a critical threshold $w>0$, and the failure time T is the first hitting time of the accumulated stress to the critical threshold. The derivation of the Birnbaum–Saunders distribution as an approximate distribution of the first hitting time of the Brownian motion can be found in [12]. In fact, the BS distribution is a two-point mixture of the inverse Gaussian (IG) and the length-biased inverse Gaussian (LB-IG) distributions with equal mixture weights. Formally, the cumulative distribution function (cdf) and probability density function (pdf) of BS distribution, respectively, are given as follows:

$$\begin{array}{c}\hfill F_{BS}(x;\alpha ,\beta )=\Phi \left({\displaystyle \frac{1}{\alpha }}\left(\sqrt{{\displaystyle \frac{x}{\beta }}}-\sqrt{{\displaystyle \frac{\beta }{x}}}\right)\right),x>0,\end{array}$$

(1)

$$\begin{array}{c}\hfill f_{BS}(x;\alpha ,\beta )={\displaystyle \frac{1}{2}}{f}_{IG}(x;\alpha ,\beta )+{\displaystyle \frac{1}{2}}{f}_{LB-IG}(x;\alpha ,\beta ),x>0,\end{array}$$

(2)

where

$$\begin{array}{ccc}\hfill f_{IG}(x;\alpha ,\beta )& =& {\displaystyle \frac{\sqrt{\beta }}{\alpha }}{x}^{-3/2}\varphi \left[{\displaystyle \frac{1}{\alpha }}\left\{{\left({\displaystyle \frac{x}{\beta }}\right)}^{1/2}-{\left({\displaystyle \frac{\beta }{x}}\right)}^{1/2}\right\}\right]\hfill \end{array}$$

(3)

$$\begin{array}{ccc}\hfill f_{LB-IG}(x;\alpha ,\beta )& =& {\displaystyle \frac{1}{\alpha \sqrt{\beta }}}{x}^{-1/2}\varphi \left[{\displaystyle \frac{1}{\alpha }}\left\{{\left({\displaystyle \frac{x}{\beta }}\right)}^{1/2}-{\left({\displaystyle \frac{\beta }{x}}\right)}^{1/2}\right\}\right],\hfill \end{array}$$

(4)

in which $\beta >0$ and $\alpha >0$ are the scale and shape parameters, respectively; $\Phi (·)$ and $\varphi (·)$ denote the cdf and pdf of the standard normal distribution, and $f_{IG}$ and $f_{LB-IG}$ denote the cdf and pdf of the two-parameter IG and the LB-IG distributions, respectively.

2.2. Gaussian Crack Distribution

The Gaussian crack distribution introduced in [9] is an extension of BS distribution by introducing the weight parameter $p\in [0,1]$ in place of the fixed weight of 1/2 in the BS distribution. Naturally, the CR distribution allows for higher flexibility compared to the classical BS distribution. The pdf of the Gaussian crack distribution is given as follows:

$$\begin{array}{ccc}\hfill f_{CR}(x;\alpha ,\beta ,p)& =& p{f}_{IG}(x;\alpha ,\beta )+q{f}_{LB-IG}(x;\alpha ,\beta )\hfill \\ & =& \left(p{\displaystyle \frac{\sqrt{\beta }}{\alpha }}{x}^{-3/2}+(1-p){\displaystyle \frac{1}{\alpha \sqrt{\beta }}}{x}^{-1/2}\right)\varphi \left[{\displaystyle \frac{1}{\alpha }}\left\{{\left({\displaystyle \frac{x}{\beta }}\right)}^{1/2}-{\left({\displaystyle \frac{\beta }{x}}\right)}^{1/2}\right\}\right],\hfill \end{array}$$

for $x>0$, where $\beta >0$, $\alpha >0$ and $0\le p\le 1$ are scale, shape, and mixture weight parameters, respectively, and $q=1-p$.

Even though the right tail of the CR distribution thickens as p decreases, the tail maintains the shape of an exponential distribution. It is easy to see that the limit of the hazard rate function of the CR distribution converges to a constant that is greater than zero, i.e.,

$$\underset{x\to \infty }{lim}{\displaystyle \frac{f_{CR}\left(x\right)}{1-F_{CR}\left(x\right)}}={\displaystyle \frac{1}{2\beta {\alpha }^2}}>0,$$

which suggests that the Gaussian crack distribution does not belong to heavy-tailed distribution class [3,18].

2.3. Generalized Crack Distribution

In [10], the Gaussian CR distribution family has been extended to a large class of generalized crack (GCR) distribution whose members depend on the specification of the base-density function $g(·)$, where $g(·)$ is a base (or auxiliary) density function which is symmetric about zero. The base-density function replaces the standard normal density function $\varphi (·)$ used for the Gaussian crack density.

Specifically, a random variable T has the generalized crack(GCR) distribution with base density g, denoted as $GCR(\alpha ,\beta ,p;g)$ where the parameters $\alpha >0$, $\beta >0$ and the mixture weight $0\le p\le 1$, if its pdf is given as

$$\begin{array}{c}\hfill f_{GCR}(x;\alpha ,\beta ,p;g)=p{f}_{IS}(x;\alpha ,\beta ;g)+q{f}_{LB-IS}(x;\alpha ,\beta ;g),\end{array}$$

where

$$\begin{array}{ccc}\hfill f_{IS}(x;\alpha ,\beta ;g)& =& {\displaystyle \frac{\sqrt{\beta }}{\alpha }}{x}^{-3/2}g\left[{\displaystyle \frac{1}{\alpha }}\left\{{\left({\displaystyle \frac{x}{\beta }}\right)}^{1/2}-{\left({\displaystyle \frac{\beta }{x}}\right)}^{1/2}\right\}\right]\hfill \\ \hfill f_{LB-IS}(x;\alpha ,\beta ;g)& =& {\displaystyle \frac{1}{\alpha \sqrt{\beta }}}{x}^{-1/2}g\left[{\displaystyle \frac{1}{\alpha }}\left\{{\left({\displaystyle \frac{x}{\beta }}\right)}^{1/2}-{\left({\displaystyle \frac{\beta }{x}}\right)}^{1/2}\right\}\right].\hfill \end{array}$$

The functions $f_{IS}$ and $f_{LB-IS}$ denotes the pdfs of the inverse symmetric (IS) and the length-biased inverse symmetric (LB-IS) distribution, respectively, on the base-density function $g(·)$. The cdf of the $GCR(\alpha ,\beta ,p;g)$ distribution is given as:

$$\begin{array}{c}\hfill F_{GCR}(x;\alpha ,\beta ,p;g)=\overline{G}\left(b\left(x\right)\right)+(2p-1){\int }_{b\left(x\right)}^{\infty }{\displaystyle \frac{s}{\sqrt{s^2+4/{\alpha }^2}}}g\left(s\right)ds\end{array}$$

where $\overline{G}\left(x\right)=1-G\left(x\right)$ is the survival function of the distribution with the cdf $G\left(x\right)={\int }_0^{x}g\left(s\right)ds$ corresponding to the base-density function $g(·)$, and $b\left(x\right)={\displaystyle \frac{1}{\alpha }}\left[\sqrt{\beta /x}-\sqrt{x/\beta }\right]$.

A general expression for the nth raw moments of the IS and LB-IS random variable, along with details on the tail behavior of the GCR distribution, is provided in [10]. Their findings demonstrated that, with an appropriate choice of the base-density function, the heavy-tailed generalized crack distribution performs better than many well-known parametric distributions, such as lognormal, Pareto type II, and Weibull distributions, which are frequently used in modeling positively skewed and heavy-tailed extreme data sets.

2.4. Type-II Birnbaum–Saunders (BS2) Distribution

The Type-II Birnbaum–Saunders distribution extends the BS distribution by introducing another shape parameter $\tau$. Applying the inverse transform method to Equation (1), it can easily be verified that the following stochastic relationship holds. Let us define a random variable T as

$$T=\beta {\left({\displaystyle \frac{\alpha Z+\sqrt{{\alpha }^2{Z}^2+4}}{2}}\right)}^2$$

where Z is a standard normal random variable. The Type-II Birnbaum–Saunders (BS2) distribution includes an additional shape parameter $\tau >0$ to the above expression, defined formally as

$$T=\beta {\left({\displaystyle \frac{\alpha Z+\sqrt{{\alpha }^2{Z}^2+4}}{2}}\right)}^{1/\tau }.$$

The pdf of the Type-II Birnbaum–Saunders (BS2) distribution is

$$\begin{array}{ccc}\hfill f_{BS2}(x;\alpha ,\beta ,\tau )& =& \left({\displaystyle \frac{2\tau }{\alpha \beta }}\right)\left[{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{\beta }{x}}\right)}^{\tau +1}+{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{x}{\beta }}\right)}^{\tau -1}\right]\varphi \left[{\displaystyle \frac{1}{\alpha }}\left\{{\left({\displaystyle \frac{x}{\beta }}\right)}^{\tau }-{\left({\displaystyle \frac{\beta }{x}}\right)}^{\tau }\right\}\right]\hfill \\ & =& {\displaystyle \frac{1}{2}}{f}_{BS2}^{\left(1\right)}(x;\alpha ,\beta ,\tau )+{\displaystyle \frac{1}{2}}{f}_{BS2}^{\left(2\right)}(x;\alpha ,\beta ,\tau ),\hfill \end{array}$$

where $\alpha >0,\beta >0,\tau >0$ are model parameters, and

$$\begin{array}{ccc}\hfill f_{BS2}^{\left(1\right)}(x;\alpha ,\beta ,\tau )& =& \left({\displaystyle \frac{2\tau }{\alpha \beta }}\right){\left({\displaystyle \frac{\beta }{x}}\right)}^{\tau +1}\varphi \left[{\displaystyle \frac{1}{\alpha }}\left\{{\left({\displaystyle \frac{x}{\beta }}\right)}^{\tau }-{\left({\displaystyle \frac{\beta }{x}}\right)}^{\tau }\right\}\right]\hfill \\ \hfill f_{BS2}^{\left(2\right)}(x;\alpha ,\beta ,\tau )& =& \left({\displaystyle \frac{2\tau }{\alpha \beta }}\right){\left({\displaystyle \frac{x}{\beta }}\right)}^{\tau -1}\varphi \left[{\displaystyle \frac{1}{\alpha }}\left\{{\left({\displaystyle \frac{x}{\beta }}\right)}^{\tau }-{\left({\displaystyle \frac{\beta }{x}}\right)}^{\tau }\right\}\right].\hfill \end{array}$$

The cdf of BS2 distribution is

$$\begin{array}{c}\hfill F_{BS2}(x;\alpha ,\beta ,\tau )=\Phi \left[{\displaystyle \frac{1}{\alpha }}\left\{{\left({\displaystyle \frac{x}{\beta }}\right)}^{\tau }-{\left({\displaystyle \frac{\beta }{x}}\right)}^{\tau }\right\}\right].\end{array}$$

(5)

Like the BS distribution, the BS2 distribution is a two-point mixture of densities $f_{BS2}^{\left(1\right)}$ and $f_{BS2}^{\left(2\right)}$ with equal weights. However, with the extra shape parameter $\tau$, the BS2 distribution provides more flexibility over the BS distribution. The BS2 distribution reduces to the Birnbaum–Saunders distribution when $\tau =1/2$.

2.5. Type-II Generalized Crack Distribution

The Type-II generalized crack distribution class introduced in [12] can be seen as a natural extension of the Type-II Birnbaum–Saunders distribution by replacing the standard normal density with a symmetric base density and including the mixture weight parameter p. Specifically, the pdf of GCR2 distribution with base density g is given as follows:

$$\begin{array}{c}\hfill f_{GCR2}(x;\alpha ,\beta ,\tau ,p;g)=p{f}_{GCR2}^{\left(1\right)}(x;\alpha ,\beta ,\tau ;g)+q{f}_{GCR2}^{\left(2\right)}(x;\alpha ,\beta ,\tau ;g)\end{array}$$

where

$$\begin{array}{ccc}\hfill f_{GCR2}^{\left(1\right)}(x;\alpha ,\beta ,\tau ;g)& =& \left({\displaystyle \frac{2\tau }{\alpha \beta }}\right){\left({\displaystyle \frac{\beta }{x}}\right)}^{\tau +1}g\left[{\displaystyle \frac{1}{\alpha }}\left\{{\left({\displaystyle \frac{x}{\beta }}\right)}^{\tau }-{\left({\displaystyle \frac{\beta }{x}}\right)}^{\tau }\right\}\right]\hfill \\ \hfill f_{GCR2}^{\left(2\right)}(x;\alpha ,\beta ,\tau ;g)& =& \left({\displaystyle \frac{2\tau }{\alpha \beta }}\right){\left({\displaystyle \frac{x}{\beta }}\right)}^{\tau -1}g\left[{\displaystyle \frac{1}{\alpha }}\left\{{\left({\displaystyle \frac{x}{\beta }}\right)}^{\tau }-{\left({\displaystyle \frac{\beta }{x}}\right)}^{\tau }\right\}\right],\hfill \end{array}$$

and the cdf is given as

$$\begin{array}{ccc}\hfill F_{GCR2}(x;\alpha ,\beta ,\tau ,p;g)& =& p{F}_{GCR2}^{\left(1\right)}(x;\alpha ,\beta ,\tau ;g)+q{F}_{GCR2}^{\left(2\right)}(x;\alpha ,\beta ,\tau ;g)\hfill \\ & =& \overline{G}\left(b\left(x\right)\right)+(2p-1){\int }_{b\left(x\right)}^{\infty }{\displaystyle \frac{s}{\sqrt{s^2+4/{\alpha }^2}}}g\left(s\right)ds\hfill \end{array}$$

(6)

where

$$\begin{array}{ccc}\hfill F_{GCR2}^{\left(1\right)}(x;\alpha ,\beta ,\tau ;g)& =& \overline{G}\left(b\left(x\right)\right)+{\int }_{b\left(x\right)}^{\infty }{\displaystyle \frac{s}{\sqrt{s^2+4/{\alpha }^2}}}g\left(s\right)ds\hfill \\ \hfill F_{GCR2}^{\left(2\right)}(x;\alpha ,\beta ,\tau ;g)& =& \overline{G}\left(b\left(x\right)\right)-{\int }_{b\left(x\right)}^{\infty }{\displaystyle \frac{s}{\sqrt{s^2+4/{\alpha }^2}}}g\left(s\right)ds,\hfill \end{array}$$

and $b\left(x\right):=b(x;\alpha ,\beta ,\tau )={\alpha }^{-1}\{{(\beta /x)}^{\tau }-{(x/\beta )}^{\tau }\}$. By the symmetry of the base density g and using L’Hopital’s rule, it is easy to verify the following asymptotic relation:

$$\begin{array}{ccc}\hfill \underset{x\to \infty }{lim}{\displaystyle \frac{{\overline{F}}_{GCR2}(x;\alpha ,\beta ,\tau ,p;g)}{G\left(b\right(x\left)\right)}}& =& \underset{x\to \infty }{lim}{\displaystyle \frac{g\left(b\left(x\right)\right)b^{\prime }\left(x\right)\left\{1+(1-2p)\frac{-b\left(x\right)}{\sqrt{{\left(b\left(x\right)\right)}^2+4/{\alpha }^2}}\right\}}{g\left(b\left(x\right)\right)b^{\prime }\left(x\right)}}\hfill \\ & =& 2(1-p).\hfill \end{array}$$

Due to the symmetry of the base density g, $G\left(x\right)=\overline{G}(-x)$, and thus, the GCR2 distribution with base density g has the following asymptotic tail:

$${\overline{F}}_{GCR2}(x;\alpha ,\beta ,\tau ,p;g)\sim 2(1-p)G\left(b\left(x\right)\right)=2(1-p)\overline{G}\left\{{\alpha }^{-1}\left[{(x/\beta )}^{\tau }-{(\beta /x)}^{\tau }\right]\right\},$$

as $x\to \infty$. Here we write $a\left(x\right)\sim b\left(x\right)$ as $x\to \infty$ if ${lim}_{x\to \infty }a\left(x\right)/b\left(x\right)=1$. Specifically when $p=0$, ${\overline{F}}_{GCR2}(x;\alpha ,\beta ,\tau ,p;g)\sim 2\overline{G}\left\{{\alpha }^{-1}\left[{(x/\beta )}^{\tau }-{(\beta /x)}^{\tau }\right]\right\}$, and, when $p=1$, the tail of the GCR2 distribution decays to zero faster than $G\left(b\right(x\left)\right)$ as $x\to \infty$.

2.6. Specific Examples of GCR2 Distributions

The specification of the base-density function on which the GCR2 distribution is built may depend on some key distributional features, such as tail characteristics that are required for each specific application. Two practical examples, the GCR2-t and the GCR2-GG distributions whose base densities are the Student’s t and the generalized Gaussian (normal) distributions, respectively, are given in [12]. In this paper, we also consider the GCR2 distribution with the logistic base density, referred to as the GCR2-LG distribution, as another member of the GCR2 distribution class.

Example 1 

(GCR2-t distribution). The Student’s t distribution has a regular varying tail, and its density function is given as

$$g_t(x;\nu )={\displaystyle \frac{\Gamma \left(\frac{\nu +1}{2}\right)}{\sqrt{\nu \pi }\Gamma \left(\frac{\nu }{2}\right)}}{\left(1+{\displaystyle \frac{x^2}{\nu }}\right)}^{-\left({\displaystyle \frac{\nu +1}{2}}\right)},\nu >0.$$

From this the density of the GCR2-t distribution is expressed as

$$\begin{array}{ccc}\hfill f_{\mathit{GCR}2}^t(x;\alpha ,\beta ,\tau ,p;\nu )& =& {\displaystyle \frac{2\tau \Gamma \left(\frac{\nu +1}{2}\right)}{\alpha \beta \sqrt{\nu \pi }\Gamma \left(\frac{\nu }{2}\right)}}\left[p{\left({\displaystyle \frac{\beta }{x}}\right)}^{\tau +1}+q{\left({\displaystyle \frac{x}{\beta }}\right)}^{\tau -1}\right]\hfill \\ & & \times {\left[1+{\displaystyle \frac{1}{{\alpha }^2\nu }}\{{\left({\displaystyle \frac{x}{\beta }}\right)}^{2\tau }+{\left({\displaystyle \frac{\beta }{x}}\right)}^{2\tau }-2\}\right]}^{-\left({\displaystyle \frac{\nu +1}{2}}\right)}.\hfill \end{array}$$

The Student’s t distribution belongs to the Maximum Domain of Attraction (MDA) of the Fréchet distribution with index ν, which means that the distribution of a properly normalized maximum of independent and identically distributed (i.i.d.) Student’s t random variables converge to a Fréchet distribution asymptotically. In [12], it is shown that the tail of the GCR2-t distribution is regularly varying with index τν, and thus, it also belongs to the MDA of the Fréchet distribution. Compared to the tail of the GCR-t distribution that has the index ${\displaystyle \frac{1}{2}}$ν, the tail of GCR2-t becomes heavier than that of the GCR-t when τ becomes smaller than $1/2$. Figure 1, Figure 2 and Figure 3 illustrate the shapes of the GCR2-t density function for a few sets of prescribed parameter values.

Example 2 

(GCR2-GG distribution). The generalized Gaussian (GG) distribution is a large distribution family that encompasses both thin to moderately heavy-tailed ones, and it can be useful for various practical applications. In particular, when $\theta <1$, the GG distribution has a subexponential tail (see [3] for details on the subexponential distribution class).

The pdf of the generalized Gaussian distribution is given as

$$g_{\mathit{GG}}(x;\theta )={\displaystyle \frac{\theta }{2\lambda \Gamma \left(\frac{1}{\theta }\right)}}exp\left(-{\left|{\displaystyle \frac{x-\mu }{\lambda }}\right|}^{\theta }\right).$$

Recall that the base-density function is required to be symmetric for the construction of the GCR2 model. For this, we set $\mu =0$, and, for the identification of parameters, we further set $\lambda =\sqrt{\Gamma \left(1/\theta \right)/\Gamma \left(3/\theta \right)}$. The pdf of the resulting GCR2-GG distribution is given as follows.

$${\displaystyle f_{\mathit{GCR}2}^{\mathit{GG}}\left(x;\alpha ,\beta ,\tau ,p;\theta \right)=\frac{\tau \theta }{\alpha \beta \lambda \Gamma \left(\frac{1}{\theta }\right)}\left[p{\left(\frac{\beta }{x}\right)}^{\tau +1}+q{\left(\frac{x}{\beta }\right)}^{\tau -1}\right]exp\left(-{\left|\frac{1}{\alpha \lambda }\left[{\left(\frac{x}{\beta }\right)}^{\tau }-{\left(\frac{\beta }{x}\right)}^{\tau }\right]\right|}^{\theta }\right).}$$

Figure 4, Figure 5 and Figure 6 illustrate the shapes of the GCR2-GG density function for a few prescribed sets of parameter values.

Example 3 

(GCR2-LG distribution). The (symmetric) logistic (LG) distribution has the pdf

$$g_{LG}\left(x\right)={\displaystyle \frac{e^{-x/s}}{s{(1+e^{-x/s})}^2}}.$$

For parameter identification, we set $s={\displaystyle \frac{\sqrt{3}}{\pi }}$, and thus, the resulting GCR2-LG distribution has one less parameters than the GCR2-t and GCR2-GG distributions. The pdf of the GCR2-LG distribution is given as

$${\displaystyle f_{\mathit{GCR}2}^{\mathit{LG}}\left(x;\alpha ,\beta ,\tau ,p\right)=\left(\frac{2\pi \tau }{\sqrt{3}\alpha \beta }\right)\left[p{\left(\frac{\beta }{x}\right)}^{\tau +1}+q{\left(\frac{x}{\beta }\right)}^{\tau -1}\right]\frac{exp\left\{-\frac{\pi }{\sqrt{3}\alpha }\left[{\left(\frac{x}{\beta }\right)}^{\tau }-{\left(\frac{\beta }{x}\right)}^{\tau }\right]\right\}}{{\left(1+exp\left\{-\frac{\pi }{\sqrt{3}\alpha }\left[{\left(\frac{x}{\beta }\right)}^{\tau }-{\left(\frac{\beta }{x}\right)}^{\tau }\right]\right\}\right)}^2}.}$$

Due to the fact that the tail of the logistic distribution is heavier than that of normal distribution, the corresponding GCR2-LG distribution can be an effective model for datasets with thin to moderately heavy tails.

Figure 7, Figure 8 and Figure 9 illustrate the shapes of the GCR2-LG density function for a few sets of prescribed parameter values.

As can be seen from Figure 1, Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9, various density shapes can be formed depending on the distributional features of the base-density function and the values of shape parameters. For each specific GCR2 distribution family, the right tail of the distribution becomes heavier as the mixture weight parameter p becomes smaller. The larger the shape parameter $\alpha$ is, the heavier the tail becomes, and the opposite is true for the shape parameter $\tau$. For details of tail properties of the GCR2 distribution, see Theorem 1 and Theorem 2 in [12].

3. Bivariante GCR2 Distribution

In this section, we introduce a bivariate distribution with GCR2 marginals and study its key theoretical properties. Formally, the bivariate Type-II generalized crack (BVGCR2) distribution is defined as follows.

Definition 1. 

A pair of random variables T = ($T_1$, $T_2$) has a bivariate Type-II generalized crack distribution with base density g and parameters $\mathsf{\alpha }$ = (${\alpha }_1$, ${\alpha }_2$), $\mathsf{\beta }$ = (${\beta }_1,{\beta }_2$), $\mathsf{\tau }$ = (${\tau }_1,{\tau }_2$) and $\mathit{p}$ = ($p_{11},p_{12},p_{21},p_{22}$), denoted as BVGCR2($\mathsf{\alpha }$, $\mathsf{\beta }$, $\mathsf{\tau }$, $\mathit{p}$; g), if and only if, its joint pdf is given as follows:

$$\begin{array}{ccc}\hfill f_{BVGCR2}(t_1,t_2;\mathsf{\alpha },\mathsf{\beta },\mathsf{\tau },\mathit{p};g)& =& p_{11}{f}_{GCR2}^{\left(1\right)}(t_1;{\alpha }_1,{\beta }_1,{\tau }_1;g)f_{GCR2}^{\left(1\right)}(t_2;{\alpha }_2,{\beta }_2,{\tau }_2;g)\hfill \\ & +& p_{12}{f}_{GCR2}^{\left(1\right)}(t_1;{\alpha }_1,{\beta }_1,{\tau }_1;g)f_{GCR2}^{\left(2\right)}(t_2;{\alpha }_2,{\beta }_2,{\tau }_2;g)\hfill \\ & +& p_{21}{f}_{GCR2}^{\left(2\right)}(t_1;{\alpha }_1,{\beta }_1,{\tau }_1;g)f_{GCR2}^{\left(1\right)}(t_2;{\alpha }_2,{\beta }_2,{\tau }_2;g)\hfill \\ & +& p_{22}{f}_{GCR2}^{\left(2\right)}(t_1;{\alpha }_1,{\beta }_1,{\tau }_1;g)f_{GCR2}^{\left(2\right)}(t_2;{\alpha }_2,{\beta }_2,{\tau }_2;g),t_1,t_2>0,\hfill \end{array}$$

where the mixture weight parameters satisfy $0\le p_j\le 1$, $j\in J:=\{11,12,21,22\}$, ${\sum }_{j\in J}{p}_j=1$, and, for each $i\in \{1,2\}$,

$$\begin{array}{ccc}\hfill f_{GCR2}^{\left(1\right)}(t_i;{\alpha }_i,{\beta }_i,{\tau }_i;g)& =& \left({\displaystyle \frac{2{\tau }_i}{{\alpha }_i{\beta }_i}}\right){\left({\displaystyle \frac{{\beta }_i}{t_i}}\right)}^{{\tau }_i+1}g\left[{\displaystyle \frac{1}{{\alpha }_i}}\left\{{\left({\displaystyle \frac{t_i}{{\beta }_i}}\right)}^{{\tau }_i}-{\left({\displaystyle \frac{{\beta }_i}{t_i}}\right)}^{{\tau }_i}\right\}\right],\hfill \end{array}$$

and

$$\begin{array}{ccc}\hfill f_{GCR2}^{\left(2\right)}(t_i;{\alpha }_i,{\beta }_i,{\tau }_i;g)& =& \left({\displaystyle \frac{2{\tau }_i}{{\alpha }_i{\beta }_i}}\right){\left({\displaystyle \frac{t_i}{{\beta }_i}}\right)}^{{\tau }_i-1}g\left[{\displaystyle \frac{1}{{\alpha }_i}}\left\{{\left({\displaystyle \frac{t_i}{{\beta }_i}}\right)}^{{\tau }_i}-{\left({\displaystyle \frac{{\beta }_i}{t_i}}\right)}^{{\tau }_i}\right\}\right].\hfill \end{array}$$

Clearly, BVGCR2($\mathsf{\alpha }$, $\mathsf{\beta }$, $\mathsf{\tau }$, $\mathit{p}$; g) is a mixture of four combinations of independent bivariate distributions. It is easy to see that the marginal distributions of $T_1$ and $T_2$ are GCR2 $({\alpha }_1,{\beta }_1,{\tau }_1,p_1=p_{11}+p_{12};g)$ and GCR2 $({\alpha }_2,{\beta }_2,{\tau }_2,p_2=p_{11}+p_{21};g)$, respectively. Please note that for simplicity, we assume that the marginal distributions are built on the same base-density function g, but any model parameters involved in g are distinct for each marginal.

3.1. Conditional Distribution

The conditional distribution is useful in simulating a pair of random variables from the BVGCR2 distribution. From the following relationship between the two mixture components of the GCR2 distribution,

$${\displaystyle \frac{f_{GCR2}^{\left(2\right)}(t;\alpha ,\beta ,\tau ;g)}{f_{GCR2}^{\left(1\right)}(t;\alpha ,\beta ,\tau ;g)}}={\left({\displaystyle \frac{t}{\beta }}\right)}^{2\tau },$$

the conditional density of $T_2$ given $T_1=t_1$, denoted as $f_{2|1}(t_2|t_1;g)$, can be expressed as

$$\begin{array}{cc}\hfill f_{2|1}(t_2|t_1;g)& :={\displaystyle \frac{p_{11}{f}^{\left(1\right)}(t_1;g)f^{\left(1\right)}(t_2;g)+p_{12}{f}^{\left(1\right)}(t_1;g)f^{\left(2\right)}(t_2;g)}{p_1{f}^{\left(1\right)}(t_1;g)+q_1{f}^{\left(2\right)}(t_1;g)}}\hfill \\ \hfill & +{\displaystyle \frac{p_{21}{f}^{\left(2\right)}(t_1;g)f^{\left(1\right)}(t_2;g)+p_{22}{f}^{\left(2\right)}(t_1;g)f^{\left(2\right)}(t_2;g)}{p_1{f}^{\left(1\right)}(t_1;g)+q_1{f}^{\left(2\right)}(t_1;g)}}\hfill \\ \hfill & ={\displaystyle \frac{p_{11}+p_{21}{\left(\frac{t_1}{{\beta }_1}\right)}^{2{\tau }_1}{f}^{\left(1\right)}(t_2;g)+p_{12}+p_{22}{\left(\frac{t_1}{{\beta }_1}\right)}^{2{\tau }_1}{f}^{\left(2\right)}(t_2;g)}{p_1+q_1{\left(\frac{t_1}{{\beta }_1}\right)}^{2{\tau }_1}}}\hfill \\ \hfill & =\left({\displaystyle \frac{p_{11}+p_{21}{\left(\frac{t_1}{{\beta }_1}\right)}^{2{\tau }_1}}{p_1+q_1{\left(\frac{t_1}{{\beta }_1}\right)}^{2{\tau }_1}}}\right)f^{\left(1\right)}(t_2;g)+\left({\displaystyle \frac{p_{12}+p_{22}{\left(\frac{t_1}{{\beta }_1}\right)}^{2{\tau }_1}}{p_1+q_1{\left(\frac{t_1}{{\beta }_1}\right)}^{2{\tau }_1}}}\right)f^{\left(2\right)}(t_2;g).\hfill \end{array}$$

Here the subscripts and parameters in the density functions have been dropped for notational convenience. From the expression, we can see that the conditional distribution of $T_2$ given $T_1=t_1$ is also GCR2 $({\alpha }_2,{\beta }_2,{\tau }_2,p_{2|1};g)$ where

$$p_{2|1}={\displaystyle \frac{p_{11}+p_{21}{(t_1/{\beta }_1)}^{2{\tau }_1}}{p_1+q_1{(t_1/{\beta }_1)}^{2{\tau }_1}}}.$$

Using the conditional distribution, one can easily simulate a pair of random variates $(T_1,T_2)$ from the BVGCR2 model; first simulate $T_1$ from the GCR2 $({\alpha }_1,{\beta }_1,{\tau }_1,p_1;g)$ using the acceptance-rejection method as given in [12], and then simulate, $T_2$ from the conditional distribution GCR2 $({\alpha }_2,{\beta }_2,{\tau }_2,p_{2|1};g)$.

3.2. Dependence Measures

In this section, we derive expressions for Spearman’s rho and Kendall’s tau of BVGCR2 random variables.

3.2.1. Spearman’s Rho

Spearman’s rho is a commonly used measure of dependence between two random variables. Due to the invariance under monotone transformations, Spearman’s rho provides a broad interpretation of the dependence structure for any bivariate distributions. With marginal densities $f_{T_1}\left(t_1\right)$ and $f_{T_2}\left(t_2\right)$ for the random variables $T_1$ and $T_2$, respectively, and the joint distribution $F(t_1,t_2)$ on $(t_1,t_2)\in {\mathbb{R}}_{+}^2$, the (population version) Spearman’s rho is defined as

$$\begin{array}{cc}\hfill {\rho }_s& ={\displaystyle \frac{E(U_1,U_2)-E\left(U_1\right)E\left(U_2\right)}{\sqrt{Var\left(U_1\right)}\sqrt{Var\left(U_2\right)}}}\hfill \\ \hfill & =12{\int }_0^{\infty }{\int }_0^{\infty }F(t_1,t_2)f_{T_1}\left(t_1\right)f_{T_2}\left(t_2\right)d{t}_{1}d{t}_2-3,\hfill \end{array}$$

(7)

where $U_1=F_{T_1}\left(T_1\right)$ and $U_2=F_{T_2}\left(T_2\right)$ are uniform random variables, i.e., Spearman’s rho is Pearson’s correlation between transformations of the original random variables into standard uniform marginals. The following provides an expression for Spearman’s rho of a pair of random variables with a BVGCR2 distribution.

Proposition 1. 

Suppose $(T_1,T_2)\sim$ BVGCR2($\mathsf{\alpha }$,$\mathsf{\beta }$,$\mathsf{\tau }$,$\mathit{p}$;g). Then, Spearman’s rho between $T_1$ and $T_2$ is expressed as

$${\rho }_s=48(p_{11}{p}_{22}-p_{12}{p}_{21}){\gamma }_1{\gamma }_2,$$

(8)

where

$${\gamma }_i={\int }_{-\infty }^{\infty }{\displaystyle \frac{t}{\sqrt{t^2+4/{\alpha }_i^2}}}{G}_i\left(t\right)g_i\left(t\right)dt,i=1,2.$$

Proof. 

See Appendix A. □

3.2.2. Kendall’s Tau

Kendall’s tau is a measure of association (concordance/discordance) between two random variables. Formally, for the random variables $T_1$ and $T_2$ with the joint distribution $F(t_1,t_2)$ on $(t_1,t_2)\in {\mathbb{R}}_{+}^2$, the (population version) Kendall’s tau is defined as

$$\begin{array}{cc}\hfill {\mathsf{\tau }}_k& =P((T_1-T_1^{\prime })(T_2-T_2^{\prime })>0)-P((T_1-T_1^{\prime })(T_2-T_2^{\prime })<0)\hfill \\ \hfill & =4{\int }_0^{\infty }{\int }_0^{\infty }F(t_1,t_2)f(t_1,t_2)d{t}_{1}d{t}_2-1,\hfill \end{array}$$

(9)

where the pair $(T_1^{\prime },T_2^{\prime })$ has the joint distribution F and is independent to $(T_1,T_2)$.

The following provides an expression for Kendall’s tau of random variables with a BVGCR2 distribution.

Proposition 2. 

Suppose $(T_1,T_2)\sim$ BVGCR2($\mathsf{\alpha }$,$\mathsf{\beta }$,$\mathsf{\tau }$,$\mathit{p}$;g). Then, Kendall’s tau between $T_1$ and $T_2$ is expressed as

$${\mathsf{\tau }}_k=32(p_{11}{p}_{22}-p_{12}{p}_{21}){\gamma }_1{\gamma }_2,$$

(10)

where

$${\gamma }_i={\int }_{-\infty }^{\infty }{\displaystyle \frac{t}{\sqrt{t^2+4/{\alpha }_i^2}}}{G}_i\left(t\right)g_i\left(t\right)dt,i=1,2.$$

Proof. 

See Appendix B. □

Remark 1. 

It is important to note that

$$0\le \gamma \le {\displaystyle \frac{1}{4}}.$$

(11)

These inequalities can be proven using the symmetry of $g\left(t\right)$, $G(-t)=1-G\left(t\right)$ and integration by parts. Specifically,

$$\begin{array}{ccc}\hfill \gamma & =& {\int }_{-\infty }^{\infty }{\displaystyle \frac{t}{\sqrt{t^2+4/{\alpha }^2}}}G\left(t\right)g\left(t\right)dt\hfill \\ & =& {\int }_{-\infty }^0{\displaystyle \frac{t}{\sqrt{t^2+4/{\alpha }^2}}}G\left(t\right)g\left(t\right)dt+{\int }_0^{\infty }{\displaystyle \frac{t}{\sqrt{t^2+4/{\alpha }^2}}}G\left(t\right)g\left(t\right)dt\hfill \\ & =& {\int }_0^{\infty }{\displaystyle \frac{-s}{\sqrt{s^2+4/{\alpha }^2}}}G(-s)g(-s)ds+{\int }_0^{\infty }{\displaystyle \frac{t}{\sqrt{t^2+4/{\alpha }^2}}}G\left(t\right)g\left(t\right)dt\hfill \\ & =& -{\int }_0^{\infty }{\displaystyle \frac{s}{\sqrt{s^2+4/{\alpha }^2}}}(1-G\left(s\right))g\left(s\right)ds+{\int }_0^{\infty }{\displaystyle \frac{t}{\sqrt{t^2+4/{\alpha }^2}}}G\left(t\right)g\left(t\right)dt\hfill \\ & =& -{\int }_0^{\infty }{\displaystyle \frac{s}{\sqrt{s^2+4/{\alpha }^2}}}g\left(s\right)ds+2{\int }_0^{\infty }{\displaystyle \frac{s}{\sqrt{s^2+4/{\alpha }^2}}}G\left(s\right)g\left(s\right)ds\hfill \\ & =& {\int }_0^{\infty }(2G\left(s\right)-1){\displaystyle \frac{s}{\sqrt{s^2+4/{\alpha }^2}}}g\left(s\right)ds.\hfill \end{array}$$

Clearly $2G\left(s\right)-1\ge 0$ and $0\le {\displaystyle \frac{s}{\sqrt{s^2+4/{\alpha }^2}}}\le 1$, for $0\le s<\infty$. Then, by the fact ${\int }_0^{\infty }G\left(s\right)g\left(s\right)ds=3/8$, we have

$$\begin{array}{ccc}\hfill 0\le {\int }_0^{\infty }(2G\left(s\right)-1){\displaystyle \frac{s}{\sqrt{s^2+4/{\alpha }^2}}}g\left(s\right)ds& \le & {\int }_0^{\infty }(2G\left(s\right)-1)g\left(s\right)ds=2\left({\displaystyle \frac{3}{8}}\right)-{\displaystyle \frac{1}{2}}={\displaystyle \frac{1}{4}}.\hfill \end{array}$$

From this and by Equations (8) and (10), we obtain the following bounds for ${\rho }_s$ and ${\mathsf{\tau }}_k$:

$$|{\rho }_s|\le {\displaystyle \frac{3}{4}},\left|{\mathsf{\tau }}_k\right|\le {\displaystyle \frac{1}{2}},$$

Note that ${\rho }_s$ and ${\mathsf{\tau }}_k$ become the maximum when $p_{11}$ = $p_{22}=1/2$ and γ attains its maximum of 1/4.

We also remark that if $p_{11}{p}_{22}=p_{12}{p}_{21}$, and thus ${\rho }_s={\mathsf{\tau }}_k=0$, the joint density of the BVGCR2 model can be expressed as a product of two GCR2 marginals, i.e., the two random variables are independent.

3.3. Tail Independence

As remarked in the previous section, the dependency measures of the proposed BVGCR2 model are bounded, and thus, the model may not be suitable for cases where extreme dependency is required, i.e., market turmoil. Here, we further investigate the tail dependence of the BVGCR2 model in terms of the upper-tail dependence, which is defined as

$$\begin{array}{ccc}\hfill {\lambda }_U=\underset{q\to 1}{lim}Pr[T_1>F_{T_1}^{-1}\left(q\right)|T_2>F_{T_2}^{-1}\left(q\right)]& =& \underset{q\to 1}{lim}{\displaystyle \frac{Pr[T_1>F_{T_1}^{-1}\left(q\right),T_2>F_{T_2}^{-1}\left(q\right)]}{1-q}}\hfill \\ & =& \underset{q\to 1}{lim}{\displaystyle \frac{1-2q+F_{T_1,T_2}(F_{T_1}^{-1}\left(q\right),F_{T_2}^{-1}\left(q\right))}{1-q}}\hfill \\ & =& 2-\underset{q\to 1}{lim}{\displaystyle \frac{d}{dq}}{F}_{T_1,T_2}(F_{T_1}^{-1}\left(q\right),F_{T_2}^{-1}\left(q\right)),\hfill \end{array}$$

(12)

where the last equality is held by L’Hopital’s rule. In fact,

$$\begin{array}{ccc}& & {\displaystyle \frac{d}{dq}}{F}_{T_1,T_2}(F_{T_1}^{-1}\left(q\right),F_{T_2}^{-1}\left(q\right))\hfill \\ & & =p_{11}{\displaystyle \frac{d}{dq}}\left\{F_{T_1}^{\left(1\right)}\left(F_{T_1}^{-1}\left(q\right)\right)F_{T_2}^{\left(1\right)}\left(F_{T_2}^{-1}\left(q\right)\right)\right\}+p_{12}{\displaystyle \frac{d}{dq}}\left\{F_{T_1}^{\left(1\right)}\left(F_{T_1}^{-1}\left(q\right)\right)F_{T_2}^{\left(2\right)}\left(F_{T_2}^{-1}\left(q\right)\right)\right\}\hfill \\ & & +p_{21}{\displaystyle \frac{d}{dq}}\left\{F_{T_1}^{\left(2\right)}\left(F_{T_1}^{-1}\left(q\right)\right)F_{T_2}^{\left(1\right)}\left(F_{T_2}^{-1}\left(q\right)\right)\right\}+p_{22}{\displaystyle \frac{d}{dq}}\left\{F_{T_1}^{\left(2\right)}\left(F_{T_1}^{-1}\left(q\right)\right)F_{T_2}^{\left(2\right)}\left(F_{T_2}^{-1}\left(q\right)\right)\right\}\hfill \\ & & =p_{11}\left\{F_{T_2}^{\left(1\right)}\left(F_{T_2}^{-1}\left(q\right)\right){\displaystyle \frac{f_{T_1}^{\left(1\right)}\left(F_{T_1}^{-1}\left(q\right)\right)}{f_{T_1}\left(F_{T_1}^{-1}\left(q\right)\right)}}+F_{T_1}^{\left(1\right)}\left(F_{T_1}^{-1}\left(q\right)\right){\displaystyle \frac{f_{T_2}^{\left(1\right)}\left(F_{T_2}^{-1}\left(q\right)\right)}{f_{T_2}\left(F_{T_2}^{-1}\left(q\right)\right)}}\right\}\hfill \\ & & +p_{12}\left\{F_{T_2}^{\left(2\right)}\left(F_{T_2}^{-1}\left(q\right)\right){\displaystyle \frac{f_{T_1}^{\left(1\right)}\left(F_{T_1}^{-1}\left(q\right)\right)}{f_{T_1}\left(F_{T_1}^{-1}\left(q\right)\right)}}+F_{T_1}^{\left(1\right)}\left(F_{T_1}^{-1}\left(q\right)\right){\displaystyle \frac{f_{T_2}^{\left(2\right)}\left(F_{T_2}^{-1}\left(q\right)\right)}{f_{T_2}\left(F_{T_2}^{-1}\left(q\right)\right)}}\right\}\hfill \\ & & +p_{21}\left\{F_{T_2}^{\left(1\right)}\left(F_{T_2}^{-1}\left(q\right)\right){\displaystyle \frac{f_{T_1}^{\left(2\right)}\left(F_{T_1}^{-1}\left(q\right)\right)}{f_{T_1}\left(F_{T_1}^{-1}\left(q\right)\right)}}+F_{T_1}^{\left(2\right)}\left(F_{T_1}^{-1}\left(q\right)\right){\displaystyle \frac{f_{T_2}^{\left(1\right)}\left(F_{T_2}^{-1}\left(q\right)\right)}{f_{T_2}\left(F_{T_2}^{-1}\left(q\right)\right)}}\right\}\hfill \\ & & +p_{22}\left\{F_{T_2}^{\left(2\right)}\left(F_{T_2}^{-1}\left(q\right)\right){\displaystyle \frac{f_{T_1}^{\left(2\right)}\left(F_{T_1}^{-1}\left(q\right)\right)}{f_{T_1}\left(F_{T_1}^{-1}\left(q\right)\right)}}+F_{T_1}^{\left(2\right)}\left(F_{T_1}^{-1}\left(q\right)\right){\displaystyle \frac{f_{T_2}^{\left(2\right)}\left(F_{T_2}^{-1}\left(q\right)\right)}{f_{T_2}\left(F_{T_2}^{-1}\left(q\right)\right)}}\right\}.\hfill \end{array}$$

Since $f_{T_1}\left(t_1\right)=(p_{11}+p_{12})f^{\left(1\right)}\left(t_1\right)+(p_{21}+p_{22})f^{\left(2\right)}\left(t_1\right)$, $f_{T_2}\left(t_2\right)=(p_{11}+p_{21})f^{\left(1\right)}\left(t_2\right)$ $+(p_{12}+p_{22})f^{\left(2\right)}\left(t_2\right)$, and $f^{\left(2\right)}\left(t_i\right)={(t_i/{\beta }_i)}^{2{\tau }_i},i=1,2$, we have

$$\begin{array}{ccc}\hfill \underset{q\to 1}{lim}{F}_{T_i}^{\left(1\right)}\left(F_{T_i}^{-1}\left(q\right)\right)& =& 1,\underset{q\to 1}{lim}{\displaystyle \frac{f_{T_i}^{\left(1\right)}\left(F_{T_i}^{-1}\left(q\right)\right)}{f_{T_i}\left(F_{T_i}^{-1}\left(q\right)\right)}}=0,i=1,2,\hfill \\ \hfill \underset{q\to 1}{lim}{\displaystyle \frac{f_{T_1}^{\left(2\right)}\left(F_{T_1}^{-1}\left(q\right)\right)}{f_{T_1}\left(F_{T_1}^{-1}\left(q\right)\right)}}& =& {\displaystyle \frac{1}{p_{21}+p_{22}}},\underset{q\to 1}{lim}{\displaystyle \frac{f_{T_2}^{\left(2\right)}\left(F_{T_2}^{-1}\left(q\right)\right)}{f_{T_2}\left(F_{T_2}^{-1}\left(q\right)\right)}}={\displaystyle \frac{1}{p_{12}+p_{22}}}.\hfill \end{array}$$

Therefore,

$$\begin{array}{c}\hfill \underset{q\to 1}{lim}{\displaystyle \frac{d}{dq}}{F}_{T_1,T_2}(F_{T_1}^{-1}\left(q\right),F_{T_2}^{-1}\left(q\right))={\displaystyle \frac{p_{12}}{p_{12}+p_{22}}}+{\displaystyle \frac{p_{21}}{p_{21}+p_{22}}}+{\displaystyle \frac{p_{22}}{p_{21}+p_{22}}}+{\displaystyle \frac{p_{22}}{p_{12}+p_{22}}}=2.\end{array}$$

Then, by Equation (12), we have ${\lambda }_U=2-{lim}_{q\to 1}{\displaystyle \frac{d}{dq}}{F}_{T_1,T_2}(F_{T_1}^{-1}\left(q\right),F_{T_2}^{-1}\left(q\right))=0$, i.e., the BVGCR2 model lacks upper-tail dependence.

4. Model Parameter Estimation

In this section, we discuss parameter estimation for the bivariate Type-II generalized crack distribution using the expectation-maximization algorithm. We briefly review the EM algorithm in a general setting and provide a specific application to the BVGCR2 models.

4.1. Maximum Likelihood Estimation

Suppose $(X_1,X_2,\dots ,X_n)$ is an independent and identically distributed random sample drawn from a density $f\left(x\right|\theta )$. Likelihood function is given as

$$L\left(\theta \right)=f\left(\underset{˜}{x}\right|\theta )=\prod _{i=1}^{n}f\left(x_i\right|\theta ).$$

The maximum likelihood estimation (MLE) aims to find the parameter estimate that maximizes the likelihood function, or equivalently, the log-likelihood function:

$${\widehat{\theta }}_{MLE}=arg\underset{\theta }{max}\{logf\left(\underset{˜}{x}\right|\theta )\}.$$

In the presence of latent (hidden/unobserved) values, however, the direct maximization method often does not provide reliable estimates; hence, other alternative methods are usually considered, and one such method is the expectation-maximization algorithm.

4.2. Expectation-Maximization Algorithm

For simple mixture models with a small number of mixture components, the direct optimization of the log-likelihood function may be used to obtain the maximum likelihood estimates of the model parameters. However, the direct optimization often fails to converge when the number of mixture components is large relative to the sample size.

The expectation-maximization (EM) algorithm [19] is an efficient iterative procedure used to compute the maximum likelihood estimates in the presence of missing or hidden data. When applied to finite mixture model settings, the expectation step renders the separation of the mixture weight parameters from other model parameters for optimization, and the maximization step gives an explicit solution for updating the mixture weights. Due to this, the algorithm is the most widely used maximum likelihood estimation technique for finite mixture models. For details on recent developments of EM-type algorithms for the Poisson mixture model and multivariate Gaussian mixture models in complex data settings, see [20,21], respectively.

Here, we briefly review the general form of the EM algorithm. Suppose we have observed data $\underset{˜}{x}=(x_1,x_2,\dots ,x_n)$ with density $p\left(\underset{˜}{x}\right|\theta )$ and some latent (hidden/unobserved data) $\underset{˜}{z}=(z_1,z_2,\dots ,z_n)$ with density $p\left(\underset{˜}{z}\right|\theta )$. The density of the complete data is denoted by $p(\underset{˜}{x},\underset{˜}{z}|\theta )$. The goal of the EM algorithm is to find the MLE, i.e., the maximum of the observed data likelihood function,

$$L\left(\theta \right)=\sum _{\underset{˜}{z}}p(\underset{˜}{x},\underset{˜}{z}|\theta ).$$

The EM algorithm proceeds by iterating between the following two steps;

• E-Step: This step calculates the expectation of the likelihood with respect to the conditional distribution of $\underset{˜}{z}$ given $\underset{˜}{x}$ and the initial parameter estimate ${\theta }^{\left(m\right)}$, i.e.,

  $$Q\left(\theta \right|{\theta }^{\left(m\right)})=E_{Z|\underset{˜}{x},{\theta }^{\left(m\right)}}\left[logp(\underset{˜}{x},Z|\theta )\right].$$
• M-Step: Choose ${\widehat{\theta }}^{(m+1)}=argma{x}_{\theta }Q\left(\theta \right|{\theta }^{\left(m\right)}).$

Lemma 1. 

The EM algorithm improves $Q\left(\theta \right|{\theta }^{\left(m\right)})$. That is, if $Q\left({\theta }^{(m+1)}\right|{\theta }^{\left(m\right)})\ge Q\left({\theta }^{\left(m\right)}\right|{\theta }^{\left(m\right)})$, then $l\left({\theta }^{(m+1)}\right)\ge l{(\theta }^{\left(m\right)})$.

Proof. 

See Appendix C. □

We now provide the EM algorithm for the estimation of the parameters in the BVGCR2 model. Let $\mathit{t}=(t_1,t_2,\dots ,t_n)$ be a random sample drawn from a BVGCR2($\mathsf{\alpha },\mathsf{\beta },\mathsf{\tau },\mathit{p}$; g) distribution, where $t_i=(t_{i1},t_{i2})$ is a pair of observations for each $i=1,\dots ,n$, and the base density g may have its own parameter(s) such as $\nu$ involved in the Student’s t density. We denote the vector of parameters involved in the base densities by $\mathsf{\theta }=({\theta }_1,{\theta }_2)$. Letting $\mathsf{\gamma }=(\mathsf{\alpha },\mathsf{\beta },\mathsf{\tau },\mathit{p},\mathsf{\theta })$ be the vector of all model parameters, the likelihood function based on the incomplete data is

$$L\left(\mathsf{\gamma }|\mathit{t}\right)=\prod _{i=1}^n\sum _{j\in J}{q}_j{f}_j(t_{i1},t_{i2};\mathsf{\alpha },\mathsf{\beta },\mathsf{\tau };\mathsf{\theta };g)$$

where $J=\{11,12,21,22\}$ is the set of indexes and

$$f_j(t_i;\mathsf{\gamma };g)=\left\{\begin{array}{cc}{f}_{GCR2}^{\left(1\right)}(t_{i1};{\alpha }_1,{\beta }_1,{\tau }_1;{\theta }_1;g)f_{GCR2}^{\left(1\right)}(t_{i2};{\alpha }_2,{\beta }_2,{\tau }_2;{\theta }_2;g),\hfill & j=11\hfill \\ f_{GCR2}^{\left(1\right)}(t_{i1};{\alpha }_1,{\beta }_1,{\tau }_1;{\theta }_1;g)f_{GCR2}^{\left(2\right)}(t_{i2};{\alpha }_2,{\beta }_2,{\tau }_2;{\theta }_2;g),\hfill & j=12\hfill \\ f_{GCR2}^{\left(2\right)}(t_{i1};{\alpha }_1,{\beta }_1,{\tau }_1;{\theta }_1;g)f_{GCR2}^{\left(1\right)}(t_{i2};{\alpha }_2,{\beta }_2,{\tau }_2;{\theta }_2;g),\hfill & j=21\hfill \\ f_{GCR2}^{\left(2\right)}(t_{i1};{\alpha }_1,{\beta }_1,{\tau }_1;{\theta }_1;g)f_{GCR2}^{\left(2\right)}(t_{i2};{\alpha }_2,{\beta }_2,{\tau }_2;{\theta }_2;g),\hfill & j=22.\hfill \end{array}\right.$$

On the other hand, letting $\mathit{Z}=(Z_1,Z_2,\dots ,Z_n)$ be a set of latent variables where Pr$(Z_i=j)=q_j,j=1,2,$ for each $i\in \{1,\dots ,n\}$, the likelihood function based on the (augmented) complete data is

$$L\left(\mathsf{\gamma }|\mathit{t},\mathit{Z}\right)=\prod _{i=1}^n\prod _j{\left[q_j{f}_j(t_{i1},t_{i2};\mathsf{\alpha },\mathsf{\beta },\mathsf{\tau };\mathsf{\theta };g)\right]}^{I(Z_i=j)}$$

where $I(·)$ denotes the indicator function.

Let ${\mathsf{\gamma }}^{\left(m\right)}=({\mathsf{\alpha }}^{\left(m\right)},{\mathsf{\beta }}^{\left(m\right)},{\mathsf{\tau }}^{\left(m\right)},{\mathit{p}}^{\left(m\right)},{\mathsf{\theta }}^{\left(m\right)})$ denote the current estimate of $\mathsf{\gamma }$ after m-th iteration of the EM algorithm. Then, by Bayes’ theorem, we have

$$p_{ij}^{\left(m\right)}:=\mathrm{Pr}(Z_i=j|t_i,{\mathsf{\gamma }}^{\left(\mathit{m}\right)})={\displaystyle \frac{q_j^{\left(m\right)}{f}_j(t_{i1},t_{i2};{\mathsf{\alpha }}^{\left(\mathit{m}\right)},{\mathsf{\beta }}^{\left(\mathit{m}\right)},{\mathsf{\tau }}^{\left(\mathit{m}\right)};{\mathsf{\theta }}^{\left(\mathit{m}\right)})}{{\sum }_{j\in J}{q}_j^{\left(m\right)}{f}_j(t_{i1},t_{i2};{\mathsf{\alpha }}^{\left(\mathit{m}\right)},{\mathsf{\beta }}^{\left(\mathit{m}\right)},{\mathsf{\tau }}^{\left(\mathit{m}\right)};{\mathsf{\theta }}^{\left(\mathit{m}\right)})}},$$

The expectation step (E-step) of the EM algorithm follows.

$$\begin{array}{ccc}\hfill E_{Z|t,{\mathsf{\gamma }}^{\left(m\right)}}\left[logL\left(\mathsf{\gamma }|\mathit{t},\mathit{Z}\right)\right]& =& \sum _{i=1}^n{E}_{Z|t,{\mathsf{\gamma }}^{\left(\mathit{m}\right)}}(log\prod _j{\left[q_j{f}_j(t_{i1},t_{i2};\mathsf{\alpha },\mathsf{\beta },\mathsf{\tau };\mathsf{\theta };g)\right]}^{I(Z_i=j)})\hfill \\ & =& \sum _{i=1}^n{E}_{Z|t,{\mathsf{\gamma }}^{\left(\mathit{m}\right)}}(\sum _{j\in J}I(Z_i=j)log\left[q_j{f}_j(t_{i1},t_{i2};\mathsf{\alpha },\mathsf{\beta },\mathsf{\tau };\mathsf{\theta };g)\right])\hfill \\ & =& \sum _{i=1}^n\sum _{j\in J}\mathrm{Pr}(Z_i=j|t_i,{\mathsf{\gamma }}^{\left(\mathit{m}\right)})log\left[q_j{f}_j(t_{i1},t_{i2};\mathsf{\alpha },\mathsf{\beta },\mathsf{\tau };\mathsf{\theta };g)\right]\hfill \\ & =& \sum _{i=1}^n\sum _{j\in J}{p}_{ij}^{\left(m\right)}log\left[q_j{f}_j(t_{i1},t_{i2};\mathsf{\alpha },\mathsf{\beta },\mathsf{\tau };\mathsf{\theta };g)\right]\hfill \\ & =& \sum _{i=1}^n\sum _{j\in j}{p}_{ij}^{\left(m\right)}log{q}_j+\sum _{i=1}^n\sum _{j\in J}{p}_{ij}^{\left(m\right)}log{f}_j(t_{i1},t_{i2};\mathsf{\alpha },\mathsf{\beta },\mathsf{\tau };\mathsf{\theta };g).\hfill \end{array}$$

(13)

The maximization step (M-step) finds the updated parameter estimates that maximize the objective function (13), which separates ${\mathit{q}}_{\mathit{j}}$ and the other parameters ($\mathsf{\alpha },\mathsf{\beta },\mathsf{\tau };\mathsf{\theta }$). The update of ${\mathit{q}}_{\mathit{j}}$ can be dealt with separately by applying the method of Lagrange multiplier. The updated estimate is

$${\mathit{q}}_{\mathit{j}}^{(\mathit{m}+\mathbf{1})}=\underset{q}{arg\; max}\left\{\sum _{i=1}^n\sum _{j\in J}{p}_{ij}^{\left(m\right)}log{q}_j\right\}={\left({\displaystyle \frac{1}{n}}\sum _{i=1}^n{p}_{i,11}^{\left(m\right)},\dots ,{\displaystyle \frac{1}{n}}\sum _{i=1}^n{p}_{i,22}^{\left(m\right)}\right)}^{\prime }.$$

The updated estimates ${\mathsf{\alpha }}^{(\mathit{m}+\mathbf{1})},{\mathsf{\beta }}^{(\mathit{m}+\mathbf{1})},{\mathsf{\tau }}^{(\mathit{m}+\mathbf{1})}$ and ${\mathsf{\theta }}^{(\mathit{m}+\mathbf{1})}$ are the maximizers of the objective function

$$\sum _{i=1}^n\sum _{j\in J}{p}_{ij}^{\left(m\right)}log{f}_j(t_{i1},t_{i2};\mathsf{\alpha },\mathsf{\beta },\mathsf{\tau };\mathsf{\theta };g)=\left(Q_1\right)+\left(Q_2\right),$$

where

$$\begin{array}{cc}\hfill \left(Q_1\right)=\sum _{i=1}^n& [(p_{i,11}^{\left(m\right)}+p_{i,12}^{\left(m\right)})log{f}_{GCR2}^{\left(1\right)}(t_{i1};{\alpha }_1,{\beta }_1,{\tau }_1;{\theta }_1;g)\hfill \\ \hfill & +(p_{i,21}^{\left(m\right)}+p_{i,22}^{\left(m\right)})log{f}_{GCR2}^{\left(2\right)}(t_{i1};{\alpha }_1,{\beta }_1,{\tau }_1;{\theta }_1;g)],\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill \left(Q_2\right)=\sum _{i=1}^n& [(p_{i,11}^{\left(m\right)}+p_{i,21}^{\left(m\right)})log{f}_{GCR2}^{\left(1\right)}(t_{i2};{\alpha }_2,{\beta }_2,{\tau }_2;{\theta }_2;g)\hfill \\ \hfill & +(p_{i,12}^{\left(m\right)}+p_{i,22}^{\left(m\right)})log{f}_{GCR2}^{\left(2\right)}(t_{i2};{\alpha }_2,{\beta }_2,{\tau }_2;{\theta }_2;g)].\hfill \end{array}$$

That is, the maximization for $({\alpha }_1,{\beta }_1,{\tau }_1;{\theta }_1)$ can proceed separately from that for $({\alpha }_2,{\beta }_2,{\tau }_2;{\theta }_2)$ and thus, the dimensionality of the optimization problem is reduced significantly.

5. Applications

In [12], the usefulness of the univariate GCR2 models for heavy-tailed data modeling has been demonstrated through an application with a real loss dataset. In this section, we fit several bivariate Type-II generalized crack distributions on a real catastrophic loss dataset compiled from the International Disaster Database (EM-DAT, www.emdat.be, accessed on 22 November 2024). Specifically, each observation in the dataset is composed of two variables: ‘Meteo’ and ‘Hydro’. Marginally, ‘Meteo’ is a quarterly time series of (estimated) losses from meteorological disasters such as storms and extreme temperatures, spanning from 1950 to 2022 in Asia, and ‘Hydro’ is a series of (estimated) losses due to hydrological disasters such as flood and landslide for the same geographical area and the time span. For bivariate model fitting, we remove the pairs with missing observations, resulting in a bivariate dataset with 166 observations. The losses are inflation-adjusted to be equivalent to the US dollar values in 2021.

Table 1. Descriptive summary statistics of the dataset (Unit: 100 million USD).

Data	n	Min	1st Qu.	Median	Mean	3rd Qu.	Max	Skewness	Kurtosis
Meteo	166	0.003	2.005	7.275	30.799	35.282	287.090	2.820	11.304
Hydro	166	0.007	1.692	12.172	45.538	53.330	621.207	4.014	23.188

presents summary statistics of losses due to meteorological and hydrological disasters in Asia. The descriptive statistics, along with the histogram and the normal Q-Q plots (Figure 10 and Figure 11), suggest that the marginal distributions of the two variables are both positively skewed and heavy-tailed. Many time series include some deterministic components such as long-term trends and seasonality. To isolate the deterministic components, each quarterly time series is decomposed under the multiplicative model assumption. The time series decompositions given in Figure 12 show the presence of strong seasonality in both time series and weak evidence of long-term trend. Since the proposed GCR2 models do not assume any deterministic seasonality, we deseasonalize each time series by dividing it by its estimated seasonal component.

Figure 13 gives the scatter plots of the two seasonally adjusted variables and their log-transformations. The figure shows some evidence of a dependent relationship between the two variables. The sample Spearman’s rho and Kendall’s tau of the variables are 0.425 and 0.291, respectively.

We apply the EM algorithm described in Section 4.2 to fit the (seasonally adjusted) data using six specific bivariate models: GCR-GG, GCR-t, GCR-LG, GCR2-GG, GCR2-t, and GCR2-LG. For each model fitting, the EM algorithm requires initialization of the parameter values. To obtain a reliable result, we first fit the marginal models separately using the estimation method given in [12], and the fitted values of the marginal parameters, i.e., $({\alpha }_i,{\beta }_i,{\tau }_i,\theta )$ for GCR models and $({\alpha }_i,{\beta }_i,{\tau }_i,\theta ),i=1,2$, for the GCR2 models, are used for the initialization of the BVGCR (or BVGCR2) model parameters. For mixture weight parameter initialization, we use a large set of possible parameter values satisfying $p_{11}+p_{12}=p_1$ and $p_{11}+p_{21}=p_2$, where $p_1$ and $p_2$ are the mixture weight parameter estimates for the marginal models. The log-likelihood function values of the EM fits under the set of initializations are compared, and the one that gives the largest log-likelihood value is selected for the final fit. To compare the performance of BVGCR2 models with some other benchmark bivariate models commonly used in loss modeling, we further implement the maximum likelihood estimation of the following seven bivariate models: BVLNorm (bivariate lognormal), Clayton-LNorm (bivariate Clayton copula with lognormal marginals), Gumbel-LNorm (bivariate Gumbel copula with lognormal marginals), Frank-LNorm (bivariate Frank copula with lognormal marginals), Clayton-Pareto (bivariate Clayton copula with two-parameter Pareto marginals), Gumbel-Pareto (bivariate Gumbel copula with two-parameter Pareto marginals), and Frank-Pareto (bivariate Frank copula with two-parameter Pareto marginals) models.

Since the number of estimated model parameters differs by model, we compare the fits of candidate models in terms of the Akaike information criterion (AIC) and the Bayesian information criterion (BIC), defined as

$$\begin{array}{ccc}\hfill AIC& =& -2\left(\mathrm{Log}\text{-}\mathrm{Likelihood}\right)+2k\hfill \\ \hfill BIC& =& -2\left(\mathrm{Log}\text{-}\mathrm{Likelihood}\right)+klogn,\hfill \end{array}$$

respectively, where k is the number of estimated parameters and n is the sample size. Among candidate models, the preferred model is the one with the smallest value of either of these criteria.

From

Table 2. Model comparison with AIC and BIC.

Model	Log-Likelihood	AIC	BIC
BVGCR-GG	−1370.332	2758.664	2786.672
BVGCR-t	−1372.381	2762.762	2790.770
BVGCR-LG	−1372.066	2758.132	2779.916
BVGCR2-GG	−1367.346	2756.692	2790.924
BVGCR2-t	−1371.782	2765.564	2799.796
BVGCR2-LG	−1371.707	2761.414	2789.422
BVLNorm	−1387.839	2785.678	2796.126
Clayton-LNorm	−1393.159	2796.318	2811.878
Gumbel-LNorm	−1385.676	2781.352	2796.912
Frank-LNorm	−1386.976	2783.952	2799.512
Clayton-Pareto	−1405.479	2820.958	2836.518
Gumbel-Pareto	−1403.067	2816.134	2831.694
Frank-Pareto	−1405.029	2820.058	2835.618

, we see that, based on the Akaike information criterion, the fitted BVGCR2-GG model outperforms all the other alternative bivariate models. However, the BVGCR-LG model is preferred in terms of the Bayesian information criterion, which heavily penalizes complex models. Comparing the BVGCR-GG to the BVGCR2-GG, one can see that the BVGCR2-GG significantly improves model fitting due to the additional shape parameter.

Table 3. Parameter estimates of marginals and Bivariate GCR2-GG models.

	${\widehat{\mathit{\alpha }}}_1$	${\widehat{\mathit{\beta }}}_1$	${\widehat{\mathit{\tau }}}_1$	${\widehat{\mathit{\theta }}}_1$	$\widehat{{\mathit{p}}_1}$	${\widehat{\mathit{\alpha }}}_2$	${\widehat{\mathit{\beta }}}_2$	${\widehat{\mathit{\tau }}}_2$	${\widehat{\mathit{\theta }}}_2$	$\widehat{{\mathit{p}}_2}$	${\widehat{\mathit{p}}}_{11}$	${\widehat{\mathit{p}}}_{12}$	${\widehat{\mathit{p}}}_{21}$	${\widehat{\mathit{p}}}_{22}$
Meteo	20.075	0.915	0.737	0.623	0.094
Hydro						11.152	1.964	0.669	0.809	0.189
Bivariate	22.644	1.494	0.791	0.477		11.442	1.964	0.675	0.798		0.124	0.040	0.072	0.764

gives parameter estimates of the fitted marginal distributions and the bivariate GCR2-GG model.

Please note that the values of $(p_{11}+p_{12})$ and $(p_{11}+p_{21})$ based on the fitted BVGCR2-GG model are larger than the corresponding estimates of $p_1$ and $p_2$. Spearman’s rho and Kendall’s tau under the fitted BVGCR2-GG model are 0.173 and 0.115, respectively, which are lower than the empirical counterparts. This may be because the empirical data contains some spurious dependence due to the deterministic trends.

To test the statistical significance of the parameter estimates of the BVGCR2-GG model, we construct the bootstrap confidence intervals by computing the 2.5th and the 97.5th percentiles of the estimates based on the 100 bootstrap samples.

Table 4. The 95% bootstrap confidence intervals for the parameters in Bivariate GCR2-GG model.

	${\widehat{\mathit{\alpha }}}_1$	${\widehat{\mathit{\beta }}}_1$	${\widehat{\mathit{\tau }}}_1$	${\widehat{\mathit{\theta }}}_1$	${\widehat{\mathit{\alpha }}}_2$	${\widehat{\mathit{\beta }}}_2$	${\widehat{\mathit{\tau }}}_2$	${\widehat{\mathit{\theta }}}_2$	${\widehat{\mathit{p}}}_{11}$	${\widehat{\mathit{p}}}_{12}$	${\widehat{\mathit{p}}}_{21}$	${\widehat{\mathit{p}}}_{22}$
Original estimate	22.644	1.494	0.791	0.477	11.442	1.964	0.675	0.798	0.124	0.040	0.072	0.764
Bootstrap mean	27.673	1.524	0.735	0.583	14.189	2.065	0.680	0.988	0.131	0.037	0.067	0.766
2.5th percentile	7.031	0.762	0.542	0.347	2.356	1.500	0.322	0.418	0.059	0.000	0.000	0.687
97.5th percentile	52.149	2.252	0.951	1.013	42.898	2.803	0.913	3.027	0.192	0.080	0.147	0.843

shows that all the parameter estimates on the original dataset fall in the corresponding 95% bootstrap confidence intervals.

Figure 14 gives the contour plots of the fitted BVGCR2-GG model, and that of the log-transformed random variables, and Figure 15 presents the scatter plots of the simulated random variables from the BVGCR2-GG model and their log-transformations. Comparing these plots with Figure 13, we can see that the fitted model explains the dependence structure of the empirical data well.

6. Concluding Remarks

In this paper, we constructed a bivariate extension of the Type-II generalized crack distribution and studied a few specific examples of the bivariate GCR2 distributions base on the generalized Gaussian, Student’s t, and logistic densities to demonstrate the applicability of the constructed model. Specifically, our main theoretical finding is that the level of dependence of the constructed BVGCR2 model in terms of Kendall’s tau and Spearman’s rho is a weak to medium association. The model fitting to catastrophic loss data showed that the fitted BVGCR2-GG model outperformed all the other alternative models based on the Akaike information criterion. Especially when compared to the BVGCR-GG model, the BVGCR2-GG model has shown a significant improvement due to the increased flexibility. With an appropriate choice of base-density function, the proposed BVGCR2 model can be effectively used for various applications that require a weak to moderate level of dependence.

With the lack of the upper-tail dependence, the bivariate GCR2 distributions may not be applicable for the situations where variables are (or assumed to be) asymptotically dependent in the upper tail, e.g., stress testing for market/credit portfolios. This limitation can be alleviated using a common parameter for both marginals and randomizing the parameter. For example, we may set $\tau ={\tau }_1={\tau }_2$, the shape parameters in the GCR2 marginals, and assume $\tau$ follows a Gamma distribution. The use of a common random parameter is expected to widen the range of dependence levels and allow for upper-tail dependence. We will pursue this approach in future research.