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 [1]. Stable distributions have been proposed as models for many types of physical and economic systems. There are solid theoretical reasons for expecting a non-Gaussian stable model, for example, reflection off a rotating mirror yielding a Cauchy distribution, hitting times for a Brownian motion yielding a Levy distribution, and the gravitational field of stars yielding the Holtsmark distribution. The generalized central limit theorem states that the only possible nontrivial 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 [2].

The lack of formulae 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 [3, 4, 5]. Rose and Smith [6] have 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. Our 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.