Add Jarque-Bera normality test (#95)

src/HypothesisTests.jl

@@ -143,4 +143,5 @@ include("anderson_darling.jl")
 include("box_test.jl")
 include("breusch_godfrey.jl")
 include("augmented_dickey_fuller.jl")
+include("jarque_bera.jl")
 end

src/jarque_bera.jl

@@ -0,0 +1,87 @@
+# jarque_bera.jl
+# Jarque-Bera goodness-of-fit test
+#
+# Copyright (C) 2017   Benjamin Born
+#
+# Permission is hereby granted, free of charge, to any person obtaining
+# a copy of this software and associated documentation files (the
+# "Software"), to deal in the Software without restriction, including
+# without limitation the rights to use, copy, modify, merge, publish,
+# distribute, sublicense, and/or sell copies of the Software, and to
+# permit persons to whom the Software is furnished to do so, subject to
+# the following conditions:
+#
+# The above copyright notice and this permission notice shall be
+# included in all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
+# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+
+export JarqueBeraTest
+
+immutable JarqueBeraTest <: HypothesisTest
+    n::Int         # number of observations
+    JB::Float64    # test statistic
+    skew::Float64  # skewness
+    kurt::Float64  # excess kurtosis
+end
+
+"""
+    JarqueBeraTest(y::AbstractVector)
+
+Compute the Jarque-Bera statistic to test the null hypothesis that a real-valued vector `y`
+is normally distributed.
+
+Note that the approximation by the Chi-squared distribution does not work well and the
+speed of convergence is slow. In small samples, the test tends to be over-sized for
+nominal levels up to about 3% and under-sized for larger nominal levels (Mantalos, 2010).
+
+# References
+
+  * Panagiotis Mantalos, 2011, "The three different measures of the sample skewness and
+    kurtosis and the effects to the Jarque-Bera test for normality", International Journal
+    of Computational Economics and Econometrics, Vol. 2, No. 1,
+    [link](http://dx.doi.org/10.1504/IJCEE.2011.040576).
+
+# External links
+
+  * [Jarque-Bera test on Wikipedia](https://en.wikipedia.org/wiki/Jarque–Bera_test)
+"""
+function JarqueBeraTest{T<:Real}(y::AbstractVector{T})
+    n = length(y)
+    M = Base.promote_op(/, T, typeof(n))
+    m1r = m2r = m3r = m4r = zero(M)
+    @inbounds for yi in y # compute raw moments
+        m1r += yi / n
+        m2r += yi^2 / n
+        m3r += yi^3 / n
+        m4r += yi^4 / n
+    end
+    # compute central moments (http://mathworld.wolfram.com/CentralMoment.html)
+    m2 = -m1r^2 + m2r
+    m3 = 2 * m1r^3 - 3 * m1r * m2r + m3r
+    m4 = -3 * m1r^4 + 6 * m1r^2 * m2r - 4 * m1r * m3r + m4r
+
+    skew = m3 / m2^(3/2)
+    kurt = m4 / m2^2
+
+    stat = n * skew^2 / 6 + n * (kurt - 3)^2 / 24
+    JarqueBeraTest(n, stat, skew, kurt)
+end
+
+testname(::JarqueBeraTest) = "Jarque-Bera normality test"
+population_param_of_interest(x::JarqueBeraTest) =
+    ("skewness and kurtosis", "0 and 3", "$(x.skew) and $(x.kurt)")
+default_tail(test::JarqueBeraTest) = :right
+
+function show_params(io::IO, x::JarqueBeraTest, ident)
+    println(io, ident, "number of observations:         ", x.n)
+    println(io, ident, "JB statistic:                   ", x.JB)
+end
+
+pvalue(x::JarqueBeraTest) = pvalue(Chisq(2), x.JB; tail=:right)

test/jarque_bera.jl

@@ -0,0 +1,58 @@
+using HypothesisTests, Base.Test
+using HypothesisTests: default_tail
+
+sim_data_h0 = 1.0 .+ [
+    0.2972879845354616, 0.3823959677906078, -0.5976344767282311, -0.01044524463737564,
+    -0.839026854388764, 0.31111133849833383, 2.2950878238373105, -2.2670863488005306,
+    0.5299655761667461, 0.43142152642291204, 0.5837082875687786, 0.9632716050381906,
+    0.45879095505371686, -0.5223367574215084, 0.40839583832475224, -0.050451229933665284,
+    -0.6936536438038856, -1.7738332534080556, 0.12076596493743928, -0.7576330377205277,
+    -1.7297561813906863, 0.7959486220046159, 0.6700619812560624, 0.5508523775007609,
+    -0.06337459242956062, 1.3369427822509223, -0.07314863333773934, -0.7454643166407556,
+    -1.2200551285346526, -0.05317733752769253, -0.1651363578028337, -2.115370248427311,
+    -0.0667679865065892, 1.2224618484680618, 0.5676954597384435, 0.7634989283233873,
+    0.37859887985335006, -0.645597056636597, -0.6646462443910314, -1.8030263242108941,
+    -1.2890476038497443, -0.3352700795181175, 0.07046760136115796, 0.34179376737982303,
+    1.7351656660543822, 1.299915964310192, 0.20636408140261675, -1.0088599555376034,
+    -0.8500556703153979, 1.1294094358843267, -1.348533456491827, 1.2402909004290597,
+    -0.055647093207324264, 0.7593902139454598, -0.030266886646762546, -0.36108266775609366,
+    -2.0073032025325626, -0.5688591273466379, -1.1496275243969976, 1.880622174631185,
+    -1.4706245413873802, -0.5352109701657319, -0.9635439598255011, -1.3851080676759138,
+    0.13430023756127, -0.6161166660088192, -1.719989356319742, 0.32076927998741056,
+    -1.4473728978815112, -0.4922711880695664, -0.4130873839831828, 1.001928512455953,
+    1.174681325975111, -1.326970345062047, -1.593214868044067, -0.43870997575879667,
+    1.225797497997867, 1.4890902780696147, 2.1594331871008534, -0.9094790180375028,
+    1.1478203006125018, -0.5050721965223731, -0.2670670352181297, 1.499082393930729,
+    0.7973041601753904, -0.17164044330709693, -0.46907960993009634, 0.21762402120860017,
+    0.35914640760234967, 0.3200247157080577, 0.2592163280556861, 0.459695985041628,
+    0.20998636004214272, 1.5231753911772645, 1.2597946036121963, -0.16528310111426064,
+    0.15619449488763584, -0.8143152943766987, -1.7098682060045913, 1.4999378761140325,
+    -0.1289156385300165, -0.4501839715598432
+]
+
+t = JarqueBeraTest(sim_data_h0)
+
+@test t.n == 102
+@test t.JB ≈ 1.0204197333021678
+@test t.skew ≈ -0.020527653857777352
+@test t.kurt ≈ 2.5117242352057993
+@test pvalue(t) ≈ 0.6003695680393418
+@test default_tail(t) == :right
+show(IOBuffer(), t)
+
+sim_data_h1 = [
+    0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1,
+    0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0,
+    1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1,
+    1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1
+]
+
+t = JarqueBeraTest(sim_data_h1)
+
+@test t.n == 102
+@test t.JB ≈ 17.00081983024691
+@test t.skew ≈ 0.11785113019775637
+@test t.kurt ≈ 1.0138888888888888
+@test pvalue(t) ≈ 0.00020338498134114293
+@test default_tail(t) == :right
+show(IOBuffer(), t)

test/runtests.jl

@@ -10,15 +10,19 @@ macro test_ci_approx(x::Expr, y::Expr)
     end
 end
 
+include("anderson_darling.jl")
+include("augmented_dickey_fuller.jl")
 include("binomial.jl")
+include("box_test.jl")
+include("breusch_godfrey.jl")
 include("circular.jl")
+include("common.jl")
 include("fisher.jl")
+include("jarque_bera.jl")
 include("kolmogorov_smirnov.jl")
 include("kruskal_wallis.jl")
 include("mann_whitney.jl")
+include("power_divergence.jl")
 include("t.jl")
-include("z.jl")
 include("wilcoxon.jl")
-include("power_divergence.jl")
-include("anderson_darling.jl")
-include("common.jl")
+include("z.jl")