Next: 1.3 Hyperbolic Distributions Up: 1. Computationally Intensive Value Previous: 1.1 Introduction

Subsections

# 1.2 Stable Distributions

It is often argued that financial asset returns are the cumulative outcome of a vast number of pieces of information and individual decisions arriving almost continuously in time ([71,87]). As such, since the pioneering work of Louis Bachelier in 1900, they have been modeled by the Gaussian distribution. The strongest statistical argument for it is based on the Central Limit Theorem, which states that the sum of a large number of independent, identically distributed variables from a finite-variance distribution will tend to be normally distributed. However, financial asset returns usually have heavier tails.

In response to the empirical evidence [64] and [36] proposed the stable distribution as an alternative model. There are at least two good reasons for modeling financial variables using stable distributions. Firstly, they are supported by the generalized Central Limit Theorem, which states that stable laws are the only possible limit distributions for properly normalized and centered sums of independent, identically distributed random variables ([58]). Secondly, stable distributions are leptokurtic. Since they can accommodate the fat tails and asymmetry, they fit empirical distributions much better.

Stable laws - also called -stable, stable Paretian or Lévy stable - were introduced by [61] during his investigations of the behavior of sums of independent random variables. A sum of two independent random variables having an -stable distribution with index  is again -stable with the same index . This invariance property, however, does not hold for different 's.

The -stable distribution requires four parameters for complete description: an index of stability also called the tail index, tail exponent or characteristic exponent, a skewness parameter , a scale parameter and a location parameter . The tail exponent  determines the rate at which the tails of the distribution taper off, see the left panel of Fig. 1.3. When , a Gaussian distribution results. When , the variance is infinite and the tails are asymptotically equivalent to a Pareto law, i.e. they exhibit a power-law behavior. More precisely, using a central limit theorem type argument it can be shown that ([48,90]):

 (1.1)

where:

When , the mean of the distribution exists and is equal to . In general, the th moment of a stable random variable is finite if and only if . When the skewness parameter  is positive, the distribution is skewed to the right, i.e. the right tail is thicker, see the right panel of Fig. 1.3. When it is negative, it is skewed to the left. When , the distribution is symmetric about . As approaches ,  loses its effect and the distribution approaches the Gaussian distribution regardless of . The last two parameters, and , are the usual scale and location parameters, i.e. determines the width and  the shift of the mode (the peak) of the distribution.

## 1.2.1 Characteristic Function Representation

From a practitioner's point of view the crucial drawback of the stable distribution is that, with the exception of three special cases, its probability density function (PDF) and cumulative distribution function (CDF) do not have closed form expressions. These exceptions include the well known Gaussian () law, whose density function is given by:

 (1.2)

and the lesser known Cauchy (, ) and Lévy ( , ) laws.

Hence, the -stable distribution can be most conveniently described by its characteristic function  - the inverse Fourier transform of the PDF. However, there are multiple parameterizations for -stable laws and much confusion has been caused by these different representations. The variety of formulas is caused by a combination of historical evolution and the numerous problems that have been analyzed using specialized forms of the stable distributions. The most popular parameterization of the characteristic function of , i.e. an -stable random variable with parameters , , and , is given by ([90,98]):

 (1.3)

Note, that the traditional scale parameter of the Gaussian distribution is not the same as  in the above representation. A comparison of formulas (1.2) and (1.3) yields the relation: .

For numerical purposes, it is often useful to use Nolan's (1997) parameterization:

 (1.4)

The representation is a variant of Zolotarev's 1986 (M)-parameterization, with the characteristic function and hence the density and the distribution function jointly continuous in all four parameters, see Fig. 1.4. In particular, percentiles and convergence to the power-law tail vary in a continuous way as  and  vary. The location parameters of the two representations are related by for and for .

## 1.2.2 Computation of Stable Density and Distribution Functions

