Value at Risk and Expected Shortfall for Normal Weighted Inverse Gaussian Distributions

Calvin B. Maina^1*, Patrick G. O. Weke², Carolyne A. Ogutu², Joseph A. M. Ottieno²

¹Department of Mathematics and Actuarial Science, Kisii University, Kisii, Kenya.
²School of Mathematics, University of Nairobi, Nairob, Kenya.

DOI: 10.4236/jmf.2022.121002                

Abstract

Value at Risk (VaR) and Expected Shortfall (ES) is commonly used measures of potential risk for losses in financial markets. In literature VaR and ES for the Normal Inverse Gaussian (NIG) distribution, a special case of Generalized Hyperbolic Distribution (GHD), is frequently used. There are however, Normal Inverse Gaussian related distributions, which are also special cases of GHD that can also be used. The objective of this paper is to calculate VaR for Normal Weighted Inverse Gaussian (NWIG) distributions. The Expectation-Maximization (EM) algorithm has been used to obtain the Maximum Likelihood (ML) estimates of the proposed models for the Range Resource Corporation (RRC) financial data. We used Kupiec likelihood ratio (LR) for backtesting of VaR. Kolmogorov-Smirnov test and Anderson-Darling test have been used for goodness of fit test. Akaike Information Creterion (AIC), Bayesian Information Creterion (BIC) and Log-likelihood have been used for model selection. The results clearly show that the NWIG distributions are good alternatives to NIG for determining VaR and ES.

Keywords

Risk Measures, Backtesting, Weighted Distribution, Normal Mixture, EM-Algorithm

1. Introduction

The most popular measures for financial risk are Value at Risk (VaR) and Expected Shortfall (ES). VaR was proposed by Till Guldimann in the late 1980s while working for J. P. Morgan as the head of global research. It is generally defined as possible maximum loss over a given holding period within a fixed confidence level. An attractive feature of VaR is the backtestability of the measure. Backtesting is a method that uses historical data to gauge accuracy and effectiveness (Zhang and Nadarajah [1] ). Backtesting VaR is used to compare the forecast/predicted losses from the actual calculated losses realised at the end of a fixed time horizon. However, the main shortcoming of VaR is that it ignores any loss beyond the value at risk level. That is, it fails to capture tail risk. It also lacks a mathematical property called subadditivity as stated by Wimmerstedt [2]. That is, VaR for two combined portfolios can be larger than VaR for the sum of the two portfolios independently. This implies that diversification could increase risk, a contradiction to standard beliefs in finance. Artzner, Delbaen, Eber and Heath ( [3] [4] ) have proposed the use of Expected Shortfall (ES) also called conditional Value at Risk (CVaR) to circumvent the problems inherent in VaR. Expected Shortfall is the conditional expectation of loss given that the loss is beyond the VaR level. However, the main problem with ES is that it lacks a mathematical property called elicitability, Gneting [5], necessary for the backtestability of the risk measure. The Basel Committee [6] proposed to replace Value at Risk with Expected Shortfall but concluded that the backtesting will still be done on VaR even though the capital would be based on Expected Shortfall. Therefore the two measures of risk still remain the most popular and useful in financial management.

Notebly, Nadarajah et al. [7] have given a detailed review of VaR and ES for various distributions. One of the distributions reviewed is the Generalized Hyperbolic Distribution (GHD) introduced by Barndorff-Nielsen [8] as a Normal Variance-Mean Mixture with the Generalized Inverse Gaussian (GIG) distribution as the mixing distribution. The GIG is a three-parameter distribution denoted as $GIG\left(\lambda \mathrm{,}\delta \mathrm{,}\gamma \right)$. It embraces a number of special and limiting cases. The GHD and its subclasses fit high frequency financial data well which are characterised by skewness, excess kurtosis and fat tail. The most common special case is Normal Inverse Gaussian (NIG) distribution introduced by Barndorff-Nielsen [9] with the Inverse Gaussian (IG) as the mixing distribution. However, there are other special cases of GHD which are related to the NIG distribution which have not been considered for VaR an ES. These special cases are Normal Weighted Inverse Gaussian Distributions.

The objective of this paper is to determine VaR and ES for Range Resource Corporation (RRC) financial data using Normal Weighted Inverse Gaussian (NWIG) distributions. In particular we consider Normal mixtures with $GIG\left(\frac{1}{2}\mathrm{,}\delta \mathrm{,}\gamma \right)$, $GIG\left(-\frac{3}{2}\mathrm{,}\delta \mathrm{,}\gamma \right)$ and $GIG\left(\frac{3}{2}\mathrm{,}\delta \mathrm{,}\gamma \right)$ as mixing distribution. We study their properties and estimate parameters using the Expectation Maximization algorithm introduced by Dempster et al. [10]. Kolmogorov-Smirnov test and Anderson-Darling test have been used for goodness of fit test. Akaike Information Creterion (AIC), Bayesian Information Creterion (BIC) and Log-likelihood have been used for model selection.

The concept of a weighted distribution was introduced by Fisher [11] and elaborated by Patil and Rao [12]. Reciprocal Inverse Gaussian and the finite mixture of Inverse Gaussian and Reciprocal Inverse Gaussian distribution are shown to be Weighted Inverse Gaussian (WIG) distributions by Akman and Gupta [13]; Gupta and Akman [14]; Gupta and Kundu [15]. The Special cases $GIG\left(-\frac{3}{2}\mathrm{,}\delta \mathrm{,}\gamma \right)$ and $GIG\left(\frac{3}{2}\mathrm{,}\delta \mathrm{,}\gamma \right)$ are also WIG distributions.

For value at Risk of these models we use the Kupiec likelihood ratio (LR) introduced by Kupiec [16]. The remainder of this paper is structured as follows: Section 2 deals with the concept of VaR and ES. Section 3 deals with Generalised Inverse Gaussian and its special cases. Weighted distribution is covered in Section 4. The concept on Generalised hyperbolic distribution is illustrated in Section 5 while its special cases of interest are in Section 6. Parameter estimation of the proposed models is performed in Section 7. Application to Range Resource Corporation is done in Section 8 and Section 8 deals with conclusion.

2. Value at Risk and Expected Shortfall: Mathematical Background

The most important risk measures despite their drawbacks are Value at Risk (VaR) and Expected Shortfall (ES). VaR was proposed by Till Guldimann in the late 1980s, and at the time he was the head of global research at J. P. Morgan.

Value at Risk is generally defined as possible maximum loss over a given holding period within a fixed confidence level. Mathematically VaR at the ( $100-\alpha$ ) percent confidence level is defined as the lower $100\alpha$ percentile of the profit-loss distribution.

In statistical terms, VaR is a quantile of distribution for financial asset returns. More formally, VaR is defined as

$P\left\{X\le -Va{R}_{1-\alpha }^X\right\}=\alpha$ (1)

where X represents the Asset’s returns. In integral form it can be expressed as

${\displaystyle {\int }_{-\infty }^{Va{R}_{\alpha }^X}}f\left(x\right)\text{d}x=\alpha$ (2)

where f(x) is the profit-loss distribution.

The concept of Expected Shortfall (ES) was first introduced in Rappoport [17]. Artzner et al. ( [3] [4] ) formally developed the concept. ES is the conditional expectation of loss given that the loss is beyond the VaR level and measures how much one can lose on average in the states beyond the VaR level.

From Equation (2.2)

$\frac{1}{\alpha }{\displaystyle {\int }_{-\infty }^{Va{R}_{\alpha }}}f\left(x\right)\text{d}x=1$ (3)

Therefore $\frac{f\left(x\right)}{\alpha }$ is a pdf for $-\infty <x<Va{R}_{\alpha }$ and we refer to it as “Tail loss distribution”.

Conditional Expectation

$E\left[X|X<Va{R}_{\alpha }\right]={\displaystyle {\int }_{-\infty }^{Va{R}_{\alpha }}}x\frac{f\left(x\right)}{\alpha }\text{d}x$ (4)

is the Expected Shortfall denoted as $E{S}_{\alpha }$. This version was used by Yamai and Yoshiba [18] to obtain the ES for a normal distribution. Equation (2.4) can be expressed in a different version as follows: Defining F(x) as the cdf of the random variable X, let

$u=F\left(x\right)\Rightarrow x=F^{-1}(\; u\; )$

$\therefore \text{d}u=f\left(x\right)\text{d}x$

when

$x=-\infty \Rightarrow u=0$

$x=Va{R}_{\alpha }\Rightarrow u=\alpha$

$\therefore E{S}_{\alpha }=\frac{1}{\alpha }{\displaystyle {\int }_0^{\alpha }}{F}^{-1}\left(u\right)\text{d}u=\frac{1}{\alpha }{\displaystyle {\int }_0^{\alpha }}Va{R}_u\text{d}u$ (5)

as presented by Zhang et al. [19].

