On Modeling the Earthquake Insurance Data via a New Member of the T-X Family

Abstract

Heavy-tailed distributions play an important role in modeling data in actuarial and financial sciences. In this article, a new method is suggested to define new distributions suitable for modeling data with a heavy right tail. The proposed method may be named as the Z-family of distributions. For illustrative purposes, a special submodel of the proposed family, called the Z-Weibull distribution, is considered in detail to model data with a heavy right tail. The method of maximum likelihood estimation is adopted to estimate the model parameters. A brief Monte Carlo simulation study for evaluating the maximum likelihood estimators is done. Furthermore, some actuarial measures such as value at risk and tail value at risk are calculated. A simulation study based on these actuarial measures is also done. An application of the Z-Weibull model to the earthquake insurance data is presented. Based on the analyses, we observed that the proposed distribution can be used quite effectively in modeling heavy-tailed data in insurance sciences and other related fields. Finally, Bayesian analysis and performance of Gibbs sampling for the earthquake data have also been carried out.

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 a number of applied areas such as finance and actuarial sciences, data sets are most often positive, and the respective distribution is unimodal hump-shaped and skewed to right having heavier tails as compared to the well-known classical distributions. These distributions are not much flexible to adequately model such types of heavy-tailed data sets. For example, (i) the Pareto distribution, which is frequently used to model financial data sets, does not provide a reasonable fit for many applications, for example, if we are interested in modeling only especially moderate-to-large losses altogether, then in such cases, the Pareto distribution may not be a suitable choice to use [1], and (ii) the Weibull model is capable of catering the behavior of small losses very closely, but, unfortunately, fails to provide an adequate fit to the large losses [2]. In such circumstances, the utilization of the heavy-tailed models may be a good choice to apply. For positive data, heavy-tailed distributions are those whose right-tail probabilities are greater than the exponential one [3], that is,$$\begin{array}{}\text{(1)}& \underset{x⟶\infty }{\text{lim}}{e}^{\text{px}}\left(1-F\left(x;\xi \right)\right)=\infty ,\hspace{1em}p>0,\end{array}$$where $F\left(x;\xi \right)$ is the cumulative distribution function (cdf) depending on the parameter vector $\xi \in ℝ$.

Due to the usefulness and flexibility of the heavy-tailed models in financial and actuarial practice, actuaries are always intended to propose new statistical distributions. Therefore, serious attempts have been made to propose new statistical models and are still growing rapidly. The new contribution is made via different approaches such as (i) transformation of variables, (ii) composition of two or more distributions, (iii) compounding of distributions, and (iv) finite mixture of distributions, see Ahmad et al. [4].

Recent investigation of Eling [5] and Adcock et al. [6] determined that skew-normal and skew Student’s t distributions are the most excellent competitors because the skewed distributions adjust right-skewness and high kurtosis; for the interested readers, one can refer to Shushi [7] and Punzo et al. [8]. However, insurance losses and monetary risks take values on the positive real line, and subsequently, these skew models may not be a suitable choice to use as they are defined on $ℝ$. In such circumstances, the transformation of variable approach, especially, the exponential transformation, has demonstrated to be considerable; for details, see Azzalini et al. [9]. Bagnato and Punzo [10] showed that the transformation method of introducing new distributions is easy to use; however, most often, the inferences become complicated.

Another useful method of proposing new versatile heavy-tailed distributions, which provide the best fit to the heavy-tailed losses, is the methodology of composition, see Paula et al. [11], Klugman et al. [12], Nadarajah and Bakar [13] and Bakar et al. [14]. However, it ought to be noted that the new statistical models introduced via this approach involve more than three parameters inflicting difficulties in estimating the parameters, and more computational efforts are needed.

Another approach of introducing new distributions to cater data modeling adequately with unimodality is compounding of distributions, see Punzo et al. [15] and Mazza and Punzo [16]. Unfortunately, the density function of the distributions obtained by this approach might not have a closed-form expression that makes the estimation more complicated as shown in Punzo et al. [15].

The method of finite mixture models is another prominent approach to obtain new, very flexible models which are able to capture, for example, multimodality of the distribution under consideration, see Bernardi et al. [17], Miljkovic and Grün [18], and Punzo et al. [15]. No doubt, the distributions obtained via this approach are much flexible, but the inferences become more complicated and computationally challenging.

Furthermore, Dutta and Perry [19] performed an empirical analysis of loss distributions, and risk was estimated by different approaches such as exploratory data analysis and other empirical approaches. These authors rejected the idea of using the exponential, gamma, and Weibull models in modeling insurance losses due to the poor results. They concluded that one would need to use a model that is flexible enough in its structure. This motivated the researchers to search for more flexible models offering greater accuracy in fitting the heavy-tailed data.

Hence, bringing flexibility to a model by introducing additional parameter(s) is a desirable feature [20–24]. In a number of recent papers, serious attempts have been made to introduce a variety of new heavy-tailed distributions, see Ahmad et al. [25] and Ahmad et al. [26]. Due to the importance of statistical distributions in financial science, a new family of distributions, called the Z-family, is introduced. The proposed family is introduced via the T-X family approach [27]. To illustrate the usefulness of the proposed method, a three-parameter submodel, called the Z-Weibull distribution, is taken and studied in detail. The proposed distribution provides a better description of the earthquake insurance data with possibly heavy tails than the available (i) two-parameter distributions such as Weibull, Burr-XII (B-XII), and generalized exponential (GE), (ii) three-parameter Weibull-claim (W-claim) and exponentiated Lomax (EL), and (iii) four-parameter beta-Weibull (BW) distributions, and possibly many others.

The rest of the article is structured in the following way. The proposed method is introduced in Section 2. The Z-Weibull model is considered in Section 3, and the shapes of its probability density function (pdf) are investigated in the same section. Estimation of parameters is discussed in Section 4. In the same section, a detailed Monte Carlo simulation study is conducted. Actuarial measures of the proposed method along with a simulation study are provided in Section 5. Distribution fit to the earthquake insurance data set is discussed in Section 6. Bayesian analysis as well as Gibbs sampling procedure for the real data set is discussed in Section 7. Future frame work is discussed in Section 8. Finally, some concluding remarks are presented in Section 9.

2. Development of the Z-Family

Let $p\left(t\right)$ be the pdf of a random variable, say T, where $T\in \left[n_1,n_2\right]$ for $-\infty \le n_1,n_2\le \infty$, and let $W\left[F\left(x;\xi \right)\right]$ be a function of $F\left(x;\xi \right)$ of a random variable, say X, depending on the vector parameter $\xi$ satisfying the conditions given in the following:

(1)

$W\left[F\left(x;\xi \right)\right]\in \left[n_1,n_2\right]$

(2)

$W\left[F\left(x;\xi \right)\right]$ is differentiable and monotonically increasing

(3)

$W\left[F\left(x;\xi \right)\right]⟶n_1$ as $x⟶-\infty$ and $W\left[F\left(x;\xi \right)\right]⟶n_2$ as $x⟶\infty$

Alzaatreh et al. [27] introduced a general method for generating new families of distributions called the T-X family which is defined by$$\begin{array}{}\text{(2)}& G\left(x\right)={\displaystyle {\int }_{n_1}^{W\left[F\left(x;\xi \right)\right]}p\left(t\right)\text{d}t},\hspace{1em}t\ge 0,\end{array}$$where $W\left[F\left(x;\xi \right)\right]$ satisfies the conditions mentioned above. The probability density function (pdf) corresponding to (2) is given by$$\begin{array}{}\text{(3)}& g\left(x\right)=\left\{\frac{d}{\text{d}x}W\left(F\left(x;\xi \right)\right)\right\}p\left(W\left[F\left(x;\xi \right)\right]\right).\end{array}$$

Deploying the T-X proposal, several new classes of distributions have been introduced in the literature [28]. Let X have the exponential distribution with the pdf given by$$\begin{array}{}\text{(4)}& p\left(t\right)=\theta e^{-\theta t},\hspace{1em}t,\theta >0.\end{array}$$

Using $\theta =1$ in (4), we get$$\begin{array}{}\text{(5)}& p\left(t\right)=e^{-t},\hspace{1em}t>0.\end{array}$$