The lack of closed form formulas for most stable densities and distribution functions has negative consequences. Numerical approximation or direct numerical integration have to be used, leading to a drastic increase in computational time and loss of accuracy. Of all the attempts to be found in the literature a few are worth mentioning. [29] developed a procedure for approximating the stable distribution function using Bergström's (1952) series expansion. Depending on the particular range of  and , [46] combined four alternative approximations to compute the stable density function. Both algorithms are computationally intensive and time consuming, making maximum likelihood estimation a nontrivial task, even for modern computers. Recently, two other techniques have been proposed.

[75] exploited the density function - characteristic function relationship and applied the fast Fourier transform (FFT). However, for data points falling between the equally spaced FFT grid nodes an interpolation technique has to be used. The authors suggested that linear interpolation suffices in most practical applications, see also [87]. Taking a larger number of grid points increases accuracy, however, at the expense of higher computational burden. Setting the number of grid points to and the grid spacing to allows to achieve comparable accuracy to the direct integration method (see below), at least for a range of 's typically found for financial data ( ). As for the computational speed, the FFT based approach is faster for large samples, whereas the direct integration method favors small data sets since it can be computed at any arbitrarily chosen point. [75] report that for the FFT based method is faster for samples exceeding  observations and slower for smaller data sets.

We must stress, however, that the FFT based approach is not as universal as the direct integration method - it is efficient only for large alpha's and only as far as the probability density function calculations are concerned. When computing the cumulative distribution function the former method must numerically integrate the density, whereas the latter takes the same amount of time in both cases.

The direct integration method, proposed by Nolan 1997, 1999) consists of a numerical integration of Zolotarev's (1986) formulas for the density or the distribution function. To save space we state only the formulas for the probability density function. Complete formulas can be also found in [12].

Set . Then the density of a standard -stable random variable in representation , i.e. , can be expressed as (note, that Zolotarev 1986, Sect. 2.2) used yet another parametrization):

• when and :

• when and :

• when and :

• when :

where

 (1.5)

and

