Add Hansen's Standardized Skewed T Distribution (#86)

* standardized t dist: consistenly use greek letter

* add Hansen's Standardized Skewed T Distribution

* incorporate code review changes for StdSkewT

* change distname for StdSkewT

* fixup! incorporate code review changes for StdSkewT

* add tests for Hansen's Standardized Skewed T Distribution

* add documentation for Hansen's Standardized Skewed T Distribution

* optimize logkernel function for StdSkewT

hand optimized the logkernel function for StdSkewT instead of using the support
{a,b,c} functions which were wasteful since calls to them were repeated. A
speed improvement of 4x has been attained on `fit(GARCH{1,1},BG96)`

* fixup! add tests for Hansen's Standardized Skewed T Distribution

docs/src/univariatetypehierarchy.md

@@ -148,6 +148,11 @@ StdT{Float64}(coefs=[3.0])
 julia> StdGED(1) # convenience constructor
 StdGED{Float64}(coefs=[1.0])
 ```
+* [`StdSkewT`](@ref), [Hansen's standardized skewed Student's ``t`` ditribution](https://en.wikipedia.org/wiki/Skewed_generalized_t_distribution#cite_note-hansen-8):
+```jldoctest TYPES
+julia> StdSkewT(3,-0.3) # convenience constructor
+StdSkewT{Float64}(coefs=[3.0,-0.3])
+```
 ### User-defined standardized distributions
 Apart from the natively supported standardized distributions, it is possible to wrap a continuous univariate distribution from the [Distributions package](https://github.com/JuliaStats/Distributions.jl) in the [`Standardized`](@ref) wrapper type. Below, we reimplement the standardized normal distribution:
 

src/ARCHModels.jl

@@ -49,7 +49,7 @@ import StatsBase: StatisticalModel, stderror, loglikelihood, nobs, fit, fit!, co
 				  informationmatrix, islinear, score, vcov, residuals, predict
 import StatsModels: TableRegressionModel
 export ARCHModel, UnivariateARCHModel, UnivariateVolatilitySpec, StandardizedDistribution, Standardized, MeanSpec,
-       simulate, simulate!, selectmodel, StdNormal, StdT, StdGED, Intercept, Regression,
+       simulate, simulate!, selectmodel, StdNormal, StdT, StdGED, StdSkewT, Intercept, Regression,
        NoIntercept, ARMA, AR, MA, BG96, volatilities, mean, quantile, VaRs, pvalue, means, VolatilitySpec,
 	   MultivariateVolatilitySpec, MultivariateStandardizedDistribution, MultivariateARCHModel, MultivariateStdNormal,
 	   EGARCH, ARCH, GARCH, TGARCH, ARCHLMTest, DQTest,

src/univariatestandardizeddistributions.jl

@@ -142,14 +142,14 @@ struct StdT{T} <: StandardizedDistribution{T}
 end
 
 """
-    StdT(v)
+    StdT(ν)
 
-Create a standardized t distribution with `v` degrees of freedom. `ν`` can be passed
+Create a standardized t distribution with `ν` degrees of freedom. `ν`` can be passed
 as a scalar or vector.
 """
 StdT(ν) = StdT([ν])
 StdT(ν::Integer) = StdT(float(ν))
-StdT(v::Vector{T}) where {T} = StdT{T}(v)
+StdT(ν::Vector{T}) where {T} = StdT{T}(ν)
 (rand(d::StdT{T})::T) where {T}  =  (ν=d.coefs[1]; tdistrand(ν)*sqrt((ν-2)/ν))
 @inline kernelinvariants(::Type{<:StdT}, coefs) = (1/ (coefs[1]-2),)
 @inline logkernel(::Type{<:StdT}, x, coefs, iv) = (-(coefs[1] + 1) / 2) * log1p(abs2(x) *iv)
@@ -184,8 +184,8 @@ function startingvals(::Type{<:StdT}, data::Array{T}) where {T}
 end
 
 function quantile(dist::StdT, q::Real)
-    v = dist.coefs[1]
-    tdistinvcdf(v, q)*sqrt((v-2)/v)
+    ν = dist.coefs[1]
+    tdistinvcdf(ν, q)*sqrt((ν-2)/ν)
 end
 
 ################################################################################
@@ -253,3 +253,68 @@ function quantile(dist::StdGED, q::Real)
     qq = 2*q-1
     return sign(qq) * (gammainvcdf(ip, 1., abs(qq)))^ip/kernelinvariants(StdGED, [p])[1]
 end
+
+################################################################################
+#Hansen's SKT-Distribution
+
+"""
+    StdSkewT{T} <: StandardizedDistribution{T}
+
+Hansen's standardized (mean zero, variance one) Skewed Student's t distribution.
+"""
+struct StdSkewT{T} <: StandardizedDistribution{T}
+    coefs::Vector{T}
+    function StdSkewT{T}(coefs::Vector) where {T}
+        length(coefs) == 2 || throw(NumParamError(2, length(coefs)))
+        new{T}(coefs)
+    end
+end
+
+"""
+    StdSkewT(v,λ)
+
+Create a standardized skewed t distribution with `v` degrees of freedom and `λ` shape parameter. `ν,λ`` can be passed
+as scalars or vectors.
+"""
+StdSkewT(ν,λ) = StdSkewT([float(ν), float(λ)])
+StdSkewT(coefs::Vector{T}) where {T} = StdSkewT{T}(coefs)
+
+(rand(d::StdSkewT{T}) where {T} = (quantile(d, rand())))
+
+@inline a(d::Type{<:StdSkewT}, coefs) =  (ν=coefs[1];λ=coefs[2]; 4λ*c(d,coefs) * ((ν-2)/(ν-1)))
+@inline b(d::Type{<:StdSkewT}, coefs) =  (ν=coefs[1];λ=coefs[2]; sqrt(1+3λ^2-a(d,coefs)^2))
+@inline c(d::Type{<:StdSkewT}, coefs) =  (ν=coefs[1];λ=coefs[2]; gamma((ν+1)/2) / (sqrt(π*(ν-2)) * gamma(ν/2)))
+
+@inline kernelinvariants(::Type{<:StdSkewT}, coefs) = (1/ (coefs[1]-2),)
+@inline function logkernel(d::Type{<:StdSkewT}, x, coefs, iv)
+    ν=coefs[1]
+    λ=coefs[2]
+    c = gamma((ν+1)/2) / (sqrt(π*(ν-2)) * gamma(ν/2))
+    a = 4λ * c * ((ν-2)/(ν-1))
+    b = sqrt(1 + 3λ^2 -a^2)
+    λsign = x < (-a/b) ? -1 : 1
+    (-(ν + 1) / 2) * log1p(1/abs2(1+λ*λsign) * abs2(b*x + a) *iv)
+end
+@inline logconst(d::Type{<:StdSkewT}, coefs)  = (log(b(d,coefs))+(log(c(d,coefs))))
+
+nparams(::Type{<:StdSkewT}) = 2
+coefnames(::Type{<:StdSkewT}) = ["ν", "λ"]
+distname(::Type{<:StdSkewT}) = "Hansen's Skewed t"
+
+function constraints(::Type{<:StdSkewT}, ::Type{T}) where {T}
+    lower = T[20/10, -one(T)]
+    upper = T[Inf,one(T)]
+    return lower, upper
+end
+
+startingvals(::Type{<:StdSkewT}, data::Array{T}) where {T<:AbstractFloat} = [startingvals(StdT, data)..., zero(T)]
+
+function quantile(d::StdSkewT{T}, q::T) where T
+    ν = d.coefs[1]
+    λ = d.coefs[2]
+    a_val = a(typeof(d),d.coefs)
+    b_val = b(typeof(d),d.coefs)
+    λconst = q < (1 - λ)/2 ? (1 - λ) : (1 + λ)
+    quant_numer = q < (1 - λ)/2 ? q : (q + λ)
+    1/b_val * ((λconst) * sqrt((ν-2)/ν) * tdistinvcdf(ν, quant_numer/λconst) - a_val)
+end

test/runtests.jl

@@ -302,6 +302,7 @@ end
     @test_throws ARCHModels.NumParamError NoIntercept([1.])
     @test_throws ARCHModels.NumParamError StdNormal([1.])
     @test_throws ARCHModels.NumParamError StdT([1., 2.])
+    @test_throws ARCHModels.NumParamError StdSkewT([2.])
     @test_throws ARCHModels.NumParamError StdGED([1., 2.])
     @test_throws ARCHModels.NumParamError Regression([1], [1 2; 3 4])
     at = zeros(10)
@@ -365,6 +366,40 @@ end
                                        3.834033118707376,
                                        2.9918481871096083], rtol=1e-4))
     end
+    @testset "HansenSkewedT" begin
+       Random.seed!(1)
+       data = rand(StdSkewT(4,-0.3), T)
+       spec = GARCH{1, 1}([1., .9, .05])
+       @test fit(StdSkewT, data).coefs[1] ≈ 4.115817082082046 rtol=1e-4
+       @test fit(StdSkewT, data).coefs[2] ≈ -0.29195923568201576 rtol=1e-4
+       @test typeof(StdSkewT(3,0)) == typeof(StdSkewT(3.,0)) == typeof(StdSkewT([3,0.0]))
+       @test coefnames(StdSkewT) == ["ν", "λ"]
+       @test ARCHModels.nparams(StdSkewT) == 2
+       @test ARCHModels.distname(StdSkewT) == "Hansen's Skewed t"
+       @test ARCHModels.constraints(StdNormal{Float64}, Float64) == (Float64[], Float64[])
+       @test quantile(StdSkewT(3,0), 0.5) == 0
+       @test quantile(StdSkewT(3,0), .05) ≈ -1.3587150125838563
+       @test ARCHModels.constraints(StdSkewT{Float64}, Float64) == (Float64[20/10, -one(Float64)], Float64[Inf,one(Float64)])
+       Random.seed!(1);
+       dataskt = simulate(spec, T; dist=StdSkewT(4,-0.3)).data
+       Random.seed!(1);
+       datam = simulate(spec, T; dist=StdSkewT(4,-0.3), meanspec=Intercept(3)).data
+       am4 = selectmodel(GARCH, dataskt; dist=StdSkewT, meanspec=NoIntercept{Float64}(), show_trace=true)
+       am5 = selectmodel(GARCH, datam; dist=StdSkewT, show_trace=true)
+       @test coefnames(am5) == ["ω", "β₁", "α₁", "ν", "λ", "μ"]
+       @test all(coeftable(am4).cols[2] .== stderror(am4))
+       @test all(isapprox(coef(am4), [1.1129628842351935,
+                                      0.8875795519570414,
+                                      0.06520263915417537,
+                                      3.886881385708421,
+                                      -0.29751925081928116], rtol=1e-4))
+       @test all(isapprox(coef(am5), [1.1108927904860744,
+                                      0.8875825749083535,
+                                      0.0651093574134582,
+                                      3.8894659913968255,
+                                      -0.29684836199362896,
+                                      3.0047117757300317], rtol=1e-4))
+    end
     @testset "GED" begin
         Random.seed!(1)
         @test typeof(StdGED(3)) == typeof(StdGED(3.)) == typeof(StdGED([3.]))