On setting $W\left[F\left(x;\xi \right)\right]=-\log \left\{\left(1-F\left(x;\xi \right)\right)/\left({\beta }^{F\left(x;\xi \right)}\right)\right\}$ and $p\left(t\right)=e^{-t}$ in (2), we define the cdf of the Z-family of distributions by$$\begin{array}{}\text{(6)}& G\left(x;\beta ,\xi \right)=1-\left\{\frac{1-F\left(x;\xi \right)}{{\beta }^{F\left(x;\xi \right)}}\right\},\hspace{1em}\beta >0,x,\xi \in ℝ,\end{array}$$where $F\left(x;\xi \right)$ is the baseline cdf. The expression in (6) represents a wide family of univariate continuous distributions. Clearly, when $\beta =1$, the cdf of the proposed family derived in (6) becomes identical to the baseline cdf. The pdf corresponding to (6) is given by$$\begin{array}{}\text{(7)}& g\left(x;\beta ,\xi \right)=f\left(x;\xi \right)\left(\frac{1+\left(\log \beta \right)\left[1-F\left(x;\xi \right)\right]}{{\beta }^{F\left(x;\xi \right)}}\right),\hspace{1em}x\in ℝ.\end{array}$$

The survival function (sf) and hazard rate function (hrf) corresponding to (6) are, respectively, given by$$\begin{array}{c}\text{(8)}& S\left(x;\beta ,\xi \right)=\frac{1-F\left(x;\xi \right)}{{\beta }^{F\left(x;\xi \right)}},\hspace{1em}x\in ℝ,h\left(x;\beta ,\xi \right)=f\left(x;\xi \right)\left\{\frac{1+\left(\log \beta \right)\left[1-F\left(x;\xi \right)\right]}{1-F\left(x;\xi \right)}\right\},\hspace{1em}x\in ℝ.\end{array}$$

Due to induction of the extra parameter, the Z-family provides greater distributional flexibility. The key motivations for using the Z-family in the practice are as follows:

(i)

A very useful and simple method of introducing an additional parameter to generalize the existing distributions

(ii)

To improve the characteristics and flexibility of the existing models

(iii)

To introduce new distributions having closed form of cdf and sf, as well as hrf

(iv)

To extend the existing distributions by introducing only one parameter, rather than adding two or more parameters

(v)

To provide the best fit to heavy-tailed insurance data sets

(vi)

To provide better fits than other modified models having same or higher number of parameters

3. The Z-Weibull Distribution

Most of the extended forms of distributions are introduced for one of the following aims: (i) an extension of the existing model to improve its characteristics, (ii) to obtain new distribution having a heavy right tail, and (iii) to introduce a model whose empirical fit is good to data. Here, we discuss the Z-Weibull distribution that can possess at least one of these aims. The Weibull random variable has the cdf and pdf given by $F\left(x;\xi \right)=1-e^{-\gamma x^{\alpha }}$ and $f\left(x;\xi \right)=\theta \gamma x^{\alpha -1}{e}^{-\gamma x^{\alpha }}$, respectively, where $\xi =\left(\alpha ,\gamma \right)$. Then, the cdf of the Z-Weibull distribution has the following form:$$\begin{array}{}\text{(9)}& G\left(x;\beta ,\xi \right)=1-\frac{e^{-\gamma x^{\alpha }}}{{\beta }^{\left(1-e^{-\gamma x^{\alpha }}\right)}},\hspace{1em}x\ge 0,\alpha ,\gamma ,\beta >0.\end{array}$$

The corresponding density is given by$$\begin{array}{}\text{(10)}& g\left(x;\beta ,\xi \right)=\frac{\alpha \gamma x^{\alpha -1}{e}^{-\gamma x^{\alpha }}}{{\beta }^{\left(1-e^{-\gamma x^{\alpha }}\right)}}\left(1+\left(\log \beta \right)e^{-\gamma x^{\alpha }}\right),\hspace{1em}x>0.\end{array}$$

The sf, hrf, and reversed hazard rate function (rhrf) of the proposed model are given by$$\begin{array}{c}\text{(11)}& S\left(x;\beta ,\xi \right)=\frac{e^{-\gamma x^{\alpha }}}{{\beta }^{\left(1-e^{-\gamma x^{\alpha }}\right)}},\hspace{1em}x>0,h\left(x;\beta ,\xi \right)=\alpha \gamma x^{\alpha -1}\left(1+\left(\log \beta \right)e^{-\gamma x^{\alpha }}\right),\hspace{1em}x>0,\\ r\left(x;\beta ,\xi \right)=\frac{\alpha \gamma x^{\alpha -1}{e}^{-\gamma x^{\alpha }}}{{\beta }^{\left(1-e^{-\gamma x^{\alpha }}\right)}-e^{-\gamma x^{\alpha }}}\left(1+\left(\log \beta \right)e^{-\gamma x^{\alpha }}\right),\hspace{1em}x>0,\end{array}$$respectively.

Different plots for the pdf of the Z-Weibull distribution for selected parameter values are given in Figure 1.

[figures omitted; refer to PDF]

4. Estimation and Monte Carlo Simulation Study

Several approaches to estimate the model parameter have been introduced in the literature, but the maximum likelihood estimation method is the most commonly employed. The maximum likelihood estimators (MLEs) enjoy several desirable properties and can be used for constructing confidence intervals and regions and also in test statistics. The normal approximation for MLEs in large samples can be easily handled either analytically or numerically. So, we estimate the parameters of the Z-family of distributions from complete samples via the maximum likelihood estimation method. Furthermore, we perform a comprehensive Monte Carlo simulation study to evaluate the performance of the MLEs.

4.1. Maximum Likelihood Estimation

In this section, we obtain the MLEs of the model parameters of the Z-family of distributions from complete samples only. Let $x_1,x_2,\dots ,x_n$ be the observed values from the Z-family of distributions with parameters $\beta$ and $\xi$. The total log-likelihood function corresponding to (7) is given by$$\begin{array}{}\text{(12)}& \ell \left(\beta ,\xi \right)={\displaystyle \sum _{i=1}^n\log \left[f\left(x_i;\xi \right)\right]}-{\displaystyle \sum _{i=1}^{n}F\left(x_i;\xi \right)\left(\log \beta \right)}+{\displaystyle \sum _{i=1}^n\log \left[1+\left(\log \beta \right)\left(1-F\left(x_i;\xi \right)\right)\right]}.\end{array}$$

The partial derivatives of (12) are$$\begin{array}{c}\text{(13)}& \frac{\partial }{\partial \beta }\ell \left(\beta ,\xi \right)=-{\displaystyle \sum _{i=1}^n\frac{F\left(x_i;\xi \right)}{\beta }}+\frac{1}{\beta }{\displaystyle \sum _{i=1}^n\frac{\left(1-F\left(x_i;\xi \right)\right)}{1+\left(\log \beta \right)\left(1-F\left(x_i;\xi \right)\right)}},\frac{\partial }{\partial \xi }\ell \left(\beta ,\xi \right)={\displaystyle \sum _{i=1}^n\frac{\partial f\left(\left(x_i;\xi \right)/\partial \xi \right)}{f\left(x_i;\xi \right)}}-{\displaystyle \sum _{i=1}^n\frac{\left(\log \beta \right)\partial F\left(\left(x_i;\xi \right)/\partial \xi \right)}{F\left(x_i;\xi \right)}}-{\displaystyle \sum _{i=1}^n\frac{\partial F\left(\left(x_i;\xi \right)/\partial \xi \right)}{1+\left(\log \beta \right)\left(1-F\left(x_i;\xi \right)\right)}}.\end{array}$$

Setting $\left(\partial /\partial \beta \right)\ell \left(\beta ,\xi \right)$ and $\left(\partial /\partial \xi \right)\ell \left(\beta ,\xi \right)$ equal to zero and solving numerically these expressions simultaneously yield the MLEs of $\left(\beta ,\xi \right)$.

4.2. Monte Carlo Simulation Study

This section offers a comprehensive simulation study to assess the behavior of the MLEs. The Z-family is easily simulated by inverting (6) as follows: if U has a uniform U (0,1) distribution, then the nonlinear equation is$$\begin{array}{}\text{(14)}& \log \left(1-u\right)+\gamma x^{\alpha }+\left(1-e^{-\gamma x^{\alpha }}\right)\left(\log \beta \right)=0.\end{array}$$

