Repository: JuliaStats/HypothesisTests.jl
Pull request #180 (2020-02-12 19:14:04)

Add Diebold-Mariano test
Closes the issue #179 .

All tests were based on the Diebold-Mariano implementation on the forecast package of R. The script is the following.

```R
library("forecast")
# Test on in-sample one-step forecasts
f1 <- ets(WWWusage)
f2 <- auto.arima(WWWusage)
accuracy(f1)
accuracy(f2)
dm.test(residuals(f1),residuals(f2),h=1)
# Test on out-of-sample one-step forecasts
f1 <- ets(WWWusage[1:80])
f2 <- auto.arima(WWWusage[1:80])
f1.out <- ets(WWWusage[81:100],model=f1)
f2.out <- Arima(WWWusage[81:100],model=f2)
accuracy(f1.out)
accuracy(f2.out)
dm.test(residuals(f1.out),residuals(f2.out),h=1)
#p-value = 0.8871
dm.test(residuals(f1.out),residuals(f2.out),h=1,alternative="less")
#p-value = 0.4435
dm.test(residuals(f1.out),residuals(f2.out),h=1,alternative="greater")
#p-value = 0.5565

dm.test(residuals(f1.out),residuals(f2.out),h=1,power=1)
#p-value = 0.7818
dm.test(residuals(f1.out),residuals(f2.out),h=1,alternative="less",power=1)
#p-value = 0.6091
dm.test(residuals(f1.out),residuals(f2.out),h=1,alternative="greater",power=1)
#p-value = 0.3909

dm.test(residuals(f1.out),residuals(f2.out),h=10)
#p-value = 0.9362
dm.test(residuals(f1.out),residuals(f2.out),h=10,alternative="less")
#p-value = 0.4681
dm.test(residuals(f1.out),residuals(f2.out),h=10,alternative="greater")
#p-value = 0.5319

f1 <- ets(WWWusage)
f2 <- auto.arima(WWWusage)
accuracy(f1)
accuracy(f2)
dm.test(residuals(f1),residuals(f2),h=1)
# Test on out-of-sample one-step forecasts
f1 <- ets(WWWusage[1:80])
f2 <- auto.arima(WWWusage[1:80])
f1.out <- ets(WWWusage[81:100],model=f1)
f2.out <- Arima(WWWusage[81:100],model=f2)
accuracy(f1.out)
accuracy(f2.out)
dm.test(residuals(f1.out),residuals(f2.out),h=1)
#p-value = 0.8871
dm.test(residuals(f1.out),residuals(f2.out),h=1,alternative="less")
#p-value = 0.4435
dm.test(residuals(f1.out),residuals(f2.out),h=1,alternative="greater")
#p-value = 0.5565

dm.test(residuals(f1.out),residuals(f2.out),h=1,power=1)
#p-value = 0.7818
dm.test(residuals(f1.out),residuals(f2.out),h=1,alternative="less",power=1)
#p-value = 0.6091
dm.test(residuals(f1.out),residuals(f2.out),h=1,alternative="greater",power=1)
#p-value = 0.3909

dm.test(residuals(f1.out),residuals(f2.out),h=10)
#p-value = 0.9362
dm.test(residuals(f1.out),residuals(f2.out),h=10,alternative="less")
#p-value = 0.4681
dm.test(residuals(f1.out),residuals(f2.out),h=10,alternative="greater")
#p-value = 0.5319
```



johnczito: Once you've computed the test statistic, this is "just" a *t*-test. So you should make `DieboldMarianoTest <: TTest`, rename the `dof` field to `df` and the `stat` field to `t`. Then you can eliminate your `pvalue` method because it's already provided for all *t*-tests here:

https://github.com/JuliaStats/HypothesisTests.jl/blob/7f21943e13c8f24bca7c171b7263f36926b41702/src/t.jl#L31

Patch:
diff --git a/Project.toml b/Project.toml
index bcd4bbb..977e3b7 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "HypothesisTests"
 uuid = "09f84164-cd44-5f33-b23f-e6b0d136a0d5"
