A New Generalized t Distribution Based on a Distribution Construction Method

Abstract

In this paper, a new generalized t (new Gt) distribution based on a distribution construction approach is proposed and proved to be suitable for fitting both the data with high kurtosis and heavy tail. The main innovation of this article consists of four parts. First of all, the main characteristics and properties of this new distribution are outined. Secondly, we derive the explicit expression for the moments of order statistics as well as its corresponding variance–covariance matrix. Thirdly, we focus on the parameter estimation of this new Gt distribution and introduce several estimation methods, such as a modified method of moments (MMOM), a maximum likelihood estimation (MLE) using the EM algorithm, a novel iterative algorithm to acquire MLE, and improved probability weighted moments (IPWM). Through simulation studies, it can be concluded that the IPWM estimation performs better than the MLE using the EM algorithm and the MMOM in general. The newly-proposed iterative algorithm has better performance than the EM algorithm when the sample kurtosis is greater than 2.7. For four parameters of the new Gt distribution, a profile maximum likelihood approach using the EM algorithm is developed to deal with the estimation problem and obtain acceptable.

1. Introduction

The generalized t (Gt) distribution was first proposed by McDonald and Newey (1988) [1] to implement a partially adaptive regression models with the following probability density function (pdf)

$$\begin{array}{c}\hfill {\displaystyle f(x;p,q)}={\displaystyle \frac{p}{2q^{1/p}B(1/p,q)}{\left(1+\frac{{\left|x\right|}^p}{q}\right)}^{-(q+1/p)},p,q>0,}\end{array}$$

(1)

where $B(a,b)$ denotes the beta function. With the addition of two shape parameters, $p,q$, the shape of the distribution becomes more flexibe so that it can fit a wider range of data. Nadarajah (2008) [2] obtained the analytic expression of the cumulative distribution function (cdf) through the incomplete beta function presented below

$$\begin{array}{c}\hfill {\displaystyle F(x;p,q)}={\displaystyle \frac{1}{2}\left[1+\mathrm{sign}\left(x\right)I_{{1-(1+|x|}^p{/q)}^{-1}}\left(\frac{1}{p},q\right)\right]},\end{array}$$

(2)

where $I_x(a,b)=\frac{1}{B(a,b)}{\int }_0^x{t}^{a-1}{(1-t)}^{b-1}\mathsf{d}t$ denotes the incomplete beta function.

With the passage of time, the generalized t distribution warrants special attention by scholars and has been widely applied in several research areas. Galbraith and Zhu (2010) [3] introduced a new class of asymmetric Student t distributions and illustrated their applications in financial econometrics. Harvey and Lange (2017) [4] applied the generalized t distribution in autoregressive conditional heteroscedasticity models. Furthermore, the generalized t distribution and its extension were involved in extensive research in the fields of robust estimation and robust statistical models, such as [5,6,7,8].

As the types of data become more and more widespread, many different generalized t distributions have been proposed and the related properties have been studied in depth. Arslan and Genç (2009) [9] obtained the skew-generalized t distribution through the scale mixture of the skew-exponential power and generalized gamma distributions. Venegas et al. (2012), [10] constructed a more flexible distribution by introducing skew coefficients $\epsilon$ to fit a wider range of data. Papastathopoulos and and Tawn (2013) [11] coined an extension of the Student’s t distribution that allows for a negative shape parameter. Acitas et al. (2015), [12] introduced an alpha-skew-generalized t distribution based on a new skewing procedure, and presented the corresponding properties as well as the estimation methods. Lak et al. (2019), [13] proposed a distribution with a variety of shape variations, namely, the alpha-beta-skew-generalized t distribution, which contains several well-known distributions and illustrated its usefulness by means of a real data analysis. For Student’s t Birnbaum-Saunders distribution, we recommend [14,15].

Although there is plenty of literature on the generalized t distribution, two aspects still need to be improved: (1) the type of data that previous distributions can be fitted to is not wide enough; (2) for the flexibility of the shape, the distribution function is often very complex, imposing limitations on the parameter estimation. To solve these two problems, the new Gt distribution we propose will fit both the data with a heavier tail and with higher kurtosis based on different shape parameters, with a more extensive application scope than the generalized t distributions introduced before. In addition, our new Gt distribution will make many classical estimation methods suitable for the parameter estimation problem of this distribution.

In this paper, we coined a new Gt distribution based on the distribution construction approach and denote it by $N\mathrm{Gt}(\mu ,\sigma ,\alpha ,\beta )$. The main innovation of this work consists of four parts. Firstly, we investigate the main properties of the new Gt distribution, including moments, skewness coefficients, kurtosis coefficients and random number generation. Secondly, we derive explicit expressions for moments of order statistics as well as its corresponding variance–covariance matrix through recurrence relations and distribution transformation techniques when the shape parameter $\frac{1}{\beta }$ is either a real non-integer or an integer. Thirdly, we established an EM-type algorithm and a novel iterative algorithm to acquire the maximum likelihood estimation (MLE) of parameters for $N\mathrm{Gt}(0,\sigma ,\alpha ,\beta )$. The explicit analytical expressions for the asymptotic variance–covariance matrix of MLE were also presented based on the EM algorithm. Finally, we introduce improved probability-weighted moments (IPWM) based on order statistics transformation and obtain the corresponding algorithm to implement this novel method. For four parameters, $N\mathrm{Gt}(\mu ,\sigma ,\alpha ,\beta )$, we propose a profile maximum likelihood approach (PLA) using the EM algorithm to deal with the estimation problem.

The rest of this article is organized as follows. In Section 2, we investigate the basic properties of this new distribution and make a comparison with the Gt distribution introduced by McDonald and Newey (1988) [1]. In Section 3, we derive explicit expressions for moments of order statistics from this distribution. In Section 4, we focus on the variance–covariance matrix of order statistics. In Section 5, several estimation methods are considered. In Section 6, we conduct several simulation studies to illustrate the performance of the proposed estimation method and compare these simulation results under different parameter values and different sample sizes. In Section 7, a real data analysis is performed to illustrate the feasibility of the new distribution we proposed. Discussions and conclusions are presented in Section 8 and proofs are given in the Appendix A.

2. Basic Theorems and Properties

In this section, some basic definitions, theorems and properties for the new Gt distribution are elucidated.

The new Gt distribution contains many common distributions, such as the normal distribution, the Laplace distribution, and the Pareto distribution. With the addition of two shape parameters, the shape of the distribution becomes more flexible, which can fit data having both heavy tails and high kurtosis.

Definition 1.

Let the random variable X with the following pdf

$$f_X(x;\mu ,\sigma ,\beta )=\frac{\beta }{2^{1+1/\beta }\Gamma (1/\beta )\sigma }\exp \left\{-\frac{{|x-\mu |}^{\beta }}{2{\sigma }^{\beta }}\right\},x\in R,$$

where $\Gamma ( )$ denotes the gamma function, $\mu \in R$ is a location parameter, $\sigma >0$ is a scale parameter, and $\beta >0$ is a shape parameter. Then the random variable X is said to have a generalized normal distribution and denoted by $X\sim GN(\mu ,\sigma ,\beta )$.

The definition of the new Gt distribution (denoted by $N\mathrm{Gt}(\mu ,\sigma ,\alpha ,\beta ))$ is stated as follows:

Definition 2.

Let the random variable $X-\mu \sim GN(0,\sigma ,\beta )$ and $Y\sim Gamma(\alpha ,\beta )$. Supposing X and Y are independent, we shall define $T=X/Y^{\frac{1}{\beta }}\sim N\mathrm{Gt}(\mu ,\sigma ,\alpha ,\beta )$, where $\mu \in R$ is a location parameter, $\sigma >0$ is a scale parameter, and $\alpha ,\beta >0$ are two shape parameters.

Theorem 1.

The pdf of $N\mathrm{Gt}(\mu ,\sigma ,\alpha ,\beta )$ can be expressed as follows

$$\begin{array}{ccc}\hfill {\displaystyle g_T(t;\mu ,\sigma ,\alpha ,\beta )}& =& {\displaystyle \frac{1}{{\left(\frac{2}{\beta }\right)}^{1+\frac{1}{\beta }}B(\alpha ,\frac{1}{\beta })\sigma }{\left[1+\frac{{\beta |t-\mu |}^{\beta }}{2{\sigma }^{\beta }}\right]}^{-(\alpha +\frac{1}{\beta })}},{\displaystyle t\in R.}\hfill \end{array}$$

Proof. 

See Appendix A.1 (Part I).   □

As is shown in the following two pictures (Figure 1 and Figure 2), two shape parameters ($\alpha$ and $\beta$) have a decisive effect on the shape of the density function. The density function has a heavier tail with a smaller value of $\alpha$ and $\beta$, but it is the opposite when the parameters $\alpha$ and $\beta$ are larger. Furthermore, the shape parameter $\beta$ plays a decisive role in the shape of the distribution and the change of parameter $\alpha$ can control the tail behavior.

Figure 3 and Figure 4 vividly describe the feature that the new Gt distribution is more suitable for fitting both the data with high kurtosis and with heavier tail than the Gt distribution coined by McDonald and Newey (1988) [1] under appropriate parameter range from two perspectives respectively.

It is well-known that the pdf with a steeper shape is suitable for fitting data with high kurtosis, while a pdf with a smoother shape can be used for fitting heavy-tailed data. Since these two distributions we compare are symmetric, we use the maximum value of the pdf to characterize the steepness of its shape. Let $K_1(\sigma ,\alpha ,\beta )$ and $K_2(\sigma ,p,q)$ denote the maximum value of the pdf of Gt distribution and new Gt distribution, respectively.

Since the value of $\sigma$ has little effect on the value of $K_1(\sigma ,\alpha ,\beta )$ and $K_2(\sigma ,p,q)$, we may assume that $\sigma =1$. On the one hand, it is not very difficult to prove that, for a fixed $\beta$, $K_1(\sigma ,\alpha ,\beta )$ decreases monotonically with respect to $\alpha$. Then we have: $\underset{\alpha ⟶\infty }{lim}{K}_1(\sigma ,\alpha ,\beta )⟶+\infty ,\underset{\alpha ⟶0^{+}}{lim}{K}_1(\sigma ,\alpha ,\beta )⟶0$. For $K_2(\sigma ,p,q)$, the interval of the shape parameters $p,q$ are (1, 10) and (0.1, 50), respectively. The reason for choosing these two ranges is that they include possible values for parameters in practical applications. From Figure 3, we can obviously find that the range of function $K_2(\sigma ,p,q)$ is never greater than 1. This means that we can only adjust the parameter $\sigma$ to make the pdf of the Gt distribution have a steeper shape. On the other hand, for fixed $\beta$ and p, when $\alpha ,q⟶0^{+}$, both $K_1(\sigma ,\alpha ,\beta )$ and $K_2(\sigma ,p,q)$ tend to 0 but $K_1(\sigma ,\alpha ,\beta )$ has a faster convergence. This means that the pdf of the new Gt distribution proposed by us has a heavier tail when the parameters are within a reasonable range. Thus, the conclusion stated in the preceding paragraph is therefore self-evident.

The pdf generated under different shape parameters can be divided into four categories according to different shapes and the tail behavior, i.e., steep in shape and with a light tail, flat in shape and with a light tail, steep in shape and with a heavy tail, flat in shape and with a heavy tail. As is depicted in Figure 4, the solid line denotes the new Gt distribution, and the dashed line represents the previous Gt distribution defined by McDonald and Newey (1988) [1]. The subfigure Figure 4a–d depicts the pdf of the Gt distribution and the new Gt distribution under the above four categories with different shape parameters, respectively. It can be seen from the figure that, if the parameter value is in a reasonable range, the shape of the new Gt distribution is more flexible and more suitable for fitting a wider range of data.

Some important properties of the new Gt distribution are stated as follows:

Proposition 1.