Expression (14) can be used to simulate any special subcase of the Z-family. Here, we consider the Z-Weibull distribution to assess the behavior of the MLEs of the proposed method. We simulate the Z-Weibull distribution for two sets of parameters (set 1: $\alpha =0.9$, $\beta =0.3$, and $\gamma =0.5$ and set 2: $\alpha =1.3$, $\beta =0.8$, and $\gamma =1$). The simulation is performed via statistical software R through the library (rootSolve) command mle. The number of Monte Carlo replications made was 750 times. For maximizing the log-likelihood function, we use the method = “L-BFGS-B” algorithm with optim(). The evaluation of the estimators was performed via the following quantities for each sample size: the empirical mean squared errors (MSEs) are calculated using the R package from the Monte Carlo replications. The MLEs are determined for each piece of simulated data, say, $\left({\widehat{\alpha }}_i,{\widehat{\gamma }}_i,{\widehat{\beta }}_i\right)$ for $i=\mathrm{1,2},\dots ,750$, and the biases and MSEs are computed, respectively, by$$\begin{array}{c}\text{(15)}& \text{bias}\left(w\right)=\left(\frac{1}{750}\right){\displaystyle \sum _{i=1}^{750}\left({\widehat{w}}_i-w\right)},\text{MSE}\left(w\right)=\left(\frac{1}{750}\right){\displaystyle \sum _{i=1}^{750}{\left({\widehat{w}}_i-w\right)}^2},\end{array}$$for $w=\alpha ,\gamma ,\beta$. We consider the sample sizes of $n=\mathrm{25,100,200,300,400,500,600,700,750}$. The empirical results are given in Tables 1 and 2. Corresponding to Tables 1 and 2, the simulation results are graphically displayed in Figures 2 and 3. Based on Tables 1 and 2 and Figures 2 and 3, the following results are concluded:

(i)

Biases for all parameters are positive

(ii)

The parameters tend to be stable

(iii)

Estimated biases decrease when the sample size n increases

(iv)

Estimated MSEs decay toward zero when the sample size n increases

Table 1

Simulation results for the Z-Weibull distribution.

Set 1: $\alpha =0.9$, $\beta =0.3$, and $\gamma =1.2$
$n$	Parameters	MLEs	MSEs	Bias
25	$\alpha$	1.1150	0.1937	0.2150
	$\beta$	1.2490	3.6872	0.9490
	$\gamma$	2.0431	4.7184	0.8431
50	$\alpha$	1.0456	0.1526	0.1456
	$\beta$	0.9536	2.4614	0.6536
	$\gamma$	1.9282	4.3933	0.7281
100	$\alpha$	0.9970	0.0783	0.0970
	$\beta$	0.7277	1.3741	0.4277
	$\gamma$	1.6722	2.3642	0.4722
200	$\alpha$	0.9605	0.0472	0.0605
	$\beta$	0.5478	0.6878	0.2478
	$\gamma$	1.4691	1.1474	0.2691
400	$\alpha$	0.9281	0.0163	0.0281
	$\beta$	0.3769	0.0901	0.0769
	$\gamma$	1.3115	0.4336	0.1115
600	$\alpha$	0.9193	0.0085	0.0193
	$\beta$	0.3471	0.0383	0.0471
	$\gamma$	1.2424	0.1812	0.0424
800	$\alpha$	0.9128	0.0067	0.0128
	$\beta$	0.3313	0.0279	0.0313
	$\gamma$	1.2405	0.1318	0.0405
1000	$\alpha$	0.9006	0.0042	0.0076
	$\beta$	0.3184	0.0119	0.0184
	$\gamma$	1.2160	0.0928	0.0360

Table 2

Simulation results for the Z-Weibull distribution.

Set 2: $\alpha =1.3$, $\beta =0.8$, and $\gamma =1.4$
$n$	Parameters	MLEs	MSEs	Bias
25	$\alpha$	1.3094	0.2002	0.2094
	$\beta$	1.9925	5.7399	1.2925
	$\gamma$	2.1937	4.3061	0.7937
50	$\alpha$	1.2340	0.1333	0.1340
	$\beta$	1.6368	4.0681	0.9368
	$\gamma$	2.0939	3.3865	0.6939
100	$\alpha$	1.1766	0.0697	0.0766
	$\beta$	1.2543	2.1560	0.5543
	$\gamma$	1.8715	2.1136	0.4715
200	$\alpha$	1.1396	0.0368	0.0396
	$\beta$	0.9740	0.9313	0.2740
	$\gamma$	1.6771	1.0320	0.2771
400	$\alpha$	1.1194	0.0153	0.0194
	$\beta$	0.8080	0.2293	0.1080
	$\gamma$	1.5204	0.3718	0.1204
600	$\alpha$	1.1171	0.0094	0.0171
	$\beta$	0.7860	0.1229	0.0860
	$\gamma$	1.4516	0.1869	0.0516
800	$\alpha$	1.1093	0.0065	0.0093
	$\beta$	0.7495	0.0833	0.0495
	$\gamma$	1.4443	0.1219	0.0443
1000	$\alpha$	1.1050	0.0049	0.0050
	$\beta$	0.7807	0.0537	0.0307
	$\gamma$	1.4147	0.0935	0.0347

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

5. Actuarial Measures

One of the most important tasks of actuarial science institutions is to evaluate the exposure to market risk in a portfolio of instruments, which arise from changes in underlying variables such as prices of equity, interest rates, or exchange rates. In this section, we calculate some important risk measures such as value at risk (VaR) and tail value at risk (TVaR) for the proposed distribution, which play a crucial role in portfolio optimization under uncertainty.

5.1. Value at Risk

In the context of actuarial sciences, the VaR is widely used by practitioners as a standard financial market risk measure. It is also known as the quantile risk measure or quantile premium principle. The VaR is always specified with a given degree of confidence, say q (typically 90%, 95%, or 99%), and represents the percentage loss in the portfolio value that will be equaled or exceeded only X percent of the time. VaR of a random variable X is the q^th quantile of its cdf, see Artzner [29]. If X follows the proposed method, then the VaR of X is$$\begin{array}{}\text{(16)}& x_q=G^{-1}\left(u\right)=F^{-1}\left(t\right),\end{array}$$where t is the solution of ${\beta }^t\left(1-q\right)+t-1=0$.

5.2. Tail Value at Risk

Another important measure is TVaR, also known as conditional tail expectation (CTE) or tail conditional expectation (TCE), which is used to quantify the expected value of the loss given that an event outside a given probability level has occurred. Let X follow the Z-Weibull distribution; then, TVaR of X is derived as$$\begin{array}{}\text{(17)}& {\text{TVaR}}_q\left(x\right)=\left(\frac{1}{1-q}\right){\displaystyle {\int }_{{\text{VaR}}_q}^{\infty }x\text{g}\left(x\right)\text{d}x}.\end{array}$$

Using (10) in (17), we get$$\begin{array}{c}\text{(18)}& {\text{TVaR}}_q\left(x\right)=\left(\frac{\alpha \gamma }{1-q}\right){\displaystyle {\int }_{{\text{VaR}}_q}^{\infty }\frac{x^{\alpha +1-1}{e}^{-\gamma x^{\alpha }}}{{\beta }^{\left(1-e^{-\gamma x^{\alpha }}\right)}}\left(1+\left(\log \beta \right)e^{-\gamma x^{\alpha }}\right)\text{d}x},{\text{TVaR}}_q\left(x\right)=\left(\frac{\alpha \gamma }{1-q}\right){\displaystyle \sum _{i=0}^{\infty }{\displaystyle \sum _{j=0}^i\frac{{\left(\log \beta \right)}^i{\left(-1\right)}^{\left(i+j\right)}}{\left(i-j\right)!j!}}{\displaystyle {\int }_{{\text{VaR}}_q}^{\infty }{x}^{\alpha +1-1}{e}^{-\gamma \left(j+1\right)x^{\alpha }}\text{d}x}}+\left(\frac{\alpha \gamma }{1-q}\right){\displaystyle \sum _{i=0}^{\infty }{\displaystyle \sum _{j=0}^i\frac{{\left(\log \beta \right)}^i{\left(-1\right)}^{i+j}}{\left(i-j\right)!j!}\left(\log \beta \right){\displaystyle {\int }_{{\text{VaR}}_q}^{\infty }{x}^{\alpha +1-1}{e}^{-\gamma \left(j+2\right)x^{\alpha }}\text{d}x}}}.\end{array}$$