-version = "0.9.0"
+version = "0.10.0"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
@@ -14,10 +14,10 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 
 [compat]
 Combinatorics = "1"
-Distributions = "0.19, 0.20, 0.21, 0.22"
+Distributions = "0.19, 0.20, 0.21, 0.22, 0.23"
 Rmath = "0.5, 0.6"
-Roots = "0.7, 0.8"
-StatsBase = "0.27, 0.28, 0.29, 0.30, 0.31, 0.32"
+Roots = "0.7, 0.8, 1.0"
+StatsBase = "0.27, 0.28, 0.29, 0.30, 0.31, 0.32, 0.33"
 julia = "1"
 
 [extras]
diff --git a/docs/src/time_series.md b/docs/src/time_series.md
index bfe77d7..e3f8e96 100644
--- a/docs/src/time_series.md
+++ b/docs/src/time_series.md
@@ -22,4 +22,9 @@ JarqueBeraTest
 ## Augmented Dickey-Fuller test
 ```@docs
 ADFTest
+```
+
+## Diebold-Mariano test
+```@docs
+DieboldMarianoTest
 ```
\ No newline at end of file
diff --git a/src/augmented_dickey_fuller.jl b/src/augmented_dickey_fuller.jl
index 61f0474..67e45e3 100644
--- a/src/augmented_dickey_fuller.jl
+++ b/src/augmented_dickey_fuller.jl
@@ -34,12 +34,12 @@ struct ADFTest <: HypothesisTest
 end
 
 """
-    ADFTest(y, deterministic, lag)
+    ADFTest(y::AbstractVector{T}, deterministic::Symbol, lag::Int) where T<:Real
 
 Compute the augmented Dickey-Fuller unit root test.
 
 `y` is the time series to be tested, `deterministic` determines the deterministic terms
-(options: `none`, `constant`, `trend`, `squared_trend`) and `lag` the number of lagged
+(options: `:none`, `:constant`, `:trend`, `:squared_trend`) and `lag` the number of lagged
 first-differences included in the test regression, respectively.
 
 Critical values and asymptotic p-values are computed based on response surface regressions
diff --git a/src/common.jl b/src/common.jl
index 67f73e9..0527841 100644
--- a/src/common.jl
+++ b/src/common.jl
@@ -29,13 +29,16 @@
 # sum(n^3-n) where n is the number of ties for each data point
 function tiedrank_adj!(ord::AbstractVector, v::AbstractArray)
     n     = length(v)
+    (n == length(ord)) ||
+        throw(DimensionMismatch("The length of ord ($(length(ord))) differs from v length ($n)"))
     place = sortperm(v)
     tieadj = 0.0
 
     i = 1
-    while i <= n
+    @inbounds while i <= n
         j = i
-        while j + 1 <= n && v[place[i]] == v[place[j + 1]]
+        vi = v[place[i]]
+        while (j + 1 <= n) && (vi == v[place[j + 1]])
             j += 1
         end
 
diff --git a/src/diebold_mariano.jl b/src/diebold_mariano.jl
new file mode 100644
index 0000000..c6bfd14
--- /dev/null
+++ b/src/diebold_mariano.jl
@@ -0,0 +1,76 @@
+# diebold_mariano.jl
+# Diebold-Mariano
+#
+# Copyright (C) 2020   Guilherme Bodin
+#
+# 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 DieboldMarianoTest
+
+struct DieboldMarianoTest <: TTest
+    n::Int         # number of observations
+    df::Int        # degrees of freedom
+    t::Float64     # test statistic
+end
+
+"""
+    DieboldMarianoTest(x::AbstractVector{<:Real}, y::AbstractVector{<:Real}; loss=abs2, h=1)
+
+Perform the modified Diebold-Mariano test proposed by Harvey, Leybourne and Newbold of the null 
+hypothesis that the two methods have the same forecast accuracy. `loss` is the loss function described
+in Diebold, F.X. and Mariano, R.S. (1995) Comparing predictive accuracy. Journal of Business and 
+Economic Statistics, 13, 253-263. and `lookahead` is the number of steps ahead of the forecast.
+
+Implements: [`pvalue`](@ref)
+
+# References
+
+  * Diebold, F.X. and Mariano, R.S. (1995) Comparing predictive accuracy. 
+    Journal of Business and Economic Statistics, 13, 253-263.
+
+  * Harvey, D., Leybourne, S., & Newbold, P. (1997). Testing the equality of prediction 
+    mean squared errors. International Journal of forecasting, 13(2), 281-291.
+  
+"""
+function DieboldMarianoTest(e1::AbstractVector{<:Real}, e2::AbstractVector{<:Real}; 
+                            loss::Function=abs2, lookahead::Integer=1)
+    length(e1) == length(e2) || throw(DimensionMismatch("inputs must have the same length"))
+    n = length(e1)
+    # Calculate the loss diferential series based on the loss function g
+    d = loss.(e1) .- loss.(e2)
+    dm_cov = autocov(d, collect(0:lookahead-1))
+    dm_var = (dm_cov[1] + 2 * sum(dm_cov[2:end]))/n
+    # Statistic from the original Diebold-Mariano test 
+    statistic_dm = mean(d)/sqrt(dm_var)
+    k = sqrt((1 + (1 - 2*lookahead + (lookahead/n)*(lookahead - 1))/n))
+    # Statistic from the modified Diebold-Mariano test proposed by Harvey, Leybourne and Newbold
+    statistic_hln = statistic_dm * k
+    return DieboldMarianoTest(n, n - 1, statistic_hln)
+end
+
+testname(::DieboldMarianoTest) = "Diebold-Mariano test"
+population_param_of_interest(x::DieboldMarianoTest) = ("mean", 0.0, x.t)
+default_tail(test::DieboldMarianoTest) = :both
+
+function show_params(io::IO, x::DieboldMarianoTest, ident)
+    println(io, ident, "number of observations: $(x.n)")
+    println(io, ident, "DM statistic:           $(x.t)")
+    println(io, ident, "degrees of freedom:     $(x.df)")
+end
diff --git a/src/mann_whitney.jl b/src/mann_whitney.jl
index 90114a5..4de0452 100644
--- a/src/mann_whitney.jl
+++ b/src/mann_whitney.jl
@@ -62,10 +62,10 @@ function mwustats(x::AbstractVector{S}, y::AbstractVector{T}) where {S<:Real,T<:
     ny = length(y)
     if nx <= ny
         (ranks, tieadj) = tiedrank_adj([x; y])
-        U = sum(ranks[1:nx]) - nx*(nx+1)/2
+        U = sum(@view ranks[1:nx]) - nx*(nx+1)/2
     else
         (ranks, tieadj) = tiedrank_adj([y; x])
-        U = nx*ny - sum(ranks[1:ny]) + ny*(ny+1)/2
+        U = nx*ny - sum(@view ranks[1:ny]) + ny*(ny+1)/2
     end
     (U, ranks, tieadj, nx, ny, median(x)-median(y))
 end
@@ -110,7 +110,7 @@ function show_params(io::IO, x::ExactMannWhitneyUTest, ident)
     println(io)
     println(io, ident, "Mann-Whitney-U statistic:             ", x.U)
     print(io, ident, "rank sums:                            ")
-    show(io, [sum(x.ranks[1:x.nx]), sum(x.ranks[x.nx+1:end])])
+    show(io, [sum(@view x.ranks[1:x.nx]), sum(@view x.ranks[x.nx+1:end])])
     println(io)
     println(io, ident, "adjustment for ties:                  ", x.tie_adjustment)
 end
@@ -120,16 +120,12 @@ end
 function mwuenumerate(x::ExactMannWhitneyUTest)
     # Get the other U if inverted by mwu_stats
     n = min(x.nx, x.ny)
-    U = x.U
-    if x.ny > x.nx
-        U = x.nx*x.ny -x.U
-    end
+    U = (x.nx >= x.ny ? x.U : x.nx * x.ny - x.U) + n*(n + 1)/2
     le = 0
     gr = 0
     tot = 0
-    k = n*(n+1)/2
     for comb in combinations(x.ranks, n)
-        Up = sum(comb) - k
+        Up = sum(comb)
         tot += 1
         le += Up <= U
         gr += Up >= U
diff --git a/src/permutation.jl b/src/permutation.jl
index edaea76..6e9e034 100644
--- a/src/permutation.jl
+++ b/src/permutation.jl
@@ -1,6 +1,6 @@
 export ExactPermutationTest, ApproximatePermutationTest
 
-function ptstats(x,y,f)
+function ptstats(x,y)
     xy = vcat(x,y)
     rx = 1:length(x)
     ry = (length(xy)-length(y)+1):length(xy)
@@ -20,7 +20,7 @@ that `f(x)` is equal to `f(y)`.  All possible permutations are sampled.
 """
 function ExactPermutationTest(x::AbstractVector{R}, y::AbstractVector{S},
                               f::Function) where {R<:Real,S<:Real}
