Fitting the parameters of a stable distribution

Question

I have a data set and I have to fit this data set with a stable distribution. The problem is that the stable distributions are known analytically only in the form of the characteristic function (Fourier transform). How can I do this?

Answer 1

As suggested in the comments, you can use fitdistr, with the density function from fBasics.

# Sample data
x <- rt(100,df=4)

# Density (I reparametrize it to remove the constraints 
# on the parameters)
library(fBasics)
f <- function(u,a,b,c,d) {
  cat(a,b,c,d,"\n")  # Some logging (it is very slow)
  dstable(u, 2*exp(a)/(1+exp(a)), 2*exp(b)/(1+exp(b))-1, exp(c), d)
}

# Fit the distribution
library(MASS)
r <- fitdistr(x, f, list(a=1, b=0, c=log(mad(x)), d=median(x)))
r

# Graphical check
plot(
  qstable(ppoints(100),
    2*exp(r$estimate[1])/(1+exp(r$estimate[1])),
    2*exp(r$estimate[2])/(1+exp(r$estimate[2]))-1,
    exp(r$estimate[3]),
    r$estimate[4]
  ),
  sort(x)
)
abline(0,1)

Answer 2

One way to fit the $\alpha$ parameter is via the Nagaev transform described by Okoneshnikov. An alternative is the 'Probability of Return' method of Mantegna and Stanley, which is considerably easier.

edit: the other 'classical' method is of Kogon & Williams (S.M. Kogon, Douglas B. Williams, "On Characteristic Function Based Stable Distribution Parameter Estimation Techniques"), see also matlab implementation of K&W

Answer 3

@Vincent's answer sounds good, but here is another approach: Since you know the Fourier transform of the distribution, take the appropriate Fourier transformation of the data, and find parameters that give the best fit in Fourier space.

I think this method should work just as well in theory, and in practice would avoid lots of numerical integration to get the form of the stable distributions. I am not coding up the test now, sorry. Anyone have any insight on this?