Recall the definition of incomplete gamma function in the form $\mathrm{\Gamma }\left(\alpha ,x\right)={\displaystyle {\int }_x^{\infty }{t}^{\alpha -1}{e}^{-t}\text{d}t}$, so from (18), we get$$\begin{array}{}\text{(19)}& {\text{TVaR}}_q\left(x\right)=\left(\frac{1}{1-q}\right){\displaystyle \sum _{i=0}^{\infty }{\displaystyle \sum _{j=0}^i\frac{{\left(\log \beta \right)}^i{\left(-1\right)}^{\left(i+j\right)}}{\left(i-j\right)!j!{\gamma }^{\left(1/\alpha \right)}}\frac{\Gamma \left(\left(1/\alpha \right)+1,\gamma \left(j+1\right){\left({\text{VaR}}_q\right)}^{\alpha }\right)}{{\left(j+1\right)}^{\left(\left(1/\alpha \right)+1\right)}}}}+\frac{1}{1-q}{\displaystyle \sum _{i=0}^{\infty }{\displaystyle \sum _{j=0}^i\frac{{\left(\log \beta \right)}^i{\left(-1\right)}^{\left(i+j\right)}}{\left(i-j\right)!j!{\gamma }^{\left(1/\alpha \right)}}\frac{\Gamma \left(\left(1/\alpha \right)+1,\gamma \left(j+2\right){\left({\text{VaR}}_q\right)}^{\alpha }\right)}{{\left(j+2\right)}^{\left(\left(1/\alpha \right)+1\right)}}}}.\end{array}$$

5.3. Numerical Study of the Risk Measures

n this section, we provide numerical study of the VaR and TVaR for the Weibull and Z-Weibull distributions for different sets of parameters. The process is described as follows:

(1)

Random samples of size n = 100 are generated from the Weibull and Z-Weibull models

(2)

The parameters have been estimated via the maximum likelihood method

(3)

1000 repetitions are made to calculate the VaR and TVaR

The numerical results of the risk measures are provided in Tables 3 and 4 and displayed graphically in Figures 4 and 5 corresponding to each table.

Table 3

Simulation results of the VaR and TVaR for the selected values of the parameters.

Dist.	Par	Level of significance	VaR	TVaR
Weibull	$\alpha =1.2$$\gamma =0.85$	0.700	4.5373	7.3911
		0.750	5.0945	7.9076
		0.800	5.7601	8.5304
		0.850	6.5967	9.3205
		0.900	7.7437	10.4136
		0.950	9.6399	12.2384
		0.990	11.4765	14.0197
		0.999	19.5352	21.9150
Z-Weibull	$\alpha =1.2$$\beta =1.4$$\gamma =0.85$	0.700	4.3778	13.2772
		0.750	5.5376	14.9465
		0.800	7.1015	17.1137
		0.850	9.3408	20.1026
		0.900	12.8999	24.6727
		0.950	20.0017	33.4043
		0.990	28.2956	43.2302
		0.999	80.4988	101.7184

Table 4

Simulation results of the VaR and TVaR for the selected values of the parameters.

Dist.	Par	Level of significance	VaR	TVaR
Weibull	$\alpha =0.65$$\gamma =0.7$	0.700	1.5273	4.4325
		0.750	1.9174	4.9762
		0.800	2.4392	5.6792
		0.850	3.1801	6.6436
		0.900	4.3463	8.1092
		0.950	6.6444	10.8853
		0.990	9.2950	13.9806
		0.999	25.5667	32.0408
Z-Weibull	$\alpha =0.65$$\beta =0.8$$\gamma =0.7$	0.700	1.7585	6.9318
		0.750	2.3222	7.9131
		0.800	3.1171	9.2175
		0.850	4.3112	11.0668
		0.900	6.3169	13.9926
		0.950	10.6148	19.8560
		0.990	16.0027	26.8019
		0.999	55.1566	73.1187

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

The simulation is performed for Weibull and Z-Weibull for the selected values of parameters. A model with higher values of the risk measures is said to have a heavier tail. The simulated results provided in Tables 3 and 4 show that the proposed Z-Weibull model has higher values of the risk measures than the traditional Weibull distribution. The simulation results are graphically displayed in Figures 4 and 5, which show that the proposed model has a heavier tail than the Weibull distribution.

6. Practical Illustration via the Earthquake Insurance Data

The main applications of the heavy tail models are the so-called extreme value theory or insurance loss phenomena. In this section, we consider heavy-tailed earthquake insurance data to illustrate the usefulness of the Z-Weibull model. The data are reported by the “National Centers for Environmental Information” available at https://ngdc.noaa.gov/hazard/earthqk.shtml. We compare the goodness-of-fit results of the Z-Weibull distribution with the other well-known heavy-tailed distributions. The distribution functions of the competitive models are as follows:

(i)

W-claim:$$\begin{array}{}\text{(20)}& G\left(x\right)=\frac{\sigma {\left(1-e^{-\gamma x^{\alpha }}\right)}^2}{1-\left(1-\sigma \right)\left(1-e^{-\gamma x^{\alpha }}\right)},\hspace{1em}x\ge ,\alpha ,\gamma ,\sigma >0.\end{array}$$

(ii)

Weibull:

$$\begin{array}{}\text{(21)}& G\left(x\right)=1-e^{-\gamma x^{\alpha }},\hspace{1em}x\ge 0,\alpha ,\gamma >0.\end{array}$$

(iii)

B-XII:$$\begin{array}{}\text{(22)}& G\left(x\right)=1-{\left(1+x^c\right)}^{-k},\hspace{1em}x\ge 0,c,k>0.\end{array}$$

(iv)

GE:$$\begin{array}{}\text{(23)}& G\left(x\right)={\left(1-e^{-\gamma x^{\alpha }}\right)}^a,\hspace{1em}x\ge 0,a,\alpha ,\gamma >0.\end{array}$$

(v)

EL:$$\begin{array}{}\text{(24)}& G\left(x\right)={\left(1-{\left(1+\gamma x\right)}^{-\alpha }\right)}^a,\hspace{1em}x\ge 0,a,\alpha ,\gamma >0.\end{array}$$

(vi)

BW:$$\begin{array}{}\text{(25)}& G\left(x\right)=I_{\left(1-e^{-\gamma x^{\alpha }}\right)}\left(a,b\right),\hspace{1em}x\ge 0,a,b,\alpha ,\gamma >0.\end{array}$$

 

Next, we consider certain analytical measures in order to verify which distribution fits better the considered data. These analytical measures include (i) discrimination measures, such as Akaike information criterion (AIC), Bayesian information criterion (BIC), Hannan–Quinn information criterion (HQIC), and consistent Akaike information criterion (CAIC), and (ii) two other goodness-of-fit measures including Anderson–Darling (AD) test statistic and $-2\ell$. The discrimination measures are given as follows.

(vii)

Akaike information criterion is given by$$\begin{array}{}\text{(26)}& \text{AIC}=2k-2\ell .\end{array}$$

(viii)

Consistent Akaike information criterion is given by$$\begin{array}{}\text{(27)}& \text{CAIC}=\frac{2\text{nk}}{n-k-1}-2\ell .\end{array}$$

(ix)

Bayesian information criterion is given by$$\begin{array}{}\text{(28)}& \text{BIC}=k\log \left(n\right)-2\ell .\end{array}$$

(x)

Hannan–Quinn information criterion is given by

$$\begin{array}{}\text{(29)}& \text{HQIC}=2k\log \left(\log \left(n\right)\right)-2\ell ,\end{array}$$

 

where $\ell$ denotes the log-likelihood function, k is the number of model parameters, and n is the sample size.

(xi)

The AD test statistic is given by

$$\begin{array}{}\text{(30)}& \text{AD}=-n-\left(\frac{1}{n}\right){\displaystyle \sum _{i=1}^n\left(2i-1\right)\left[\log G\left(x_i\right)+\log \left\{1-G\left(x_{n-i+1}\right)\right\}\right]},\end{array}$$where n is the sample size and $x_i$ is the $i^{\text{th}}$ observation in the sample, calculated when the data are sorted in the ascending order.

All the computations are carried out using the optim() R-function with the argument method = “BFGS” (see Appendix). A model with lowest values for these measures could be chosen as the best model to fit the data. The values of MLEs of the parameters along with standard errors in parenthesis are presented in Table 5, whereas the discrimination measures are displayed in Table 6. The AD statistic and $-2\ell$ are provided in Table 7. Based on the considered data set, we have observed that the Z-Weibull distribution is the best fitted model among the above considered models.

Table 5

Estimated values of the proposed and other competitive models.

