Add Diebold-Mariano test

[guilhermebodin]

Closes the issue #179 .

All tests were based on the Diebold-Mariano implementation on the forecast package of R. The script is the following.

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 TTest subtypes here:

HypothesisTests.jl/src/t.jl

Line 31 in 7f21943

pvalue(x::TTest; tail=:both) = pvalue(TDist(x.df), x.t; tail=tail)

[PaulSoderlind]

I would recommend adding the Clark-West test (for nested models), where the test is based on

2*residuals(f1.out).*(residuals(f1.out)-residuals(f2.out))

assuming that f1 is the smaller model and that residuals(fx.out) is a vector of residuals. With this, the usual tests can be applied. See, for instance, equation (5.40) of the Empirical Asset Pricing notes at my homepage.

EDIT: I believe it is well recognized that applying a DM test for nested models is a non-trivial mistake and that either (a) a bootstrap; (b) moving data window (rather than an expanding ones) or (c) a CW test are the possible ways to handle this case. In case it's too late for introducing a CW test in this PR, please let me know so I'll do it in a new PR.

[PaulSoderlind]

please see my updated comment (above)

[guilhermebodin]

@PaulSoderlind I think it is better to make the Clark-West test in another PR.

[PaulSoderlind]

@guilhermebodin

how do we move on with the Clark-West? I am happy to contribute, but I must admit that I am not really an expert on the API of this package.

One idea would be to just have it as an option in the DM test. This would be easy to both implement and make a PR of. What do you think?

[guilhermebodin]

@PaulSoderlind I think the easiest way to implement this is basically to copy this file in a new one named clark_west.jl and change the expression for d for the one you commented. Moving into a new struct will be more clear to the user and will have more clear documentation.

After making the new file you should add some tests and add it to the time series tests section on the documentation.

[PaulSoderlind]

@guilhermebodin OK, sounds right. Maybe quicker if you take the lead? (I could review.)

[guilhermebodin]

@ararslan is it good to merge?

[mschauer]

Looks quite good to me. @PaulSoderlind, I am unsure how to read your comment, do you think this is good to go as well?

[PaulSoderlind]

My comment referred to the idea of adding another PR for a Clark-West test. Anyhow I could review this PR. Just a stupid question before a try to do that: how do I see the latest version of the PR (after all the updates)? It's a bit of a hassle to wade through all the diffs.

…

On Wed, 20 May 2020 at 22:26, Moritz Schauer ***@***.***> wrote: Looks quite good to me. @PaulSoderlind <https://github.com/PaulSoderlind>, I am unsure how to read your comment, do you think this is good to go as well? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#180 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACHRTVL57YUHO5HD3VJDCZDRSQ4H5ANCNFSM4KUBXAIA> .

[mschauer]

The "files changed" tab above shows the accumulated changes,

https://github.com/JuliaStats/HypothesisTests.jl/pull/180/files

[PaulSoderlind]

On the src/diebold_mariano.jl: the DM test looks fine (and gives the same result as I get from my own home-brewed code). I don't understand the TTest type so I cannot comment on that.