Special Cases of AlphaStable

[azev77]

According to Wiki:
For α = 2 the distribution reduces to a Gaussian distribution with variance σ2 = 2c2 and mean μ; the skewness parameter β has no effect.
For α = 1 and β = 0 the distribution reduces to a Cauchy distribution with scale parameter c and shift parameter μ.
For α = 1/2 and β = 1 the distribution reduces to a Lévy distribution with scale parameter c and shift parameter μ.

I get:

julia> using Distributions, AlphaStableDistributions
julia> d_train = rand(Normal(3,4), 100_000);
julia> fit(AlphaStable, d_train)
AlphaStable{Float64}(α=1.9914574535808538, β=0.0, scale=2.8151451718957867, location=2.9906397632987503)
julia> d_train = rand(Cauchy(3,4), 100_000);
julia> fit(AlphaStable, d_train)
AlphaStable{Float64}(α=0.9998014163853479, β=0.0, scale=4.0005882600294305, location=2.9901037116713165)
julia> d_train = rand(Levy(3,4), 100_000);
julia> fit(AlphaStable, d_train)
AlphaStable{Float64}(α=0.5845238495639795, β=0.0, scale=16.770767994549328, location=15.212645543834048)

There are some discrepancies. Not sure if this is due to a different parametrization or something else.

[azev77]

I resolved part of the issue at this link.
Updated code:

using Distributions, AlphaStableDistributions;
μ=3.0;σ=4.0;
#
#==#Normal(μ,σ) == AlphaStable(α=2.0, β=β, scale=σ/√2≈2.83, location=μ)
d_train = rand(Normal(μ,σ), 100_000);
D̂=fit(AlphaStable, d_train);
isapprox.((D̂.α,D̂.β,D̂.scale,D̂.location),(2.0,0.0,σ/√2,μ), atol=.25)  # All true, except β=0
#
#==#Cauchy(μ,σ) == AlphaStable(α=1.0, β=0.0, scale=σ, location=μ)
d_train = rand(Cauchy(μ,σ), 100_000);
D̂=fit(AlphaStable, d_train);
isapprox.((D̂.α,D̂.β,D̂.scale,D̂.location),(1.0,0.0,σ,μ), atol=.25)  # All true
#
#==#Levy(μ,σ) == AlphaStable(α=1/2, β=1.0, scale=σ, location=μ)
d_train = rand(Levy(μ,σ), 100_000);
D̂=fit(AlphaStable, d_train);
isapprox.((D̂.α,D̂.β,D̂.scale,D̂.location),(0.5,1.0,σ,μ), atol=.25) 
AlphaStable{Float64}(α=0.58, β=0.0, scale=16.55, location=15.15) # only α approx correct
# α=0.58, β=0.0, scale=16.55, location=15.15 

Looks like part of the issue is it currently fits symmetric (β=0.0) alpha stable distributions. Thus you have β=zero(α).

I wonder if this explains why the location & scale parameters are off for the Levy distribution?

[ymtoo]

You're right. The current implementation fits Symmetric Alpha Stable Distributions. The location and scale parameter estimations are based on Fama & Roll (1971) method which is restricted to the symmetrical case β=0, and α values in the range of [1, 2].

[mchitre]

@ymtoo it might make sense to have 2 distributions exported, one for the symmetric case, and one for the asymmetric case? The symmetric version will provide better estimates for parameters where symmetry can be inferred from the problem at hand.

[azev77]

@mchitre great idea.
Maybe call one: SymmetricAlphaStable()
Other: AlphaStable()

[azev77]

@ymtoo works great now except for Levy(), which I guess is bc its α=.5, but currently the package needs α>.6.
Hopefully, it can be expanded to α in [0,2].

[ymtoo]

@azev77 thanks for the feedback! We'll definitely consider the more general case if we come across any efficient methods.

[azev77]

What about Koutrouvelis (1980,1981)?
It's slightly slower than McCulloch (1986) but more accurate.

Koutrouvelis (1980)
"Regression-Type Estimation of the Paramters of Stable Laws. JASA, Vol 75, No. 372
Koutrouvelis (1981) "An Iterative Procedure for the estimation of the Parameters of Stable Laws"Commun. Stat. - Simul. Comput.
McCulloch (1986) "Simple Consistent Estimators of Stable Distribution Parameters" Cummun. Stat. Simul. Comput.

[ymtoo]

Thanks for the references. Will take a look when I've some free time. In the meantime, PRs are always welcome!