Skip to content

Commit 220eab8

Browse files
committed
Add SymmetricAlphaStable
1 parent fb74898 commit 220eab8

File tree

1 file changed

+61
-56
lines changed

1 file changed

+61
-56
lines changed

src/AlphaStableDistributions.jl

Lines changed: 61 additions & 56 deletions
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,7 @@ using StatsBase, Distributions, StaticArrays
55
using MAT, SpecialFunctions, ToeplitzMatrices
66
using Interpolations
77

8-
export AlphaStable, AlphaSubGaussian, fit
8+
export AlphaStable, SymmetricAlphaStable, AlphaSubGaussian, fit
99

1010
Base.@kwdef struct AlphaStable{T} <: Distributions.ContinuousUnivariateDistribution
1111
α::T = 1.5
@@ -33,27 +33,25 @@ Statistics.var(d::AlphaStable) = d.α == 2 ? 2d.scale^2 : Inf
3333
# mgf(d::AlphaStable, ::Any) = error("Not implemented")
3434
# cf(d::AlphaStable, ::Any) = error("Not implemented")
3535

36-
# # lookup table from McCulloch (1986)
37-
# const _ena = [
38-
# 2.4388
39-
# 2.5120
40-
# 2.6080
41-
# 2.7369
42-
# 2.9115
43-
# 3.1480
44-
# 3.4635
45-
# 3.8824
46-
# 4.4468
47-
# 5.2172
48-
# 6.3140
49-
# 7.9098
50-
# 10.4480
51-
# 14.8378
52-
# 23.4831
53-
# 44.2813
54-
# ]
55-
5636
# lookup tables from McCulloch (1986)
37+
const _ena = [
38+
2.4388
39+
2.5120
40+
2.6080
41+
2.7369
42+
2.9115
43+
3.1480
44+
3.4635
45+
3.8824
46+
4.4468
47+
5.2172
48+
6.3140
49+
7.9098
50+
10.4480
51+
14.8378
52+
23.4831
53+
44.2813
54+
]
5755
const _να = [
5856
2.439
5957
2.500
@@ -183,8 +181,8 @@ Fit an α stable distribution to data.
183181
184182
returns `AlphaStable`
185183
186-
α, β, c and δ are the characteristic exponent, skewness parameter, scale parameter
187-
(dispersion^1/α) and location parameter respectively.
184+
α∈[0.6,2.0], β∈[-1,1] , c∈[0,∞] and δ∈[-∞,∞] are the characteristic exponent,
185+
skewness parameter, scale parameter (dispersion^1/α) and location parameter respectively.
188186
189187
α, β, c and δ are computed based on McCulloch (1986) fractile.
190188
"""
@@ -214,39 +212,46 @@ function Distributions.fit(::Type{<:AlphaStable}, x)
214212
return AlphaStable=α, β=β, scale=c, location=oftype(α, δ))
215213
end
216214

217-
# """
218-
# Fit a symmetric α stable distribution to data.
219-
220-
# :param x: data
221-
# :returns: (α, c, δ)
222-
223-
# α, c and δ are the characteristic exponent, scale parameter
224-
# (dispersion^1/α) and location parameter respectively.
225-
226-
# α is computed based on McCulloch (1986) fractile.
227-
# c is computed based on Fama & Roll (1971) fractile.
228-
# δ is the 50% trimmed mean of the sample.
229-
# """
230-
# function Distributions.fit(d::Type{<:AlphaStable}, x)
231-
# δ = mean(StatsBase.trim(x,prop=0.25))
232-
# p = quantile.(Ref(sort(x)), (0.05, 0.25, 0.28, 0.72, 0.75, 0.95), sorted=true)
233-
# c = (p[4]-p[3])/1.654
234-
# an = (p[6]-p[1])/(p[5]-p[2])
235-
# if an < 2.4388
236-
# α = 2.
237-
# else
238-
# α = 0.
239-
# j = findfirst(>=(an), _ena) # _np.where(an <= _ena[:,0])[0]
240-
# (j === nothing || j == length(_ena)) && (j = length(_ena))
241-
# t = (an-_ena[j-1])/(_ena[j]-_ena[j-1])
242-
# α = (22-j-t)/10
243-
244-
# end
245-
# if α < 0.5
246-
# α = 0.5
247-
# end
248-
# return AlphaStable(α=α, β=zero(α), scale=c, location=oftype(α, δ))
249-
# end
215+
Base.@kwdef struct SymmetricAlphaStable{T} <: Distributions.ContinuousUnivariateDistribution
216+
α::T = 1.5
217+
scale::T = one(α)
218+
location::T = zero(α)
219+
end
220+
221+
"""
222+
fit(d::Type{<:SymmetricAlphaStable}, x; alg=QuickSort)
223+
224+
Fit a symmetric α stable distribution to data.
225+
226+
returns `SymmetricAlphaStable`
227+
228+
α∈[1,2], c∈[0,∞] and δ∈[-∞,∞] are the characteristic exponent, scale parameter
229+
(dispersion^1/α) and location parameter respectively.
230+
231+
α is computed based on McCulloch (1986) fractile.
232+
scale is computed based on Fama & Roll (1971) fractile.
233+
location is the 50% trimmed mean of the sample.
234+
"""
235+
function Distributions.fit(::Type{<:SymmetricAlphaStable}, x)
236+
δ = mean(StatsBase.trim(x,prop=0.25))
237+
p = quantile.(Ref(sort(x)), (0.05, 0.25, 0.28, 0.72, 0.75, 0.95), sorted=true)
238+
c = (p[4]-p[3])/1.654
239+
an = (p[6]-p[1])/(p[5]-p[2])
240+
if an < 2.4388
241+
α = 2.
242+
else
243+
α = 0.
244+
j = findfirst(>=(an), _ena) # _np.where(an <= _ena[:,0])[0]
245+
(j === nothing || j == length(_ena)) && (j = length(_ena))
246+
t = (an-_ena[j-1])/(_ena[j]-_ena[j-1])
247+
α = (22-j-t)/10
248+
249+
end
250+
if α < 0.5
251+
α = 0.5
252+
end
253+
return AlphaStable=α, β=zero(α), scale=c, location=oftype(α, δ))
254+
end
250255

251256
"""
252257
Generate independent stable random numbers.

0 commit comments

Comments
 (0)