If $Y|U=u\sim GN(\mu ,\sigma u^{-\frac{1}{\beta }},\beta )$, and $U\sim Gamma(\alpha ,\beta )$, then $Y\sim N\mathrm{Gt}(\mu ,\sigma ,\alpha ,\beta )$. $U|Y=y\sim Gamma(\alpha +\frac{1}{\beta },\frac{1}{\beta }+\frac{{|y-\mu |}^{\beta }}{2{\sigma }^{\beta }})$.

Proof. 

See Appendix A.1 (Part I).   □

Proposition 2.

If $Y|U=u\sim GN(\mu ,\sigma {\left(\frac{\beta }{2}u\right)}^{-\frac{1}{\beta }},\beta )$ and $U\sim {\chi }^2\left(2\alpha \right)$, then $Y\sim N\mathrm{Gt}(\mu ,\sigma ,\alpha ,\beta )$. $U|Y=y\sim Gamma(\alpha +\frac{1}{\beta },\frac{1}{2}+\frac{{\beta |y-\mu |}^{\beta }}{4{\sigma }^{\beta }})$.

Proposition 3.

If $X\sim Beta(\frac{1}{\beta },\alpha )$, let $Y=\mu +\sigma ·W·{\left(2Z\right)}^{\frac{1}{\beta }}$ and $Z=\frac{X}{\beta (1-X)}$, then $Y\sim NGt(\mu ,\sigma ,\alpha ,\beta )$, where W is a random variable that only takes $\pm 1$ and satisfies $P(W=1)=P(W=-1)=\frac{1}{2}$.

From the above proposition, we can follow the following steps to generate random samples from $N\mathrm{Gt}(\mu ,\sigma ,\alpha ,\beta )$:

1.

Generate $X\sim Beta(1/\beta ,\alpha )$ and set $Z=\frac{X}{\beta (1-X)}$.

2.

Generate $U\sim U(0,1)$.

if $U\le \frac{1}{2}$, set $W=1$.

else set $W=-1$.

3.

Set $Y=\mu +\sigma W{\left(2Z\right)}^{\frac{1}{\beta }}$ and return Y.

Proposition 4.

If $T\sim N\mathrm{Gt}(\sigma ,\alpha ,\beta )$, we have

• the $2k$-th order moment of T is given by

  $$E\left(X^{2k}\right)={\left(\frac{2}{\beta }\right)}^{\frac{2k}{\beta }}\frac{\Gamma (\alpha -2k/\beta )\Gamma \left[\right(2k+1)/\beta ]}{\Gamma \left(\alpha \right)\Gamma (1/\beta )}{\sigma }^{2k},k=1,2,\cdots ,$$

  and the above equation exists for $\alpha \beta >2k$. Furthermore, from the symmetry of the distribution, we can conclude that $E\left(X^{2k-1}\right)=0$, $k=1,2,\dots ,$
• the kurtosis coefficient of T is given by

  $$G_2(\alpha ,\beta )=\frac{\Gamma \left(\alpha \right)\Gamma (\alpha -4/\beta )\Gamma (1+1/\beta )\Gamma (5/\beta )}{{\Gamma }^2(\alpha -2/\beta ){\Gamma }^2(3/\beta )}-3\alpha \beta >4,$$
• the mean deviation of T is given by

  $$\delta =\frac{2\sigma B(2/\beta ,\alpha -1/\beta )}{{(2/\beta )}^{1-1/\beta }B(\alpha ,1/\beta )\beta },\alpha \beta >1,$$
• the shannon entropy is given by

  $$E[-log\left(t\right)]=\left(1+\frac{1}{\beta }\right)log\left(\frac{2}{\beta }\right)+logB\left(\alpha ,\frac{1}{\beta }\right)+log\sigma -\left(\alpha +\frac{1}{\beta }\right)\left[\phi \left(\alpha \right)-\phi (\alpha +1/\beta )\right],$$
• the Rényi entropy is given by

  $${{\rm Y}}_R\left(\eta \right)=\frac{2B(1/\beta ,\eta \alpha +(\eta -1)/\beta )}{(1-\eta ){(2/\beta )}^{\eta +(\eta -1)/\beta }{\left[B(\alpha ,1/\beta )\right]}^{\eta }{\sigma }^{\eta -1}\beta },$$

  where $\phi ( )$ denote the digamma function.

3. Moments of Order Statistics

In this section, we derive the explicit expression for moments of order statistics from this distribution under the independent identically-distributed (IID) case and independent not identically-distributed (INID) case when the shape parameters $\frac{1}{\beta }$ and $\alpha$ are both real non-integer or $\frac{1}{\beta }$ is an integer.

It is well-known that order statistics has always been a significant field in statistical research, especially in the study of distribution theory. In addition, it is necessary to study the related properties of order statistics, since many classical parameter estimation methods are based on the application of order statistics.

3.1. The IID Case

In this subsection, we obtain explicit expression for moments of order statistics from $N\mathrm{Gt}(\mu ,\sigma ,\alpha ,\beta )$ under IID case based on distribution transformation approach and recurrence relations when the shape parameter $\frac{1}{\beta }$ is either an integer or a real non-integer.

3.1.1. Shape Parameter $\frac{1}{\beta }$ Integer Case

Govindarajulu (1963) [16] pointed out that the k-th order moment of order statistics from a symmetric distribution can be linearly expressed by the moments of the order statistics from the corresponding folded distribution

$$\begin{array}{c}\hfill {\displaystyle 2^n{\delta }_{r:n}^{\left(k\right)}}={\displaystyle \sum _{m=0}^{r-1}\left(\genfrac{}{}{0pt}{}{n}{m}\right){\gamma }_{r-m:n-m}^{\left(k\right)}+{(-1)}^k\sum _{m=r}^n\left(\genfrac{}{}{0pt}{}{n}{m}\right){\gamma }_{m-r+1:m}^{\left(k\right)}},\end{array}$$

where ${\delta }_{r:n}^{\left(k\right)}$ denotes the k-th moments of order statistics from symmetric distribution, while ${\gamma }_{r:n}^{\left(k\right)}$ denotes the k-th moments of order statistics from the corresponding folded distribution.

Without loss of generality, let $\mu =0,\sigma =1$. Making transformation $\frac{{\left|T\right|}^{\beta }}{2}=Z$, we can derive the pdf of the random variable Z as follows

$$\begin{array}{ccc}\hfill {\displaystyle f_Z(z;\alpha ,\beta )}& =& {\displaystyle \frac{1}{{\left(\frac{1}{\beta }\right)}^{\frac{1}{\beta }}B(\alpha ,\frac{1}{\beta })}{z}^{\frac{1}{\beta }-1}{(1+\beta z)}^{-(\alpha +\frac{1}{\beta })}},{\displaystyle z>0}\hfill \end{array}$$

(3)

Under the above transformation, the recurrence relations introduced by Govindarajulu (1963) [16] becomes

Lemma 1.

$$\begin{array}{c}\hfill {\displaystyle 2^{n-\frac{k}{\beta }}{\mu }_{r:n}^{\left(k\right)}}={\displaystyle \sum _{m=0}^{r-1}\left(\genfrac{}{}{0pt}{}{n}{m}\right){\nu }_{r-m:n-m}^{\left(k\right)}+{(-1)}^k\sum _{m=r}^n\left(\genfrac{}{}{0pt}{}{n}{m}\right){\nu }_{m-r+1:m}^{\left(k\right)}},\end{array}$$

(4)

where ${\mu }_{r:n}^{\left(k\right)}$ denotes the k-th moments of order statistics from $N\mathrm{Gt}(0,1,\alpha ,\beta )$ and ${\nu }_{r:n}^{\left(k\right)}$ denotes the k-th moments of order statistics from the corresponding folded distribution $F_Z\left(z\right)$.

Proof. 

See Appendix A.2 (Part II).   □

Let ${\alpha }_{r:n}^{\left(k\right)}$ denote $E\left(Z_{r:n}^k\right)$, then ${\nu }_{r:n}^{\left(k\right)}={\alpha }_{r:n}^{(k/\beta )}$. The cdf of the random variable Z can be written as

$$\begin{array}{ccc}\hfill {\displaystyle F_Z(z;\alpha ,\beta )}& =& {\displaystyle \frac{1}{{\left(\frac{1}{\beta }\right)}^{\frac{1}{\beta }}B\left(\alpha ,\frac{1}{\beta }\right)}{\int }_0^z{y}^{\frac{1}{\beta }-1}{(1+\beta y)}^{-(\alpha +\frac{1}{\beta })}\mathsf{d}y}\hfill \\ & =& {\displaystyle \frac{1}{B\left(\alpha ,\frac{1}{\beta }\right)}{\int }_0^{1-\frac{1}{1+\beta z}}{y}^{\frac{1}{\beta }-1}{(1-y)}^{\alpha -1}\mathsf{d}y},\hfill \end{array}$$

by binomial expansion we have

$$\begin{array}{c}\hfill 1-F_Z\left(z\right)=\frac{1}{B\left(\alpha ,\frac{1}{\beta }\right)}\sum _{i=0}^{\frac{1}{\beta }-1}{C}_i\frac{1}{{(1+\beta z)}^{\alpha +i}},\end{array}$$

(5)

where $C_i={(-1)}^i\left(\genfrac{}{}{0pt}{}{1/\beta -1}{i}\right)\frac{1}{\alpha +i}$.

Let $H_v(\frac{1}{\beta },p)$ denote the coefficient in the expansion of ${\left[{\displaystyle \sum _{i=0}^{\frac{1}{\beta }-1}}{C}_i\frac{1}{{(1+\beta y)}^i}\right]}^p$. Obviously, the coefficients of this polynomial satisfy the following recursive relation

$$\begin{array}{c}\hfill {\displaystyle H_v\left(\frac{1}{\beta },p\right)}={\displaystyle \sum _{i=0}^v{C}_i{H}_{v-i}\left(\frac{1}{\beta },p-1\right)},\end{array}$$

and with the following constraints

$$H_v\left(\frac{1}{\beta },1\right)=C_v,H_v\left(\frac{1}{\beta },0\right)=1,$$

we can cauculate every coefficient striaghtly. Then ${[1-F_Z\left(z\right)]}^p$ can be expressed as

$$\begin{array}{c}\hfill {[1-F_Z\left(z\right)]}^p=\frac{1}{{\left[B\left(\alpha ,\frac{1}{\beta }\right)\right]}^p}{\displaystyle \sum _{v=0}^{(\frac{1}{\beta }-1)p}}{H}_v\left(\frac{1}{\beta },p\right)\frac{1}{{(1+\beta z)}^{\alpha p+v}}.\end{array}$$

(6)

It is noted that Equations (5) and (6) have similar forms, so the distribution $F_Z\left(z\right)$ is said to have a minimum domains of attraction. With this conclusion, $E\left(Z_{1:n}^k\right)$ is derived first. Through the traditional way, we obtain

$$\begin{array}{ccc}\hfill {\displaystyle E\left(Z_{1:n}^k\right)}& =& {\displaystyle n{\int }_0^{+\infty }{z}^k{[1-F_Z\left(z\right)]}^{n-1}f\left(z\right)\mathsf{d}z}\hfill \\ & =& {\displaystyle \frac{n{\left(\frac{1}{\beta }\right)}^k}{{\left[B\left(\alpha ,\frac{1}{\beta }\right)\right]}^n}\sum _{v=0}^{(\frac{1}{\beta }-1)(n-1)}{H}_v\left(\frac{1}{\beta },n-1\right)\times B\left(\frac{1}{\beta }+k,\alpha n+v-k\right)},\hfill \end{array}$$

and the above expression exists for $k<\alpha n$.

After that, using the well-known recurrence relations from David and Nagaraja (2003) [17] presented below

