Stable Distributions in Mathematica

Robert H. Rimmer MD
HC 4, Box 3-C
Payson, AZ 85541-9571
rrimmer@npgcable . com

John P. Nolan PhD							
Math/Stat Department
American University
4400 Massachusetts Avenue, NW
Washington, DC 20016-8050
jpnolan@american.edu

Introduction

Stable distributions are a rich class of probability distributions that allow skewness and heavy tails and have many intriguing mathematical properties.  The class was characterized by Paul Levy in his study of sums of independent identically distributed terms in the 1920s.  Stable distributions have been proposed as a model for many types of physical and economic systems.  There are solid theoretical reasons for expecting a non-Gaussian stable model, e.g. reflection off a rotating mirror yielding a Cauchy distribution, hitting times for a Brownian motion yielding a Levy distribution, the gravitational field of stars yielding the Holtsmark distribution.  The Generalized Central Limit Theorem states that the only possible non-trivial limit of normalized sums of independent identically distributed terms is stable.  It is argued that some observed quantities are the sum of many small terms; hence a stable model should be used to describe such systems.   For a detailed exposition of stable distributions and their applications, see Nolan (2005).  

The lack of formulas for densities and distribution functions for the full parameter range,  has limited the use of stable distributions by practitioners.  Recently high speed processors have made numerical calculations of stable functions feasible, and interest in their use has been increasing, especially in the study of financial markets based on Mandelbrot's work.  Rose (2002) has used Mathematica to generate special functions to calculate stable densities and cumulative distributions, but these special functions exist only for certain rational values of stable parameters.  This program takes advantage of the Zolotarev integral representation of the stable density and the capabilities of the function, NIntegrate, to accurately calculate stable distributions for nearly the whole parameter range.  Its availability in Mathematica should make stable distributions accessible to a wide audience across all commonly used computing environments.

Parameterizations of Stable Laws

A general stable distribution requires four parameters to describe it:

	index of stability or characteristic exponent α ∈ (0, 2],
	skewness parameter β ∈ [-1, 1], 
	scale parameter γ > 0,
	a location parameter δ ∈ ℝ. 

There are multiple parameterizations for stable laws and much confusion has been caused by these different parameterizations.  This package includes two parameterizations.  Below are the characteristic functions in the S(α, β, γ, δ; 0) and S(α, β, γ, δ; 1) parameterizations.

S(α, β, γ, δ; 0):           	 \(\(\(e\^\(i\ u\ δ - γ\^α\ Abs[u]\^α\ \((1 + i\ β\ \((\(-1\) + Abs[u\ γ]\^\(1 - α\))\)\ Sign[
                      u]\ Tan[\(π\ α\)\/2])\)\)\)\(\ \ \ \ \)\)\)α ≠ 1
		
				\(e\^\(i\ u\ δ - \(γ\ Abs[u]\ \((π + 2\ i\ β\ Log[Abs[u\ γ]]\ Sign[u])\)\)\/π\)\)     		 α == 1
and

S(α, β, γ, δ; 1):   		\(e\^\(i\ u\ δ - γ\^α\ Abs[u]\^α\ \((1 - i\ β\ Sign[u]\ Tan[\(π\ α\)\/2])\)\)\)             		 α ≠ 1
		
				\(e\^\(i\ u\ δ - γ\ Abs[u]\ \((1 + \(2\ i\ β\ Log[Abs[u]]\ Sign[u]\)\/π)\)\)\)          		α == 1

The S(α, β, γ, δ; 0) parameterization is the best choice for numerical computation of stable distributions, since it has the simplest form for the characteristic function that is continuous in all four parameters.  The S(α, β, γ, δ; 1) has been commonly used in economic literature and has the property that the location parameter,  δ, is the mean when α > 1.  The graphics array below demonstrates these two parameterizations.  The scale parameter, γ, has been set to 1 and the location parameter, δ is 0.  Each graph shows varying values of β {0, 1/3, 2/3, 1} in {black, red, green, blue} colors respectively at the value of α and type listed on the graph.  Negative values of β would give mirror image graphs.  When α is near 1, the mode in the S(α, β, γ, δ; 1) parameterization moves dramatically; the S(α, β, γ, δ; 0) parameterization progressively reduces this motion as alpha approaches one.

<< StableDistribution`SMath`
SetOptions[NIntegrate, 
        GaussPoints ->  30];