Remarks: Equation (2.3) is the mean of the loss distribution. Equation (2.4) represents the average of the VaR between 0 and $\alpha$. The loss distribution, $\frac{f\left(x\right)}{\alpha }$, $-\infty <x<Va{R}_{\alpha }$ gives the tail distribution.

For the purpose of VaR and ES analysis, a model for the return distribution is important because it describes the potential behaviour of a financial security in the future Bams and Wilhouwer [20]. A Normal distribution supposedly underestimates the tail and hence VaR. Recently alternative distributions have been proposed that focus more on tail behaviour of the returns. One such candidate is the Normal Inverse Gaussian (NIG) distribution. We consider extensions of NIG distribution as Normal Weighted Inverse Gaussian (NWIG) distributions. In the next few sections we give a detailed illustration on their construction, properties and parameter estimation via EM-algorithm.

3. Generalised Inverse Gaussian Distribution

The Generalised Inverse Gaussian (GIG) Distribution is based on modified Bessel function of the third kind. Modified Bessel function of the third kind of order $\lambda$ evaluated at $\omega$ denoted by $K_{\lambda }\left(\omega \right)$ is defined as

$K_{\lambda }\left(\omega \right)=\frac{1}{2}{\displaystyle {\int }_0^{\infty }}{x}^{\lambda -1}{\text{e}}^{-\frac{\omega }{2}\left(x+\frac{1}{x}\right)}\text{d}x$ (6)

with the following properties

a) $K_{\frac{1}{2}}\left(\omega \right)=K_{-\frac{1}{2}}\left(\omega \right)=\sqrt{\frac{\pi }{2\omega }}{\text{e}}^{-\omega }$ (7)

b) $K_{\frac{3}{2}}\left(\omega \right)=K_{-\frac{3}{2}}\left(\omega \right)=\sqrt{\frac{\pi }{2\omega }}{\text{e}}^{-\omega }\left(1+\frac{1}{\omega }\right)$ (8)

c) $K_{\frac{5}{2}}\left(\omega \right)=K_{-\frac{5}{2}}\left(\omega \right)=\sqrt{\frac{\pi }{2\omega }}{\text{e}}^{-\omega }\left(1+\frac{3}{\omega }+\frac{3}{{\omega }^2}\right)$ (9)

d) $K_{\frac{7}{2}}\left(\omega \right)=K_{-\frac{7}{2}}\left(\omega \right)=\sqrt{\frac{\pi }{2\omega }}{\text{e}}^{-\omega }\left(1+\frac{6}{\omega }+\frac{15}{{\omega }^2}+\frac{15}{{\omega }^3}\right)$ (10)

e) $K_{\frac{9}{2}}\left(\omega \right)=K_{-\frac{9}{2}}\left(\omega \right)=\sqrt{\frac{\pi }{2\omega }}{\text{e}}^{-\omega }\left(1+\frac{10}{\omega }+\frac{45}{{\omega }^2}+\frac{105}{{\omega }^3}+\frac{105}{{\omega }^4}\right)$ (11)

f) $K_{\frac{11}{2}}\left(\omega \right)=K_{-\frac{11}{2}}\left(\omega \right)=\sqrt{\frac{\pi }{2\omega }}{\text{e}}^{-\omega }\left(1+\frac{15}{\omega }+\frac{105}{{\omega }^2}+\frac{420}{{\omega }^3}+\frac{945}{{\omega }^4}+\frac{945}{{\omega }^5}\right)$ (12)

which are necessary in deriving the properties and estimates of the proposed models. For more definition and properties see Abramowitz and Stegun [21].

Using Parametrization $\omega =\delta \gamma$ and transformation $x=\frac{\gamma }{\delta }z$ then formula (2.1) becomes

$K_{\lambda }\left(\delta \gamma \right)=\frac{1}{2}{\displaystyle {\int }_0^{\infty }}{\left(\frac{\gamma }{\delta }\right)}^{\lambda }{z}^{\lambda -1}\exp \left\{-\frac{1}{2}\left(\frac{{\delta }^2}{z}+{\gamma }^{2}z\right)\right\}\text{d}z$

Hence

$g\left(z\right)={\left(\frac{\gamma }{\delta }\right)}^{\lambda }\frac{z^{\lambda -1}}{2K_{\lambda }\left(\delta \gamma \right)}\exp \left\{-\frac{1}{2}\left(\frac{{\delta }^2}{z}+{\gamma }^{2}z\right)\right\}$ (13)

$z>0;-\infty <\lambda <\infty ,\delta >0,\gamma >0$

which is a Generalized Inverse Gaussian (GIG) distribution with parameter $\lambda \mathrm{,}\delta \mathrm{,}\gamma$.

Thus

$Z\sim GIG\left(\lambda \mathrm{,}\delta \mathrm{,}\gamma \right)$

with

$E\left(Z^r\right)={\left(\frac{\delta }{\gamma }\right)}^r\frac{K_{\lambda +r}\left(\delta \gamma \right)}{K_{\lambda }(\; \delta \; \gamma \; )}$

where r can be positive or negative integers.

4. Weighted Inverse Gaussian Distribution

Let X be a random variable with pdf $f\left(x\right)$. A function of X, $w\left(X\right)$ is also a random variable with expectation

$E\left[w\left(X\right)\right]={\displaystyle {\int }_{-\infty }^{\infty }}w\left(x\right)f\left(x\right)\text{d}x$

$\therefore 1={\displaystyle {\int }_{-\infty }^{\infty }}\frac{w\left(x\right)}{E\left[w\left(X\right)\right]}f\left(x\right)\text{d}x$

Thus

$f_W\left(x\right)=\frac{w\left(x\right)}{E\left[w\left(X\right)\right]}f\left(x\right),-\infty <x<\infty$ (14)

is a weighted distribution. It was introduced by Fisher [11] and elaborated by Patil and Rao [12].

From Equation (3.8) When $\lambda =-\frac{1}{2}$ we have $GIG\left(-\frac{1}{2}\mathrm{,}\delta \mathrm{,}\gamma \right)$ with pdf

$g_1\left(z\right)=\frac{\delta {\text{e}}^{\delta \gamma }}{\sqrt{2\pi }}{z}^{-\frac{3}{2}}\exp \left\{-\frac{1}{2}\left(\frac{{\delta }^2}{z}+{\gamma }^{2}z\right)\right\}$ (15)

This is called Inverse Gaussian (IG) distribution.

The following special cases of GIG distribution can be expressed in terms of IG distribution. They are weighted Inverse Gaussian distributions.

Example 1: $GIG\left(\frac{1}{2}\mathrm{,}\delta \mathrm{,}\gamma \right)$

When $\lambda =\frac{1}{2}$, then the pdf of $GIG\left(\lambda \mathrm{,}\delta \mathrm{,}\gamma \right)$ becomes

$g_2\left(z\right)=\frac{\gamma {\text{e}}^{\delta \gamma }}{\sqrt{2\pi }}{z}^{-\frac{1}{2}}\exp \left\{-\frac{1}{2}\left(\frac{{\delta }^2}{z}+{\gamma }^{2}z\right)\right\}$

This is a Reciprocal Inverse Gaussian (RIG) distribution. It can be written as

$g_2\left(z\right)=\frac{\gamma }{\delta }z\left[\frac{\delta {\text{e}}^{\delta \gamma }}{\sqrt{2\pi }}{z}^{-\frac{3}{2}}\exp \left\{-\frac{1}{2}\left(\frac{{\delta }^2}{z}+{\gamma }^{2}z\right)\right\}\right]$

$\therefore g_2\left(z\right)=\left(\frac{\gamma }{\delta }z\right)g_1(\; z\; )$

i.e.,

$GIG\left(\frac{1}{2}\mathrm{,}\delta \mathrm{,}\gamma \right)=\left(\frac{\gamma }{\delta }z\right)GIG\left(-\frac{1}{2}\mathrm{,}\delta \mathrm{,}\gamma \right)$ (16)

Thus a $GIG\left(\frac{1}{2}\mathrm{,}\delta \mathrm{,}\gamma \right)$ is a weighted Inverse Gaussian distribution with weights

$w\left(Z\right)=Z$ (17)

where

$Z\sim GIG\left(-\frac{1}{2}\mathrm{,}\delta \mathrm{,}\gamma \right)$

$\therefore E\left[w\left(Z\right)\right]=E\left(Z\right)=\frac{\delta }{\gamma }$

RIG distribution is called Length Biased Inverse Gaussian distribution.

Example 2: $GIG\left(-\frac{3}{2}\mathrm{,}\delta \mathrm{,}\gamma \right)$

The pdf is

$g_3\left(z\right)=\frac{{\delta }^3{\text{e}}^{\delta \gamma }}{\sqrt{2\pi }\left(1+\delta \gamma \right)}{z}^{-\frac{5}{2}}{\text{e}}^{-\frac{1}{2}\left(\frac{{\delta }^2}{z}+{\gamma }^{2}z\right)}=\frac{{\delta }^2}{1+\delta \gamma }{z}^{-1}{g}_1(\; z\; )$

