Description Usage Arguments Details Value Tail Behavior Note Author(s) References See Also Examples
A collection and description of functions to compute density, distribution and quantile function and to generate random variates of the stable distribution.
The four functions are:
[dpqr]stable  the (skewed) stable distribution. 
1 2 3 4 5 6 7 8 9 10 11 12  dstable(x, alpha, beta, gamma = 1, delta = 0, pm = 0,
log = FALSE,
tol = 64*.Machine$double.eps, zeta.tol = NULL,
subdivisions = 1000)
pstable(q, alpha, beta, gamma = 1, delta = 0, pm = 0,
lower.tail = TRUE, log.p = FALSE, silent = FALSE,
tol = 64*.Machine$double.eps, subdivisions = 1000)
qstable(p, alpha, beta, gamma = 1, delta = 0, pm = 0,
lower.tail = TRUE, log.p = FALSE,
tol = .Machine$double.eps^0.25, maxiter = 1000, trace = 0,
integ.tol = 1e7, subdivisions = 200)
rstable(n, alpha, beta, gamma = 1, delta = 0, pm = 0)

alpha, beta, gamma, delta 
value of the index parameter 
n 
sample size (integer). 
p 
numeric vector of probabilities. 
pm 
parameterization, an integer in 
x, q 
numeric vector of quantiles. 
log, log.p 
logical; if TRUE, probabilities p are given as log(p). 
lower.tail 
logical; if TRUE (default), probabilities are P[X ≤ x] otherwise, P[X > x]. 
silent 
logical indicating that e.g., warnings should be
suppressed when 
integ.tol 
positive number, the tolerance used for numerical
integration, see 
tol 
numerical tolerance,

zeta.tol 
( 
subdivisions 
maximal number of intervals for integration, see

maxiter, trace 
maximal number of iterations and verboseness in

Skew Stable Distribution:
The function uses the approach of J.P. Nolan for general stable
distributions. Nolan (1997) derived expressions in form of integrals
based on the characteristic function for standardized stable random
variables. For dstable
and pstable
, these integrals
are numerically evaluated using R's integrate()
function.
“S0” parameterization [pm=0]: based on the (M) representation
of Zolotarev for an alpha stable distribution with skewness
beta. Unlike the Zolotarev (M) parameterization, gamma and
delta are straightforward scale and shift parameters. This
representation is continuous in all 4 parameters, and gives
an intuitive meaning to gamma and delta that is lacking in
other parameterizations.
Switching the sign of beta
mirrors the distribution at
the vertical axis x = delta, i.e.,
f(x, α, β, γ, δ, 0) = f(2δx, α, +β, γ, δ, 0),
see the graphical example below.
“S” or “S1” parameterization [pm=1]: the parameterization used
by Samorodnitsky and Taqqu in the book Stable NonGaussian
Random Processes. It is a slight modification of Zolotarev's
(A) parameterization.
“S*” or “S2” parameterization [pm=2]: a modification of the S0
parameterization which is defined so that (i) the scale gamma
agrees with the Gaussian scale (standard dev.) when alpha=2
and the Cauchy scale when alpha=1, (ii) the mode is exactly at
delta. For this parametrization,
stableMode(alpha,beta)
is needed.
“S3” parameterization [pm=3]: an internal parameterization,
currently not available for these functions. The scale is the same
as the “S2” parameterization, the shift is
β*g(α), where g(α) is defined in Nolan(1999).
All values for the *stable
functions
are numeric vectors:
d*
returns the density,
p*
returns the distribution function,
q*
returns the quantile function, and
r*
generates random deviates.
The asymptotic behavior for large x, aka “tail behavior” for the cumulative F(x) = P(X <= x) is (for x > Inf)
1  F(x) ~ (1+b) C_a x^a,
where a=alpha, b=beta, C_a = Gamma(a)/pi * sin(a*pi/2); hence also
F(x) ~ (1+b) C_a x^a.
Differentiating F() above gives
f(x) ~ a(1+b) C_a x^(1+a).
In the case β = 1, the distributions are “maximally
skewed to the right” or simply “extremal stable”
(Zolotarev). In that case, the package FMStable provides
dpq*
functions which are faster and more accurate than ours (if
accuracy higher than about 6 digits is needed), see,
pEstable
.
When alpha
is close to 1 or close to 0 (“close”,
e.g., meaning distance d < 0.01), the computations typically are
numerically considerably more challenging, and the results may not be
accurate.
As we
plan to improve on this, and
as it is unknown when exactly the numerical difficulties arise, we
currently only do warn here (in the documentation), but not by
giving explicit warning()
s.
Diethelm Wuertz for the original Rmetrics Rport. Many numerical improvements by Martin Maechler.
Chambers J.M., Mallows, C.L. and Stuck, B.W. (1976) A Method for Simulating Stable Random Variables, J. Amer. Statist. Assoc. 71, 340–344.
John P. Nolan (2012) Stable Distributions  Models for Heavy Tailed Data Birkhauser, Boston; in progress, chapter 1 online at http://academic2.american.edu/~jpnolan/stable/chap1.pdf
Nolan J.P. (1997)
Numerical calculation of stable densities and distribution functions.
Stochastic Models 13(4), 759–774.
Also available as ‘density.ps’ from Nolan's web page.
Samoridnitsky G., Taqqu M.S. (1994); Stable NonGaussian Random Processes, Stochastic Models with Infinite Variance, Chapman and Hall, New York, 632 pages.
Weron, A., Weron R. (1999); Computer Simulation of Levy alphaStable Variables and Processes, Preprint Technical University of Wroclaw, 13 pages.
the stableSlider()
function from package
fBasics for displaying densities and probabilities of these
distributions, for educational purposes.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37  ## stable 
## Plot stable random number series
set.seed(1953)
r < rstable(n = 1000, alpha = 1.9, beta = 0.3)
plot(r, type = "l", main = "stable: alpha=1.9 beta=0.3",
col = "steelblue")
grid()
## Plot empirical density and compare with true density:
hist(r, n = 25, probability = TRUE, border = "white",
col = "steelblue")
x < seq(5, 5, 0.25)
lines(x, dstable(x, alpha = 1.9, beta = 0.3, tol= 1e3), lwd = 2)
## Plot df and compare with true df:
plot(ecdf(r), do.points=TRUE, col = "steelblue",
main = "Probabilities: ecdf(rstable(1000,..)) and true cdf F()")
rug(r)
lines(x, pstable(q = x, alpha = 1.9, beta = 0.3),
col="#0000FF88", lwd= 2.5)
## Switching sign(beta) <==> Mirror the distribution around x == delta:
curve(dstable(x, alpha=1.2, beta = .8, gamma = 3, delta = 2), 10, 10)
curve(dstable(x, alpha=1.2, beta = .8, gamma = 3, delta = 2),
add=TRUE, col=2)
## or the same
curve(dstable(2*2x, alpha=1.2, beta = +.8, gamma = 3, delta = 2),
add=TRUE, col=adjustcolor("gray",0.2), lwd=5)
abline(v = 2, col = "gray", lty=2, lwd=2)
axis(1, at = 2, label = expression(delta == 2))
## Compute quantiles:
x. < 4:4
px < pstable(x., alpha = 1.9, beta = 0.3)
(qs < qstable(px, alpha = 1.9, beta = 0.3))
stopifnot(all.equal(as.vector(qs), x., tol = 1e5))

[1] 4.000020e+00 3.000005e+00 2.000000e+00 1.000002e+00 1.173636e08
[6] 1.000000e+00 1.999998e+00 3.000002e+00 4.000000e+00
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.