Dist.	$\widehat{\alpha }$	$\widehat{\gamma }$	$\widehat{\beta }$	$\widehat{a}$	$\widehat{b}$	$\widehat{c}$	$\widehat{k}$	$\widehat{\sigma }$
Z-Weibull	2.346 (0.020)	0.564 (0.029)	2.295 (0.248)
W-claim	2.023 (0.076)	0.709 (0.056)						1.783 (0.197)
Weibull	2.110 (0.091)	0.889 (0.086)
B-XII						1.480 (0.030)	2.851 (0.107)
GE		2.283 (0.015)		4.308 (0.053)
EL	4.804 (4.910)	3.879 (5.400)		4.335 (0.054)
BW	1.609 (0.055)	0.644 (0.304)		1.616 (0.087)	2.237 (1.228)

Table 6

Discrimination measures of the proposed and other competitive models.

Dist.	AIC	BIC	CAIC	HQIC
Z-Weibull	23286.95	23310.57	23289.92	23294.69
W-claim	23345.86	23364.09	23355.04	23354.73
Weibull	23379.64	23395.38	23386.84	23384.80
B-XII	24197.63	24213.37	24197.63	24202.78
GE	23846.71	23862.45	23853.57	23851.86
EL	23863.99	23887.60	23871.99	23871.72
BW	23338.24	23369.73	23343.17	23348.56

Table 7

Analytical measures of the fitted models.

Dist.	$-2\ell$	AD
Z-Weibull	11640.48	0.840
W-claim	11667.98	0.986
Weibull	11687.82	1.071
B-XII	12096.81	5.139
GE	11921.35	4.965
EL	11928.99	5.012
BW	11665.12	0.952

As we can see, the results show that the Z-Weibull distribution provides a better fit than the other competitors. Hence, the proposed model can be used as a best candidate model for modeling insurance data sets. Furthermore, in support of Tables 6 and 7, the estimated cdf and pdf are plotted in Figure 6. The Kaplan–Meier survival plot and PP plot are sketched in Figure 7, whereas the QQ plot of the proposed distribution and box plot of the earthquake data are presented in Figure 8. From the estimated pdf in Figure 6, it is clear that the proposed distribution provides an adequate fit to the heavy-tailed earthquake data. From Figures 6 and 7, we can easily detect that the proposed distribution fits the estimated cdf and Kaplan–Meier survival plots very closely. The PP and QQ plots which serve as a tool for graphical display of the analytical measures show that the Z-Weibull distribution provides the best fit to real data. Finally, the box plot (Figure 8) of the data is graphical evidence that the data possess tail skewed to the right.

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

Furthermore, using the earthquake insurance data, we obtained the values of the Kolmogorov–Smirnov (KS) statistic of the proposed and other competing models. Then, we applied the parametric bootstrap technique [30] and bootstrapped the $p$ value for all the competing models. The KS statistic and the corresponding bootstrapped $p$ value are provided in Table 8. Based on the results provided in Table 8, we conclude that the proposed model is the best candidate model among the competing distributions for modeling the insurance claim data.

Table 8

KS and the corresponding bootstrapped $p$ value of the proposed and other models.

Dist.	KS	Bootstrapped $p$ value
Z-Weibull	0.205	0.976
W-claim	0.367	0.806
Weibull	0.476	0.704
B-XII	0.954	0.408
GE	0.590	0.605
EL	0.869	0.502
BW	0.406	0.775

7. Bayesian Estimation

Bayesian inference procedure has been taken into consideration by many statistical researchers, especially those in the field of survival analysis and reliability engineering. In this section, complete sample data are analyzed through the Bayesian point of view. We assume that the parameters $\alpha ,\gamma$, and $\beta$ of the Z-Weibull distribution have independent prior distributions as follows:$$\begin{array}{c}\text{(31)}& \alpha \sim \text{Gamma}\left(a,b\right),\gamma \sim \text{Gamma}\left(c,d\right),\\ \beta \sim \text{Gamma}\left(e,f\right),\end{array}$$where $a,b,c,d,e$, and $f$ are positive. More about choosing gamma priors, refer Kundu and Howlader [31], S. Dey and T. Dey [32], Dey et al. [33], and Dey et al. [34]. Hence, the joint prior density function is formulated as follows:$$\begin{array}{}\text{(32)}& \pi \left(\alpha ,\gamma ,\beta \right)=\frac{b^{\alpha }{d}^c{f}^e}{\left(\Gamma a\right)\left(\Gamma c\right)\left(\Gamma e\right)}{\alpha }^{a-1}{\gamma }^{c-1}{\gamma }^{e-1}\exp \left\{-\left(b\alpha +d\gamma +f\beta \right)\right\}.\end{array}$$

In the Bayesian estimation, the actual value of the parameter may be adversely affected by the loss when choosing an estimator. This loss can be measured by a function of the parameter and the corresponding estimator. Five well-known loss functions and associated Bayesian estimators and corresponding posterior risk are presented in Table 9.

Table 9

Bayes estimator and posterior risk under different loss functions.

Loss function	Bayes estimator	Posterior risk
L1 = SELF = ${\left(\theta -d\right)}^2$	$E\left(\theta /x\right)$	$V\left(\theta /x\right)$
L2 = WSELF = $\left(\left({\left(\theta -d\right)}^2\right)/\theta \right)$	${\left(E\left({\theta }^{-1}/x\right)\right)}^{-1}$	$E\left(\theta /x\right)-{\left(E\left({\theta }^{-1}/x\right)\right)}^{-1}$
L3 = MSELF = ${\left(1-\left(d/\theta \right)\right)}^2$	$\left(E\left({\theta }^{-1}/x\right)\right)/\left(E\left({\theta }^{-2}/x\right)\right)$	$1-\left(\left(E{\left({\theta }^{-1}/x\right)}^2\right)/\left(E{\left({\theta }^{-2}/x\right)}^2\right)\right)$
L4 = PLF = $\left(\left({\left(\theta -d\right)}^2\right)/d\right)$	$\sqrt{E\left({\theta }^2/x\right)}$	$2\left(\sqrt{E\left({\theta }^2/x\right)}-E\left(\theta /x\right)\right)$
L5 = KLF = $\left(\sqrt{\left(d/\theta \right)-\sqrt{\left(d/\theta \right)}}\right)$	$\sqrt{\left(E\left(\theta /x\right)\right)/\left(E\left({\theta }^{-1}/x\right)\right)}$	$2\left(\sqrt{E\left(\theta /x\right)E\left({\theta }^{-1}/x\right)}-1\right)$

For more details, see Calabria and Pulcini [35] and Dey et al. [36]. Next, we provide the posterior probability distribution for a complete data set. We define the function $\phi$ as$$\begin{array}{}\text{(33)}& \phi \left(\alpha ,\gamma ,\beta \right)={\alpha }^{a-1}{\gamma }^{c-1}{\beta }^{e-1}\exp \left\{-\left(b\alpha +d\gamma +f\beta \right)\right\},\hspace{1em}\alpha >0,\gamma >0,\beta >0,\beta \ne 1.\end{array}$$

The joint posterior distribution in terms of a given likelihood function L(data) and joint prior distribution $\pi \left(\alpha ,\gamma ,\beta \right)$ is defined as$$\begin{array}{}\text{(34)}& {\pi }^{\ast }\left(\left.\alpha ,\gamma ,\beta \right|\text{data}\right)\alpha \pi \left(\alpha ,\gamma ,\beta \right)L\left(\text{data}\right).\end{array}$$

Hence, the joint posterior density of parameters $\alpha ,\gamma$, and $\beta$ for complete sample data is obtained by combining the likelihood function and joint prior density (32). Therefore, the joint posterior density function is given by$$\begin{array}{}\text{(35)}& {\pi }^{\ast }\left(\left.\alpha ,\gamma ,\beta \right|\underset{¯}{x}\right)=K\phi \left(\left.\alpha ,\gamma ,\beta \right|\underset{¯}{x}\right){\displaystyle \prod _{i=1}^n\frac{\alpha \gamma x_i^{\alpha -1}{e}^{-\gamma x_i^{\alpha }}\left(1+\log \left(\beta \right)\right)e^{-\gamma x_i^{\alpha }}}{{\beta }^{1-e^{-\gamma x_i^{\alpha }}}}},\end{array}$$where$$\begin{array}{}\text{(36)}& K^{-1}={\displaystyle {\int }_0^{\infty }{\displaystyle {\int }_0^{\infty }{\displaystyle {\int }_0^{\infty }\phi \left(\alpha ,\gamma ,\beta \right){\displaystyle \prod _{i=1}^n\frac{\alpha \gamma x_i^{\alpha -1}{e}^{-\gamma x_i^{\alpha }}\left(1+\log \left(\beta \right)\right)e^{-\gamma x_i^{\alpha }}}{{\beta }^{1-e^{-\gamma x_i^{\alpha }}}}\text{d}\alpha \text{d}\gamma \text{d}\beta }}}}.\end{array}$$

