Next: 1.4 Value at Risk, Up: 1. Computationally Intensive Value Previous: 1.2 Stable Distributions

Subsections

# 1.3 Hyperbolic Distributions

In response to remarkable regularities discovered by geomorphologists in the 1940s, [4] introduced the hyperbolic law for modeling the grain size distribution of windblown sand. Excellent fits were also obtained for the log-size distribution of diamonds from a large mining area in South West Africa. Almost twenty years later the hyperbolic law was found to provide a very good model for the distributions of daily stock returns from a number of leading German enterprises ([31,57]), giving way to its today's use in stock price modeling ([9]) and market risk measurement ([32]). The name of the distribution is derived from the fact that its log-density forms a hyperbola, see Fig. 1.8. Recall that the log-density of the normal distribution is a parabola. Hence the hyperbolic distribution provides the possibility of modeling heavier tails.

The hyperbolic distribution is defined as a normal variance-mean mixture where the mixing distribution is the generalized inverse Gaussian (GIG) law with parameter , i.e. it is conditionally Gaussian, see [4] and [6]. More precisely, a random variable has the hyperbolic distribution if:

 N (1.19)

where is a generalized inverse Gaussian GIG random variable and N denotes the Gaussian distribution with mean  and variance . The GIG law is a very versatile positive domain three parameter distribution. It arises in the context of the first passage time of a diffusion process, when the drift and variance of displacement per unit time are dependent upon the current position of the particle. The probability density function of a GIG variable is given by:

 (1.20)

with the parameters taking values in one of the ranges: (1)  if , (2)  if or (3)  if . The generalized inverse Gaussian law has a number of interesting properties that we will use later in this section. The distribution of the inverse of a GIG variable is again GIG but with a different , namely if:

 GIG   thenGIG (1.21)

A GIG variable can be also reparameterized by setting and , and defining , where:

 GIG (1.22)

The normalizing constant K in formula (1.20) is the modified Bessel function of the third kind with index , also known as the MacDonald function. It is defined as:

 K (1.23)

In the context of hyperbolic distributions, the Bessel functions are thoroughly discussed in [6]. Here we recall only two properties that will be used later. Namely, (1)  K is symmetric with respect to , i.e. K   K, and (2) for it can be written in a simpler form:

 K (1.24)

For other values of numerical approximations of the integral in (1.23) have to be used, see e.g. [19], [86] or [96].

Relation (1.19) implies that a hyperbolic random variable    H can be represented in the form:

 N

with the characteristic function:

 (1.25)

Here denotes the distribution function of a generalized inverse Gaussian random variable with parameter , see (1.20). Hence, the hyperbolic PDF is given by:

 (1.26)

Sometimes another parameterization of the hyperbolic distribution with and is used. Then the probability density function of the hyperbolic H law can be written as:

 (1.27)

where is the scale parameter, is the location parameter and . The latter two parameters -  and  - determine the shape, with  being responsible for the steepness and  for the skewness. In XploRe the hyperbolic density and distribution functions are implemented in the pdfhyp and cdfhyp quantlets, respectively. The calculation of the PDF is straightforward, however, the CDF has to be numerically integrated from (1.27).

The hyperbolic law is a member of a more general class of generalized hyperbolic distributions. The generalized hyperbolic law can be represented as a normal variance-mean mixture where the mixing distribution is the generalized inverse Gaussian (GIG) law with any . Hence, the generalized hyperbolic distribution is described by five parameters . Its probability density function is given by:

 (1.28)

where:

 (1.29)

For its moment generating function takes the form:

 (1.30)

Note, that is smooth, i.e. infinitely many times differentiable, near 0 and hence every moment exists. If we set then the first two moments lead to the following formulas for the mean and variance of a generalized hyperbolic random variable:

 (1.31) Var (1.32)

The normal-inverse Gaussian (NIG) distributions were introduced by [5] as a subclass of the generalized hyperbolic laws obtained for . The density of the normal-inverse Gaussian distribution is given by:

 (1.33)

