The positive stable distribution in R

Question

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?

Thank you in advance!

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

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

Answer 2

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.