# R/MDM.test.R In multDM: Multivariate Version of the Diebold-Mariano Test

#### Documented in d_tlossMDM.selectionMDM.testprint.MDM

```### computes multivariate DM test for equal predictive accuracy (EPA)
### of two or more non-nested forecasting models,
### as in R.S. Mariano and D. Preve (2012)
### Statistical tests for multiple forecast comparison,
### Journal of Econometrics 169, 123-130
###
### d - kxP matrix with observations of d_t
###       k+1 - number of considered models
###       P - lenght of time-series
###       d_t = g(y.hat_1t - y.real_t ) - g(y.hat_2t - y.real_t), i.e.,
###       d_jt = g(e_jt) - g(e_j+1,t), j = 1,...,k,
###       where g is the loss function
###
### q - order of the VMA representation of d_t,
###     i.e., a lag length beyond which we are willing
###     to assume that the correlation between d_t and d_t-h
###     is essentially zero
###
### statistic - "S" for the basic version,
###             "Sc"  for the finite-sample correction
###
### the null hypothesis of EPA is E(d_t)=0

.in.MDM.test <- function(d,q,statistic)
{
n <- deparse(substitute(d))

G <- function(d,h)
{
SCM <- matrix(0,nrow(d),nrow(d))
for (t in (h+1):ncol(d))
{
SCM <- SCM + (d[,t,drop=FALSE] - dbar) %*% t(d[,t-h,drop=FALSE] - dbar)
}
SCM <- SCM / ncol(d)

return(SCM)
}

O <- function(d,q)
{
SLRV <- G(d,0)
if (q>0)
{
for (h in 1:q)
{
TEMP <- G(d,h)
SLRV <- SLRV + TEMP + t(TEMP)
}
}

return(SLRV)
}

dbar <- rowMeans(d)
c <- 1 - (1 + 2 * q) / ncol(d) + q * (q+1) / (ncol(d))^2
Om <- O(d,q)
s <- as.vector(dbar) / sqrt(as.vector(abs(diag(Om))) / ncol(d))
S <- ncol(d) * t(dbar) %*% solve(Om) %*% dbar
if (statistic=="Sc")
{
S <- c * S
}

pval <- pchisq(q=S,df=nrow(d),lower.tail=FALSE)

names(S) <- "statistic"
names(q) <- "lag length"

ret <- list(S,q,"Equal predictive accuracy does not hold.",pval,"multivariate Diebold-Mariano test",n)
names(ret) <- c("statistic","parameter","alternative","p.value","method","data.name")
class(ret) <-"htest"
ret <- list(s,ret)
return(ret)
}

MDM.test <- function(realized,evaluated,q,statistic="Sc",loss.type="SE")
{
n <- paste(deparse(substitute(realized)),deparse(substitute(evaluated)),sep=" and ")
l <- loss(realized=realized,evaluated=evaluated,loss.type=loss.type)
d <- d_t(l)
out <- .in.MDM.test(d=d,q=q,statistic=statistic)[[2]]
out\$data.name <- n
return(out)
}

##################################################################

### computes loss function for a set of forecasts
###
### realized - vector of realized values
###
### evaluated - (k+1)xP matrix of forecasts
###               P - length of time-series
###               k+1 - number of models
###
### loss.type - "SE" for squared errors,
###             "AE" for absolute errors,
###              "SPE" for squared proportional error
###                (useful if errors are heteroskedastic)
###                see S. J. Taylor, 2005. Asset Price Dynamics, Volatility, and Prediction,
###                Princeton University Press,
###             "ASE" for absolute scaled error (R. J. Hyndman, A. B. Koehler, 2006,
###               Another Look at Measures of Forecast Accuracy,
###               International Journal of Forecasting volume 22, 679-688,
###             positive numeric value for loss function of type
###              exp(loss*errors)-1-loss*errors
###              (useful when it is more costly to underpredict y than to overpredict)
###              see U. Triacca,  Comparing Predictive Accuracy of Two Forecasts,
###              http://www.phdeconomics.sssup.it/documents/Lesson19.pdf

loss <- function(realized,evaluated,loss.type)
{
e <- evaluated
for (i in 1:nrow(evaluated))
{
e[i,] <- realized - as.vector(evaluated[i,])
}

if (loss.type=="SE")
{
e <- e^2
}
if (loss.type=="AE")
{
e <- abs(e)
}
if (loss.type=="SPE")
{
for (i in 1:nrow(e))
{
e[i,] <- (as.vector(e[i,]) / as.vector(evaluated[i,]))^2
}
}
if (loss.type=="ASE")
{
for (i in 1:nrow(e))
{
e[i,] <- abs(as.vector(e[i,])) / mean(abs(realized-c(NA,realized[-length(realized)]))[-1])
}
e <- e[,-1,drop=FALSE]
}
if (is.numeric(loss.type))
{
e <- exp(loss.type*e)-1-loss.type*e
}

return(e)
}

### computes d_t based on loss function

d_t <- function(e)
{
for (i in 1:(nrow(e)-1))
{
e[i,] <- (as.vector(e[i,]) - as.vector(e[i+1,]))
}
e <- e[-nrow(e),,drop=FALSE]

return(e)
}

MDM.selection <- function(realized,evaluated,q,alpha,statistic="Sc",loss.type="SE")
{
p <- 0
e <- evaluated
models <- (1:nrow(e))
if (is.null(rownames(e))) { rownames(e) <- seq(1:nrow(e)) }
n.mods <- nrow(e)

while(p<alpha && length(models)>1)
{
d <- loss(realized=realized,evaluated=e,loss.type=loss.type)
d <- d_t(d)
mdm <- .in.MDM.test(d=d,q=q,statistic=statistic)
p <- mdm[[2]]\$p.value
if (p<alpha)
{
j <- which.max(abs(mdm[[1]]))
if (mdm[[1]][j]>0)
{
models <- models[-j]
j.drop <- j
}
else
{
models <- models[-(j+1)]
j.drop <- j+1
}
e <- e[-j.drop,,drop=FALSE]
models <- (1:nrow(e))
}
}

n.mods <- n.mods - nrow(e)
ret <- matrix(NA,ncol=3,nrow=nrow(e))
rownames(ret) <- rownames(e)

ret[-nrow(ret),2] <- mdm[[1]]
ret[-nrow(ret),1] <- rank(abs(ret[-nrow(ret),2]))
ret[,3] <- rowMeans(loss(realized=realized,evaluated=e,loss.type=loss.type))

colnames(ret) <- c("Rank","s","Mean loss")

ret <- list(ret,as.numeric(p),alpha,n.mods)
names(ret) <- c("outcomes","p.value","alpha","eliminated")
class(ret) <- "MDM"
return(ret)
}

print.MDM <- function(x, ...)
{
if (x[[2]]>x[[3]])
{
cat("####################################################")
cat("\n")
cat("Models with outstanding predictive ability:")
cat("\n")
cat("\n")
print(round(x[[1]],digits=4),na.print="")
cat("\n")
cat("p-value: ")
cat(round(as.numeric(x[[2]]),digits=4))
cat("\n")
cat("\n")
cat("Number of eliminated models: ",x[[4]])
cat("\n")
cat("####################################################")
}
else
{
cat("No models with outstanding predictive ability were found.")
}
}
```

## Try the multDM package in your browser

Any scripts or data that you put into this service are public.

multDM documentation built on June 9, 2022, 5:06 p.m.