In XploRe the NIG density and distribution functions are implemented in the pdfnig and cdfnig quantlets, respectively. Like for the hyperbolic distribution the calculation of the PDF is straightforward, but the CDF has to be numerically integrated from (1.33).

At the ''expense'' of four parameters, the NIG distribution is able to model symmetric and asymmetric distributions with possibly long tails in both directions. Its tail behavior is often classified as ''semi-heavy'', i.e. the tails are lighter than those of non-Gaussian stable laws, but much heavier than Gaussian. Interestingly, if we let  tend to zero the NIG distribution converges to the Cauchy distribution (with location parameter  and scale parameter ), which exhibits extremely heavy tails. The tail behavior of the NIG density is characterized by the following asymptotic relation:

 (1.34)

In fact, this is a special case of a more general relation with the exponent of  being equal to (instead of ), which is valid for all generalized hyperbolic laws ([6]). Obviously, the NIG distribution may not be adequate to deal with cases of extremely heavy tails such as those of Pareto or non-Gaussian stable laws. However, empirical experience suggests an excellent fit of the NIG law to financial data ([52,62,89,97]). Moreover, the class of normal-inverse Gaussian distributions possesses an appealing feature that the class of hyperbolic laws does not have. Namely, it is closed under convolution, i.e. a sum of two independent NIG random variables is again NIG ([5]). In particular, if  and  are independent normal inverse Gaussian random variables with common parameters  and  but having different scale and location parameters  and , respectively, then is NIG. This feature is especially useful in time scaling of risks, e.g. in deriving -day risks from daily risks. Only two subclasses of the generalized hyperbolic distributions are closed under convolution. The other class with this important property is the class of variance-gamma (VG) distributions, which is a limiting case obtained for . The variance-gamma distributions (with ) were introduced to the financial literature by [63].

## 1.3.1 Simulation of Generalized Hyperbolic Variables

The most natural way of simulating generalized hyperbolic variables stems from the fact that they can be represented as normal variance-mean mixtures. Since the mixing distribution is the generalized inverse Gaussian law, the resulting algorithm reads as follows:

1. simulate a random variable    GIG   GIG;

2. simulate a standard normal random variable , e.g. using the Box-Muller algorithm, see Sect. 1.2.3;

3. return
The algorithm is fast and efficient if we have a handy way of simulating generalized inverse Gaussian variates. For , i.e. when sampling from the so-called inverse Gaussian (IG) distribution, there exists an efficient procedure that utilizes a transformation yielding two roots. It starts with the observation that if we let then the GIG   IG density, see (1.20), of  can be written as:

Now, following [92] we may write:

 (1.35)

i.e. is distributed as a chi-square random variable with one degree of freedom. As such it can be simply generated by taking a square of a standard normal random number. Unfortunately, the value of  is not uniquely determined by (1.35). Solving this equation for  yields two roots:

 and

The difficulty in generating observations with the desired distribution now lies in choosing between the two roots. [74] showed that  can be simulated by choosing  with probability . So for each random observation  from a  distribution the smaller root  has to be calculated. Then an auxiliary Bernoulli trial is performed with probability . If the trial results in a ''success'', is chosen; otherwise, the larger root is selected. The rndnig quantlet of XploRe, as well as the rnig function of the Rmetrics collection of software packages for S-plus/R (see also Sect. 1.2.2 where Rmetrics was briefly described), utilize this routine.

In the general case, the GIG distribution - as well as the (generalized) hyperbolic law - can be simulated via the rejection algorithm. An adaptive version of this technique is used to obtain hyperbolic random numbers in the rhyp function of Rmetrics. Rejection is also implemented in the HyperbolicDist package for S-plus/R developed by David Scott, see the R-project home page http://cran.r-project.org/. The package utilizes a version of the algorithm proposed by [3], i.e. rejection coupled either with a two (''GIG algorithm'' for any admissible value of ) or a three part envelope (''GIGLT1 algorithm'' for ). Envelopes, also called hat or majorizing functions, provide an upper limit for the PDF of the sampled distribution. The proper choice of such functions can substantially increase the speed of computations, see Chap. II.2. As [3] shows, once the parameter values for these envelopes have been determined, the algorithm efficiency is reasonable for most values of the parameter space. However, finding the appropriate parameters requires optimization and makes the technique burdensome.