-    xy, rx, ry = ptstats(x,y,f)
+    xy, rx, ry = ptstats(x,y)
     P = permutations(xy)
     samples = [f(view(p,rx)) - f(view(p,ry)) for p in P]
     PermutationTest(f(x) - f(y), samples)
@@ -35,14 +35,14 @@ permutations are sampled at random.
 """
 function ApproximatePermutationTest(x::AbstractVector{R}, y::AbstractVector{S},
                                     f::Function, n::Int) where {R<:Real,S<:Real}
-    xy, rx, ry = ptstats(x,y,f)
+    xy, rx, ry = ptstats(x,y)
     samples = [(shuffle!(xy); f(view(xy,rx)) - f(view(xy,ry))) for i = 1:n]
     PermutationTest(f(x) - f(y), samples)
 end
 
 function pvalue(apt::PermutationTest; tail=:both)
     if tail == :both
-        count = sum(abs(apt.observation) >= abs(x) for x in apt.samples)
+        count = sum(abs(apt.observation) <= abs(x) for x in apt.samples)
     elseif tail == :left
         count = sum(apt.observation >= x for x in apt.samples)
     elseif tail == :right


Test patch:
diff --git a/src/HypothesisTests.jl b/src/HypothesisTests.jl
index 9ddff68..4a5cf22 100644
--- a/src/HypothesisTests.jl
+++ b/src/HypothesisTests.jl
@@ -190,4 +190,6 @@ include("bartlett.jl")
 include("wald_wolfowitz.jl")
 include("f.jl")
 include("correlation.jl")
+include("diebold_mariano.jl")
+
 end
diff --git a/test/diebold_mariano.jl b/test/diebold_mariano.jl
new file mode 100644
index 0000000..d13834c
--- /dev/null
+++ b/test/diebold_mariano.jl
@@ -0,0 +1,45 @@
+using HypothesisTests, Test
+using HypothesisTests: default_tail
+
+@testset "Diebold-Mariano tests" begin
+    e_ets = [0.002241270, -0.000661906, -1.200002613, -4.000218178, 3.799373533,
+             3.401005111, -1.199692125, -0.800490048, -1.200048999, 3.399847229,
+             7.000708948,  2.200989241, -4.800165569, -1.201043375, 3.400167556, 
+             3.000709085,  0.400267950, -6.800169207, -2.401259050, 1.200112198]
+
+    e_arima = [0.1209998,  8.0622269, -0.9189429, -1.7952436, 5.3431277, 2.0459823,
+               -0.6073261, 0.5516816, -0.7633238, 4.1229472, 5.8785333, 1.2828331,
+               -3.4854233, 1.1971370, 3.1709314, 2.3408457, 0.4962286, -5.9841362,
+                0.1537114, 0.4207951]
+
+    atol = 1e-3
+    dm_test = DieboldMarianoTest(e_ets, e_arima)
+    @test pvalue(dm_test) ≈ 0.8871 atol=atol
+    @test pvalue(dm_test, tail=:right) ≈ 0.5565 atol=atol
+    @test pvalue(dm_test, tail=:left) ≈ 0.4435 atol=atol
+
+    dm_test = DieboldMarianoTest(e_ets, e_arima; lookahead=10)
+    @test pvalue(dm_test) ≈ 0.9362 atol=atol
+    @test pvalue(dm_test, tail=:right) ≈ 0.5319 atol=atol
+    @test pvalue(dm_test, tail=:left) ≈ 0.4681 atol=atol
+
+    dm_test = DieboldMarianoTest(e_ets, e_arima; loss=abs)
+    @test pvalue(dm_test) ≈ 0.7818 atol=atol
+    @test pvalue(dm_test, tail=:right) ≈ 0.3909 atol=atol
+    @test pvalue(dm_test, tail=:left) ≈ 0.6091 atol=atol
+
+    # Random forecast must be worse then ets and arima
+    Random.seed!(1)
+    e_rand = rand(20)
+    dm_test = DieboldMarianoTest(e_ets, e_rand)
+    @test pvalue(dm_test) ≈ 0.0045 atol=atol
+    @test pvalue(dm_test, tail=:right) ≈ 0.0024 atol=atol
+    @test pvalue(dm_test, tail=:left) ≈ 0.9975 atol=atol
+
+    dm_test = DieboldMarianoTest(e_arima, e_rand)
+    @test pvalue(dm_test) ≈ 0.012 atol=atol
+    @test pvalue(dm_test, tail=:right) ≈ 0.0062 atol=atol
+    @test pvalue(dm_test, tail=:left) ≈ 0.9937 atol=atol
+
+    @test_throws DimensionMismatch DieboldMarianoTest(rand(3), rand(4))
+end
diff --git a/test/permutation.jl b/test/permutation.jl
index c232c46..fbc503e 100644
--- a/test/permutation.jl
+++ b/test/permutation.jl
@@ -1,14 +1,14 @@
 using HypothesisTests, Test
 
 @testset "Permutation" begin
-@test isapprox(pvalue(ExactPermutationTest([1,2,3],[4,5,6],mean), tail=:both), 1.0)
-@test isapprox(pvalue(ExactPermutationTest([4,5,6],[1,2,3],mean), tail=:both), 1.0)
+@test isapprox(pvalue(ExactPermutationTest([1,2,3],[4,5,6],mean), tail=:both), 0.1)
+@test isapprox(pvalue(ExactPermutationTest([4,5,6],[1,2,3],mean), tail=:both), 0.1)
 @test isapprox(pvalue(ExactPermutationTest([1,2,3],[4,5,6],mean), tail=:left), 0.05)
 @test isapprox(pvalue(ExactPermutationTest([4,5,6],[1,2,3],mean), tail=:right), 0.05)
 
 Random.seed!(12345)
-@test isapprox(pvalue(ApproximatePermutationTest([1,2,3],[4,5,6],mean,100), tail=:both), 1.0)
-@test isapprox(pvalue(ApproximatePermutationTest([4,5,6],[1,2,3],mean,100), tail=:both), 1.0)
+@test isapprox(pvalue(ApproximatePermutationTest([1,2,3],[4,5,6],mean,100), tail=:both), 0.09)
+@test isapprox(pvalue(ApproximatePermutationTest([4,5,6],[1,2,3],mean,100), tail=:both), 0.04)
 @test isapprox(pvalue(ApproximatePermutationTest([1,2,3],[4,5,6],mean,100), tail=:left), 0.08)
 @test isapprox(pvalue(ApproximatePermutationTest([4,5,6],[1,2,3],mean,100), tail=:right), 0.07)
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 629a098..a2f91a4 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -36,3 +36,4 @@ include("wald_wolfowitz.jl")
 include("wilcoxon.jl")
 include("z.jl")
 include("f.jl")
+include("diebold_mariano.jl")