Moreover, the marginal posterior density of $\alpha$, $\gamma$, and $\beta$, assuming that $\mathrm{\Theta }=\left(\alpha ,\gamma ,\beta \right)$, is given by$$\begin{array}{}\text{(37)}& \pi \left(\left.{\Theta }_i\right|\underset{¯}{x}\right)={\displaystyle {\int }_0^{\infty }{\displaystyle {\int }_0^{\infty }{\pi }^{\ast }\left(\left.{\Theta }_i\right|\underset{¯}{x}\right)\text{d}{\Theta }_j\text{d}{\Theta }_k}},\end{array}$$where $i,j,k=\mathrm{1,2,3},i\ne j\ne k$, and also ${\mathrm{\Theta }}_i$ is the $i$-th member of vector $\mathrm{\Theta }$.

From (35) and (37), it is clear that there is no closed form for the Bayesian estimators under the five loss functions described in Table 8. Therefore, we use the MCMC procedure based on 10,000 replicates to compute Bayesian estimators.

Because of intractable integrals associated with joint posterior and marginal posterior distributions, one needs to use numerical software to solve the integral equations numerically. The two most popular MCMC methods are the Metropolis–Hastings algorithm [37, 38] and the Gibbs sampling [39]. Gibbs sampling is a special case of the Metropolis–Hastings algorithm which generates a Markov chain by sampling from the full set of conditional distributions. Often, Bayesian inference requires computing intractable integrals to generate posterior samples. In practice, simulations related to Gibbs sampling are conducted through special software WinBUGS. WinBUGS software was developed in 1997 to simulate data of complex posterior distributions, where analytical or numerical integration techniques cannot be applied. One may also use OpenBUGS software, which is an open-source version of WinBUGS. Using Gibbs sampling, we obtain samples from the joint posterior distribution and then use OpenBUGS software to carry out the Bayesian analysis.

The process is described as follows:

(i)

Gibbs sampling technique is used to generate posterior samples

(ii)

10,000 replicates are made to compute the Bayesian estimators via OpenBUGS software

(iii)

The idea of Congdon [40] was to implement and choose $a=b=c=d=e=f=0.0001$ as we do not have any prior information about hyperparameters

The corresponding Bayesian estimates and posterior risk are provided in Table 10. Table 11 provides 95% credible and HPD intervals for each parameter of the Z-Weibull distribution. Moreover, we provide the posterior summary plots in Figures 9–11. These plots confirm that the sampling process is of the prime quality, and the convergence does occur.

Table 10

Bayesian estimates and their posterior risks of the parameters under different loss functions based on the earthquake insurance data.

Bayes	$\widehat{\alpha }$		$\widehat{\gamma }$		$\widehat{\beta }$
Loss functions	Estimate	Risk	Estimate	Risk	Estimate	Risk
SELF	2.34808	0.00048	0.56078	0.00113	2.35742	0.08661
WSELF	2.34787	0.00020	0.55876	0.00201	2.32162	0.03580
MSELF	2.34766	8.90e − 05	0.55675	0.00359	2.28686	0.01497
PLF	2.34819	0.00020	0.56179	0.00202	2.37572	0.03659
KLF	2.34798	8.879e − 05	0.55977	0.00360	2.33945	0.01536

Table 11

Credible and HPD intervals of the parameters for the earthquake insurance data.

Parameters	Credible interval	HPD interval
$\alpha$	(2.334, 2.364)	(2.303, 2.388)
$\gamma$	(0.536, 0.585)	(0.410, 0.625)
$\beta$	(2.139, 2.562)	(1.827, 2.883)

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

8. Discussion and Future Framework

Statistical decision theory addresses the state of uncertainty and provides a rational framework for dealing with problems of actuarial and financial decision-making. The insurance data sets are generally skewed to the right and heavy-tailed. The traditional distributions are not flexible enough to counter complex forms of data such as insurance science data.

Due to the importance of statistical distributions in actuarial sciences, a number of papers have been appeared in the literature aiming to improve the characteristics of the existing distributions. Although this has been achieved, unfortunately, the numbers of parameters have been increased, and the estimation of parameters and derivation of mathematical problems become complicated.

To provide a better description of the insurance science data, therefore, in this study, an attempt has been made to introduce a new family of statistical distributions aiming to increase the flexibility of the existing distributions. A special submodel of the proposed family offers the best fitting to the heavy-tailed insurance science data. The maximum likelihood method is adopted to estimate the model parameters, and a comprehensive Monte Carlo simulation study is done to evaluate the behavior of the estimators.

To show the usefulness of the proposed method in insurance sciences, a real-life application of the earthquake insurance data is discussed. Analyzing the data set, it showed that the proposed model performs much better than the other competitive distributions.

From the above discussion, it is obvious that the researchers are always in search of new flexible distributions. Therefore, to bring further flexibility in the proposed model, we suggest to introduce its extended versions. The proposed method can further be extended by introducing a shape parameter to the model.

(i)

A random variable X is said to follow the extended version of the Z-family if its cdf is given by$$\begin{array}{}\text{(38)}& G\left(x;a,\beta ,\xi \right)=1-\left\{\frac{1-F{\left(x;\xi \right)}^a}{{\beta }^{F\left(x;\xi \right)}}\right\},\hspace{1em}a,\beta >0,x,\xi \in ℝ,\end{array}$$

 

where $a$ is the additional shape parameter. For $a=1$, expression (38) reduces to (6). The new proposal may be named as the exponentiated Z-family. For the illustrative purposes, one may consider its special subcase, called the exponentiated Z-Weibull (EZ-Weibull) distribution, defined by the cdf$$\begin{array}{}\text{(39)}& G\left(x;a,\beta ,\xi \right)=1-\left\{\frac{1-{\left(1-e^{-\gamma x^{\alpha }}\right)}^a}{{\beta }^{\left(1-e^{-\gamma x^{\alpha }}\right)}}\right\},\hspace{1em}a,\alpha ,\gamma ,\beta >0,x\in ℝ.\end{array}$$

 

Due to the introduction of the additional shape parameter, the suggested extension may be much flexible in modeling data in insurance sciences and other related fields.

(ii)

Another extension of the Z-family is given by$$\begin{array}{}\text{(40)}& G\left(x;\theta ,\beta ,\xi \right)=1-{\left\{\frac{1-F\left(x;\xi \right)}{{\beta }^{F\left(x;\xi \right)}}\right\}}^{\theta },\hspace{1em}\theta ,\beta >0,x,\xi \in ℝ,\end{array}$$

 

where $\theta$ is the additional shape parameter. For $\theta =1$, expression (40) reduces to (6). The model defined in (40) may be named as the extended Z-family.

(iii)

Another generalized version of the Z-family can be introduced via

$$\begin{array}{}\text{(41)}& G\left(x;a,\theta ,\beta ,\xi \right)=1-{\left\{\frac{1-F{\left(x;\xi \right)}^a}{{\beta }^{F\left(x;\xi \right)}}\right\}}^{\theta },\hspace{1em}a,\theta ,\beta >0,x,\xi \in ℝ,\end{array}$$

 

where $a$ and $\theta$ are the additional shape parameters. Clearly, for $a=1$, expression (41) reduces to (40), and for $\theta =1$, expression (41) reduces to (38). However, for $a=\theta =1$, expression (41) reduces to (6). The model introduced in (41) may be named as the extended exponentiated Z-family.

(iv)

Another generalized version of the new extended Z-family can be introduced via$$\begin{array}{}\text{(42)}& G\left(x;a,\theta ,\beta ,\xi \right)=1-{\left\{\frac{1-F{\left(x;\xi \right)}^a}{{\beta }^{F\left(x;\xi \right)}}\right\}}^{\theta },\hspace{1em}a,\theta ,\beta >0,x,\xi \in ℝ.\end{array}$$

9. Concluding Remarks