This difficulty led to a search for a short algorithm which would give comparable efficiencies but without the drawback of extensive numerical optimizations. A solution, based on the ''ratio-of-uniforms'' method, was provided a few years later by [25]. First, recalling properties (1.21) and (1.22), observe that we only need to find a method to simulate    GIG variables and only for . Next, define the relocated variable , where is the mode of the density of  . Then the relocated variable can be generated by taking , where the pair is uniformly distributed over the region with:

 for

Since this region is irregularly shaped, it is more convenient to generate the pair uniformly over a minimal enclosing rectangle . Finally, the variate is accepted if . The efficiency of the algorithm depends on the method of deriving and the actual choice of and . Further, for and there is no need for the shift at mode . Such a version of the algorithm is implemented in the *gigru* functions of UNU.RAN, a library of C functions for non-uniform random number generation developed at the Department for Statistics, Vienna University of Economics, see http://statistik.wu-wien.ac.at/unuran/. It is also implemented in the gigru function of the SSC library (a Stochastic Simulation library in C developed originally by Pierre L'Ecuyer, see http://www.iro.umontreal.ca/~lecuyer and Chap. II.2) and in the rndghd quantlet of XploRe.

## 1.3.2 Estimation of Parameters

### 1.3.2.1 Maximum Likelihood Method

The parameter estimation of generalized hyperbolic distributions can be performed by the maximum likelihood method, since there exist closed-form formulas (although, involving special functions) for the densities of these laws. The computational burden is not as heavy as for -stable laws, but it still is considerable.

In general, the maximum likelihood estimation algorithm is as follows. For a vector of observations , the ML estimate of the parameter vector is obtained by maximizing the log-likelihood function:

 K (1.36)

where is defined by (1.29). Obviously, for hyperbolic ( ) distributions the algorithm uses simpler expressions of the log-likelihood function due to relation (1.24).

The routines proposed in the literature differ in the choice of the optimization scheme. The first software product that allowed statistical inference with hyperbolic distributions - the HYP program - used a gradient search technique, see [10]. In a large simulation study [84] utilized the bracketing method. The XploRe quantlets mlehyp and mlenig use yet another technique - the downhill simplex method of [77], with slight modifications due to parameter restrictions.

The main factor for the speed of the estimation is the number of modified Bessel functions to compute. Note, that for (i.e. the hyperbolic distribution) this function appears only in the constant . For a data set with  independent observations we need to evaluate  and Bessel functions for NIG and generalized hyperbolic distributions, respectively, whereas only one for the hyperbolic. This leads to a considerable reduction in the time necessary to calculate the likelihood function in the hyperbolic case. [84] reported a reduction of ca. , however, the efficiency results are highly sample and implementation dependent. For example, limited simulation studies performed in XploRe revealed a , and reduction in CPU time for samples of size , and , respectively.

We also have to say that the optimization is challenging. Some of the parameters are hard to separate since a flat-tailed generalized hyperbolic distribution with a large scale parameter is hard to distinguish from a fat-tailed distribution with a small scale parameter, see [6] who observed such a behavior already for the hyperbolic law. The likelihood function with respect to these parameters then becomes very flat, and may have local mimima. In the case of NIG distributions [97] proposed simple estimates of  and  that can be used as staring values for the ML scheme. Starting from relation (1.34) for the tails of the NIG density they derived the following approximation:

where is the -th population quantile, see Sect. 1.2.4. After the choice of a suitable value for , [97] used , the ''tail estimates'' of  and  are obtained by replacing the quantiles and expectations by their sample values in the above relations.

