##============================================================================##
# Jianhao Lin, Jiacheng Fan, Yifan Zhang, Liangyuan Chen, "Real-time 
#   Macroeconomic Projection Using Narrative Central Bank Communication", 
#   Journal of Applied Econometrics, forthcoming
##----------------------------------------------------------------------------##
  
# Forecast PPI in real time using PBC communication indices (PCIs) and 
# macro factors.

# Note: 
#   1) Estimate numerous individual models with different lags of predictors.
#   2) Select the best individual model with the smallest out-of-sample (OOS) RMSE.
#   3) Use the best individual model to obtain the combination model.
#   4) Use the combination model to generate the final OOS forecasts.

##----------------------------------------------------------------------------##
# Jiacheng Fan, 2022 (fanjch7@mail2.sysu.edu.cn)
##============================================================================##
library(readxl)
library(midasr)
library(openxlsx)
rm(list = ls())

# Adjust this path to where you stored the files
path = "~/lin-fan-zhang-chen-files"

# Set up paths for reading data and storing results
rawdatpath = paste0(path, "/Data/Raw_Data/")
tmpdatpath = paste0(path, "/Data/Temp_Data/")

setwd(path)


#=============================   Step (1) & (2)  ==============================#
#===============================================================================
# Nowcasting
#===============================================================================
#-------------------------------------
# Load Data
#-------------------------------------
# Target
tag <- na.omit(read.csv(paste0(rawdatpath, 'Target_PPI.csv'))[,2])

# PCIs: GDP
hdmr_gdp <- read.csv(paste0(tmpdatpath, 'PCIs_GDP_Realtime_Nowcast.csv'), header = FALSE)
zp1 <- hdmr_gdp[, seq(1, (length(hdmr_gdp)-3), 4)]
zo1 <- hdmr_gdp[, seq(2, (length(hdmr_gdp)-2), 4)]

# PCIs: CPI
hdmr_cpi <- read.csv(paste0(tmpdatpath, 'PCIs_CPI_Realtime_Nowcast.csv'), header = FALSE)
zp2 <- hdmr_cpi[, seq(1, (length(hdmr_cpi)-3), 4)]
zo2 <- hdmr_cpi[, seq(2, (length(hdmr_cpi)-2), 4)]

# PCIs: IVA
hdmr_iva <- read.csv(paste0(tmpdatpath, 'PCIs_IVA_Realtime_Nowcast.csv'), header = FALSE)
zp3 <- hdmr_iva[, seq(1, (length(hdmr_iva)-3), 4)]
zo3 <- hdmr_iva[, seq(2, (length(hdmr_iva)-2), 4)]

# PCIs: FAI
hdmr_fai <- read.csv(paste0(tmpdatpath, 'PCIs_FAI_Realtime_Nowcast.csv'), header = FALSE)
zp4 <- hdmr_fai[, seq(1, (length(hdmr_fai)-3), 4)]
zo4 <- hdmr_fai[, seq(2, (length(hdmr_fai)-2), 4)]

# PCIs: M2
hdmr_m2 <- read.csv(paste0(tmpdatpath, 'PCIs_M2_Realtime_Nowcast.csv'), header = FALSE)
zp5 <- hdmr_m2[, seq(1, (length(hdmr_m2)-3), 4)]
zo5 <- hdmr_m2[, seq(2, (length(hdmr_m2)-2), 4)]

# PCIs: PPI
hdmr_ppi <- read.csv(paste0(tmpdatpath, 'PCIs_PPI_Realtime_Nowcast.csv'), header = FALSE)
zp6 <- hdmr_ppi[, seq(1, (length(hdmr_ppi)-3), 4)]
zo6 <- hdmr_ppi[, seq(2, (length(hdmr_ppi)-2), 4)]

# Macro Factors
factor <- read.csv(paste0(tmpdatpath, 'MFactors_Realtime.csv'), header = F)
df1 <- data.frame(factor[, seq(1, (length(factor)-4) ,5)])
df2 <- data.frame(factor[, seq(2, (length(factor)-3), 5)])
df3 <- data.frame(factor[, seq(3, (length(factor)-2), 5)])
df4 <- data.frame(factor[, seq(4, (length(factor)-1), 5)])
df5 <- data.frame(factor[, seq(5, (length(factor)), 5)])


#-------------------------------------
# Set Estimating Windows & Variables
#-------------------------------------
len <- 7            # Initial window size = 7 years
leny <- 4 * len
lenx <- 12 * len