A variety of methods for proposing new heavy-tailed distributions have been developed to model data related to financial and actuarial sciences. We carried out this area of research further and introduced a new heavy-tailed distribution family. Some distributional properties are derived, and the method of maximum likelihood estimation is discussed to estimate the model parameters. In addition to distributional properties, some actuarial properties are also derived. Based on the actuarial measures, a comprehensive simulation study is conducted. We focused our concentration on a three-parameter special model called the Z-Weibull distribution. To prove the potential and usefulness of the Z-Weibull distribution, earthquake insurance data are analyzed, and its comparison is made with the other well-known distributions. While analyzing the earthquake insurance data, it is observed that the proposed model performs better than the other competitive models. Bayesian analysis using the earthquake data is also provided. Finally, some new extensions are also suggested which may further improve the characteristics of the proposed family.

Disclosure

This article is drafted from the PhD work of the first author (Zubair Ahmad).

Appendix

Appendix

R-Code for Analysis

Note: in the following R-code, a is used for $\alpha$, $g$ is used for $\gamma$, b is used for $\beta$, and pm is used for the proposed model.

 

data ← read.csv(file.choose(), header = TRUE)

 

data = data[1]

 

data = data[!is.na(data)]

 

data = data/100

 

data

 

hist(data)

 

###—Probability density function

 

pdf_pm ← function(par,x)

 

