diebold.mariano.test <- function(x, alternative = c("two.sided", "less", "greater"), k)
{
    if (NCOL(x) > 1) 
        stop("x is not a vector or univariate time series")
    if (any(is.na(x))) 
        stop("NAs in x")
    alternative <- match.arg(alternative)
    DNAME <- deparse(substitute(x))
    n <- NROW(x)
    cv <- acf(x, lag.max=k, type="covariance", plot=FALSE)$acf[,,1]
    eps <- 1.0e-8
    vr <- max(eps, sum(c(cv[1], 2*cv[-1])) / n)
    STATISTIC <- mean(x) / sqrt(vr)
    names(STATISTIC) <- "Standard Normal"
    METHOD <- "Diebold-Mariano Test"
    if (alternative == "two.sided") 
        PVAL <- 2 * pnorm(-abs(STATISTIC))
    else if (alternative == "less") 
        PVAL <- pnorm(STATISTIC)
    else if (alternative == "greater") 
        PVAL <- pnorm(STATISTIC, lower.tail = FALSE)
    PARAMETER <- k
    names(PARAMETER) <- "Truncation lag"
    structure(list(statistic = STATISTIC, parameter = PARAMETER, alternative = alternative, 
                   p.value = PVAL, method = METHOD, data.name = DNAME), 
              class = "htest")
}


g <- function(x)
{
    abs(x)
}

e1 <- rnorm(500)
e2 <- rnorm(500)

diebold.mariano.test(g(e1)-g(e2), k = 3)

e1 <- rnorm(500, sd=1)
e2 <- rnorm(500, sd=1.3)

diebold.mariano.test(g(e1)-g(e2), k = 3)

