## Out of sample forecast of forward rate using ese model
## 1971Q4 - 2018Q1 
# setwd("D:/trends_replication/oos/US2")
# setwd("D:/HM/JMCB/BR_oos/oos/US2")
setwd("D:/HM/JMCB/US4/oostest/US2")
library(KFAS)
library(MCMCpack)
source("R/dtsm_fns.R")
source("R/data_fns.R")
source("R/util_fns.R")

# load("myese_616.RData")
# load("ese_616_GSW_Quartermean.RData")
# load("results/myese_616.RData")
load("results/myeseLW2_616.RData")
LW.qend <- read.csv("data/LW_forward_quarterend.csv",header = F)
LW.qmean <- read.csv("data/LW_forward_quartermean.csv",header = F)
LW.qend[,2:18] <- LW.qend[,2:18]/100
LW.qmean[,2:18] <- LW.qmean[,2:18]/100


N <- 3
ltri <- lower.tri(matrix(0, N, N), diag=TRUE)
T <- nrow(Y)

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")
}

ssm <- function(Y, pars) {
  ## state vector: 5 state variables
  ##  (1) constant, (2) i*, (3-5) Ptilde(t)
  N <- length(pars$lamQ)
  J <- ncol(Y)
  ## measurement equation
  Z <- cbind(t(pars$AcP + crossprod(pars$Pbar, pars$BcP)), # constant
             t(crossprod(pars$gamma, pars$BcP)),         # tau
             t(pars$BcP))                                # Ptilde(t)
  H <- diag(J)*pars$sige2
  ## transition equation
  T <- diag(2+N)
  T[3:5,3:5] <- pars$Phi
  R <- matrix(0, 2+N, N+1); R[2:(N+2),] <- diag(N+1)
  Q <- matrix(0, 1+N, 1+N); Q[1,1] <- pars$sigtau2; Q[-1,-1] <- pars$Sigma.tilde
  ## initial conditions
  a1 <- c(1, 0.02, rep(0, 3))
  P1 <- diag(c(0, rep(1, 4)))
  ## model
  # Y[zlb.ind==1] <- NA  #  yields that are constrained at ZLB treated as unobserved
  SSModel(Y ~ -1 + SSMcustom(Z=Z, T=T, R=R, Q=Q, a1=a1, P1=P1), H=H)
}

Sigma_temp <- matrix(NA,3,3)
Sigma_temp[ltri] <- colMeans(Sigma)
Sigma_temp[1,2] <- Sigma_temp[2,1]
Sigma_temp[1,3] <- Sigma_temp[3,1]
Sigma_temp[2,3] <- Sigma_temp[3,2]

Sigmatilde_temp <- matrix(NA,3,3)
Sigmatilde_temp[ltri] <- colMeans(Sigma.tilde)
Sigmatilde_temp[1,2] <- Sigmatilde_temp[2,1]
Sigmatilde_temp[1,3] <- Sigmatilde_temp[3,1]
Sigmatilde_temp[2,3] <- Sigmatilde_temp[3,2]

Phi_temp <- colMeans(Phi)
Phi_temp <- matrix(Phi_temp,3,3)


pars.mean <- within(list(),{
  kinfQ <- mean(kinfQ)
  lamQ <- colMeans(lamQ)
  Sigma <- Sigma_temp
  gamma <- colMeans(gamma)
  Pbar <- colMeans(Pbar)
  Phi <- Phi_temp
  sige2 <- mean(sige2)
  sigtau2 <- mean(sigtau2)
  Sigma.tilde <- Sigmatilde_temp
})

# loadings for 10y and 11y yields
# mats <- c(0.25, 0.5, 1:15)
# W <- 1*cbind(mats==0.25, mats==2, mats==10)
# WN <- t(W)
loads.yields <- jsz.loadings(WN, diag(pars.mean$lamQ-1), pars.mean$kinfQ, pars.mean$Sigma, mats, dt)
pars.mean$AcP <- loads.yields$AcP
pars.mean$BcP <- loads.yields$BcP
# Loading for 10y and 11y log-price
A10 <- -10*loads.yields$AcP[12] 
A11 <- -11*loads.yields$AcP[13]
B10 <- -10*loads.yields$BcP[,12]
B11 <- -11*loads.yields$BcP[,13]