Therefore $GIG\left(-\frac{3}{2}\mathrm{,}\delta \mathrm{,}\gamma \right)$ is a weighted Inverse Gaussian distribution with weights

$w\left(Z\right)=Z^{-1}$ (18)

and

$E\left[w\left(Z\right)\right]=\frac{1+\delta \gamma }{{\delta }^2}$

Example 3: $GIG\left(\frac{3}{2}\mathrm{,}\delta \mathrm{,}\gamma \right)$

The pdf of $GIG\left(\frac{3}{2}\mathrm{,}\delta \mathrm{,}\gamma \right)$

$\begin{array}{c}{g}_4\left(z\right)=\frac{{\gamma }^3}{1+\delta \gamma }\frac{{\text{e}}^{\delta \gamma }}{\sqrt{2\pi }}{z}^{\frac{1}{2}}{\text{e}}^{-\frac{1}{2}\left(\frac{{\delta }^2}{z}+{\gamma }^{2}z\right)}\\ =\frac{{\gamma }^3}{\delta \left(1+\delta \gamma \right)}{z}^2\left[\frac{\delta {\text{e}}^{\delta \gamma }}{\sqrt{2\pi }}{z}^{-\frac{3}{2}}\exp \left\{-\frac{1}{2}\left(\frac{{\delta }^2}{z}+{\gamma }^{2}z\right)\right\}\right]\\ =\frac{{\gamma }^3}{\delta \left(1+\delta \gamma \right)}{z}^2{g}_1(\; z\; )\end{array}$

Thus, $GIG\left(\frac{3}{2}\mathrm{,}\delta \mathrm{,}\gamma \right)$ is a weighted Inverse Gaussian distribution with weights

$w\left(Z\right)=Z^2$ (19)

where

$Z\sim GIG\left(-\frac{1}{2}\mathrm{,}\delta \mathrm{,}\gamma \right)$

$\therefore E\left[w\left(Z\right)\right]=\frac{\delta \left(1+\delta \gamma \right)}{{\gamma }^3}$

5. Generalized Hyperbolic Distribution

A stochastic representation of a Normal Variance-Mean mixture is given by letting Let

$X=\mu +\beta Z+\sqrt{Z}Y$

where

$Y\sim N(\; 0,1\; )$

and Z, independent of Y, is a positive random variable.

If F(x) is a cdf of X, then

$\begin{array}{c}F\left(x\right)=prob\left\{X\le x\right\}\\ =\left\{Y\le \frac{x-\mu -\beta z}{\sqrt{z}},0<z<\infty \right\}\\ ={\displaystyle {\int }_0^{\infty }{\displaystyle {\int }_{-\infty }^{\frac{x-\mu -\beta z}{\sqrt{z}}}\varphi \left(y\right)g\left(z\right)\text{d}y\text{d}z}}\\ ={\displaystyle {\int }_0^{\infty }\Phi \left(\frac{x-\mu -\beta z}{\sqrt{z}}\right)g\left(z\right)\text{d}z}\end{array}$

where $\varphi (\cdot )$ and $\Phi (\cdot )$ are pdf and cdf of a standard normal distribution, respectively.

$\therefore f\left(x\right)={\displaystyle {\int }_0^{\infty }}\frac{1}{\sqrt{z}}\varphi \left(\frac{x-\mu -\beta z}{\sqrt{z}}\right)g\left(z\right)\text{d}z={\displaystyle {\int }_0^{\infty }}\frac{1}{\sqrt{2\pi z}}{\text{e}}^{-\frac{{\left[x-\left(\mu +\beta z\right)\right]}^2}{2z}}g\left(z\right)\text{d}z$ (20)

Thus we have a hierarchical representation as

$X/Z=z\sim N\left(\mu +\beta z\mathrm{,}z\right)$ (21)

being the conditional pdf and g(z) the mixing distribution.

If

$Z\sim GIG\left(\lambda \mathrm{,}\delta \mathrm{,}\gamma \right)$ (22)

then

$f\left(x\right)=\frac{{\left(\frac{\gamma }{\delta }\right)}^{\lambda }{\text{e}}^{\beta \left(x-\mu \right)}}{\sqrt{2\pi }{K}_{\lambda }\left(\delta \gamma \right)}\frac{1}{2}{\displaystyle {\int }_0^{\infty }}{z}^{\lambda -\frac{3}{2}}{\text{e}}^{-\frac{1}{2}\left({\alpha }^{2}z+\frac{{\delta }^2}{z}\varphi \left(x\right)\right)}\text{d}z$

where

${\alpha }^2={\gamma }^2+{\beta }^2\Rightarrow \gamma =\sqrt{{\alpha }^2-{\beta }^2}$

and

$\varphi \left(x\right)=1+\frac{{\left(x-\mu \right)}^2}{{\delta }^2}$

This is an integral form of a GHD

for $-\infty <x<\infty ;-\infty <\mu <\infty ,\delta >0,\alpha >0,\beta >0$.

Using the stochastic representation of NVM mixture, the properties are expressed in terms of the mixing distribution as shown in the following

Proposition 1

$M_X\left(t\right)={\text{e}}^{\mu t}{M}_Z\left(\beta t+\frac{t^2}{2}\right)$

$E\left(X\right)=\mu +\beta E(\; Z\; )$

$Var\left(X\right)=E\left(Z\right)+{\beta }^{2}Var(\; Z\; )$

${\mu }_3\left(X\right)=3\beta var\left(Z\right)+{\beta }^3{\mu }_3(\; Z\; )$

${\mu }_4\left(X\right)={\beta }^4{\mu }_4\left(Z\right)+6{\beta }^2{\mu }_3\left(Z\right)+6{\beta }^{2}E\left[Z\right]var\left(Z\right)+3E\left[Z^2\right]$

6. Special Cases of Interest

When $\lambda =-\frac{1}{2}$, the GHD becomes NIG whose pdf can be expressed in the integral form as

$f_1\left(x\right)=\frac{\delta {\text{e}}^{\delta \gamma }{\text{e}}^{\beta \left(x-\mu \right)}}{2\pi }{\displaystyle {\int }_0^{\infty }}{z}^{-2}{\text{e}}^{-\frac{1}{2}\left\{\frac{{\delta }^2{\varphi }_x}{z}+{\alpha }^{2}z\right\}}\text{d}z$ (23)

Therefore the Normal Weighted Inverse Gaussian formulation can be represented as

$\therefore f\left(x\right)=\frac{\delta {\text{e}}^{\delta \gamma }{\text{e}}^{\beta \left(x-\mu \right)}}{2\pi }{\displaystyle {\int }_0^{\infty }}\frac{w\left(z\right)}{E\left[w\left(Z\right)\right]}{z}^{-2}{\text{e}}^{-\frac{1}{2}\left\{\frac{{\delta }^2{\varphi }_x}{z}+{\alpha }^{2}z\right\}}\text{d}z$ (24)

Hence

$\lambda =-\frac{1}{2}\Rightarrow \frac{w\left(z\right)}{E\left(w\left(Z\right)\right)}=1$ (25)

$\therefore f_1\left(x\right)=\frac{\alpha \delta {\text{e}}^{\delta \gamma +\beta \left(x-\mu \right)}{K}_1\left(\alpha \sqrt{{\delta }^2+{\left(x-\mu \right)}^2}\right)}{\pi \sqrt{{\delta }^2+{\left(x-\mu \right)}^2}}$

which is the Normal Inverse Gaussian (NIG) distribution with properties given in Table 1 below.

when

$\lambda =\frac{1}{2}\Rightarrow \frac{w\left(z\right)}{E\left[w\left(Z\right)\right]}=\frac{\gamma }{\delta }z$ (26)

$f_2\left(x\right)=\frac{\gamma {\text{e}}^{\delta \gamma +\beta \left(x-\mu \right)}{K}_0\left(\alpha \sqrt{{\delta }^2+{\left(x-\mu \right)}^2}\right)}{\pi }$

Which is the Normal Reciprocal Inverse Gaussian (NRIG) distribution with properties given in Table 2 below.

Similarly, when

$\lambda =-\frac{3}{2}\Rightarrow \frac{w\left(z\right)}{E\left[w\left(Z\right)\right]}=\frac{{\delta }^2}{\left(1+\delta \gamma \right)z}$ (27)

$\therefore f_3\left(x\right)=\frac{{\alpha }^2{\delta }^3{\text{e}}^{\delta \gamma +\beta \left(x-\mu \right)}{K}_2\left(\alpha \sqrt{{\delta }^2+{\left(x-\mu \right)}^2}\right)}{\pi \left(1+\delta \gamma \right)\left[{\delta }^2+{\left(x-\mu \right)}^2\right]}$

which is the Normal- $GIG\left(\lambda =-\frac{3}{2},\delta ,\gamma \right)$ with properties given in Table 3 below.