XploRe offers the direct integration method through the cdfstab and pdfstab quantlets, see [22] for a thorough exposition of quantlets related to stable distributions. On a PC equipped with a Centrino GHz processor the calculation of the stable distribution or density function at  points takes about  seconds. As default, the integrals found in the above formulas are numerically integrated using  subintervals. These computational times can be improved when using a numerically more efficient environment. For example, the program STABLE (downloadable from John Nolan's web page: http://academic2.american.edu/~jpnolan/stable/stable.html) needs about  seconds for performing corresponding calculations. It was written in Fortran and calls several external IMSL routines, see [79] for details. Apart from speed, the STABLE program also exhibits higher relative accuracy (for default tolerance settings in both programs): about compared to for values used in typical financial applications (like approximating asset return distributions). Naturally, the accuracy of both programs can be increased at the cost of computational time.

It is interesting to note, that currently no other statistical computing environment offers the computation of stable density and distribution functions in its standard release. Users have to rely on third-party libraries or commercial products. A few are worth mentioning. John Nolan offers the STABLE program in library form through Robust Analysis Inc., see http://www.robustanalysis.com. This library (in C, Visual Basic or Fortran) provides interfaces to Matlab, S-plus (or its GNU version - R) and Mathematica. Diethelm Würtz has developed Rmetrics, an open source collection of software packages for S-plus/R, which may be useful for teaching computational finance, see http://www.itp.phys.ethz.ch/econophysics/R/. Stable PDF and CDF calculations are performed using the direct integration method, with the integrals being computed by R's function integrate. Interestingly, for symmetric stable distributions Rmetrics utilizes McCulloch's (1998) approximation, which was obtained by interpolating between the complements of the Cauchy and Gaussian CDFs in a transformed space. For the absolute precision of the stable PDF and CDF approximation is . The FFT based approach is utilized in Cognity, a commercial risk management platform that offers derivatives pricing and portfolio optimization based on the assumption of stably distributed returns, see http://www.finanalytica.com.

## 1.2.3 Simulation of -stable Variables

The complexity of the problem of simulating sequences of -stable random variables stems from the fact that there are no analytic expressions for the inverse nor the cumulative distribution function . All standard approaches like the rejection or the inversion methods would require tedious computations. See Chap. II.2 for a review of non-uniform random number generation techniques.

A much more elegant and efficient solution was proposed by [21]. They noticed that a certain integral formula derived by [100] yielded the following algorithm:

• generate a random variable uniformly distributed on and an independent exponential random variable  with mean ;

• for compute:

 (1.6)

• for compute:

 (1.7)

where is given by (1.5). This algorithm yields a random variable , in representation (1.3). For a detailed proof see [98].

Given the formulas for simulation of a standard -stable random variable, we can easily simulate a stable random variable for all admissible values of the parameters , , and  using the following property: if then

is . It is interesting to note that for (and ) the Chambers-Mallows-Stuck method reduces to the well known [14] algorithm for generating Gaussian random variables ([49]).

Many other approaches have been proposed in the literature, including application of [8] and LePage ([60]) series expansions, see [65] and [47], respectively. However, this method is regarded as the fastest and the most accurate. In XploRe the algorithm is implemented in the rndstab quantlet. On a PC equipped with a Centrino GHz processor one million variables are generated in about  seconds, compared to about  seconds for one million standard normal random variables obtained via the Box-Muller algorithm (normal2). Because of its unquestioned superiority and relative simplicity, the Chambers-Mallows-Stuck method is implemented in some statistical computing environments (e.g. the rstable function in S-plus/R) even if no other routines related to stable distributions are provided.

## 1.2.4 Estimation of Parameters

The estimation of stable law parameters is in general severely hampered by the lack of known closed-form density functions for all but a few members of the stable family. Numerical approximation or direct numerical integration are nontrivial and burdensome from a computational point of view. As a consequence, the maximum likelihood (ML) estimation algorithm based on such approximations is difficult to implement and time consuming for samples encountered in modern finance. However, there are also other numerical methods that have been found useful in practice and are discussed in this section.

All presented methods work quite well assuming that the sample under consideration is indeed -stable. Since there are no formal tests for assessing the -stability of a data set we suggest to first apply the ''visual inspection'' and tail exponent estimators to see whether the empirical densities resemble those of -stable laws ([12,99]).

Given a sample of independent and identically distributed (i.i.d.) observations, in what follows, we provide estimates , , and of all four stable law parameters. We start the discussion with the simplest, fastest and  least accurate quantile methods, then develop the slower, yet much more accurate sample characteristic function methods and, finally, conclude with the slowest but most accurate maximum likelihood approach.

### 1.2.4.1 Sample Quantile Methods

[37] provided very simple estimates for parameters of symmetric ( ) stable laws with . They proposed to estimate  by:

 (1.8)

where denotes the -th population quantile, so that . As [70] noticed, Fama and Roll based their estimator of on the fortuitous observation that lies within of for all , when . This enabled them to estimate by (1.8) with less than asymptotic bias without first knowing . However, when , the search for an invariant range such as the one they found becomes futile.

The characteristic exponent , on the other hand, can be estimated from the tail behavior of the distribution. Fama and Roll suggested to take  satisfying:

 (1.9)

They found that , , worked best for estimating . This method unnecessarily compounds the small asymptotic bias in the estimator of  into the estimator of .

For , the stable distribution has finite mean. Hence, the sample mean is a consistent estimate of the location parameter . However, a more robust estimate is the truncated sample mean - the arithmetic mean of the middle  percent of the ranked observations. The truncated mean is often suggested in the literature when the range of  is unknown.

Fama and Roll's (1971) method is simple but suffers from a small asymptotic bias in and and restrictions on  and . [70] generalized and improved the quantile method. He analyzed stable law quantiles and provided consistent estimators of all four stable parameters, with the restriction , while retaining the computational simplicity of Fama and Roll's method. After McCulloch define:

 and (1.10)

Statistics and are functions of and  only, i.e. they are independent of both  and . This relationship may be inverted and the parameters  and  may be viewed as functions of and . Substituting and by their sample values and applying linear interpolation between values found in tables given in [70] yields estimators and  .

Scale and location parameters, and , can be estimated in a similar way. However, due to the discontinuity of the characteristic function for and in representation (1.3), this procedure is more complicated. We refer the interested reader to the original work of [70]. This estimation technique is implemented in XploRe in the stabcull quantlet.

### 1.2.4.2 Sample Characteristic Function Methods

Given an i.i.d. random sample of size , define the sample characteristic function by:

Since is bounded by unity all moments of are finite and, for any fixed , it is the sample average of i.i.d. random variables . Hence, by the law of large numbers, is a consistent estimator of the characteristic function .

[85] proposed a simple estimation method, called the method of moments, based on transformations of the characteristic function. From (1.3) we have for all :

 (1.11)

Hence, . Now, assuming , choose two nonzero values of , say . Then for we have:

 (1.12)

Solving these two equations for and , and substituting for yields:

 (1.13)

and

 (1.14)

In order to estimate and we have to choose two nonzero values of , say , and apply a similar trick to . The estimators are consistent since they are based upon estimators of , and , which are known to be consistent. However, convergence to the population values depends on the choice of . The optimal selection of these values is problematic and still is an open question. The XploRe implementation of the method of moments (the stabmom quantlet) uses , , , and as proposed by [56] in his simulation study.

In the same paper Koutrouvelis presented a much more accurate regression-type method which starts with an initial estimate of the parameters and proceeds iteratively until some prespecified convergence criterion is satisfied. Each iteration consists of two weighted regression runs. The number of points to be used in these regressions depends on the sample size and starting values of . Typically no more than two or three iterations are needed. The speed of the convergence, however, depends on the initial estimates and the convergence criterion.

The regression method is based on the following observations concerning the characteristic function . First, from (1.3) we can easily derive:

 (1.15)

The real and imaginary parts of are for given by:

and

The last two equations lead, apart from considerations of principal values, to:

 sign (1.16)

Equation (1.15) depends only on and and suggests that we estimate these parameters by regressing on in the model:

 (1.17)

where is an appropriate set of real numbers, , and denotes an error term. [56] proposed to use ; with ranging between  and  for different values of  and sample sizes.

Once and have been obtained and and have been fixed at these values, estimates of  and  can be obtained using (1.16). Next, the regressions are repeated with , , and as the initial parameters. The iterations continue until a prespecified convergence criterion is satisfied. Koutrouvelis proposed to use the Fama-Roll estimator (1.8) and the truncated mean for initial estimates of  and , respectively.

[54] eliminated this iteration procedure and simplified the regression method. For initial estimation they applied McCulloch's method, worked with the continuous representation (1.4) of the characteristic function instead of the classical one (1.3) and used a fixed set of only 10 equally spaced frequency points . In terms of computational speed their method compares favorably to the original method of Koutrouvelis, see Table 1.1. It has a significantly better performance near and due to the elimination of discontinuity of the characteristic function. However, it returns slightly worse results for other values of . In XploRe both regression algorithms are implemented in the stabreg quantlet. An optional parameter lets the user choose between the original Koutrouvelis code and the Kogon-Williams modification.

 Method CPU time McCulloch s ( ) ( ) ( ) ( ) Moments s ( ) ( ) ( ) ( ) Koutrouvelis s ( ) ( ) ( ) ( ) Kogon-Williams s ( ) ( ) ( ) ( )

A typical performance of the described estimators is summarized in Table 1.1, see also Fig. 1.5. McCulloch's quantile technique, the method of moments, the regression approach of Koutrouvelis and the method of Kogon and Williams were applied to 100 simulated samples of two thousand random numbers each. The method of moments yielded the worst estimates, clearly outside any admissible error range. McCulloch's method came in next with acceptable results and computational time significantly lower than the regression approaches. On the other hand, both the Koutrouvelis and the Kogon-Williams implementations yielded good estimators with the latter performing considerably faster, but slightly less accurate. We have to say, though, that all methods had problems with estimating . Like it or not, our search for the optimal estimation technique is not over yet. We have no other choice but turn to the last resort - the maximum likelihood method.

### 1.2.4.3 Maximum Likelihood Method

The maximum likelihood (ML) estimation scheme for -stable distributions does not differ from that for other laws, at least as far as the theory is concerned. For a vector of observations , the ML estimate of the parameter vector is obtained by maximizing the log-likelihood function:

 (1.18)

where is the stable density function. The tilde denotes the fact that, in general, we do not know the explicit form of the stable PDF and have to approximate it numerically. The ML methods proposed in the literature differ in the choice of the approximating algorithm. However, all of them have an appealing common feature - under certain regularity conditions the maximum likelihood estimator is asymptotically normal with the variance specified by the Fischer information matrix ([30]). The latter can be approximated either by using the Hessian matrix arising in maximization or, as in [81], by numerical integration.

Because of computational complexity there are only a few documented attempts of estimating stable law parameters via maximum likelihood. [29] developed an approximate ML method, which was based on grouping the data set into bins and using a combination of means to compute the density (FFT for the central values of  and series expansions for the tails) to compute an approximate log-likelihood function. This function was then numerically maximized.

Applying Zolotarev's (1964) integral formulas, [17] formulated another approximate ML method, however, only for symmetric stable random variables. To avoid the discontinuity and nondifferentiability of the symmetric stable density function at , the tail index  was restricted to be greater than one.

Much better, in terms of accuracy and computational time, are more recent maximum likelihood estimation techniques. [76] utilized the FFT approach for approximating the stable density function, whereas [81] used the direct integration method. Both approaches are comparable in terms of efficiency. The differences in performance are the result of different approximation algorithms, see Sect. 1.2.2.

As [83] observes, the ML estimates are almost always the most accurate, closely followed by the regression-type estimates, McCulloch's quantile method, and finally the method of moments. However, as we have already said in the introduction to this section, maximum likelihood estimation techniques are certainly the slowest of all the discussed methods. For example, ML estimation for a sample of  observations using a gradient search routine which utilizes the direct integration method needs  seconds or about  minutes! The calculations were performed on a PC equipped with a Centrino GHz processor and running STABLE ver.  (see also Sect. 1.2.2 where the program was briefly described). For comparison, the STABLE implementation of the Kogon-Williams algorithm performs the same calculations in only seconds (the XploRe quantlet stabreg needs roughly four times more time, see Table 1.1). Clearly, the higher accuracy does not justify the application of ML estimation in many real life problems, especially when calculations are to be performed on-line. For this reason the program STABLE also offers an alternative - a fast quasi ML technique. It quickly approximates stable densities using a -dimensional spline interpolation based on pre-computed values of the standardized stable density on a grid of values. At the cost of a large array of coefficients, the interpolation is highly accurate over most values of the parameter space and relatively fast -  seconds for a sample of  observations.

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

## 1.2.5 Are Asset Returns -stable?

In this paragraph we want to apply the discussed techniques to financial data. Due to limited space we have chosen only one estimation method - the regression approach of [56], as it offers high accuracy at moderate computational time. We start the empirical analysis with the most prominent example - the Dow Jones Industrial Average (DJIA) index. The data set covers the period January 2, 1985 - November 30, 1992 and comprises  returns. Recall, that this period includes the largest crash in Wall Street history - the Black Monday of October 19, 1987. Clearly the -stable law offers a much better fit to the DJIA returns than the Gaussian distribution, see Table 1.2. Its superiority, especially in the tails of the distribution, is even better visible in Fig.  1.6. In this figure we also plotted vertical lines representing the -stable, 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 stable VaR estimates are almost twice closer to the empirical estimates than the Gaussian ones, see Table 1.2.

Recall that calculating the VaR number reduces to finding the quantile of a given distribution or equivalently to evaluating the inverse of the distribution function at . Unfortunately no simple algorithms for inverting the stable CDF are known. The qfstab quantlet of XploRe utilizes a simple binary search routine for large 's and values near the mode of the distribution. In the extreme tails the approximate range of the quantile values is first estimated via the power law formula (1.1), then a binary search is conducted.

To make our statistical analysis more sound, we also compare both fits through Anderson-Darling and Kolmogorov test statistics ([24]). The former may be treated as a weighted Kolmogorov statistics which puts more weight to the differences in the tails of the distributions. Although no asymptotic results are known for the stable laws, approximate critical values for these goodness-of-fit tests can be obtained via the bootstrap technique ([12,94]). In this chapter, though, we will not perform hypothesis testing and just compare the test values. Naturally, the lower the values the better the fit. The stable law seems to be tailor-cut for the DJIA index returns. But does it fit other asset returns as well?

The second analyzed data set comprises  returns of the Deutsche Aktienindex (DAX) index from the period January 2, 1995 - December 5, 2002. Also in this case the -stable law offers a much better fit than the Gaussian, see Table 1.3. However, the test statistics suggest that the fit is not as good as for the DJIA returns (observe that both data sets are of the same size and the test values in both cases can be compared). This can be also seen in Fig. 1.7. The left tail seems to drop off at some point and the very tail is largely overestimated by the stable distribution. At the same time it is better approximated by the Gaussian law. This results in a surprisingly good fit of the daily VaR by the Gaussian distribtion, see Table 1.3, and an overestimation of the daily VaR estimate at the confidence level by the -stable distribtion. In fact, the latter is a rather typical situation. For a risk manager who likes to play safe this may not be a bad idea, as the stable law overestimates the risks and thus provides an upper limit of losses.

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

This example clearly shows that the -stable distribution is not a panacea. Although it gives a very good fit to a number of empirical data sets, there surely are distributions that recover the characteristics of other data sets better. We devote the rest of this chapter to such alternative heavy tailed distributions. We start with a modification of the stable law and in Sect. 1.3 concentrate on the class of generalized hyperbolic distributions.

## 1.2.6 Truncated Stable Distributions

Mandelbrot's (1963) pioneering work on applying -stable distributions to asset returns gained support in the first few years after its publication ([36,82,95]). Subsequent works, however, have questioned the stable distribution hypothesis ([1,11]). By the definition of the stability property, the sum of i.i.d. stable random variables is also stable. Thus, if short term asset returns are distributed according to a stable law, longer term returns should retain the same functional form. However, from the empirical data it is evident that as the time interval between price observations grows longer, the distribution of returns deviates from the short term heavy tailed distribution, and converges to the Gaussian law. This indicates that the returns probably are not -stable (but it could mean as well that the returns are just not independent). Over the next few years, the stable distribution temporarily lost favor and alternative processes were suggested as mechanisms generating stock returns.

In mid 1990s the stable distribution hypothesis has made a dramatic comeback. Several authors have found a very good agreement of high-frequency returns with a stable distribution up to six standard deviations away from the mean ([23,67]). For more extreme observations, however, the distribution they have found falls off approximately exponentially. To cope with such observations the truncated Lévy distributions (TLD) were introduced by [66]. The original definition postulated a sharp truncation of the -stable probability density function at some arbitrary point. However, later an exponential smoothing was proposed by [55].

For the characteristic function of a symmetric TLD random variable is given by:

where is the tail exponent, is the scale parameter and  is the truncation coefficient. Clearly the TLD reduces to the symmetric -stable distribution ( ) when . The TLD distribution exhibits the following behavior: for small and intermediate returns it behaves like a stable distribution, but for extreme returns the truncation causes the distribution to converge to a Gaussian distribution. Thus the observation that the asset returns distribution is a TLD explains both the short-term -stable behavior and the long run convergence to the normal distribution.

Despite these interesting features the truncated Lévy distributions have not been applied extensively to date. The most probable reason for this being the complicated definition of the TLD law. Like for -stable distributions, only the characteristic function is known. No closed form formulas exist for the density or the distribution function. Since no integral formulas, like Zolotarev's (1986) for the -stable laws, have been discovered as yet, statistical inference is, in general, limited to maximum likelihood utilizing the FFT technique for approximating the PDF. Moreover, compared to the stable distribution, the TLD introduces one more parameter (asymmetric TLD laws have also been considered in the literature, see e.g. [15] and [55]) making the estimation procedure even more complicated. Other parameter fitting techniques proposed so far comprise a combination of ad hoc approaches and moment matching ([15,69]). Better techniques have to be discovered before TLDs become a common tool in finance.

Next: 1.3 Hyperbolic Distributions Up: 1. Computationally Intensive Value Previous: 1.1 Introduction