$$\begin{array}{c}\hfill {\displaystyle (n-r){\nu }_{r:n}^{\left(k\right)}+r{\nu }_{r+1:n}^{\left(k\right)}}={\displaystyle n{\nu }_{r:n-1}^{\left(k\right)}},\end{array}$$

(7)

where $r=1,\cdots ,n-1,$ and $k=1,2,\cdots .$

We acquire the k-th order moment of order statistics from the distribution $F_Z(z;\alpha ,\beta )$. When the sample size is n, the method of obtaining all the k-th moments of order statistics from $N\mathrm{Gt}(0,1,\alpha ,\beta )$ when $\frac{1}{\beta }$ is an integer can be briefly described as follows

1.

Derive ${\nu }_{i,j}^{\left(k\right)}$ ($1\le i\le n,1\le j\le n$) by using the above approach.

2.

For $1\le r\le \left[\frac{n}{2}\right]$, derive ${\mu }_{r:n}^{\left(k\right)}$ by using the recurrence relations in Equation (4).

3.

By using the symmetry of Gt distribution, the k-th moments of the remaining order statistics can be calculated through ${\mu }_{n-r+1:n}^{\left(k\right)}=-{\mu }_{r:n}^{\left(k\right)}$.

Since this method makes full use of the characteristics of the distribution function and recursive relations, it is very efficient and can avoid a lot of unnecessary calculations.

3.1.2. Shape Parameter $\frac{1}{\beta }$ or $\alpha$ Non-Integer Case

The generalized Kampe de Feriet function from Exton (1978) [18] is defined by

$$\begin{array}{cc}\hfill {\displaystyle F_{C:D}^{A:B}(\left(a\right):\left(b_1\right);\cdots ;\left(b_n\right);\left(c\right):\left(d_1\right);\cdots ,\left(d_n\right);x_1,\cdots ,x_n)}=& \sum _{m_1=0}^{\infty }\cdots \sum _{m_n=0}^{\infty }\frac{{\left(\left(a\right)\right)}_{m_1+\cdots +m_n}{\left(\left(b_1\right)\right)}_{m_1}\cdots {\left(\left(b_n\right)\right)}_{m_n}}{{\left(\left(c\right)\right)}_{m_1+\cdots +m_n}{\left(\left(d_1\right)\right)}_{m_1}\cdots {\left(\left(d_n\right)\right)}_{m_n}}\hfill \\ \hfill & \times \frac{x_1^{m_1}\cdots x_n^{m_n}}{m_1!\cdots m_n!},\hfill \end{array}$$

where $a=(a_1,a_2,\cdots ,a_A)$, $b_i=(b_{i,1},b_{i,2},\cdots ,b_{i,B})$ for $i=1,2,\cdots ,n$, $c=(c_1,\cdots ,c_C)$, $d_i=(d_{i,1},\cdots ,d_{i,D})$, or $i=1,2,\cdots ,n$ and ${\left(\left(f\right)\right)}_k={\left((f_1,f_2,\cdots ,f_p)\right)}_k={\left(f_1\right)}_k\cdots {\left(f_p\right)}_k$, ${\left(f_i\right)}_k=f_i(f_i+1)\cdots (f_i+k-1)$.

By using the above special function, we derive the following theorem:

Theorem 2.

If the shape parameter α is a real non-integer, the k-th order moment of order statistics $X_{r:n}$ from $N\mathrm{Gt}(0,1,\alpha ,\beta )$ can be calculated by the following convergent expression

$$E\left(X_{r:n}^k\right)=C_{r,n}\left[I_1(k,r-1,n-r)+I_2(k,r-1,n-r)\right],$$

where

$$I_1(k,r-1,n-r)=\frac{2^{\frac{k}{\beta }-n}{\left(\frac{1}{\beta }\right)}^{\frac{k}{\beta }}}{B\left(\frac{1}{\beta },\alpha \right)}\sum _{i=0}^{r-1}\sum _{j=0}^{n-r}{(-1)}^j\left(\genfrac{}{}{0pt}{}{r-1}{i}\right)\left(\genfrac{}{}{0pt}{}{n-r}{j}\right)A\left(\frac{k+1}{\beta },\alpha -\frac{k}{\beta }-1,i+j\right),$$

$$I_2(k,r-1,n-r)={(-1)}^k\frac{2^{\frac{k}{\beta }-n}{\left(\frac{1}{\beta }\right)}^{\frac{k}{\beta }}}{B\left(\frac{1}{\beta },\alpha \right)}\sum _{i=0}^{n-r}\sum _{j=0}^{r-1}{(-1)}^j\left(\genfrac{}{}{0pt}{}{n-r}{i}\right)\left(\genfrac{}{}{0pt}{}{r-1}{j}\right)A\left(\frac{k+1}{\beta },\alpha -\frac{k}{\beta }-1,i+j\right),$$

$$\begin{array}{cc}\hfill A\left(\frac{k+1}{\beta },\alpha -\frac{k}{\beta }-1,i+j\right)=& \frac{{\beta }^{i+j}B\left(\frac{k+i+j+1}{\beta },\alpha -\frac{k}{\beta }\right)}{{\left[B\left(\alpha ,\frac{1}{\beta }\right)\right]}^{i+j}}\times F_{1:1}^{1:2}(\left(\frac{k+i+j+1}{\beta }\right):(1-\alpha ,\frac{1}{\beta });\cdots ;\hfill \\ \hfill & (1-\alpha ,\frac{1}{\beta }):(\alpha +\frac{i+j+1}{\beta }):(\frac{1}{\beta }+1);\cdots ;(\frac{1}{\beta }+1):1;\cdots ;1\left)\right),\hfill \end{array}$$

$$C_{r:n}=\frac{n!}{(r-1)!(n-r)!}.$$

Proof. 

See Appendix A.2 (Part II).   □

Theorem 3.

If the shape parameter $\frac{1}{\beta }$ is a real non-integer, $A\left(\frac{k+1}{\beta },\alpha -\frac{k}{\beta }-1,i+j\right)$ in Theorem 2 can be calculated by

$$\begin{array}{ccc}\hfill {\displaystyle A(\frac{k+1}{\beta },\alpha -\frac{k}{\beta }-1,i+j)}& =& {\displaystyle \frac{{\left(\frac{1}{\beta }\right)}^{\frac{k+1}{\beta }}{\alpha }^{-(i+j)}}{{\left[B\left(\alpha ,\frac{1}{\beta }\right)\right]}^{i+j}}B\left(\frac{k+1}{\beta },\alpha (i+j+1)-\frac{k}{\beta }\right)\times F_{1:1}^{1:2}((\alpha (i+j+1)-\frac{k}{\beta }):}\hfill \\ & & {\displaystyle (1-\frac{1}{\beta },\alpha );\cdots ;(1-\frac{1}{\beta },\alpha ):(\alpha (i+j+1)+\frac{1}{\beta }):(\alpha +1);\cdots ;(\alpha +1):}\hfill \\ & & {\displaystyle 1;\cdots 1)},\hfill \end{array}$$

Proof. 

See Appendix A.2 (Part II).   □

3.2. The INID Case

Suppose that $X_1,\cdots ,X_n$ are independent but not identical random samples from the following pdf with the same shape parameter $\beta$ but the rest of the parameters are different

$$\begin{array}{c}\hfill {\displaystyle g_i\left(x\right)}={\displaystyle \frac{1}{{\left(\frac{2}{\beta }\right)}^{1+\frac{1}{\beta }}B\left({\alpha }_i,\frac{1}{\beta }\right){\sigma }_i}{\left[1+\frac{\beta |t-{\mu }_i{|}^{\beta }}{2{\sigma }_i^{\beta }}\right]}^{-({\alpha }_i+\frac{1}{\beta })}}.\end{array}$$

Theorem 4.

If the shape parameter $\frac{1}{\beta }$ is a real non-integer, the m-th order moment of order statistics $X_{r:n}$ from $N\mathrm{Gt}(0,1,{\alpha }_i,\beta )$ can be calculated by the following convergent expression

$${\mu }_{r:n}^m=\left[\sum _{j=n-r+1}^n{(-1)}^{j-(n-r+1)}\left(\genfrac{}{}{0pt}{}{j-1}{n-r}\right)+\sum _{j=r}^n{(-1)}^{m+j-r}\left(\genfrac{}{}{0pt}{}{j-1}{r-1}\right)\right]I_j^{+}\left(m\right),$$

where

$$\begin{array}{ccc}\hfill {\displaystyle I_j^{+}\left(m\right)}& =& {\displaystyle \frac{m{2}^{\frac{m}{\beta }-1-j}}{{\displaystyle \prod _{t=1}^j}{\alpha }_{i_t}B({\alpha }_{i_t},\frac{1}{\beta })}{\left(\frac{1}{\beta }\right)}^{\frac{m}{\beta }-1}B\left(\frac{m}{\beta },\sum _{t=1}^j{\alpha }_{i_t}-\frac{m}{\beta }\right)\times F_{1:1}^{1:2}((\sum _{t=1}^j{\alpha }_{i_t}-\frac{m}{\beta }):(1-\frac{1}{\beta },{\alpha }_{i_1});\cdots ;}\hfill \\ & & {\displaystyle (1-\frac{1}{\beta },{\alpha }_{i_j}):\left(\sum _{t=1}^j{\alpha }_{i_t}\right):({\alpha }_{i_1}+1);\cdots ;({\alpha }_{i_j}+1):1;\cdots ;1),{\alpha }_{i_t}>\frac{m}{\beta }.}\hfill \end{array}$$

(8)

Proof. 

See Appendix A.2 (Part II).   □

4. Product Moments and Variance–Covariance Matrix of Order Statistics

In this section, we obtain explicit expressions for product moments and the corresponding variance–covariance matrix of order statistics from this distribution based on recurrence relations and distribution transformation when the shape parameter $\frac{1}{\beta }$ is an integer.

The covariance matrix of order statistics plays a big role in a number of estimation methods based on order statistics, such as the least square estimation and the BLUE for location-scale parameters. In addition, obtaining a covariance matrix of order statistics is usually not an easy task and requires some related skills.

Under transformation $\frac{{\left|T\right|}^{\beta }}{2}=Z$ introduced in the previous section, the recursive relations for the product moment of order statistics from a symmetric distribution between its corresponding folded distribution becomes

$$\begin{array}{ccc}\hfill {\displaystyle 2^{n-\frac{k_1+k_2}{\beta }}{\mu }_{r,s:n}^{(k_1,k_2)}}& =& {\displaystyle \sum _{m=0}^{r-1}\left(\genfrac{}{}{0pt}{}{n}{m}\right){\nu }_{r-m,s-m:n-m}^{(k_1,k_2)}+{(-1)}^{k_1}\sum _{m=r}^{s-1}\left(\genfrac{}{}{0pt}{}{n}{m}\right){\nu }_{m-r+1:m}^{\left(k_1\right)}{\nu }_{s-m:n-m}^{\left(k_2\right)}}\hfill \\ & & {\displaystyle +{(-1)}^{k_1+k_2}\sum _{m=s}^n\left(\genfrac{}{}{0pt}{}{n}{m}\right){\nu }_{m+1-s,m+1-r:m}^{(k_1,k_2)}}.\hfill \end{array}$$

(9)

The cdf of the random variable Z is denoted by $F_Z\left(z\right)$ and suppose $U,V\sim F_Z\left(z\right)$, then let ${\alpha }_{r,s:n}^{(k_1,k_2)}$ denote $E\left(U_{r:n}^{k_1}{V}_{s:n}^{k_2}\right) (r\le s)$. Obsiviously, we can write that ${\nu }_{r,s:n}^{(k_1,k_2)}={\alpha }_{r,s:n}^{(k_1/\beta ,k_2/\beta )}$.

First of all, we decide to derive ${\alpha }_{r,s:n}^{(k_1,k_2)}$. Taking the well-known expression of the product moment of order statistics and then using the binomial expansion, we obtain

