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

 A: The short answer is that your $\delta$ is fine, but your $\gamma$ is wrong. In order to get the positive stable distribution given by your formula in R, you need to set
$$
\gamma = |1 - i \tan \left(\pi \alpha / 2\right)|^{-1/\alpha}.
$$
The earliest example I could find of the formula you gave was in (Feller, 1971), but I've only found that book in physical form. However (Hougaard, 1986) gives the same formula, along with the Laplace transform
$$
\mathrm{L}(s) = \mathrm{E}\left[\exp(-sX)\right] = \exp\left(-s^\alpha\right).
$$
From the stabledist manual (stabledist is used in fBasics), the pm=1 parameterization is from (Samorodnitsky and Taqqu, 1994), another resource whose online reproduction has eluded me. However (Weron, 2001) gives the characteristic function in Samorodnitsky and Taqqu's parameterization for $\alpha \neq 1$ to be
$$
\varphi(t) = \mathrm{E}\left[\exp(i t X) \right] = \exp\left[i \delta t - \gamma^\alpha |t|^\alpha \left(1 - i \beta \mathrm{sign}(t) \tan{\frac{\pi \alpha}{2}} \right) \right].
$$
I've renamed some parameters from Weron's paper to coinside with the notation we're using. He uses $\mu$ for $\delta$ and $\sigma$ for $\gamma$. In any case, plugging in $\beta=1$ and $\delta=0$, we get
$$
\varphi(t) = \exp\left[-\gamma^\alpha |t|^\alpha \left(1 - i \mathrm{sign}(t) \tan \frac{\pi \alpha}{2} \right) \right].
$$
Note that $(1 - i \tan (\pi \alpha / 2)) / |1 - i \tan(\pi \alpha / 2)| = \exp(-i \pi \alpha / 2)$ for $\alpha \in (0, 1)$ and that $i^\alpha = \exp(i \pi \alpha / 2)$. Formally, $\mathrm{L}(s)=\varphi(is)$, so by setting $\gamma = |1 - i \tan \left(\pi \alpha / 2\right)|^{-1/\alpha}$ in $\varphi(t)$ we get
$$
\varphi(is) = \exp\left(-s^\alpha\right) = \mathrm{L}(s).
$$
One interesting point to note is that the $\gamma$ that corresponds to $\alpha=1/2$ is also $1/2$, so if you were to try $\gamma=\alpha$ or $\gamma=1-\alpha$, which is actually not a bad approximation, you end up exactly correct for $\alpha=1/2$.
Here's an example in R to check correctness:
library(stabledist)

# Series representation of the density
PSf <- function(x, alpha, K) {
  k <- 1:K
  return(
    -1 / (pi * x) * sum(
      gamma(k * alpha + 1) / factorial(k) * 
        (-x ^ (-alpha)) ^ k * sin(alpha * k * pi)
    )
  )
}

# Derived expression for gamma
g <- function(a) {
  iu <- complex(real=0, imaginary=1)
  return(abs(1 - iu * tan(pi * a / 2)) ^ (-1 / a))
}

x=(1:100)/100
plot(0, xlim=c(0, 1), ylim=c(0, 2), pch='', 
     xlab='x', ylab='f(x)', main="Density Comparison")
legend('topright', legend=c('Series', 'gamma=g(alpha)'),
       lty=c(1, 2), col=c('gray', 'black'),
       lwd=c(5, 2))
text(x=c(0.1, 0.25, 0.7), y=c(1.4, 1.1, 0.7), 
     labels=c(expression(paste(alpha, " = 0.4")),
              expression(paste(alpha, " = 0.5")),
              expression(paste(alpha, " = 0.6"))))

for(a in seq(0.4, 0.6, by=0.1)) {
  y <- vapply(x, PSf, FUN.VALUE=1, alpha=a, K=100)
  lines(x, y, col="gray", lwd=5, lty=1)
  lines(x, dstable(x, alpha=a, beta=1, gamma=g(a), delta=0, pm=1), 
        col="black", lwd=2, lty=2)
}

$\hskip1in$



*

*Feller, W. (1971). An Introduction to Probability Theory and Its Applications, 2, 2nd ed. New York: Wiley.

*Hougaard, P. (1986). Survival Models for Heterogeneous Populations Derived from Stable Distributions, Biometrika 73, 387-396.

*Samorodnitsky, G., Taqqu, M.S. (1994). Stable Non-Gaussian Random Processes, Chapman & Hall, New York, 1994.

*Weron, R. (2001). Levy-stable distributions revisited: tail index > 2 does not exclude the Levy-stable regime, International Journal of Modern Physics C, 2001, 12(2), 209-223.

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