splot[a_, type_] := 
        Module[{style}, 
style  =  {{RGBColor[0, 0, 0], 
                Thickness[0.01]}, {RGBColor[1, 0, 0], 
                Thickness[0.01]}, {RGBColor[0, 1, 0], 
                Thickness[0.01]}, {RGBColor[0, 0, 1], 
                Thickness[0.01]}}; 
Plot[{SPDF[x, {a, 0}, 
                type], SPDF[x, {a, 1/3}, type], SPDF[x, {a, 2/3}, type], 
              SPDF[x, {a, 1}, type]}, {x, -5, 5}, 
            PlotStyle -> style, PlotRange ->  All, Frame -> True, 
            Axes -> False, 
            PlotLabel  -> 
              StringJoin["S(", ToString[a], ", β, 1, 0; ", 
                ToString[type], ")"], 
            DisplayFunction -> Identity]];
Show[
        GraphicsArray[
          Table[splot[a, t], {a, 0.8, 1.6, 0.4}, {t, 0, 
              1}]]];
Quit

Using the Package

The StableDistribution package provides several basic functions needed for the calculation of stable distributions.  They are all designed with identical input, output, and default formats.  The package comes with a help browser that explains in more detail how each function works; in addition there are a few more functions that are not described in this paper.  The package allows both the 0 and 1 parameterizations for all the functions listed below, but in the examples and discussion only S(α, β, γ, δ; 1) will be used.

<< StableDistribution`SMath`

Print[StringTake[SPDF::usage, 114]];
Print[StringTake[SCDF::usage, 139]]; 
    Print[StringTake[SQuantile::usage, 124]]; 
    Print[StringTake[SRandom::usage, 90]]; 
    Print[StringTake[fMLFit::usage, 123]]; 
    Print[StringTake[D1SPDF::usage, 142]];
Print[StringTake[SMode::usage, 108]];

The SPDF function calculates the stable density using a variation of the Zolotarev integral transformation, (Nolan 1997).  This transformation permits accurate calculations of both the density and distribution functions across nearly the whole domain of stable parameters.  Attempts to simply numerically integrate the inverse Fourier transform of the characteristic function run into difficulties when alpha is less than one or when x is large.

Plot[SPDF[x, {1.3, 0.5, 1, 0.5}], {x, -20, 20}];

Chambers, Mallows, and Stuck (1976) published an algorithm for the generation of  stable random variables.  It is implemented in this package as the function, SRandom.  The example below generates a stable random sequence of  10 numbers in the where p1 is a parameter set, {α, β, γ, δ}.

p1  =  {1.3, 0.5, 1, 0.5};
SRandom[10, p1]

The SRandom generator makes it easy to experience the behavior of stable distributions.  Below three stable random samples are generated.  Note how extreme outliers may be for lower values of alpha. The graphs are at the same scale to illustrate the differences in dispersion, but at this scale some of the extreme points in the case of α = 0.9 are not shown.  A stable sample with α = 2 has the same dispersion as the normal distribution.  α = 1.7 is typical of stock price returns.

s1 = SRandom[10000, {0.9, 0, 1, 0}, 1];
g1 = 
        ListPlot[s1, PlotRange -> {All, {-100, 100}}, 
          PlotStyle -> Hue[1], PlotLabel ->  "α = 0.9", 
          Frame -> True, FrameTicks -> {None, Automatic}, 
          DisplayFunction -> Identity];
s2 = SRandom[10000, {1.7, 0, 1, 0}, 1];
g2 = 
        ListPlot[s2, PlotRange -> {All, {-100, 100}}, 
          PlotStyle -> Hue[1], PlotLabel ->  "α = 1.7", 
          Frame -> True, FrameTicks -> {None, Automatic}, 
          DisplayFunction -> Identity];
s3 = SRandom[10000, {2, 0, 1, 0}, 1];
g3 = 
        ListPlot[s3, PlotRange -> {All, {-100, 100}}, 
          PlotStyle -> Hue[1], 
          PlotLabel ->  "α = 2 (Normal)", 
          Frame -> True, FrameTicks -> {None, Automatic}, 
          DisplayFunction -> Identity];
Show[GraphicsArray[{g1, g2, g3}]];

The example below simulates what a series of daily stock market prices might look like if generated by a stable return model with a starting price of 50 and an expected annual return of 7% (there are about 252 trading days in a year); notice the sudden large jumps.