$$\begin{array}{ccc}\hfill {\displaystyle {\alpha }_{r,r+1:s}^{(k_1,k_2)}}& =& {\displaystyle C_{r,r+1:s}\underset{0<x<y<+\infty }{\int \int }{x}^{k_1}{y}^{k_2}{\left[F\left(x\right)\right]}^{r-1}{[1-F\left(y\right)]}^{s-r-1}f\left(x\right)f\left(y\right)\mathsf{d}x\mathsf{d}y}\hfill \\ & =& {\displaystyle C_{r,r+1:s}\sum _{i=0}^{r-1}{(-1)}^i\left(\genfrac{}{}{0pt}{}{r-1}{i}\right){\int }_0^{+\infty }{x}^{k_1}{[1-F\left(x\right)]}^{i}f\left(x\right)\left({\int }_x^{+\infty }{y}^{k_2}{[1-F\left(y\right)]}^{s-r-1}f\left(y\right)\mathsf{d}y\right)\mathsf{d}x}\hfill \\ & =& {\displaystyle C_{r,r+1:s}\sum _{i=0}^{r-1}{(-1)}^i\left(\genfrac{}{}{0pt}{}{r-1}{i}\right)I(k_1,k_2,i,s-r-1)},\hfill \end{array}$$

(10)

where

$$C_{r,s:n}=\frac{s!}{(r-1)!(s-r-1)!}.$$

After that, we attempt to derive the double integral $I(k_1,k_2,i,s-r-1)$. By the multinomial theorem we have

$$\begin{array}{c}\hfill {\displaystyle {[1-F\left(y\right)]}^{s-r-1}}={\displaystyle \frac{1}{{\left[B\left(\alpha ,\frac{1}{\beta }\right)\right]}^{s-r-1}}\sum _{v=0}^{(s-r-1)(\frac{1}{\beta }-1)}{H}_v\left(\frac{1}{\beta },s-r-1\right)\times {(1+\beta y)}^{-\left[\alpha \right(s-r-1)+v]}}.\end{array}$$

We denote the inner integral ${\int }_x^{+\infty }{y}^{k_2}{[1-F\left(y\right)]}^{s-r-1}f\left(y\right)\mathsf{d}y$ of the above double integral $I(k_1,k_2,i,s-r-1)$ by $T_x(k_2,s-r-1)$ and it can be calculated as follows

$$\begin{array}{cc}\hfill T_x(k_2,s-r-1)& =\frac{{\left(\frac{1}{\beta }\right)}^{-\frac{1}{\beta }}}{{\left[B\left(\alpha ,\frac{1}{\beta }\right)\right]}^{s-r}}\sum _{v=0}^{(s-r-1)(\frac{1}{\beta }-1)}{H}_v\left(\frac{1}{\beta }-1,s-r-1\right)\times {\int }_x^{+\infty }{y}^{k_2+\frac{1}{\beta }-1}{(1+\beta y)}^{-[\alpha (s-r)+v+\frac{1}{\beta }]}\mathsf{d}y\hfill \\ \hfill & =\frac{{\left(\frac{1}{\beta }\right)}^{k_2}}{{\left[B\left(\alpha ,\frac{1}{\beta }\right)\right]}^{s-r}}\sum _{v=0}^{(s-r-1)(\frac{1}{\beta }-1)}{H}_v\left(\frac{1}{\beta },s-r-1\right)\times {\int }_{\frac{\beta x}{1+\beta x}}^1{w}^{k_2+\frac{1}{\beta }-1}{(1-w)}^{\alpha (s-r)+v-k_2-1}\mathsf{d}w.\hfill \end{array}$$

Notice that $(k_2+\frac{1}{\beta }-1)$ is still a positive integer, then expanding the integrand on the right hand side of the above expression in terms of $(1-w)$, we acquire

$$\begin{array}{ccc}\hfill {\displaystyle T_x(k_2,s-r-1)}& =& {\displaystyle \frac{{\left(\frac{1}{\beta }\right)}^{k_2}}{{\left[B\left(\alpha ,\frac{1}{\beta }\right)\right]}^{s-r}}\sum _{v=0}^{(s-r-1)(\frac{1}{\beta }-1)}{H}_v\left(\frac{1}{\beta },s-r-1\right)\sum _{m=0}^{k_2+\frac{1}{\beta }-1}{(-1)}^m}\hfill \\ & & {\displaystyle \times \left(\genfrac{}{}{0pt}{}{k_2+\frac{1}{\beta }-1}{m}\right)\frac{{(1+\beta x)}^{-[\alpha (s-r)+m+v-k_2]}}{\alpha (s-r)+m+v-k_2}}.\hfill \end{array}$$

(11)

Using the fact that

$$\begin{array}{c}\hfill {\displaystyle {[1-F\left(x\right)]}^i}={\displaystyle \frac{1}{{\left[B\left(\alpha ,\frac{1}{\beta }\right)\right]}^i}\sum _{u=0}^{i(\frac{1}{\beta }-1)}{H}_u\left(\frac{1}{\beta },i\right)\times {(1+\beta y)}^{-(\alpha i+u)}},\end{array}$$

(12)

by putting Equations (11) and (12) into $I(k_1,k_2,i,s-r-1)$, we obtain

$$\begin{array}{ccc}\hfill {\displaystyle I(k_1,k_2,i,s-r-1)}& =& {\displaystyle \frac{{\left(\frac{1}{\beta }\right)}^{k_1+k_2}}{{\left[B\left(\alpha ,\frac{1}{\beta }\right)\right]}^{s-r+i+1}}\sum _{u=0}^{(\frac{1}{\beta }-1)i}\sum _{v=0}^{(s-r-1)(\frac{1}{\beta }-1)}{H}_u\left(\frac{1}{\beta },i\right)H_v\left(\frac{1}{\beta },s-r-1\right)}\hfill \\ & & {\displaystyle \times \sum _{m=0}^{k_2+\frac{1}{\beta }-1}\left(\genfrac{}{}{0pt}{}{k_2+\frac{1}{\beta }-1}{m}\right)\frac{{(-1)}^m}{\alpha (s-r)+m+v-k_2}}\hfill \\ & & {\displaystyle \times B\left[k_1+\frac{1}{\beta },\alpha (s-r+i+1)+m+u+v-k_1-k_2\right]}.\hfill \end{array}$$

(13)

The above equation will exists for $\alpha >(k_1+k_2)/2$. Finally, using recurrence relations from Arnold et al. (2008) [19]:

For any arbitray distribution and $1\le p<q\le n$

$$\begin{array}{c}\hfill {\displaystyle {\alpha }_{p,q:n}^{(k_1,k_2)}}={\displaystyle \sum _{r=p}^{q-1}\sum _{s=n-q+r+1}^n{(-1)}^{s+n-p-q+1}\left(\genfrac{}{}{0pt}{}{r-1}{p-1}\right)\left(\genfrac{}{}{0pt}{}{s-r-1}{n-q}\right)\left(\genfrac{}{}{0pt}{}{n}{s}\right){\alpha }_{r,r+1:s}^{(k_1,k_2)}},\end{array}$$

(14)

we can derive a closed-form expression for every product moment and the corresponding variance–covariance matrix of order statistics from $F_Z\left(z\right)$.

When the sample size is n, the method of obtaining all the product moments and variance–covariance matrix for order statistics from $N\mathrm{Gt}(\mu ,\sigma ,\alpha ,\beta )$ can be briefly described as:

1.

Derive ${\nu }_{i,j:n}^{(k_1,k_2)}$ ($1\le i\le j\le n$) by using the above approach.

2.

For $1\le r\le \left[\frac{n}{2}\right],r\le s\le n-r+1$, obtain ${\mu }_{r,s:n}^{(k_1,k_2)}$ by using recurrence relations in Equation (9) and obtain the covariance ${\sigma }_{r,s:n}$ by ${\sigma }_{r,s:n}={\mu }_{r,s:n}-{\mu }_{r:n}{\mu }_{s:n}$.

3.

By using the symmetry of Gt distribution and the variance–covariance matrix, the remaining covariance of order statistics can be calculated by ${\sigma }_{r,s:n}={\sigma }_{n-s+1,n-r+1:n}$ and ${\sigma }_{r,s:n}={\sigma }_{s,r:n}$.

5. Parameter Estimation

In this section, we propose several parameter estimation methods for the new Gt distribution.

5.1. Modfied Method of Moments

In this subsection, a modified method of moments (MMOM) for the parameters of $N\mathrm{Gt}(0,\sigma ,\alpha ,\beta )$ is proposed.

Let $X_1,\dots ,X_n$ be a random sample from $N\mathrm{Gt}(\sigma ,\alpha ,\beta )$. Making transformation $Y={\left|X\right|}^2$, the k-th theoretical central moments of Y are given by

$$E\left(Y^k\right)={\left(\frac{2}{\beta }\right)}^{\frac{2k}{\beta }}\frac{\Gamma (\alpha -2k/\beta )\Gamma \left[\right(2k+1)/\beta ]}{\Gamma \left(\alpha \right)\Gamma (1/\beta )}{\sigma }^{2k},k=1,2,\cdots ,$$

The MMEs of the new distribution after transformation can be obtained by equating the first three theoretical central moments of Y with the corresponding central sample moments. Eliminating the scale parameter $\sigma$, we acquire

$$\begin{array}{ccc}\hfill {\displaystyle \frac{[\Gamma (\alpha \left)\right]\Gamma (\alpha -4/\beta )}{{[\Gamma (\alpha -2/\beta )]}^2}\times \frac{[\Gamma (1/\beta \left)\right]\Gamma (5/\beta )}{{[\Gamma (3/\beta )]}^2}}& =& {\displaystyle \frac{n{\sum }_{i=1}^n{Y}_i^2}{{\left({\sum }_{i=1}^n{Y}_i\right)}^2}},\hfill \end{array}$$

(15)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{{[\Gamma \left(\alpha \right)]}^2\Gamma (\alpha -6/\beta )}{{[\Gamma (\alpha -2/\beta )]}^3}\times \frac{{[\Gamma (1/\beta )]}^2\Gamma (7/\beta )}{{[\Gamma (3/\beta )]}^3}}& =& {\displaystyle \frac{n^2{\sum }_{i=1}^n{Y}_i^3}{{\left({\sum }_{i=1}^n{Y}_i\right)}^3}}.\hfill \end{array}$$

(16)

In order to establish an effective algorithm to obtain the numerical solution of the above equations, we propose the following theorem:

Theorem 5.

Let

$$G_k(\alpha ,\beta )=\frac{{[\Gamma \left(\alpha \right)]}^{k-1}\Gamma (\alpha -2k/\beta )}{{[\Gamma (\alpha -2/\beta )]}^k}\times \frac{{[\Gamma (1/\beta )]}^{k-1}\Gamma [(2k+1)/\beta ]}{{[\Gamma (3/\beta )]}^k}k>2,$$

then we assert that $G_k(\alpha ,\beta )$ is a monotonically-decreasing function of the parameters α and β. In addition, the function has the following limiting properties

$$Forfixed\alpha :\underset{\beta ⟶2k/{\alpha }^{+}}{\lim }{G}_k(\alpha ,\beta )=\infty ,\underset{\beta ⟶+\infty }{\lim }{G}_k(\alpha ,\beta )=\frac{3^{k/2}}{\sqrt{2k+1}}\times \frac{{\displaystyle \prod _{i=1}^k}\Gamma [i/(2k+1)]}{{[\Gamma (1/3)]}^k{[\Gamma (2/3)]}^k}.$$

$$Forfixed\beta :\underset{\alpha ⟶2k/{\beta }^{+}}{\lim }{G}_k(\alpha ,\beta )=\infty ,\underset{\alpha ⟶+\infty }{\lim }{G}_k(\alpha ,\beta )=\frac{{[\Gamma (1/\beta )]}^{k-1}\Gamma [(2k+1)/\beta ]}{{[\Gamma (3/\beta )]}^k}.$$