# Loading for 1y and 2y log-price
A1 <- -1*loads.yields$AcP[3] 
A2 <- -2*loads.yields$AcP[4]
B1 <- -1*loads.yields$BcP[,3]
B2 <- -2*loads.yields$BcP[,4]

#############################################################
# We conduct 1, 2, 4, 8-Quarter forecast
# Forecast starts at 1976Q3 when 20 observations are available
# ends at 2016Q1. There are 159 observations in total

ep <- T-8
fhat <- matrix(NA,ep,1)

## 1-Quarter Forecast
fhat <- matrix(NA,T,1)
f <- matrix(NA,T,1)
err <- matrix(NA,T,1)
for (i in 20:ep){
  sim <- simulateSSM(ssm(Y[1:i,], pars.mean), type="states") ## (N+1) x i
  sim.tau <- sim[,2,1]  # tau_exp = sim.tau
  sim.Ptilde <- sim[,3:5,1]
  Ptilde_exp <- sim.Ptilde %*% t(pars.mean$Phi)
  cp_exp <- rep(1, i) %*% rbind(pars.mean$Pbar) + sim.tau %*% rbind(pars.mean$gamma) + Ptilde_exp
  f[i+1] <- 11*Y[i+1,13] - 10*Y[i+1,12]
  fhat[i+1] <- A10 - A11 +  cp_exp[i,] %*% (B10-B11)
  err[i+1] <- f[i+1] - fhat[i+1]
}
err1 <- err[21:(ep+1)]*100
rmse1 <- sqrt((t(err1)%*%err1)/(ep+1-20))


## 2-Quarter Forecast
fhat <- matrix(NA,T,1)
f <- matrix(NA,T,1)
err <- matrix(NA,T,1)
for (i in 20:ep){
  sim <- simulateSSM(ssm(Y[1:i,], pars.mean), type="states") ## (N+1) x i
  sim.tau <- sim[,2,1]  # tau_exp = sim.tau
  sim.Ptilde <- sim[,3:5,1]
  Ptilde_exp <- sim.Ptilde
  for (j in 1:2){
    Ptilde_exp <- Ptilde_exp %*% t(pars.mean$Phi)
  }
  cp_exp <- rep(1, i) %*% rbind(pars.mean$Pbar) + sim.tau %*% rbind(pars.mean$gamma) + Ptilde_exp
  f[i+2] <- 11*Y[i+2,13] - 10*Y[i+2,12]
  fhat[i+2] <- A10 - A11 +  cp_exp[i,] %*% (B10-B11)
  err[i+2] <- f[i+2] - fhat[i+2]
}
err2 <- err[22:(ep+2)]*100
rmse2 <- sqrt((t(err2)%*%err2)/length(err2))


## 4-Quarter Forecast
fhat <- matrix(NA,T,1)
f <- matrix(NA,T,1)
err <- matrix(NA,T,1)
for (i in 20:ep){
  sim <- simulateSSM(ssm(Y[1:i,], pars.mean), type="states") ## (N+1) x i
  sim.tau <- sim[,2,1]  # tau_exp = sim.tau
  sim.Ptilde <- sim[,3:5,1]
  Ptilde_exp <- sim.Ptilde
  for (j in 1:4){
    Ptilde_exp <- Ptilde_exp %*% t(pars.mean$Phi)
  }
  cp_exp <- rep(1, i) %*% rbind(pars.mean$Pbar) + sim.tau %*% rbind(pars.mean$gamma) + Ptilde_exp
  f[i+4] <- 11*Y[i+4,13] - 10*Y[i+4,12]
  fhat[i+4] <- A10 - A11 +  cp_exp[i,] %*% (B10-B11)
  err[i+4] <- f[i+4] - fhat[i+4]
}
err4 <- err[24:(ep+4)]*100
rmse4 <- sqrt((t(err4)%*%err4)/length(err4))


