From bd3e84508ccb34722366d093a8c76b28d6f33d02 Mon Sep 17 00:00:00 2001 From: ymtoo Date: Wed, 14 Oct 2020 00:39:10 +0800 Subject: [PATCH 1/9] Implement McCulloch's quantile method --- Project.toml | 1 + src/AlphaStableDistributions.jl | 254 ++++++++++++++++++++++++++------ 2 files changed, 207 insertions(+), 48 deletions(-) diff --git a/Project.toml b/Project.toml index 6283855..70f08c6 100644 --- a/Project.toml +++ b/Project.toml @@ -5,6 +5,7 @@ version = "1.0.0" [deps] Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MAT = "23992714-dd62-5051-b70f-ba57cb901cac" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" diff --git a/src/AlphaStableDistributions.jl b/src/AlphaStableDistributions.jl index 26c75ae..1948c67 100644 --- a/src/AlphaStableDistributions.jl +++ b/src/AlphaStableDistributions.jl @@ -3,7 +3,7 @@ module AlphaStableDistributions using LinearAlgebra, Statistics, Random using StatsBase, Distributions, StaticArrays using MAT, SpecialFunctions, ToeplitzMatrices - +using Interpolations export AlphaStable, AlphaSubGaussian, fit @@ -33,62 +33,220 @@ Statistics.var(d::AlphaStable) = d.α == 2 ? 2d.scale^2 : Inf # 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 +# # 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 +# ] + +# lookup tables from McCulloch (1986) +const _να = [ + 2.439 + 2.500 + 2.600 + 2.700 + 2.800 + 3.000 + 3.200 + 3.500 + 4.000 + 5.000 + 6.000 + 8.000 + 10.000 + 15.000 + 25.000 +] +const _νβ = [ + 0.0 + 0.1 + 0.2 + 0.3 + 0.5 + 0.7 + 1.0 +] +const ψ₁ = [ + 2.000 2.000 2.000 2.000 2.000 2.000 2.000 + 1.916 1.924 1.924 1.924 1.924 1.924 1.924 + 1.808 1.813 1.829 1.829 1.829 1.829 1.829 + 1.729 1.730 1.737 1.745 1.745 1.745 1.745 + 1.664 1.663 1.663 1.668 1.676 1.676 1.676 + 1.563 1.560 1.553 1.548 1.547 1.547 1.547 + 1.484 1.480 1.471 1.460 1.448 1.438 1.438 + 1.391 1.386 1.378 1.364 1.337 1.318 1.318 + 1.279 1.273 1.266 1.250 1.210 1.184 1.150 + 1.128 1.121 1.114 1.101 1.067 1.027 0.973 + 1.029 1.021 1.014 1.004 0.974 0.935 0.874 + 0.896 0.892 0.887 0.883 0.855 0.823 0.769 + 0.818 0.812 0.806 0.801 0.780 0.756 0.691 + 0.698 0.695 0.692 0.689 0.676 0.656 0.595 + 0.593 0.590 0.588 0.586 0.579 0.563 0.513 +] +const ψ₂ = [ + 0.0 2.160 1.000 1.000 1.000 1.000 1.000 + 0.0 1.592 3.390 1.000 1.000 1.000 1.000 + 0.0 0.759 1.800 1.000 1.000 1.000 1.000 + 0.0 0.482 1.048 1.694 2.229 1.000 1.000 + 0.0 0.360 0.760 1.232 2.229 1.000 1.000 + 0.0 0.253 0.518 0.823 1.575 1.000 1.000 + 0.0 0.203 0.410 0.632 1.244 1.906 1.000 + 0.0 0.165 0.332 0.499 0.943 1.560 1.000 + 0.0 0.136 0.271 0.404 0.689 1.230 2.195 + 0.0 0.109 0.216 0.323 0.539 0.827 1.917 + 0.0 0.096 0.190 0.284 0.472 0.693 1.759 + 0.0 0.082 0.163 0.243 0.412 0.601 1.596 + 0.0 0.074 0.147 0.220 0.377 0.546 1.482 + 0.0 0.064 0.128 0.191 0.330 0.478 1.362 + 0.0 0.056 0.112 0.167 0.285 0.428 1.274 +] +const _α =[ + 0.5 + 0.6 + 0.7 + 0.8 + 0.9 + 1.0 + 1.1 + 1.2 + 1.3 + 1.4 + 1.5 + 1.6 + 1.7 + 1.8 + 1.9 + 2.0 +] +const _β = [ + 0.0 + 0.25 + 0.50 + 0.75 + 1.00 +] +const ϕ₃ = [ + 2.588 3.073 4.534 6.636 9.144 + 2.337 2.635 3.542 4.808 6.247 + 2.189 2.392 3.004 3.844 4.775 + 2.098 2.244 2.676 3.265 3.912 + 2.040 2.149 2.461 2.886 3.356 + 2.000 2.085 2.311 2.624 2.973 + 1.980 2.040 2.205 2.435 2.696 + 1.965 2.007 2.125 2.294 2.491 + 1.955 1.984 2.067 2.188 2.333 + 1.946 1.967 2.022 2.106 2.211 + 1.939 1.952 1.988 2.045 2.116 + 1.933 1.940 1.962 1.997 2.043 + 1.927 1.930 1.943 1.961 1.987 + 1.921 1.922 1.927 1.936 1.947 + 1.914 1.915 1.916 1.918 1.921 + 1.908 1.908 1.908 1.908 1.908 +] +const ϕ₅ = [ + 0.0 0.0 0.0 0.0 0.0 + 0.0 -0.017 -0.032 -0.049 -0.064 + 0.0 -0.030 -0.061 -0.092 -0.123 + 0.0 -0.043 -0.088 -0.132 -0.179 + 0.0 -0.056 -0.111 -0.170 -0.232 + 0.0 -0.066 -0.134 -0.206 -0.283 + 0.0 -0.075 -0.154 -0.241 -0.335 + 0.0 -0.084 -0.173 -0.276 -0.390 + 0.0 -0.090 -0.192 -0.310 -0.447 + 0.0 -0.095 -0.208 -0.346 -0.508 + 0.0 -0.098 -0.223 -0.383 -0.576 + 0.0 -0.099 -0.237 -0.424 -0.652 + 0.0 -0.096 -0.250 -0.469 -0.742 + 0.0 -0.089 -0.262 -0.520 -0.853 + 0.0 -0.078 -0.272 -0.581 -0.997 + 0.0 -0.061 -0.279 -0.659 -1.198 ] """ - fit(d::Type{<:AlphaStable}, x; alg=QuickSort) - -Fit a symmetric α stable distribution to data. +Fit an α stable distribution to data. -returns `AlphaStable` +:param x: data +:returns: (α, β, c, δ) -α, c and δ are the characteristic exponent, scale parameter +α, β, c and δ are the characteristic exponent, skewness parameter, 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. +α, β, c and δ are computed based on McCulloch (1986) fractile. """ -function Distributions.fit(d::Type{<:AlphaStable}, x; alg=QuickSort) - sx = sort(x, alg=alg) - δ = mean(@view(sx[end÷4:(3*end)÷4])) - p = quantile.(Ref(sx), (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. +function Distributions.fit(::Type{<:AlphaStable}, x) + p = quantile.(Ref(sort(x)), (0.05, 0.25, 0.28, 0.5, 0.72, 0.75, 0.95), sorted=true) + να = (p[7]-p[1]) / (p[6]-p[2]) + νβ = (p[7]+p[1]-2p[4]) / (p[7]-p[1]) + (να < _να[1]) && (να = _να[1]) + (να > _να[end]) && (να = _να[end]) + + itp₁ = interpolate((_να, _νβ), ψ₁, Gridded(Linear())) + α = itp₁(να, abs(νβ)) + itp₂ = interpolate((_να, _νβ), ψ₂, Gridded(Linear())) + β = sign(νβ) * itp₂(να, abs(νβ)) + + (β > 1.0) && (β = 1.0) + (β < -1.0) && (β = -1.0) + itp₃ = interpolate((_α, _β), ϕ₃, Gridded(Linear())) + c = (p[6]-p[2]) / itp₃(α, abs(β)) + itp₄ = interpolate((_α, _β), ϕ₅, Gridded(Linear())) + ζ = p[4] + c * sign(β) * itp₄(α, abs(β)) + if abs(α - 1.0) < 0.1 + δ = ζ else - α = 0. - j = findfirst(>=(an), _ena) # _np.where(an <= _ena[:,0])[0] - (j === nothing || j == length(_ena)) && (j = length(_ena)) - t = (an-_ena[j-1])/(_ena[j]-_ena[j-1]) - α = (22-j-t)/10 - + δ = ζ - β * c * tan(π*α/2) end - if α < 0.5 - α = 0.5 - end - return AlphaStable(α=α, β=zero(α), scale=c, location=oftype(α, δ)) + return AlphaStable(α=α, β=β, scale=c, location=oftype(α, δ)) end +# """ +# 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 Distributions.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] +# (j === nothing || j == length(_ena)) && (j = length(_ena)) +# t = (an-_ena[j-1])/(_ena[j]-_ena[j-1]) +# α = (22-j-t)/10 + +# end +# if α < 0.5 +# α = 0.5 +# end +# return AlphaStable(α=α, β=zero(α), scale=c, location=oftype(α, δ)) +# end + """ Generate independent stable random numbers. @@ -278,7 +436,7 @@ function Distributions.fit(d::Type{<:AlphaSubGaussian}, x::AbstractVector{T}, m: α = d1.α; scale=d1.scale cov = zeros(T, m+1, m+1) xlen = length(x) - c = ((sum(x->abs(x)^p, x)/xlen)^(1/p))/scale + c = ((sum(abs.(x).^p)/xlen)^(1/p))/scale for i in 1:m tempxlen = xlen-mod(xlen, i) xtemp = reshape(x[1:end-mod(xlen, i)], i, tempxlen÷i) @@ -287,7 +445,7 @@ function Distributions.fit(d::Type{<:AlphaSubGaussian}, x::AbstractVector{T}, m: tempxlen = size(xtemp, 1)*size(xtemp, 2) end xtemp = reshape(xtemp', 2, tempxlen÷2) - @views r = (2/(c^p))*(scale^(2-p))*(xtemp[1, :]'*((sign.(xtemp[2, :]).*(abs.(xtemp[2, :]).^(p-1)))))/(tempxlen/2) + r = (2/(c^p))*(scale^(2-p))*(xtemp[1, :]'*((sign.(xtemp[2, :]).*(abs.(xtemp[2, :]).^(p-1)))))/(tempxlen/2) cov[diagind(cov, i)] .+= r end cov = (cov+cov')+2*(scale^2)*I(m+1) @@ -295,4 +453,4 @@ function Distributions.fit(d::Type{<:AlphaSubGaussian}, x::AbstractVector{T}, m: AlphaSubGaussian(α=α, R=cov, n=length(x)) end -end # module +end # module \ No newline at end of file From dc2dfadb4962216cb2a1de1f478986047602ea2a Mon Sep 17 00:00:00 2001 From: ymtoo Date: Wed, 14 Oct 2020 00:39:22 +0800 Subject: [PATCH 2/9] Add tests for special cases --- test/runtests.jl | 34 +++++++++++++++++++++------------- 1 file changed, 21 insertions(+), 13 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 0ecf0f2..10d5419 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,20 +1,34 @@ using AlphaStableDistributions -using Test, Random +using Test, Random, Distributions @testset "AlphaStableDistributions.jl" begin - for α in 0.5:0.1:2 + for α in 0.6:0.1:2 d1 = AlphaStable(α=α) s = rand(d1, 100000) d2 = fit(AlphaStable, s) @test d1.α ≈ d2.α rtol=0.1 - @test d1.β ≈ d2.β rtol=0.1 + @test d1.β ≈ d2.β atol=0.2 @test d1.scale ≈ d2.scale rtol=0.1 - @test d1.location ≈ d2.location atol=0.03 + @test d1.location ≈ d2.location atol=0.2 end + xnormal = rand(Normal(3.0, 4.0), 96000) + d = fit(AlphaStable, xnormal) + @test d.α ≈ 2 rtol=0.2 + @test d.β ≈ 0 atol=0.2 + @test d.scale ≈ 4/√2 rtol=0.2 + @test d.location ≈ 3 rtol=0.1 + + xcauchy = rand(Cauchy(3.0, 4.0), 96000) + d = fit(AlphaStable, xcauchy) + @test d.α ≈ 1 rtol=0.2 + @test d.β ≈ 0 atol=0.2 + @test d.scale ≈ 4 rtol=0.2 + @test d.location ≈ 3 rtol=0.1 + for α in 1.1:0.1:1.9 d = AlphaSubGaussian(n=96000, α=α) @@ -25,9 +39,9 @@ using Test, Random d3 = fit(AlphaStable, x) @test d3.α ≈ α rtol=0.2 - @test d3.β == 0 + @test d3.β ≈ 0 atol=0.2 @test d3.scale ≈ 1 rtol=0.2 - @test d3.location ≈ 0 atol=0.03 + @test d3.location ≈ 0 atol=0.1 end d4 = AlphaSubGaussian(n=96000) @@ -85,10 +99,4 @@ end # # using Cthulhu # d = AlphaSubGaussian(n=96) -# @descend_code_warntype rand(Random.GLOBAL_RNG, d) -# -# using AlphaStableDistributions -# d1 = AlphaStable(α=1.5) -# s = rand(d1, 100000) -# using ThreadsX -# @btime fit($AlphaStable, $s, $ThreadsX.MergeSort) +# @descend_code_warntype rand(Random.GLOBAL_RNG, d) \ No newline at end of file From 76f917876dc9dea74152cc87284ed3e0ce4812b1 Mon Sep 17 00:00:00 2001 From: ymtoo Date: Wed, 14 Oct 2020 10:59:11 +0800 Subject: [PATCH 3/9] remove changes --- src/AlphaStableDistributions.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/AlphaStableDistributions.jl b/src/AlphaStableDistributions.jl index 1948c67..5fa4ab3 100644 --- a/src/AlphaStableDistributions.jl +++ b/src/AlphaStableDistributions.jl @@ -436,7 +436,7 @@ function Distributions.fit(d::Type{<:AlphaSubGaussian}, x::AbstractVector{T}, m: α = d1.α; scale=d1.scale cov = zeros(T, m+1, m+1) xlen = length(x) - c = ((sum(abs.(x).^p)/xlen)^(1/p))/scale + c = ((sum(x->abs(x)^p, x)/xlen)^(1/p))/scale for i in 1:m tempxlen = xlen-mod(xlen, i) xtemp = reshape(x[1:end-mod(xlen, i)], i, tempxlen÷i) @@ -445,7 +445,7 @@ function Distributions.fit(d::Type{<:AlphaSubGaussian}, x::AbstractVector{T}, m: tempxlen = size(xtemp, 1)*size(xtemp, 2) end xtemp = reshape(xtemp', 2, tempxlen÷2) - r = (2/(c^p))*(scale^(2-p))*(xtemp[1, :]'*((sign.(xtemp[2, :]).*(abs.(xtemp[2, :]).^(p-1)))))/(tempxlen/2) + @views r = (2/(c^p))*(scale^(2-p))*(xtemp[1, :]'*((sign.(xtemp[2, :]).*(abs.(xtemp[2, :]).^(p-1)))))/(tempxlen/2) cov[diagind(cov, i)] .+= r end cov = (cov+cov')+2*(scale^2)*I(m+1) From 9a3effdb4b44c86c76b8a866240872ee79844f22 Mon Sep 17 00:00:00 2001 From: ymtoo Date: Wed, 14 Oct 2020 11:06:13 +0800 Subject: [PATCH 4/9] Reduce atol --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 10d5419..29f0571 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -12,7 +12,7 @@ using Test, Random, Distributions @test d1.α ≈ d2.α rtol=0.1 @test d1.β ≈ d2.β atol=0.2 @test d1.scale ≈ d2.scale rtol=0.1 - @test d1.location ≈ d2.location atol=0.2 + @test d1.location ≈ d2.location atol=0.1 end xnormal = rand(Normal(3.0, 4.0), 96000) From 189f75fa2fb3051ec186ecbb536505e2762c5816 Mon Sep 17 00:00:00 2001 From: ymtoo Date: Wed, 14 Oct 2020 11:11:25 +0800 Subject: [PATCH 5/9] Retain test code --- test/runtests.jl | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 29f0571..135dfed 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -99,4 +99,10 @@ end # # using Cthulhu # d = AlphaSubGaussian(n=96) -# @descend_code_warntype rand(Random.GLOBAL_RNG, d) \ No newline at end of file +# @descend_code_warntype rand(Random.GLOBAL_RNG, d) +# +# using AlphaStableDistributions +# d1 = AlphaStable(α=1.5) +# s = rand(d1, 100000) +# using ThreadsX +# @btime fit($AlphaStable, $s, $ThreadsX.MergeSort) \ No newline at end of file From fb74898c9a7613190ac3acceb14f42ff5682bf04 Mon Sep 17 00:00:00 2001 From: ymtoo Date: Wed, 14 Oct 2020 11:16:19 +0800 Subject: [PATCH 6/9] Update docstring --- src/AlphaStableDistributions.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/AlphaStableDistributions.jl b/src/AlphaStableDistributions.jl index 5fa4ab3..124af76 100644 --- a/src/AlphaStableDistributions.jl +++ b/src/AlphaStableDistributions.jl @@ -177,10 +177,11 @@ const ϕ₅ = [ ] """ + fit(d::Type{<:AlphaStable}, x; alg=QuickSort) + Fit an α stable distribution to data. -:param x: data -:returns: (α, β, c, δ) +returns `AlphaStable` α, β, c and δ are the characteristic exponent, skewness parameter, scale parameter (dispersion^1/α) and location parameter respectively. From 220eab8d3510ec9f79539f194c4fa8f287445308 Mon Sep 17 00:00:00 2001 From: ymtoo Date: Thu, 15 Oct 2020 10:30:41 +0800 Subject: [PATCH 7/9] Add SymmetricAlphaStable --- src/AlphaStableDistributions.jl | 117 +++++++++++++++++--------------- 1 file changed, 61 insertions(+), 56 deletions(-) diff --git a/src/AlphaStableDistributions.jl b/src/AlphaStableDistributions.jl index 124af76..99ebe1a 100644 --- a/src/AlphaStableDistributions.jl +++ b/src/AlphaStableDistributions.jl @@ -5,7 +5,7 @@ using StatsBase, Distributions, StaticArrays using MAT, SpecialFunctions, ToeplitzMatrices using Interpolations -export AlphaStable, AlphaSubGaussian, fit +export AlphaStable, SymmetricAlphaStable, AlphaSubGaussian, fit Base.@kwdef struct AlphaStable{T} <: Distributions.ContinuousUnivariateDistribution α::T = 1.5 @@ -33,27 +33,25 @@ Statistics.var(d::AlphaStable) = d.α == 2 ? 2d.scale^2 : Inf # 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 -# ] - # lookup tables 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 +] const _να = [ 2.439 2.500 @@ -183,8 +181,8 @@ Fit an α stable distribution to data. returns `AlphaStable` -α, β, c and δ are the characteristic exponent, skewness parameter, scale parameter -(dispersion^1/α) and location parameter respectively. +α∈[0.6,2.0], β∈[-1,1] , c∈[0,∞] and δ∈[-∞,∞] are the characteristic exponent, +skewness parameter, scale parameter (dispersion^1/α) and location parameter respectively. α, β, c and δ are computed based on McCulloch (1986) fractile. """ @@ -214,39 +212,46 @@ function Distributions.fit(::Type{<:AlphaStable}, x) return AlphaStable(α=α, β=β, scale=c, location=oftype(α, δ)) end -# """ -# 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 Distributions.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] -# (j === nothing || j == length(_ena)) && (j = length(_ena)) -# t = (an-_ena[j-1])/(_ena[j]-_ena[j-1]) -# α = (22-j-t)/10 - -# end -# if α < 0.5 -# α = 0.5 -# end -# return AlphaStable(α=α, β=zero(α), scale=c, location=oftype(α, δ)) -# end +Base.@kwdef struct SymmetricAlphaStable{T} <: Distributions.ContinuousUnivariateDistribution + α::T = 1.5 + scale::T = one(α) + location::T = zero(α) +end + +""" + fit(d::Type{<:SymmetricAlphaStable}, x; alg=QuickSort) + +Fit a symmetric α stable distribution to data. + +returns `SymmetricAlphaStable` + +α∈[1,2], c∈[0,∞] and δ∈[-∞,∞] are the characteristic exponent, scale parameter +(dispersion^1/α) and location parameter respectively. + +α is computed based on McCulloch (1986) fractile. +scale is computed based on Fama & Roll (1971) fractile. +location is the 50% trimmed mean of the sample. +""" +function Distributions.fit(::Type{<:SymmetricAlphaStable}, 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] + (j === nothing || j == length(_ena)) && (j = length(_ena)) + t = (an-_ena[j-1])/(_ena[j]-_ena[j-1]) + α = (22-j-t)/10 + + end + if α < 0.5 + α = 0.5 + end + return AlphaStable(α=α, β=zero(α), scale=c, location=oftype(α, δ)) +end """ Generate independent stable random numbers. From b39b9889d893157947e54aa9b530425175d87247 Mon Sep 17 00:00:00 2001 From: ymtoo Date: Thu, 15 Oct 2020 10:37:39 +0800 Subject: [PATCH 8/9] Add more tests --- test/runtests.jl | 47 +++++++++++++++++++++++++---------------------- 1 file changed, 25 insertions(+), 22 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 135dfed..7eaf0e7 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,34 +3,37 @@ using Test, Random, Distributions @testset "AlphaStableDistributions.jl" begin - for α in 0.6:0.1:2 - d1 = AlphaStable(α=α) - s = rand(d1, 100000) + stabletypes = [AlphaStable,SymmetricAlphaStable] + αs = [0.6:0.1:2,1:0.1:2] + for (i, stabletype) in enumerate(stabletypes) + for α in αs[i] + d1 = AlphaStable(α=α) + s = rand(d1, 100000) - d2 = fit(AlphaStable, s) + d2 = fit(stabletype, s) - @test d1.α ≈ d2.α rtol=0.1 - @test d1.β ≈ d2.β atol=0.2 - @test d1.scale ≈ d2.scale rtol=0.1 - @test d1.location ≈ d2.location atol=0.1 - end + @test d1.α ≈ d2.α rtol=0.1 + stabletype != SymmetricAlphaStable && @test d1.β ≈ d2.β atol=0.2 + @test d1.scale ≈ d2.scale rtol=0.1 + @test d1.location ≈ d2.location atol=0.1 + end - xnormal = rand(Normal(3.0, 4.0), 96000) - d = fit(AlphaStable, xnormal) - @test d.α ≈ 2 rtol=0.2 - @test d.β ≈ 0 atol=0.2 - @test d.scale ≈ 4/√2 rtol=0.2 - @test d.location ≈ 3 rtol=0.1 + xnormal = rand(Normal(3.0, 4.0), 96000) + d = fit(stabletype, xnormal) + @test d.α ≈ 2 rtol=0.2 + stabletype != SymmetricAlphaStable && @test d.β ≈ 0 atol=0.2 + @test d.scale ≈ 4/√2 rtol=0.2 + @test d.location ≈ 3 rtol=0.1 - xcauchy = rand(Cauchy(3.0, 4.0), 96000) - d = fit(AlphaStable, xcauchy) - @test d.α ≈ 1 rtol=0.2 - @test d.β ≈ 0 atol=0.2 - @test d.scale ≈ 4 rtol=0.2 - @test d.location ≈ 3 rtol=0.1 + xcauchy = rand(Cauchy(3.0, 4.0), 96000) + d = fit(stabletype, xcauchy) + @test d.α ≈ 1 rtol=0.2 + @test d.β ≈ 0 atol=0.2 + @test d.scale ≈ 4 rtol=0.2 + @test d.location ≈ 3 rtol=0.1 + end for α in 1.1:0.1:1.9 - d = AlphaSubGaussian(n=96000, α=α) x = rand(d) x2 = copy(x) From 6bbda71e36faf3d6355926d017206a598b46b50b Mon Sep 17 00:00:00 2001 From: ymtoo Date: Thu, 15 Oct 2020 10:52:13 +0800 Subject: [PATCH 9/9] Fix bugs --- src/AlphaStableDistributions.jl | 2 +- test/runtests.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/AlphaStableDistributions.jl b/src/AlphaStableDistributions.jl index 99ebe1a..b6acba3 100644 --- a/src/AlphaStableDistributions.jl +++ b/src/AlphaStableDistributions.jl @@ -250,7 +250,7 @@ function Distributions.fit(::Type{<:SymmetricAlphaStable}, x) if α < 0.5 α = 0.5 end - return AlphaStable(α=α, β=zero(α), scale=c, location=oftype(α, δ)) + return SymmetricAlphaStable(α=α, scale=c, location=oftype(α, δ)) end """ diff --git a/test/runtests.jl b/test/runtests.jl index 7eaf0e7..d02bed5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -28,7 +28,7 @@ using Test, Random, Distributions xcauchy = rand(Cauchy(3.0, 4.0), 96000) d = fit(stabletype, xcauchy) @test d.α ≈ 1 rtol=0.2 - @test d.β ≈ 0 atol=0.2 + stabletype != SymmetricAlphaStable && @test d.β ≈ 0 atol=0.2 @test d.scale ≈ 4 rtol=0.2 @test d.location ≈ 3 rtol=0.1 end