Table 1. Properties of NIG.

Table 2. Properties of NRIG.

Table 3. Properties for Normal- $GIG\left(-\frac{3}{2}\mathrm{,}\delta \mathrm{,}\gamma \right)$.

Finally, when

$\lambda =\frac{3}{2}\Rightarrow \frac{w\left(z\right)}{E\left[w\left(Z\right)\right]}=\frac{{\gamma }^3{z}^2}{\delta \left(1+\delta \gamma \right)}$ (28)

$\therefore f_4\left(x\right)=\frac{{\gamma }^3{\text{e}}^{\delta \gamma +\beta \left(x-\mu \right)}\left[\sqrt{{\delta }^2+{\left(x-\mu \right)}^2}\right]K_1\left(\alpha \sqrt{{\delta }^2+{\left(x-\mu \right)}^2}\right)}{\alpha \pi \left(1+\delta \gamma \right)}$

which is the Normal- $GIG\left(\lambda =\frac{3}{2},\delta ,\gamma \right)$ with properties given in Table 4 below

7. Parameter Estimation via EM Algorithm

7.1. Theory/Concept

EM algorithm is a powerful technique for maximum likelihood estimation for data containing missing values or data that can be considered as containing missing values. It was introduced by Dempster et al. [10].

Table 4. Properties for Normal- $GIG\left(\frac{3}{2}\mathrm{,}\delta \mathrm{,}\gamma \right)$.

Assume that the true data are made of an observed part X and unobserved part Z. This then ensures the log likelihood of the complete data $\left(x_i\mathrm{,}{z}_i\right)$ for $i=1,2,3,\cdots ,n$ factorizes into two parts, Kostas [22].

This implies that the joint density of X and Z is given by

$f\left(x\mathrm{,}z\right)=f\left(x/z\right)g(\; z\; )$

The likelihood function is

$L={\displaystyle \underset{i=1}{\overset{n}{\prod }}}f\left(x_i/z_i\right)g\left(z_i\right)={\displaystyle \underset{i=1}{\overset{n}{\prod }}}f\left(x_i/z_i\right){\displaystyle \underset{i=1}{\overset{n}{\prod }}}g(\; z\; i\; )$