## 8-Quarter Forecast
fhat <- matrix(NA,T,1)
f <- matrix(NA,T,1)
err <- matrix(NA,T,1)
for (i in 20:ep){
  sim <- simulateSSM(ssm(Y[1:i,], pars.mean), type="states") ## (N+1) x i
  sim.tau <- sim[,2,1]  # tau_exp = sim.tau
  sim.Ptilde <- sim[,3:5,1]
  Ptilde_exp <- sim.Ptilde
  for (j in 1:8){
    Ptilde_exp <- Ptilde_exp %*% t(pars.mean$Phi)
  }
  cp_exp <- rep(1, i) %*% rbind(pars.mean$Pbar) + sim.tau %*% rbind(pars.mean$gamma) + Ptilde_exp
  f[i+8] <- 11*Y[i+8,13] - 10*Y[i+8,12]
  fhat[i+8] <- A10 - A11 +  cp_exp[i,] %*% (B10-B11)
  err[i+8] <- f[i+8] - fhat[i+8]
}
err8 <- err[28:(ep+8)]*100
rmse8 <- sqrt((t(err8)%*%err8)/length(err8))


## Compare SDTSM and ESE
sdtsm.Q1 <- as.numeric(unlist(read.csv('sdtsm_Q1_full.csv', header = F)))
dmtest_Q1 <- myDMtest(sdtsm.Q1, err1, alternative="less", h=1, power=2)

sdtsm.Q2 <- as.numeric(unlist(read.csv('sdtsm_Q2_full.csv', header = F)))
dmtest_Q2 <- myDMtest(sdtsm.Q2, err2, alternative="less", h=2, power=2)

sdtsm.Q4 <- as.numeric(unlist(read.csv('sdtsm_Q4_full.csv', header = F)))
dmtest_Q4 <- myDMtest(sdtsm.Q4, err4, alternative="less", h=4, power=2)

sdtsm.Q8 <- as.numeric(unlist(read.csv('sdtsm_Q8_full.csv', header = F)))
dmtest_Q8 <- myDMtest(sdtsm.Q8, err8, alternative="less", h=8, power=2)

save.image("./results/myoos_Quarterend_short.RData")


#############################################################
# We conduct 4, 10, 20, 30, 40-Quarter forecast
# Forecast starts at 1976Q3 when 20 observations are available
# ends at 2008Q1. There are 127 observations in total
set.seed(616)
ep <- T-40
# num_forecast <- 127

# Y <- as.matrix(Y_GSW2)/100
# Y <- Y[1:186,] 


## 4-Quarter Forecast
fhat <- matrix(NA,T,1)
f <- matrix(NA,T,1)
err <- matrix(NA,T,1)

for (i in 20:ep){
  sim <- simulateSSM(ssm(Y[1:i,], pars.mean), type="states") ## (N+1) x i
  sim.tau <- sim[,2,1]  # tau_exp = sim.tau
  sim.Ptilde <- sim[,3:5,1]
  Ptilde_exp <- sim.Ptilde
  for (j in 1:4){
    Ptilde_exp <- Ptilde_exp %*% t(pars.mean$Phi)
  }
  cp_exp <- rep(1, i) %*% rbind(pars.mean$Pbar) + sim.tau %*% rbind(pars.mean$gamma) + Ptilde_exp
  # f[i+4] <- 11*Y[i+4,13] - 10*Y[i+4,12]
  f[i+4] <- LW.qmean[i+4,13]
  fhat[i+4] <- A10 - A11 +  cp_exp[i,] %*% (B10-B11)
  err[i+4] <- f[i+4] - fhat[i+4]
}

# err4 <- err[24:(ep+4)]*100
err4 <- err[25:(ep+4)]*100   # 1977Q4 - 2009Q1
rmse4 <- sqrt((t(err4)%*%err4)/length(err4))


## 10-Quarter Forecast
fhat <- matrix(NA,T,1)
f <- matrix(NA,T,1)
err <- matrix(NA,T,1)

for (i in 20:ep){
  sim <- simulateSSM(ssm(Y[1:i,], pars.mean), type="states") ## (N+1) x i
  sim.tau <- sim[,2,1]  # tau_exp = sim.tau
  sim.Ptilde <- sim[,3:5,1]
  Ptilde_exp <- sim.Ptilde
  for (j in 1:10){
    Ptilde_exp <- Ptilde_exp %*% t(pars.mean$Phi)
  }
  cp_exp <- rep(1, i) %*% rbind(pars.mean$Pbar) + sim.tau %*% rbind(pars.mean$gamma) + Ptilde_exp
 # f[i+10] <- 11*Y[i+10,13] - 10*Y[i+10,12]
  f[i+10] <- LW.qmean[i+10,13]
  fhat[i+10] <- A10 - A11 +  cp_exp[i,] %*% (B10-B11)
  err[i+10] <- f[i+10] - fhat[i+10]
}

