For the efficient generation of random variates, we use the following useful fact (see e.g. Theorem 19.1 in Haubold, Mathai, and Saxena (2011)): A standard α-Mittag-Leffler random variable Y has the representation:
$$Y \stackrel{d}{=} X^{1/\alpha} Z$$
where X is standard exponentially distributed, Z is α-stable with Laplace Transform E[exp (−sZ)] = exp (−sα), X and Z are independent, and $\stackrel{d}{=}$ means equality in distribution.
To generate such random variates Z, we use
a <- 0.8
sigma <- (cos(pi*a/2))^(1/a)
z <- stabledist::rstable(n = n, alpha = a, beta = 1, gamma = sigma, delta = 0, pm = 1)Below are the details of the calculation. We use the parametrization of the stable distribution by Samorodnitsky and Taqqu (1994) as it has become standard. For α ∈ (0, 1) and α ∈ (1, 2),
$$\mathbf E[\exp(it Z)] = \exp\left\lbrace -\sigma^\alpha |t|^\alpha \left[1 - i \beta {\rm sgn}t \tan \frac{\pi \alpha}{2}\right] + i a t\right\rbrace$$
As in Meerschaert and Scheffler (2001), Equation (7.28), set
$$\sigma^\alpha = C \Gamma(1-\alpha) \cos \frac{\pi\alpha}{2},$$
for some constant C > 0, set β = 1, set a = 0, and the log-characteristic function becomes
Setting t = is recovers the Laplace transform, and to match the Laplace transform exp (−sα) of Z, it is necessary that CΓ(1 − α) = 1. But then σα = cos (πα/2), and we see that
Z ∼ S(α, β, σ, a) = S(α, 1, cos (πα/2)1/α, 0)