Proof. 

See Appendix A.3 (Part III).   □

With the above theorem, the steps of the algorithm are presented by the following forms

1.

Assign an initial values ${\alpha }_{\left(1\right)}$, and then substitute ${\alpha }_{\left(1\right)}$ into Equation (15). Since $G_k(\alpha ,\beta )$ is monotonically-decreasing with respect to parameter $\beta$, $uniroot$ in the R function can be used to derive ${\widehat{\beta }}_{\left(1\right)}$.

2.

Substitute ${\widehat{\beta }}_{\left(1\right)}$ into Equation (16), since $G_k(\alpha ,\beta )$ is monotonically-decreasing with respect to parameter $\alpha$, using $uniroot$ in R function then we obtain ${\widehat{\alpha }}_{\left(2\right)}$.

3.

If the norm of vector $({\widehat{\alpha }}_{(i+1)}-{\widehat{\alpha }}_{\left(i\right)},{\widehat{\beta }}_{(i+1)}-{\widehat{\beta }}_{\left(i\right)})$ exceed the prescribed value $\epsilon >0$, then 2–3 are repeated. Otherwise, the iteration is broken up and the values generated in the last iteration are determined to be estimates of the shape parameters $\alpha$ and $\beta$.

5.2. MLE Using the EM Algorithm

In this subsection, an EM-type algorithm to determine the MLEs for the parameters of $N\mathrm{Gt}(0,\sigma ,\alpha ,\beta )$ is established.

Through the method of distribution construction, the pdf of the new Gt distribution is treated as a marginal pdf of a bivariate distribution and this bivariate distribution function can be expressed as the product of two common distributions. After that, we treat one random variable in the bivariate distribution as an observed variable, while the other as an unobserved variable, so that the EM algorithm can be applied to obtain the MLE. At the same time, this idea is also widely used to derive the MLE of other distributions. See [15,20,21] for more details about recent works in this topic.

First proposed by Dempster (1977) [22], the EM algorithm is a powerful tool to deal with estimation problem under missing data situation. In this subsection, we are going to use the EM algorithm to acquire the MLEs of $N\mathrm{Gt}(\sigma ,\alpha ,\beta )$ based on Proposition 2. Let $Y={(y_1,\dots ,y_n)}^T$ denote the observed data, $U={(u_1,\dots ,u_n)}^T$ on behalf of the unobserved data. Combining Y and U together we derive the complete data denoted by $Z=(Y,U)$. Besides that, let $\theta ={(\sigma ,\alpha ,\beta )}^T$ represent the vector of parameters. Invoking Proposition 2, the complete log-likelihood function can be expressed as follows

$$\begin{array}{ccc}\hfill {\displaystyle lnL\left(\theta \right|Z)}& =& {\displaystyle n\left[\left(1+\frac{1}{\beta }\right)ln\beta -\left(1+\frac{2}{\beta }+\alpha \right)ln2-ln\Gamma \left(\frac{1}{\beta }\right)-ln\Gamma \left(\alpha \right)-ln\sigma \right]}\hfill \\ & & {\displaystyle +\left(\alpha +\frac{1}{\beta }-1\right)\sum _{i=1}^{n}ln{u}_i-\frac{\beta }{4{\sigma }^{\beta }}\sum _{i=1}^n{\left|y_i\right|}^{\beta }{u}_i-\frac{1}{2}\sum _{i=1}^n{u}_i}.\hfill \end{array}$$

(17)

It is well-known that the EM framework is an iterative algorithm consisting of two steps; namely, the expectation step (E-step) and the maximization step (M-step).

E-step: Given the observed data set $Y=(y_1,\dots ,y_n)$ and parameter estimation values ${\theta }_{\left(h\right)}$ for the h-th iteration, the conditional expectation of the complete log-likelihood function can be calculated as follows

$$\begin{array}{ccc}\hfill {\displaystyle Q\left(\theta \right|{\theta }_{\left(h\right)})}={\displaystyle E[lnL\left(\theta \right|Z)|Y,{\theta }_{\left(h\right)}]}& =& {\displaystyle n\left[\left(1+\frac{1}{\beta }\right)ln\beta -\left(1+\frac{2}{\beta }+\alpha \right)ln2-ln\Gamma \left(\frac{1}{\beta }\right)-ln\Gamma \left(\alpha \right)-ln\sigma \right]}\hfill \\ & & {\displaystyle +\left(\alpha +\frac{1}{\beta }-1\right)\sum _{i=1}^{n}E[ln{u}_i|y_i,{\theta }_{\left(h\right)}]-\frac{\beta }{4{\sigma }^{\beta }}\sum _{i=1}^n{\left|y_i\right|}^{\beta }E\left[u_i\right|y_i,{\theta }_{\left(h\right)}]}\hfill \\ & & {\displaystyle -\frac{1}{2}\sum _{i=1}^{n}E\left[u_i\right|y_i,{\theta }_{\left(h\right)}]},\hfill \end{array}$$

where

$$\begin{array}{ccc}\hfill {\displaystyle E\left[u_i\right|y_i,{\theta }_{\left(h\right)}]}& =& {\displaystyle \frac{4{\sigma }_{\left(h\right)}^{{\beta }_{\left(h\right)}}({\alpha }_{\left(h\right)}+1/{\beta }_{\left(h\right)})}{2{\sigma }_{\left(h\right)}^{{\beta }_{\left(h\right)}}+{\beta }_{\left(h\right)}{\left|y_i\right|}^{{\beta }_{\left(h\right)}}}},\hfill \\ \hfill {\displaystyle E[ln{u}_i|y_i,{\theta }_{\left(h\right)}]}& =& {\displaystyle \phi \left({\alpha }_{\left(h\right)}+\frac{1}{{\beta }_{\left(h\right)}}\right)-ln\left(\frac{1}{2}+\frac{{\beta }_{\left(h\right)}{\left|y_i\right|}^{{\beta }_{\left(h\right)}}}{4{\sigma }_{\left(h\right)}^{{\beta }_{\left(h\right)}}}\right)}.\hfill \end{array}$$

In the EM algorithm, the M-step needs to maximize the conditional expectation obtained by E-step and the suggested framework can be briefly described as follows

1.

Update ${\alpha }_{(h+1)}$ through the following nonlinear equation and the corresponding numerical procedures is recommended for $uniroot$ in R program.

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial Q\left(\theta \right|{\theta }_{\left(h\right)})}{\partial \alpha }}={\displaystyle n[\phi \left(\alpha \right)+ln2]-\sum _{i=1}^{n}E[ln{u}_i|y_i,{\theta }_{\left(h\right)}]}.\end{array}$$

2.

Update ${\sigma }_{(h+1)}$ through the following explicit expression

$$\begin{array}{c}\hfill {\displaystyle {\sigma }_{(h+1)}}={\displaystyle {\left(\frac{{\beta }_{\left(h\right)}^2}{4n}\right)}^{\frac{1}{{\beta }_{\left(h\right)}}}{\left(\sum _{i=1}^n{\left|y_i\right|}^{{\beta }_{\left(h\right)}}E\left[u_i\right|y_i,{\sigma }_{\left(h\right)},{\alpha }_{(h+1)},{\beta }_{\left(h\right)}]\right)}^{\frac{1}{{\beta }_{\left(h\right)}}}}.\end{array}$$

3.

Update ${\beta }_{(h+1)}$ by maximizing the observed data log-likelihood function

$$\begin{array}{c}\hfill {\displaystyle {\beta }_{(h+1)}}={\displaystyle \underset{\beta }{\arg \max }\sum _{i=1}^{n}Q\left(\theta \right|{\sigma }_{(h+1)},{\alpha }_{(h+1)},\beta ).}\end{array}$$

For the initial values ${\theta }_0=({\sigma }_0,{\alpha }_0,{\beta }_0)$, we recommend conducting numerical procedures such as $optim$ in R program for maximizing the corresponding likelihood function or using the estimation values from the method of moments.

5.3. Explicit Expressions for Asymptotic Variances and Covariances of the Estimates under the EM Algorithm

In this subsection, we obtain closed-form expression for the asymptotic variances and covariances of the estimates under the EM algorithm. Louis (1982) [23] coined the missing information principle to acquire an information matrix in a situation containing latent variables. This principle is given by

$$\mathrm{Observed}\mathrm{information}=\mathrm{Complete}\mathrm{information}-\mathrm{Missing}\mathrm{information}.$$

In the following paper, we will denote the complete, observed and missing information by $I^Z\left(\theta \right)$, $I^Y\left(\theta \right)$ and $I^{\left(U\right|Y)}\left(\theta \right)$, respectively. Then the observed information $I^Y\left(\theta \right)$ is calculated by:

$$I^Y\left(\theta \right)=I^Z\left(\theta \right)-I^{\left(U\right|Y)}\left(\theta \right).$$

5.3.1. Complete Information Matrix $I^Z\left(\theta \right)$

Taking the second-order partial derivatives of the complete log-likelihood function in Equation (17) and applying the definition of the fisher information, then after some calculation, we obtain

$$\begin{array}{c}\hfill {\displaystyle I_{11}^Z}={\displaystyle \frac{n\beta }{{\sigma }^2}},{\displaystyle I_{22}^Z=n{\phi }^{\prime }\left(\alpha \right)},{\displaystyle I_{12}^Z=I_{21}^Z=I_{23}^Z=I_{32}^Z=0},\end{array}$$

$$\begin{array}{c}\hfill {\displaystyle I_{13}^Z=I_{31}^Z=-\frac{n}{\sigma \beta }+\frac{{\beta }^2}{4\sigma }n{T}_1(\alpha ,\beta )+\frac{n{\beta }^{\frac{1}{\beta }}}{2^{\frac{1}{\beta }}B(\alpha ,1/\beta )}},\end{array}$$

$$\begin{array}{ccc}\hfill {\displaystyle I_{33}^Z}& =& {\displaystyle \frac{2n}{{\beta }^3}\left[-ln\beta +(1+\beta )+2ln2+\phi \left(\frac{1}{\beta }\right)\right]+\frac{n}{{\beta }^2}\left[\frac{1}{\beta }-1+\frac{\phi \left(\frac{1}{\beta }\right)}{{\beta }^2}\right]-\frac{2}{{\beta }^3}n\left[\phi \left(\alpha \right)+\frac{ln2}{2^{\alpha }}\right]}\hfill \\ & & {\displaystyle +\frac{n}{2}{T}_1(\alpha ,\beta )+\frac{\beta n}{4}{T}_2(\alpha ,\beta )},\hfill \end{array}$$

where

$$\begin{array}{c}\hfill {\displaystyle T_1(\alpha ,\beta )}={\displaystyle \frac{4}{{\beta }^3}\left[\phi (1+\frac{1}{\beta })-ln\frac{\beta }{2}-\phi \left(\alpha \right)\right]},\end{array}$$

$$\begin{array}{ccc}\hfill {\displaystyle T_2(\alpha ,\beta )}& =& {\displaystyle \frac{2}{{\beta }^3}\left\{{\left[\phi \left(1+\frac{1}{\beta }\right)\right]}^2+\frac{1}{{\beta }^2}{\phi }^{\prime }\left(1+\frac{1}{\beta }\right)\right\}+\frac{8}{{\beta }^3}\left[\phi \left(1+\frac{1}{\beta }\right)ln\frac{\beta }{2}+\phi \left(1+\frac{1}{\beta }\right)\phi \left(\alpha \right)\right]\}}\hfill \\ & & {\displaystyle +\frac{4}{{\beta }^3}\left\{{\left(\frac{\beta }{2}\right)}^2+{\left[\phi \left(\alpha \right)\right]}^2+{\phi }^{\prime }\left(\alpha \right)+2\phi \left(\alpha \right)ln\frac{\beta }{2}\right\}}.\hfill \end{array}$$