$\begin{array}{l}\therefore \log L=\log {\displaystyle \underset{i=1}{\overset{n}{\prod }}}f\left(x_i/z_i\right)+\log {\displaystyle \underset{i=1}{\overset{n}{\prod }}}g\left(z_i\right)\\ ={\displaystyle \underset{i=1}{\overset{n}{\sum }}}\log f\left(x_i/z_i\right)+{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\log g\left(z_i\right)\\ =l_1+l_2\end{array}$

where

$l_1={\displaystyle \underset{i=1}{\overset{n}{\sum }}}\log f\left(x_i/z_i\right)$

and

$l_2={\displaystyle \underset{i=1}{\overset{n}{\sum }}}\log g(\; z\; i\; )$

Karlis [23] applied EM algorithm to mixtures which he considered to consist of two parts; the conditional pdf is for observed data and the mixing distribution is based on an unobserved data, the missing values.

7.2. M-Step for Conditional pdf

In this study

$l_1=-\frac{n}{2}\log \left(2\pi \right)-\frac{1}{2}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\log z_i-{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\frac{{\left(x_i-\mu -\beta z_i\right)}^2}{2z_i}$

Therefore

$\frac{\partial }{\partial \beta }{l}_1=0\Rightarrow {\displaystyle \underset{i=1}{\overset{n}{\sum }}}\left(x_i-\hat{\mu }-\hat{\beta }{z}_i\right)=0$

i.e., ${\displaystyle \underset{i=1}{\overset{n}{\sum }}}{x}_i-n\hat{\mu }-\hat{\beta }{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{z}_i=0$

$\therefore \hat{\mu }=\bar{x}-\hat{\beta }\bar{z}$

where $\bar{x}={\displaystyle {\sum }_{i=1}^n}\frac{x_i}{n}$ and $\bar{z}={\displaystyle {\sum }_{i=1}^n}\frac{z_i}{n}$.

Similarly,

$\frac{\partial }{\partial \mu }{l}_1=0\Rightarrow {\displaystyle \underset{i=1}{\overset{n}{\sum }}}\frac{x_i}{z_i}-\hat{\mu }{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\frac{1}{z_i}-n\hat{\beta }=0$

$\therefore {\displaystyle \underset{i=1}{\overset{n}{\sum }}}\frac{x_i}{z_i}-\bar{x}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\frac{1}{z_i}+\hat{\beta }\bar{z}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\frac{1}{z_i}-n\hat{\beta }=0$

$\therefore \hat{\beta }=\frac{{\displaystyle {\sum }_{i=1}^n}\frac{x_i}{z_i}-\bar{x}{\displaystyle {\sum }_{i=1}^n}\frac{1}{z_i}}{n-\bar{z}{\displaystyle {\sum }_{i=1}^n}\frac{1}{z_i}}$

7.3. E-Step

Values of random variables $Z_i$ and $\frac{1}{Z_i}$ are not known. So we estiamte them by considering posterior expectations $E\left(Z_i/X_i\right)$ and $E\left(\frac{1}{Z_i}\right)$.

The posterior distribution is defined as

$f\left(z/x\right)=\frac{f\left(x/z\right)g\left(z\right)}{{\displaystyle {\int }_0^{\infty }}f\left(x/z\right)g\left(z\right)\text{d}z}$

Therefore

$E\left(Z/X\right)=\frac{{\displaystyle {\int }_0^{\infty }}zf\left(x/z\right)g\left(z\right)\text{d}z}{{\displaystyle {\int }_0^{\infty }}f\left(x/z\right)g\left(z\right)\text{d}z}$

$E\left(Z/X\right)=\frac{{\displaystyle {\int }_0^{\infty }}\frac{1}{z}f\left(x/z\right)g\left(z\right)\text{d}z}{{\displaystyle {\int }_0^{\infty }}f\left(x/z\right)g\left(z\right)\text{d}z}$

7.4. For Inverse Gaussian Mixing Distribution

7.4.1. M-Step

$l_2=-\frac{n}{2}\log \left(2\pi \right)+n\log \delta +n\delta \gamma -\frac{3}{2}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\log z_i-\frac{{\delta }^2}{2}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\frac{1}{z_i}-\frac{{\gamma }^2}{2}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{z}_i$

$\therefore \frac{\partial }{\partial \gamma }{l}_2=n\delta -\gamma {\displaystyle \underset{i=1}{\overset{n}{\sum }}}{z}_i$

$\therefore \frac{\partial }{\partial \gamma }{l}_2=0\Rightarrow \hat{\gamma }=\frac{\hat{\delta }}{\bar{z}}$

and

$\frac{\partial }{\partial \delta }{l}_2=0\Rightarrow \frac{n}{\delta }+n\gamma -\delta {\displaystyle \underset{i=1}{\overset{n}{\sum }}}\frac{1}{z_i}=0$

$\therefore \hat{\delta }=\sqrt{\frac{n}{{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\frac{1}{z_i}-\frac{n}{\bar{z}}}}$

7.4.2. Posterior Expectations

For NIG

$E\left(Z/X\right)=\frac{\delta \sqrt{\varphi \left(x\right)}}{\alpha }\frac{K_0\left(\alpha \delta \sqrt{\varphi \left(x\right)}\right)}{K_1\left(\alpha \delta \sqrt{\varphi (\; x\; )}\right)}$

and

$E\left(\frac{1}{Z}/X\right)=\frac{\alpha }{\delta \sqrt{\varphi \left(x\right)}}\frac{K_2\left(\alpha \delta \sqrt{\varphi \left(x\right)}\right)}{K_1\left(\alpha \delta \sqrt{\varphi (\; x\; )}\right)}$

7.4.3. Iterations

Let

$s_i=E\left(Z_i/X_i\right)=\frac{\delta \sqrt{\varphi \left(x_i\right)}}{\alpha }\frac{K_0\left(\alpha \delta \sqrt{\varphi \left(x_i\right)}\right)}{K_1\left(\alpha \delta \sqrt{\varphi \left(x_i\right)}\right)}$

where

$\varphi \left(x_i\right)=1+\frac{{\left(x_i-\mu \right)}^2}{{\delta }^2}$

and

$\alpha =\sqrt{{\gamma }^2+{\beta }^2}$

Therefore $s_i$ is a function of $\alpha \mathrm{,}\beta \mathrm{,}\delta \mathrm{,}\gamma \mathrm{,}\mu$

Next, let

$w_i=E\left(\frac{1}{Z_i}/X_i\right)=\frac{\alpha }{\delta \sqrt{\varphi \left(x_i\right)}}\frac{K_2\left(\alpha \delta \sqrt{\varphi \left(x_i\right)}\right)}{K_1\left(\alpha \delta \sqrt{\varphi \left(x_i\right)}\right)}$

which is also a function of $\alpha \mathrm{,}\beta \mathrm{,}\delta \mathrm{,}\gamma \mathrm{,}\mu$

For computation, we have the following in the k-th iteration

$s_i^{\left(k\right)}=\frac{{\delta }^{\left(k\right)}\sqrt{{\varphi }^{\left(k\right)}\left(x_i\right)}}{{\alpha }^{\left(k\right)}}\frac{K_0\left({\alpha }^{\left(k\right)}{\delta }^{\left(k\right)}\sqrt{{\varphi }^{\left(k\right)}\left(x_i\right)}\right)}{K_1\left({\alpha }^{\left(k\right)}{\delta }^{\left(k\right)}\sqrt{{\varphi }^{\left(k\right)}\left(x_i\right)}\right)}$ (29)

$w_i^{\left(k\right)}=\frac{{\alpha }^{\left(k\right)}}{{\delta }^{\left(k\right)}\sqrt{{\varphi }^{\left(k\right)}\left(x_i\right)}}\frac{K_2\left({\alpha }^{\left(k\right)}{\delta }^{\left(k\right)}\sqrt{{\varphi }^{\left(k\right)}\left(x_i\right)}\right)}{K_1\left({\alpha }^{\left(k\right)}{\delta }^{\left(k\right)}\sqrt{{\varphi }^{\left(k\right)}\left(x_i\right)}\right)}$ (30)

${\hat{\delta }}^{\left(k+1\right)}=\sqrt{\frac{n}{{\displaystyle {\sum }_{i=1}^n}{w}_i^{\left(k\right)}-\frac{n^2}{{\displaystyle {\sum }_{i=1}^n}{s}_i^{\left(k\right)}}}}$ (31)

${\hat{\gamma }}^{\left(k+1\right)}=\frac{{\hat{\delta }}^{\left(k+1\right)}}{\frac{1}{n}{\displaystyle {\sum }_{i=1}^n}{s}_i^{\left(k\right)}}$ (32)

The log-likelihood function of NIG distribution is given by

$\begin{array}{c}\log L={\displaystyle \underset{i=1}{\overset{n}{\sum }}}\log \left\{\frac{\alpha \delta {\text{e}}^{\delta \gamma +\beta \left(x-\mu \right)}{K}_1\left(\alpha \sqrt{{\delta }^2+{\left(x-\mu \right)}^2}\right)}{\pi \sqrt{{\delta }^2+{\left(x-\mu \right)}^2}}\right\}\\ =n\log \alpha +n\left(\delta \gamma -\beta \mu \right)+\beta {\displaystyle \underset{i=1}{\overset{n}{\sum }}}{x}_i-n\log \pi \\ -\frac{1}{2}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\varphi \left(x_i\right)+{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\log K_1\left(\alpha \delta \varphi {\left(x_i\right)}^{\frac{1}{2}}\right)\end{array}$ (33)

The k-th iteration of the loglikelihood function of the NIG distribution is

$\begin{array}{c}{l}^{\left(k\right)}=n\log {\alpha }^{\left(k\right)}+n\left({\delta }^{\left(k\right)}{\gamma }^{\left(k\right)}-{\beta }^{\left(k\right)}{\mu }^{\left(k\right)}\right)+{\beta }^{\left(k\right)}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{x}_i-n\log \pi \\ -\frac{1}{2}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\varphi }^{\left(k\right)}\left(x_i\right)+{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\log K_1\left({\alpha }^{\left(k\right)}{\delta }^{\left(k\right)}{\varphi }^{\left(k\right)}{\left(x_i\right)}^{\frac{1}{2}}\right)\end{array}$ (34)

7.5. For Length Biased (Reciprocal) Inverse Gaussian Distribution

7.5.1. M-Step

$l_2\left(\delta ,\gamma \right)=n\log \gamma -\frac{n}{2}\log \left(2\pi \right)+n\delta \gamma -\frac{{\delta }^2}{2}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\frac{1}{z_i}-\frac{1}{2}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\log z_i$

$\therefore \frac{\partial }{\partial \delta }{l}_2=0\Rightarrow \hat{\delta }=\frac{n\hat{\gamma }}{{\displaystyle {\sum }_{i=1}^n}\frac{1}{z_i}}$

$\frac{\partial }{\partial \gamma }{l}_2=0\Rightarrow \hat{\gamma }=\sqrt{\frac{n}{{\displaystyle {\sum }_{i=1}^n}{z}_i-n^2{\left({\displaystyle {\sum }_{i=1}^n}\frac{1}{z_i}\right)}^{-1}}}$

7.5.2. Posterior Expectation

$E\left(Z/X\right)=\frac{\delta \sqrt{\varphi \left(x\right)}}{\alpha }\frac{K_1\left(\alpha \delta \sqrt{\varphi \left(x\right)}\right)}{K_0\left(\alpha \delta \sqrt{\varphi \left(x\right)}\right)}$ (35)

and

$E\left(\frac{1}{Z}/X\right)=\frac{\alpha }{\delta \sqrt{\varphi \left(x\right)}}\frac{K_1\left(\alpha \delta \sqrt{\varphi \left(x\right)}\right)}{K_0\left(\alpha \delta \sqrt{\varphi \left(x\right)}\right)}$ (36)

7.5.3. Iterations

Let

$s_i=E\left(Z_i/X_i\right)=\frac{\delta \sqrt{\varphi \left(x_i\right)}}{\alpha }\frac{K_1\left(\alpha \delta \sqrt{\varphi \left(x_i\right)}\right)}{K_0\left(\alpha \delta \sqrt{\varphi \left(x_i\right)}\right)}$

and

$w_i=E\left(\frac{1}{Z_i}/X_i\right)=\frac{\alpha }{\delta \sqrt{\varphi \left(x_i\right)}}\frac{K_1\left(\alpha \delta \sqrt{\varphi \left(x_i\right)}\right)}{K_0\left(\alpha \delta \sqrt{\varphi \left(x_i\right)}\right)}$

For iterations we have

$s_i^{\left(k\right)}=\frac{{\delta }^{\left(k\right)}\sqrt{{\varphi }^{\left(k\right)}\left(x_i\right)}}{{\alpha }^{\left(k\right)}}\frac{K_1\left({\alpha }^{\left(k\right)}{\delta }^{\left(k\right)}\sqrt{{\varphi }^{\left(k\right)}\left(x_i\right)}\right)}{K_0\left({\alpha }^{\left(k\right)}{\delta }^{\left(k\right)}\sqrt{{\varphi }^{\left(k\right)}\left(x_i\right)}\right)}$ (37)

$w_i^{\left(k\right)}=E\left(\frac{1}{Z_i}/X_i\right)=\frac{{\alpha }^{\left(k\right)}}{{\delta }^{\left(k\right)}\sqrt{{\varphi }^{\left(k\right)}\left(x_i\right)}}\frac{K_1\left({\alpha }^{\left(k\right)}{\delta }^{\left(k\right)}\sqrt{{\varphi }^{\left(k\right)}\left(x_i\right)}\right)}{K_0\left({\alpha }^{\left(k\right)}{\delta }^{\left(k\right)}\sqrt{{\varphi }^{\left(k\right)}\left(x_i\right)}\right)}$ (38)

${\hat{\gamma }}^{\left(k+1\right)}=\sqrt{\frac{n}{{\displaystyle {\sum }_{i=1}^n}{s}_i^{\left(k\right)}-\frac{n^2}{{\displaystyle {\sum }_{i=1}^n}{w}_i^{\left(k\right)}}}}$ (39)

${\hat{\delta }}^{\left(k+1\right)}=\frac{{\gamma }^{\left(k+1\right)}}{\frac{1}{n}{\displaystyle {\sum }_{i=1}^n}{w}_i^{\left(k\right)}}$ (40)

The loglikelihood function of NRIG distribution is given by

$\log L=n\log \gamma -n\log \pi +n\left(\delta \gamma -\beta \mu \right)+\beta {\displaystyle \underset{i=1}{\overset{n}{\sum }}}{x}_i+{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\log K_0\left(\alpha \delta \varphi {\left(x_i\right)}^{\frac{1}{2}}\right)$ (41)

The k-th iteration becomes

$\begin{array}{c}{l}^{\left(k\right)}=n\log {\gamma }^{\left(k\right)}-n\log \pi +n\left({\delta }^{\left(k\right)}{\gamma }^{\left(k\right)}-{\beta }^{\left(k\right)}{\mu }^{\left(k\right)}\right)+{\beta }^{\left(k\right)}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{x}_i\\ +{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\log K_0\left({\alpha }^{\left(k\right)}{\delta }^{\left(k\right)}{\varphi }^{\left(k\right)}{\left(x_i\right)}^{\frac{1}{2}}\right)\end{array}$ (42)

7.6. Special Case When the Index Parameter Is −3/2

$l_2=3n\log \delta -\frac{n}{2}\log \left(2\pi \right)-\frac{5}{2}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\log z_i-n\log \left(1+\delta \gamma \right)+n\delta \gamma -\frac{1}{2}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\left(\frac{{\delta }^2}{z_i}+{\gamma }^2{z}_i\right)$

$\therefore \frac{\partial }{\partial \gamma }{l}_2=0\Rightarrow \hat{\gamma }=\frac{{\delta }^2-\bar{z}}{\delta \bar{z}}$

$\frac{\partial }{\partial \delta }{l}_2=0\Rightarrow \frac{3n}{\delta }-n\left(\frac{n{\delta }^2-{\displaystyle {\sum }_{i=1}^n}{z}_i}{\delta {\displaystyle {\sum }_{i=1}^n}{z}_i}\right){\left(1+\delta \gamma \right)}^{-1}+n\left(\frac{n{\delta }^2-{\displaystyle {\sum }_{i=1}^n}{z}_i}{\delta {\displaystyle {\sum }_{i=1}^n}{z}_i}\right)-\delta {\displaystyle \underset{i=1}{\overset{n}{\sum }}}\frac{1}{z_i}=0$

$\therefore {\delta }^4\left(n^2-{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{z}_i{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\frac{1}{z_i}\right)+n{\delta }^2{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{z}_i+{\left({\displaystyle \underset{i=1}{\overset{n}{\sum }}}{z}_i\right)}^2=0$

$\therefore \left(n^2-{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{z}_i{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\frac{1}{z_i}\right)t^2+\left(n{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{z}_i\right)t+{\left({\displaystyle \underset{i=1}{\overset{n}{\sum }}}{z}_i\right)}^2=0$

i.e., $\left[{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\frac{1}{z_i}-\frac{n}{\bar{z}}\right]t^2-nt-n\bar{z}=0$

where

$t={\delta }^2$

Therefore

$\hat{\delta }={\left[\frac{n{\displaystyle {\sum }_{i=1}^n}{z}_i+{\displaystyle {\sum }_{i=1}^n}{z}_i{\left(4{\displaystyle {\sum }_{i=1}^n}{z}_i{\displaystyle {\sum }_{i=1}^n}\frac{1}{z_i}-3n^2\right)}^{\frac{1}{2}}}{2\left[{\displaystyle {\sum }_{i=1}^n}{z}_i{\displaystyle {\sum }_{i=1}^n}\frac{1}{z_i}-n^2\right]}\right]}^{\frac{1}{2}}$

7.6.1. Posterior Expectation

$E\left(Z/X\right)=\frac{\delta \sqrt{\varphi \left(x\right)}}{\alpha }\frac{K_1\left(\alpha \delta \sqrt{\varphi \left(x\right)}\right)}{K_2\left(\alpha \delta \sqrt{\varphi \left(x\right)}\right)}$ (43)

and

$E\left(\frac{1}{Z}/X\right)=\frac{\alpha }{\delta \sqrt{\varphi \left(x\right)}}\frac{K_3\left(\alpha \delta \sqrt{\varphi \left(x\right)}\right)}{K_2\left(\alpha \delta \sqrt{\varphi \left(x\right)}\right)}$ (44)

7.6.2. Iterations

Let

$s_i=E\left(Z_i/X_i\right)=\frac{\delta \sqrt{\varphi \left(x_i\right)}}{\alpha }\frac{K_1\left(\alpha \delta \sqrt{\varphi \left(x_i\right)}\right)}{K_2\left(\alpha \delta \sqrt{\varphi \left(x_i\right)}\right)}$ (45)

and

$w_i=E\left(\frac{1}{Z_i}/X_i\right)=\frac{\alpha }{\delta \sqrt{\varphi \left(x_i\right)}}\frac{K_3\left(\alpha \delta \sqrt{\varphi \left(x_i\right)}\right)}{K_2\left(\alpha \delta \sqrt{\varphi \left(x_i\right)}\right)}$ (46)

For iterations we have

$s_i^{\left(k\right)}=\frac{{\delta }^{\left(k\right)}\sqrt{{\varphi }^{\left(k\right)}\left(x_i\right)}}{{\alpha }^{\left(k\right)}}\frac{K_1\left({\alpha }^{\left(k\right)}{\delta }^{\left(k\right)}\sqrt{{\varphi }^{\left(k\right)}\left(x_i\right)}\right)}{K_2\left({\alpha }^{\left(k\right)}{\delta }^{\left(k\right)}\sqrt{{\varphi }^{\left(k\right)}\left(x_i\right)}\right)}$ (47)

$w_i^{\left(k\right)}=E\left(\frac{1}{Z_i}/X_i\right)=\frac{{\alpha }^{\left(k\right)}}{{\delta }^{\left(k\right)}\sqrt{{\varphi }^{\left(k\right)}\left(x_i\right)}}\frac{K_3\left({\alpha }^{\left(k\right)}{\delta }^{\left(k\right)}\sqrt{{\varphi }^{\left(k\right)}\left(x_i\right)}\right)}{K_2\left({\alpha }^{\left(k\right)}{\delta }^{\left(k\right)}\sqrt{{\varphi }^{\left(k\right)}\left(x_i\right)}\right)}$ (48)

${\hat{\delta }}^{\left(k+1\right)}={\left[\frac{n{\displaystyle {\sum }_{i=1}^n}{s}_i^{\left(k\right)}+{\displaystyle {\sum }_{i=1}^n}{s}_i^{\left(k\right)}{\left(4{\displaystyle {\sum }_{i=1}^n}{s}_i^{\left(k\right)}{\displaystyle {\sum }_{i=1}^n}{w}_i^{\left(k\right)}-3n^2\right)}^{\frac{1}{2}}}{2\left[{\displaystyle {\sum }_{i=1}^n}{s}_i^{\left(k\right)}{\displaystyle {\sum }_{i=1}^n}{w}_i^{\left(k\right)}-n^2\right]}\right]}^{\frac{1}{2}}$ (49)

${\hat{\gamma }}^{\left(k+1\right)}=\frac{{\delta }^{\left(k+1\right)}}{\frac{1}{n}{\displaystyle {\sum }_{i=1}^n}{s}_i^{\left(k\right)}}-\frac{1}{{\delta }^{\left(k+1\right)}}$ (50)

The log-likelihood function of Normal- $GIG\left(-\frac{3}{2}\mathrm{,}\delta \mathrm{,}\gamma \right)$ is given by

$\begin{array}{c}\log L=-n\log \pi -n\log \left(1+\delta \gamma \right)+2n\log \alpha +n\log \delta +n\left(\delta \gamma -\beta \mu \right)\\ +\beta {\displaystyle \underset{i=1}{\overset{n}{\sum }}}{x}_i-{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\varphi \left(x_i\right)+{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\log K_2\left(\alpha \delta \varphi {\left(x_i\right)}^{\frac{1}{2}}\right)\end{array}$ (51)

The k-th iteration becomes

$\begin{array}{c}{l}^{\left(k\right)}=-n\log \pi -n\log \left(1+{\delta }^{\left(k\right)}{\gamma }^{\left(k\right)}\right)+2n\log {\alpha }^{\left(k\right)}+n\log {\delta }^{\left(k\right)}+n({\delta }^{\left(k\right)}{\gamma }^{\left(k\right)}\\ -{\beta }^{\left(k\right)}{\mu }^{\left(k\right)})+{\beta }^{\left(k\right)}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{x}_i-{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\varphi }^{\left(k\right)}\left(x_i\right)+{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\log K_2\left({\alpha }^{\left(k\right)}{\delta }^{\left(k\right)}{\varphi }^{\left(k\right)}{\left(x_i\right)}^{\frac{1}{2}}\right)\end{array}$ (52)

7.7. Special Case When the Index Parameter Is 3/2

7.7.1. M-Step

$\begin{array}{c}{l}_2={\displaystyle \underset{i=1}{\overset{n}{\sum }}}\log g_4\left(z_i\right)\\ =-\frac{n}{2}-n\log \left(1+\delta \gamma \right)+3n\log \gamma +n\delta \gamma +\frac{1}{2}\log {\displaystyle \underset{i=1}{\overset{n}{\sum }}}\left(z_i\right)\\ -\frac{{\delta }^2}{2}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\frac{1}{z_i}-\frac{{\gamma }^2}{2}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{z}_i\end{array}$

$\therefore \frac{\partial }{\partial \delta }{l}_2=0\Rightarrow \frac{n{\gamma }^2\delta }{1+\delta \gamma }-\delta {\displaystyle \underset{i=1}{\overset{n}{\sum }}}\frac{1}{z_i}=0$

$\therefore \hat{\delta }=\frac{n{\hat{\gamma }}^2-{\displaystyle {\sum }_{i=1}^n}\frac{1}{z_i}}{\hat{\gamma }{\displaystyle {\sum }_{i=1}^n}\frac{1}{z_i}}$

Similarly

$\frac{\partial }{\partial \gamma }{l}_2=0\Rightarrow \frac{{\delta }^2\gamma }{1+\delta \gamma }+\frac{3}{\gamma }-\gamma \bar{z}=0$

$\therefore {\gamma }^2{\delta }^2+\left(3\gamma -{\gamma }^3\bar{z}\right)\delta -{\gamma }^2\bar{z}+3=0$

Substituting for $\delta$ and letting ${\gamma }^2=t$ we obtain

$\left(n^2-{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{z}_i{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\frac{1}{z_i}\right)t^2+\left(n{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\frac{1}{z_i}\right)t+{\left({\displaystyle \underset{i=1}{\overset{n}{\sum }}}\frac{1}{z_i}\right)}^2=0$

Therefore

$t=\frac{-b\pm \sqrt{b^2-4ac}}{2a}$

where

$a=n^2-{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{z}_i{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\frac{1}{z_i}$

$b=n{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\frac{1}{z_i}$

$c={\left({\displaystyle \underset{i=1}{\overset{n}{\sum }}}\frac{1}{z_i}\right)}^2$

$\therefore \hat{\gamma }={\left[\frac{-b+\sqrt{b^2-4ac}}{2a}\right]}^{\frac{1}{2}}$

7.7.2. E-Step

Posterior Expectation

$E\left(Z/X\right)=\frac{\delta \sqrt{\varphi \left(x\right)}}{\alpha }\frac{K_2\left(\alpha \delta \sqrt{\varphi \left(x\right)}\right)}{K_1\left(\alpha \delta \sqrt{\varphi \left(x\right)}\right)}$ (53)

$E\left(\frac{1}{Z}/X\right)=\frac{\alpha K_0\left(\alpha \delta \sqrt{\varphi \left(x\right)}\right)}{\delta \sqrt{\varphi \left(x\right)}{K}_1\left(\alpha \delta \sqrt{\varphi \left(x\right)}\right)}$ (54)

7.8. Iterations

Let

$s_i=E\left(Z_i/X_i\right)=\frac{\delta \sqrt{\varphi \left(x_i\right)}}{\alpha }\frac{K_2\left(\alpha \delta \sqrt{\varphi \left(x_i\right)}\right)}{K_1\left(\alpha \delta \sqrt{\varphi \left(x_i\right)}\right)}$

and

$w_i=E\left[\frac{1}{Z_i}/X_i\right]=\frac{\alpha K_0\left(\alpha \delta \sqrt{\varphi \left(x_i\right)}\right)}{\delta \sqrt{\varphi \left(x\right)}{K}_1\left(\alpha \delta \sqrt{\varphi (\; x\; )}\right)}$

For iterations we have

$s_i^{\left(k\right)}=\frac{{\delta }^{\left(k\right)}}{{\alpha }^{\left(k\right)}}\sqrt{{\varphi }^{\left(k\right)}\left(x_i\right)}\frac{K_2\left({\alpha }^{\left(k\right)}{\delta }^{\left(k\right)}\sqrt{{\varphi }^{\left(k\right)}\left(x_i\right)}\right)}{K_1\left({\alpha }^{\left(k\right)}{\delta }^{\left(k\right)}\sqrt{{\varphi }^{\left(k\right)}\left(x_i\right)}\right)}$ (55)

$w_i^{\left(k\right)}=\frac{{\alpha }^{\left(k\right)}{K}_0\left({\alpha }^{\left(k\right)}{\delta }^{\left(k\right)}\sqrt{{\varphi }^{\left(k\right)}\left(x_i\right)}\right)}{{\delta }^{\left(k\right)}\sqrt{{\varphi }^{\left(k\right)}\left(x_i\right)}{K}_1\left({\alpha }^{\left(k\right)}{\delta }^{\left(k\right)}\sqrt{{\varphi }^{\left(k\right)}\left(x_i\right)}\right)}$ (56)

${\hat{\gamma }}^{\left(k+1\right)}=\sqrt{t^{\left(k\right)}}$ (57)

${\hat{\delta }}^{\left(k+1\right)}=\frac{n{\left({\hat{\gamma }}^{\left(k+1\right)}\right)}^2-{\displaystyle {\sum }_{i=1}^n}{w}_i^{\left(k\right)}}{{\hat{\gamma }}^{\left(k+1\right)}{\displaystyle {\sum }_{i=1}^n}{w}_i^{\left(k\right)}}$ (58)

The log-likelihood function for the Normal- $GIG\left(\frac{3}{2}\mathrm{,}\delta \mathrm{,}\gamma \right)$ distribution is given by

$\begin{array}{l}l=3n\log \gamma +n\left(\delta \gamma -\beta \mu \right)-n\log \left(\alpha \pi \left(1+\delta \gamma \right)\right)+n\log \delta \\ +\beta {\displaystyle \underset{i=1}{\overset{n}{\sum }}}{x}_i+\frac{1}{2}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\log \varphi \left(x_i\right)+{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\log K_1\left(\alpha \delta \sqrt{\varphi \left(x\right)}\right)\end{array}$ (59)

For the k-th iteration we have

$\begin{array}{c}{l}^{\left(k\right)}=3n\log {\gamma }^{\left(k\right)}+n\left({\delta }^{\left(k\right)}{\gamma }^{\left(k\right)}-{\beta }^{\left(k\right)}{\mu }^{\left(k\right)}\right)-n\log \left({\alpha }^{\left(k\right)}\pi \left(1+{\delta }^{\left(k\right)}{\gamma }^{\left(k\right)}\right)\right)\\ +n\log {\delta }^{\left(k\right)}+{\beta }^{\left(k\right)}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{x}_i+\frac{1}{2}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\log {\varphi }^{\left(k\right)}\left(x_i\right)+{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\log K_1\left({\alpha }^{\left(k\right)}{\delta }^{\left(k\right)}\sqrt{{\varphi }^{\left(k\right)}\left(x\right)}\right)\end{array}$ (60)

Remarks:

1) For all the proposed models, the $\beta$, $\mu$ and $\alpha$ parameters of the conditional distribution are updated as follows:

${\hat{\beta }}^{\left(k+1\right)}=\frac{{\displaystyle {\sum }_{i=1}^n}\left(x_i-\bar{x}\right)w_i^{\left(k\right)}}{n-\frac{1}{n}{\displaystyle {\sum }_{i=1}^n}{s}_i^{\left(k\right)}{\displaystyle {\sum }_{i=1}^n}{w}_i^{(\; k\; )}}$

${\hat{\mu }}^{\left(k+1\right)}=\bar{x}-{\hat{\beta }}^{\left(k+1\right)}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\frac{s_i^{\left(k\right)}}{n}$

${\hat{\alpha }}^{\left(k+1\right)}={\left[{\left({\gamma }^{\left(k+1\right)}\right)}^2+{\left({\beta }^{\left(k+1\right)}\right)}^2\right]}^{\frac{1}{2}}$

2) The stopping criterion is when

$\frac{l^{\left(k\right)}-l^{\left(k-1\right)}}{l^{\left(k\right)}}<tol$

where tol is the tolerance level chosen; e.g 10⁻⁶.

3) Initial values used are moment estimates of NIG ditsribution as proposed by Karlis [23].

8. Application

8.1. Fitting of the Proposed Models

The data used in this research is the Range Resource Corporation weekly returns for the period 3/01/2000 to 1/07/2013 with 702 observations. The histogram for the weekly log-returns in Figure 1 shows that the data is negatively skewed and exhibiting heavy tails.The Q-Q plot shows that the normal distribution is not a good fit for the data especially at the tails.

Table 5 provides descriptive statistics for the return series in consideration. We observe that the excess kurtosis of 2.768252 indicates the leptokurtic behaviour of the returns. The log-returns has a distributions with relatively heavier tails than the normal distribution. We observe skewness of −0.1886714 which indicates that the two tails of the returns behave slightly differently.

The proposed models are now fitted to RRC weekly log-returns. Using the sample estimates and the NIG estimators to the RRC data we obtain the following estimates as initial values for the EM algorithm (see Karlis [23] ).

$\hat{\alpha }=0.3722511,\hat{\beta }=-0.02456226,\hat{\delta }=2.950864,\hat{\mu }=0.4284473$.

Figure 1. Fitting Model 1 to RRC weekly returns.

Table 5. Summary Statistics for RRC weekly log-returns.

The initial values were used in all the proposed models to obtain the maximum likelihood estimates as shown in Table 6 below

The parameter estimates from Table 2 are now fitted to RRC weekly log-returns. Figures 2-5 show the histogram and Q-Q plots of the RRC returns fitted with the proposed models. Figure 2-5 show that the proposed model fit the data well.

Table 7 present results of Kolmogorv-Smirnov and Anderson-Darling test performed on the models. All models produce high p-values, a strong evidence that we can not reject the null hypothesis that the returns data follow the proposed models.

Table 8 presents values of Alkaike Information Criterion (AIC), Bayesian Information Creterion (BIC) and Log-likelihood. The values illustrate that the models are alternative to each other.

Model 3 has the lowest AIC and BIC with the highest log-likelihood. It is the best fit for the data.

Figure 2. Fitting Model 1 to RRC weekly returns.

Figure 3. Fitting Model 2 to RRC weekly returns.

Figure 4. Fitting Model 3 to RRC weekly returns.

Figure 5. Fitting Model 4 to RRC weekly returns.

Table 6. Estimates for the proposed models.

Table 7. Results from Kolmogorov-Smirnov Test and Anderson-Darling Test.

Table 8. AIC, BIC and Log-likelihood Values.

8.2. Risk Estimation and Backtesting

We use the parameter estimates for our proposed model to determine the VaR (Table 9) and ES (Table 10) at levels $\alpha \in \left\{0.001,0.01,0.05,0.95,0.99,0.999\right\}$. The first three level are used to measure the risk of long position, while the last three levels are used to measure the risk of short positions. We apply the Kupiec Likelihood Ratio (LR) test given by Kupiec [16] which test the hypothesis that the expected proposition of violations is equal to $\alpha$. The method consist of calculating $\tau \left(\alpha \right)$ the number of times the observed returns, $x_t$ falls below (for long position) or above (for short position) the $Va{R}_{\alpha }$ estimates at level $\alpha$ ; i.e., $x_t<Va{R}_{\alpha }$ or $x_t<Va{R}_{1-\alpha }$, and compare the corresponding failure rate to $\alpha$.

The likelihood ratio statistic is given by

$2\log {\left(\frac{\tau \left(\alpha \right)}{n}\right)}^{\tau \left(\alpha \right)}{\left(1-\frac{\tau \left(\alpha \right)}{n}\right)}^{n-\tau \left(\alpha \right)}-2\log \left({\alpha }^{\tau \left(\alpha \right)}-{\left(1-\alpha \right)}^{n-\tau \left(\alpha \right)}\right)$ (61)

where $\tau \left(\alpha \right)$ is the number of violations. Under the null hypothesis this statistic is distributed as ${\chi }^2$ distribution with one degree of freedom. Table 11 gives the number of violations for the models. Table 12 gives the p-values based on Kupiec test.

Model 3 has the highest VaR and ES value indicating that it perform well than the other models at the tails.

Table 9. VaR Values of RRC log-returns based on normal and proposed models.

Table 10. ES Values of RRC log-returns based on normal and proposed models.

Table 11. Number of violations of VaR for each distribution at different levels.

Table 12. P-value for the kupiec test for each distribution at different levels.

Remark: At 5 percent level of significant, the Normal distribution is rejected at levels: 0.001, 0.95 and 0.999. In addition it is also rejected at level 0.99 at 10 percent level of significant. The Normal weighted Inverse Gaussian distributions were all effective and well specified on all levels of VaR. It can be noted model 3 (Normal- $GIG\left(-\frac{3}{2}\mathrm{,}\delta \mathrm{,}\gamma \right)$ ) outperforms the other models at level 0.99.

9. Conclusions

In this paper we obtained VaR using NWIG distributions. We first constructed Normal Variance-Mean Mixture when the mixing distributions are $GIG\left(-\frac{1}{2}\mathrm{,}\delta \mathrm{,}\gamma \right)$, $GIG\left(\frac{1}{2}\mathrm{,}\delta \mathrm{,}\gamma \right)$, $GIG\left(-\frac{3}{2}\mathrm{,}\delta \mathrm{,}\gamma \right)$ and $GIG\left(\frac{3}{2}\mathrm{,}\delta \mathrm{,}\gamma \right)$. We have shown that these mixing distribution are WIG distributions.

The parameters of the mixed models are estimated using EM-algorithm. The iterative schemes used are based on explicit solutions of normal equations. We used method of moment estimates of NIG as initial values and obtained monotonic convergence.

We used AIC, BIC and loglikelihood for model selection. Normal- $GIG\left(-\frac{3}{2}\mathrm{,}\delta \mathrm{,}\gamma \right)$ was found to be the best model. The results show that the three NWIG distributions are as good as NIG for VaR computation.

Further work can be done on Normal Mixtures when the mixing distributions are Finite mixtures of WIG distributions.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Zhang, Y. and Nadarajah, S. (2017) A Review of Backtesting for Values at Risk. In: Communications in Statistics—Theory and Methods, 47, 3616-3639.
[2]	Wimmerstedt, L. (2015) Backtesting Expected Shortfall: The Design and Implementation of Different Backtest. Master’s Thesis in Mathematical Statistics, Royal Institute of Technology, Stockholm.
[3]	Artzner, P., Delbaen, F., Eber, J.-M. and Heath, D. (1997) Thinking Coherently. Risk Journal, 10, 68-71.
[4]	Artzner, P., Delbaen, F., Eber, J.-M. and Heath, D. (1997) Coherent Measures of Risk. Mathematical Finance, 9, 203-228. https://doi.org/10.1111/1467-9965.00068
[5]	Gneiting, T. (2011) Making and Evaluating Point Forecasts. Journal of the American Statistical Association, 106, 746-762. https://doi.org/10.1198/jasa.2011.r10138
[6]	Basel Committee (2013) Fundamental Review of the Trading Book: A Revised Market Risk Framework.
[7]	Nadarajah, S., Zhang, B. and Chan, S. (2014) Estimation Methods for Expected Shortfall. Quantitative Finance, 14, 271-291. https://doi.org/10.1080/14697688.2013.816767
[8]	Barndorff-Nielsen, O.E. (1977) Exponentially Decreasing Distributions for the Logarithm of Particle Size. Proceedings of the Royal Society of London. Series A, 353, 409-419. https://doi.org/10.1098/rspa.1977.0041
[9]	Barndorff-Nielsen, O.E. (1997) Normal Inverse Gaussian Distribution and Stochastic Volatility Modelling. Scandinavian Journal of Statistics, 24, 1-13. https://doi.org/10.1111/1467-9469.00045
[10]	Dempster, A.P., Laird, N.M. and Rubin, D. (1977) Maximum Likelihood from Incomplete Data Using the EM Algorithm. Journal of the Royal Statistical Society: Series B, 39, 1-38. https://doi.org/10.1111/j.2517-6161.1977.tb01600.x
[11]	Fisher, R.A. (1934) The Effect of Methods of Ascertainment upon the Estimation of Frequencies. The Annals of Human Genetics, 6, 13-25. https://doi.org/10.1111/j.1469-1809.1934.tb02105.x
[12]	Patil, G.P. and Rao, C.R. (1978) Weighted Distributions and Size-Biased Sampling with Applications to Wildlife Populations and Human Families. Biometrics, 34, 179-189. https://doi.org/10.2307/2530008
[13]	Akman, H.O. and Gupta, R.C. (1992) A Comparison of Various Estimators of the Mean of an Inverse Gaussian Distribution. Journal of Statistical Computation and Simulation, 40, 71-81. https://doi.org/10.1080/00949659208811366
[14]	Gupta, R.C. and Akman, O. (1995) On Reliability Studies of a Weighted Inverse Gaussian Distribution Model. Statistical Papers, 38, 445-452. https://doi.org/10.1007/BF02925999
[15]	Gupta, R.C. and Kundu, D. (2011) Weighted Inverse Gaussian—A Versatile Lifetime Model. Journal of Applied Statistics, 38, 2695-2708. https://doi.org/10.1080/02664763.2011.567251
[16]	Kupiec, P. (1995) Techniques for Verifying the Accuracy of Risk Measurement Models. https://doi.org/10.3905/jod.1995.407942
[17]	Rappoport, P. (1993) A New Approach: Average Shortfall. Technical Report, J.P. Morgan.
[18]	Yamai, Y. and Yoshiba, T. (2002) On the Validity of Value at Risk Comparative Analyses with Expected Shortfall. Monetary and Economic Studies.
[19]	Zhang, Y., Chu, J., Chan, S. and Chan, B. (2019) The Generalised Hyperbolic Distribution and Its Subclass in Analysis of a New Era of Cryptocurrencies: Ethereum and Its Financial Risk. Physica A, 526, 1-12. https://doi.org/10.1016/j.physa.2019.04.136
[20]	Bams, D. and Wilhouwer, J.L. (2001) Empirical Issues in Value at Risk. ASTIN Bulletin, 31, 299-315. https://doi.org/10.2143/AST.31.2.1007
[21]	Abramowitz, M and Stegun, I.A. (1972) Handbook of Mathematical Function. Dover, New York.
[22]	Kostas, F. (2007) Tests of Fit for Symmetric Variance Gamma Distributions. UADPhilEcon, National and Kapodistrian University of Athens, Athens.
[23]	Karlis, D. (2002) An EM Type Algorithm for Maximum Likelihood Estimation of the Normal-Inverse Gaussian Distribution. Statistics and Probability Letters, 57, 43-52. https://doi.org/10.1016/S0167-7152(02)00040-8