The Generalized Hyperbolic Distribution

The Generalized Hyperbolic Distribution#

In this section we briefly review the basic properties of the Generalized Hyperbolic (GH) distribution.

Definition as Normal Mixture#

Let \(Y\) be a GIG random variable with parameters \((p, a, b)\), and \(Z\) be an independent Gaussian random vector with zero mean and covariance \(\Sigma\). Then the random vector:

(1)#\[X \stackrel{d}{=} \mu + \gamma Y + \sqrt{Y} Z\]

follows the generalized hyperbolic distribution with parameters \((\mu, \gamma, \Sigma, p, a, b)\), where \(\mu, \gamma \in \mathbb{R}^d\) and \(\Sigma \in \mathbb{R}^{d \times d}\) is a positive definite matrix.

Parameter Interpretation#

  • \(\mu\): location parameter

  • \(\gamma\): skewness parameter

  • \(\Sigma\): models the dependency structure of the multivariate distribution

  • \((p, a, b)\): GIG parameters controlling the heavy-tailedness

In general, many multivariate heavy-tailed distributions can be defined by (1) given some non-negative random variable \(Y\). These distributions are usually called normal mixture or Gaussian mixture distributions.

Note

In many references, “normal mixture” refers to a discrete mixture of normal densities. In normix, we use “normal mixture” to refer to any random variable that can be expressed by (1).

Joint GH Distribution#

The joint distribution of \(X\) and \(Y\) is crucial in analyzing the GH distribution. We call the distribution of \((X, Y)\) the joint-GH distribution. Its density function is:

(2)#\[\begin{split}f(x, y | \mu, \gamma, \Sigma, p, a, b) &= \frac{1}{\sqrt{(2\pi)^d |\Sigma|}} \frac{(a/b)^{p/2}}{2 K_p(\sqrt{ab})} y^{p - 1 - d/2} \\ &\quad \times \exp\left( -\frac{1}{2}(x - \mu - \gamma y)^\top \Sigma^{-1} (x - \mu - \gamma y) y^{-1} - \frac{1}{2}(b y^{-1} + a y) \right),\end{split}\]

for \(y > 0\).

Marginal GH Density#

Integrating out \(y\) from (2), we obtain the marginal distribution of \(x\), which is the GH density function:

(3)#\[\begin{split}f(x | \mu, \gamma, \Sigma, p, a, b) &= \int_0^\infty f(x, y | \mu, \gamma, \Sigma, p, a, b) \, dy \\ &= c \frac{K_{p - d/2}\left(\sqrt{(b + q(x))(a + \gamma^\top \Sigma^{-1} \gamma)}\right)} {\left(\sqrt{(b + q(x))(a + \gamma^\top \Sigma^{-1} \gamma)}\right)^{d/2 - p}} e^{(x - \mu)^\top \Sigma^{-1} \gamma},\end{split}\]

where \(q(x) = (x - \mu)^\top \Sigma^{-1} (x - \mu)\) is the Mahalanobis distance, and

\[c = \frac{(a/b)^{p/2} (a + \gamma^\top \Sigma^{-1} \gamma)^{d/2 - p}} {(2\pi)^{d/2} |\Sigma|^{1/2} K_p(\sqrt{ab})}.\]

Alternative Parameterization#

Similar to the GIG distribution, one can use another parameterization with \(\delta = \sqrt{b/a}\) and \(\eta = \sqrt{ab}\):

\[f(x | \mu, \gamma, \Sigma, p, \delta, \eta) = c \frac{K_{p - d/2}\left(\sqrt{(\eta + q_\delta(x))(\eta + \tilde{q}_\delta)}\right)} {\left(\sqrt{(\eta + q_\delta(x))(\eta + \tilde{q}_\delta)}\right)^{d/2 - p}} e^{(x - \mu)^\top (\delta \Sigma)^{-1} \delta \gamma},\]

where \(q_\delta(x) = (x - \mu)^\top (\delta \Sigma)^{-1} (x - \mu)\), \(\tilde{q}_\delta = (\delta \gamma)^\top (\delta \Sigma)^{-1} (\delta \gamma)\), and

\[c = \frac{(\eta + \tilde{q}_\delta)^{d/2 - p}}{(2\pi)^{d/2} |\delta \Sigma|^{1/2} K_p(\eta)}.\]

Model Identifiability#

From the above representation, one can observe that the GH model is not regular since the parameter sets \((\mu, \gamma/c, \Sigma/c, p, c\delta, \eta)\) give the same distribution for any \(c > 0\). Therefore, the Fisher information matrix of the GH distribution is singular.

There are several ways to regularize the GH family:

  1. Set \(\delta = 1\) (simplest approach)

  2. Set \(b = 1\) in the EM algorithm [Protassov2004]

  3. Fix \(b\) when \(p > -1\) and fix \(a\) when \(p < 1\) [Hu2005]

  4. Fix the determinant \(|\Sigma| = 1\) [McNeil2010]