The second-order partial derivatives of the complete likelihood function and the calculation method of the complete information are shown in the Appendix A.3 (Part III).

5.3.2. Missing Information Matrix $I^{\left(U\right|Y)}\left(\theta \right)$

From Proposition 2, the logarithm of the joint density function of latent variable U given Y can be written as

$$ln{f}_{U|Y}\left(u\right|y)=\sum _{i=1}^n\left(\alpha +\frac{1}{\beta }\right)ln\left(\frac{1}{2}+\frac{\beta |y_i|}{4{\sigma }^{\beta }}\right)-n\Gamma (\alpha +1/\beta )+\left(\alpha +\frac{1}{\beta }-1\right)\sum _{i=1}^{n}ln{U}_i-\sum _{i=1}^n\left(\frac{1}{2}+\frac{\beta |y_i|}{4{\sigma }^{\beta }}\right)U_i.$$

After taking the second-order partial derivatives with respect to $\sigma$, $\alpha$ and $\beta$, and utilizing the results in Section 5.2, the elements of the missing information matrix $I^{\left(U\right|Y)}\left(\theta \right)$ are presented as follows

$$\begin{array}{ccc}\hfill {\displaystyle I_{11}^{\left(U\right|Y)}}& =& {\displaystyle -\beta (\alpha +1)\sum _{i=1}^n\left\{\frac{|y_i\left|\right[2(\beta +1){\sigma }^{\beta }+\beta |y_i{|}^2]}{2{\sigma }^{\beta +1}+\beta \sigma {\left|y_i\right|}^2}-\frac{|y_i|(\beta +1)}{{\sigma }^2(2{\sigma }^{\beta }+\beta |y_i{|}^{\beta })}\right\}},\hfill \\ \hfill {\displaystyle I_{12}^{\left(U\right|Y)}}& =& {\displaystyle I_{21}^{\left(U\right|Y)}}={\displaystyle -\sum _{i=1}^n\frac{{\beta }^2\left|y_i\right|}{2{\sigma }^{\beta +1}+\beta \sigma \left|y_i\right|}},\hfill \\ \hfill {\displaystyle I_{13}^{\left(U\right|Y)}}& =& {\displaystyle I_{31}^{\left(U\right|Y)}}={\displaystyle -\sum _{i=1}^n\frac{|y_i|}{2{\sigma }^{\beta +1}+\beta \sigma \left|y_i\right|}+(\alpha +\frac{1}{\beta })\sum _{i=1}^n\frac{2\beta |y_i\left|\right(2{\sigma }^{\beta +1}+\beta \sigma |y_i\left|\right)-{\beta }^2|y_i\left|\right(2{\sigma }^{\beta +1}ln\sigma +\sigma |y_i\left|\right)}{(2{\sigma }^{\beta +1}+\beta \sigma |y_i{\left|\right)}^2}}\hfill \\ & & {\displaystyle +\frac{2\beta |y_i|-{\beta }^2\sigma \left|y_i\right|ln\sigma }{4{\sigma }^{\beta +1}}},\hfill \\ \hfill {\displaystyle I_{22}^{\left(U\right|Y)}}& =& {\displaystyle n\frac{{\partial }^2\Gamma (\alpha +1/\beta )}{\partial {\alpha }^2}},\hfill \\ \hfill {\displaystyle I_{23}^{\left(U\right|Y)}}& =& {\displaystyle I_{32}^{\left(U\right|Y)}}={\displaystyle -\frac{|y_i|(1-\beta ln\sigma )}{2{\sigma }^{\beta }+\beta \left|y_i\right|}+\frac{n}{{\beta }^4}{\phi }^{\prime }(\alpha +1/\beta )},\hfill \\ \hfill {\displaystyle I_{33}^{\left(U\right|Y)}}& =& {\displaystyle -\frac{2}{{\beta }^2}\sum _{i=1}^{n}ln\left(\frac{1}{2}+\frac{\beta |y_i|}{4{\sigma }^{\beta }}\right)+\sum _{i=1}^n\frac{|y_i|}{2{\sigma }^{\beta +1}+\beta \sigma \left|y_i\right|}+\frac{1}{{\beta }^2}\sum _{i=1}^n\frac{|y_i|(1-\beta ln\sigma )}{2{\sigma }^{\beta }+\beta \left|y_i\right|}+\left(\frac{2n}{{\beta }^3}+\frac{n}{{\beta }^4}\right)\frac{\partial \Gamma (\alpha +1/\beta )}{\partial \beta }}\hfill \\ & & {\displaystyle +\left(\alpha +\frac{1}{\beta }\right)\sum _{i=1}^n\left[\frac{|y_i\left|\right[(2{\sigma }^{\beta }+\beta |y_i\left|\right)ln\sigma +(2{\sigma }^{\beta }ln\sigma +|y_i\left|\right)|y_i\left|(1-\beta ln\sigma )\right]}{(2{\sigma }^{\beta }+\beta |y_i{\left|\right)}^2}+\frac{|y_i|(2-\beta ln\sigma )}{2{\sigma }^{\beta }+\beta {\left|y_i\right|}^{\beta }}\right]}\hfill \\ & & {\displaystyle -\frac{2n}{{\beta }^3}\left[\phi \left(\alpha +\frac{1}{\beta }\right)-ln\left(\frac{1}{2}+\frac{\beta |y_i|}{4{\sigma }^{\beta }}\right)\right]}.\hfill \end{array}$$

5.4. MLE via a New Iterative Algorithm

In this subsection, we propose a new iterative algorithm to obtain the MLEs for the parameters of $N\mathrm{Gt}(0,\sigma ,\alpha ,\beta )$ because of the difficulty of sloving the estimation equations directly.

Let $X_1,\cdots ,X_n$ be a random sample from $N\mathrm{Gt}(0,\sigma ,\alpha ,\beta )$, and the log-likelihood function is given by

$$lnL(X;\sigma ,\alpha ,\beta )=-n(1+1/\beta )ln(2/\beta )-nlnB(\alpha ,1/\beta )-nln\sigma -\left(\alpha +\frac{1}{\beta }\right)\sum _{i=1}^{n}ln\left[1+\frac{\beta |x_i{|}^{\beta }}{2{\sigma }^{\beta }}\right].$$

With partial derivatives respect to $(\sigma ,\alpha ,\sigma )$, we obtain

