# The positive stable distribution in R

Positive stable distributions are described by four parameters: the skewness parameter $\beta\in[-1,1]$, the scale parameter $\sigma>0$, the location parameter $\mu\in(-\infty,\infty)$, and the so-called index parameter $\alpha\in(0,2]$. When $\beta$ is zero the distribution is symmetric around $\mu$, when it is positive (resp. negative) the distribution is skewed to the right (resp. to the left). Stable distributions allow fat tails when $\alpha$ decreases.

When $\alpha$ is strictly less than one and $\beta=1$ the support of the distribution restricts to $(\mu,\infty)$.

The density function only has a closed-form expression for some particular combinations of values for the parameters. When $\mu=0$, $\alpha<1$, $\beta=1$, and $\sigma=\alpha$ it is (see formula (4.4) here):

$f(y) = -\frac{1}{\pi y} \sum_{k=1}^{\infty} \frac{\Gamma(k\alpha+1)}{k!} (-y^{-\alpha})^k \sin(\alpha k \pi)$

It has infinite mean and variance.

Question

I would like to use that density in R. I use

> alpha <- ...
> dstable(y, alpha=alpha, beta=1, gamma=alpha, delta=0, pm=1)


where the dstable function comes with the fBasics package.

Can you confirm this is the right way to compute that density in R?

EDIT

One reason why I am suspicious is that, in the output, the value of delta is different from that in the input. Example:

> library(fBasics)
> alpha <- 0.4
> dstable(4, alpha=alpha, beta=1, gamma=alpha, delta=0, pm=1)
[1] 0.02700602
attr(,"control")
dist alpha beta gamma    delta pm
stable   0.4    1   0.4 0.290617  1

-

What I think is happening is that in the output delta may be reporting an internal location value, while in the input delta is describing the shift. [There seems to be a similar issue with gamma when pm=2.] So if you try increasing the shift to 2

> dstable(4, alpha=0.4, beta=1, gamma=0.4, delta=2, pm=1)
[1] 0.06569375
attr(,"control")
dist alpha beta gamma    delta pm
stable   0.4    1   0.4 2.290617  1


then you add 2 to the location value.

With beta=1 and pm=1 you have a positive random variable with a distribution lower bound at 0.

> min(rstable(100000, alpha=0.4, beta=1, gamma=0.4, delta=0, pm=1))
[1] 0.002666507


Shift by 2 and the lower bound rises by the same amount

> min(rstable(100000, alpha=0.4, beta=1, gamma=0.4, delta=2, pm=1))
[1] 2.003286


But if you want the delta input to be the internal location value rather than the shift or lower bound, then you need to use a different specification for the parameters. For example if you try the following (with pm=3 and trying delta=0 and the delta=0.290617 you found earlier), you seem to get the same delta in and out. With pm=3 and delta=0.290617 you get the same density of 0.02700602 you found earlier and a lower bound at 0. With pm=3 and delta=0 you get a negative lower bound (in fact -0.290617).

> dstable(4, alpha=0.4, beta=1, gamma=0.4, delta=0, pm=3)
[1] 0.02464434
attr(,"control")
dist alpha beta gamma delta pm
stable   0.4    1   0.4     0  3
> dstable(4, alpha=0.4, beta=1, gamma=0.4, delta=0.290617, pm=3)
[1] 0.02700602
attr(,"control")
dist alpha beta gamma    delta pm
stable   0.4    1   0.4 0.290617  3
> min(rstable(100000, alpha=0.4, beta=1, gamma=0.4, delta=0, pm=3))
[1] -0.2876658
> min(rstable(100000, alpha=0.4, beta=1, gamma=0.4, delta=0.290617, pm=3))
[1] 0.004303485


You may find it easier simply to ignore delta in the output, and so long as you keep beta=1 then using pm=1 means delta in the input is the distribution lower bound, which it seems you want to be 0.

-
 Thank you very much for these explanations/experiments! If I could vote twice, I would do! – ocram Mar 20 '11 at 16:36 +1 for this lengthy one, deleted mine wrong guess with negative default pm = 0 parametrization. – Dmitrij Celov Mar 20 '11 at 17:40

Also of note: Martin Maechler just refactored the code for the stable distributed and added some improvements.

His new package stabledist will be used by fBasics as well, so you may want to give this a look as well.

-
 Yes, I will have a look at it! Thx!! – ocram Mar 21 '11 at 5:59