## Out of sample forecast of 10-year forward rate
## Diebold-Mariano Test

setwd("D:/HM/JMCB/UK2/UK_oostest")

myDMtest <- function (e1, e2, alternative = c("two.sided", "less", "greater"),
                      h = 1, power = 2, lrvar = c("HH", "NW")) {
  ## Diebold-Mariano test
  ## based on dm.test from forecast package
  alternative <- match.arg(alternative)
  d <- c(abs(e1))^power - c(abs(e2))^power
  lrvar <- match.arg(lrvar)
  if (lrvar == "HH") {
    ## Hansen-Hodrick
    d.cov <- acf(d, na.action = na.omit, lag.max = h - 1, type = "covariance",
                 plot = FALSE)$acf[, , 1]
    dv <- sum(c(d.cov[1], 2 * d.cov[-1]))/length(d)
  } else {
    ## Newey-West
    L <- floor(1.5*(h-1)) # need to use 1.5 x more lags, instead of L = h-1
    d.cov <- acf(d, na.action = na.omit, lag.max = L, type = "covariance",
                 plot = FALSE)$acf[, , 1]
    dv <- sum(c(d.cov[1], 2*(1 - (1:L)/(L+1)) * d.cov[-1]))/length(d)
    ## this is identical to NeweyWest(lm(d ~ 1), lag=L, prewhite=FALSE)
  }
  if (dv > 0)
    STATISTIC <- mean(d, na.rm = TRUE)/sqrt(dv)
  else if (h == 1)
    stop("Variance of DM statistic is zero")
  else {
    cat("Variance for h =", h, "is negative, using horizon h=", h-1, "\n")
    return(myDMtest(e1, e2, alternative, h = h-1, power))
    ## warning("Variance is negative, using NW")
    ## return(myDMtest(e1, e2, alternative, h, power, lrvar="NW"))
  }
  n <- length(d)
  k <- ((n + 1 - 2 * h + (h/n) * (h - 1))/n)^(1/2)
  STATISTIC <- STATISTIC * k
  names(STATISTIC) <- "DM"
  if (alternative == "two.sided")
    PVAL <- 2 * pt(-abs(STATISTIC), df = n - 1)
  else if (alternative == "less")
    PVAL <- pt(STATISTIC, df = n - 1)
  else if (alternative == "greater")
    PVAL <- pt(STATISTIC, df = n - 1, lower.tail = FALSE)
  PARAMETER <- c(h, power)
  names(PARAMETER) <- c("Forecast horizon", "Loss function power")
  structure(list(statistic = STATISTIC, parameter = PARAMETER,
                 alternative = alternative, p.value = PVAL, method = "Diebold-Mariano Test",
                 data.name = c(deparse(substitute(e1)), deparse(substitute(e2)))),
            class = "htest")
}


###########################################################
############## Monthly Forecast Comparison ################
###########################################################

# 6-month forecast
# sdtsm.M6 <- as.numeric(unlist(read.csv('eseM6.csv', header = F)))
sdtsm.M6 <- as.numeric(unlist(read.csv('eseM6_rt.csv', header = F)))
rw.M6 <- as.numeric(unlist(read.csv('rw6.csv', header = F)))
jsz.M6 <- as.numeric(unlist(read.csv('jsz6.csv', header = F)))
dmtest.M6 <- c(NA,NA)
dmtest.M6[1] <-  myDMtest(sdtsm.M6, rw.M6, alternative="less", h=6, power=2)$p.value
dmtest.M6[2] <-  myDMtest(sdtsm.M6, jsz.M6, alternative="less", h=6, power=2)$p.value

sqrt((t(rw.M6)%*%rw.M6)/length(rw.M6))
sqrt((t(jsz.M6)%*%jsz.M6)/length(jsz.M6))
sqrt((t(sdtsm.M6)%*%sdtsm.M6)/length(sdtsm.M6))


