An implementation of an alpha stable distribution in Julia

AlphaStable.jl

using StatsBase, Distributions

Base.@kwdef struct AlphaStable{T} <: Distributions.ContinuousUnivariateDistribution
    α::T = 1.5
    β::T = 0.0
    scale::T = 1.0
    location::T = 0.0
end


# sampler(d::AlphaStable) = error("Not implemented")
# pdf(d::AlphaStable, x::Real) = error("Not implemented")
# logpdf(d::AlphaStable, x::Real) = error("Not implemented")
# cdf(d::AlphaStable, x::Real) = error("Not implemented")
# quantile(d::AlphaStable, q::Real) = error("Not implemented")
# minimum(d::AlphaStable) = error("Not implemented")
# maximum(d::AlphaStable) = error("Not implemented")
# insupport(d::AlphaStable, x::Real) = error("Not implemented")
Statistics.mean(d::AlphaStable) = d.α > 1 ? d.location : error("Not defined")
Statistics.var(d::AlphaStable) = d.α == 2 ? 2d.scale^2 : Inf
# modes(d::AlphaStable) = error("Not implemented")
# mode(d::AlphaStable) = error("Not implemented")
# skewness(d::AlphaStable) = error("Not implemented")
# kurtosis(d::Distribution, ::Bool) = error("Not implemented")
# entropy(d::AlphaStable, ::Real) = error("Not implemented")
# mgf(d::AlphaStable, ::Any) = error("Not implemented")
# cf(d::AlphaStable, ::Any) = error("Not implemented")

# lookup table from McCulloch (1986)
const _ena = [
           2.4388
           2.5120
           2.6080
           2.7369
           2.9115
           3.1480
           3.4635
           3.8824
           4.4468
           5.2172
           6.3140
           7.9098
          10.4480
          14.8378
          23.4831
          44.2813
]

"""
Fit a symmetric α stable distribution to data.

:param x: data
:returns: (α, c, δ)

α, c and δ are the characteristic exponent, scale parameter
(dispersion^1/α) and location parameter respectively.

α is computed based on McCulloch (1986) fractile.
c is computed based on Fama & Roll (1971) fractile.
δ is the 50% trimmed mean of the sample.
"""
function fit(d::Type{<:AlphaStable}, x)
    δ = mean(StatsBase.trim(x,prop=0.25))
    p = quantile.(Ref(sort(x)), (0.05, 0.25, 0.28, 0.72, 0.75, 0.95), sorted=true)
    c = (p[4]-p[3])/1.654
    an = (p[6]-p[1])/(p[5]-p[2])
    if an < 2.4388
        α = 2.
    else
        α = 0.
        j = findfirst(>=(an), _ena) # _np.where(an <= _ena[:,0])[0]
        if j !== nothing
            t = (an-_ena[j,1])/(_ena[j+1]-_ena[j])
            α = (21-j-t)/10
        end
    end
    if α < 0.5
        α = 0.5
    end
    return α, c, δ
end

"""
Generate independent stable random numbers.

:param α: characteristic exponent (0.1 to 2.0)
:param β: skew (-1 to +1)
:param scale: scale parameter
:param loc: location parameter (mean for α > 1, median/mode when β=0)


This implementation is based on the method in J.M. Chambers, C.L. Mallows
and B.W. Stuck, "A Method for Simulating Stable Random Variables," JASA 71 (1976): 340-4.
McCulloch's MATLAB implementation (1996) served as a reference in developing this code.
"""
function Base.rand(rng::AbstractRNG, d::AlphaStable)
    α=d.α; β=d.β; scale=d.scale; loc=d.location
    (α < 0.1 || α > 2) && throw(DomainError(α, "α must be in the range 0.1 to 2"))
    abs(β) > 1 && throw(DomainError(β, "β must be in the range -1 to 1"))
    ϕ = (rand(rng) - 0.5) * π
    if α == 1 && β == 0
        return loc + scale*tan(ϕ)
    end
    w = -log(rand(rng))
    α == 2 && (return loc + 2*scale*sqrt(w)*sin(ϕ))
    β == 0 && (return loc + scale * ((cos((1-α)*ϕ) / w)^(1.0/α - 1) * sin(α * ϕ) / cos(ϕ)^(1.0/α)))
    cosϕ = cos(ϕ)
    if abs(α-1) > 1e-8
        ζ = β * tan(π*α/2)
        aϕ = α * ϕ
        a1ϕ = (1-α) * ϕ
        return loc + scale * (( (sin(aϕ)+ζ*cos(aϕ))/cosϕ * ((cos(a1ϕ)+ζ*sin(a1ϕ))) / ((w*cosϕ)^((1-α)/α)) ))
    end
    bϕ = π/2 + β*ϕ
    x = 2/π * (bϕ*tan(ϕ) - β*log(π/2*w*cosϕ/bϕ))
    α == 1 || (x += β * tan(π*α/2))

    return loc + scale*x
end

# @btime rand(AlphaStable())
#
# s = [rand(AlphaStable()) for _ in 1:100000]
# histogram(s)
# @btime fit(AlphaStable, $s)