{

 

a = par[1]

 

$g$ = par[2]

 

b = par[3]

 

(a ∗ $g$ ∗ (x^(a − 1)) ∗ (exp($-g$ ∗ x^a)) ∗ (b^(exp($-g$ ∗ x^a) − 1))) ∗ (1 + log(b) ∗ exp($-g$ ∗ x^a))

 

###—Cumulative distribution function.

 

cdf_pm ← function(par,x)

 

{

 

a = par[1]

 

$g$ = par[2]

 

b = par[3]

 

1 − (exp($-g$ ∗ x^a)/(b^(1 − exp($-g$ ∗ x^a))))

 

}

 

set.seed(0)

 

goodness.fit(pdf = pdf_pm,

 

cdf = cdf_pm,

 

starts = c(0.5,0.5,0.5), data = data,

 

method = “BFGS,” domain = c(0,Inf),mle = NULL)

 

%%%%---------------------------------------------------

 

% Estimated pdf

 

%%%%---------------------------------------------------

 

x ← read.csv(file.choose(), header = TRUE)

 

x = x[,1]

 

x = x[!is.na(x)]

 

x = x/100

 

x

 

###—Parameter values

 

a = 2.346

 

$g$ = 0.564

 

b = 2.295

 

pdf = (a ∗ $g$ ∗ (x^(a − 1)) ∗ (exp($-g$ ∗ x^a)) ∗ (b^(exp($-g$ ∗ x^a) − 1))) ∗ (1 + log(b) ∗ exp($-g$ ∗ x^a))

 

f = pdf

 

x = sort(x)

 

yrange = c(0,1)

 

xrange = c(min(x),max(x))

 

hist(x, freq = FALSE, breaks = 15, xlim = xrange, ylim = yrange,

 

ylab = “Estimated pdf,”xlab = “x,” main = ””)

 

par(new = TRUE)

 

lines(x,f, xlim = xrange, lty = 1, ylab = ,” ylim = yrange, lwd = 3,

 

col = “blue,” xlab = ””)

 

par(new = TRUE)

 

%%%%---------------------------------------------------

 

%Estimated cdf

 

%%%%---------------------------------------------------

 

x ← read.csv(file.choose(), header = TRUE)

 

x = x[,1]

 

x = x[!is.na(x)]

 

x = x/100

 

x

 

###—Parameter values

 

a = 2.346

 

$g$ = 0.564

 

b = 2.295

 

x ← sort(x)

 

F1 ← ecdf(x)

 

ecdf ← F1(c(x))

 

zweibullcdf ← 1 − ((exp($-g$ ∗ x^a)/(b^(1 − exp($-g$ ∗ x^a)))))

 

plot(x, ecdf, lty = 1, lwd = 2.5, type = “s,” xlab = “x,”

 

ylab = “Estimated cdf,” ylim = c(0,1), xlim = c(min(x),

 

max(x)), col = “black”)

 

par(new = TRUE)

 

plot(x, zweibullcdf, lty = 1, lwd = 3, type = “l,” xlab = “x,”

 

ylab = “Estimated cdf,” ylim = c(0,1), xlim = c(min(x),

 

max(x)), col = “blue”)

 

par(new = TRUE)

 

%%%%---------------------------------------------------

 

%Kaplan-Meier survival plot

 

%%%%---------------------------------------------------

 

x ← read.csv(file.choose(), header = TRUE)

 

x = x[,1]

 

x = x[!is.na(x)]

 

x = x/100

 

x

 

###—Parameter values

 

library(survival)

 

delta = rep(1,length(x))

 

x ← sort(x)

 

km = survfit(Surv(x,delta)∼1)

 

plot(km,conf.int = FALSE, ylab = “Kaplan-Meier survival plot,”xlab = “x”)

 

a = 2.346

 

$g$ = 0.564

 

b = 2.295

 

ss ← function(x)

 

{

 

(exp($-g$ ∗ x^a)/(b^(1 − exp($-g$ ∗ x^a))))

 

}

 

lines(seq(0, 3.5, length.out = 100),ss(seq(0, 3.5, length.out = 100)),

 

col = “blue,” lwd = 3)

 

%%%%---------------------------------------------------

 

%PP plot

 

%%%%---------------------------------------------------

 

x ← read.csv(file.choose(), header = TRUE)

 

x = x[,1]

 

x = x[!is.na(x)]

 

x = x/100

 

x

 

cdfLD1 = function(x, a, $g$, b)

 

{

 

1 − (exp($-g$ ∗ x^a)/(b^(1 − exp($-g$ ∗ x^a))))

 

}

 

x = sort(x)

 

n = length(x)

 

###—Empirical distribution function

 

Fn = seq(1,n)/n

 

plot(Fn, cdfLD1(x, 2.346, 0.564, 2.295), xlab = “x,”

 

ylab = “PP Plot,” pch = 21, col = “blue,” bg = “blue”)

 

abline(0,1)

 

%%%%---------------------------------------------------

 

%QQ plot

 

%%%%---------------------------------------------------

 

x ← read.csv(file.choose(), header = TRUE)

 

x = x[,1]

 

x = x[!is.na(x)]

 

x = x/100

 

x

 

###—Parameter values

 

a = 2.346

 

$g$ = 0.564

 

b = 2.295

 

x ← sort(x)

 

F1 ← ecdf(x)

 

ecdf ← F1(c(x))

 

zweibullcdf ← 1 − (exp($-g$ ∗ x^a)/(b^(1 − exp($-g$ ∗ x^a))))

 

qqnorm(x, pch = “@,” col = “black,” main = ””)

 

qqline(x, pch = “@,” col = “blue,” ylab = “X,” lty = 1, lwd = 3)

 

%%%%---------------------------------------------------

 

%Box plot

 

%%%%---------------------------------------------------

 

x < read.csv(file.choose(), header = TRUE)

 

x = x[,1]

 

x = x[!is.na(x)]

 

x = x/100

 

x

 

boxplot(x, main = ,” col = “blue,” ylab = “x,” xlab = “Box Plot”)

 

%%%%---------------------------------------------------

References

[1] M. Guillen, F. Prieto, J. M. Sarabia, "Modelling losses and locating the tail with the Pareto positive stable distribution," Insurance: Mathematics and Economics, vol. 49 no. 3, pp. 454-461, DOI: 10.1016/j.insmatheco.2011.07.004, 2011.

[2] D. Bhati, S. Ravi, "On generalized log-Moyal distribution: a new heavy tailed size distribution," Insurance: Mathematics and Economics, vol. 79, pp. 247-259, DOI: 10.1016/j.insmatheco.2018.02.002, 2018.

[3] J. Beirlant, G. Matthys, G. Dierckx, "Heavy-tailed distributions and rating," ASTIN Bulletin, vol. 31 no. 1, pp. 37-58, DOI: 10.2143/ast.31.1.993, 2001.

[4] Z. Ahmad, E. Mahmoudi, G. G. Hamedani, O. Kharazmi, "New methods to define heavy-tailed distributions with applications to insurance data," Journal of Taibah University for Science, vol. 14 no. 1, pp. 359-382, DOI: 10.1080/16583655.2020.1741942, 2020.

[5] M. Eling, "Fitting insurance claims to skewed distributions: are the skew-normal and skew-student good models?," Insurance: Mathematics and Economics, vol. 51 no. 2, pp. 239-248, DOI: 10.1016/j.insmatheco.2012.04.001, 2012.

[6] C. Adcock, M. Eling, N. Loperfido, "Skewed distributions in finance and actuarial science: a review," The European Journal of Finance, vol. 21 no. 13-14, pp. 1253-1281, DOI: 10.1080/1351847x.2012.720269, 2015.

[7] T. Shushi, "Skew-elliptical distributions with applications in risk theory," European Actuarial Journal, vol. 7 no. 1, pp. 277-296, DOI: 10.1007/s13385-016-0144-9, 2017.

[8] A. Punzo, A. Mazza, A. Maruotti, "Fitting insurance and economic data with outliers: a flexible approach based on finite mixtures of contaminated gamma distributions," Journal of Applied Statistics, vol. 45 no. 14, pp. 2563-2584, DOI: 10.1080/02664763.2018.1428288, 2018.

[9] A. Azzalini, T. Del Cappello, S. Kotz, "Log-skew-normal and log-skew- t distributions as models for family income data," Journal of Income Distribution, vol. 11 no. 3-4, 2002.

[10] L. Bagnato, A. Punzo, "Finite mixtures of unimodal beta and gamma densities and the k -bumps algorithm," Computational Statistics, vol. 28 no. 4, pp. 1571-1597, DOI: 10.1007/s00180-012-0367-4, 2013.

[11] G. A. Paula, V. Leiva, M. Barros, S. Liu, "Robust statistical modeling using the Birnbaum-Saunders- t distribution applied to insurance," Applied Stochastic Models in Business and Industry, vol. 28 no. 1, pp. 16-34, DOI: 10.1002/asmb.887, 2012.

[12] S. A. Klugman, H. H. Panjer, G. E. Willmot, Loss Models: from Data to Decisions, vol. 715, 2012.

[13] S. Nadarajah, S. A. A. Bakar, "New composite models for the Danish fire insurance data," Scandinavian Actuarial Journal, vol. 2014 no. 2, pp. 180-187, DOI: 10.1080/03461238.2012.695748, 2014.

[14] S. A. Bakar, N. A. Hamzah, M. Maghsoudi, S. Nadarajah, "Modeling loss data using composite models," Insurance: Mathematics and Economics, vol. 61, pp. 146-154, DOI: 10.1016/j.insmatheco.2014.08.008, 2015.

[15] A. Punzo, L. Bagnato, A. Maruotti, "Compound unimodal distributions for insurance losses," Insurance: Mathematics and Economics, vol. 81, pp. 95-107, DOI: 10.1016/j.insmatheco.2017.10.007, 2018.

[16] A. Mazza, A. Punzo, "Modeling household income with contaminated unimodal distributions," New Statistical Developments in Data Science, pp. 373-391, DOI: 10.1007/978-3-030-21158-5_28, 2019.

[17] M. Bernardi, A. Maruotti, L. Petrella, "Skew mixture models for loss distributions: a Bayesian approach," Insurance: Mathematics and Economics, vol. 51 no. 3, pp. 617-623, DOI: 10.1016/j.insmatheco.2012.08.002, 2012.

[18] T. Miljkovic, B. Grün, "Modeling loss data using mixtures of distributions," Insurance: Mathematics and Economics, vol. 70, pp. 387-396, DOI: 10.1016/j.insmatheco.2016.06.019, 2016.

[19] K. Dutta, J. Perry, "A tale of tails: an empirical analysis of loss distribution models for estimating operational risk capital," 2006. Working Papers, No. 06-13

[20] S. Nasiru, "Extended odd Fréchet-G family of distributions," Journal of Probability and Statistics, vol. 2018,DOI: 10.1155/2018/2931326, 2018.

[21] M. A. Cortés, D. Elal-Olivero, J. F. Olivares-Pacheco, "A new class of distributions generated by the extended bimodal-normal distribution," Journal of Probability and Statistics, vol. 2018,DOI: 10.1155/2018/9753439, 2018.

[22] H. Al-Mofleh, "On generating a new family of distributions using the tangent function," Pakistan Journal of Statistics and Operation Research, vol. 14 no. 3, pp. 471-499, DOI: 10.18187/pjsor.v14i3.1472, 2018.

[23] M. E. Mead, G. M. Cordeiro, A. Z. Afify, H. A. Mofleh, "The alpha power transformation family: properties and applications," Pakistan Journal of Statistics and Operation Research, vol. 15, pp. 525-545, DOI: 10.18187/pjsor.v15i3.2969, 2019.

[24] W. He, Z. Ahmad, A. Z. Afify, H. Goual, "The arcsine exponentiated- X family: validation and insurance application," Complexity, vol. 2020,DOI: 10.1155/2020/8394815, 2020.

[25] Z. Ahmad, E. Mahmoudi, S. Dey, S. K. Khosa, "Modeling vehicle insurance loss data using a new member of T-X family of distribution," Journal of Statistical Theory & Practice, vol. 19 no. 2, pp. 133-147, DOI: 10.2991/jsta.d.200421.001, 2020.

[26] Z. Ahmad, E. Mahmoudi, S. Dey, "A new family of heavy-tailed distribution with an application to the heavy-tailed insurance loss data," Communication in Statistics: Simulation & Computation, vol. 49,DOI: 10.1080/03610918.2020.1741623, 2020.

[27] A. Alzaatreh, C. Lee, F. Famoye, "A new method for generating families of continuous distributions," Metron, vol. 71 no. 1, pp. 63-79, DOI: 10.1007/s40300-013-0007-y, 2013.

[28] Z. Ahmad, G. G. Hamedani, N. S. Butt, "Recent developments in distribution theory: a brief survey and some new generalized classes of distributions," Pakistan Journal of Statistics and Operation Research, vol. 15 no. 1, pp. 87-110, DOI: 10.18187/pjsor.v15i1.2803, 2019.

[29] P. Artzner, "Application of coherent risk measures to capital requirements in insurance," North American Actuarial Journal, vol. 3 no. 2, pp. 11-25, DOI: 10.1080/10920277.1999.10595795, 1999.

[30] W. Stute, W. G. Manteiga, M. P. Quindimil, "Bootstrap based goodness-of-fit-tests," Metrika, vol. 40 no. 1, pp. 243-256, DOI: 10.1007/bf02613687, 1993.

[31] D. Kundu, H. Howlader, "Bayesian inference and prediction of the inverse Weibull distribution for type-II censored data," Computational Statistics & Data Analysis, vol. 54 no. 6, pp. 1547-1558, DOI: 10.1016/j.csda.2010.01.003, 2010.

[32] S. Dey, T. Dey, "On progressively censored generalized inverted exponential distribution," Journal of Applied Statistics, vol. 41 no. 12, pp. 2557-2576, DOI: 10.1080/02664763.2014.922165, 2014.

[33] S. Dey, S. Ali, C. Park, "Weighted exponential distribution: properties and different methods of estimation," Journal of Statistical Computation and Simulation, vol. 85 no. 18, pp. 3641-3661, DOI: 10.1080/00949655.2014.992346, 2015.

[34] S. Dey, T. Dey, S. Ali, M. S. Mulekar, "Two-parameter Maxwell distribution: properties and different methods of estimation," Journal of Statistical Theory and Practice, vol. 10 no. 2, pp. 291-310, DOI: 10.1080/15598608.2015.1135090, 2016.

[35] R. Calabria, G. Pulcini, "Point estimation under asymmetric loss functions for left truncated exponential samples," Communications in Statistics-Theory and Methods, vol. 25 no. 3, pp. 585-600, DOI: 10.1080/02664763.2018.1542668, 1996.

[36] S. Dey, S. Singh, Y. M. Tripathi, A. Asgharzadeh, "Estimation and prediction for a progressively censored generalized inverted exponential distribution," Statistical Methodology, vol. 32, pp. 185-202, DOI: 10.1016/j.stamet.2016.05.007, 2016.

[37] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, E. Teller, "Equation of state calculations by fast computing machines," The Journal of Chemical Physics, vol. 21 no. 6, pp. 1087-1092, DOI: 10.1063/1.1699114, 1953.

[38] W. K. Hastings, "Monte Carlo sampling methods using Markov chains and their applications," Biometrika, vol. 57 no. 1, pp. 97-109, DOI: 10.1093/biomet/57.1.97, 1970.

[39] S. Geman, D. Geman, "Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images," IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 6 no. 6, pp. 721-741, DOI: 10.1109/tpami.1984.4767596, 1984.

[40] P. Congdon, Bayesian Statistical Modelling, 2001.