StableDistribution: Stable Distribution Function

Description

A collection and description of functions to compute density, distribution and quantile function and to generate random variates of the stable distribution.

The functions are: ll{ [dpqr]stable the (skewed) stable distribution. }

Usage

dstable(x, alpha, beta, gamma = 1, delta = 0, pm = 0,
        log = FALSE,
        tol = 64*.Machine$double.eps, zeta.tol = 1e-5,
        subdivisions = 1000)
pstable(q, alpha, beta, gamma = 1, delta = 0, pm = 0,
        lower.tail = TRUE, log.p = 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,
        integ.tol = 1e-7, subdivisions = 200)
rstable(n, alpha, beta, gamma = 1, delta = 0, pm = 0)

Arguments

alpha, beta, gamma, delta

value of the index parameter alpha with alpha = (0,2]; skewness parameter beta, in the range [-1, 1]; scale parameter gamma; and shift parameter delta.

sample size (integer).

numeric vector of probabilities.

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 \le x]$ otherwise, $P[X > x]$.

integ.tol

positive number, the tolerance used for numerical integration, see integrate.

tol

numerical tolerance, [object Object],[object Object]

zeta.tol

(dstable) numerical tolerance for checking if x is close to $\zeta(\alpha,\beta)$. As this is experimental and may well go away in the future, its use is not recommended in production code. Rather e-mail th

subdivisions

maximal number of intervals for integration, see integrate.

maxiter

maximal number of iterations in uniroot.

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.

concept

Stable Distribution

Tail Behavior

The asymptotic behavior for large $x$, aka tail behavior for the cumulative $F(x) = P(X \le x)$ is (for $x\to\infty$) $$1 - F(x) \sim (1+\beta) C_\alpha x^{-\alpha},$$ where $C_\alpha = \Gamma(\alpha)/\pi \sin(\alpha\pi/2)$; hence also $$F(-x) \sim (1-\beta) C_\alpha x^{-\alpha}.$$

Differentiating $F()$ above gives $$f(x) \sim \alpha(1+\beta) C_\alpha x^{-(1+\alpha)}.$$

Details

Skew Stable Distribution: The function uses the approach of J.P. Nolan for general stable distributions. Nolan derived expressions in form of integrals based on the characteristic function for standardized stable random variables. These integrals are numerically evaluated using R's integrate() function, following Nolan (1997). 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, \alpha, -\beta, \gamma, \delta, 0) = f(2\delta-x, \alpha, +\beta, \gamma, \delta, 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. 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 $-beta*g(alpha)$, where $g(alpha)$ is defined in Nolan(1999).

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. (1999); Stable Distributions, Preprint, University Washington DC, 30 pages.

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 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.

Examples

Run this code

## 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))

Run the code above in your browser using DataLab