Проблема в том, что gdk возвращает отрицательное число, а sqrt с отрицательным числом равно NaN. Если мы добавим несколько print операторов в функцию и запустим, станет ясно
library(multDM)
DM.test1(forecast1,forecast2,ts,loss.type="SE",h=4,c=TRUE)
[1] "gdk : 885835705046.250000"
[1] "gdk/T : 68141208080.480766"
[1] "sqrt(gdk/T): 261038.709927"
Diebold-Mariano test
data: forecast1 and forecast2 and ts
statistic = -44.562, forecast horizon = 4, p-value = 1.066e-14
alternative hypothesis: Forecast f1 and f2 have different accuracy.
с h= 13
DM.test1(forecast1,forecast2,ts,loss.type="SE",h=13,c=TRUE)
[1] "gdk : 317749091646343.500000"
[1] "gdk/T : 24442237818949.500000"
[1] "sqrt(gdk/T): 4943909.163703"
Diebold-Mariano test
data: forecast1 and forecast2 and ts
statistic = -6.655, forecast horizon = 13, p-value = 2.343e-05
alternative hypothesis: Forecast f1 and f2 have different accuracy.
с h = 5
DM.test1(forecast1,forecast2,ts,loss.type="SE",h=5,c=TRUE)
[1] "gdk : -201824682703276.000000"
[1] "gdk/T : -15524975592559.691406"
[1] "sqrt(gdk/T): NaN"
Diebold-Mariano test
data: forecast1 and forecast2 and ts
statistic = NaN, forecast horizon = 5, p-value = NA
alternative hypothesis: Forecast f1 and f2 have different accuracy.
Warning messages:
1: In sqrt(gdk/T) : NaNs produced
2: In sqrt(gdk/T) : NaNs produced
-функция с оператором печати
DM.test1 <- function (f1, f2, y, loss.type = "SE", h = 1, c = FALSE, H1 = "same")
{
n <- paste(deparse(substitute(f1)), deparse(substitute(f2)),
deparse(substitute(y)), sep = " and ")
n1 <- deparse(substitute(f1))
n2 <- deparse(substitute(f2))
e1 <- f1 - y
e2 <- f2 - y
if (loss.type == "SE") {
g1 <- e1^2
g2 <- e2^2
}
if (loss.type == "AE") {
g1 <- abs(e1)
g2 <- abs(e2)
}
if (loss.type == "SPE") {
g1 <- ((y - f1)/f1)^2
g2 <- ((y - f2)/f2)^2
}
if (loss.type == "ASE") {
g1 <- abs(e1[-1])/mean(abs((y - c(NA, y[-length(y)]))[-1]))
g2 <- abs(e2[-1])/mean(abs((y - c(NA, y[-length(y)]))[-1]))
}
if (is.numeric(loss.type)) {
g1 <- exp(loss.type * e1) - 1 - loss.type * e1
g2 <- exp(loss.type * e2) - 1 - loss.type * e2
}
d <- g1 - g2
T <- length(d)
dbar <- mean(d)
gammahat <- function(k) {
temp1 <- d - dbar
temp2 <- rep(NA, abs(k))
temp2 <- c(temp2, temp1)
temp2 <- temp2[1:T]
temp2 <- temp2 - dbar
temp <- temp1 * temp2
temp <- temp[(1 + abs(k)):T]
temp <- sum(temp)/T
return(temp)
}
if (h > 1) {
gdk <- lapply(seq(from = 1, to = h - 1, by = 1), gammahat)
gdk <- unlist(gdk)
}
else {
gdk <- 0
}
gdk <- gammahat(0) + 2 * sum(gdk)
print(sprintf("gdk : %f", gdk))
print(sprintf("gdk/T : %f", gdk/T))
print(sprintf("sqrt(gdk/T): %f",sqrt(gdk/T) ))
DM <- dbar/sqrt(gdk/T)
if (H1 == "same") {
pval <- 2 * min(pnorm(DM, lower.tail = FALSE), 1 - pnorm(DM,
lower.tail = FALSE))
}
if (H1 == "less") {
pval <- pnorm(DM, lower.tail = FALSE)
}
if (H1 == "more") {
pval <- 1 - pnorm(DM, lower.tail = FALSE)
}
if (c) {
DM <- DM * sqrt((T + 1 - 2 * h + h * (h - 1))/T)
if (H1 == "same") {
pval <- 2 * min(pt(q = DM, df = T - 1, lower.tail = FALSE),
1 - pt(q = DM, df = T - 1, lower.tail = FALSE))
}
if (H1 == "less") {
pval <- pt(q = DM, df = T - 1, lower.tail = FALSE)
}
if (H1 == "more") {
pval <- 1 - pt(q = DM, df = T - 1, lower.tail = FALSE)
}
}
names(DM) <- "statistic"
names(h) <- "forecast horizon"
if (H1 == "same") {
alt <- "Forecast f1 and f2 have different accuracy."
}
if (H1 == "less") {
alt <- "Forecast f1 is less accurate than f2."
}
if (H1 == "more") {
alt <- "Forecast f1 is more accurate than f2."
}
ret <- list(DM, h, paste(alt), pval, "Diebold-Mariano test",
n)
names(ret) <- c("statistic", "parameter", "alternative",
"p.value", "method", "data.name")
class(ret) <- "htest"
return(ret)
}
data
ts <- c(14500, 8300, 10900, 11400, 5600, 16500, 11000, 18500, 11800,
17900, 12400, 9100, 10600)
forecast1 <- c(11965.53, 11918.98, 11761.95, 11925.01, 11751.1, 11671.07,
11688.99, 11790.13, 11888.66, 11894.13, 11850.57, 11681.65, 11488.7
)
forecast2 <- c(12100, 10900, 15300, 13100, 16600, 24900, 11100, 17400, 9600,
11400, 8100, 11600, 11000)