Another method of providing the starting values for the ML scheme was suggested by [84]. He estimated the parameters of a symmetric ( ) generalized hyperbolic law with a reasonable kurtosis (i.e. with ) that had the variance equal to that of the empirical distribution.

### 1.3.2.2 Other Methods

Besides the ML approach other estimation methods have been proposed in the literature. [84] tested different estimation techniques by replacing the log-likelihood function with other score functions, like the Anderson-Darling and Kolmogorov statistics or -norms. But the results were disappointing. [62] made use of the Markov chain Monte Carlo technique (see Chap. II.3), however, again the results obtained were not impressive. [52] described an EM type algorithm (see Chap. II.5) for maximum likelihood estimation of the normal inverse Gaussian distribution. The algorithm can be programmed in any statistical package supporting Bessel functions and it has all the properties of the standard EM algorithm, like sure, but slow, convergence, parameters in the admissible range, etc. The EM scheme can be also generalized to the family of generalized hyperbolic distributions.

## 1.3.3 Are Asset Returns NIG Distributed?

It is always necessary to find a reasonable tradeoff between the introduction of additional parameters and the possible improvement of the fit. [6] mentioned the flatness of the likelihood function for the hyperbolic distribution. The variation in the likelihood function of the generalized hyperbolic distribution is even smaller for a wide range of parameters. Consequently, the generalized hyperbolic distribution applied as a model for financial data leads to overfitting ([84]). In the empirical analysis that follows we will thus concentrate only on NIG distributions. They possess nice analytic properties and have been reported to fit financial data better than hyperbolic laws ([52,62,97]).

 Parameters or NIG fit () Gaussian fit () Test values Anderson-Darling Kolmogorov NIG fit Gaussian fit INF VaR estimates ( ) Empirical NIG fit ( ) ( ) Gaussian fit ( ) ( )

Now, we can return to the empirical analysis. This time we want to check whether DJIA and/or DAX returns can be approximated by the NIG distribution. We estimate the parameters using the maximum likelihood approach. As can be seen in Fig. 1.9 the fitted NIG distribution ''misses'' the very extreme DJIA returns. However, it seems to give a better fit to the central part of the empirical distribution than the -stable law. This is confirmed by a lower value of the Kolmogorov statistics, compare Tables 1.2 and 1.4. Surprisingly, also the Anderson-Darling statistics returns a lower value, implying a better fit in the tails of the distribution as well.

In the right panel of Fig. 1.9 we also plotted vertical lines representing the NIG, Gaussian and empirical daily VaR estimates at the and confidence levels. These estimates correspond to a one day VaR of a virtual portfolio consiting of one long position in the DJIA index. The NIG VaR estimate matches the empirical VaR almost perfectly and the NIG VaR estimate also yields a smaller difference than the stable estimate, compare Tables 1.2 and 1.4. However, if we were interested in very high confidence levels (i.e. very low quantiles) then the NIG fit would be less favorable than the stable one. Like in the stable case, no simple algorithms for inverting the NIG CDF are known but finding the right quantile could be performed through a binary search routine. For some members of the generalized hyperbolic family specialized inversion techniques have been developed. For example, [59] showed that the approximate inverse of the hyperbolic CDF can be computed as the solution of a first-order differential equation.

The second analyzed data set comprises  returns of the Deutsche Aktienindex (DAX) index. In this case the NIG distribution offers an indisputably better fit than the Gaussian or even the -stable law, see Table 1.5 and compare with Table 1.3. This can be also seen in Fig. 1.10. The ''drop off'' in the left tail of the empirical distribution is nicely caught by the NIG distribution. The empirical VaR estimates are also ''caught'' almost perfectly.

 Parameters or NIG fit () Gaussian fit () Test values Anderson-Darling Kolmogorov NIG fit Gaussian fit VaR estimates ( ) Empirical NIG fit ( ) ( ) Gaussian fit ( ) ( )

Next: 1.4 Value at Risk, Up: 1. Computationally Intensive Value Previous: 1.2 Stable Distributions