# err10 <- err[30:(ep+10)]*100
err10 <- err[31:(ep+10)]*100
rmse10 <- sqrt((t(err10)%*%err10)/length(err10))


## 20-Quarter Forecast
fhat <- matrix(NA,T,1)
f <- matrix(NA,T,1)
err <- matrix(NA,T,1)

for (i in 20:ep){
  sim <- simulateSSM(ssm(Y[1:i,], pars.mean), type="states") ## (N+1) x i
  sim.tau <- sim[,2,1]  # tau_exp = sim.tau
  sim.Ptilde <- sim[,3:5,1]
  Ptilde_exp <- sim.Ptilde
  for (j in 1:20){
    Ptilde_exp <- Ptilde_exp %*% t(pars.mean$Phi)
  }
  cp_exp <- rep(1, i) %*% rbind(pars.mean$Pbar) + sim.tau %*% rbind(pars.mean$gamma) + Ptilde_exp
  # f[i+20] <- 11*Y[i+20,13] - 10*Y[i+20,12]
  f[i+20] <- LW.qmean[i+20,13]
  fhat[i+20] <- A10 - A11 +  cp_exp[i,] %*% (B10-B11)
  err[i+20] <- f[i+20] - fhat[i+20]
}

# err20 <- err[40:(ep+20)]*100
err20 <- err[41:(ep+20)]*100
rmse20 <- sqrt((t(err20)%*%err20)/length(err20))


## 30-Quarter Forecast
fhat <- matrix(NA,T,1)
f <- matrix(NA,T,1)
err <- matrix(NA,T,1)

for (i in 20:ep){
  sim <- simulateSSM(ssm(Y[1:i,], pars.mean), type="states") ## (N+1) x i
  sim.tau <- sim[,2,1]  # tau_exp = sim.tau
  sim.Ptilde <- sim[,3:5,1]
  Ptilde_exp <- sim.Ptilde
  for (j in 1:30){
    Ptilde_exp <- Ptilde_exp %*% t(pars.mean$Phi)
  }
  cp_exp <- rep(1, i) %*% rbind(pars.mean$Pbar) + sim.tau %*% rbind(pars.mean$gamma) + Ptilde_exp
  # f[i+30] <- 11*Y[i+30,13] - 10*Y[i+30,12]
  f[i+30] <- LW.qmean[i+30,13]
  fhat[i+30] <- A10 - A11 +  cp_exp[i,] %*% (B10-B11)
  err[i+30] <- f[i+30] - fhat[i+30]
}

# err30 <- err[50:(ep+30)]*100
err30 <- err[51:(ep+30)]*100
rmse30 <- sqrt((t(err30)%*%err30)/length(err30))


## 40-Quarter Forecast
fhat <- matrix(NA,T,1)
f <- matrix(NA,T,1)
err <- matrix(NA,T,1)

for (i in 20:ep){
  sim <- simulateSSM(ssm(Y[1:i,], pars.mean), type="states") ## (N+1) x i
  sim.tau <- sim[,2,1]  # tau_exp = sim.tau
  sim.Ptilde <- sim[,3:5,1]
  Ptilde_exp <- sim.Ptilde
  for (j in 1:40){
    Ptilde_exp <- Ptilde_exp %*% t(pars.mean$Phi)
  }
  cp_exp <- rep(1, i) %*% rbind(pars.mean$Pbar) + sim.tau %*% rbind(pars.mean$gamma) + Ptilde_exp
  # f[i+40] <- 11*Y[i+40,13] - 10*Y[i+40,12]
  f[i+40] <- LW.qmean[i+40,13]
  fhat[i+40] <- A10 - A11 +  cp_exp[i,] %*% (B10-B11)
  err[i+40] <- f[i+40] - fhat[i+40]
}