Note

When the dimension \(d\) is high, fixing \(|\Sigma| = 1\) is recommended. Since \(|\Sigma/c| = |\Sigma|/c^d\), any small perturbation of the matrix scale will make \(|\Sigma|\) change dramatically when \(d\) is large. The matrix inversion would be intractable if \(|\Sigma|\) is too large or too small.

Exponential Family Form#

The joint-GH distribution (2) belongs to the exponential family with density:

\[f(x, y|\theta) = h(x, y) \exp\left(\theta^\top t(x, y) - \psi(\theta)\right)\]

Sufficient Statistics:

\[\begin{split}t(x, y) = \begin{pmatrix} \log y \\ y^{-1} \\ y \\ x \\ x y^{-1} \\ x x^\top y^{-1} \end{pmatrix}\end{split}\]

where \(t_1, t_2, t_3 \in \mathbb{R}\), \(t_4, t_5 \in \mathbb{R}^d\), and \(t_6 \in \mathbb{R}^{d \times d}\).

Natural Parameters:

The natural parameters are derived from the classical parameters \((\mu, \gamma, \Sigma, p, a, b)\). By expanding the exponent in (2):

\[\begin{split}&-\frac{1}{2}(x - \mu - \gamma y)^\top \Sigma^{-1} (x - \mu - \gamma y) y^{-1} - \frac{1}{2}(b y^{-1} + a y) \\ &= -\frac{1}{2} x^\top \Sigma^{-1} x \, y^{-1} + x^\top \Sigma^{-1} \mu \, y^{-1} + x^\top \Sigma^{-1} \gamma \\ &\quad - \frac{1}{2} \mu^\top \Sigma^{-1} \mu \, y^{-1} - \mu^\top \Sigma^{-1} \gamma - \frac{1}{2} \gamma^\top \Sigma^{-1} \gamma \, y \\ &\quad - \frac{1}{2} b \, y^{-1} - \frac{1}{2} a \, y\end{split}\]

we identify the natural parameters:

(4)#\[\begin{split}\theta_1 &= p - 1 - \frac{d}{2} \quad \text{(coefficient of } \log y \text{)} \\ \theta_2 &= -\frac{1}{2}\left(b + \mu^\top \Sigma^{-1} \mu\right) \quad \text{(coefficient of } y^{-1} \text{)} \\ \theta_3 &= -\frac{1}{2}\left(a + \gamma^\top \Sigma^{-1} \gamma\right) \quad \text{(coefficient of } y \text{)} \\ \theta_4 &= \Sigma^{-1} \gamma \quad \text{(coefficient of } x \text{)} \\ \theta_5 &= \Sigma^{-1} \mu \quad \text{(coefficient of } x y^{-1} \text{)} \\ \theta_6 &= -\frac{1}{2} \Sigma^{-1} \quad \text{(coefficient of } x x^\top y^{-1} \text{)}\end{split}\]

Base Measure:

\[h(x, y) = \frac{1}{(2\pi)^{d/2}} \mathbf{1}_{y > 0}\]

Log Partition Function:

(5)#\[\psi(\theta) = \frac{1}{2} \log|\Sigma| + \log 2 + \log K_p(\sqrt{ab}) + \frac{p}{2} \log\left(\frac{b}{a}\right) + \mu^\top \Sigma^{-1} \gamma\]

Expectation Parameters:

The expectation parameters \(\eta = \nabla\psi(\theta) = E[t(X, Y)]\) are:

(6)#\[\begin{split}\eta_1 &:= E[\log Y] = \left.\frac{\partial}{\partial \alpha} \left(\sqrt{\frac{b}{a}}\right)^\alpha \frac{K_{p+\alpha}(\sqrt{ab})}{K_p(\sqrt{ab})} \right|_{\alpha=0} \\ \eta_2 &:= E[Y^{-1}] = \sqrt{\frac{a}{b}} \frac{K_{p-1}(\sqrt{ab})}{K_p(\sqrt{ab})} \\ \eta_3 &:= E[Y] = \sqrt{\frac{b}{a}} \frac{K_{p+1}(\sqrt{ab})}{K_p(\sqrt{ab})} \\ \eta_4 &:= E[X] = \mu + \gamma \eta_3 \\ \eta_5 &:= E[X Y^{-1}] = \mu \eta_2 + \gamma \\ \eta_6 &:= E[X X^\top Y^{-1}] = \Sigma + \mu \mu^\top \eta_2 + \gamma \gamma^\top \eta_3 + \mu \gamma^\top + \gamma \mu^\top\end{split}\]

where \(\eta_1, \eta_2, \eta_3 \in \mathbb{R}\), \(\eta_4, \eta_5 \in \mathbb{R}^d\), and \(\eta_6 \in \mathbb{R}^{d \times d}\).

Note that \(\eta_1, \eta_2, \eta_3\) are exactly the expectation parameters of the GIG random variable \(Y\).

