A collection and description of functions to compute
density, distribution and quantile function and
to generate random variates of the stable distribution.
value of the index parameter alpha with alpha = (0,2];
skewness parameter beta, in the range [-1, 1];
scale parameter gamma; and location (or ‘shift’)
parameter delta.
n
sample size (integer).
p
numeric vector of probabilities.
pm
parameterization, an integer in 0, 1, 2; by default pm=0,
the ‘S0’ parameterization.
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].
integ.tol
positive number, the tolerance used for numerical
integration, see integrate.
tol
numerical tolerance,
dstable(), pstable():
used for numerical integration, see
integ.tol above. Note that earlier versions had tighter
tolerances – which seem too tight as default values.
qstable():
used for rootfinding, see uniroot.
zeta.tol
(dstable) numerical tolerance for checking if
x is close to ζ(α,β). The default,
NULL depends itself on (α,β).
As it is experimental and not guaranteed to remain in the
future, its use is not recommended in production code. Rather e-mail the
package maintainer about it.
subdivisions
maximal number of intervals for integration, see
integrate.
maxiter, trace
maximal number of iterations and verboseness in
uniroot, see there.
Details
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 betamirrors 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 Non-Gaussian
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).
Value
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.
Tail Behavior
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).
Note
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.
Author(s)
Diethelm Wuertz for the original Rmetrics R-port.
Many numerical improvements by Martin Maechler.
References
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.
Nolan J.P. (1997)
Numerical calculation of stable densities and distribution functions.
Stochastic Models13(4), 759–774.
Also available as ‘density.ps’ from Nolan's web page.
Samoridnitsky G., Taqqu M.S. (1994);
Stable Non-Gaussian Random Processes, Stochastic Models
with Infinite Variance,
Chapman and Hall, New York, 632 pages.
Weron, A., Weron R. (1999);
Computer Simulation of Levy alpha-Stable Variables and
Processes,
Preprint Technical University of Wroclaw, 13 pages.
See Also
the stableSlider() function from package
fBasics for displaying densities and probabilities of these
distributions, for educational purposes.
Examples
## 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= 1e-3), 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*2-x, 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 = 1e-5))