yy = ts(c(tag[1:leny], tag[leny+1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx1 = ts(c(zp1[1:lenx, 1], zp1[lenx+1, 1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx2 = ts(c(zo1[1:lenx, 1], zo1[lenx+1, 1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx3 = ts(c(zp2[1:lenx, 1], zp2[lenx+1, 1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx4 = ts(c(zo2[1:lenx, 1], zo2[lenx+1, 1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx5 = ts(c(zp3[1:lenx, 1], zp3[lenx+1, 1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx6 = ts(c(zo3[1:lenx, 1], zo3[lenx+1, 1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx7 = ts(c(zp4[1:lenx, 1], zp4[lenx+1, 1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx8 = ts(c(zo4[1:lenx, 1], zo4[lenx+1, 1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx9 = ts(c(zp5[1:lenx, 1], zp5[lenx+1, 1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx10 = ts(c(zo5[1:lenx, 1], zo5[lenx+1, 1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx11 = ts(c(zp6[1:lenx, 1], zp6[lenx+1, 1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx12 = ts(c(zo6[1:lenx, 1], zo6[lenx+1, 1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
zz1 = ts(c(df1[1:lenx, 1], df1[lenx+1, 2]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
zz2 = ts(c(df2[1:lenx, 1], df2[lenx+1, 2]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
zz3 = ts(c(df3[1:lenx, 1], df3[lenx+1, 2]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
zz4 = ts(c(df4[1:lenx, 1], df4[lenx+1, 2]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
zz5 = ts(c(df5[1:lenx, 1], df5[lenx+1, 2]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)


#-------------------------------------
# Estimate and Nowcast (h=0)
#-------------------------------------
Weightfunction <- c()
Ylag <- c()
Xlag <- c()

result_rec_now <- list(Ylag = c(),
                       Xlag = c(),
                       Recu_df1 = c(),
                       Recu_df2 = c(),
                       Recu_df3 = c(),
                       Recu_df4 = c(),
                       Recu_df5 = c(),
                       Recu_zpos_gdp = c(),
                       Recu_zzero_gdp = c(),
                       Recu_zpos_cpi = c(),
                       Recu_zzero_cpi = c(),
                       Recu_zpos_iva = c(),
                       Recu_zzero_iva = c(),
                       Recu_zpos_fai = c(),
                       Recu_zzero_fai = c(),
                       Recu_zpos_m2 = c(),
                       Recu_zzero_m2 = c(),
                       Recu_zpos_ppi = c(),
                       Recu_zzero_ppi = c())

Ylag.max <- 12
Xlag.max <- 12

for (ylag in seq(1, Ylag.max, 1)) {
  for (xlag in seq(1, Xlag.max, 1)) {
    result_rec_now$Ylag <- append(result_rec_now$Ylag, ylag)
    result_rec_now$Xlag <- append(result_rec_now$Xlag, xlag)
    
    # Setting
    lead.xx <- 0
    ##lagx.xx <- lead.xx + xlag
    lead.zz <- 1
    ##lagx.zz <- lead.zz + xlag
    
    # Training
    ## MF 
    mdz1 <- midas_r(yy ~ mls(yy, 1:ylag, 1) + mls(zz1, lead.zz:xlag, 1), start = NULL) # DF1
    mdz2 <- midas_r(yy ~ mls(yy, 1:ylag, 1) + mls(zz2, lead.zz:xlag, 1), start = NULL) # DF2
    mdz3 <- midas_r(yy ~ mls(yy, 1:ylag, 1) + mls(zz3, lead.zz:xlag, 1), start = NULL) # DF3
    mdz4 <- midas_r(yy ~ mls(yy, 1:ylag, 1) + mls(zz4, lead.zz:xlag, 1), start = NULL) # DF4
    mdz5 <- midas_r(yy ~ mls(yy, 1:ylag, 1) + mls(zz5, lead.zz:xlag, 1), start = NULL) # DF5
    
    ## HDMR: GDP
    mdx1 <- midas_r(yy ~ mls(yy, 1:ylag, 1) + mls(xx1, lead.xx:xlag, 1), start = NULL) # zp1
    mdx2 <- midas_r(yy ~ mls(yy, 1:ylag, 1) + mls(xx2, lead.xx:xlag, 1), start = NULL) # zo1
    ## HDMR: CPI
    mdx3 <- midas_r(yy ~ mls(yy, 1:ylag, 1) + mls(xx3, lead.xx:xlag, 1), start = NULL) # zp2
    mdx4 <- midas_r(yy ~ mls(yy, 1:ylag, 1) + mls(xx4, lead.xx:xlag, 1), start = NULL) # zo2
    ## HDMR: IVA
    mdx5 <- midas_r(yy ~ mls(yy, 1:ylag, 1) + mls(xx5, lead.xx:xlag, 1), start = NULL) # zp3
    mdx6 <- midas_r(yy ~ mls(yy, 1:ylag, 1) + mls(xx6, lead.xx:xlag, 1), start = NULL) # zo3
    ## HDMR: FAI
    mdx7 <- midas_r(yy ~ mls(yy, 1:ylag, 1) + mls(xx7, lead.xx:xlag, 1), start = NULL) # zp4
    mdx8 <- midas_r(yy ~ mls(yy, 1:ylag, 1) + mls(xx8, lead.xx:xlag, 1), start = NULL) # zo4
    ## HDMR: M2
    mdx9 <- midas_r(yy ~ mls(yy, 1:ylag, 1) + mls(xx9, lead.xx:xlag, 1), start = NULL) # zp5
    mdx10 <- midas_r(yy ~ mls(yy, 1:ylag, 1) + mls(xx10, lead.xx:xlag, 1), start = NULL) # zo5
    ## HDMR: PPI
    mdx11 <- midas_r(yy ~ mls(yy, 1:ylag, 1) + mls(xx11, lead.xx:xlag, 1), start = NULL) # zp6
    mdx12 <- midas_r(yy ~ mls(yy, 1:ylag, 1) + mls(xx12, lead.xx:xlag, 1), start = NULL) # zo6
    
    # Predicting
    Out_recu_df_1 <- c()
    Out_recu_df_2 <- c()
    Out_recu_df_3 <- c()
    Out_recu_df_4 <- c()
    Out_recu_df_5 <- c()
    Out_recu_zp1 <- c()
    Out_recu_zo1 <- c()
    Out_recu_zp2 <- c()
    Out_recu_zo2 <- c()
    Out_recu_zp3 <- c()
    Out_recu_zo3 <- c()
    Out_recu_zp4 <- c()
    Out_recu_zo4 <- c()
    Out_recu_zp5 <- c()
    Out_recu_zo5 <- c()
    Out_recu_zp6 <- c()
    Out_recu_zo6 <- c()
    
    for (i in seq(1, length(zp2), 1)) {
      # Sampling
      i1 = i + 1
      m1 = i + leny - 1
      m2 = m1 + 1
      
      # Start & End
      start = c(2003,1)
      end = start + c(len, i-1)
      fulldata <- list(
        yy = ts(c(tag[1:m1], tag[m2]),start, end, 12),
        xx1 = ts(c(zp1[1:m1, i], zp1[m2, i]), start, end, 12),
        xx2 = ts(c(zo1[1:m1, i], zo1[m2, i]), start, end, 12),
        xx3 = ts(c(zp2[1:m1, i], zp2[m2, i]), start, end, 12),
        xx4 = ts(c(zo2[1:m1, i], zo2[m2, i]), start, end, 12),
        xx5 = ts(c(zp3[1:m1, i], zp3[m2, i]), start, end, 12),
        xx6 = ts(c(zo3[1:m1, i], zo3[m2, i]), start, end, 12),
        xx7 = ts(c(zp4[1:m1, i], zp4[m2, i]), start, end, 12),
        xx8 = ts(c(zo4[1:m1, i], zo4[m2, i]), start, end, 12),
        xx9 = ts(c(zp5[1:m1, i], zp5[m2, i]), start, end, 12),
        xx10 = ts(c(zo5[1:m1, i], zo5[m2, i]), start, end, 12),
        xx11 = ts(c(zp6[1:m1, i], zp6[m2, i]), start, end, 12),
        xx12 = ts(c(zo6[1:m1, i], zo6[m2, i]), start, end, 12),
        zz1 = ts(c(df1[1:m1, i], df1[m1, i]), start, end, 12),
        zz2 = ts(c(df2[1:m1, i], df2[m1, i]), start, end, 12),
        zz3 = ts(c(df3[1:m1, i], df3[m1, i]), start, end, 12),
        zz4 = ts(c(df4[1:m1, i], df4[m1, i]), start, end, 12),
        zz5 = ts(c(df5[1:m1, i], df5[m1, i]), start, end, 12)
      )
      insample <- 1:m1
      outsample <- (1:length(fulldata$yy))[-insample]
      
      ## Combine MF & CB Recuing window
      recu_combmdxz <- try({
        average_forecast(list(mdz1, mdz2, mdz3, mdz4, mdz5, 
                              mdx1, mdx2, mdx3, mdx4, mdx5, mdx6,
                              mdx7, mdx8, mdx9, mdx10, mdx11, mdx12), 
                         data = fulldata, insample = insample, outsample = outsample, 
                         type = "recursive",
                         measures = c("MSE"),
                         fweights = c("EW", "BICW", "MSFE", "DMSFE"), 
                         show_progress = FALSE)}, silent = T) 
      if ('try-error' %in% class(recu_combmdxz)) {
        recu_forexz <- rep(NA, 17)
      } else {
        recu_forexz <- recu_combmdxz[["forecast"]]
      }
      
      # Collecting
      Out_recu_df_1 <- append(Out_recu_df_1, recu_forexz[1])
      Out_recu_df_2 <- append(Out_recu_df_2, recu_forexz[2])
      Out_recu_df_3 <- append(Out_recu_df_3, recu_forexz[3])
      Out_recu_df_4 <- append(Out_recu_df_4, recu_forexz[4])
      Out_recu_df_5 <- append(Out_recu_df_5, recu_forexz[5])
      Out_recu_zp1 <- append(Out_recu_zp1, recu_forexz[6])
      Out_recu_zo1 <- append(Out_recu_zo1, recu_forexz[7])
      Out_recu_zp2 <- append(Out_recu_zp2, recu_forexz[8])
      Out_recu_zo2 <- append(Out_recu_zo2, recu_forexz[9])
      Out_recu_zp3 <- append(Out_recu_zp3, recu_forexz[10])
      Out_recu_zo3 <- append(Out_recu_zo3, recu_forexz[11])
      Out_recu_zp4 <- append(Out_recu_zp4, recu_forexz[12])
      Out_recu_zo4 <- append(Out_recu_zo4, recu_forexz[13])
      Out_recu_zp5 <- append(Out_recu_zp5, recu_forexz[14])
      Out_recu_zo5 <- append(Out_recu_zo5, recu_forexz[15])
      Out_recu_zp6 <- append(Out_recu_zp6, recu_forexz[16])
      Out_recu_zo6 <- append(Out_recu_zo6, recu_forexz[17])
    }
    
    Rmse_df1 <- sqrt(sum((Out_recu_df_1 - tag[-c(1:leny)])^2) / length(Out_recu_df_1))
    Rmse_df2 <- sqrt(sum((Out_recu_df_2 - tag[-c(1:leny)])^2) / length(Out_recu_df_2))
    Rmse_df3 <- sqrt(sum((Out_recu_df_3 - tag[-c(1:leny)])^2) / length(Out_recu_df_3))
    Rmse_df4 <- sqrt(sum((Out_recu_df_4 - tag[-c(1:leny)])^2) / length(Out_recu_df_4))
    Rmse_df5 <- sqrt(sum((Out_recu_df_5 - tag[-c(1:leny)])^2) / length(Out_recu_df_5))
    Rmse_zpos_gdp <- sqrt(sum((Out_recu_zp1 - tag[-c(1:leny)])^2) / length(Out_recu_zp1))
    Rmse_zzero_gdp <- sqrt(sum((Out_recu_zo1 - tag[-c(1:leny)])^2) / length(Out_recu_zo1))
    Rmse_zpos_cpi <- sqrt(sum((Out_recu_zp2 - tag[-c(1:leny)])^2) / length(Out_recu_zp2))
    Rmse_zzero_cpi <- sqrt(sum((Out_recu_zo2 - tag[-c(1:leny)])^2) / length(Out_recu_zo2))
    Rmse_zpos_iva <- sqrt(sum((Out_recu_zp3 - tag[-c(1:leny)])^2) / length(Out_recu_zp3))
    Rmse_zzero_iva <- sqrt(sum((Out_recu_zo3 - tag[-c(1:leny)])^2) / length(Out_recu_zo3))
    Rmse_zpos_fai <- sqrt(sum((Out_recu_zp4 - tag[-c(1:leny)])^2) / length(Out_recu_zp4))
    Rmse_zzero_fai <- sqrt(sum((Out_recu_zo4 - tag[-c(1:leny)])^2) / length(Out_recu_zo4))
    Rmse_zpos_m2 <- sqrt(sum((Out_recu_zp5 - tag[-c(1:leny)])^2) / length(Out_recu_zp5))
    Rmse_zzero_m2 <- sqrt(sum((Out_recu_zo5 - tag[-c(1:leny)])^2) / length(Out_recu_zo5))
    Rmse_zpos_ppi <- sqrt(sum((Out_recu_zp6 - tag[-c(1:leny)])^2) / length(Out_recu_zp6))
    Rmse_zzero_ppi <- sqrt(sum((Out_recu_zo6 - tag[-c(1:leny)])^2) / length(Out_recu_zo6))
    
    result_rec_now$Recu_df1 <- append(result_rec_now$Recu_df1, Rmse_df1)
    result_rec_now$Recu_df2 <- append(result_rec_now$Recu_df2, Rmse_df2)
    result_rec_now$Recu_df3 <- append(result_rec_now$Recu_df3, Rmse_df3)
    result_rec_now$Recu_df4 <- append(result_rec_now$Recu_df4, Rmse_df4)
    result_rec_now$Recu_df5 <- append(result_rec_now$Recu_df5, Rmse_df5)
    result_rec_now$Recu_zpos_gdp <- append(result_rec_now$Recu_zpos_gdp, Rmse_zpos_gdp)
    result_rec_now$Recu_zzero_gdp <- append(result_rec_now$Recu_zzero_gdp, Rmse_zzero_gdp)
    result_rec_now$Recu_zpos_cpi <- append(result_rec_now$Recu_zpos_cpi, Rmse_zpos_cpi)
    result_rec_now$Recu_zzero_cpi <- append(result_rec_now$Recu_zzero_cpi, Rmse_zzero_cpi)
    result_rec_now$Recu_zpos_iva <- append(result_rec_now$Recu_zpos_iva, Rmse_zpos_iva)
    result_rec_now$Recu_zzero_iva <- append(result_rec_now$Recu_zzero_iva, Rmse_zzero_iva)
    result_rec_now$Recu_zpos_fai <- append(result_rec_now$Recu_zpos_fai, Rmse_zpos_fai)
    result_rec_now$Recu_zzero_fai <- append(result_rec_now$Recu_zzero_fai, Rmse_zzero_fai)
    result_rec_now$Recu_zpos_m2 <- append(result_rec_now$Recu_zpos_m2, Rmse_zpos_m2)
    result_rec_now$Recu_zzero_m2 <- append(result_rec_now$Recu_zzero_m2, Rmse_zzero_m2)
    result_rec_now$Recu_zpos_ppi <- append(result_rec_now$Recu_zpos_ppi, Rmse_zpos_ppi)
    result_rec_now$Recu_zzero_ppi <- append(result_rec_now$Recu_zzero_ppi, Rmse_zzero_ppi)
  }
}


#===============================================================================
# Forecasting
#===============================================================================
#-------------------------------------
# Load Data
#-------------------------------------
# Target
tag <- na.omit(read.csv(paste0(rawdatpath, 'Target_PPI.csv'))[,2])

# PCIs: GDP
hdmr_gdp <- read.csv(paste0(tmpdatpath, 'PCIs_GDP_Realtime.csv'), header = FALSE)
zp1 <- hdmr_gdp[, seq(1, (length(hdmr_gdp)-3), 4)]
zo1 <- hdmr_gdp[, seq(2, (length(hdmr_gdp)-2), 4)]

# PCIs: CPI
hdmr_cpi <- read.csv(paste0(tmpdatpath, 'PCIs_CPI_Realtime.csv'), header = FALSE)
zp2 <- hdmr_cpi[, seq(1, (length(hdmr_cpi)-3), 4)]
zo2 <- hdmr_cpi[, seq(2, (length(hdmr_cpi)-2), 4)]

# PCIs: IVA
hdmr_iva <- read.csv(paste0(tmpdatpath, 'PCIs_IVA_Realtime.csv'), header = FALSE)
zp3 <- hdmr_iva[, seq(1, (length(hdmr_iva)-3), 4)]
zo3 <- hdmr_iva[, seq(2, (length(hdmr_iva)-2), 4)]

# PCIs: FAI
hdmr_fai <- read.csv(paste0(tmpdatpath, 'PCIs_FAI_Realtime.csv'), header = FALSE)
zp4 <- hdmr_fai[, seq(1, (length(hdmr_fai)-3), 4)]
zo4 <- hdmr_fai[, seq(2, (length(hdmr_fai)-2), 4)]

# PCIs: M2
hdmr_m2 <- read.csv(paste0(tmpdatpath, 'PCIs_M2_Realtime.csv'), header = FALSE)
zp5 <- hdmr_m2[, seq(1, (length(hdmr_m2)-3), 4)]
zo5 <- hdmr_m2[, seq(2, (length(hdmr_m2)-2), 4)]

# PCIs: PPI
hdmr_ppi <- read.csv(paste0(tmpdatpath, 'PCIs_PPI_Realtime.csv'), header = FALSE)
zp6 <- hdmr_ppi[, seq(1, (length(hdmr_ppi)-3), 4)]
zo6 <- hdmr_ppi[, seq(2, (length(hdmr_ppi)-2), 4)]

# Macro Factors
factor <- read.csv(paste0(tmpdatpath, 'MFactors_Realtime.csv'), header = F)
df1 <- data.frame(factor[, seq(1, (length(factor)-4) ,5)])
df2 <- data.frame(factor[, seq(2, (length(factor)-3), 5)])
df3 <- data.frame(factor[, seq(3, (length(factor)-2), 5)])
df4 <- data.frame(factor[, seq(4, (length(factor)-1), 5)])
df5 <- data.frame(factor[, seq(5, (length(factor)), 5)])


#-------------------------------------
# Set Estimating Windows & Variables
#-------------------------------------
len <- 7            # Initial window size = 7 years
leny <- 4 * len
lenx <- 12 * len

yy = ts(c(tag[1:leny], tag[leny+1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx1 = ts(c(zp1[1:lenx, 1], zp1[(lenx+1), 2]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx2 = ts(c(zo1[1:lenx, 1], zo1[(lenx+1), 2]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx3 = ts(c(zp2[1:lenx, 1], zp2[(lenx+1), 2]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx4 = ts(c(zo2[1:lenx, 1], zo2[(lenx+1), 2]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx5 = ts(c(zp3[1:lenx, 1], zp3[(lenx+1), 2]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx6 = ts(c(zo3[1:lenx, 1], zo3[(lenx+1), 2]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx7 = ts(c(zp4[1:lenx, 1], zp4[(lenx+1), 2]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx8 = ts(c(zo4[1:lenx, 1], zo4[(lenx+1), 2]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx9 = ts(c(zp5[1:lenx, 1], zp5[(lenx+1), 2]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx10 = ts(c(zo5[1:lenx, 1], zo5[(lenx+1), 2]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx11 = ts(c(zp6[1:lenx, 1], zp6[(lenx+1), 2]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx12 = ts(c(zo6[1:lenx, 1], zo6[(lenx+1), 2]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
zz1 = ts(c(df1[1:lenx, 1], df1[lenx, 1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
zz2 = ts(c(df2[1:lenx, 1], df2[lenx, 1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
zz3 = ts(c(df3[1:lenx, 1], df3[lenx, 1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
zz4 = ts(c(df4[1:lenx, 1], df4[lenx, 1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
zz5 = ts(c(df5[1:lenx, 1], df5[lenx, 1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)


#-------------------------------------
# Estimate and Forecast (h=1)
#-------------------------------------
Weightfunction <- c()
Ylag <- c()
Xlag <- c()

result_rec_for1 <- list(Ylag = c(),
                        Xlag = c(),
                        Tau = c(),
                        Recu_df1 = c(),
                        Recu_df2 = c(),
                        Recu_df3 = c(),
                        Recu_df4 = c(),
                        Recu_df5 = c(),
                        Recu_zpos_gdp = c(),
                        Recu_zzero_gdp = c(),
                        Recu_zpos_cpi = c(),
                        Recu_zzero_cpi = c(),
                        Recu_zpos_iva = c(),
                        Recu_zzero_iva = c(),
                        Recu_zpos_fai = c(),
                        Recu_zzero_fai = c(),
                        Recu_zpos_m2 = c(),
                        Recu_zzero_m2 = c(),
                        Recu_zpos_ppi = c(),
                        Recu_zzero_ppi = c())

Ylag.max <- 12
Xlag.max <- 12

for (ylag in seq(1, Ylag.max, 1)) {
  for (xlag in seq(1, Xlag.max, 1)) {
    result_rec_for1$Ylag <- append(result_rec_for1$Ylag, ylag)
    result_rec_for1$Xlag <- append(result_rec_for1$Xlag, xlag)
    result_rec_for1$Tau <- append(result_rec_for1$Tau, 1)
    
    # Training
    ## MF 
    mdz1 <- midas_r(yy ~ mls(yy, 1:ylag, 1) + mls(zz1, 1:xlag, 1), start = NULL) # DF1
    mdz2 <- midas_r(yy ~ mls(yy, 1:ylag, 1) + mls(zz2, 1:xlag, 1), start = NULL) # DF2
    mdz3 <- midas_r(yy ~ mls(yy, 1:ylag, 1) + mls(zz3, 1:xlag, 1), start = NULL) # DF3
    mdz4 <- midas_r(yy ~ mls(yy, 1:ylag, 1) + mls(zz4, 1:xlag, 1), start = NULL) # DF4
    mdz5 <- midas_r(yy ~ mls(yy, 1:ylag, 1) + mls(zz5, 1:xlag, 1), start = NULL) # DF5
    
    ## HDMR: GDP
    mdx1 <- midas_r(yy ~ mls(yy, 1:ylag, 1) + mls(xx1, 1:xlag, 1), start = NULL) # zp1
    mdx2 <- midas_r(yy ~ mls(yy, 1:ylag, 1) + mls(xx2, 1:xlag, 1), start = NULL) # zo1
    ## HDMR: CPI
    mdx3 <- midas_r(yy ~ mls(yy, 1:ylag, 1) + mls(xx3, 1:xlag, 1), start = NULL) # zp2
    mdx4 <- midas_r(yy ~ mls(yy, 1:ylag, 1) + mls(xx4, 1:xlag, 1), start = NULL) # zo2
    ## HDMR: IVA
    mdx5 <- midas_r(yy ~ mls(yy, 1:ylag, 1) + mls(xx5, 1:xlag, 1), start = NULL) # zp3
    mdx6 <- midas_r(yy ~ mls(yy, 1:ylag, 1) + mls(xx6, 1:xlag, 1), start = NULL) # zo3
    ## HDMR: FAI
    mdx7 <- midas_r(yy ~ mls(yy, 1:ylag, 1) + mls(xx7, 1:xlag, 1), start = NULL) # zp4
    mdx8 <- midas_r(yy ~ mls(yy, 1:ylag, 1) + mls(xx8, 1:xlag, 1), start = NULL) # zo4
    ## HDMR: M2
    mdx9 <- midas_r(yy ~ mls(yy, 1:ylag, 1) + mls(xx9, 1:xlag, 1), start = NULL) # zp5
    mdx10 <- midas_r(yy ~ mls(yy, 1:ylag, 1) + mls(xx10, 1:xlag, 1), start = NULL) # zo5
    ## HDMR: PPI
    mdx11 <- midas_r(yy ~ mls(yy, 1:ylag, 1) + mls(xx11, 1:xlag, 1), start = NULL) # zp6
    mdx12 <- midas_r(yy ~ mls(yy, 1:ylag, 1) + mls(xx12, 1:xlag, 1), start = NULL) # zo6
    
    # Predicting
    Out_recu_df_1 <- c()
    Out_recu_df_2 <- c()
    Out_recu_df_3 <- c()
    Out_recu_df_4 <- c()
    Out_recu_df_5 <- c()
    Out_recu_zp1 <- c()
    Out_recu_zo1 <- c()
    Out_recu_zp2 <- c()
    Out_recu_zo2 <- c()
    Out_recu_zp3 <- c()
    Out_recu_zo3 <- c()
    Out_recu_zp4 <- c()
    Out_recu_zo4 <- c()
    Out_recu_zp5 <- c()
    Out_recu_zo5 <- c()
    Out_recu_zp6 <- c()
    Out_recu_zo6 <- c()
    
    for (i in seq(1, (length(zp2)-1), 1)) {
      # Sampling
      i1 = i + 1
      m1 = i + leny - 1
      m2 = m1 + 1
      
      # Start & End
      start = c(2003,1)
      end = start + c(len, i-1)
      fulldata <- list(
        yy = ts(c(tag[1:m1], tag[m2]),start, end, 12),
        xx1 = ts(c(zp1[1:m1, i], zp1[m2, i1]), start, end, 12),
        xx2 = ts(c(zo1[1:m1, i], zo1[m2, i1]), start, end, 12),
        xx3 = ts(c(zp2[1:m1, i], zp2[m2, i1]), start, end, 12),
        xx4 = ts(c(zo2[1:m1, i], zo2[m2, i1]), start, end, 12),
        xx5 = ts(c(zp3[1:m1, i], zp3[m2, i1]), start, end, 12),
        xx6 = ts(c(zo3[1:m1, i], zo3[m2, i1]), start, end, 12),
        xx7 = ts(c(zp4[1:m1, i], zp4[m2, i1]), start, end, 12),
        xx8 = ts(c(zo4[1:m1, i], zo4[m2, i1]), start, end, 12),
        xx9 = ts(c(zp5[1:m1, i], zp5[m2, i1]), start, end, 12),
        xx10 = ts(c(zo5[1:m1, i], zo5[m2, i1]), start, end, 12),
        xx11 = ts(c(zp6[1:m1, i], zp6[m2, i1]), start, end, 12),
        xx12 = ts(c(zo6[1:m1, i], zo6[m2, i1]), start, end, 12),
        zz1 = ts(c(df1[1:m1, i], df1[m1, i]), start, end, 12),
        zz2 = ts(c(df2[1:m1, i], df2[m1, i]), start, end, 12),
        zz3 = ts(c(df3[1:m1, i], df3[m1, i]), start, end, 12),
        zz4 = ts(c(df4[1:m1, i], df4[m1, i]), start, end, 12),
        zz5 = ts(c(df5[1:m1, i], df5[m1, i]), start, end, 12)
      )
      insample <- 1:m1
      outsample <- (1:length(fulldata$yy))[-insample]
      
      ## Combine MF & CB Recuing window
      recu_combmdxz <- try({
        average_forecast(list(mdz1, mdz2, mdz3, mdz4, mdz5, 
                              mdx1, mdx2, mdx3, mdx4, mdx5, mdx6,
                              mdx7, mdx8, mdx9, mdx10, mdx11, mdx12), 
                         data = fulldata, insample = insample, outsample = outsample, 
                         type = "recursive",
                         measures = c("MSE"),
                         fweights = c("EW", "BICW", "MSFE", "DMSFE"), 
                         show_progress = FALSE)}, silent = T) 
      if ('try-error' %in% class(recu_combmdxz)) {
        recu_forexz <- rep(NA, 17)
      } else {
        recu_forexz <- recu_combmdxz[["forecast"]]
      }
      
      # Collecting
      Out_recu_df_1 <- append(Out_recu_df_1, recu_forexz[1])
      Out_recu_df_2 <- append(Out_recu_df_2, recu_forexz[2])
      Out_recu_df_3 <- append(Out_recu_df_3, recu_forexz[3])
      Out_recu_df_4 <- append(Out_recu_df_4, recu_forexz[4])
      Out_recu_df_5 <- append(Out_recu_df_5, recu_forexz[5])
      Out_recu_zp1 <- append(Out_recu_zp1, recu_forexz[6])
      Out_recu_zo1 <- append(Out_recu_zo1, recu_forexz[7])
      Out_recu_zp2 <- append(Out_recu_zp2, recu_forexz[8])
      Out_recu_zo2 <- append(Out_recu_zo2, recu_forexz[9])
      Out_recu_zp3 <- append(Out_recu_zp3, recu_forexz[10])
      Out_recu_zo3 <- append(Out_recu_zo3, recu_forexz[11])
      Out_recu_zp4 <- append(Out_recu_zp4, recu_forexz[12])
      Out_recu_zo4 <- append(Out_recu_zo4, recu_forexz[13])
      Out_recu_zp5 <- append(Out_recu_zp5, recu_forexz[14])
      Out_recu_zo5 <- append(Out_recu_zo5, recu_forexz[15])
      Out_recu_zp6 <- append(Out_recu_zp6, recu_forexz[16])
      Out_recu_zo6 <- append(Out_recu_zo6, recu_forexz[17])
    }
    
    Rmse_df1 <- sqrt(sum((Out_recu_df_1 - tag[-c(1:leny)])^2) / length(Out_recu_df_1))
    Rmse_df2 <- sqrt(sum((Out_recu_df_2 - tag[-c(1:leny)])^2) / length(Out_recu_df_2))
    Rmse_df3 <- sqrt(sum((Out_recu_df_3 - tag[-c(1:leny)])^2) / length(Out_recu_df_3))
    Rmse_df4 <- sqrt(sum((Out_recu_df_4 - tag[-c(1:leny)])^2) / length(Out_recu_df_4))
    Rmse_df5 <- sqrt(sum((Out_recu_df_5 - tag[-c(1:leny)])^2) / length(Out_recu_df_5))
    Rmse_zpos_gdp <- sqrt(sum((Out_recu_zp1 - tag[-c(1:leny)])^2) / length(Out_recu_zp1))
    Rmse_zzero_gdp <- sqrt(sum((Out_recu_zo1 - tag[-c(1:leny)])^2) / length(Out_recu_zo1))
    Rmse_zpos_cpi <- sqrt(sum((Out_recu_zp2 - tag[-c(1:leny)])^2) / length(Out_recu_zp2))
    Rmse_zzero_cpi <- sqrt(sum((Out_recu_zo2 - tag[-c(1:leny)])^2) / length(Out_recu_zo2))
    Rmse_zpos_iva <- sqrt(sum((Out_recu_zp3 - tag[-c(1:leny)])^2) / length(Out_recu_zp3))
    Rmse_zzero_iva <- sqrt(sum((Out_recu_zo3 - tag[-c(1:leny)])^2) / length(Out_recu_zo3))
    Rmse_zpos_fai <- sqrt(sum((Out_recu_zp4 - tag[-c(1:leny)])^2) / length(Out_recu_zp4))
    Rmse_zzero_fai <- sqrt(sum((Out_recu_zo4 - tag[-c(1:leny)])^2) / length(Out_recu_zo4))
    Rmse_zpos_m2 <- sqrt(sum((Out_recu_zp5 - tag[-c(1:leny)])^2) / length(Out_recu_zp5))
    Rmse_zzero_m2 <- sqrt(sum((Out_recu_zo5 - tag[-c(1:leny)])^2) / length(Out_recu_zo5))
    Rmse_zpos_ppi <- sqrt(sum((Out_recu_zp6 - tag[-c(1:leny)])^2) / length(Out_recu_zp6))
    Rmse_zzero_ppi <- sqrt(sum((Out_recu_zo6 - tag[-c(1:leny)])^2) / length(Out_recu_zo6))
    
    result_rec_for1$Recu_df1 <- append(result_rec_for1$Recu_df1, Rmse_df1)
    result_rec_for1$Recu_df2 <- append(result_rec_for1$Recu_df2, Rmse_df2)
    result_rec_for1$Recu_df3 <- append(result_rec_for1$Recu_df3, Rmse_df3)
    result_rec_for1$Recu_df4 <- append(result_rec_for1$Recu_df4, Rmse_df4)
    result_rec_for1$Recu_df5 <- append(result_rec_for1$Recu_df5, Rmse_df5)
    result_rec_for1$Recu_zpos_gdp <- append(result_rec_for1$Recu_zpos_gdp, Rmse_zpos_gdp)
    result_rec_for1$Recu_zzero_gdp <- append(result_rec_for1$Recu_zzero_gdp, Rmse_zzero_gdp)
    result_rec_for1$Recu_zpos_cpi <- append(result_rec_for1$Recu_zpos_cpi, Rmse_zpos_cpi)
    result_rec_for1$Recu_zzero_cpi <- append(result_rec_for1$Recu_zzero_cpi, Rmse_zzero_cpi)
    result_rec_for1$Recu_zpos_iva <- append(result_rec_for1$Recu_zpos_iva, Rmse_zpos_iva)
    result_rec_for1$Recu_zzero_iva <- append(result_rec_for1$Recu_zzero_iva, Rmse_zzero_iva)
    result_rec_for1$Recu_zpos_fai <- append(result_rec_for1$Recu_zpos_fai, Rmse_zpos_fai)
    result_rec_for1$Recu_zzero_fai <- append(result_rec_for1$Recu_zzero_fai, Rmse_zzero_fai)
    result_rec_for1$Recu_zpos_m2 <- append(result_rec_for1$Recu_zpos_m2, Rmse_zpos_m2)
    result_rec_for1$Recu_zzero_m2 <- append(result_rec_for1$Recu_zzero_m2, Rmse_zzero_m2)
    result_rec_for1$Recu_zpos_ppi <- append(result_rec_for1$Recu_zpos_ppi, Rmse_zpos_ppi)
    result_rec_for1$Recu_zzero_ppi <- append(result_rec_for1$Recu_zzero_ppi, Rmse_zzero_ppi)
  }
}


#-------------------------------------
# Estimate and Forecast (h=3)
#-------------------------------------
Weightfunction <- c()
Ylag <- c()
Xlag <- c()

result_rec_for3 <- list(Ylag = c(),
                        Xlag = c(),
                        Tau = c(),
                        Recu_df1 = c(),
                        Recu_df2 = c(),
                        Recu_df3 = c(),
                        Recu_df4 = c(),
                        Recu_df5 = c(),
                        Recu_zpos_gdp = c(),
                        Recu_zzero_gdp = c(),
                        Recu_zpos_cpi = c(),
                        Recu_zzero_cpi = c(),
                        Recu_zpos_iva = c(),
                        Recu_zzero_iva = c(),
                        Recu_zpos_fai = c(),
                        Recu_zzero_fai = c(),
                        Recu_zpos_m2 = c(),
                        Recu_zzero_m2 = c(),
                        Recu_zpos_ppi = c(),
                        Recu_zzero_ppi = c())

Ylag.max <- 12
Xlag.max <- 12

for (ylag in seq(1, Ylag.max, 1)) {
  for (xlag in seq(1, Xlag.max, 1)) {
    result_rec_for3$Ylag <- append(result_rec_for3$Ylag, ylag)
    result_rec_for3$Xlag <- append(result_rec_for3$Xlag, xlag)
    result_rec_for3$Tau <- append(result_rec_for3$Tau, 3)
    
    # Setting
    lagy <- 2 + ylag
    lagx <- 2 + xlag
    
    # Training
    ## MF 
    mdz1 <- midas_r(yy ~ mls(yy, 3:lagy, 1) + mls(zz1, 3:lagx, 1), start = NULL) # DF1
    mdz2 <- midas_r(yy ~ mls(yy, 3:lagy, 1) + mls(zz2, 3:lagx, 1), start = NULL) # DF2
    mdz3 <- midas_r(yy ~ mls(yy, 3:lagy, 1) + mls(zz3, 3:lagx, 1), start = NULL) # DF3
    mdz4 <- midas_r(yy ~ mls(yy, 3:lagy, 1) + mls(zz4, 3:lagx, 1), start = NULL) # DF4
    mdz5 <- midas_r(yy ~ mls(yy, 3:lagy, 1) + mls(zz5, 3:lagx, 1), start = NULL) # DF5
    
    ## HDMR: GDP
    mdx1 <- midas_r(yy ~ mls(yy, 3:lagy, 1) + mls(xx1, 3:lagx, 1), start = NULL) # zp1
    mdx2 <- midas_r(yy ~ mls(yy, 3:lagy, 1) + mls(xx2, 3:lagx, 1), start = NULL) # zo1
    ## HDMR: CPI
    mdx3 <- midas_r(yy ~ mls(yy, 3:lagy, 1) + mls(xx3, 3:lagx, 1), start = NULL) # zp2
    mdx4 <- midas_r(yy ~ mls(yy, 3:lagy, 1) + mls(xx4, 3:lagx, 1), start = NULL) # zo2
    ## HDMR: IVA
    mdx5 <- midas_r(yy ~ mls(yy, 3:lagy, 1) + mls(xx5, 3:lagx, 1), start = NULL) # zp3
    mdx6 <- midas_r(yy ~ mls(yy, 3:lagy, 1) + mls(xx6, 3:lagx, 1), start = NULL) # zo3
    ## HDMR: FAI
    mdx7 <- midas_r(yy ~ mls(yy, 3:lagy, 1) + mls(xx7, 3:lagx, 1), start = NULL) # zp4
    mdx8 <- midas_r(yy ~ mls(yy, 3:lagy, 1) + mls(xx8, 3:lagx, 1), start = NULL) # zo4
    ## HDMR: M2
    mdx9 <- midas_r(yy ~ mls(yy, 3:lagy, 1) + mls(xx9, 3:lagx, 1), start = NULL) # zp5
    mdx10 <- midas_r(yy ~ mls(yy, 3:lagy, 1) + mls(xx10, 3:lagx, 1), start = NULL) # zo5
    ## HDMR: PPI
    mdx11 <- midas_r(yy ~ mls(yy, 3:lagy, 1) + mls(xx11, 3:lagx, 1), start = NULL) # zp6
    mdx12 <- midas_r(yy ~ mls(yy, 3:lagy, 1) + mls(xx12, 3:lagx, 1), start = NULL) # zo6
    
    # Predicting
    Out_recu_df_1 <- c()
    Out_recu_df_2 <- c()
    Out_recu_df_3 <- c()
    Out_recu_df_4 <- c()
    Out_recu_df_5 <- c()
    Out_recu_zp1 <- c()
    Out_recu_zo1 <- c()
    Out_recu_zp2 <- c()
    Out_recu_zo2 <- c()
    Out_recu_zp3 <- c()
    Out_recu_zo3 <- c()
    Out_recu_zp4 <- c()
    Out_recu_zo4 <- c()
    Out_recu_zp5 <- c()
    Out_recu_zo5 <- c()
    Out_recu_zp6 <- c()
    Out_recu_zo6 <- c()
    
    for (i in seq(1, (length(zp2)-1), 1)) {
      # Sampling
      i1 = i + 1
      m1 = i + leny - 1
      m2 = m1 + 1
      
      # Start & End
      start = c(2003,1)
      end = start + c(len, i-1)
      fulldata <- list(
        yy = ts(c(tag[1:m1], tag[m2]),start, end, 12),
        xx1 = ts(c(zp1[1:m1, i], zp1[m2, i1]), start, end, 12),
        xx2 = ts(c(zo1[1:m1, i], zo1[m2, i1]), start, end, 12),
        xx3 = ts(c(zp2[1:m1, i], zp2[m2, i1]), start, end, 12),
        xx4 = ts(c(zo2[1:m1, i], zo2[m2, i1]), start, end, 12),
        xx5 = ts(c(zp3[1:m1, i], zp3[m2, i1]), start, end, 12),
        xx6 = ts(c(zo3[1:m1, i], zo3[m2, i1]), start, end, 12),
        xx7 = ts(c(zp4[1:m1, i], zp4[m2, i1]), start, end, 12),
        xx8 = ts(c(zo4[1:m1, i], zo4[m2, i1]), start, end, 12),
        xx9 = ts(c(zp5[1:m1, i], zp5[m2, i1]), start, end, 12),
        xx10 = ts(c(zo5[1:m1, i], zo5[m2, i1]), start, end, 12),
        xx11 = ts(c(zp6[1:m1, i], zp6[m2, i1]), start, end, 12),
        xx12 = ts(c(zo6[1:m1, i], zo6[m2, i1]), start, end, 12),
        zz1 = ts(c(df1[1:m1, i], df1[m1, i]), start, end, 12),
        zz2 = ts(c(df2[1:m1, i], df2[m1, i]), start, end, 12),
        zz3 = ts(c(df3[1:m1, i], df3[m1, i]), start, end, 12),
        zz4 = ts(c(df4[1:m1, i], df4[m1, i]), start, end, 12),
        zz5 = ts(c(df5[1:m1, i], df5[m1, i]), start, end, 12)
      )
      insample <- 1:m1
      outsample <- (1:length(fulldata$yy))[-insample]
      
      ## Combine MF & CB Recuing window
      recu_combmdxz <- try({
        average_forecast(list(mdz1, mdz2, mdz3, mdz4, mdz5, 
                              mdx1, mdx2, mdx3, mdx4, mdx5, mdx6,
                              mdx7, mdx8, mdx9, mdx10, mdx11, mdx12), 
                         data = fulldata, insample = insample, outsample = outsample, 
                         type = "recursive",
                         measures = c("MSE"),
                         fweights = c("EW", "BICW", "MSFE", "DMSFE"), 
                         show_progress = FALSE)}, silent = T) 
      if ('try-error' %in% class(recu_combmdxz)) {
        recu_forexz <- rep(NA, 17)
      } else {
        recu_forexz <- recu_combmdxz[["forecast"]]
      }
      
      # Collecting
      Out_recu_df_1 <- append(Out_recu_df_1, recu_forexz[1])
      Out_recu_df_2 <- append(Out_recu_df_2, recu_forexz[2])
      Out_recu_df_3 <- append(Out_recu_df_3, recu_forexz[3])
      Out_recu_df_4 <- append(Out_recu_df_4, recu_forexz[4])
      Out_recu_df_5 <- append(Out_recu_df_5, recu_forexz[5])
      Out_recu_zp1 <- append(Out_recu_zp1, recu_forexz[6])
      Out_recu_zo1 <- append(Out_recu_zo1, recu_forexz[7])
      Out_recu_zp2 <- append(Out_recu_zp2, recu_forexz[8])
      Out_recu_zo2 <- append(Out_recu_zo2, recu_forexz[9])
      Out_recu_zp3 <- append(Out_recu_zp3, recu_forexz[10])
      Out_recu_zo3 <- append(Out_recu_zo3, recu_forexz[11])
      Out_recu_zp4 <- append(Out_recu_zp4, recu_forexz[12])
      Out_recu_zo4 <- append(Out_recu_zo4, recu_forexz[13])
      Out_recu_zp5 <- append(Out_recu_zp5, recu_forexz[14])
      Out_recu_zo5 <- append(Out_recu_zo5, recu_forexz[15])
      Out_recu_zp6 <- append(Out_recu_zp6, recu_forexz[16])
      Out_recu_zo6 <- append(Out_recu_zo6, recu_forexz[17])
    }
    
    Rmse_df1 <- sqrt(sum((Out_recu_df_1 - tag[-c(1:leny)])^2) / length(Out_recu_df_1))
    Rmse_df2 <- sqrt(sum((Out_recu_df_2 - tag[-c(1:leny)])^2) / length(Out_recu_df_2))
    Rmse_df3 <- sqrt(sum((Out_recu_df_3 - tag[-c(1:leny)])^2) / length(Out_recu_df_3))
    Rmse_df4 <- sqrt(sum((Out_recu_df_4 - tag[-c(1:leny)])^2) / length(Out_recu_df_4))
    Rmse_df5 <- sqrt(sum((Out_recu_df_5 - tag[-c(1:leny)])^2) / length(Out_recu_df_5))
    Rmse_zpos_gdp <- sqrt(sum((Out_recu_zp1 - tag[-c(1:leny)])^2) / length(Out_recu_zp1))
    Rmse_zzero_gdp <- sqrt(sum((Out_recu_zo1 - tag[-c(1:leny)])^2) / length(Out_recu_zo1))
    Rmse_zpos_cpi <- sqrt(sum((Out_recu_zp2 - tag[-c(1:leny)])^2) / length(Out_recu_zp2))
    Rmse_zzero_cpi <- sqrt(sum((Out_recu_zo2 - tag[-c(1:leny)])^2) / length(Out_recu_zo2))
    Rmse_zpos_iva <- sqrt(sum((Out_recu_zp3 - tag[-c(1:leny)])^2) / length(Out_recu_zp3))
    Rmse_zzero_iva <- sqrt(sum((Out_recu_zo3 - tag[-c(1:leny)])^2) / length(Out_recu_zo3))
    Rmse_zpos_fai <- sqrt(sum((Out_recu_zp4 - tag[-c(1:leny)])^2) / length(Out_recu_zp4))
    Rmse_zzero_fai <- sqrt(sum((Out_recu_zo4 - tag[-c(1:leny)])^2) / length(Out_recu_zo4))
    Rmse_zpos_m2 <- sqrt(sum((Out_recu_zp5 - tag[-c(1:leny)])^2) / length(Out_recu_zp5))
    Rmse_zzero_m2 <- sqrt(sum((Out_recu_zo5 - tag[-c(1:leny)])^2) / length(Out_recu_zo5))
    Rmse_zpos_ppi <- sqrt(sum((Out_recu_zp6 - tag[-c(1:leny)])^2) / length(Out_recu_zp6))
    Rmse_zzero_ppi <- sqrt(sum((Out_recu_zo6 - tag[-c(1:leny)])^2) / length(Out_recu_zo6))
    
    result_rec_for3$Recu_df1 <- append(result_rec_for3$Recu_df1, Rmse_df1)
    result_rec_for3$Recu_df2 <- append(result_rec_for3$Recu_df2, Rmse_df2)
    result_rec_for3$Recu_df3 <- append(result_rec_for3$Recu_df3, Rmse_df3)
    result_rec_for3$Recu_df4 <- append(result_rec_for3$Recu_df4, Rmse_df4)
    result_rec_for3$Recu_df5 <- append(result_rec_for3$Recu_df5, Rmse_df5)
    result_rec_for3$Recu_zpos_gdp <- append(result_rec_for3$Recu_zpos_gdp, Rmse_zpos_gdp)
    result_rec_for3$Recu_zzero_gdp <- append(result_rec_for3$Recu_zzero_gdp, Rmse_zzero_gdp)
    result_rec_for3$Recu_zpos_cpi <- append(result_rec_for3$Recu_zpos_cpi, Rmse_zpos_cpi)
    result_rec_for3$Recu_zzero_cpi <- append(result_rec_for3$Recu_zzero_cpi, Rmse_zzero_cpi)
    result_rec_for3$Recu_zpos_iva <- append(result_rec_for3$Recu_zpos_iva, Rmse_zpos_iva)
    result_rec_for3$Recu_zzero_iva <- append(result_rec_for3$Recu_zzero_iva, Rmse_zzero_iva)
    result_rec_for3$Recu_zpos_fai <- append(result_rec_for3$Recu_zpos_fai, Rmse_zpos_fai)
    result_rec_for3$Recu_zzero_fai <- append(result_rec_for3$Recu_zzero_fai, Rmse_zzero_fai)
    result_rec_for3$Recu_zpos_m2 <- append(result_rec_for3$Recu_zpos_m2, Rmse_zpos_m2)
    result_rec_for3$Recu_zzero_m2 <- append(result_rec_for3$Recu_zzero_m2, Rmse_zzero_m2)
    result_rec_for3$Recu_zpos_ppi <- append(result_rec_for3$Recu_zpos_ppi, Rmse_zpos_ppi)
    result_rec_for3$Recu_zzero_ppi <- append(result_rec_for3$Recu_zzero_ppi, Rmse_zzero_ppi)
  }
}


#-------------------------------------
# Estimate and Forecast (h=6)
#-------------------------------------
Weightfunction <- c()
Ylag <- c()
Xlag <- c()

result_rec_for6 <- list(Ylag = c(),
                        Xlag = c(),
                        Tau = c(),
                        Recu_df1 = c(),
                        Recu_df2 = c(),
                        Recu_df3 = c(),
                        Recu_df4 = c(),
                        Recu_df5 = c(),
                        Recu_zpos_gdp = c(),
                        Recu_zzero_gdp = c(),
                        Recu_zpos_cpi = c(),
                        Recu_zzero_cpi = c(),
                        Recu_zpos_iva = c(),
                        Recu_zzero_iva = c(),
                        Recu_zpos_fai = c(),
                        Recu_zzero_fai = c(),
                        Recu_zpos_m2 = c(),
                        Recu_zzero_m2 = c(),
                        Recu_zpos_ppi = c(),
                        Recu_zzero_ppi = c())

Ylag.max <- 12
Xlag.max <- 12

for (ylag in seq(1, Ylag.max, 1)) {
  for (xlag in seq(1, Xlag.max, 1)) {
    result_rec_for6$Ylag <- append(result_rec_for6$Ylag, ylag)
    result_rec_for6$Xlag <- append(result_rec_for6$Xlag, xlag)
    result_rec_for6$Tau <- append(result_rec_for6$Tau, 6)
    
    # Setting
    lagy <- 5 + ylag
    lagx <- 5 + xlag
    
    # Training
    ## MF 
    mdz1 <- midas_r(yy ~ mls(yy, 6:lagy, 1) + mls(zz1, 6:lagx, 1), start = NULL) # DF1
    mdz2 <- midas_r(yy ~ mls(yy, 6:lagy, 1) + mls(zz2, 6:lagx, 1), start = NULL) # DF2
    mdz3 <- midas_r(yy ~ mls(yy, 6:lagy, 1) + mls(zz3, 6:lagx, 1), start = NULL) # DF3
    mdz4 <- midas_r(yy ~ mls(yy, 6:lagy, 1) + mls(zz4, 6:lagx, 1), start = NULL) # DF4
    mdz5 <- midas_r(yy ~ mls(yy, 6:lagy, 1) + mls(zz5, 6:lagx, 1), start = NULL) # DF5
    
    ## HDMR: GDP
    mdx1 <- midas_r(yy ~ mls(yy, 6:lagy, 1) + mls(xx1, 6:lagx, 1), start = NULL) # zp1
    mdx2 <- midas_r(yy ~ mls(yy, 6:lagy, 1) + mls(xx2, 6:lagx, 1), start = NULL) # zo1
    ## HDMR: CPI
    mdx3 <- midas_r(yy ~ mls(yy, 6:lagy, 1) + mls(xx3, 6:lagx, 1), start = NULL) # zp2
    mdx4 <- midas_r(yy ~ mls(yy, 6:lagy, 1) + mls(xx4, 6:lagx, 1), start = NULL) # zo2
    ## HDMR: IVA
    mdx5 <- midas_r(yy ~ mls(yy, 6:lagy, 1) + mls(xx5, 6:lagx, 1), start = NULL) # zp3
    mdx6 <- midas_r(yy ~ mls(yy, 6:lagy, 1) + mls(xx6, 6:lagx, 1), start = NULL) # zo3
    ## HDMR: FAI
    mdx7 <- midas_r(yy ~ mls(yy, 6:lagy, 1) + mls(xx7, 6:lagx, 1), start = NULL) # zp4
    mdx8 <- midas_r(yy ~ mls(yy, 6:lagy, 1) + mls(xx8, 6:lagx, 1), start = NULL) # zo4
    ## HDMR: M2
    mdx9 <- midas_r(yy ~ mls(yy, 6:lagy, 1) + mls(xx9, 6:lagx, 1), start = NULL) # zp5
    mdx10 <- midas_r(yy ~ mls(yy, 6:lagy, 1) + mls(xx10, 6:lagx, 1), start = NULL) # zo5
    ## HDMR: PPI
    mdx11 <- midas_r(yy ~ mls(yy, 6:lagy, 1) + mls(xx11, 6:lagx, 1), start = NULL) # zp6
    mdx12 <- midas_r(yy ~ mls(yy, 6:lagy, 1) + mls(xx12, 6:lagx, 1), start = NULL) # zo6
    
    # Predicting
    Out_recu_df_1 <- c()
    Out_recu_df_2 <- c()
    Out_recu_df_3 <- c()
    Out_recu_df_4 <- c()
    Out_recu_df_5 <- c()
    Out_recu_zp1 <- c()
    Out_recu_zo1 <- c()
    Out_recu_zp2 <- c()
    Out_recu_zo2 <- c()
    Out_recu_zp3 <- c()
    Out_recu_zo3 <- c()
    Out_recu_zp4 <- c()
    Out_recu_zo4 <- c()
    Out_recu_zp5 <- c()
    Out_recu_zo5 <- c()
    Out_recu_zp6 <- c()
    Out_recu_zo6 <- c()
    
    for (i in seq(1, (length(zp2)-1), 1)) {
      # Sampling
      i1 = i + 1
      m1 = i + leny - 1
      m2 = m1 + 1
      
      # Start & End
      start = c(2003,1)
      end = start + c(len, i-1)
      fulldata <- list(
        yy = ts(c(tag[1:m1], tag[m2]),start, end, 12),
        xx1 = ts(c(zp1[1:m1, i], zp1[m2, i1]), start, end, 12),
        xx2 = ts(c(zo1[1:m1, i], zo1[m2, i1]), start, end, 12),
        xx3 = ts(c(zp2[1:m1, i], zp2[m2, i1]), start, end, 12),
        xx4 = ts(c(zo2[1:m1, i], zo2[m2, i1]), start, end, 12),
        xx5 = ts(c(zp3[1:m1, i], zp3[m2, i1]), start, end, 12),
        xx6 = ts(c(zo3[1:m1, i], zo3[m2, i1]), start, end, 12),
        xx7 = ts(c(zp4[1:m1, i], zp4[m2, i1]), start, end, 12),
        xx8 = ts(c(zo4[1:m1, i], zo4[m2, i1]), start, end, 12),
        xx9 = ts(c(zp5[1:m1, i], zp5[m2, i1]), start, end, 12),
        xx10 = ts(c(zo5[1:m1, i], zo5[m2, i1]), start, end, 12),
        xx11 = ts(c(zp6[1:m1, i], zp6[m2, i1]), start, end, 12),
        xx12 = ts(c(zo6[1:m1, i], zo6[m2, i1]), start, end, 12),
        zz1 = ts(c(df1[1:m1, i], df1[m1, i]), start, end, 12),
        zz2 = ts(c(df2[1:m1, i], df2[m1, i]), start, end, 12),
        zz3 = ts(c(df3[1:m1, i], df3[m1, i]), start, end, 12),
        zz4 = ts(c(df4[1:m1, i], df4[m1, i]), start, end, 12),
        zz5 = ts(c(df5[1:m1, i], df5[m1, i]), start, end, 12)
      )
      insample <- 1:m1
      outsample <- (1:length(fulldata$yy))[-insample]
      
      ## Combine MF & CB Recuing window
      recu_combmdxz <- try({
        average_forecast(list(mdz1, mdz2, mdz3, mdz4, mdz5, 
                              mdx1, mdx2, mdx3, mdx4, mdx5, mdx6,
                              mdx7, mdx8, mdx9, mdx10, mdx11, mdx12), 
                         data = fulldata, insample = insample, outsample = outsample, 
                         type = "recursive",
                         measures = c("MSE"),
                         fweights = c("EW", "BICW", "MSFE", "DMSFE"), 
                         show_progress = FALSE)}, silent = T) 
      if ('try-error' %in% class(recu_combmdxz)) {
        recu_forexz <- rep(NA, 17)
      } else {
        recu_forexz <- recu_combmdxz[["forecast"]]
      }
      
      # Collecting
      Out_recu_df_1 <- append(Out_recu_df_1, recu_forexz[1])
      Out_recu_df_2 <- append(Out_recu_df_2, recu_forexz[2])
      Out_recu_df_3 <- append(Out_recu_df_3, recu_forexz[3])
      Out_recu_df_4 <- append(Out_recu_df_4, recu_forexz[4])
      Out_recu_df_5 <- append(Out_recu_df_5, recu_forexz[5])
      Out_recu_zp1 <- append(Out_recu_zp1, recu_forexz[6])
      Out_recu_zo1 <- append(Out_recu_zo1, recu_forexz[7])
      Out_recu_zp2 <- append(Out_recu_zp2, recu_forexz[8])
      Out_recu_zo2 <- append(Out_recu_zo2, recu_forexz[9])
      Out_recu_zp3 <- append(Out_recu_zp3, recu_forexz[10])
      Out_recu_zo3 <- append(Out_recu_zo3, recu_forexz[11])
      Out_recu_zp4 <- append(Out_recu_zp4, recu_forexz[12])
      Out_recu_zo4 <- append(Out_recu_zo4, recu_forexz[13])
      Out_recu_zp5 <- append(Out_recu_zp5, recu_forexz[14])
      Out_recu_zo5 <- append(Out_recu_zo5, recu_forexz[15])
      Out_recu_zp6 <- append(Out_recu_zp6, recu_forexz[16])
      Out_recu_zo6 <- append(Out_recu_zo6, recu_forexz[17])
    }
    
    Rmse_df1 <- sqrt(sum((Out_recu_df_1 - tag[-c(1:leny)])^2) / length(Out_recu_df_1))
    Rmse_df2 <- sqrt(sum((Out_recu_df_2 - tag[-c(1:leny)])^2) / length(Out_recu_df_2))
    Rmse_df3 <- sqrt(sum((Out_recu_df_3 - tag[-c(1:leny)])^2) / length(Out_recu_df_3))
    Rmse_df4 <- sqrt(sum((Out_recu_df_4 - tag[-c(1:leny)])^2) / length(Out_recu_df_4))
    Rmse_df5 <- sqrt(sum((Out_recu_df_5 - tag[-c(1:leny)])^2) / length(Out_recu_df_5))
    Rmse_zpos_gdp <- sqrt(sum((Out_recu_zp1 - tag[-c(1:leny)])^2) / length(Out_recu_zp1))
    Rmse_zzero_gdp <- sqrt(sum((Out_recu_zo1 - tag[-c(1:leny)])^2) / length(Out_recu_zo1))
    Rmse_zpos_cpi <- sqrt(sum((Out_recu_zp2 - tag[-c(1:leny)])^2) / length(Out_recu_zp2))
    Rmse_zzero_cpi <- sqrt(sum((Out_recu_zo2 - tag[-c(1:leny)])^2) / length(Out_recu_zo2))
    Rmse_zpos_iva <- sqrt(sum((Out_recu_zp3 - tag[-c(1:leny)])^2) / length(Out_recu_zp3))
    Rmse_zzero_iva <- sqrt(sum((Out_recu_zo3 - tag[-c(1:leny)])^2) / length(Out_recu_zo3))
    Rmse_zpos_fai <- sqrt(sum((Out_recu_zp4 - tag[-c(1:leny)])^2) / length(Out_recu_zp4))
    Rmse_zzero_fai <- sqrt(sum((Out_recu_zo4 - tag[-c(1:leny)])^2) / length(Out_recu_zo4))
    Rmse_zpos_m2 <- sqrt(sum((Out_recu_zp5 - tag[-c(1:leny)])^2) / length(Out_recu_zp5))
    Rmse_zzero_m2 <- sqrt(sum((Out_recu_zo5 - tag[-c(1:leny)])^2) / length(Out_recu_zo5))
    Rmse_zpos_ppi <- sqrt(sum((Out_recu_zp6 - tag[-c(1:leny)])^2) / length(Out_recu_zp6))
    Rmse_zzero_ppi <- sqrt(sum((Out_recu_zo6 - tag[-c(1:leny)])^2) / length(Out_recu_zo6))
    
    result_rec_for6$Recu_df1 <- append(result_rec_for6$Recu_df1, Rmse_df1)
    result_rec_for6$Recu_df2 <- append(result_rec_for6$Recu_df2, Rmse_df2)
    result_rec_for6$Recu_df3 <- append(result_rec_for6$Recu_df3, Rmse_df3)
    result_rec_for6$Recu_df4 <- append(result_rec_for6$Recu_df4, Rmse_df4)
    result_rec_for6$Recu_df5 <- append(result_rec_for6$Recu_df5, Rmse_df5)
    result_rec_for6$Recu_zpos_gdp <- append(result_rec_for6$Recu_zpos_gdp, Rmse_zpos_gdp)
    result_rec_for6$Recu_zzero_gdp <- append(result_rec_for6$Recu_zzero_gdp, Rmse_zzero_gdp)
    result_rec_for6$Recu_zpos_cpi <- append(result_rec_for6$Recu_zpos_cpi, Rmse_zpos_cpi)
    result_rec_for6$Recu_zzero_cpi <- append(result_rec_for6$Recu_zzero_cpi, Rmse_zzero_cpi)
    result_rec_for6$Recu_zpos_iva <- append(result_rec_for6$Recu_zpos_iva, Rmse_zpos_iva)
    result_rec_for6$Recu_zzero_iva <- append(result_rec_for6$Recu_zzero_iva, Rmse_zzero_iva)
    result_rec_for6$Recu_zpos_fai <- append(result_rec_for6$Recu_zpos_fai, Rmse_zpos_fai)
    result_rec_for6$Recu_zzero_fai <- append(result_rec_for6$Recu_zzero_fai, Rmse_zzero_fai)
    result_rec_for6$Recu_zpos_m2 <- append(result_rec_for6$Recu_zpos_m2, Rmse_zpos_m2)
    result_rec_for6$Recu_zzero_m2 <- append(result_rec_for6$Recu_zzero_m2, Rmse_zzero_m2)
    result_rec_for6$Recu_zpos_ppi <- append(result_rec_for6$Recu_zpos_ppi, Rmse_zpos_ppi)
    result_rec_for6$Recu_zzero_ppi <- append(result_rec_for6$Recu_zzero_ppi, Rmse_zzero_ppi)
  }
}


#-------------------------------------
# Estimate and Forecast (h=12)
#-------------------------------------
Weightfunction <- c()
Ylag <- c()
Xlag <- c()

result_rec_for12 <- list(Ylag = c(),
                         Xlag = c(),
                         Tau = c(),
                         Recu_df1 = c(),
                         Recu_df2 = c(),
                         Recu_df3 = c(),
                         Recu_df4 = c(),
                         Recu_df5 = c(),
                         Recu_zpos_gdp = c(),
                         Recu_zzero_gdp = c(),
                         Recu_zpos_cpi = c(),
                         Recu_zzero_cpi = c(),
                         Recu_zpos_iva = c(),
                         Recu_zzero_iva = c(),
                         Recu_zpos_fai = c(),
                         Recu_zzero_fai = c(),
                         Recu_zpos_m2 = c(),
                         Recu_zzero_m2 = c(),
                         Recu_zpos_ppi = c(),
                         Recu_zzero_ppi = c())

Ylag.max <- 12
Xlag.max <- 12

for (ylag in seq(1, Ylag.max, 1)) {
  for (xlag in seq(1, Xlag.max, 1)) {
    result_rec_for12$Ylag <- append(result_rec_for12$Ylag, ylag)
    result_rec_for12$Xlag <- append(result_rec_for12$Xlag, xlag)
    result_rec_for12$Tau <- append(result_rec_for12$Tau, 12)
    
    # Setting
    lagy <- 11 + ylag
    lagx <- 11 + xlag    
    
    # Training
    ## MF 
    mdz1 <- midas_r(yy ~ mls(yy, 12:lagy, 1) + mls(zz1, 12:lagx, 1), start = NULL) # DF1
    mdz2 <- midas_r(yy ~ mls(yy, 12:lagy, 1) + mls(zz2, 12:lagx, 1), start = NULL) # DF2
    mdz3 <- midas_r(yy ~ mls(yy, 12:lagy, 1) + mls(zz3, 12:lagx, 1), start = NULL) # DF3
    mdz4 <- midas_r(yy ~ mls(yy, 12:lagy, 1) + mls(zz4, 12:lagx, 1), start = NULL) # DF4
    mdz5 <- midas_r(yy ~ mls(yy, 12:lagy, 1) + mls(zz5, 12:lagx, 1), start = NULL) # DF5
    
    ## HDMR: GDP
    mdx1 <- midas_r(yy ~ mls(yy, 12:lagy, 1) + mls(xx1, 12:lagx, 1), start = NULL) # zp1
    mdx2 <- midas_r(yy ~ mls(yy, 12:lagy, 1) + mls(xx2, 12:lagx, 1), start = NULL) # zo1
    ## HDMR: CPI
    mdx3 <- midas_r(yy ~ mls(yy, 12:lagy, 1) + mls(xx3, 12:lagx, 1), start = NULL) # zp2
    mdx4 <- midas_r(yy ~ mls(yy, 12:lagy, 1) + mls(xx4, 12:lagx, 1), start = NULL) # zo2
    ## HDMR: IVA
    mdx5 <- midas_r(yy ~ mls(yy, 12:lagy, 1) + mls(xx5, 12:lagx, 1), start = NULL) # zp3
    mdx6 <- midas_r(yy ~ mls(yy, 12:lagy, 1) + mls(xx6, 12:lagx, 1), start = NULL) # zo3
    ## HDMR: FAI
    mdx7 <- midas_r(yy ~ mls(yy, 12:lagy, 1) + mls(xx7, 12:lagx, 1), start = NULL) # zp4
    mdx8 <- midas_r(yy ~ mls(yy, 12:lagy, 1) + mls(xx8, 12:lagx, 1), start = NULL) # zo4
    ## HDMR: M2
    mdx9 <- midas_r(yy ~ mls(yy, 12:lagy, 1) + mls(xx9, 12:lagx, 1), start = NULL) # zp5
    mdx10 <- midas_r(yy ~ mls(yy, 12:lagy, 1) + mls(xx10, 12:lagx, 1), start = NULL) # zo5
    ## HDMR: PPI
    mdx11 <- midas_r(yy ~ mls(yy, 12:lagy, 1) + mls(xx11, 12:lagx, 1), start = NULL) # zp6
    mdx12 <- midas_r(yy ~ mls(yy, 12:ylag, 1) + mls(xx12, 12:lagx, 1), start = NULL) # zo6
    
    # Predicting
    Out_recu_df_1 <- c()
    Out_recu_df_2 <- c()
    Out_recu_df_3 <- c()
    Out_recu_df_4 <- c()
    Out_recu_df_5 <- c()
    Out_recu_zp1 <- c()
    Out_recu_zo1 <- c()
    Out_recu_zp2 <- c()
    Out_recu_zo2 <- c()
    Out_recu_zp3 <- c()
    Out_recu_zo3 <- c()
    Out_recu_zp4 <- c()
    Out_recu_zo4 <- c()
    Out_recu_zp5 <- c()
    Out_recu_zo5 <- c()
    Out_recu_zp6 <- c()
    Out_recu_zo6 <- c()
    
    for (i in seq(1, (length(zp2)-1), 1)) {
      # Sampling
      i1 = i + 1
      m1 = i + leny - 1
      m2 = m1 + 1
      
      # Start & End
      start = c(2003,1)
      end = start + c(len, i-1)
      fulldata <- list(
        yy = ts(c(tag[1:m1], tag[m2]),start, end, 12),
        xx1 = ts(c(zp1[1:m1, i], zp1[m2, i1]), start, end, 12),
        xx2 = ts(c(zo1[1:m1, i], zo1[m2, i1]), start, end, 12),
        xx3 = ts(c(zp2[1:m1, i], zp2[m2, i1]), start, end, 12),
        xx4 = ts(c(zo2[1:m1, i], zo2[m2, i1]), start, end, 12),
        xx5 = ts(c(zp3[1:m1, i], zp3[m2, i1]), start, end, 12),
        xx6 = ts(c(zo3[1:m1, i], zo3[m2, i1]), start, end, 12),
        xx7 = ts(c(zp4[1:m1, i], zp4[m2, i1]), start, end, 12),
        xx8 = ts(c(zo4[1:m1, i], zo4[m2, i1]), start, end, 12),
        xx9 = ts(c(zp5[1:m1, i], zp5[m2, i1]), start, end, 12),
        xx10 = ts(c(zo5[1:m1, i], zo5[m2, i1]), start, end, 12),
        xx11 = ts(c(zp6[1:m1, i], zp6[m2, i1]), start, end, 12),
        xx12 = ts(c(zo6[1:m1, i], zo6[m2, i1]), start, end, 12),
        zz1 = ts(c(df1[1:m1, i], df1[m1, i]), start, end, 12),
        zz2 = ts(c(df2[1:m1, i], df2[m1, i]), start, end, 12),
        zz3 = ts(c(df3[1:m1, i], df3[m1, i]), start, end, 12),
        zz4 = ts(c(df4[1:m1, i], df4[m1, i]), start, end, 12),
        zz5 = ts(c(df5[1:m1, i], df5[m1, i]), start, end, 12)
      )
      insample <- 1:m1
      outsample <- (1:length(fulldata$yy))[-insample]
      
      ## Combine MF & CB Recuing window
      recu_combmdxz <- try({
        average_forecast(list(mdz1, mdz2, mdz3, mdz4, mdz5, 
                              mdx1, mdx2, mdx3, mdx4, mdx5, mdx6,
                              mdx7, mdx8, mdx9, mdx10, mdx11, mdx12), 
                         data = fulldata, insample = insample, outsample = outsample, 
                         type = "recursive",
                         measures = c("MSE"),
                         fweights = c("EW", "BICW", "MSFE", "DMSFE"), 
                         show_progress = FALSE)}, silent = T) 
      if ('try-error' %in% class(recu_combmdxz)) {
        recu_forexz <- rep(NA, 17)
      } else {
        recu_forexz <- recu_combmdxz[["forecast"]]
      }
      
      # Collecting
      Out_recu_df_1 <- append(Out_recu_df_1, recu_forexz[1])
      Out_recu_df_2 <- append(Out_recu_df_2, recu_forexz[2])
      Out_recu_df_3 <- append(Out_recu_df_3, recu_forexz[3])
      Out_recu_df_4 <- append(Out_recu_df_4, recu_forexz[4])
      Out_recu_df_5 <- append(Out_recu_df_5, recu_forexz[5])
      Out_recu_zp1 <- append(Out_recu_zp1, recu_forexz[6])
      Out_recu_zo1 <- append(Out_recu_zo1, recu_forexz[7])
      Out_recu_zp2 <- append(Out_recu_zp2, recu_forexz[8])
      Out_recu_zo2 <- append(Out_recu_zo2, recu_forexz[9])
      Out_recu_zp3 <- append(Out_recu_zp3, recu_forexz[10])
      Out_recu_zo3 <- append(Out_recu_zo3, recu_forexz[11])
      Out_recu_zp4 <- append(Out_recu_zp4, recu_forexz[12])
      Out_recu_zo4 <- append(Out_recu_zo4, recu_forexz[13])
      Out_recu_zp5 <- append(Out_recu_zp5, recu_forexz[14])
      Out_recu_zo5 <- append(Out_recu_zo5, recu_forexz[15])
      Out_recu_zp6 <- append(Out_recu_zp6, recu_forexz[16])
      Out_recu_zo6 <- append(Out_recu_zo6, recu_forexz[17])
    }
    
    Rmse_df1 <- sqrt(sum((Out_recu_df_1 - tag[-c(1:leny)])^2) / length(Out_recu_df_1))
    Rmse_df2 <- sqrt(sum((Out_recu_df_2 - tag[-c(1:leny)])^2) / length(Out_recu_df_2))
    Rmse_df3 <- sqrt(sum((Out_recu_df_3 - tag[-c(1:leny)])^2) / length(Out_recu_df_3))
    Rmse_df4 <- sqrt(sum((Out_recu_df_4 - tag[-c(1:leny)])^2) / length(Out_recu_df_4))
    Rmse_df5 <- sqrt(sum((Out_recu_df_5 - tag[-c(1:leny)])^2) / length(Out_recu_df_5))
    Rmse_zpos_gdp <- sqrt(sum((Out_recu_zp1 - tag[-c(1:leny)])^2) / length(Out_recu_zp1))
    Rmse_zzero_gdp <- sqrt(sum((Out_recu_zo1 - tag[-c(1:leny)])^2) / length(Out_recu_zo1))
    Rmse_zpos_cpi <- sqrt(sum((Out_recu_zp2 - tag[-c(1:leny)])^2) / length(Out_recu_zp2))
    Rmse_zzero_cpi <- sqrt(sum((Out_recu_zo2 - tag[-c(1:leny)])^2) / length(Out_recu_zo2))
    Rmse_zpos_iva <- sqrt(sum((Out_recu_zp3 - tag[-c(1:leny)])^2) / length(Out_recu_zp3))
    Rmse_zzero_iva <- sqrt(sum((Out_recu_zo3 - tag[-c(1:leny)])^2) / length(Out_recu_zo3))
    Rmse_zpos_fai <- sqrt(sum((Out_recu_zp4 - tag[-c(1:leny)])^2) / length(Out_recu_zp4))
    Rmse_zzero_fai <- sqrt(sum((Out_recu_zo4 - tag[-c(1:leny)])^2) / length(Out_recu_zo4))
    Rmse_zpos_m2 <- sqrt(sum((Out_recu_zp5 - tag[-c(1:leny)])^2) / length(Out_recu_zp5))
    Rmse_zzero_m2 <- sqrt(sum((Out_recu_zo5 - tag[-c(1:leny)])^2) / length(Out_recu_zo5))
    Rmse_zpos_ppi <- sqrt(sum((Out_recu_zp6 - tag[-c(1:leny)])^2) / length(Out_recu_zp6))
    Rmse_zzero_ppi <- sqrt(sum((Out_recu_zo6 - tag[-c(1:leny)])^2) / length(Out_recu_zo6))
    
    result_rec_for12$Recu_df1 <- append(result_rec_for12$Recu_df1, Rmse_df1)
    result_rec_for12$Recu_df2 <- append(result_rec_for12$Recu_df2, Rmse_df2)
    result_rec_for12$Recu_df3 <- append(result_rec_for12$Recu_df3, Rmse_df3)
    result_rec_for12$Recu_df4 <- append(result_rec_for12$Recu_df4, Rmse_df4)
    result_rec_for12$Recu_df5 <- append(result_rec_for12$Recu_df5, Rmse_df5)
    result_rec_for12$Recu_zpos_gdp <- append(result_rec_for12$Recu_zpos_gdp, Rmse_zpos_gdp)
    result_rec_for12$Recu_zzero_gdp <- append(result_rec_for12$Recu_zzero_gdp, Rmse_zzero_gdp)
    result_rec_for12$Recu_zpos_cpi <- append(result_rec_for12$Recu_zpos_cpi, Rmse_zpos_cpi)
    result_rec_for12$Recu_zzero_cpi <- append(result_rec_for12$Recu_zzero_cpi, Rmse_zzero_cpi)
    result_rec_for12$Recu_zpos_iva <- append(result_rec_for12$Recu_zpos_iva, Rmse_zpos_iva)
    result_rec_for12$Recu_zzero_iva <- append(result_rec_for12$Recu_zzero_iva, Rmse_zzero_iva)
    result_rec_for12$Recu_zpos_fai <- append(result_rec_for12$Recu_zpos_fai, Rmse_zpos_fai)
    result_rec_for12$Recu_zzero_fai <- append(result_rec_for12$Recu_zzero_fai, Rmse_zzero_fai)
    result_rec_for12$Recu_zpos_m2 <- append(result_rec_for12$Recu_zpos_m2, Rmse_zpos_m2)
    result_rec_for12$Recu_zzero_m2 <- append(result_rec_for12$Recu_zzero_m2, Rmse_zzero_m2)
    result_rec_for12$Recu_zpos_ppi <- append(result_rec_for12$Recu_zpos_ppi, Rmse_zpos_ppi)
    result_rec_for12$Recu_zzero_ppi <- append(result_rec_for12$Recu_zzero_ppi, Rmse_zzero_ppi)
  }
}


#-------------------------------------
# Collect Results
#-------------------------------------
result_rec_now  <- data.frame(result_rec_now)
result_rec_for1 <- data.frame(result_rec_for1)
result_rec_for2 <- data.frame(result_rec_for3)
result_rec_for3 <- data.frame(result_rec_for6)
result_rec_for4 <- data.frame(result_rec_for12)
sheets <- list("sheet1" = result_rec_now, "sheet2" = result_rec_for1, "sheet3" = result_rec_for3,
               "sheet4" = result_rec_for6, "sheet5" = result_rec_for12)
write.xlsx(sheets, paste0(tmpdatpath, "Forecast_PPI_Realtime.xlsx"), sheetName=c("rec_now", "rec_for1", "rec_for3", "rec_for6", "rec_for12"))


#================================   Step (3)  =================================#
#===============================================================================
# Nowcasting
#===============================================================================
#-------------------------------------
# Load Data
#-------------------------------------
# Target
tag <- na.omit(read.csv(paste0(rawdatpath, 'Target_PPI.csv'))[,2])

# PCIs: GDP
hdmr_gdp <- read.csv(paste0(tmpdatpath, 'PCIs_GDP_Realtime_Nowcast.csv'), header = FALSE)
zp1 <- hdmr_gdp[, seq(1, (length(hdmr_gdp)-3), 4)]
zo1 <- hdmr_gdp[, seq(2, (length(hdmr_gdp)-2), 4)]

# PCIs: CPI
hdmr_cpi <- read.csv(paste0(tmpdatpath, 'PCIs_CPI_Realtime_Nowcast.csv'), header = FALSE)
zp2 <- hdmr_cpi[, seq(1, (length(hdmr_cpi)-3), 4)]
zo2 <- hdmr_cpi[, seq(2, (length(hdmr_cpi)-2), 4)]

# PCIs: IVA
hdmr_iva <- read.csv(paste0(tmpdatpath, 'PCIs_IVA_Realtime_Nowcast.csv'), header = FALSE)
zp3 <- hdmr_iva[, seq(1, (length(hdmr_iva)-3), 4)]
zo3 <- hdmr_iva[, seq(2, (length(hdmr_iva)-2), 4)]

# PCIs: FAI
hdmr_fai <- read.csv(paste0(tmpdatpath, 'PCIs_FAI_Realtime_Nowcast.csv'), header = FALSE)
zp4 <- hdmr_fai[, seq(1, (length(hdmr_fai)-3), 4)]
zo4 <- hdmr_fai[, seq(2, (length(hdmr_fai)-2), 4)]

# PCIs: M2
hdmr_m2 <- read.csv(paste0(tmpdatpath, 'PCIs_M2_Realtime_Nowcast.csv'), header = FALSE)
zp5 <- hdmr_m2[, seq(1, (length(hdmr_m2)-3), 4)]
zo5 <- hdmr_m2[, seq(2, (length(hdmr_m2)-2), 4)]

# PCIs: PPI
hdmr_ppi <- read.csv(paste0(tmpdatpath, 'PCIs_PPI_Realtime_Nowcast.csv'), header = FALSE)
zp6 <- hdmr_ppi[, seq(1, (length(hdmr_ppi)-3), 4)]
zo6 <- hdmr_ppi[, seq(2, (length(hdmr_ppi)-2), 4)]

# Macro Factors
factor <- read.csv(paste0(tmpdatpath, 'MFactors_Realtime.csv'), header = F)
df1 <- data.frame(factor[, seq(1, (length(factor)-4) ,5)])
df2 <- data.frame(factor[, seq(2, (length(factor)-3), 5)])
df3 <- data.frame(factor[, seq(3, (length(factor)-2), 5)])
df4 <- data.frame(factor[, seq(4, (length(factor)-1), 5)])
df5 <- data.frame(factor[, seq(5, (length(factor)), 5)])


#-------------------------------------
# Set Estimating Windows & Variables
#-------------------------------------
len <- 7            # Initial window size = 7 years
leny <- 4 * len
lenx <- 12 * len
mdate <- format(seq.Date(from = as.Date("2003/01/01",format = "%Y/%m/%d"), by = "month", length.out = length(tag)), format = "%Y/%m")

yy = ts(c(tag[1:leny], tag[leny+1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx1 = ts(c(zp1[1:lenx, 1], zp1[lenx+1, 1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx2 = ts(c(zo1[1:lenx, 1], zo1[lenx+1, 1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx3 = ts(c(zp2[1:lenx, 1], zp2[lenx+1, 1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx4 = ts(c(zo2[1:lenx, 1], zo2[lenx+1, 1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx5 = ts(c(zp3[1:lenx, 1], zp3[lenx+1, 1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx6 = ts(c(zo3[1:lenx, 1], zo3[lenx+1, 1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx7 = ts(c(zp4[1:lenx, 1], zp4[lenx+1, 1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx8 = ts(c(zo4[1:lenx, 1], zo4[lenx+1, 1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx9 = ts(c(zp5[1:lenx, 1], zp5[lenx+1, 1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx10 = ts(c(zo5[1:lenx, 1], zo5[lenx+1, 1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx11 = ts(c(zp6[1:lenx, 1], zp6[lenx+1, 1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx12 = ts(c(zo6[1:lenx, 1], zo6[lenx+1, 1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
zz1 = ts(c(df1[1:lenx, 1], df1[lenx+1, 2]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
zz2 = ts(c(df2[1:lenx, 1], df2[lenx+1, 2]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
zz3 = ts(c(df3[1:lenx, 1], df3[lenx+1, 2]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
zz4 = ts(c(df4[1:lenx, 1], df4[lenx+1, 2]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
zz5 = ts(c(df5[1:lenx, 1], df5[lenx+1, 2]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)


#-------------------------------------
# Model Combination (h=0)
#-------------------------------------
result_rec_now <- na.omit(read_xlsx(paste0(tmpdatpath, "Forecast_PPI_Realtime.xlsx"), sheet = 1, col_names = T))

# Select best individual MIDAS model
best_recu_df1 <- result_rec_now[which(result_rec_now$Recu_df1 == min(result_rec_now$Recu_df1)), ]
best_recu_df2 <- result_rec_now[which(result_rec_now$Recu_df2 == min(result_rec_now$Recu_df2)), ]
best_recu_df3 <- result_rec_now[which(result_rec_now$Recu_df3 == min(result_rec_now$Recu_df3)), ]
best_recu_df4 <- result_rec_now[which(result_rec_now$Recu_df4 == min(result_rec_now$Recu_df4)), ]
best_recu_df5 <- result_rec_now[which(result_rec_now$Recu_df5 == min(result_rec_now$Recu_df5)), ]
best_recu_z_pos_gdp <- result_rec_now[which(result_rec_now$Recu_zpos_gdp == min(result_rec_now$Recu_zpos_gdp)), ][1,]
best_recu_z_zero_gdp <- result_rec_now[which(result_rec_now$Recu_zzero_gdp == min(result_rec_now$Recu_zzero_gdp)), ][1,]
best_recu_z_pos_cpi <- result_rec_now[which(result_rec_now$Recu_zpos_cpi == min(result_rec_now$Recu_zpos_cpi)), ]
best_recu_z_zero_cpi <- result_rec_now[which(result_rec_now$Recu_zzero_cpi == min(result_rec_now$Recu_zzero_cpi)), ]
best_recu_z_pos_iva <- result_rec_now[which(result_rec_now$Recu_zpos_iva == min(result_rec_now$Recu_zpos_iva)), ]
best_recu_z_zero_iva <- result_rec_now[which(result_rec_now$Recu_zzero_iva == min(result_rec_now$Recu_zzero_iva)), ]
best_recu_z_pos_fai <- result_rec_now[which(result_rec_now$Recu_zpos_fai == min(result_rec_now$Recu_zpos_fai)), ]
best_recu_z_zero_fai <- result_rec_now[which(result_rec_now$Recu_zzero_fai == min(result_rec_now$Recu_zzero_fai)), ]
best_recu_z_pos_m2 <- result_rec_now[which(result_rec_now$Recu_zpos_m2 == min(result_rec_now$Recu_zpos_m2)), ]
best_recu_z_zero_m2 <- result_rec_now[which(result_rec_now$Recu_zzero_m2 == min(result_rec_now$Recu_zzero_m2)), ]
best_recu_z_pos_ppi <- result_rec_now[which(result_rec_now$Recu_zpos_ppi == min(result_rec_now$Recu_zpos_ppi)), ]
best_recu_z_zero_ppi <- result_rec_now[which(result_rec_now$Recu_zzero_ppi == min(result_rec_now$Recu_zzero_ppi)), ]
## DF1
best_recu_df1_ylag <- best_recu_df1$Ylag
best_recu_df1_xlag <- best_recu_df1$Xlag
## DF2
best_recu_df2_ylag <- best_recu_df2$Ylag
best_recu_df2_xlag <- best_recu_df2$Xlag
## DF3
best_recu_df3_ylag <- best_recu_df3$Ylag
best_recu_df3_xlag <- best_recu_df3$Xlag
## DF4
best_recu_df4_ylag <- best_recu_df4$Ylag
best_recu_df4_xlag <- best_recu_df4$Xlag
## DF5
best_recu_df5_ylag <- best_recu_df5$Ylag
best_recu_df5_xlag <- best_recu_df5$Xlag
## z_pos_gdp
best_recu_zpos_gdp_ylag <- best_recu_z_pos_gdp$Ylag
best_recu_zpos_gdp_xlag <- best_recu_z_pos_gdp$Xlag
## z_zero_gdp
best_recu_zzero_gdp_ylag <- best_recu_z_zero_gdp$Ylag
best_recu_zzero_gdp_xlag <- best_recu_z_zero_gdp$Xlag
## z_pos_cpi
best_recu_zpos_cpi_ylag <- best_recu_z_pos_cpi$Ylag
best_recu_zpos_cpi_xlag <- best_recu_z_pos_cpi$Xlag
## z_zero_cpi
best_recu_zzero_cpi_ylag <- best_recu_z_zero_cpi$Ylag
best_recu_zzero_cpi_xlag <- best_recu_z_zero_cpi$Xlag
## z_pos_iva
best_recu_zpos_iva_ylag <- best_recu_z_pos_iva$Ylag
best_recu_zpos_iva_xlag <- best_recu_z_pos_iva$Xlag
## z_zero_iva
best_recu_zzero_iva_ylag <- best_recu_z_zero_iva$Ylag
best_recu_zzero_iva_xlag <- best_recu_z_zero_iva$Xlag
## z_pos_fai
best_recu_zpos_fai_ylag <- best_recu_z_pos_fai$Ylag
best_recu_zpos_fai_xlag <- best_recu_z_pos_fai$Xlag
## z_zero_fai
best_recu_zzero_fai_ylag <- best_recu_z_zero_fai$Ylag
best_recu_zzero_fai_xlag <- best_recu_z_zero_fai$Xlag
## z_pos_m2
best_recu_zpos_m2_ylag <- best_recu_z_pos_m2$Ylag
best_recu_zpos_m2_xlag <- best_recu_z_pos_m2$Xlag
## z_zero_m2
best_recu_zzero_m2_ylag <- best_recu_z_zero_m2$Ylag
best_recu_zzero_m2_xlag <- best_recu_z_zero_m2$Xlag
## z_pos_ppi
best_recu_zpos_ppi_ylag <- best_recu_z_pos_ppi$Ylag
best_recu_zpos_ppi_xlag <- best_recu_z_pos_ppi$Xlag
## z_zero_ppi
best_recu_zzero_ppi_ylag <- best_recu_z_zero_ppi$Ylag
best_recu_zzero_ppi_xlag <- best_recu_z_zero_ppi$Xlag

# Training and combine best individual models to get combination MIDAS
## MF
best_md1 <- midas_r(yy ~ mls(yy, 1:best_recu_df1_ylag, 1) + mls(zz1, 1:best_recu_df1_xlag, 1), start = NULL) # DF1
best_md2 <- midas_r(yy ~ mls(yy, 1:best_recu_df2_ylag, 1) + mls(zz2, 1:best_recu_df2_xlag, 1), start = NULL) # DF2
best_md3 <- midas_r(yy ~ mls(yy, 1:best_recu_df3_ylag, 1) + mls(zz3, 1:best_recu_df3_xlag, 1), start = NULL) # DF3
best_md4 <- midas_r(yy ~ mls(yy, 1:best_recu_df4_ylag, 1) + mls(zz4, 1:best_recu_df4_xlag, 1), start = NULL) # DF4
best_md5 <- midas_r(yy ~ mls(yy, 1:best_recu_df5_ylag, 1) + mls(zz5, 1:best_recu_df5_xlag, 1), start = NULL) # DF5
## MF+CB_2
best_md6 <- midas_r(yy ~ mls(yy, 1:best_recu_zpos_gdp_ylag, 1) + mls(xx1, 0:best_recu_zpos_gdp_xlag, 1), start = NULL) # z_pos_gdp
best_md7 <- midas_r(yy ~ mls(yy, 1:best_recu_zzero_gdp_ylag, 1) + mls(xx2, 0:best_recu_zzero_gdp_xlag, 1), start = NULL) # z_zero_gdp
best_md8 <- midas_r(yy ~ mls(yy, 1:best_recu_zpos_cpi_ylag, 1) + mls(xx3, 0:best_recu_zpos_cpi_xlag, 1), start = NULL) # z_pos_cpi
best_md9 <- midas_r(yy ~ mls(yy, 1:best_recu_zzero_cpi_ylag, 1) + mls(xx4, 0:best_recu_zzero_cpi_xlag, 1), start = NULL) # z_zero_cpi
best_md10 <- midas_r(yy ~ mls(yy, 1:best_recu_zpos_iva_ylag, 1) + mls(xx5, 0:best_recu_zpos_iva_xlag, 1), start = NULL) # z_pos_iva
best_md11 <- midas_r(yy ~ mls(yy, 1:best_recu_zzero_iva_ylag, 1) + mls(xx6, 0:best_recu_zzero_iva_xlag, 1), start = NULL) # z_zero_iva
best_md12 <- midas_r(yy ~ mls(yy, 1:best_recu_zpos_fai_ylag, 1) + mls(xx7, 0:best_recu_zpos_fai_xlag, 1), start = NULL) # z_pos_fai
best_md13 <- midas_r(yy ~ mls(yy, 1:best_recu_zzero_fai_ylag, 1) + mls(xx8, 0:best_recu_zzero_fai_xlag, 1), start = NULL) # z_zero_fai
best_md14 <- midas_r(yy ~ mls(yy, 1:best_recu_zpos_m2_ylag, 1) + mls(xx9, 0:best_recu_zpos_m2_xlag, 1), start = NULL) # z_pos_m2
best_md15 <- midas_r(yy ~ mls(yy, 1:best_recu_zzero_m2_ylag, 1) + mls(xx10, 0:best_recu_zzero_m2_xlag, 1), start = NULL) # z_zero_m2
best_md16 <- midas_r(yy ~ mls(yy, 1:best_recu_zpos_ppi_ylag, 1) + mls(xx11, 0:best_recu_zpos_ppi_xlag, 1), start = NULL) # z_pos_ppi
best_md17 <- midas_r(yy ~ mls(yy, 1:best_recu_zzero_ppi_ylag, 1) + mls(xx12, 0:best_recu_zzero_ppi_xlag, 1), start = NULL) # z_zero_ppi

# Estimating & Forecasting
Recu_mf_ew = c()
Recu_mf_bicw = c()
Recu_mf_msfew = c()
Recu_mf_dmsfew = c()
Recu_zp_ew = c()
Recu_zp_bicw = c()
Recu_zp_msfew = c()
Recu_zp_dmsfew = c()
Recu_zz_ew = c()
Recu_zz_bicw = c()
Recu_zz_msfew = c()
Recu_zz_dmsfew = c()
Recu_cb_ew = c()
Recu_cb_bicw = c()
Recu_cb_msfew = c()
Recu_cb_dmsfew = c()
Recu_mfcb_ew = c()
Recu_mfcb_bicw = c()
Recu_mfcb_msfew = c()
Recu_mfcb_dmsfew = c()

for (i in seq(1, length(zp2), 1)) {
  # Sampling
  m1 = i + leny - 1
  m2 = m1 + 1
  
  # Start & End
  start = c(2003,1)
  end = start + c(len, i-1)
  fulldata <- list(
    yy = ts(c(tag[1:m1], tag[m2]),start, end, 12),
    xx1 = ts(c(zp1[1:m1, i], zp1[m2, i]), start, end, 12),
    xx2 = ts(c(zo1[1:m1, i], zo1[m2, i]), start, end, 12),
    xx3 = ts(c(zp2[1:m1, i], zp2[m2, i]), start, end, 12),
    xx4 = ts(c(zo2[1:m1, i], zo2[m2, i]), start, end, 12),
    xx5 = ts(c(zp3[1:m1, i], zp3[m2, i]), start, end, 12),
    xx6 = ts(c(zo3[1:m1, i], zo3[m2, i]), start, end, 12),
    xx7 = ts(c(zp4[1:m1, i], zp4[m2, i]), start, end, 12),
    xx8 = ts(c(zo4[1:m1, i], zo4[m2, i]), start, end, 12),
    xx9 = ts(c(zp5[1:m1, i], zp5[m2, i]), start, end, 12),
    xx10 = ts(c(zo5[1:m1, i], zo5[m2, i]), start, end, 12),
    xx11 = ts(c(zp6[1:m1, i], zp6[m2, i]), start, end, 12),
    xx12 = ts(c(zo6[1:m1, i], zo6[m2, i]), start, end, 12),
    zz1 = ts(c(df1[1:m1, i], df1[m1, i]), start, end, 12),
    zz2 = ts(c(df2[1:m1, i], df2[m1, i]), start, end, 12),
    zz3 = ts(c(df3[1:m1, i], df3[m1, i]), start, end, 12),
    zz4 = ts(c(df4[1:m1, i], df4[m1, i]), start, end, 12),
    zz5 = ts(c(df5[1:m1, i], df5[m1, i]), start, end, 12)
  )
  insample <- 1:m1
  outsample <- (1:length(fulldata$yy))[-insample]
  
  # Combination Model
  ## MF
  best_nowmdzz <- average_forecast(list(best_md1, best_md2, best_md3, best_md4, best_md5),
                                   data = fulldata, insample = insample, outsample = outsample,
                                   type = "recursive", 
                                   measures = c("MSE"),
                                   fweights = c("EW", "BICW", "MSFE", "DMSFE"),
                                   show_progress = FALSE)
  Recu_mf_ew <- append(Recu_mf_ew, best_nowmdzz[["avgforecast"]][1])
  Recu_mf_bicw <- append(Recu_mf_bicw, best_nowmdzz[["avgforecast"]][2])
  Recu_mf_msfew <- append(Recu_mf_msfew, best_nowmdzz[["avgforecast"]][3])
  Recu_mf_dmsfew <- append(Recu_mf_dmsfew, best_nowmdzz[["avgforecast"]][4])
  ## MF+CB_+
  best_nowmdxzp <- average_forecast(list(best_md1, best_md2, best_md3, best_md4, best_md5, best_md12),
                                    data = fulldata, insample = insample, outsample = outsample,
                                    type = "recursive", 
                                    measures = c("MSE"),
                                    fweights = c("EW", "BICW", "MSFE", "DMSFE"),
                                    show_progress = FALSE)
  Recu_zp_ew <- append(Recu_zp_ew, best_nowmdxzp[["avgforecast"]][1])
  Recu_zp_bicw <- append(Recu_zp_bicw, best_nowmdxzp[["avgforecast"]][2])
  Recu_zp_msfew <- append(Recu_zp_msfew, best_nowmdxzp[["avgforecast"]][3])
  Recu_zp_dmsfew <- append(Recu_zp_dmsfew, best_nowmdxzp[["avgforecast"]][4])  
  ## MF+CB_0
  best_nowmdxzz <- average_forecast(list(best_md1, best_md2, best_md3, best_md4, best_md5, best_md13),
                                    data = fulldata, insample = insample, outsample = outsample,
                                    type = "recursive", 
                                    measures = c("MSE"),
                                    fweights = c("EW", "BICW", "MSFE", "DMSFE"),
                                    show_progress = FALSE)
  Recu_zz_ew <- append(Recu_zz_ew, best_nowmdxzz[["avgforecast"]][1])
  Recu_zz_bicw <- append(Recu_zz_bicw, best_nowmdxzz[["avgforecast"]][2])
  Recu_zz_msfew <- append(Recu_zz_msfew, best_nowmdxzz[["avgforecast"]][3])
  Recu_zz_dmsfew <- append(Recu_zz_dmsfew, best_nowmdxzz[["avgforecast"]][4])      
  ## MF+CB_2
  best_nowmdxz_2 <- average_forecast(list(best_md1, best_md2, best_md3, best_md4,                                                                best_md5, best_md12, best_md13),
                                     data = fulldata, insample = insample, outsample = outsample,
                                     type = "recursive", 
                                     measures = c("MSE"),
                                     fweights = c("EW", "BICW", "MSFE", "DMSFE"),
                                     show_progress = FALSE)
  Recu_cb_ew <- append(Recu_cb_ew, best_nowmdxz_2[["avgforecast"]][1])
  Recu_cb_bicw <- append(Recu_cb_bicw, best_nowmdxz_2[["avgforecast"]][2])
  Recu_cb_msfew <- append(Recu_cb_msfew, best_nowmdxz_2[["avgforecast"]][3])
  Recu_cb_dmsfew <- append(Recu_cb_dmsfew, best_nowmdxz_2[["avgforecast"]][4])  
  ## MF+CB_12
  best_nowmdxz_12 <- average_forecast(list(best_md1, best_md2, best_md3, best_md4,
                                           best_md5, best_md6, best_md7, best_md8, best_md9,
                                           best_md10, best_md11, best_md12, best_md13,
                                           best_md14, best_md15, best_md16, best_md17),
                                      data = fulldata, insample = insample, outsample = outsample,
                                      type = "recursive", 
                                      measures = c("MSE"),
                                      fweights = c("EW", "BICW", "MSFE", "DMSFE"),
                                      show_progress = FALSE)
  Recu_mfcb_ew <- append(Recu_mfcb_ew, best_nowmdxz_12[["avgforecast"]][1])
  Recu_mfcb_bicw <- append(Recu_mfcb_bicw, best_nowmdxz_12[["avgforecast"]][2])
  Recu_mfcb_msfew <- append(Recu_mfcb_msfew, best_nowmdxz_12[["avgforecast"]][3])
  Recu_mfcb_dmsfew <- append(Recu_mfcb_dmsfew, best_nowmdxz_12[["avgforecast"]][4])
}

# Saving predicting results
best_result_rec_now <- data.frame(mdate = mdate,
                                  real = tag,
                                  mf_ew = append(rep(0, leny), Recu_mf_ew),
                                  mf_bicw = append(rep(0, leny), Recu_mf_bicw),
                                  mf_msfew = append(rep(0, leny), Recu_mf_msfew),
                                  mf_dmsfew = append(rep(0, leny), Recu_mf_dmsfew),
                                  zp_ew = append(rep(0, leny), Recu_zp_ew),
                                  zp_bicw = append(rep(0, leny), Recu_zp_bicw),
                                  zp_msfew = append(rep(0, leny), Recu_zp_msfew),
                                  zp_dmsfew = append(rep(0, leny), Recu_zp_dmsfew),
                                  zz_ew = append(rep(0, leny), Recu_zz_ew),
                                  zz_bicw = append(rep(0, leny), Recu_zz_bicw),
                                  zz_msfew = append(rep(0, leny), Recu_zz_msfew),
                                  zz_dmsfew = append(rep(0, leny), Recu_zz_dmsfew),                                       
                                  cb_ew = append(rep(0, leny), Recu_cb_ew),
                                  cb_bicw = append(rep(0, leny), Recu_cb_bicw),
                                  cb_msfew = append(rep(0, leny), Recu_cb_msfew),
                                  cb_dmsfew = append(rep(0, leny), Recu_cb_dmsfew),
                                  mfcb_ew = append(rep(0, leny), Recu_mfcb_ew),
                                  mfcb_bicw = append(rep(0, leny), Recu_mfcb_bicw),
                                  mfcb_msfew = append(rep(0, leny), Recu_mfcb_msfew),
                                  mfcb_dmsfew = append(rep(0, leny), Recu_mfcb_dmsfew))


#===============================================================================
# Forecasting
#===============================================================================
#-------------------------------------
# Load Data
#-------------------------------------
# Target
tag <- na.omit(read.csv(paste0(rawdatpath, 'Target_PPI.csv'))[,2])

# PCIs: GDP
hdmr_gdp <- read.csv(paste0(tmpdatpath, 'PCIs_GDP_Realtime.csv'), header = FALSE)
zp1 <- hdmr_gdp[, seq(1, (length(hdmr_gdp)-3), 4)]
zo1 <- hdmr_gdp[, seq(2, (length(hdmr_gdp)-2), 4)]

# PCIs: CPI
hdmr_cpi <- read.csv(paste0(tmpdatpath, 'PCIs_CPI_Realtime.csv'), header = FALSE)
zp2 <- hdmr_cpi[, seq(1, (length(hdmr_cpi)-3), 4)]
zo2 <- hdmr_cpi[, seq(2, (length(hdmr_cpi)-2), 4)]

# PCIs: IVA
hdmr_iva <- read.csv(paste0(tmpdatpath, 'PCIs_IVA_Realtime.csv'), header = FALSE)
zp3 <- hdmr_iva[, seq(1, (length(hdmr_iva)-3), 4)]
zo3 <- hdmr_iva[, seq(2, (length(hdmr_iva)-2), 4)]

# PCIs: FAI
hdmr_fai <- read.csv(paste0(tmpdatpath, 'PCIs_FAI_Realtime.csv'), header = FALSE)
zp4 <- hdmr_fai[, seq(1, (length(hdmr_fai)-3), 4)]
zo4 <- hdmr_fai[, seq(2, (length(hdmr_fai)-2), 4)]

# PCIs: M2
hdmr_m2 <- read.csv(paste0(tmpdatpath, 'PCIs_M2_Realtime.csv'), header = FALSE)
zp5 <- hdmr_m2[, seq(1, (length(hdmr_m2)-3), 4)]
zo5 <- hdmr_m2[, seq(2, (length(hdmr_m2)-2), 4)]

# PCIs: PPI
hdmr_ppi <- read.csv(paste0(tmpdatpath, 'PCIs_PPI_Realtime.csv'), header = FALSE)
zp6 <- hdmr_ppi[, seq(1, (length(hdmr_ppi)-3), 4)]
zo6 <- hdmr_ppi[, seq(2, (length(hdmr_ppi)-2), 4)]

# Macro Factors
factor <- read.csv(paste0(tmpdatpath, 'MFactors_Realtime.csv'), header = F)
df1 <- data.frame(factor[, seq(1, (length(factor)-4) ,5)])
df2 <- data.frame(factor[, seq(2, (length(factor)-3), 5)])
df3 <- data.frame(factor[, seq(3, (length(factor)-2), 5)])
df4 <- data.frame(factor[, seq(4, (length(factor)-1), 5)])
df5 <- data.frame(factor[, seq(5, (length(factor)), 5)])


#-------------------------------------
# Set Estimating Windows & Variables
#-------------------------------------
len <- 7            # Initial window size = 7 years
leny <- 4 * len
lenx <- 12 * len
mdate <- format(seq.Date(from = as.Date("2003/01/01",format = "%Y/%m/%d"), by = "month", length.out = length(tag)), format = "%Y/%m")

yy = ts(c(tag[1:leny], tag[leny+1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx1 = ts(c(zp1[1:lenx, 1], zp1[(lenx+1), 2]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx2 = ts(c(zo1[1:lenx, 1], zo1[(lenx+1), 2]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx3 = ts(c(zp2[1:lenx, 1], zp2[(lenx+1), 2]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx4 = ts(c(zo2[1:lenx, 1], zo2[(lenx+1), 2]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx5 = ts(c(zp3[1:lenx, 1], zp3[(lenx+1), 2]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx6 = ts(c(zo3[1:lenx, 1], zo3[(lenx+1), 2]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx7 = ts(c(zp4[1:lenx, 1], zp4[(lenx+1), 2]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx8 = ts(c(zo4[1:lenx, 1], zo4[(lenx+1), 2]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx9 = ts(c(zp5[1:lenx, 1], zp5[(lenx+1), 2]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx10 = ts(c(zo5[1:lenx, 1], zo5[(lenx+1), 2]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx11 = ts(c(zp6[1:lenx, 1], zp6[(lenx+1), 2]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
xx12 = ts(c(zo6[1:lenx, 1], zo6[(lenx+1), 2]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
zz1 = ts(c(df1[1:lenx, 1], df1[lenx, 1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
zz2 = ts(c(df2[1:lenx, 1], df2[lenx, 1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
zz3 = ts(c(df3[1:lenx, 1], df3[lenx, 1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
zz4 = ts(c(df4[1:lenx, 1], df4[lenx, 1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)
zz5 = ts(c(df5[1:lenx, 1], df5[lenx, 1]), start = c(2003,1), end = c(2003,1)+c(len,0), frequency = 12)


#-------------------------------------
# Model Combination (h=1)
#-------------------------------------
result_rec_for1 <- na.omit(read_xlsx(paste0(tmpdatpath, "Forecast_CPI_Realtime.xlsx"), sheet = 2, col_names = T))

# Select best individual MIDAS model
best_recu_df1 <- result_rec_for1[which(result_rec_for1$Recu_df1 == min(result_rec_for1$Recu_df1)), ]
best_recu_df2 <- result_rec_for1[which(result_rec_for1$Recu_df2 == min(result_rec_for1$Recu_df2)), ]
best_recu_df3 <- result_rec_for1[which(result_rec_for1$Recu_df3 == min(result_rec_for1$Recu_df3)), ]
best_recu_df4 <- result_rec_for1[which(result_rec_for1$Recu_df4 == min(result_rec_for1$Recu_df4)), ]
best_recu_df5 <- result_rec_for1[which(result_rec_for1$Recu_df5 == min(result_rec_for1$Recu_df5)), ]
best_recu_z_pos_gdp <- result_rec_for1[which(result_rec_for1$Recu_zpos_gdp == min(result_rec_for1$Recu_zpos_gdp)), ][1,]
best_recu_z_zero_gdp <- result_rec_for1[which(result_rec_for1$Recu_zzero_gdp == min(result_rec_for1$Recu_zzero_gdp)), ][1,]
best_recu_z_pos_cpi <- result_rec_for1[which(result_rec_for1$Recu_zpos_cpi == min(result_rec_for1$Recu_zpos_cpi)), ]
best_recu_z_zero_cpi <- result_rec_for1[which(result_rec_for1$Recu_zzero_cpi == min(result_rec_for1$Recu_zzero_cpi)), ]
best_recu_z_pos_iva <- result_rec_for1[which(result_rec_for1$Recu_zpos_iva == min(result_rec_for1$Recu_zpos_iva)), ]
best_recu_z_zero_iva <- result_rec_for1[which(result_rec_for1$Recu_zzero_iva == min(result_rec_for1$Recu_zzero_iva)), ]
best_recu_z_pos_fai <- result_rec_for1[which(result_rec_for1$Recu_zpos_fai == min(result_rec_for1$Recu_zpos_fai)), ]
best_recu_z_zero_fai <- result_rec_for1[which(result_rec_for1$Recu_zzero_fai == min(result_rec_for1$Recu_zzero_fai)), ]
best_recu_z_pos_m2 <- result_rec_for1[which(result_rec_for1$Recu_zpos_m2 == min(result_rec_for1$Recu_zpos_m2)), ]
best_recu_z_zero_m2 <- result_rec_for1[which(result_rec_for1$Recu_zzero_m2 == min(result_rec_for1$Recu_zzero_m2)), ]
best_recu_z_pos_ppi <- result_rec_for1[which(result_rec_for1$Recu_zpos_ppi == min(result_rec_for1$Recu_zpos_ppi)), ]
best_recu_z_zero_ppi <- result_rec_for1[which(result_rec_for1$Recu_zzero_ppi == min(result_rec_for1$Recu_zzero_ppi)), ]
## DF1
best_recu_df1_ylag <- best_recu_df1$Ylag
best_recu_df1_xlag <- best_recu_df1$Xlag
## DF2
best_recu_df2_ylag <- best_recu_df2$Ylag
best_recu_df2_xlag <- best_recu_df2$Xlag
## DF3
best_recu_df3_ylag <- best_recu_df3$Ylag
best_recu_df3_xlag <- best_recu_df3$Xlag
## DF4
best_recu_df4_ylag <- best_recu_df4$Ylag
best_recu_df4_xlag <- best_recu_df4$Xlag
## DF5
best_recu_df5_ylag <- best_recu_df5$Ylag
best_recu_df5_xlag <- best_recu_df5$Xlag
## z_pos_gdp
best_recu_zpos_gdp_ylag <- best_recu_z_pos_gdp$Ylag
best_recu_zpos_gdp_xlag <- best_recu_z_pos_gdp$Xlag
## z_zero_gdp
best_recu_zzero_gdp_ylag <- best_recu_z_zero_gdp$Ylag
best_recu_zzero_gdp_xlag <- best_recu_z_zero_gdp$Xlag
## z_pos_cpi
best_recu_zpos_cpi_ylag <- best_recu_z_pos_cpi$Ylag
best_recu_zpos_cpi_xlag <- best_recu_z_pos_cpi$Xlag
## z_zero_cpi
best_recu_zzero_cpi_ylag <- best_recu_z_zero_cpi$Ylag
best_recu_zzero_cpi_xlag <- best_recu_z_zero_cpi$Xlag
## z_pos_iva
best_recu_zpos_iva_ylag <- best_recu_z_pos_iva$Ylag
best_recu_zpos_iva_xlag <- best_recu_z_pos_iva$Xlag
## z_zero_iva
best_recu_zzero_iva_ylag <- best_recu_z_zero_iva$Ylag
best_recu_zzero_iva_xlag <- best_recu_z_zero_iva$Xlag
## z_pos_fai
best_recu_zpos_fai_ylag <- best_recu_z_pos_fai$Ylag
best_recu_zpos_fai_xlag <- best_recu_z_pos_fai$Xlag
## z_zero_fai
best_recu_zzero_fai_ylag <- best_recu_z_zero_fai$Ylag
best_recu_zzero_fai_xlag <- best_recu_z_zero_fai$Xlag
## z_pos_m2
best_recu_zpos_m2_ylag <- best_recu_z_pos_m2$Ylag
best_recu_zpos_m2_xlag <- best_recu_z_pos_m2$Xlag
## z_zero_m2
best_recu_zzero_m2_ylag <- best_recu_z_zero_m2$Ylag
best_recu_zzero_m2_xlag <- best_recu_z_zero_m2$Xlag
## z_pos_ppi
best_recu_zpos_ppi_ylag <- best_recu_z_pos_ppi$Ylag
best_recu_zpos_ppi_xlag <- best_recu_z_pos_ppi$Xlag
## z_zero_ppi
best_recu_zzero_ppi_ylag <- best_recu_z_zero_ppi$Ylag
best_recu_zzero_ppi_xlag <- best_recu_z_zero_ppi$Xlag

# Training and combine best individual models to get combination MIDAS
## MF
best_md1 <- midas_r(yy ~ mls(yy, 1:best_recu_df1_ylag, 1) + mls(zz1, 1:best_recu_df1_xlag, 1), start = NULL) # DF1
best_md2 <- midas_r(yy ~ mls(yy, 1:best_recu_df2_ylag, 1) + mls(zz2, 1:best_recu_df2_xlag, 1), start = NULL) # DF2
best_md3 <- midas_r(yy ~ mls(yy, 1:best_recu_df3_ylag, 1) + mls(zz3, 1:best_recu_df3_xlag, 1), start = NULL) # DF3
best_md4 <- midas_r(yy ~ mls(yy, 1:best_recu_df4_ylag, 1) + mls(zz4, 1:best_recu_df4_xlag, 1), start = NULL) # DF4
best_md5 <- midas_r(yy ~ mls(yy, 1:best_recu_df5_ylag, 1) + mls(zz5, 1:best_recu_df5_xlag, 1), start = NULL) # DF5
## MF+CB_2
best_md6 <- midas_r(yy ~ mls(yy, 1:best_recu_zpos_gdp_ylag, 1) + mls(xx1, 1:best_recu_zpos_gdp_xlag, 1), start = NULL) # z_pos_gdp
best_md7 <- midas_r(yy ~ mls(yy, 1:best_recu_zzero_gdp_ylag, 1) + mls(xx2, 1:best_recu_zzero_gdp_xlag, 1), start = NULL) # z_zero_gdp
best_md8 <- midas_r(yy ~ mls(yy, 1:best_recu_zpos_cpi_ylag, 1) + mls(xx3, 1:best_recu_zpos_cpi_xlag, 1), start = NULL) # z_pos_cpi
best_md9 <- midas_r(yy ~ mls(yy, 1:best_recu_zzero_cpi_ylag, 1) + mls(xx4, 1:best_recu_zzero_cpi_xlag, 1), start = NULL) # z_zero_cpi
best_md10 <- midas_r(yy ~ mls(yy, 1:best_recu_zpos_iva_ylag, 1) + mls(xx5, 1:best_recu_zpos_iva_xlag, 1), start = NULL) # z_pos_iva
best_md11 <- midas_r(yy ~ mls(yy, 1:best_recu_zzero_iva_ylag, 1) + mls(xx6, 1:best_recu_zzero_iva_xlag, 1), start = NULL) # z_zero_iva
best_md12 <- midas_r(yy ~ mls(yy, 1:best_recu_zpos_fai_ylag, 1) + mls(xx7, 1:best_recu_zpos_fai_xlag, 1), start = NULL) # z_pos_fai
best_md13 <- midas_r(yy ~ mls(yy, 1:best_recu_zzero_fai_ylag, 1) + mls(xx8, 1:best_recu_zzero_fai_xlag, 1), start = NULL) # z_zero_fai
best_md14 <- midas_r(yy ~ mls(yy, 1:best_recu_zpos_m2_ylag, 1) + mls(xx9, 1:best_recu_zpos_m2_xlag, 1), start = NULL) # z_pos_m2
best_md15 <- midas_r(yy ~ mls(yy, 1:best_recu_zzero_m2_ylag, 1) + mls(xx10, 1:best_recu_zzero_m2_xlag, 1), start = NULL) # z_zero_m2
best_md16 <- midas_r(yy ~ mls(yy, 1:best_recu_zpos_ppi_ylag, 1) + mls(xx11, 1:best_recu_zpos_ppi_xlag, 1), start = NULL) # z_pos_ppi
best_md17 <- midas_r(yy ~ mls(yy, 1:best_recu_zzero_ppi_ylag, 1) + mls(xx12, 1:best_recu_zzero_ppi_xlag, 1), start = NULL) # z_zero_ppi

# Estimating & Forecasting
Recu_mf_ew = c()
Recu_mf_bicw = c()
Recu_mf_msfew = c()
Recu_mf_dmsfew = c()
Recu_zp_ew = c()
Recu_zp_bicw = c()
Recu_zp_msfew = c()
Recu_zp_dmsfew = c()
Recu_zz_ew = c()
Recu_zz_bicw = c()
Recu_zz_msfew = c()
Recu_zz_dmsfew = c()
Recu_cb_ew = c()
Recu_cb_bicw = c()
Recu_cb_msfew = c()
Recu_cb_dmsfew = c()
Recu_mfcb_ew = c()
Recu_mfcb_bicw = c()
Recu_mfcb_msfew = c()
Recu_mfcb_dmsfew = c()

for (i in seq(1, (length(zp2)-1), 1)) {
  # Sampling
  i1 = i + 1
  m1 = i + leny - 1
  m2 = m1 + 1
  
  # Start & End
  start = c(2003,1)
  end = start + c(len, i-1)
  fulldata <- list(
    yy = ts(c(tag[1:m1], tag[m2]),start, end, 12),
    xx1 = ts(c(zp1[1:m1, i], zp1[m2, i1]), start, end, 12),
    xx2 = ts(c(zo1[1:m1, i], zo1[m2, i1]), start, end, 12),
    xx3 = ts(c(zp2[1:m1, i], zp2[m2, i1]), start, end, 12),
    xx4 = ts(c(zo2[1:m1, i], zo2[m2, i1]), start, end, 12),
    xx5 = ts(c(zp3[1:m1, i], zp3[m2, i1]), start, end, 12),
    xx6 = ts(c(zo3[1:m1, i], zo3[m2, i1]), start, end, 12),
    xx7 = ts(c(zp4[1:m1, i], zp4[m2, i1]), start, end, 12),
    xx8 = ts(c(zo4[1:m1, i], zo4[m2, i1]), start, end, 12),
    xx9 = ts(c(zp5[1:m1, i], zp5[m2, i1]), start, end, 12),
    xx10 = ts(c(zo5[1:m1, i], zo5[m2, i1]), start, end, 12),
    xx11 = ts(c(zp6[1:m1, i], zp6[m2, i1]), start, end, 12),
    xx12 = ts(c(zo6[1:m1, i], zo6[m2, i1]), start, end, 12),
    zz1 = ts(c(df1[1:m1, i], df1[m1, i]), start, end, 12),
    zz2 = ts(c(df2[1:m1, i], df2[m1, i]), start, end, 12),
    zz3 = ts(c(df3[1:m1, i], df3[m1, i]), start, end, 12),
    zz4 = ts(c(df4[1:m1, i], df4[m1, i]), start, end, 12),
    zz5 = ts(c(df5[1:m1, i], df5[m1, i]), start, end, 12)
  )
  insample <- 1:m1
  outsample <- (1:length(fulldata$yy))[-insample]
  
  # Combination Model
  ## MF
  best_nowmdzz <- average_forecast(list(best_md1, best_md2, best_md3, best_md4, best_md5),
                                   data = fulldata, insample = insample, outsample = outsample,
                                   type = "recursive", 
                                   measures = c("MSE"),
                                   fweights = c("EW", "BICW", "MSFE", "DMSFE"),
                                   show_progress = FALSE)
  Recu_mf_ew <- append(Recu_mf_ew, best_nowmdzz[["avgforecast"]][1])
  Recu_mf_bicw <- append(Recu_mf_bicw, best_nowmdzz[["avgforecast"]][2])
  Recu_mf_msfew <- append(Recu_mf_msfew, best_nowmdzz[["avgforecast"]][3])
  Recu_mf_dmsfew <- append(Recu_mf_dmsfew, best_nowmdzz[["avgforecast"]][4])
  ## MF+CB_+
  best_nowmdxzp <- average_forecast(list(best_md1, best_md2, best_md3, best_md4, best_md5, best_md12),
                                    data = fulldata, insample = insample, outsample = outsample,
                                    type = "recursive", 
                                    measures = c("MSE"),
                                    fweights = c("EW", "BICW", "MSFE", "DMSFE"),
                                    show_progress = FALSE)
  Recu_zp_ew <- append(Recu_zp_ew, best_nowmdxzp[["avgforecast"]][1])
  Recu_zp_bicw <- append(Recu_zp_bicw, best_nowmdxzp[["avgforecast"]][2])
  Recu_zp_msfew <- append(Recu_zp_msfew, best_nowmdxzp[["avgforecast"]][3])
  Recu_zp_dmsfew <- append(Recu_zp_dmsfew, best_nowmdxzp[["avgforecast"]][4])  
  ## MF+CB_0
  best_nowmdxzz <- average_forecast(list(best_md1, best_md2, best_md3, best_md4, best_md5, best_md13),
                                    data = fulldata, insample = insample, outsample = outsample,
                                    type = "recursive", 
                                    measures = c("MSE"),
                                    fweights = c("EW", "BICW", "MSFE", "DMSFE"),
                                    show_progress = FALSE)
  Recu_zz_ew <- append(Recu_zz_ew, best_nowmdxzz[["avgforecast"]][1])
  Recu_zz_bicw <- append(Recu_zz_bicw, best_nowmdxzz[["avgforecast"]][2])
  Recu_zz_msfew <- append(Recu_zz_msfew, best_nowmdxzz[["avgforecast"]][3])
  Recu_zz_dmsfew <- append(Recu_zz_dmsfew, best_nowmdxzz[["avgforecast"]][4])        
  ## MF+CB_2
  best_nowmdxz_2 <- average_forecast(list(best_md1, best_md2, best_md3, best_md4,                                                                best_md5, best_md12, best_md13),
                                     data = fulldata, insample = insample, outsample = outsample,
                                     type = "recursive", 
                                     measures = c("MSE"),
                                     fweights = c("EW", "BICW", "MSFE", "DMSFE"),
                                     show_progress = FALSE)
  Recu_cb_ew <- append(Recu_cb_ew, best_nowmdxz_2[["avgforecast"]][1])
  Recu_cb_bicw <- append(Recu_cb_bicw, best_nowmdxz_2[["avgforecast"]][2])
  Recu_cb_msfew <- append(Recu_cb_msfew, best_nowmdxz_2[["avgforecast"]][3])
  Recu_cb_dmsfew <- append(Recu_cb_dmsfew, best_nowmdxz_2[["avgforecast"]][4])  
  ## MF+CB_12
  best_nowmdxz_12 <- average_forecast(list(best_md1, best_md2, best_md3, best_md4,
                                           best_md5, best_md6, best_md7, best_md8, best_md9,
                                           best_md10, best_md11, best_md12, best_md13,
                                           best_md14, best_md15, best_md16, best_md17),
                                      data = fulldata, insample = insample, outsample = outsample,
                                      type = "recursive", 
                                      measures = c("MSE"),
                                      fweights = c("EW", "BICW", "MSFE", "DMSFE"),
                                      show_progress = FALSE)
  Recu_mfcb_ew <- append(Recu_mfcb_ew, best_nowmdxz_12[["avgforecast"]][1])
  Recu_mfcb_bicw <- append(Recu_mfcb_bicw, best_nowmdxz_12[["avgforecast"]][2])
  Recu_mfcb_msfew <- append(Recu_mfcb_msfew, best_nowmdxz_12[["avgforecast"]][3])
  Recu_mfcb_dmsfew <- append(Recu_mfcb_dmsfew, best_nowmdxz_12[["avgforecast"]][4])
}

# Saving predicting results
best_result_rec_for1 <- data.frame(mdate = mdate,
                                   real = tag,
                                   mf_ew = append(rep(0, leny), Recu_mf_ew),
                                   mf_bicw = append(rep(0, leny), Recu_mf_bicw),
                                   mf_msfew = append(rep(0, leny), Recu_mf_msfew),
                                   mf_dmsfew = append(rep(0, leny), Recu_mf_dmsfew),
                                   zp_ew = append(rep(0, leny), Recu_zp_ew),
                                   zp_bicw = append(rep(0, leny), Recu_zp_bicw),
                                   zp_msfew = append(rep(0, leny), Recu_zp_msfew),
                                   zp_dmsfew = append(rep(0, leny), Recu_zp_dmsfew),
                                   zz_ew = append(rep(0, leny), Recu_zz_ew),
                                   zz_bicw = append(rep(0, leny), Recu_zz_bicw),
                                   zz_msfew = append(rep(0, leny), Recu_zz_msfew),
                                   zz_dmsfew = append(rep(0, leny), Recu_zz_dmsfew),                                   
                                   cb_ew = append(rep(0, leny), Recu_cb_ew),
                                   cb_bicw = append(rep(0, leny), Recu_cb_bicw),
                                   cb_msfew = append(rep(0, leny), Recu_cb_msfew),
                                   cb_dmsfew = append(rep(0, leny), Recu_cb_dmsfew),
                                   mfcb_ew = append(rep(0, leny), Recu_mfcb_ew),
                                   mfcb_bicw = append(rep(0, leny), Recu_mfcb_bicw),
                                   mfcb_msfew = append(rep(0, leny), Recu_mfcb_msfew),
                                   mfcb_dmsfew = append(rep(0, leny), Recu_mfcb_dmsfew))


#-------------------------------------
# Model Combination (h=3)
#-------------------------------------
result_rec_for2 <- na.omit(read_xlsx(paste0(tmpdatpath, "Forecast_PPI_Realtime.xlsx"), sheet = 3, col_names = T))

# Select best individual MIDAS model
best_recu_df1 <- result_rec_for3[which(result_rec_for3$Recu_df1 == min(result_rec_for3$Recu_df1)), ]
best_recu_df2 <- result_rec_for3[which(result_rec_for3$Recu_df2 == min(result_rec_for3$Recu_df2)), ]
best_recu_df3 <- result_rec_for3[which(result_rec_for3$Recu_df3 == min(result_rec_for3$Recu_df3)), ]
best_recu_df4 <- result_rec_for3[which(result_rec_for3$Recu_df4 == min(result_rec_for3$Recu_df4)), ]
best_recu_df5 <- result_rec_for3[which(result_rec_for3$Recu_df5 == min(result_rec_for3$Recu_df5)), ]
best_recu_z_pos_gdp <- result_rec_for3[which(result_rec_for3$Recu_zpos_gdp == min(result_rec_for3$Recu_zpos_gdp)), ][1,]
best_recu_z_zero_gdp <- result_rec_for3[which(result_rec_for3$Recu_zzero_gdp == min(result_rec_for3$Recu_zzero_gdp)), ][1,]
best_recu_z_pos_cpi <- result_rec_for3[which(result_rec_for3$Recu_zpos_cpi == min(result_rec_for3$Recu_zpos_cpi)), ]
best_recu_z_zero_cpi <- result_rec_for3[which(result_rec_for3$Recu_zzero_cpi == min(result_rec_for3$Recu_zzero_cpi)), ]
best_recu_z_pos_iva <- result_rec_for3[which(result_rec_for3$Recu_zpos_iva == min(result_rec_for3$Recu_zpos_iva)), ]
best_recu_z_zero_iva <- result_rec_for3[which(result_rec_for3$Recu_zzero_iva == min(result_rec_for3$Recu_zzero_iva)), ]
best_recu_z_pos_fai <- result_rec_for3[which(result_rec_for3$Recu_zpos_fai == min(result_rec_for3$Recu_zpos_fai)), ]
best_recu_z_zero_fai <- result_rec_for3[which(result_rec_for3$Recu_zzero_fai == min(result_rec_for3$Recu_zzero_fai)), ]
best_recu_z_pos_m2 <- result_rec_for3[which(result_rec_for3$Recu_zpos_m2 == min(result_rec_for3$Recu_zpos_m2)), ]
best_recu_z_zero_m2 <- result_rec_for3[which(result_rec_for3$Recu_zzero_m2 == min(result_rec_for3$Recu_zzero_m2)), ]
best_recu_z_pos_ppi <- result_rec_for3[which(result_rec_for3$Recu_zpos_ppi == min(result_rec_for3$Recu_zpos_ppi)), ]
best_recu_z_zero_ppi <- result_rec_for3[which(result_rec_for3$Recu_zzero_ppi == min(result_rec_for3$Recu_zzero_ppi)), ]
## DF1
best_recu_df1_ylag <- 2 + best_recu_df1$Ylag
best_recu_df1_xlag <- 2 + best_recu_df1$Xlag
## DF2
best_recu_df2_ylag <- 2 + best_recu_df2$Ylag
best_recu_df2_xlag <- 2 + best_recu_df2$Xlag
## DF3
best_recu_df3_ylag <- 2 + best_recu_df3$Ylag
best_recu_df3_xlag <- 2 + best_recu_df3$Xlag
## DF4
best_recu_df4_ylag <- 2 + best_recu_df4$Ylag
best_recu_df4_xlag <- 2 + best_recu_df4$Xlag
## DF5
best_recu_df5_ylag <- 2 + best_recu_df5$Ylag
best_recu_df5_xlag <- 2 + best_recu_df5$Xlag
## z_pos_gdp
best_recu_zpos_gdp_ylag <- 2 + best_recu_z_pos_gdp$Ylag
best_recu_zpos_gdp_xlag <- 2 + best_recu_z_pos_gdp$Xlag
## z_zero_gdp
best_recu_zzero_gdp_ylag <- 2 + best_recu_z_zero_gdp$Ylag
best_recu_zzero_gdp_xlag <- 2 + best_recu_z_zero_gdp$Xlag
## z_pos_cpi
best_recu_zpos_cpi_ylag <- 2 + best_recu_z_pos_cpi$Ylag
best_recu_zpos_cpi_xlag <- 2 + best_recu_z_pos_cpi$Xlag
## z_zero_cpi
best_recu_zzero_cpi_ylag <- 2 + best_recu_z_zero_cpi$Ylag
best_recu_zzero_cpi_xlag <- 2 + best_recu_z_zero_cpi$Xlag
## z_pos_iva
best_recu_zpos_iva_ylag <- 2 + best_recu_z_pos_iva$Ylag
best_recu_zpos_iva_xlag <- 2 + best_recu_z_pos_iva$Xlag
## z_zero_iva
best_recu_zzero_iva_ylag <- 2 + best_recu_z_zero_iva$Ylag
best_recu_zzero_iva_xlag <- 2 + best_recu_z_zero_iva$Xlag
## z_pos_fai
best_recu_zpos_fai_ylag <- 2 + best_recu_z_pos_fai$Ylag
best_recu_zpos_fai_xlag <- 2 + best_recu_z_pos_fai$Xlag
## z_zero_fai
best_recu_zzero_fai_ylag <- 2 + best_recu_z_zero_fai$Ylag
best_recu_zzero_fai_xlag <- 2 + best_recu_z_zero_fai$Xlag
## z_pos_m2
best_recu_zpos_m2_ylag <- 2 + best_recu_z_pos_m2$Ylag
best_recu_zpos_m2_xlag <- 2 + best_recu_z_pos_m2$Xlag
## z_zero_m2
best_recu_zzero_m2_ylag <- 2 + best_recu_z_zero_m2$Ylag
best_recu_zzero_m2_xlag <- 2 + best_recu_z_zero_m2$Xlag
## z_pos_ppi
best_recu_zpos_ppi_ylag <- 2 + best_recu_z_pos_ppi$Ylag
best_recu_zpos_ppi_xlag <- 2 + best_recu_z_pos_ppi$Xlag
## z_zero_ppi
best_recu_zzero_ppi_ylag <- 2 + best_recu_z_zero_ppi$Ylag
best_recu_zzero_ppi_xlag <- 2 + best_recu_z_zero_ppi$Xlag

# Training and combine best individual models to get combination MIDAS
## MF
best_md1 <- midas_r(yy ~ mls(yy, 3:best_recu_df1_ylag, 1) + mls(zz1, 3:best_recu_df1_xlag, 1), start = NULL) # DF1
best_md2 <- midas_r(yy ~ mls(yy, 3:best_recu_df2_ylag, 1) + mls(zz2, 3:best_recu_df2_xlag, 1), start = NULL) # DF2
best_md3 <- midas_r(yy ~ mls(yy, 3:best_recu_df3_ylag, 1) + mls(zz3, 3:best_recu_df3_xlag, 1), start = NULL) # DF3
best_md4 <- midas_r(yy ~ mls(yy, 3:best_recu_df4_ylag, 1) + mls(zz4, 3:best_recu_df4_xlag, 1), start = NULL) # DF4
best_md5 <- midas_r(yy ~ mls(yy, 3:best_recu_df5_ylag, 1) + mls(zz5, 3:best_recu_df5_xlag, 1), start = NULL) # DF5
## MF+CB_2
best_md6 <- midas_r(yy ~ mls(yy, 3:best_recu_zpos_gdp_ylag, 1) + mls(xx1, 3:best_recu_zpos_gdp_xlag, 1), start = NULL) # z_pos_gdp
best_md7 <- midas_r(yy ~ mls(yy, 3:best_recu_zzero_gdp_ylag, 1) + mls(xx2, 3:best_recu_zzero_gdp_xlag, 1), start = NULL) # z_zero_gdp
best_md8 <- midas_r(yy ~ mls(yy, 3:best_recu_zpos_cpi_ylag, 1) + mls(xx3, 3:best_recu_zpos_cpi_xlag, 1), start = NULL) # z_pos_cpi
best_md9 <- midas_r(yy ~ mls(yy, 3:best_recu_zzero_cpi_ylag, 1) + mls(xx4, 3:best_recu_zzero_cpi_xlag, 1), start = NULL) # z_zero_cpi
best_md10 <- midas_r(yy ~ mls(yy, 3:best_recu_zpos_iva_ylag, 1) + mls(xx5, 3:best_recu_zpos_iva_xlag, 1), start = NULL) # z_pos_iva
best_md11 <- midas_r(yy ~ mls(yy, 3:best_recu_zzero_iva_ylag, 1) + mls(xx6, 3:best_recu_zzero_iva_xlag, 1), start = NULL) # z_zero_iva
best_md12 <- midas_r(yy ~ mls(yy, 3:best_recu_zpos_fai_ylag, 1) + mls(xx7, 3:best_recu_zpos_fai_xlag, 1), start = NULL) # z_pos_fai
best_md13 <- midas_r(yy ~ mls(yy, 3:best_recu_zzero_fai_ylag, 1) + mls(xx8, 3:best_recu_zzero_fai_xlag, 1), start = NULL) # z_zero_fai
best_md14 <- midas_r(yy ~ mls(yy, 3:best_recu_zpos_m2_ylag, 1) + mls(xx9, 3:best_recu_zpos_m2_xlag, 1), start = NULL) # z_pos_m2
best_md15 <- midas_r(yy ~ mls(yy, 3:best_recu_zzero_m2_ylag, 1) + mls(xx10, 3:best_recu_zzero_m2_xlag, 1), start = NULL) # z_zero_m2
best_md16 <- midas_r(yy ~ mls(yy, 3:best_recu_zpos_ppi_ylag, 1) + mls(xx11, 3:best_recu_zpos_ppi_xlag, 1), start = NULL) # z_pos_ppi
best_md17 <- midas_r(yy ~ mls(yy, 3:best_recu_zzero_ppi_ylag, 1) + mls(xx12, 3:best_recu_zzero_ppi_xlag, 1), start = NULL) # z_zero_ppi

# Estimating & Forecasting
Recu_mf_ew = c()
Recu_mf_bicw = c()
Recu_mf_msfew = c()
Recu_mf_dmsfew = c()
Recu_zp_ew = c()
Recu_zp_bicw = c()
Recu_zp_msfew = c()
Recu_zp_dmsfew = c()
Recu_zz_ew = c()
Recu_zz_bicw = c()
Recu_zz_msfew = c()
Recu_zz_dmsfew = c()
Recu_cb_ew = c()
Recu_cb_bicw = c()
Recu_cb_msfew = c()
Recu_cb_dmsfew = c()
Recu_mfcb_ew = c()
Recu_mfcb_bicw = c()
Recu_mfcb_msfew = c()
Recu_mfcb_dmsfew = c()

for (i in seq(1, (length(zp2)-1), 1)) {
  # Sampling
  i1 = i + 1
  m1 = i + leny - 1
  m2 = m1 + 1
  
  # Start & End
  start = c(2003,1)
  end = start + c(len, i-1)
  fulldata <- list(
    yy = ts(c(tag[1:m1], tag[m2]),start, end, 12),
    xx1 = ts(c(zp1[1:m1, i], zp1[m2, i1]), start, end, 12),
    xx2 = ts(c(zo1[1:m1, i], zo1[m2, i1]), start, end, 12),
    xx3 = ts(c(zp2[1:m1, i], zp2[m2, i1]), start, end, 12),
    xx4 = ts(c(zo2[1:m1, i], zo2[m2, i1]), start, end, 12),
    xx5 = ts(c(zp3[1:m1, i], zp3[m2, i1]), start, end, 12),
    xx6 = ts(c(zo3[1:m1, i], zo3[m2, i1]), start, end, 12),
    xx7 = ts(c(zp4[1:m1, i], zp4[m2, i1]), start, end, 12),
    xx8 = ts(c(zo4[1:m1, i], zo4[m2, i1]), start, end, 12),
    xx9 = ts(c(zp5[1:m1, i], zp5[m2, i1]), start, end, 12),
    xx10 = ts(c(zo5[1:m1, i], zo5[m2, i1]), start, end, 12),
    xx11 = ts(c(zp6[1:m1, i], zp6[m2, i1]), start, end, 12),
    xx12 = ts(c(zo6[1:m1, i], zo6[m2, i1]), start, end, 12),
    zz1 = ts(c(df1[1:m1, i], df1[m1, i]), start, end, 12),
    zz2 = ts(c(df2[1:m1, i], df2[m1, i]), start, end, 12),
    zz3 = ts(c(df3[1:m1, i], df3[m1, i]), start, end, 12),
    zz4 = ts(c(df4[1:m1, i], df4[m1, i]), start, end, 12),
    zz5 = ts(c(df5[1:m1, i], df5[m1, i]), start, end, 12)
  )
  insample <- 1:m1
  outsample <- (1:length(fulldata$yy))[-insample]
  
  # Combination Model
  ## MF
  best_nowmdzz <- average_forecast(list(best_md1, best_md2, best_md3, best_md4, best_md5),
                                   data = fulldata, insample = insample, outsample = outsample,
                                   type = "recursive", 
                                   measures = c("MSE"),
                                   fweights = c("EW", "BICW", "MSFE", "DMSFE"),
                                   show_progress = FALSE)
  Recu_mf_ew <- append(Recu_mf_ew, best_nowmdzz[["avgforecast"]][1])
  Recu_mf_bicw <- append(Recu_mf_bicw, best_nowmdzz[["avgforecast"]][2])
  Recu_mf_msfew <- append(Recu_mf_msfew, best_nowmdzz[["avgforecast"]][3])
  Recu_mf_dmsfew <- append(Recu_mf_dmsfew, best_nowmdzz[["avgforecast"]][4])
  ## MF+CB_+
  best_nowmdxzp <- average_forecast(list(best_md1, best_md2, best_md3, best_md4, best_md5, best_md12),
                                    data = fulldata, insample = insample, outsample = outsample,
                                    type = "recursive", 
                                    measures = c("MSE"),
                                    fweights = c("EW", "BICW", "MSFE", "DMSFE"),
                                    show_progress = FALSE)
  Recu_zp_ew <- append(Recu_zp_ew, best_nowmdxzp[["avgforecast"]][1])
  Recu_zp_bicw <- append(Recu_zp_bicw, best_nowmdxzp[["avgforecast"]][2])
  Recu_zp_msfew <- append(Recu_zp_msfew, best_nowmdxzp[["avgforecast"]][3])
  Recu_zp_dmsfew <- append(Recu_zp_dmsfew, best_nowmdxzp[["avgforecast"]][4])  
  ## MF+CB_0
  best_nowmdxzz <- average_forecast(list(best_md1, best_md2, best_md3, best_md4, best_md5, best_md13),
                                    data = fulldata, insample = insample, outsample = outsample,
                                    type = "recursive", 
                                    measures = c("MSE"),
                                    fweights = c("EW", "BICW", "MSFE", "DMSFE"),
                                    show_progress = FALSE)
  Recu_zz_ew <- append(Recu_zz_ew, best_nowmdxzz[["avgforecast"]][1])
  Recu_zz_bicw <- append(Recu_zz_bicw, best_nowmdxzz[["avgforecast"]][2])
  Recu_zz_msfew <- append(Recu_zz_msfew, best_nowmdxzz[["avgforecast"]][3])
  Recu_zz_dmsfew <- append(Recu_zz_dmsfew, best_nowmdxzz[["avgforecast"]][4])          
  ## MF+CB_2
  best_nowmdxz_2 <- average_forecast(list(best_md1, best_md2, best_md3, best_md4,                                                                best_md5, best_md12, best_md13),
                                     data = fulldata, insample = insample, outsample = outsample,
                                     type = "recursive", 
                                     measures = c("MSE"),
                                     fweights = c("EW", "BICW", "MSFE", "DMSFE"),
                                     show_progress = FALSE)
  Recu_cb_ew <- append(Recu_cb_ew, best_nowmdxz_2[["avgforecast"]][1])
  Recu_cb_bicw <- append(Recu_cb_bicw, best_nowmdxz_2[["avgforecast"]][2])
  Recu_cb_msfew <- append(Recu_cb_msfew, best_nowmdxz_2[["avgforecast"]][3])
  Recu_cb_dmsfew <- append(Recu_cb_dmsfew, best_nowmdxz_2[["avgforecast"]][4])  
  ## MF+CB_12
  best_nowmdxz_12 <- average_forecast(list(best_md1, best_md2, best_md3, best_md4,
                                           best_md5, best_md6, best_md7, best_md8, best_md9,
                                           best_md10, best_md11, best_md12, best_md13,
                                           best_md14, best_md15, best_md16, best_md17),
                                      data = fulldata, insample = insample, outsample = outsample,
                                      type = "recursive", 
                                      measures = c("MSE"),
                                      fweights = c("EW", "BICW", "MSFE", "DMSFE"),
                                      show_progress = FALSE)
  Recu_mfcb_ew <- append(Recu_mfcb_ew, best_nowmdxz_12[["avgforecast"]][1])
  Recu_mfcb_bicw <- append(Recu_mfcb_bicw, best_nowmdxz_12[["avgforecast"]][2])
  Recu_mfcb_msfew <- append(Recu_mfcb_msfew, best_nowmdxz_12[["avgforecast"]][3])
  Recu_mfcb_dmsfew <- append(Recu_mfcb_dmsfew, best_nowmdxz_12[["avgforecast"]][4])
}

# Saving predicting results
best_result_rec_for3 <- data.frame(mdate = mdate,
                                   real = tag,
                                   mf_ew = append(rep(0, leny), Recu_mf_ew),
                                   mf_bicw = append(rep(0, leny), Recu_mf_bicw),
                                   mf_msfew = append(rep(0, leny), Recu_mf_msfew),
                                   mf_dmsfew = append(rep(0, leny), Recu_mf_dmsfew),
                                   zp_ew = append(rep(0, leny), Recu_zp_ew),
                                   zp_bicw = append(rep(0, leny), Recu_zp_bicw),
                                   zp_msfew = append(rep(0, leny), Recu_zp_msfew),
                                   zp_dmsfew = append(rep(0, leny), Recu_zp_dmsfew),
                                   zz_ew = append(rep(0, leny), Recu_zz_ew),
                                   zz_bicw = append(rep(0, leny), Recu_zz_bicw),
                                   zz_msfew = append(rep(0, leny), Recu_zz_msfew),
                                   zz_dmsfew = append(rep(0, leny), Recu_zz_dmsfew),   
                                   cb_ew = append(rep(0, leny), Recu_cb_ew),
                                   cb_bicw = append(rep(0, leny), Recu_cb_bicw),
                                   cb_msfew = append(rep(0, leny), Recu_cb_msfew),
                                   cb_dmsfew = append(rep(0, leny), Recu_cb_dmsfew),
                                   mfcb_ew = append(rep(0, leny), Recu_mfcb_ew),
                                   mfcb_bicw = append(rep(0, leny), Recu_mfcb_bicw),
                                   mfcb_msfew = append(rep(0, leny), Recu_mfcb_msfew),
                                   mfcb_dmsfew = append(rep(0, leny), Recu_mfcb_dmsfew))


#-------------------------------------
# Model Combination (h=6)
#-------------------------------------
result_rec_for3 <- na.omit(read_xlsx(paste0(tmpdatpath, "Forecast_PPI_Realtime.xlsx"), sheet = 4, col_names = T))

# Select best individual MIDAS model
best_recu_df1 <- result_rec_for6[which(result_rec_for6$Recu_df1 == min(result_rec_for6$Recu_df1)), ]
best_recu_df2 <- result_rec_for6[which(result_rec_for6$Recu_df2 == min(result_rec_for6$Recu_df2)), ]
best_recu_df3 <- result_rec_for6[which(result_rec_for6$Recu_df3 == min(result_rec_for6$Recu_df3)), ]
best_recu_df4 <- result_rec_for6[which(result_rec_for6$Recu_df4 == min(result_rec_for6$Recu_df4)), ]
best_recu_df5 <- result_rec_for6[which(result_rec_for6$Recu_df5 == min(result_rec_for6$Recu_df5)), ]
best_recu_z_pos_gdp <- result_rec_for6[which(result_rec_for6$Recu_zpos_gdp == min(result_rec_for6$Recu_zpos_gdp)), ][1,]
best_recu_z_zero_gdp <- result_rec_for6[which(result_rec_for6$Recu_zzero_gdp == min(result_rec_for6$Recu_zzero_gdp)), ][1,]
best_recu_z_pos_cpi <- result_rec_for6[which(result_rec_for6$Recu_zpos_cpi == min(result_rec_for6$Recu_zpos_cpi)), ]
best_recu_z_zero_cpi <- result_rec_for6[which(result_rec_for6$Recu_zzero_cpi == min(result_rec_for6$Recu_zzero_cpi)), ]
best_recu_z_pos_iva <- result_rec_for6[which(result_rec_for6$Recu_zpos_iva == min(result_rec_for6$Recu_zpos_iva)), ]
best_recu_z_zero_iva <- result_rec_for6[which(result_rec_for6$Recu_zzero_iva == min(result_rec_for6$Recu_zzero_iva)), ]
best_recu_z_pos_fai <- result_rec_for6[which(result_rec_for6$Recu_zpos_fai == min(result_rec_for6$Recu_zpos_fai)), ]
best_recu_z_zero_fai <- result_rec_for6[which(result_rec_for6$Recu_zzero_fai == min(result_rec_for6$Recu_zzero_fai)), ]
best_recu_z_pos_m2 <- result_rec_for6[which(result_rec_for6$Recu_zpos_m2 == min(result_rec_for6$Recu_zpos_m2)), ]
best_recu_z_zero_m2 <- result_rec_for6[which(result_rec_for6$Recu_zzero_m2 == min(result_rec_for6$Recu_zzero_m2)), ]
best_recu_z_pos_ppi <- result_rec_for6[which(result_rec_for6$Recu_zpos_ppi == min(result_rec_for6$Recu_zpos_ppi)), ]
best_recu_z_zero_ppi <- result_rec_for6[which(result_rec_for6$Recu_zzero_ppi == min(result_rec_for6$Recu_zzero_ppi)), ]
## DF1
best_recu_df1_ylag <- 5 + best_recu_df1$Ylag
best_recu_df1_xlag <- 5 + best_recu_df1$Xlag
## DF2
best_recu_df2_ylag <- 5 + best_recu_df2$Ylag
best_recu_df2_xlag <- 5 + best_recu_df2$Xlag
## DF3
best_recu_df3_ylag <- 5 + best_recu_df3$Ylag
best_recu_df3_xlag <- 5 + best_recu_df3$Xlag
## DF4
best_recu_df4_ylag <- 5 + best_recu_df4$Ylag
best_recu_df4_xlag <- 5 + best_recu_df4$Xlag
## DF5
best_recu_df5_ylag <- 5 + best_recu_df5$Ylag
best_recu_df5_xlag <- 5 + best_recu_df5$Xlag
## z_pos_gdp
best_recu_zpos_gdp_ylag <- 5 + best_recu_z_pos_gdp$Ylag
best_recu_zpos_gdp_xlag <- 5 + best_recu_z_pos_gdp$Xlag
## z_zero_gdp
best_recu_zzero_gdp_ylag <- 5 + best_recu_z_zero_gdp$Ylag
best_recu_zzero_gdp_xlag <- 5 + best_recu_z_zero_gdp$Xlag
## z_pos_cpi
best_recu_zpos_cpi_ylag <- 5 + best_recu_z_pos_cpi$Ylag
best_recu_zpos_cpi_xlag <- 5 + best_recu_z_pos_cpi$Xlag
## z_zero_cpi
best_recu_zzero_cpi_ylag <- 5 + best_recu_z_zero_cpi$Ylag
best_recu_zzero_cpi_xlag <- 5 + best_recu_z_zero_cpi$Xlag
## z_pos_iva
best_recu_zpos_iva_ylag <- 5 + best_recu_z_pos_iva$Ylag
best_recu_zpos_iva_xlag <- 5 + best_recu_z_pos_iva$Xlag
## z_zero_iva
best_recu_zzero_iva_ylag <- 5 + best_recu_z_zero_iva$Ylag
best_recu_zzero_iva_xlag <- 5 + best_recu_z_zero_iva$Xlag
## z_pos_fai
best_recu_zpos_fai_ylag <- 5 + best_recu_z_pos_fai$Ylag
best_recu_zpos_fai_xlag <- 5 + best_recu_z_pos_fai$Xlag
## z_zero_fai
best_recu_zzero_fai_ylag <- 5 + best_recu_z_zero_fai$Ylag
best_recu_zzero_fai_xlag <- 5 + best_recu_z_zero_fai$Xlag
## z_pos_m2
best_recu_zpos_m2_ylag <- 5 + best_recu_z_pos_m2$Ylag
best_recu_zpos_m2_xlag <- 5 + best_recu_z_pos_m2$Xlag
## z_zero_m2
best_recu_zzero_m2_ylag <- 5 + best_recu_z_zero_m2$Ylag
best_recu_zzero_m2_xlag <- 5 + best_recu_z_zero_m2$Xlag
## z_pos_ppi
best_recu_zpos_ppi_ylag <- 5 + best_recu_z_pos_ppi$Ylag
best_recu_zpos_ppi_xlag <- 5 + best_recu_z_pos_ppi$Xlag
## z_zero_ppi
best_recu_zzero_ppi_ylag <- 5 + best_recu_z_zero_ppi$Ylag
best_recu_zzero_ppi_xlag <- 5 + best_recu_z_zero_ppi$Xlag

# Training and combine best individual models to get combination MIDAS
## MF
best_md1 <- midas_r(yy ~ mls(yy, 6:best_recu_df1_ylag, 1) + mls(zz1, 6:best_recu_df1_xlag, 1), start = NULL) # DF1
best_md2 <- midas_r(yy ~ mls(yy, 6:best_recu_df2_ylag, 1) + mls(zz2, 6:best_recu_df2_xlag, 1), start = NULL) # DF2
best_md3 <- midas_r(yy ~ mls(yy, 6:best_recu_df3_ylag, 1) + mls(zz3, 6:best_recu_df3_xlag, 1), start = NULL) # DF3
best_md4 <- midas_r(yy ~ mls(yy, 6:best_recu_df4_ylag, 1) + mls(zz4, 6:best_recu_df4_xlag, 1), start = NULL) # DF4
best_md5 <- midas_r(yy ~ mls(yy, 6:best_recu_df5_ylag, 1) + mls(zz5, 6:best_recu_df5_xlag, 1), start = NULL) # DF5
## MF+CB_2
best_md6 <- midas_r(yy ~ mls(yy, 6:best_recu_zpos_gdp_ylag, 1) + mls(xx1, 6:best_recu_zpos_gdp_xlag, 1), start = NULL) # z_pos_gdp
best_md7 <- midas_r(yy ~ mls(yy, 6:best_recu_zzero_gdp_ylag, 1) + mls(xx2, 6:best_recu_zzero_gdp_xlag, 1), start = NULL) # z_zero_gdp
best_md8 <- midas_r(yy ~ mls(yy, 6:best_recu_zpos_cpi_ylag, 1) + mls(xx3, 6:best_recu_zpos_cpi_xlag, 1), start = NULL) # z_pos_cpi
best_md9 <- midas_r(yy ~ mls(yy, 6:best_recu_zzero_cpi_ylag, 1) + mls(xx4, 6:best_recu_zzero_cpi_xlag, 1), start = NULL) # z_zero_cpi
best_md10 <- midas_r(yy ~ mls(yy, 6:best_recu_zpos_iva_ylag, 1) + mls(xx5, 6:best_recu_zpos_iva_xlag, 1), start = NULL) # z_pos_iva
best_md11 <- midas_r(yy ~ mls(yy, 6:best_recu_zzero_iva_ylag, 1) + mls(xx6, 6:best_recu_zzero_iva_xlag, 1), start = NULL) # z_zero_iva
best_md12 <- midas_r(yy ~ mls(yy, 6:best_recu_zpos_fai_ylag, 1) + mls(xx7, 6:best_recu_zpos_fai_xlag, 1), start = NULL) # z_pos_fai
best_md13 <- midas_r(yy ~ mls(yy, 6:best_recu_zzero_fai_ylag, 1) + mls(xx8, 6:best_recu_zzero_fai_xlag, 1), start = NULL) # z_zero_fai
best_md14 <- midas_r(yy ~ mls(yy, 6:best_recu_zpos_m2_ylag, 1) + mls(xx9, 6:best_recu_zpos_m2_xlag, 1), start = NULL) # z_pos_m2
best_md15 <- midas_r(yy ~ mls(yy, 6:best_recu_zzero_m2_ylag, 1) + mls(xx10, 6:best_recu_zzero_m2_xlag, 1), start = NULL) # z_zero_m2
best_md16 <- midas_r(yy ~ mls(yy, 6:best_recu_zpos_ppi_ylag, 1) + mls(xx11, 6:best_recu_zpos_ppi_xlag, 1), start = NULL) # z_pos_ppi
best_md17 <- midas_r(yy ~ mls(yy, 6:best_recu_zzero_ppi_ylag, 1) + mls(xx12, 6:best_recu_zzero_ppi_xlag, 1), start = NULL) # z_zero_ppi

# Estimating & Forecasting
Recu_mf_ew = c()
Recu_mf_bicw = c()
Recu_mf_msfew = c()
Recu_mf_dmsfew = c()
Recu_zp_ew = c()
Recu_zp_bicw = c()
Recu_zp_msfew = c()
Recu_zp_dmsfew = c()
Recu_zz_ew = c()
Recu_zz_bicw = c()
Recu_zz_msfew = c()
Recu_zz_dmsfew = c()
Recu_cb_ew = c()
Recu_cb_bicw = c()
Recu_cb_msfew = c()
Recu_cb_dmsfew = c()
Recu_mfcb_ew = c()
Recu_mfcb_bicw = c()
Recu_mfcb_msfew = c()
Recu_mfcb_dmsfew = c()

for (i in seq(1, (length(zp2)-1), 1)) {
  # Sampling
  i1 = i + 1
  m1 = i + leny - 1
  m2 = m1 + 1
  
  # Start & End
  start = c(2003,1)
  end = start + c(len, i-1)
  fulldata <- list(
    yy = ts(c(tag[1:m1], tag[m2]),start, end, 12),
    xx1 = ts(c(zp1[1:m1, i], zp1[m2, i1]), start, end, 12),
    xx2 = ts(c(zo1[1:m1, i], zo1[m2, i1]), start, end, 12),
    xx3 = ts(c(zp2[1:m1, i], zp2[m2, i1]), start, end, 12),
    xx4 = ts(c(zo2[1:m1, i], zo2[m2, i1]), start, end, 12),
    xx5 = ts(c(zp3[1:m1, i], zp3[m2, i1]), start, end, 12),
    xx6 = ts(c(zo3[1:m1, i], zo3[m2, i1]), start, end, 12),
    xx7 = ts(c(zp4[1:m1, i], zp4[m2, i1]), start, end, 12),
    xx8 = ts(c(zo4[1:m1, i], zo4[m2, i1]), start, end, 12),
    xx9 = ts(c(zp5[1:m1, i], zp5[m2, i1]), start, end, 12),
    xx10 = ts(c(zo5[1:m1, i], zo5[m2, i1]), start, end, 12),
    xx11 = ts(c(zp6[1:m1, i], zp6[m2, i1]), start, end, 12),
    xx12 = ts(c(zo6[1:m1, i], zo6[m2, i1]), start, end, 12),
    zz1 = ts(c(df1[1:m1, i], df1[m1, i]), start, end, 12),
    zz2 = ts(c(df2[1:m1, i], df2[m1, i]), start, end, 12),
    zz3 = ts(c(df3[1:m1, i], df3[m1, i]), start, end, 12),
    zz4 = ts(c(df4[1:m1, i], df4[m1, i]), start, end, 12),
    zz5 = ts(c(df5[1:m1, i], df5[m1, i]), start, end, 12)
  )
  insample <- 1:m1
  outsample <- (1:length(fulldata$yy))[-insample]
  
  # Combination Model
  ## MF
  best_nowmdzz <- average_forecast(list(best_md1, best_md2, best_md3, best_md4, best_md5),
                                   data = fulldata, insample = insample, outsample = outsample,
                                   type = "recursive", 
                                   measures = c("MSE"),
                                   fweights = c("EW", "BICW", "MSFE", "DMSFE"),
                                   show_progress = FALSE)
  Recu_mf_ew <- append(Recu_mf_ew, best_nowmdzz[["avgforecast"]][1])
  Recu_mf_bicw <- append(Recu_mf_bicw, best_nowmdzz[["avgforecast"]][2])
  Recu_mf_msfew <- append(Recu_mf_msfew, best_nowmdzz[["avgforecast"]][3])
  Recu_mf_dmsfew <- append(Recu_mf_dmsfew, best_nowmdzz[["avgforecast"]][4])
  ## MF+CB_+
  best_nowmdxzp <- average_forecast(list(best_md1, best_md2, best_md3, best_md4, best_md5, best_md12),
                                    data = fulldata, insample = insample, outsample = outsample,
                                    type = "recursive", 
                                    measures = c("MSE"),
                                    fweights = c("EW", "BICW", "MSFE", "DMSFE"),
                                    show_progress = FALSE)
  Recu_zp_ew <- append(Recu_zp_ew, best_nowmdxzp[["avgforecast"]][1])
  Recu_zp_bicw <- append(Recu_zp_bicw, best_nowmdxzp[["avgforecast"]][2])
  Recu_zp_msfew <- append(Recu_zp_msfew, best_nowmdxzp[["avgforecast"]][3])
  Recu_zp_dmsfew <- append(Recu_zp_dmsfew, best_nowmdxzp[["avgforecast"]][4])  
  ## MF+CB_0
  best_nowmdxzz <- average_forecast(list(best_md1, best_md2, best_md3, best_md4, best_md5, best_md13),
                                    data = fulldata, insample = insample, outsample = outsample,
                                    type = "recursive", 
                                    measures = c("MSE"),
                                    fweights = c("EW", "BICW", "MSFE", "DMSFE"),
                                    show_progress = FALSE)
  Recu_zz_ew <- append(Recu_zz_ew, best_nowmdxzz[["avgforecast"]][1])
  Recu_zz_bicw <- append(Recu_zz_bicw, best_nowmdxzz[["avgforecast"]][2])
  Recu_zz_msfew <- append(Recu_zz_msfew, best_nowmdxzz[["avgforecast"]][3])
  Recu_zz_dmsfew <- append(Recu_zz_dmsfew, best_nowmdxzz[["avgforecast"]][4])          
  ## MF+CB_2
  best_nowmdxz_2 <- average_forecast(list(best_md1, best_md2, best_md3, best_md4,                                                                best_md5, best_md12, best_md13),
                                     data = fulldata, insample = insample, outsample = outsample,
                                     type = "recursive", 
                                     measures = c("MSE"),
                                     fweights = c("EW", "BICW", "MSFE", "DMSFE"),
                                     show_progress = FALSE)
  Recu_cb_ew <- append(Recu_cb_ew, best_nowmdxz_2[["avgforecast"]][1])
  Recu_cb_bicw <- append(Recu_cb_bicw, best_nowmdxz_2[["avgforecast"]][2])
  Recu_cb_msfew <- append(Recu_cb_msfew, best_nowmdxz_2[["avgforecast"]][3])
  Recu_cb_dmsfew <- append(Recu_cb_dmsfew, best_nowmdxz_2[["avgforecast"]][4])  
  ## MF+CB_12
  best_nowmdxz_12 <- average_forecast(list(best_md1, best_md2, best_md3, best_md4,
                                           best_md5, best_md6, best_md7, best_md8, best_md9,
                                           best_md10, best_md11, best_md12, best_md13,
                                           best_md14, best_md15, best_md16, best_md17),
                                      data = fulldata, insample = insample, outsample = outsample,
                                      type = "recursive", 
                                      measures = c("MSE"),
                                      fweights = c("EW", "BICW", "MSFE", "DMSFE"),
                                      show_progress = FALSE)
  Recu_mfcb_ew <- append(Recu_mfcb_ew, best_nowmdxz_12[["avgforecast"]][1])
  Recu_mfcb_bicw <- append(Recu_mfcb_bicw, best_nowmdxz_12[["avgforecast"]][2])
  Recu_mfcb_msfew <- append(Recu_mfcb_msfew, best_nowmdxz_12[["avgforecast"]][3])
  Recu_mfcb_dmsfew <- append(Recu_mfcb_dmsfew, best_nowmdxz_12[["avgforecast"]][4])
}

# Saving predicting results
best_result_rec_for6 <- data.frame(mdate = mdate,
                                   real = tag,
                                   mf_ew = append(rep(0, leny), Recu_mf_ew),
                                   mf_bicw = append(rep(0, leny), Recu_mf_bicw),
                                   mf_msfew = append(rep(0, leny), Recu_mf_msfew),
                                   mf_dmsfew = append(rep(0, leny), Recu_mf_dmsfew),
                                   zp_ew = append(rep(0, leny), Recu_zp_ew),
                                   zp_bicw = append(rep(0, leny), Recu_zp_bicw),
                                   zp_msfew = append(rep(0, leny), Recu_zp_msfew),
                                   zp_dmsfew = append(rep(0, leny), Recu_zp_dmsfew),
                                   zz_ew = append(rep(0, leny), Recu_zz_ew),
                                   zz_bicw = append(rep(0, leny), Recu_zz_bicw),
                                   zz_msfew = append(rep(0, leny), Recu_zz_msfew),
                                   zz_dmsfew = append(rep(0, leny), Recu_zz_dmsfew),                                   
                                   cb_ew = append(rep(0, leny), Recu_cb_ew),
                                   cb_bicw = append(rep(0, leny), Recu_cb_bicw),
                                   cb_msfew = append(rep(0, leny), Recu_cb_msfew),
                                   cb_dmsfew = append(rep(0, leny), Recu_cb_dmsfew),
                                   mfcb_ew = append(rep(0, leny), Recu_mfcb_ew),
                                   mfcb_bicw = append(rep(0, leny), Recu_mfcb_bicw),
                                   mfcb_msfew = append(rep(0, leny), Recu_mfcb_msfew),
                                   mfcb_dmsfew = append(rep(0, leny), Recu_mfcb_dmsfew))


#-------------------------------------
# Model Combination (h=12)
#-------------------------------------
result_rec_for4 <- na.omit(read_xlsx(paste0(tmpdatpath, "Forecast_PPI_Realtime.xlsx"), sheet = 5, col_names = T))

# Select best individual MIDAS model
best_recu_df1 <- result_rec_for12[which(result_rec_for12$Recu_df1 == min(result_rec_for12$Recu_df1)), ]
best_recu_df2 <- result_rec_for12[which(result_rec_for12$Recu_df2 == min(result_rec_for12$Recu_df2)), ]
best_recu_df3 <- result_rec_for12[which(result_rec_for12$Recu_df3 == min(result_rec_for12$Recu_df3)), ]
best_recu_df4 <- result_rec_for12[which(result_rec_for12$Recu_df4 == min(result_rec_for12$Recu_df4)), ]
best_recu_df5 <- result_rec_for12[which(result_rec_for12$Recu_df5 == min(result_rec_for12$Recu_df5)), ]
best_recu_z_pos_gdp <- result_rec_for12[which(result_rec_for12$Recu_zpos_gdp == min(result_rec_for12$Recu_zpos_gdp)), ][1,]
best_recu_z_zero_gdp <- result_rec_for12[which(result_rec_for12$Recu_zzero_gdp == min(result_rec_for12$Recu_zzero_gdp)), ][1,]
best_recu_z_pos_cpi <- result_rec_for12[which(result_rec_for12$Recu_zpos_cpi == min(result_rec_for12$Recu_zpos_cpi)), ]
best_recu_z_zero_cpi <- result_rec_for12[which(result_rec_for12$Recu_zzero_cpi == min(result_rec_for12$Recu_zzero_cpi)), ]
best_recu_z_pos_iva <- result_rec_for12[which(result_rec_for12$Recu_zpos_iva == min(result_rec_for12$Recu_zpos_iva)), ]
best_recu_z_zero_iva <- result_rec_for12[which(result_rec_for12$Recu_zzero_iva == min(result_rec_for12$Recu_zzero_iva)), ]
best_recu_z_pos_fai <- result_rec_for12[which(result_rec_for12$Recu_zpos_fai == min(result_rec_for12$Recu_zpos_fai)), ]
best_recu_z_zero_fai <- result_rec_for12[which(result_rec_for12$Recu_zzero_fai == min(result_rec_for12$Recu_zzero_fai)), ]
best_recu_z_pos_m2 <- result_rec_for12[which(result_rec_for12$Recu_zpos_m2 == min(result_rec_for12$Recu_zpos_m2)), ]
best_recu_z_zero_m2 <- result_rec_for12[which(result_rec_for12$Recu_zzero_m2 == min(result_rec_for12$Recu_zzero_m2)), ]
best_recu_z_pos_ppi <- result_rec_for12[which(result_rec_for12$Recu_zpos_ppi == min(result_rec_for12$Recu_zpos_ppi)), ]
best_recu_z_zero_ppi <- result_rec_for12[which(result_rec_for12$Recu_zzero_ppi == min(result_rec_for12$Recu_zzero_ppi)), ]
## DF1
best_recu_df1_ylag <- 11 + best_recu_df1$Ylag
best_recu_df1_xlag <- 11 + best_recu_df1$Xlag
## DF2
best_recu_df2_ylag <- 11 + best_recu_df2$Ylag
best_recu_df2_xlag <- 11 + best_recu_df2$Xlag
## DF3
best_recu_df3_ylag <- 11 + best_recu_df3$Ylag
best_recu_df3_xlag <- 11 + best_recu_df3$Xlag
## DF4
best_recu_df4_ylag <- 11 + best_recu_df4$Ylag
best_recu_df4_xlag <- 11 + best_recu_df4$Xlag
## DF5
best_recu_df5_ylag <- 11 + best_recu_df5$Ylag
best_recu_df5_xlag <- 11 + best_recu_df5$Xlag
## z_pos_gdp
best_recu_zpos_gdp_ylag <- 11 + best_recu_z_pos_gdp$Ylag
best_recu_zpos_gdp_xlag <- 11 + best_recu_z_pos_gdp$Xlag
## z_zero_gdp
best_recu_zzero_gdp_ylag <- 11 + best_recu_z_zero_gdp$Ylag
best_recu_zzero_gdp_xlag <- 11 + best_recu_z_zero_gdp$Xlag
## z_pos_cpi
best_recu_zpos_cpi_ylag <- 11 + best_recu_z_pos_cpi$Ylag
best_recu_zpos_cpi_xlag <- 11 + best_recu_z_pos_cpi$Xlag
## z_zero_cpi
best_recu_zzero_cpi_ylag <- 11 + best_recu_z_zero_cpi$Ylag
best_recu_zzero_cpi_xlag <- 11 + best_recu_z_zero_cpi$Xlag
## z_pos_iva
best_recu_zpos_iva_ylag <- 11 + best_recu_z_pos_iva$Ylag
best_recu_zpos_iva_xlag <- 11 + best_recu_z_pos_iva$Xlag
## z_zero_iva
best_recu_zzero_iva_ylag <- 11 + best_recu_z_zero_iva$Ylag
best_recu_zzero_iva_xlag <- 11 + best_recu_z_zero_iva$Xlag
## z_pos_fai
best_recu_zpos_fai_ylag <- 11 + best_recu_z_pos_fai$Ylag
best_recu_zpos_fai_xlag <- 11 + best_recu_z_pos_fai$Xlag
## z_zero_fai
best_recu_zzero_fai_ylag <- 11 + best_recu_z_zero_fai$Ylag
best_recu_zzero_fai_xlag <- 11 + best_recu_z_zero_fai$Xlag
## z_pos_m2
best_recu_zpos_m2_ylag <- 11 + best_recu_z_pos_m2$Ylag
best_recu_zpos_m2_xlag <- 11 + best_recu_z_pos_m2$Xlag
## z_zero_m2
best_recu_zzero_m2_ylag <- 11 + best_recu_z_zero_m2$Ylag
best_recu_zzero_m2_xlag <- 11 + best_recu_z_zero_m2$Xlag
## z_pos_ppi
best_recu_zpos_ppi_ylag <- 11 + best_recu_z_pos_ppi$Ylag
best_recu_zpos_ppi_xlag <- 11 + best_recu_z_pos_ppi$Xlag
## z_zero_ppi
best_recu_zzero_ppi_ylag <- 11 + best_recu_z_zero_ppi$Ylag
best_recu_zzero_ppi_xlag <- 11 + best_recu_z_zero_ppi$Xlag

# Training and combine best individual models to get combination MIDAS
## MF
best_md1 <- midas_r(yy ~ mls(yy, 12:best_recu_df1_ylag, 1) + mls(zz1, 12:best_recu_df1_xlag, 1), start = NULL) # DF1
best_md2 <- midas_r(yy ~ mls(yy, 12:best_recu_df2_ylag, 1) + mls(zz2, 12:best_recu_df2_xlag, 1), start = NULL) # DF2
best_md3 <- midas_r(yy ~ mls(yy, 12:best_recu_df3_ylag, 1) + mls(zz3, 12:best_recu_df3_xlag, 1), start = NULL) # DF3
best_md4 <- midas_r(yy ~ mls(yy, 12:best_recu_df4_ylag, 1) + mls(zz4, 12:best_recu_df4_xlag, 1), start = NULL) # DF4
best_md5 <- midas_r(yy ~ mls(yy, 12:best_recu_df5_ylag, 1) + mls(zz5, 12:best_recu_df5_xlag, 1), start = NULL) # DF5
## MF+CB_2
best_md6 <- midas_r(yy ~ mls(yy, 12:best_recu_zpos_gdp_ylag, 1) + mls(xx1, 12:best_recu_zpos_gdp_xlag, 1), start = NULL) # z_pos_gdp
best_md7 <- midas_r(yy ~ mls(yy, 12:best_recu_zzero_gdp_ylag, 1) + mls(xx2, 12:best_recu_zzero_gdp_xlag, 1), start = NULL) # z_zero_gdp
best_md8 <- midas_r(yy ~ mls(yy, 12:best_recu_zpos_cpi_ylag, 1) + mls(xx3, 12:best_recu_zpos_cpi_xlag, 1), start = NULL) # z_pos_cpi
best_md9 <- midas_r(yy ~ mls(yy, 12:best_recu_zzero_cpi_ylag, 1) + mls(xx4, 12:best_recu_zzero_cpi_xlag, 1), start = NULL) # z_zero_cpi
best_md10 <- midas_r(yy ~ mls(yy, 12:best_recu_zpos_iva_ylag, 1) + mls(xx5, 12:best_recu_zpos_iva_xlag, 1), start = NULL) # z_pos_iva
best_md11 <- midas_r(yy ~ mls(yy, 12:best_recu_zzero_iva_ylag, 1) + mls(xx6, 12:best_recu_zzero_iva_xlag, 1), start = NULL) # z_zero_iva
best_md12 <- midas_r(yy ~ mls(yy, 12:best_recu_zpos_fai_ylag, 1) + mls(xx7, 12:best_recu_zpos_fai_xlag, 1), start = NULL) # z_pos_fai
best_md13 <- midas_r(yy ~ mls(yy, 12:best_recu_zzero_fai_ylag, 1) + mls(xx8, 12:best_recu_zzero_fai_xlag, 1), start = NULL) # z_zero_fai
best_md14 <- midas_r(yy ~ mls(yy, 12:best_recu_zpos_m2_ylag, 1) + mls(xx9, 12:best_recu_zpos_m2_xlag, 1), start = NULL) # z_pos_m2
best_md15 <- midas_r(yy ~ mls(yy, 12:best_recu_zzero_m2_ylag, 1) + mls(xx10, 12:best_recu_zzero_m2_xlag, 1), start = NULL) # z_zero_m2
best_md16 <- midas_r(yy ~ mls(yy, 12:best_recu_zpos_ppi_ylag, 1) + mls(xx11, 12:best_recu_zpos_ppi_xlag, 1), start = NULL) # z_pos_ppi
best_md17 <- midas_r(yy ~ mls(yy, 12:best_recu_zzero_ppi_ylag, 1) + mls(xx12, 12:best_recu_zzero_ppi_xlag, 1), start = NULL) # z_zero_ppi

# Estimating & Forecasting
Recu_mf_ew = c()
Recu_mf_bicw = c()
Recu_mf_msfew = c()
Recu_mf_dmsfew = c()
Recu_zp_ew = c()
Recu_zp_bicw = c()
Recu_zp_msfew = c()
Recu_zp_dmsfew = c()
Recu_zz_ew = c()
Recu_zz_bicw = c()
Recu_zz_msfew = c()
Recu_zz_dmsfew = c()
Recu_cb_ew = c()
Recu_cb_bicw = c()
Recu_cb_msfew = c()
Recu_cb_dmsfew = c()
Recu_mfcb_ew = c()
Recu_mfcb_bicw = c()
Recu_mfcb_msfew = c()
Recu_mfcb_dmsfew = c()

for (i in seq(1, (length(zp2)-1), 1)) {
  # Sampling
  i1 = i + 1
  m1 = i + leny - 1
  m2 = m1 + 1
  
  # Start & End
  start = c(2003,1)
  end = start + c(len, i-1)
  fulldata <- list(
    yy = ts(c(tag[1:m1], tag[m2]),start, end, 12),
    xx1 = ts(c(zp1[1:m1, i], zp1[m2, i1]), start, end, 12),
    xx2 = ts(c(zo1[1:m1, i], zo1[m2, i1]), start, end, 12),
    xx3 = ts(c(zp2[1:m1, i], zp2[m2, i1]), start, end, 12),
    xx4 = ts(c(zo2[1:m1, i], zo2[m2, i1]), start, end, 12),
    xx5 = ts(c(zp3[1:m1, i], zp3[m2, i1]), start, end, 12),
    xx6 = ts(c(zo3[1:m1, i], zo3[m2, i1]), start, end, 12),
    xx7 = ts(c(zp4[1:m1, i], zp4[m2, i1]), start, end, 12),
    xx8 = ts(c(zo4[1:m1, i], zo4[m2, i1]), start, end, 12),
    xx9 = ts(c(zp5[1:m1, i], zp5[m2, i1]), start, end, 12),
    xx10 = ts(c(zo5[1:m1, i], zo5[m2, i1]), start, end, 12),
    xx11 = ts(c(zp6[1:m1, i], zp6[m2, i1]), start, end, 12),
    xx12 = ts(c(zo6[1:m1, i], zo6[m2, i1]), start, end, 12),
    zz1 = ts(c(df1[1:m1, i], df1[m1, i]), start, end, 12),
    zz2 = ts(c(df2[1:m1, i], df2[m1, i]), start, end, 12),
    zz3 = ts(c(df3[1:m1, i], df3[m1, i]), start, end, 12),
    zz4 = ts(c(df4[1:m1, i], df4[m1, i]), start, end, 12),
    zz5 = ts(c(df5[1:m1, i], df5[m1, i]), start, end, 12)
  )
  insample <- 1:m1
  outsample <- (1:length(fulldata$yy))[-insample]
  
  # Combination Model
  ## MF
  best_nowmdzz <- average_forecast(list(best_md1, best_md2, best_md3, best_md4, best_md5),
                                   data = fulldata, insample = insample, outsample = outsample,
                                   type = "recursive", 
                                   measures = c("MSE"),
                                   fweights = c("EW", "BICW", "MSFE", "DMSFE"),
                                   show_progress = FALSE)
  Recu_mf_ew <- append(Recu_mf_ew, best_nowmdzz[["avgforecast"]][1])
  Recu_mf_bicw <- append(Recu_mf_bicw, best_nowmdzz[["avgforecast"]][2])
  Recu_mf_msfew <- append(Recu_mf_msfew, best_nowmdzz[["avgforecast"]][3])
  Recu_mf_dmsfew <- append(Recu_mf_dmsfew, best_nowmdzz[["avgforecast"]][4])
  ## MF+CB_+
  best_nowmdxzp <- average_forecast(list(best_md1, best_md2, best_md3, best_md4, best_md5, best_md12),
                                    data = fulldata, insample = insample, outsample = outsample,
                                    type = "recursive", 
                                    measures = c("MSE"),
                                    fweights = c("EW", "BICW", "MSFE", "DMSFE"),
                                    show_progress = FALSE)
  Recu_zp_ew <- append(Recu_zp_ew, best_nowmdxzp[["avgforecast"]][1])
  Recu_zp_bicw <- append(Recu_zp_bicw, best_nowmdxzp[["avgforecast"]][2])
  Recu_zp_msfew <- append(Recu_zp_msfew, best_nowmdxzp[["avgforecast"]][3])
  Recu_zp_dmsfew <- append(Recu_zp_dmsfew, best_nowmdxzp[["avgforecast"]][4])  
  ## MF+CB_0
  best_nowmdxzz <- average_forecast(list(best_md1, best_md2, best_md3, best_md4, best_md5, best_md13),
                                    data = fulldata, insample = insample, outsample = outsample,
                                    type = "recursive", 
                                    measures = c("MSE"),
                                    fweights = c("EW", "BICW", "MSFE", "DMSFE"),
                                    show_progress = FALSE)
  Recu_zz_ew <- append(Recu_zz_ew, best_nowmdxzz[["avgforecast"]][1])
  Recu_zz_bicw <- append(Recu_zz_bicw, best_nowmdxzz[["avgforecast"]][2])
  Recu_zz_msfew <- append(Recu_zz_msfew, best_nowmdxzz[["avgforecast"]][3])
  Recu_zz_dmsfew <- append(Recu_zz_dmsfew, best_nowmdxzz[["avgforecast"]][4])         
  ## MF+CB_2
  best_nowmdxz_2 <- average_forecast(list(best_md1, best_md2, best_md3, best_md4,                                                                best_md5, best_md12, best_md13),
                                     data = fulldata, insample = insample, outsample = outsample,
                                     type = "recursive", 
                                     measures = c("MSE"),
                                     fweights = c("EW", "BICW", "MSFE", "DMSFE"),
                                     show_progress = FALSE)
  Recu_cb_ew <- append(Recu_cb_ew, best_nowmdxz_2[["avgforecast"]][1])
  Recu_cb_bicw <- append(Recu_cb_bicw, best_nowmdxz_2[["avgforecast"]][2])
  Recu_cb_msfew <- append(Recu_cb_msfew, best_nowmdxz_2[["avgforecast"]][3])
  Recu_cb_dmsfew <- append(Recu_cb_dmsfew, best_nowmdxz_2[["avgforecast"]][4])  
  ## MF+CB_12
  best_nowmdxz_12 <- average_forecast(list(best_md1, best_md2, best_md3, best_md4,
                                           best_md5, best_md6, best_md7, best_md8, best_md9,
                                           best_md10, best_md11, best_md12, best_md13,
                                           best_md14, best_md15, best_md16, best_md17),
                                      data = fulldata, insample = insample, outsample = outsample,
                                      type = "recursive", 
                                      measures = c("MSE"),
                                      fweights = c("EW", "BICW", "MSFE", "DMSFE"),
                                      show_progress = FALSE)
  Recu_mfcb_ew <- append(Recu_mfcb_ew, best_nowmdxz_12[["avgforecast"]][1])
  Recu_mfcb_bicw <- append(Recu_mfcb_bicw, best_nowmdxz_12[["avgforecast"]][2])
  Recu_mfcb_msfew <- append(Recu_mfcb_msfew, best_nowmdxz_12[["avgforecast"]][3])
  Recu_mfcb_dmsfew <- append(Recu_mfcb_dmsfew, best_nowmdxz_12[["avgforecast"]][4])
}

# Saving predicting results
best_result_rec_for12 <- data.frame(mdate = mdate,
                                    real = tag,
                                    mf_ew = append(rep(0, leny), Recu_mf_ew),
                                    mf_bicw = append(rep(0, leny), Recu_mf_bicw),
                                    mf_msfew = append(rep(0, leny), Recu_mf_msfew),
                                    mf_dmsfew = append(rep(0, leny), Recu_mf_dmsfew),
                                    zp_ew = append(rep(0, leny), Recu_zp_ew),
                                    zp_bicw = append(rep(0, leny), Recu_zp_bicw),
                                    zp_msfew = append(rep(0, leny), Recu_zp_msfew),
                                    zp_dmsfew = append(rep(0, leny), Recu_zp_dmsfew),
                                    zz_ew = append(rep(0, leny), Recu_zz_ew),
                                    zz_bicw = append(rep(0, leny), Recu_zz_bicw),
                                    zz_msfew = append(rep(0, leny), Recu_zz_msfew),
                                    zz_dmsfew = append(rep(0, leny), Recu_zz_dmsfew),                                     
                                    cb_ew = append(rep(0, leny), Recu_cb_ew),
                                    cb_bicw = append(rep(0, leny), Recu_cb_bicw),
                                    cb_msfew = append(rep(0, leny), Recu_cb_msfew),
                                    cb_dmsfew = append(rep(0, leny), Recu_cb_dmsfew),
                                    mfcb_ew = append(rep(0, leny), Recu_mfcb_ew),
                                    mfcb_bicw = append(rep(0, leny), Recu_mfcb_bicw),
                                    mfcb_msfew = append(rep(0, leny), Recu_mfcb_msfew),
                                    mfcb_dmsfew = append(rep(0, leny), Recu_mfcb_dmsfew))


#-------------------------------------
# Save Results
#-------------------------------------
sheets <- list("sheet1" = best_result_rec_now,
               "sheet2" = best_result_rec_for1,
               "sheet3" = best_result_rec_for3,
               "sheet4" = best_result_rec_for6,
               "sheet5" = best_result_rec_for12)

write.xlsx(sheets, paste0(path, "/Data/Forecast_PPI_Best_Realtime.xlsx"), 
           sheetName=c("rec_now", "rec_for1", "rec_for3", "rec_for6", "rec_for12"))