ListPlot[
        50 Exp[FoldList[Plus, 0, 
              SRandom[1000, {1.7, -0.1, 0.012, Log[1.07]/252}]]], 
        PlotJoined ->  True];

The StableDistribution package also has data analysis capability and can perform a fit to data.  The fMLFit function provides a maximum likelihood fit to the data.  The algorithm uses a fast Fourier transform of the characteristic function (Rachev 2000) to approximate the stable probability density function rapidly.

p1 = {1.4, 0.9, 1, -0.5}
s4  =  SRandom[10000, p1];
p2 = fMLFit[s4]

The graph below compares the fit to the empirical density of the sample.

<< Graphics`Graphics`
empiricdens  =  
        Histogram[s4, HistogramScale -> 1, 
          HistogramRange -> {-8, 8}, 
          DisplayFunction -> Identity];
calcdens  =  
        Plot[SPDF[x, p2], {x, -8, 8}, 
          PlotStyle ->  {Hue[0.67], Thickness[0.008]}, 
          DisplayFunction ->  Identity];
Show[empiricdens, calcdens, 
        PlotLabel ->  "Empirical Density Fit", 
        DisplayFunction -> $DisplayFunction];

The plot of the distribution function is shown below.

Plot[SCDF[x, p2], {x, -8, 8}, 
          PlotStyle ->  {Hue[0.67], Thickness[0.008]}];

The SQuantile function uses FindRoot to solve for x at a value of the SCDF.  It can be used to print a table of cumulative probabilities for the data set in the above example.

Table[{p, SQuantile[p, p2]},  {p, 0, 1, 0.1}] // TableForm

The density and distribution functions use NIntegrate to make the calculations.  A direct integral form for the first derivative of the density function is also provided, but further derivatives can also be calculated using the derivative function in Mathematica.

stl = Table[Hue[i], {i, 0.6, 1, 0.1}];
f2[x_, a_] := D[SPDF[x, a], {x, 2}];
f3[x_, a_] := D[SPDF[x, a], {x, 3}];
Plot[
        Evaluate[{SCDF[x, p2], SPDF[x, p2], D1SPDF[x, p2], f2[x, p2], 
            f3[x, p2]}], {x, -5, 5}, PlotStyle -> stl, 
        PlotRange -> All];

While the package is able to calculate most all of the range of stable distributions, it can be at times be slow and may require adjustment of the input parameters to NIntegrate to obtain accurate results.  With such adjustments, however, it is possible to make extremely precise calculations.  The example below shows high precision density calculations for values on the light tail of a maximally skewed stable distribution.

Table[{2^i, 
            SPDF[2^i, {13/10, -1},  GaussPoints -> 30, 
              MaxRecursion -> 10, WorkingPrecision -> 64]}, {i, 
            6}] // TableForm

The SMode function can be used to locate the mode of a stable distribution, this is particularly useful when α is close to 1, where the mode may be quite distant from the origin.  The function uses FindRoot to determine where the first derivative of the stable density is equal to zero.

p3 =  {1.001, 0.5};
md = SMode[p3]
Plot[SPDF[x, p3], {x, md - 5, md + 5}];

Program Accuracy, Limitations, and Versatility

The functions illustrated above are designed to have the look and feel of true mathematical functions.  SPDF uses the Zolotarev transformation of the inverse Fourier transform of the characteristic function in the S(α, β, γ, δ; 0) parameterization.  The integrand is represented in the program symbolically so that transformation to the S(α, β, γ, δ; 1) parameterization is exact.  The SCDF and D1SPDF integrands are created by taking the integral and derivative of the SPDF integrand with respect to x before proceeding with transform integration.  The function, NIntegrate, gives a numerical approximation to the integral.  It is important to remember that NIntegrate may not always come up with the correct answer.  This is particularly true if the integrand has a sharp narrow spike.  The integrand for the stable density can be problematic.  The integrand is continuous and does not oscillate over the interval of integration; it begins with a value of 0, rises to a peak,  (1/e), and falls to 0 at the other end of the interval.  Unfortunately the peak for some values of x, α, and β can be quite narrow and it is possible for NIntegrate to miss the peak with the default settings, resulting in an inaccurate approximation.  Since the integrand is otherwise well behaved it appears to be almost always possible to overcome this limitation and produce calculations to extreme accuracy.  When the warning messages for NIntegrate are turned on, the first message will generally give a hint as to how to proceed.  If the goal is high precision the first strategy should be to increase the option, WorkingPrecision.  If machine precision is acceptable, but the problem is that the peak is missed (most common with PDF calculations) the best strategy is to increase the option, GaussPoints.  Even though the integrand does not contain a singularity, sometimes it will be necessary to increase the option, SingularityDepth to catch the peak.