Recovering Parameters from Expectations#

Given all the expectation parameters, we can recover the original parameters as follows:

(7)#\[\begin{split}\mu &= \frac{\eta_4 - \eta_3 \eta_5}{1 - \eta_2 \eta_3}, \\ \gamma &= \frac{\eta_5 - \eta_2 \eta_4}{1 - \eta_2 \eta_3}, \\ \Sigma &= \eta_6 - \eta_5 \mu^\top - \mu \eta_5^\top + \eta_2 \mu \mu^\top - \eta_3 \gamma \gamma^\top, \\ (p, a, b) &= \arg\max_{p,a,b} L_{\text{GIG}}(p, a, b | \eta_1, \eta_2, \eta_3),\end{split}\]

where \(L_{\text{GIG}}\) is the GIG log-likelihood function given in (8).

These equations form the M-step in the EM algorithm, where the expectations in (6) are replaced by conditional expectations.

Hellinger Distance#

While there is no analytical formulation of the Hellinger distance between two GH distributions, the Hellinger distance of the joint-GH distributions is tractable.

Proposition. Let \(\theta_1 = (\mu_1, \gamma_1, \Sigma_1, p_1, a_1, b_1)\) and \(\theta_2 = (\mu_2, \gamma_2, \Sigma_2, p_2, a_2, b_2)\) be the parameters of two joint-GH distributions. The squared Hellinger distance between the two distributions is:

\[H_{\text{JGH}}^2(\theta_1 \| \theta_2) = 1 - \frac{|\Sigma_1 \Sigma_2|^{1/4}}{|\bar{\Sigma}|^{1/2}} \frac{(a_1/b_1)^{p_1/4} (a_2/b_2)^{p_2/4}} {\sqrt{K_{p_1}(\sqrt{a_1 b_1}) K_{p_2}(\sqrt{a_2 b_2})}} \frac{K_{\bar{p}}(\sqrt{\bar{b}' \bar{a}'})}{(\bar{a}'/\bar{b}')^{\bar{p}/2}} e^{-\frac{1}{4} \Delta\mu^\top \bar{\Sigma}^{-1} \Delta\gamma},\]

where:

  • \(\Delta\mu = \mu_1 - \mu_2\), \(\Delta\gamma = \gamma_1 - \gamma_2\)

  • \(\bar{\Sigma} = (\Sigma_1 + \Sigma_2)/2\)

  • \(\bar{p} = (p_1 + p_2)/2\), \(\bar{a} = (a_1 + a_2)/2\), \(\bar{b} = (b_1 + b_2)/2\)

  • \(\bar{b}' = \bar{b} + \frac{1}{4} \Delta\mu^\top \bar{\Sigma}^{-1} \Delta\mu\)

  • \(\bar{a}' = \bar{a} + \frac{1}{4} \Delta\gamma^\top \bar{\Sigma}^{-1} \Delta\gamma\)

If \(\mu_1 = \mu_2\), \(\gamma_1 = \gamma_2\), and \(\Sigma_1 = \Sigma_2\), then \(H_{\text{JGH}}(\theta_1 \| \theta_2) = H_{\text{GIG}}(p_1, a_1, b_1 \| p_2, a_2, b_2)\).

Although \(H_{\text{JGH}}\) differs from the Hellinger distance of the marginal GH, it provides an upper bound, so we can use it to measure how close two GH distributions are.

Numerical Stability#

The computation of (7) is not stable in terms of relative error under certain conditions. However, it is relatively stable in terms of the Hellinger distance.

Numerical experiments show that:

  • The relative errors of \(\mu\) and \(\gamma\) are around machine epsilon

  • The relative error of some GIG parameters (especially when \(|p|\) is large) can be large

  • However, the Hellinger distance between true and estimated parameters remains small

This behavior is consistent with the ill-conditioning of the GIG optimization problem discussed in the The Generalized Inverse Gaussian Distribution section.

Special Cases#

Several important distributions are special cases of the GH family:

  • Normal-Inverse Gaussian (NIG): \(p = -1/2\)

  • Variance-Gamma (VG): \(b \to 0\) (Gamma mixing)

  • Normal-Inverse Gamma (NInvG): \(a \to 0\) (Inverse-Gamma mixing)

  • Student-t: \(p < 0\), \(a = 0\), \(\gamma = 0\)

  • Hyperbolic: \(p = 1\)

These are implemented as separate classes in normix:

  • NormalInverseGaussian

  • VarianceGamma

  • NormalInverseGamma

References#

[Protassov2004]

Protassov, R. S. (2004). EM-based maximum likelihood parameter estimation for multivariate generalized hyperbolic distributions.

[Hu2005]

Hu, W. (2005). Calibration of multivariate generalized hyperbolic distributions using the EM algorithm.

[McNeil2010]

McNeil, A. J., Frey, R., & Embrechts, P. (2010). Quantitative Risk Management. Princeton University Press.