Add Box-Pierce and Ljung-Box Autocorrelation Tests (#87)

* Add Box-Pierce and Ljung-Box tests * Add more informative output print * Change indent * Fix lag<dof error * Reformat data matrices for tests
JuliaStats · Mar 20, 2017 · d3597b2 · d3597b2
1 parent c44172b
commit d3597b2
Show file tree

Hide file tree

Showing 3 changed files with 215 additions and 0 deletions.
diff --git a/src/HypothesisTests.jl b/src/HypothesisTests.jl
@@ -142,4 +142,5 @@ include("z.jl")
 include("wilcoxon.jl")
 include("power_divergence.jl")
 include("anderson_darling.jl")
+include("box_test.jl")
 end
diff --git a/src/box_test.jl b/src/box_test.jl
@@ -0,0 +1,121 @@
+# box_test.jl
+# Box-Pierce and Ljung-Box tests for autocorrelation
+#
+# 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 BoxPierceTest, LjungBoxTest
+
+# Box-Pierce test
+
+immutable BoxPierceTest <: HypothesisTest
+    y::Vector{Float64}  # time series vector
+    n::Int              # number of observations
+    lag::Int            # number of lags in test statistic
+    dof::Int            # degrees-of-freedom correction (if y is residual)
+    Q::Float64          # test statistic
+end
+
+"""
+    BoxPierceTest(y, lag, dof=0)
+
+Compute the Box-Pierce `Q` statistic to test the null hypothesis of
+independence in a time series `y`.
+
+`lag` specifies the number of lags used in the construction of `Q`. When
+testing the residuals of an estimated model, `dof` has to be set to the number
+of estimated parameters. E.g., when testing the residuals of an ARIMA(p,0,q)
+model, set `dof=p+q`.
+
+External links
+
+* [Box-Pierce test on Wikipedia](https://en.wikipedia.org/wiki/Ljung–Box_test#Box-Pierce_test)
+"""
+function BoxPierceTest{T<:Real}(y::AbstractVector{T}, lag::Int, dof::Int=0)
+    if dof>=lag
+        throw(ArgumentError("Number of lags has to be larger than degrees of" *
+        " freedom correction"))
+    end
+    n = size(y,1)
+    Q = n * first(sum(k -> autocor(y,k:k).^2, 1:lag))
+    BoxPierceTest(y,n,lag,dof,Q)
+end
+
+testname(::BoxPierceTest) = "Box-Pierce autocorrelation test"
+population_param_of_interest(x::BoxPierceTest) = ("autocorrelations up to lag k",
+    "all zero", NaN)
+
+function show_params(io::IO, x::BoxPierceTest, ident)
+    println(io, ident, "number of observations:         ", x.n)
+    println(io, ident, "number of lags:                 ", x.lag)
+    println(io, ident, "degrees of freedom correction:  ", x.dof)
+    println(io, ident, "Q statistic:                    ", x.Q)
+end
+
+pvalue(x::BoxPierceTest) = pvalue(Chisq(x.lag-x.dof), x.Q; tail=:right)
+
+#Ljung-Box test
+
+immutable LjungBoxTest <: HypothesisTest
+    y::Vector{Float64}  # time series vector
+    n::Int              # number of observations
+    lag::Int            # number of lags in test statistic
+    dof::Int            # degrees-of-freedom correction (if y is residual)
+    Q::Float64          # test statistic
+end
+
+"""
+    LjungBoxTest(y, lag, dof=0)
+
+Compute the Ljung-Box `Q` statistic to test the null hypothesis of
+independence in a time series `y`.
+
+`lag` specifies the number of lags used in the construction of `Q`. When
+testing the residuals of an estimated model, `dof` has to be set to the number
+of estimated parameters. E.g., when testing the residuals of an ARIMA(p,0,q)
+model, set `dof=p+q`.
+
+External links
+
+* [Ljung-Box test on Wikipedia](https://en.wikipedia.org/wiki/Ljung–Box_test)
+"""
+function LjungBoxTest{T<:Real}(y::AbstractVector{T}, lag::Int, dof::Int=0)
+    if dof>=lag
+        throw(ArgumentError("Number of lags has to be larger than degrees of" *
+        " freedom correction"))
+    end
+    n = size(y,1)
+    Q = n*(n+2)* first(sum(k -> autocor(y,k:k).^2/(n-k), 1:lag))
+    LjungBoxTest(y,n,lag,dof,Q)
+end
+
+testname(::LjungBoxTest) = "Ljung-Box autocorrelation test"
+population_param_of_interest(x::LjungBoxTest) = ("autocorrelations up to lag k",
+    "all zero", NaN)
+
+function show_params(io::IO, x::LjungBoxTest, ident)
+    println(io, ident, "number of observations:         ", x.n)
+    println(io, ident, "number of lags:                 ", x.lag)
+    println(io, ident, "degrees of freedom correction:  ", x.dof)
+    println(io, ident, "Q statistic:                    ", x.Q)
+end
+
+pvalue(x::LjungBoxTest) = pvalue(Chisq(x.lag-x.dof), x.Q; tail=:right)
diff --git a/test/box_test.jl b/test/box_test.jl
@@ -0,0 +1,93 @@
+using HypothesisTests, Base.Test
+
+sim_data_h0=[
+    0.297287984535462;0.382395967790608;-0.597634476728231;-0.0104452446373756;
+    -0.839026854388764;0.311111338498334;2.29508782383731;-2.26708634880053;
+    0.529965576166746;0.431421526422912;0.583708287568779;0.963271605038191;
+    0.458790955053717;-0.522336757421508;0.408395838324752;-0.0504512299336653;
+    -0.693653643803886;-1.77383325340806;0.120765964937439;-0.757633037720528;
+    -1.72975618139069;0.795948622004616;0.670061981256062;0.550852377500761;
+    -0.0633745924295606;1.33694278225092;-0.0731486333377393;-0.745464316640756;
+    -1.22005512853465;-0.0531773375276925;-0.165136357802834;-2.11537024842731;
+    -0.0667679865065892;1.22246184846806;0.567695459738444;0.763498928323387;
+    0.378598879853350;-0.645597056636597;-0.664646244391031;-1.80302632421089;
+    -1.28904760384974;-0.335270079518118;0.0704676013611580;0.341793767379823;
+    1.73516566605438;1.29991596431019;0.206364081402617;-1.00885995553760;
+    -0.850055670315398;1.12940943588433;-1.34853345649183;1.24029090042906;
+    -0.0556470932073243;0.759390213945460;-0.0302668866467625;-0.361082667756094;
+    -2.00730320253256;-0.568859127346638;-1.14962752439700;1.88062217463119;
+    -1.47062454138738;-0.535210970165732;-0.963543959825501;-1.38510806767591;
+    0.134300237561270;-0.616116666008819;-1.71998935631974;0.320769279987411;
+    -1.44737289788151;-0.492271188069566;-0.413087383983183;1.00192851245595;
+    1.17468132597511;-1.32697034506205;-1.59321486804407;-0.438709975758797;
+    1.22579749799787;1.48909027806961;2.15943318710085;-0.909479018037503;
+    1.14782030061250;-0.505072196522373;-0.267067035218130;1.49908239393073;
+    0.797304160175390;-0.171640443307097;-0.469079609930096;0.217624021208600;
+    0.359146407602350;0.320024715708058;0.259216328055686;0.459695985041628;
+    0.209986360042143;1.52317539117726;1.25979460361220;-0.165283101114261;
+    0.156194494887636;-0.814315294376699
+    ]
+
+t = HypothesisTests.BoxPierceTest(sim_data_h0,2,1)
+
+@test t.n == 98
+@test t.lag == 2
+@test t.dof == 1
+@test t.Q ≈ 1.233942980734545
+@test pvalue(t) ≈ 0.2666415904008932
+show(IOBuffer(), t)
+
+t = HypothesisTests.LjungBoxTest(sim_data_h0,5,2)
+
+@test t.n == 98
+@test t.lag == 5
+@test t.dof == 2
+@test t.Q ≈ 3.2090126519163626
+@test pvalue(t) ≈ 0.36050846449240337
+show(IOBuffer(), t)
+
+sim_data_h1 = [
+    0.297287984535462;0.739141549233162;0.200148986990924;0.00799107498178517;
+    -0.889482260507899;-0.758664696605681;1.65153486606286;-0.0576451005433905;
+    -0.0346690043041814;0.407112251420912;1.08264369056513;2.14031035829007;
+    2.70237027783226;2.07841446849019;2.09178211716330;1.83616297011523;
+    0.882207285185406;-1.26603340222014;-1.66313630328235;-2.37358658099330;
+    -4.07911918759795;-3.38691842881493;-2.17050437704247;-1.03767734630572;
+    -0.657436094883687;0.859322672282215;1.15526940186602;0.383062163913809;
+    -1.10696135239789;-1.49644960957930;-1.62878748357863;-3.62098034584788;
+    -3.92330815645045;-2.39921383551812;-1.13436869594816;0.122020643841029;
+    0.865334261247033;0.356197863707534;-0.496809086316101;-2.50605658690248;
+    -4.14727278223788;-4.56018044213284;-4.15756709452688;-3.27923261341258;
+    -0.952643341682653;1.14051373831478;1.86077356988515;0.881914206830145;
+    -0.349990693084770;0.444846342133559;-0.709720638006125;0.255172232181642;
+    0.463475776812484;1.23900947646595;1.31750175206863;0.848216591786478;
+    -1.38469381800938;-2.48495668649384;-3.71616740278679;-1.83329170276481;
+    -2.55572436386911;-3.05209269597923;-3.85933788583984;-5.10068572188995;
+    -4.82872126295472;-4.88037646498750;-6.12782473541832;-5.56850746301833;
+    -6.29123443287801;-6.37120026861768;-6.17115737646099;-4.49210025871193;
+    -2.36449177154091;-2.81673039329756;-4.26394380853887;-4.71042342801617;
+    -3.14752747305988;-0.874815661197385;2.05391263558195;1.81766084302006;
+    2.71283952156198;2.20503697644599;1.56512548004846;2.71572187705509;
+    3.58663276862696;3.31760231592873;2.43605333859629;2.14560733274553;
+    2.20305920531810;2.32001356226611;2.38231484117959;2.62246972577731;
+    2.64225557862103;3.90714116778931;5.15568733137306;4.84939934619662;
+    4.42876751091166;3.04538591485831
+    ]
+
+t = HypothesisTests.BoxPierceTest(sim_data_h1,3)
+
+@test t.n == 98
+@test t.lag == 3
+@test t.dof == 0
+@test t.Q ≈ 176.16899390632153
+@test pvalue(t) ≈ 5.925682677971866e-38
+show(IOBuffer(), t)
+
+t = HypothesisTests.LjungBoxTest(sim_data_h1,12)
+
+@test t.n == 98
+@test t.lag == 12
+@test t.dof == 0
+@test t.Q ≈ 271.8894095341075
+@test pvalue(t) ≈ 3.6622231247462687e-51
+show(IOBuffer(), t)