# err40 <- err[60:(ep+40)]*100
err40 <- err[61:(ep+40)]*100
rmse40 <- sqrt((t(err40)%*%err40)/length(err40))


## Compare SDTSM and ESE
sdtsm.Q4 <- as.numeric(unlist(read.csv('eseQ4.csv', header = F)))
# sdtsm.Q4 <- as.numeric(unlist(read.csv('data/eseQ4.csv', header = F)))
# sdtsm.Q4 <- as.numeric(unlist(read.csv('eseQ4_start.csv', header = F)))
# sdtsm.Q4 <- as.numeric(unlist(read.csv('eseQ4_end.csv', header = F)))
dmtest4 <- myDMtest(sdtsm.Q4, err4, alternative="less", h=4, power=2)

sdtsm.Q10 <- as.numeric(unlist(read.csv('eseQ10.csv', header = F)))
# sdtsm.Q10 <- as.numeric(unlist(read.csv('data/eseQ10.csv', header = F)))
# sdtsm.Q10 <- as.numeric(unlist(read.csv('eseQ10_start.csv', header = F)))
# sdtsm.Q10 <- as.numeric(unlist(read.csv('eseQ10_end.csv', header = F)))
dmtest_Q10 <- myDMtest(sdtsm.Q10, err10, alternative="less", h=10, power=2)

sdtsm.Q20 <- as.numeric(unlist(read.csv('eseQ20.csv', header = F)))
# sdtsm.Q20 <- as.numeric(unlist(read.csv('data/eseQ20.csv', header = F)))
# sdtsm.Q20 <- as.numeric(unlist(read.csv('eseQ20_start.csv', header = F)))
# sdtsm.Q20 <- as.numeric(unlist(read.csv('eseQ20_end.csv', header = F)))
dmtest_Q20 <- myDMtest(sdtsm.Q20, err20, alternative="less", h=20, power=2)

sdtsm.Q30 <- as.numeric(unlist(read.csv('eseQ30.csv', header = F)))
# sdtsm.Q30 <- as.numeric(unlist(read.csv('data/eseQ30.csv', header = F)))
# sdtsm.Q30 <- as.numeric(unlist(read.csv('eseQ30_start.csv', header = F)))
# sdtsm.Q30 <- as.numeric(unlist(read.csv('eseQ30_end.csv', header = F)))
dmtest_Q30 <- myDMtest(sdtsm.Q30, err30, alternative="less", h=30, power=2)

sdtsm.Q40 <- as.numeric(unlist(read.csv('eseQ40.csv', header = F)))
# sdtsm.Q40 <- as.numeric(unlist(read.csv('data/eseQ40.csv', header = F)))
# sdtsm.Q40 <- as.numeric(unlist(read.csv('eseQ40_start.csv', header = F)))
# sdtsm.Q40 <- as.numeric(unlist(read.csv('eseQ40_end.csv', header = F)))
dmtest_Q40 <- myDMtest(sdtsm.Q40, err40, alternative="less", h=40, power=2)

save.image("./results/myoos.RData")



###################################
# 40-quarter forecast of 1-year forward rate
##########################################

ep <- T-40

fhat <- matrix(NA,T,1)
f <- matrix(NA,T,1)
err <- matrix(NA,T,1)

for (i in 20:ep){
  sim <- simulateSSM(ssm(Y[1:i,], pars.mean), type="states") ## (N+1) x i
  sim.tau <- sim[,2,1]  # tau_exp = sim.tau
  sim.Ptilde <- sim[,3:5,1]
  Ptilde_exp <- sim.Ptilde
  for (j in 1:40){
    Ptilde_exp <- Ptilde_exp %*% t(pars.mean$Phi)
  }
  cp_exp <- rep(1, i) %*% rbind(pars.mean$Pbar) + sim.tau %*% rbind(pars.mean$gamma) + Ptilde_exp
  f[i+40] <- 2*Y[i+40,4] - Y[i+40,3]
  fhat[i+40] <- A1 - A2 +  cp_exp[i,] %*% (B1-B2)
  err[i+40] <- f[i+40] - fhat[i+40]
}

err40 <- err[60:(ep+40)]*100
rmse40 <- sqrt((t(err40)%*%err40)/length(err40))



###########################################################
############## 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


# 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


# 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


# 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


# 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