# 10-month forecast
# sdtsm.M10 <- as.numeric(unlist(read.csv('eseM10.csv', header = F)))
sdtsm.M10 <- as.numeric(unlist(read.csv('eseM10_rt.csv', header = F)))
rw.M10 <- as.numeric(unlist(read.csv('rw10.csv', header = F)))
jsz.M10 <- as.numeric(unlist(read.csv('jsz10.csv', header = F)))
dmtest.M10 <- c(NA,NA)
dmtest.M10[1] <-  myDMtest(sdtsm.M10, rw.M10, alternative="less", h=10, power=2)$p.value
dmtest.M10[2] <-  myDMtest(sdtsm.M10, jsz.M10, alternative="less", h=10, power=2)$p.value

sqrt((t(rw.M10)%*%rw.M10)/length(rw.M10))
sqrt((t(jsz.M10)%*%jsz.M10)/length(jsz.M10))
sqrt((t(sdtsm.M10)%*%sdtsm.M10)/length(sdtsm.M10))


# 20-month forecast
# sdtsm.M20 <- as.numeric(unlist(read.csv('eseM20.csv', header = F)))
sdtsm.M20 <- as.numeric(unlist(read.csv('eseM20_rt.csv', header = F)))
rw.M20 <- as.numeric(unlist(read.csv('rw20.csv', header = F)))
jsz.M20 <- as.numeric(unlist(read.csv('jsz20.csv', header = F)))
dmtest.M20 <- c(NA,NA)
dmtest.M20[1] <-  myDMtest(sdtsm.M20, rw.M20, alternative="less", h=20, power=2)$p.value
dmtest.M20[2] <-  myDMtest(sdtsm.M20, jsz.M20, alternative="less", h=20, power=2)$p.value

sqrt((t(rw.M20)%*%rw.M20)/length(rw.M20))
sqrt((t(jsz.M20)%*%jsz.M20)/length(jsz.M20))
sqrt((t(sdtsm.M20)%*%sdtsm.M20)/length(sdtsm.M20))


# 30-month forecast
# sdtsm.M30 <- as.numeric(unlist(read.csv('eseM30.csv', header = F)))
sdtsm.M30 <- as.numeric(unlist(read.csv('eseM30_rt.csv', header = F)))
rw.M30 <- as.numeric(unlist(read.csv('rw30.csv', header = F)))
jsz.M30 <- as.numeric(unlist(read.csv('jsz30.csv', header = F)))
dmtest.M30 <- c(NA,NA)
dmtest.M30[1] <-  myDMtest(sdtsm.M30, rw.M30, alternative="less", h=30, power=2)$p.value
dmtest.M30[2] <-  myDMtest(sdtsm.M30, jsz.M30, alternative="less", h=30, power=2)$p.value

sqrt((t(rw.M30)%*%rw.M30)/length(rw.M30))
sqrt((t(jsz.M30)%*%jsz.M30)/length(jsz.M30))
sqrt((t(sdtsm.M30)%*%sdtsm.M30)/length(sdtsm.M30))


# 40-month forecast
# sdtsm.M40 <- as.numeric(unlist(read.csv('eseM40.csv', header = F)))
sdtsm.M40 <- as.numeric(unlist(read.csv('eseM40_rt.csv', header = F)))
rw.M40 <- as.numeric(unlist(read.csv('rw40.csv', header = F)))
jsz.M40 <- as.numeric(unlist(read.csv('jsz40.csv', header = F)))
dmtest.M40 <- c(NA,NA)
dmtest.M40[1] <-  myDMtest(sdtsm.M40, rw.M40, alternative="less", h=40, power=2)$p.value
dmtest.M40[2] <-  myDMtest(sdtsm.M40, jsz.M40, alternative="less", h=40, power=2)$p.value

sqrt((t(rw.M40)%*%rw.M40)/length(rw.M40))
sqrt((t(jsz.M40)%*%jsz.M40)/length(jsz.M40))
sqrt((t(sdtsm.M40)%*%sdtsm.M40)/length(sdtsm.M40))