OnNIntMsg;

SPDF[0.001, {1.3, 0}]

Here one could increase the recursion limit as suggested by the error message, but for these functions it is often a better idea to start with GaussPoints which sets a higher number of initial points, thus requiring less recursive sub-division.  The difference between the calculations is small.

SPDF[0.001, {1.3, 0},  GaussPoints ->  15]
% - %%

For extreme precision start with working precision.  As the output precision will be limited by the specification of the precision of the input values, it is simplest to use rational numbers.  GaussPoints and/or MaxRecursion will also likely need to be increased.   The default setting for PrecisionGoal in NIntegrate is WorkingPrecision minus 10 digits.

Table[
        SPDF[1/1000, {13/10, 0},  GaussPoints ->  30, MaxRecursion ->  10, 
          WorkingPrecision ->  p], {p, 16, 64, 4}] // TableForm

The package help browser gives additional examples to help handle error messages generated by the other functions.  It contains a section that demonstrates other methods of calculating stable distributions using Mathematica, and it includes the derivation of the Zolotarev integral.  There is also an algorithm to set up an interpolation method for approximating the density function for maximum likelihood calculations.

MathLink Implementation

The direct use of Mathematica for data analysis can be slow.  For in depth analysis of data there is an excellent MathLink package that includes many more functions and parameterizations. The MathLink program consists of compiled Fortran routines that are passed to Mathematica with wrapper functions that are compatible with the functions in the StableDistribution package.  It solves the problem of the shape of the Zolotarev integrand by dividing it into intervals.  The peak is found numerically and a narrow interval on either side of the peak is handled with a more involved routine, than the other parts.  It's complied code can thus use simpler and faster numerical integration routines for greater speed than possible with Mathematica's very versatile but complicated NIntegrate function.  It is also bullet-proof, at the cost of not allowing calculations when Abs[α-1] <0.01 but rounding to α = 1, and not permitting calculations when α <  0.1.  It is limited to double precision calculations.  The interface can perform maximum likelihood calculations many times faster than the Mathematica package described above, using a very fast interpolation method.  These fits are usually a minimum of 50 times faster than possible with the package described in this paper.  The MathLink interface has been used with Mathematica's JLink to download stock market data from internet data bases, and automatically analyze the time series of thousands of stocks.  The MathLink interface is available (currently for Windows only) for $200 from  http://www.robustanalysis.com.  The package described in this paper is designed so that both programs can be run simultaneously turning Mathematica into a powerful stable work station.

References

Chambers, J.M., Mallows, C.L. and Stuck, B.W. (1976). A Method for Simulating Stable Random Variables, Journal of the American Statistical Association (71), 1340-1344.

Levy, P. (1925). Calcul des Probabilites. Paris: Gauthier-Villars.

Mandelbrot, B. (1963). The Variation of Certain Speculative Prices, J.Business, 394-419.

Mandelbrot, B.(1997). Fractals and Scaling in Finance, Discontinuity Concentration, Risk,  Springer, 1997.

Nolan, J.P., American University Stable Web Site: http://academic2.american.edu/~jpnolan/stable/stable.html

Nolan, J.P.,  Stable MathLink Package: http://www.robustanalysis.com

Nolan, J.P.,  Stable Distributions, Models for Heavy Tailed Data, Birkhauser, (2005).
Chapter 1 provides a good introduction and is available on the American University Stable Web Site.

Nolan, J.P. (1997). Numerical calculation of stable densities and distribution functions. Commun. Statist. -Stochastic Models 13, 759–774.

Rachev, S. and Mittnik, S.(2000)  Stable Paretian Models in Finance. John Wiley & Sons, 2000.

Rimmer, R.H.,  Stable Probability Distribution Calculations: http://www.npgcable.com/~rrimmer/StableCalculation
An updated version of the Mathematica StableDistribution package will be maintained here.

Zolotarev, V.M. (1986), One-dimensional Stable Distributions, Amer. Math. Soc. Transl. of Math. Monographs, Vol. 65. Amer. Math. Soc., Providence, RI. (Transl. of the original 1983 Russian)