$$\left\{\begin{array}{ccc}\hfill {\displaystyle \frac{\partial lnL(X;\sigma ,\alpha ,\beta )}{\partial \sigma }}& =& {\displaystyle n-(\alpha \beta +1)\sum _{i=1}^n\frac{\beta |x_i{|}^{\beta }}{2{\sigma }^{\beta }+\beta {\left|x_i\right|}^{\beta }}},\hfill \\ \hfill {\displaystyle \frac{\partial lnL(X;\sigma ,\alpha ,\beta )}{\partial \alpha }}& =& {\displaystyle n\left[\phi (\alpha +1/\beta )-\phi \left(\alpha \right)\right]-\sum _{i=1}^{n}ln\left[1+\frac{\beta |x_i{|}^{\beta }}{2{\sigma }^{\beta }}\right]},\hfill \\ \hfill {\displaystyle \frac{\partial lnL(X;\sigma ,\alpha ,\beta )}{\partial \beta }}& =& {\displaystyle nln(2/\beta )+n\beta (1+1/\beta )-n\left[\phi (\alpha +1/\beta )-\phi (1/\beta )\right]+\sum _{i=1}^{n}ln\left[1+\frac{\beta |x_i{|}^{\beta }}{2{\sigma }^{\beta }}\right]}\hfill \\ \hfill & & {\displaystyle -\beta (\alpha \beta +1)\sum _{i=1}^n\frac{\left(|x_i{|}^{\beta }+\beta |x_i{|}^{\beta }ln\left(\right|x_i|/\sigma )\right)}{2{\sigma }^{\beta }+\beta {\left|x_i\right|}^{\beta }}}.\hfill \end{array}\right.\left(18\right)$$

(18)

Instead of sloving all three preceding equations simultaneously, we propose a new algorithm because of the complexity of the third equation of Equation (18). In order to establish the new approach, the following lemma is necessary

Lemma 2.

For any fixed β, the following equations

$$\left\{\begin{array}{c}{\displaystyle n-(\alpha \beta +1)\sum _{i=1}^n\frac{\beta |x_i{|}^{\beta }}{2{\sigma }^{\beta }+\beta {\left|x_i\right|}^{\beta }}}=0,\hfill \\ {\displaystyle n\left[\phi (\alpha +1/\beta )-\phi \left(\alpha \right)\right]-\sum _{i=1}^{n}ln\left[1+\frac{\beta |x_i{|}^{\beta }}{2{\sigma }^{\beta }}\right]}=0,\hfill \end{array}\right.\left(19\right)$$

(19)

always possess a convergent solution of parameter $(\sigma ,\alpha )$ and is independent of the initial values $({\sigma }_0,{\alpha }_0)$.

Proof. 

See Appendix A.3 (Part III).   □

Moreover, the algorithm for obtaining the solution of the equation in Lemma 2 is shown below

1.

Given arbitrary initial values $({\sigma }_{\left(0\right)},{\alpha }_{\left(0\right)})$ and ${\sigma }_{\left(0\right)}>0,{\alpha }_{\left(0\right)}>0$.

2.

Substitute ${\alpha }_{\left(0\right)}$ into the first equation of Equation (19), and we obtain ${\sigma }_{\left(1\right)}$.

3.

Substitute ${\sigma }_{\left(1\right)}$ into the second equation of Equation (19), and we obtain ${\alpha }_{\left(1\right)}$.

4.

If the norm of vector $({\sigma }_{(h+1)}-{\sigma }_{\left(h\right)},{\alpha }_{(h+1)}-{\alpha }_{\left(h\right)})$ exceed the prescribed value $\epsilon >0$, then 2–3 are repeated. Otherwise, the iteration is broken up and the values generated in the last iteration are determined to be solution of this equation.

In what follows, the new algorithm to acquire the solution of Equation (18) can be briefly described as

1.

Give an appropriate interval for the shape parameter $\beta$, and then generate a set of values of this parameter at an appropriate distance within the interval and represent the j-th element of this vector by ${\beta }^{\left(j\right)}$.

2.

Assigning ${\beta }^{\left(j\right)}$ into the first two equations of Equation (18)and applying Lemma 2, we can obtain their convergence solution denoted by $({\sigma }^{\left(j\right)},{\alpha }^{\left(j\right)})$.

3.

Repeat this procedure and observe the sign-changing in the left-hand side of the third equation of Equation (18). Then, by interpolating on $\beta$ values, we derive the final estimate ${\widehat{\beta }}_{MLE}$, which makes the third equation equals to 0 or sufficiently close to 0. Finally, by substituting ${\widehat{\beta }}_{MLE}$ into the first two equations of Equation (18), we derive the final estimates denoted by $({\widehat{\sigma }}_{MLE},{\widehat{\alpha }}_{MLE})$.

It is noteworthy that this algorithm is applicable when the sample kurtosis is greater than 2.7. Comparing with other numerical methods, our new algorithm requires less precision of initial values $({\sigma }_0,{\alpha }_0)$ and can obtain convergent solutions.

5.5. An Improved PWM Estimation

In this subsection, we propose an improved PWM (IPWM) estimation based on order statistics transformation for $N\mathrm{Gt}(0,\sigma ,\alpha ,\beta )$ and develop the corresponding algorithm to implement this novel method.

PWM estimation is a very effective and commonly used method for parameter estimation of distributions with heavy tails. Among many parameter estimation methods proposed so far, PWM and MLE are the two most commonly-used methods that have similar priority.

Let $X_1,\cdots ,X_n$ be a random sample from $N\mathrm{Gt}(0,\sigma ,\alpha ,\beta )$ and $X_{1:n}\le \cdots \le X_{n:n}$ behalf of the corresponding order statistics. The traditional theoretical PWM was introduced by Greenwood et al. [24] as an alternative to the method of moments given in the following equation

$$\begin{array}{ccc}\hfill {\displaystyle M_{p,r,s}\left(x\right)}& =& {\displaystyle E\left\{X^p{\left[G\left(x\right)\right]}^r{[1-G\left(x\right)]}^s\right\},r=0,1,\cdots ,s=0,1,\dots .}\hfill \end{array}$$

Taking $p=1,s=0$, the moments become $M_{1,r,0}\left(x\right)$ denoted by ${\beta }_r$. Wang (1990) [25] obtained an unbiased estimator of ${\beta }_r$ through linear combination of order statistics as follows

$$b_r=\frac{1}{n}{\left(\genfrac{}{}{0pt}{}{n-1}{r}\right)}^{-1}\sum _{i=1}^n\left(\genfrac{}{}{0pt}{}{i-1}{r}\right)x_{i:n},r=0,1,\cdots .$$

However, applying the traditional PWM directly can cause many problems due to the complexity of the distribution function. For example, intricacy expressions of the theoretical PWM will make the estimation equation tedious and difficult to be numerically optimized, thus impeding the use of this estimation method. As a result, we decide to make isotone transformation $Y=\sigma X$. After this monotone transformation, the estimation equation becomes simple and easy to solve by a numerical approach. The improved theoretical PWM (denoted by $M_{1,r,0}^Y$) under this transformation is given by

$$\begin{array}{c}\hfill {\displaystyle {\beta }_r^Y=M_{1,r,0}^Y}={\displaystyle E\left\{Y{\left[G\left(y\right)\right]}^r\right\},r=0,1,\cdots .}\end{array}$$

(20)

The corresponding improved sample PWM (denoted by $b_r^Y$) is defined as follows

$$\begin{array}{c}\hfill {\displaystyle b_r^Y}={\displaystyle \frac{1}{n}{\left(\genfrac{}{}{0pt}{}{n-1}{r}\right)}^{-1}\sum _{i=1}^n\left(\genfrac{}{}{0pt}{}{i-1}{r}\right)y_{i:n},r=0,1,\cdots ,}\end{array}$$

(21)

where $y_{i:n}=\tilde{\sigma }{x}_{i:n}$ and $\tilde{\sigma }$ is an initial guess of the scale parameter $\sigma$.

From Equation (20), we can derive the 2–5th improved theoretical PWM as follows

$$\begin{array}{ccc}\hfill {\displaystyle {\beta }_1^Y}& =& {\displaystyle \frac{{\left(\frac{2}{\beta }\right)}^{\frac{1}{\beta }}}{2B(\alpha ,\frac{1}{\beta })}{\int }_0^1{z}^{\frac{2}{\beta }-1}{(1-z)}^{\alpha -\frac{1}{\beta }-1}{I}_z\left(\frac{1}{\beta },\alpha \right)\mathsf{d}z},\hfill \\ \hfill {\displaystyle {\beta }_2^Y}& =& {\displaystyle {\beta }_1^Y},\hfill \\ \hfill {\displaystyle {\beta }_3^Y}& =& {\displaystyle \frac{2^{1/\beta -3}{\left(\frac{1}{\beta }\right)}^{\frac{1}{\beta }}}{B(\alpha ,\frac{1}{\beta })}\left\{{\int }_0^1{z}^{\frac{2}{\beta }-1}{(1-z)}^{\alpha -\frac{1}{\beta }-1}{I}_z\left(\frac{1}{\beta },\alpha \right)\mathsf{d}z+3{\int }_0^1{z}^{\frac{2}{\beta }-1}{(1-z)}^{\alpha -\frac{1}{\beta }-1}{\left[I_z\left(\frac{1}{\beta },\alpha \right)\right]}^3\mathsf{d}z\right\}},\hfill \end{array}$$

$$\begin{array}{ccc}\hfill {\displaystyle {\beta }_4^Y}& =& {\displaystyle \frac{2^{1/\beta -2}{\left(\frac{1}{\beta }\right)}^{\frac{1}{\beta }}}{B(\alpha ,\frac{1}{\beta })}\left\{{\int }_0^1{z}^{\frac{2}{\beta }-1}{(1-z)}^{\alpha -\frac{1}{\beta }-1}{I}_z\left(\frac{1}{\beta },\alpha \right)\mathsf{d}z+{\int }_0^1{z}^{\frac{2}{\beta }-1}{(1-z)}^{\alpha -\frac{1}{\beta }-1}{\left[I_z\left(\frac{1}{\beta },\alpha \right)\right]}^3\mathsf{d}z\right\}}.\hfill \end{array}$$

According to the traditional idea of the method of moments, let the improved theoretical PWM ${\beta }_3^Y,{\beta }_4^Y$ equal to the corresponding improved sample PWM. Then we can obtain an equation set whose solution is the estimation of the unknown parameters. Due to the complexity of these equations, we transform the search for the solution of this equation set into an optimization problem and acquire the simultaneous estimation of the shape parameters denoted by ${\widehat{\alpha }}_{PWM},{\widehat{\beta }}_{PWM}$ through minimizing the following expression to avoid failures during optimization

$$\begin{array}{c}\hfill ({\widehat{\alpha }}_{PWM},{\widehat{\beta }}_{PWM})={\displaystyle \underset{(\alpha ,\beta )}{\arg \min }{\left({\beta }_3^Y+{\beta }_4^Y-b_3^Y-b_4^Y\right)}^2}.\end{array}$$

(22)

Numerical techniques are recommended, such as the $optim$ function in R program and the method will be chosen for $L-BFGS-B$.

Finally, the estimation of the scale parameter ${\widehat{\sigma }}_{PWM}$ can be obtained through the first original sample PWM of the random variable X with analytical expression given by the following equation

$$\begin{array}{c}\hfill {\displaystyle {\widehat{\sigma }}_{PWM}}={\displaystyle \frac{2b_{1}B({\widehat{\alpha }}_{PWM},1/{\widehat{\beta }}_{PWM})}{{\left(\frac{2}{{\widehat{\beta }}_{PWM}}\right)}^{1/{\widehat{\beta }}_{PWM}}{\int }_0^1{u}^{2/{\widehat{\beta }}_{PWM}-1}{(1-u)}^{{\widehat{\alpha }}_{PWM}-{\widehat{\beta }}_{PWM}-1}{I}_u(1/{\widehat{\beta }}_{PWM},{\widehat{\alpha }}_{PWM})\mathsf{d}u}}.\end{array}$$

(23)

In the study, we found that the estimation results are sensitive to the choice of the initial guess $\tilde{\sigma }$. Therefore, we introduce the following two methods to discuss the choice of this parameter.

Method 1 (A two step estimation method)

1.

Take a consistent estimate of parameter $\sigma$ as the initial guess of it, such as MLE.

2.

Substitute this estimate value $\widehat{\sigma }$ into Equation (22) then we derive the improved PWM-estimation of the two shape parameters. The estimation for the scale paramter $\sigma$ can be obtained through Equation (23).

Through this method, considering that it is relatively fixed to select the value of parameter $\tilde{\sigma }$, and the final estimations are largely dependent on this parameter, sometimes it is inevitable that the fixed parameter $\tilde{\sigma }$ does not apply to all samples. As a result, we establish the following algorithm, which can flexibly select the parameter $\tilde{\sigma }$ and the suggested framework can be briefly described as:

Method 2 (A selection algorithm)

1.

Give an appropriate interval for the parameter $\tilde{\sigma }$, generate a set of $\tilde{\sigma }$ at appropriate distance within the interval and represent the h-th element of the vector by ${\tilde{\sigma }}_{\left(h\right)}$.

2.

Substitute each ${\tilde{\sigma }}_{\left(h\right)}$ into Equation (22) to acquire simultaneously estimation for two shape parameters denoted by ${\widehat{\alpha }}_{PWM}^{\left(h\right)}$ and ${\widehat{\beta }}_{PWM}^{\left(h\right)}$. After that, cauculate the corresponding original 6th theroetical PWM-kuritosis through the following equation

$$T_k({\widehat{\alpha }}_{PWM}^{\left(h\right)},{\widehat{\beta }}_{PWM}^{\left(h\right)})={\int }_0^1{u}^{\frac{2}{{\widehat{\beta }}_{PWM}^{\left(h\right)}}-1}{(1-u)}^{{\widehat{\alpha }}_{PWM}^{\left(h\right)}-\frac{1}{{\widehat{\beta }}_{PWM}^{\left(h\right)}}-1}{\left[I_u\left(\frac{1}{{\widehat{\beta }}_{PWM}^{\left(h\right)}},{\widehat{\alpha }}_{PWM}^{\left(h\right)}\right)\right]}^k\mathsf{d}u,$$

$${\tau }_6^{\left(h\right)PWM}=\frac{T_1({\widehat{\alpha }}_{PWM}^{\left(h\right)},{\widehat{\beta }}_{PWM}^{\left(h\right)})+2T_3({\widehat{\alpha }}_{PWM}^{\left(h\right)},{\widehat{\beta }}_{PWM}^{\left(h\right)})+T_5({\widehat{\alpha }}_{PWM}^{\left(h\right)},{\widehat{\beta }}_{PWM}^{\left(h\right)})}{T_1({\widehat{\alpha }}_{PWM}^{\left(h\right)},{\widehat{\beta }}_{PWM}^{\left(h\right)})}.$$

3.

Repeate 2 and cauculate $q^{PWM}{\left(\tilde{\sigma }\right)}_{\left(h\right)}=|{\tau }_6^{\left(h\right)PWM}-{\widehat{M}}_{1,0,5}^X/{\widehat{M}}_{1,0,1}^X|$ (${\widehat{M}}_{1,0,5}^X$ and ${\widehat{M}}_{1,0,1}^X$ denote the original 6th and 2th sample PWM of the random variable X, respectively). Then, interpolate on $\tilde{\sigma }$ values and derive $\breve{\sigma }$ that minimize $q^{PWM}\left(\tilde{\sigma }\right)$.

4.

Substituting $\breve{\sigma }$ into Equation (22) we derive the final estimation for two shape parameters and final estimation for $\sigma$ can be obtained form Equation (23).

5.6. Parameter Estimation of $N\mathrm{Gt}(\mu ,\sigma ,\alpha ,\beta )$ $(\mu \ne 0)$

In this subsection, we propose a profile maximum likelihood approach (PLA) using the EM algorithm to cope with parameter estimation for the new Gt distribution with four parameters.

In fitting problems with the distribution of data with better symmetry, three factors are usually comprehensively scrutinized, namely location information, scale parameters, and shape parameters. Therefore, with the addition of the location parameter, we focus on parameter estimation problem of $N\mathrm{Gt}(\mu ,\sigma ,\alpha ,\beta )$.

Balakrishnan (2019) [20] applied this approach to determine the MLEs for the parameters of Student’s t Birnbaum-Saunders distribution. As Balakrishnan points out, this method is an example of an exhaustive search. Due to the addition of location parameter $\mu$, the estimation equation becomes very complicated. To overcome this difficulty, $\mu$ can be assumed as a natural number through profile likelihood approach, reducing the dimension of the problem with four parameters to that with three. Then the EM algorithm for the three-parameter case described in the preceding chapter comes in handy here.

Let $Y_1,\dots ,Y_n$ be a random sample from $N\mathrm{Gt}(\mu ,\sigma ,\alpha ,\beta )$. In this case, the log-likelihood function is given by

$$lnL\left(y\right|\mu ,\sigma ,\alpha ,\beta )=-n(1+1/\beta )ln(2/\beta )-nlnB(\alpha ,1/\beta )-nln\sigma -\left(\alpha +\frac{1}{\beta }\right)\sum _{i=1}^{n}ln\left[1+\frac{\beta |y_i{-\mu |}^{\beta }}{2{\sigma }^{\beta }}\right].$$

Since $\frac{\beta |y_i{-\mu |}^{\beta }}{2{\sigma }^{\beta }+\beta {|y_i-\mu |}^{\beta }}<ln(1+y_i)<\frac{\beta |y_i{-\mu |}^{\beta }}{2{\sigma }^{\beta }}$, we can assert that the likelihood function $L\left(\mu \right)$ is bounded when the parameter $(\sigma ,\alpha ,\beta )$ are fixed.

A suggested framework for the PLA using the EM algorithm is briefly described as follows:

1.

Give an appropriate interval for the location parameter $\mu$. Then generate a set of $\mu$ at appropriate distance within the interval and represent the j-th element of the vector by ${\mu }^{\left(j\right)}$.

2.

Substitute each ${\mu }^{\left(j\right)}$ in the log-likelihood function and implement the EM algorithm to derive the MLEs of $(\sigma ,\alpha ,\beta )$.

3.

Repeating this procedure, and then interpolating on $\mu$ values, we obtain a curve of values of the log-likelihood function and choose the value of $\mu$ that maximizes the curve as ${\widehat{\mu }}_{MLE}$.

4.

Substitute ${\widehat{\mu }}_{MLE}$ in the log-likelihood function and implement the EM algorithm to determine the final MLEs denoted by $({\widehat{\sigma }}_{MLE},{\widehat{\alpha }}_{MLE},{\widehat{\beta }}_{MLE})$.

The implementation for the EM algorithm is similar to Section 5.2, just replace all the $|y_i|$ in the partial derivative of $\sigma ,\alpha ,\beta$ and the integral with $|y_i-{\widehat{\mu }}_{MLE}|$. The trace of the log-likelihood function is depicted in Figure 5. Through simulation research, we find that this algorithm performs well for $\alpha >\beta$ and is more accurate and stable than using the sample mean as the estimation of the location parameter $\mu$ directly under a small sample size.

6. Simulation Study

In this section, we conduct a simulation study to illustrate the performance of the proposed estimation methods and compare these simulation results under different parameter values and different sample sizes.

According to the actual demand, the sample size will be chosen for 25, 50, 100 and 200. We perform Monte Carlo simulations and evaluate the performance of the considered methods through relative bias (RB) and relative mean-square errors (RMSE). These two indicators will be computed with 5000 replications at the same time and they are defined as follows

$$RB\left(\widehat{\theta }\right)=\frac{1}{N}\sum _{i=1}^N\frac{{\widehat{\theta }}_i}{{\theta }_i},$$

$$RMSE\left(\widehat{\theta }\right)=\frac{1}{N}\sum _{i=1}^N\left[\frac{{\widehat{\theta }}_i}{{\theta }_i}-RB\left(\widehat{\theta }\right)\right].$$

This section will be detailed in three parts. First of all, we conduct a comparative study for three different estimation methods (MMOM, MLE using the EM algorithm and IPWM) of $N\mathrm{Gt}(0,\sigma ,\alpha ,\beta )$. We mainly study the influence of shape parameter $\beta$ on the performance of these methods. As a result, we fix $\sigma =1$, $\alpha =3$, and $\beta$ will be chosen for $2,3,4$, respectively. To show the comparison more vividly, we draw the following pictures (Figure 6, Figure 7 and Figure 8) and the estimated values of these three parameters are shown in Appendix A.4 (Part IV).

It can be seen that the two indicators RB and RMSE of these three estimation methods gradually decrease with the increase of sample size. As the value of the shape parameter, $\beta$, increases, the estimation accuracy of the IPWM and the MLE using the EM algorithm for shape parameter $\alpha$ gradually becomes better. Meanwhile, we can conclude that the IPWM performs best for both accuracy and stability, especially for the distribution has a heavier tail, while the MLE using the EM algorithm is the second on the whole. The performance of these two methods is comparable under a large sample size or high theoretical kurtosis. Furthermore, the IPWM is suitable for distributions with heavy tails and the MLE using the EM algorithm can become an alternative estimation method for steep-shaped distributions or under large sample size.

Second, a comparative study of the MLE using the EM algorithm and the MLE via a new iterative algorithm with sample kurtosis more than 2.7 is considered. The numerical results are presented in

Table 1. RB, RMSE and estimated values (in parentheses) of $(\sigma ,\alpha ,\beta )$.

Sample Size	MLE Using the EM Algorithm			MLE via a New Iterative Algorithm
	RB			RB
	$\mathbf{\sigma }$	$\mathbf{\alpha }$	$\mathbf{\beta }$	$\mathbf{\sigma }$	$\mathbf{\alpha }$	$\mathbf{\beta }$
$n=25$	0.8788 (0.879)	1.3800 (4.14)	1.2524 (3.757)	0.7235 (0.724)	0.7337 (2.201)	1.1098 (3.323)
$n=50$	0.9371 (0.937)	1.1867 (3.56)	1.2078 (3.624)	0.8101 (0.81)	0.8392 (2.518)	1.0586 (3.176)
$n=100$	0.9482 (0.948)	1.1041 (3.312)	1.1861 (3.558)	0.9144 (0.914)	0.9572 (2.872)	0.9989 (2.997)
$n=200$	0.9678 (0.968)	0.9555 (2.867)	1.1188 (3.356)	0.9373 (0.937)	1.0052 (3.016)	0.9837(2.951)
	RMSE			RMSE
	$\sigma$	$\alpha$	$\beta$	$\sigma$	$\alpha$	$\beta$
$n=25$	0.2348	1.2119	0.6866	0.2199	0.9791	0.1478
$n=50$	0.1632	0.9105	0.8034	0.1783	0.8141	0.09916
$n=100$	0.0874	0.7222	0.7311	0.1269	0.6223	0.0509
$n=200$	0.0401	0.2789	0.2175	0.0753	0.3971	0.0211

for fixed $\mu =0$, $\sigma =1$, $\alpha =3$, $\beta =3$.

Table 1 suggests that the new iterative algorithm performs better than the EM algorithm under sample kurtosis that exceeds 2.7 in general, especially for the estimation accuracy and stability of shape parameter $\beta$. For shape parameter $\alpha$, the iterative algorithm provides higher estimation accuracy while the EM algorithm is more stable.

At last, a simulation study for the PLA using the EM algorithm for $N\mathrm{Gt}(\mu ,\sigma ,\alpha ,\beta )$ ($\mu \ne 0$) considered in Section 5.6 will be performed. We mainly focus on the influence of shape parameter $\alpha$ and sample size on the estimation performance of this method. Therefore, we fix $\mu =1,\sigma =1,\beta =2$ and $\alpha$ will be chosen for $2,2.5,3$, respectively. Line graphs of the relationship between RB, RMSE and sample sizes with different value of shape parameter $\alpha$ were shown below (Figure 9) and the estimated value for these parameters are given in

Table A4. Estimation value for PLA using the EM algorithm under different $\alpha$ (fixed $\mu =1$, $\sigma =1$, $\beta =2$).

Value of $\mathit{\alpha }$	Sample Size	$\widehat{\mathit{\mu }}$	$\widehat{\mathit{\sigma }}$	$\widehat{\mathit{\alpha }}$	$\widehat{\mathit{\beta }}$
	n = 25	1.0082	1.4411	4.3878	4.8281
	n = 50	0.9986	1.3755	3.7841	3.3574
$\alpha =2$	n = 100	0.9997	1.2325	3.0712	2.6304
	n = 200	0.9999	1.1233	2.5676	2.2731
	n = 25	1.0036	1.3598	4.9423	4.6686
	n = 50	1.0029	1.2715	4.0851	3.4920
$\alpha =2.5$	n = 100	0.9994	1.1692	3.4725	2.6501
	n = 200	1.0006	1.0876	2.962	2.3084
	n = 25	0.9741	1.2494	5.0373	4.6552
	n = 50	0.9896	1.1827	4.3119	3.4574
$\alpha =3$	n = 100	1.0005	1.1058	3.7419	2.6380
	n = 200	0.954	1.0006	3.3182	2.2532

.

Figure 9 reveals that as the value of the shape parameter $\alpha$ increases, the estimation results become better and better. It is concluded that with the increase of the theoretical kurtosis, the accuracy and stability of the MLEs are gradually improved. In addition, the estimation effect of location parameter, $\mu$, and shape parameter, $\beta$, is less affected by the change of shape parameter $\alpha$, while that two indicators for the scale parameter, $\sigma$, and shape parameter, $\alpha$, are more sensitive to the change of it. This means that it is not very difficult to estimate the location parameter for a distribution with good symmetry, while the shape parameters of a heavy tail distribution are more difficult to estimate.

7. Real Data Analysis

In this section, a real data analysis is performed to illustrate the feasibility of the new Gt distribution and a comparison with the Gt distribution is also presented.

These data concerned U.S. stock market and were introduced by Shen et al. (2020), [26]. The basic description of the statistical characteristics of the data set are shown in

Table 2. Descriptive statistics for stock data.

n	$\overline{\mathit{Y}}$	S	$\sqrt{{\mathit{b}}_1}$	${\mathit{b}}_2$	min (Y)	max (Y)
4278	−2.295 × 10${}^{-3}$	0.01143	0.09466	5.28822	−0.1039	0.08113

, which containing the sample size n, the mean value $\overline{Y}$, the standard deviation S, the asymmetry coefficient $\sqrt{b_1}$, the kurtosis coefficient $b_2$, the minimum value min$\left(Y\right)$, and the maximum value max$\left(Y\right)$. It can be seen from the table that the data set has good symmetry and higher kurtosis.

To obtain the MLEs of the unknown parameters, the PLA depicted in Figure 10 is recommended to avoid possible numerical optimization problems. The fitting results of these two distributions and the estimated parameter values are presented in

Table 3. Estimation value and the fitting index.

Parameter Estimates	New Gt Distribution	Gt Distribution
$\sigma$	0.02143	0.01025
$\alpha$	5.46063	—
$\beta$	1.22367	—
p	—	1.92009
q	—	1.72641
log-likelihood value	13,563.13	13,543.77
AIC	−27,118.26	−27,079.54

. Based on the AIC criterion, we can infer from it that the new Gt distribution is slightly more suitable for fitting this data set than the Gt distribution. The histogram of the considered data with the two density functions we fitted is drawn in Figure 11. It is clearly seen that the pdf of the new Gt distribution is steeper in shape with a lighter tail, while the Gt distribution is the opposite.

8. Conclusions

In this paper, we coin a new Gt distribution that is suitable for fitting both the data with high kurtosis and a heavy tail based on a construction approach. Firstly, we have investigated the main properties of the new distribution including moments, skewness coefficients, kurtosis coefficients and random number generation. Secondly, we have derived the explicit expression for the moments of order statistics as well as its corresponding variance–covariance matrix through recurrence relations and the distribution transformation technique. The so obtained method has high efficiency and can greatly reduce the computation time. After that, we have focused on the parameter estimation of this new Gt distribution. Several estimation methods including MMOM, MLE using the EM algorithm, MLE using a new iterative algorithm and IPWM have been introduced. Among all these estimation methods, the IPWM performs best on the whole and this novel method makes the parameter estimation method of this distribution not limited to MLE and MMOM. Furthermore, the new iterative algorithm to acquire MLE is more suitable than the EM algorithm when the sample kurtosis is more than 2.7. For four-parameter new Gt distribution, we have established an EM-type algorithm through the profile maximum likelihood approach and discovered that the variation of the shape parameter $\alpha$ has a significant effect on the estimation performance of scale parameter $\sigma$ and shape parameter $\alpha$. However, there are still some limitations and areas to be improved in our research. The parameter estimation method is limited to PLA using the EM algorithm when the distribution has a heavier tail. Therefore, proposing more efficient and accurate estimation methods will be the focus of our future research. Besides, the new Gt distribution is suitable for fitting data with good symmetry. In future work, it can be applied for the asymmetric situation by adding a skew parameter.