################################################################################ # R code for Empirical Application (Section 6 of paper) # File: New-test-Combining-2022.r [18-07-2022] # Coded by JDG # # INPUT FILE: SP195101_201612_30cols.csv (792 * 30) matrix # [effective use: 792-1=791 obs., 1-step ahead forecast] # NO TECHNICAL VARIABLES # # Reference: # Jan G. De Gooijer (2022) # "Penalized Averaging of Quantile Forecasts from GARCH Models with # Many Exogenous Predictors", # Computational Economics. DOI: https://doi.org/10.1007/s10614-022-10289-9 # ################################################################################# # rm(list = ls()) # empty memory library(MASS) # Functions and datasets to support Venables and Ripley, # ``Modern Applied Statistics with S'' (4th edition, 2002). library(zoo) # Infrastructure for Regular and Irregular Time Series library(dynlm) # for dynamic linear regression models library(rugarch) # for modeling univariate (G)ARCH-type models library(quantreg) # Quantile Regression library(rqPen) # Penalized quantile regression for LASSO, SCAD, and MCP pfs library(mboost) # for gamboost, Model-Based Boosting # library(FHDQR) # for ADMM penenalized adaptive LASSO library(cqrReg) # Quantile Regression (qr) with Adaptive Lasso Penalty (lasso) library(sandwich) # used for Newey-West HAC cov matrix estimation, when H>1 library(xts) # eXtensible Time Series library(stats) library(matrixStats) # for colSds(X) each column Stand.Dev of a matrix library(forecast) # for tail and dm.test ########################################### ######## Supporting functions ############ # quantile loss function Q.loss <- function(u, tt) u * ( tt - (u < 0)* 1) ########################################### Set.lam <- function(x, y, tt) { # We use this function to help us # to select the optimal lambda values. # Penalty for Q.BIC is log(error) + sm * (log(n)/ 2n) * Cn # In order of severity: # 1: where Cn = 0.5 --- curiosity # 2: where Cn = 1 --- the usual BIC -- meant for d << n # 3: where Cn = log(d) * log(log(n)) / log(n) # 4: where Cn = log(d)/2 # 5: where Cn = log(d) # # nfolds = 10 CV (default) ################################################################ pd <- dim(x)[2] n <- dim(x)[1] options(warn=-1) bic_model <- cv.rq.pen(x,y,tau=tt,penalty="LASSO",nlambda=60,criteria="BIC") lm.vec <- bic_model$cv[,1] kk <- length(lm.vec) # QBIC.c.1 <- rep(0,kk) QBIC.c.3 <- rep(0,kk) # QBIC.c.4 <- rep(0,kk) # QBIC.c.5 <- rep(0,kk) for (j in 1:kk) { fit_model <- rq.lasso.fit(x,y,tau=tt,lambda=lm.vec[j]) # ,intercept=FALSE) # rq.nc.fit = Non-convex penalized quantile regression s.m <- length(which(coefficients(fit_model)[2:(pd+1)] != 0)) uu <- residuals(fit_model) # penl.c.1 <- s.m * log(pd) / (4 * n) # Cn = 0.5 -- drastic penl.c.3 <- s.m * log(pd) * log(log(n)) / (2 * n) # Sherwood # and Wang (2016), The Annals of Statistics,44, 288-317. # penl.c.4 <- s.m * log(pd) * log(n) / (4 * n) # penl.c.5 <- s.m * log(pd) * log(n) / (2 * n) # stringent # QBIC.c.1[j] <- log(sum(Q.loss(uu,tt))) + penl.c.1 QBIC.c.3[j] <- log(sum(Q.loss(uu,tt))) + penl.c.3 # QBIC.c.4[j] <- log(sum(Q.loss(uu,tt))) + penl.c.4 # QBIC.c.5[j] <- log(sum(Q.loss(uu,tt))) + penl.c.5 } # pp2 <- which.min(bic_model$cv[,2]) # pp1 <- which.min(QBIC.c.1) pp3 <- which.min(QBIC.c.3) # pp4 <- which.min(QBIC.c.4) # pp5 <- which.min(QBIC.c.5) # Optimal lambda vector based on each penalty coefficient Cn out <- lm.vec[c(pp3)] # [c(pp1,pp2,pp3,pp4,pp5)] return(out) } lagmatrix <- function(x,max.lag) embed(c(rep(NA,max.lag), x), max.lag+1) ################################################################ ## Set.ENC(Y.est.np,EW.forecast[,1],Y.M10,alpha,N.total) Set.ENC <- function(Y.est.np,Y1,Y2,alpha,N.total){ ### Encompassing (ENC) Wald type test statistic ############### ### Reference: Fuertes and Olmo (IJF 29, 28-42, 2013) ######### ################################################################ m <- 3 # no. lambda's in ENC Wald type test statistic, so 3 p-values N <- 300 fit1 <- rq(Y.est.np[N:(N+N.total-1)] ~ Y1 + Y2, tau = alpha) thetaopt <- fit1$coefficients h <- N.total^(-1/4) # choice of bandwidth parameter q.mat <- cbind(rep(1,N.total), Y1,Y2) VaR <- q.mat %*% thetaopt # 100 values e <- Y.est.np[N:(N+N.total-1)]-VaR[1:(100-5),1] # 300 - 399 = 100 values e = e[!is.na(e)] # remove NaNs phi = t(alpha-(e<0)) %*% e ff = mean(phi) Omega0 <- matrix(0,m,m) Omega1 <- matrix(0,m,m) # initializing variance matrices for (t in 1:length(e)){ Omega01 <- (q.mat[t, ]) %*% t(q.mat[t,]) Omega0 <- Omega0+ Omega01/length(e) Omega1 <- Omega1+Omega01*(abs(e[t])crit) pvalue <- pchisq(ENC, df=m, lower.tail=FALSE) return(list(ENC=ENC,pvalue=pvalue,ff=ff,Sigma=Sigmadia,y=VaR)) } ####################################### #' Test for Forecast Encompassing #' #' The Test for Forecast Encompassing compares the forecasts #' of two forecast methods. #' #' The null hypothesis is that the method A forecast encompasses #' of method B, i.e., method B have the same information #' of method A. The alternative hypothesis is that method B #' have additional information. #' #' @param fA Forecast from method A. #' @param fB Forecast from method B. #' @param y The observed time series. #' @param h Length ahead of the forecast. #' #' @return An object of class "coeftest" which is #' essentially a coefficient matrix with columns #' containing the estimates, associated standard #' errors, test statistics and p values. #' #' @import lmtest sandwich stats #' @export comp.test <- function(y, fA, fB, h){ yh <- y - glag(y, -h) fA <- fA - glag(y, -h) fB <- fB - glag(y, -h) fit <- stats::lm(yh ~ fA + fB) return(lmtest::coeftest(fit, vcov. = sandwich::NeweyWest(fit))) } #' Generalized Lag #' #' Compute a lagged version of a vector like time series #' #' @param x A vector or univariate time series #' @param k The number of lags #' #' @return A vector object with the same length of x. #' #' @import stats #' #' @export #' glag <- function(x, k = -1){ h <- abs(k) if(k>0){ xh <- c(x, rep(NA, h)) ans <- embed(xh, h+1)[,1] } if(k<0){ xh <- c(rep(NA, h), x) ans <- embed(xh, h+1)[,(h+1)] } if(k==0){ ans <- x } return(ans) } ###################################################### ########### The main program ######################### library(readr) H1 <- 6 # Forecast horizon H <- H1 # Forecast horizon ## CHANGE drives !!! Returns <- read_csv("C:/Risk-Premium/Neely-Rapach-Risk-Premium-Matl/SP195101_201612_30cols.csv", col_names = FALSE) Rf.sp <- read_csv("C:/Risk-Premium/Analysis-2018/Economic-profits/Risk_free.csv", col_names = FALSE) ## CHANGE DRIVES !!!! ## Returns <- read_csv("D:/Simul-Hybrid-paper-tau08-corr05/SP195101_201612_30cols.csv", ## col_names = FALSE) ## Rf.sp <- read_csv("D:/Simul-Hybrid-paper-tau08-corr05/Risk_free.csv", ## col_names = FALSE) # Date-C1 EquityPremium-C2 DP-C3 DY-C4 EP-C5 DE-C6 RVOL-C7 BM-C8 # NTIS-C9 TBL-C10 LTY-C11 # LTR-C12 TMS-C13 DFY-C14 DFR-C15 INFL-C16 MA(1,9)-C17 MA(1,12)-C18 # MA(2,9)-C19 MA(2,12)-C20 MA(3,9)-C21 # MA(3,12)-C22 MOM(9)-C23 MOM(12)-C24 VOL(1,9)-C25 VOL(1,12)-C26 # VOL(2,9)-C27 VOL(2,12)-C28 VOL(3,9)-C29 # VOL(3,12)-C30 # 30 Vars = C1 (Date) + C2 (Eq Prem) + (C3 - C16) # (Macroeconomics) + (C17 - C30) (14 Technical variables, not used in analysis) sp <- as.zoo(ts(Returns[,2],start=c(1951,1),frequency=12)) # S&P500 returns Y <- as.ts(sp) # H=1 792 obs., 1951:1 - 2016:12 S&P500 Rf.s <- as.zoo(ts(Rf.sp[,2],start=c(1951,1),frequency=12)) # Risk free Rf <- as.ts(Rf.s) # risk free Xdata <- as.ts(Returns[,3:30],frequency=1) # [792 * 28] exogenous covariates dates <- seq(as.Date("1951-01-01"), length = 792, by = "months") inputs.Y <- xts(x = sp, order.by = dates) # NOT standardized ncol <- ncol(Xdata)-14 # NO 14 TECHNICAL VARIABLES nrow <- nrow(Xdata) N <- 168 # length first sample # N <- 792 # CHANGE !! ## QQ <- c(0.05, 0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9, 0.95) # quantiles considered tt <- 0.1 # QQ[1] # CHANGE FOR DIFFERENT QUANTILE LEVELS !! ## nr.q <- length(QQ) n <- length(Y) ########################################################## Y.total <- window(Y, start=c(1951,1),end=c(1964,12)) # 168 obser. # Rf.total <- window(Rf,start=c(1951,1),end=c(1964,12)) # 168 obser. Est.For <- window(Y,start=c(1964,12),end=c(2016,11)) # 624 obser. N.total <- length(Est.For) ################### number of recursions ################################# ## N.total=624 # N=791(last; 1951:1-2016:11)-N=168(first)+1=no.recursions ########################################################################### P.err.1 <- matrix(0,N.total,ncol) # PEs,[ 4(G)ARCH models], ncol=14 predictors P.err.2 <- matrix(0,N.total,ncol) P.err.3 <- matrix(0,N.total,ncol) P.err.4 <- matrix(0,N.total,ncol) PE.EW <- matrix(0,N.total,4) # PEs, old benchmark EW.forecast <- matrix(0,N.total,4) # Equal weighting forecasts # (mean of cond.quantiles for each of the 4 Parametric models) P.err.min.all <- matrix(0,N.total,4) # minimum 4 Parametric GARCH models P.err.loc.all <- matrix(0,N.total,4) # location of minimum %% NOTE: Method 5 = Method 2 in paper, Method 6 = Method 3 in paper, etc. ### Method 5 (PA-L): penalized averaging of all parametric models P.err.Meth.5.1 <- matrix(0,N.total,1) # PE linear fit, LASSO, Method 5-1 P.err.Meth.5.2 <- matrix(0,N.total,1) # PE linear fit, LASSO, Method 5-2 P.err.Meth.5.3 <- matrix(0,N.total,1) # PE linear fit, LASSO, Method 5-3 P.err.Meth.5.4 <- matrix(0,N.total,1) # PE linear fit, LASSO, Method 5-4 S.Meth.5.1 <- matrix(0,N.total,2*ncol) S.Meth.5.2 <- matrix(0,N.total,2*ncol) S.Meth.5.3 <- matrix(0,N.total,2*ncol) S.Meth.5.4 <- matrix(0,N.total,2*ncol) ### Method 6 (PA-AL): penalized averaging of all NP marginal quantiles P.err.Meth.6 <- matrix(0,N.total,1) S.Meth.6 <- matrix(0,N.total,2*ncol) S.sel.6 <- matrix(0,N.total,2*ncol) ### Method 7 (P-Lin-QR): penalized linear model P.err.Meth.7 <- matrix(0,N.total,1) S.Meth.7 <- matrix(0,N.total,ncol) ### Method 8 (Ad-QR): like Method (6) but uses selected variables P.err.Meth.8 <- matrix(0,N.total,1) S.Meth.8 <- S.Meth.6 ### Method 9: Hybrid, NOT included in paper QQ.Meth.9.1 <- matrix(0,N.total,2*ncol) QQ.Meth.9.2 <- matrix(0,N.total,2*ncol) QQ.Meth.9.3 <- matrix(0,N.total,2*ncol) QQ.Meth.9.4 <- matrix(0,N.total,2*ncol) S.Meth.9.1 <- matrix(0,N.total,2*ncol) S.Meth.9.2 <- matrix(0,N.total,2*ncol) S.Meth.9.3 <- matrix(0,N.total,2*ncol) S.Meth.9.4 <- matrix(0,N.total,2*ncol) P.err.Meth.9.1 <- matrix(0,N.total) P.err.Meth.9.2 <- matrix(0,N.total) P.err.Meth.9.3 <- matrix(0,N.total) P.err.Meth.9.4 <- matrix(0,N.total) #### Method 10 (H-aL): Hybrid with adaptive-LASSO P.err.Meth.10.1 <- matrix(0,N.total) P.err.Meth.10.2 <- matrix(0,N.total) P.err.Meth.10.3 <- matrix(0,N.total) P.err.Meth.10.4 <- matrix(0,N.total) S.Meth.10.1 <- matrix(0,N.total,2*ncol) S.Meth.10.2 <- matrix(0,N.total,2*ncol) S.Meth.10.3 <- matrix(0,N.total,2*ncol) S.Meth.10.4 <- matrix(0,N.total,2*ncol) ######################################## ######################################## W.l.lasso <- matrix(0,N.total,1) W.np.lasso <- matrix(0,N.total,1) EW.1 <- matrix(0,N.total,1) EW.2 <- matrix(0,N.total,1) EW.3 <- matrix(0,N.total,1) EW.4 <- matrix(0,N.total,1) R.l.lasso <- matrix(0,N.total,1) R.np.lasso <- matrix(0,N.total,1) Y.total.series <- ts(Y,freq=1) S.l.ad.sel <- array(0,dim=c(N.total,140)) # Selected vars, linear, additive LASSO S.np.ad.sel <- array(0,dim=c(N.total,140)) # Selected vars, nonp., additive LASSO fit.l.lasso = array(0,dim=c(N.total,141)) fit.np.lasso = array(0,dim=c(N.total,141)) ################################################################### # Start 1st BIG loop for (u in 1:N.total){ # 478:N.total){ # N.total = 624; # start BIG loop for estimation & forecasting # u=1 nH <- (N+u-1)-H+1 # nH <- (N+u-1)-H+1 # H = forecast horizon, here H=1 Y.H <- Y.total.series[1:(N+u-H+1)] # NOT standardized # u=1, Y.H[1:N+1)=[1:169]; u=2, Y.H[1:(N+2)];... # u=624=N.total, Y.H[1:(N+N.total-1+1)=792] X.H <- Xdata[1:(N+u-H+1),] # u=1 [(N+1) * ncol.X] = [169 * 28] [169 * 14] # u=N.total=624 [(N+N.total-1+1) * ncol] = [792 * 28] [792 * 14] N.u <- (N+u-1)-H+1 cat("u =", u,"; (N+u-1)-H+1 =", (N+u-1)-H+1,"; nH =" , nH,"\n\n") # Y.est <- as.ts(Y.H[1:nH]) # estimation part, increases for each u ## New part: ### Y.est <- as.ts(Y.H[1:(N+u-H+1)])*100 # in percentages X.H <- Xdata[1:(N+u-H+1),] ## Y.est <- as.ts(Y[1:168])*100 # Y[1:732]=1951:1 - 2011:12 in Neely in % ## X.H <- Xdata[1:168,] # ,,, PX <- matrix(0,(N+u-H+1),14) # standardized macro-economic vars. p.mean <- colMeans(X.H) p.sd <- colSds(X.H) for (i in 1:14){ PX[,i] <- t((X.H[,i]-p.mean[i])/p.sd[i]) } ## NO TECHNICAL VARIABLES # X.H.jit <- jitter(X.H[,15:28]) # jittering technical indicators mean.Y <- mean(Y.est) sd.Y <- sd(Y.est) m.Y <- mean(Y.est/100) s.Y <- sd(Y.est/100) Y.est <- (Y.est-mean.Y)/sd.Y # standardized Y=S&P500 (Y.est=Y.H*100) # X.H <- cbind(PX,X.H.jit) # stand. macro + [0-1] techn. vars X.H <- as.ts(X.H) # NO TECHNICAL VARIABLES (28 macro vars only) X.H <- unname(X.H) # Remove the names or dimnames attribute of an R object. ## STAGES 1 and 2: ####### PARAMETRIC conditional marginal fits from linear (G)ARCH ######## ####### Model.1 [=Model (2) of Cenesizoglu & Timmermann, ####### J. Banking & Finance, 2012] qf.1 <- array(0,dim=c((N+u-1-H+1),ncol)) dfm.coef.1 <- matrix(0,ncol,1) for (i in 1:ncol){ x.cov <- as.ts(X.H[,i]) # Note: [1 - 14 (Macro), 15-28 (Technical) dfm <- dynlm(Y.est ~ L(x.cov, 1)) # Y.est, L=lag = 1 dfm.coef.1[i] <- dfm$coefficients[2] ## Conditional quantiles for dynamic linear regression with normal CDF qf.1.e <- dfm$fitted.values+sd(dfm$residuals)*rep(1,length(dfm$residuals))*qnorm(tt) qf.1[,i] <- qf.1.e # one-step ahead conditional quantiles P.err.1[u,i] <- Y.est[N+u-H+1]-qf.1[(N+u-1-H+1),i] } # end i =1:nol=28 #### Model.2 [=Model (4) of Cenesizoglu & Timmermann, #### J. Banking & Finance, 2012] qf.2 <- array(0,dim=c((N+u-1-H+1),ncol)) for (i in 1:ncol){ inputs.X <- xts(x = Xdata[,i],order.by=dates) inputs.X.std <- (inputs.X[1:(N+u-H+1)]-mean(inputs.X[1:(N+u-H+1)]))/sd(inputs.X[1:(N+u-H+1)]) inputs.Y.std <- (inputs.Y[1:(N+u-H+1)]-mean(inputs.Y[1:(N+u-H+1)]))/sd(inputs.Y[1:(N+u-H+1)]) inputs <- cbind(inputs.Y.std,inputs.X.std) # inputs consists of 2 variables model.2.e <- ugarchspec( variance.model = list(model="eGARCH",garchOrder=c(1,1),external.regressors=inputs[1:(N+u-1-H+1),2]), mean.model = list(armaOrder=c(0,0),include.mean=TRUE), distribution.model="norm") model.2.e.fit <- ugarchfit(spec=model.2.e,data=inputs[2:(N+u-H+1),1],solver.control=list()) qf.2.e <- as.matrix(quantile(model.2.e.fit,tt)) # one-step ahead cond. quantiles qf.2[,i] <- qf.2.e P.err.2[u,i] <- Y.est[N+u-H+1]-qf.2[(N+u-1-H+1),i] } ##### Model.3 [=Model (5) of Cenesizoglu & Timmermann, J. Banking & Finance, 2012] qf.3 <- array(0,dim=c((N+u-1-H+1),ncol)) for (i in 1:ncol){ inputs.X <- xts(x = Xdata[,i],order.by=dates) inputs.X.std <- (inputs.X[1:(N+u-H+1)]-mean(inputs.X[1:(N+u-H+1)]))/sd(inputs.X[1:(N+u-H+1)]) inputs.Y.std <- (inputs.Y[1:(N+u-H+1)]-mean(inputs.Y[1:(N+u-H+1)]))/sd(inputs.Y[1:(N+u-H+1)]) inputs <- cbind(inputs.Y.std,inputs.X.std) # inputs consists of 2 variables model.3.e <- ugarchspec( variance.model = list(model="eGARCH",garchOrder=c(1,1)), mean.model = list(armaOrder=c(0,0),external.regressors=inputs[1:(N+u-1-H+1),2], include.mean = TRUE),distribution.model="norm") model.3.e.fit <- ugarchfit(spec=model.3.e,data=inputs[2:(N+u-H+1),1],solver.control=list()) qf.3.e <- as.matrix(quantile(model.3.e.fit,tt)) # one-step ahead cond. quantiles qf.3[,i] <- qf.3.e P.err.3[u,i] <- Y.est[N+u-H+1]-qf.3[(N+u-1-H+1),i] } ##### Model.4 [=Model (6) of Cenesizoglu & Timmermann, J. Banking & Finance, 2012] qf.4 <- array(0,dim=c((N+u-1-H+1),ncol)) for (i in 1:ncol){ inputs.X <- xts(x = Xdata[,i],order.by=dates) inputs.X.std <- (inputs.X[1:(N+u-H+1)]-mean(inputs.X[1:(N+u-H+1)]))/sd(inputs.X[1:(N+u-H+1)]) inputs.Y.std <- (inputs.Y[1:(N+u-H+1)]-mean(inputs.Y[1:(N+u-H+1)]))/sd(inputs.Y[1:(N+u-H+1)]) inputs <- cbind(inputs.Y.std,inputs.X.std) # inputs consists of 2 variables model.4.e <- ugarchspec(variance.model=list(model='eGARCH',garchOrder=c(1,1), external.regressors = inputs[1:(N+u-H+1),2],variance.targeting=T), mean.model = list(armaOrder=c(0,0),external.regressors = inputs[1:(N+u-1-H+1),2], include.mean=TRUE),distribution.model="norm") model.4.e.fit <- ugarchfit(spec=model.4.e,data=inputs[2:(N+u-H+1),1],solver.control=list()) qf.4.e <- as.matrix(quantile(model.4.e.fit,tt)) # one-step ahead cond. quantiles qf.4[,i] <- qf.4.e P.err.4[u,i] <- Y.est[N+u-H+1]-qf.4[(N+u-1-H+1),i] } #### Evaluation: Models 1 - 4, Prediction errors #### #### 4 NEW benchmark models #### P.err.min.all[u,1] = min(P.err.1[u, ]) # minimum value PE P.err.loc.all[u,1] = which.min(P.err.1[u, ]) # location predictor P.err.min.all[u,2] = min(P.err.2[u, ]) P.err.loc.all[u,2] = which.min(P.err.2[u, ]) P.err.min.all[u,3] = min(P.err.3[u, ]) P.err.loc.all[u,3] = which.min(P.err.3[u, ]) P.err.min.all[u,4] = min(P.err.4[u, ]) P.err.loc.all[u,4] = which.min(P.err.4[u, ]) ### EW conditional quantile forecast across 28 variables: "OLD" benchmark ## EW.forecast[u,1] = mean(qf.1[(N+u-1-H+1),]) EW.forecast[u,2] = mean(qf.2[(N+u-1-H+1),]) EW.forecast[u,3] = mean(qf.3[(N+u-1-H+1),]) EW.forecast[u,4] = mean(qf.4[(N+u-1-H+1),]) PE.EW[u,1] <- Y.est[N+u-H+1]-EW.forecast[u,1] PE.EW[u,2] <- Y.est[N+u-H+1]-EW.forecast[u,2] PE.EW[u,3] <- Y.est[N+u-H+1]-EW.forecast[u,3] PE.EW[u,4] <- Y.est[N+u-H+1]-EW.forecast[u,4] } # END FIRST BIG LOOP, u=[1:N.total=624], for ONE quantile level # EXECUTE THE ABOVE PART FIRST ################################################################# ################################################################# ###### SUMMARY MEASURES 1st four methods: skip THIS PART ###### ###### 4 Frequency Tables with selected benchmark models ################################################################# library(plyr) # Tools for Splitting, Applying and Combining Data # Model 1 new benchmark (best selected models) S.1.c <- c(P.err.loc.all[ ,1]) # convert matrix in column vector S.1.nz <- S.1.c[S.1.c !=0] # remove zeros x.S1 <- S.1.nz # frequency table cat("Model 1 new benchmark:") cat("Series ID ", "Freq.", " Cumulative", " Rel. Freq =100*Cum/total(Cum)") cbind(Freq=table(x.S1),Cumul=cumsum(table(x.S1)), Relative=prop.table(table(x.S1))*100) # Model 2 new benchmark (best selected models) S.2.c <- c(P.err.loc.all[ ,2]) # convert matrix in column vector S.2.nz <- S.2.c[S.2.c !=0] # remove zeros x.S2 <- S.2.nz # frequency table # tt.id <- as.matrix(count(x.S2)$x) # IDs of series cat("Model 2 new benchmark:") cat("Series ID ", "Freq.", " Cumulative", " Rel. Freq =100*Cum/total(Cum)") cbind(Freq=table(x.S2),Cumul=cumsum(table(x.S2)), Relative=prop.table(table(x.S2))*100) # Model 3 new benchmark (best selected models) S.3.c <- c(P.err.loc.all[ ,3]) # convert matrix in column vector S.3.nz <- S.3.c[S.3.c !=0] # remove zeros x.S3 <- S.3.nz # frequency table # tt.id <- as.matrix(count(x.S3)$x) # IDs of series cat("Model 3 new benchmark:") cat("Series ID ", "Freq.", " Cumulative", " Rel. Freq =100*Cum/total(Cum)") cbind(Freq=table(x.S3),Cumul=cumsum(table(x.S3)), Relative=prop.table(table(x.S3))*100) # Model 4 new benchmark (best selected models) S.4.c <- c(P.err.loc.all[ ,4]) # convert matrix in column vector S.4.nz <- S.4.c[S.4.c !=0] # remove zeros x.S4 <- S.4.nz # frequency table # tt.id <- as.matrix(count(x.S4)$x) # IDs of series cat("Model 4 new benchmark:") cat("Series ID ", "Freq.", " Cumulative", " Rel. Freq =100*Cum/total(Cum)") cbind(Freq=table(x.S4),Cumul=cumsum(table(x.S4)), Relative=prop.table(table(x.S4))*100) ############################################# ########### start WITH NEXT PART ############ ## NEW BENCHMARK MODELS Ti.1 = mean((tt-(P.err.min.all[ ,1]<0))*P.err.min.all[ ,1]) Ti.2 = mean((tt-(P.err.min.all[ ,2]<0))*P.err.min.all[ ,2]) Ti.3 = mean((tt-(P.err.min.all[ ,3]<0))*P.err.min.all[ ,3]) Ti.4 = mean((tt-(P.err.min.all[ ,4]<0))*P.err.min.all[ ,4]) ### "OLD EW " Ti.naive.1 = mean((tt-(PE.EW[ ,1]<0))*PE.EW[ ,1]) Ti.naive.2 = mean((tt-(PE.EW[ ,2]<0))*PE.EW[ ,2]) Ti.naive.3 = mean((tt-(PE.EW[ ,3]<0))*PE.EW[ ,3]) Ti.naive.4 = mean((tt-(PE.EW[ ,4]<0))*PE.EW[ ,4]) ## Ratios of Tick functions; "OLD=naive" models M1, ... , M4 ## vs. new benchmark models # M1 R1.1 = Ti.naive.1/Ti.1 # >1 = Ti.naive is WORSE than Ti.1 R1.2 = Ti.naive.1/Ti.2 R1.3 = Ti.naive.1/Ti.3 R1.4 = Ti.naive.1/Ti.4 Ratio.naive1.new = cbind(R1.1,R1.2,R1.3,R1.4) # 1st index = naive, 2nd "new" ##### Ratios of Tick functions vs naive models M1 ... M4 # M2 R2.1 = Ti.naive.2/Ti.1 # >1 = Ti.naive is worse than Ti.1 R2.2 = Ti.naive.2/Ti.2 R2.3 = Ti.naive.2/Ti.3 R2.4 = Ti.naive.2/Ti.4 Ratio.naive2.new = cbind(R2.1,R2.2,R2.3,R2.4) ## Ratios of Tick functions; "OLD=naive" models M1 ... M4 ## vs. new benchmark models # M3 R3.1 = Ti.naive.3/Ti.1 # >1 = Ti.naive is worse than Ti.1 R3.2 = Ti.naive.3/Ti.2 R3.3 = Ti.naive.3/Ti.3 R3.4 = Ti.naive.3/Ti.4 Ratio.naive3.new = cbind(R3.1,R3.2,R3.3,R3.4) # 1st index = naive, # 2nd "new" ##### Ratios of Tick functions vs naive models M1, ... ,M4 # M4 R4.1 = Ti.naive.4/Ti.1 # >1 = Ti.naive.4 is worse than Ti.1 R4.2 = Ti.naive.4/Ti.2 R4.3 = Ti.naive.4/Ti.3 R4.4 = Ti.naive.4/Ti.4 Ratio.naive4.new = cbind(R4.1,R4.2,R4.3,R4.4) ############################################## ############# Methods 5 - 10 ################# for (u in 1:N.total){ # u = 1 cat("u =", u,"; (N+u-1)-H+1 =", (N+u-1)-H+1,"\n\n") ## Use 4 parametric (P) new benchmark models, one at a time ###### X.par.1 <- qf.1[1:(N+u-1-H+1), ] X.par.2 <- qf.2[1:(N+u-1-H+1), ] X.par.3 <- qf.3[1:(N+u-1-H+1), ] X.par.4 <- qf.4[1:(N+u-1-H+1), ] X.par.1 <- jitter(X.par.1[,15:28]) X.par.2 <- jitter(X.par.2[,15:28]) X.par.3 <- jitter(X.par.3[,15:28]) X.par.4 <- jitter(X.par.4[,15:28]) constant <- matrix(1,(N+u-1-H+1),1) QQ.l.1 = cbind(constant,X.par.1) QQ.l.2 = cbind(constant,X.par.2) QQ.l.3 = cbind(constant,X.par.3) QQ.l.4 = cbind(constant,X.par.4) # FIVE METHODS IN ALL [1st FOUR METHODS are "OLD" benchmarks] Y.est.np = Y.est[2:(N+u-H+1)] # Method (5) LINEAR FIT # LASSO Penalized Quantile Averaging - general via # linear parametric conditional # marginal quantiles fits QQ.l.1, ..., QQ.l.4 (with constant) # Lambda[1,]= drastic, Lambda[2,]= CV, Lambda[3,]= Sherwood & Wang, # Lambda{4,]= standard, Lambda[5,]= stringent penalty # Method 5.1 lam.3 <- Set.lam(X.par.1,Y.est.np,tt) # selecting penalty: S&W choice fit_model <- rq.lasso.fit(X.par.1,Y.est.np,tau=tt,lambda=lam.3) fit.5.1 <- as.ts(fit_model$coefficients) X.pred <- QQ.l.1 Y.pred.l.lasso <- X.pred %*% fit.5.1 # 1-step ahead prediction P.err.Meth.5.1[u] <- Y.est[N+u-H+1]-Y.pred.l.lasso[(N+u-1-H+1)] fit.5.1.nz <- which(fit.5.1 != 0) # remove zeros nr.fit.5.1 <- length(fit.5.1.nz)-1 # length - constant if (nr.fit.5.1 >= 1){ S.Meth.5.1[u,1:nr.fit.5.1] = fit.5.1.nz[2:(nr.fit.5.1+1)]-1 # exclude constant, and shift all indices 1 lower (-1 here) } # Method 5.2 lam.3 <- Set.lam(X.par.2,Y.est.np,tt) # selecting penalty: S&W choice fit_model <- rq.lasso.fit(X.par.2,Y.est.np,tau=tt,lambda=lam.3) fit.5.2 <- as.ts(fit_model$coefficients) X.pred <- QQ.l.2 Y.pred.l.lasso <- X.pred %*% fit.5.2 # 1-step ahead prediction P.err.Meth.5.2[u] <- Y.est[N+u-H+1]-Y.pred.l.lasso[(N+u-1-H+1)] fit.5.2.nz <- which(fit.5.2 != 0) # remove zeros nr.fit.5.2 <- length(fit.5.2.nz)-1 # length - constant if (nr.fit.5.2 >= 1){ S.Meth.5.2[u,1:nr.fit.5.2] = fit.5.2.nz[2:(nr.fit.5.2+1)]-1 # exclude constant, and shift all indices 1 lower (-1 here) } # Method 5.3 lam.3 <- Set.lam(X.par.3,Y.est.np,tt) # selecting penalty: S&W choice fit_model <- rq.lasso.fit(X.par.3,Y.est.np,tau=tt,lambda=lam.3) fit.5.3 <- as.ts(fit_model$coefficients) X.pred <- QQ.l.3 Y.pred.l.lasso <- X.pred %*% fit.5.3 # 1-step ahead prediction P.err.Meth.5.3[u] <- Y.est[N+u-H+1]-Y.pred.l.lasso[(N+u-1-H+1)] fit.5.3.nz <- which(fit.5.3 != 0) # remove zeros nr.fit.5.3 <- length(fit.5.3.nz)-1 # length - constant if (nr.fit.5.3 >= 1){ S.Meth.5.3[u,1:nr.fit.5.3] = fit.5.3.nz[2:(nr.fit.5.3+1)]-1 # exclude constant, and shift all indices 1 lower (-1 here) } # Method 5.4 lam.3 <- Set.lam(X.par.4,Y.est.np,tt) # selecting penalty: S&W choice fit_model <- rq.lasso.fit(X.par.4,Y.est.np,tau=tt,lambda=lam.3) fit.5.4 <- as.ts(fit_model$coefficients) X.pred <- QQ.l.4 Y.pred.l.lasso <- X.pred %*% fit.5.4 # 1-step ahead prediction P.err.Meth.5.4[u] <- Y.est[N+u-H+1]-Y.pred.l.lasso[(N+u-1-H+1)] fit.5.4.nz <- which(fit.5.4 != 0) # remove zeros nr.fit.5.4 <- length(fit.5.4.nz)-1 # length - constant if (nr.fit.5.4 >= 1){ S.Meth.5.4[u,1:nr.fit.5.4] = fit.5.4.nz[2:(nr.fit.5.4+1)] # exclude constant, and shift all indices 1 lower (-1 here) } #################################################################### # Method (6) = like Method (5) (LASSO penalty function), using NP # marginal quantiles [DeGooijer - Zerom, IJF, 2019 method] ##### NONPARAMETRIC (np) conditional marginal fits ############ options(warn=-1) q.l <- array(0,dim=c((N+u-1-H+1),ncol)) # nonp. linear fit for (i in 1:ncol){ X.cov.np <- X.H[1:(N+u-1-H+1),i] # marginal fits q.l.fit <- fitted(rq(Y.est.np ~ X.cov.np,tau=tt)) # linear (l) fit q.l[,i] <- q.l.fit } QQ.l.1 <- q.l QQ.l.1.new <- jitter(QQ.l.1,amount=2) lam.3 <- Set.lam(QQ.l.1.new,Y.est.np,tt) # selecting penalty: S&W choice fit_model.6 <- rq.lasso.fit(QQ.l.1.new,Y.est.np,tau=tt,lambda=lam.3) fit.6 <- as.ts(fit_model.6$coefficients) X.pred <- cbind(constant,q.l) Y.pred.l.lasso <- X.pred %*% fit.6 # 1-step ahead prediction P.err.Meth.6[u] <- Y.est[N+u-H+1]-Y.pred.l.lasso[(N+u-1-H+1)] ### To use with Method (8) below: fit.6.nz <- which(fit.6 != 0) # remove zeros nr.fit.6 <- length(fit.6.nz)-1 # length - constant S.sel.6 <- matrix(0,(N+u-1-H+1),nr.fit.6) if (nr.fit.6 >= 1){ for (i in 1:nr.fit.6){ S.sel.6[1:(N+u-1-H+1),i] = QQ.l.1.new[1:(N+u-1-H+1),fit.6.nz[i]] } } ## if (nr.fit.6 == 0){ ## S.sel.6[1:(N+u-1-H+1),1] = QQ.l.1[1:(N+u-2-H+1),fit.6.nz[i]] ## } if (nr.fit.6 >= 1){ S.Meth.6[u,1:nr.fit.6] = fit.6.nz[2:(nr.fit.6+1)]-1 # exclude constant, and shift all indices 1 lower (-1 here) } #################################################################### # Method (7) = LASSO penalized linear regression QQ.Meth.7 <- X.H[1:(N+u-1-H+1), ] lam.3.7 <- Set.lam(QQ.Meth.7,Y.est.np,tt) # selecting penalty: S&W choice fit_model.7 <- rq.lasso.fit(QQ.Meth.7,Y.est.np,tau=tt,lambda=lam.3.7) fit.7 <- as.ts(fit_model.7$coefficients) X.pred.7 <- cbind(constant,QQ.Meth.7) Y.pred.pen.7.lasso <- X.pred.7 %*% fit.7 # 1-step ahead prediction P.err.Meth.7[u] <- Y.est[N+u-H+1]-Y.pred.pen.7.lasso[(N+u-1-H+1)] fit.7.nz <- which(fit.7 != 0) # remove zeros nr.fit.7 <- length(fit.7.nz)-1 # length - constant S.sel.7 <- matrix(0,(N+u-1-H+1),nr.fit.7) ## if (nr.fit.7 >= 1){ ## S.Meth.7[u,1:nr.fit.7] = fit.7.nz[2:(nr.fit.7+1)]-1 ## # exclude constant, and shift all indices 1 lower (-1 here) ## } #################################################################### # Method (8) = Method (6) with nonzero variables selected # from averaging # based on fit.6.nz, mstop = 5000 so hardly # any penalization QQ.Meth.8 = S.sel.6 if (nr.fit.6 >= 1){ fit.8 <- predict(gamboost(Y.est.np ~ bols(QQ.Meth.8),control = boost_control(mstop =5000), family = QuantReg(tau=tt) ) ) P.err.Meth.8[u] <- Y.est[N+u-H+1]-fit.8[(N+u-1-H+1)] } else { P.err.Meth.8[u] <- 1111} S.Meth.8 = S.Meth.6 ################################################################################ # Method (9) = Hybrid: 14 averages from NP and 14 from Parametric (P) model # 28 potential predictors, [4 times = number of models] QQ.Meth.9.1 = cbind(qf.1[1:(N+u-1-H+1),],q.l[1:(N+u-1-H+1),]) # merge P and NP predictors QQ.Meth.9.2 = cbind(qf.2[1:(N+u-1-H+1),],q.l[1:(N+u-1-H+1),]) # merge P and NP predictors QQ.Meth.9.3 = cbind(qf.3[1:(N+u-1-H+1),],q.l[1:(N+u-1-H+1),]) # merge P and NP predictors QQ.Meth.9.4 = cbind(qf.4[1:(N+u-1-H+1),],q.l[1:(N+u-1-H+1),]) # merge P and NP predictors ##################################################################### # Method (10) = Hybrid: 14 averages from NP and 14 from P model # 28 potential predictors, [4 times = number of models] ## Method 10.1 f.10.1 <- rq(Y.est.np ~ QQ.Meth.9.1, method="lasso",lambda=30,tau=tt) Meth.10.1 <- qrfit.lasso(QQ.Meth.9.1,Y.est.np,tau=tt,method="admm",beta=f.10.1$coefficients) beta.10.1 <- Meth.10.1$beta b0.10.1 <- Meth.10.1$b fit.10.1 <- cbind(b0.10.1,t(beta.10.1)) # all coefficients new.10.1 <- beta.10.1[beta.10.1 >=0.01] fit.10.1.nz <- which(beta.10.1 >= 0.01) # indices X.pred <- cbind(matrix(1,(N+u-1-H+1),1),QQ.Meth.9.1) Y.10.1 <- X.pred %*% t(fit.10.1) # 1-step ahead prediction P.err.Meth.10.1[u] <- Y.est[N+u-H+1]-Y.10.1[(N+u-1-H+1)] nr.fit.10.1 <- length(fit.10.1.nz) if (nr.fit.10.1 >= 1){ S.Meth.10.1[u,1:nr.fit.10.1] <- fit.10.1.nz } ## Method 10.2 f.10.2 <- rq(Y.est.np ~ QQ.Meth.9.2, method="lasso",lambda=30,tau=tt) Meth.10.2 <- qrfit.lasso(QQ.Meth.9.2,Y.est.np,tau=tt,method="admm",beta=f.10.2$coefficients) beta.10.2 <- Meth.10.2$beta b0.10.2 <- Meth.10.2$b fit.10.2 <- cbind(b0.10.2,t(beta.10.2)) # all coefficients new.10.2 <- beta.10.2[beta.10.2 >=0.01] fit.10.2.nz <- which(beta.10.2 >= 0.01) # indices X.pred <- cbind(matrix(1,(N+u-1-H+1),1),QQ.Meth.9.2) Y.10.2 <- X.pred %*% t(fit.10.2) # 1-step ahead prediction P.err.Meth.10.2[u] <- Y.est[N+u-H+1]-Y.10.2[(N+u-1-H+1)] nr.fit.10.2 <- length(fit.10.2.nz) if (nr.fit.10.2 >= 1){ S.Meth.10.2[u,1:nr.fit.10.2] <- fit.10.2.nz } ## Method 10.3 f.10.3 <- rq(Y.est.np ~ QQ.Meth.9.3, method="lasso",lambda=30,tau=tt) Meth.10.3 <- qrfit.lasso(QQ.Meth.9.3,Y.est.np,tau=tt,method="admm",beta=f.10.3$coefficients) beta.10.3 <- Meth.10.3$beta b0.10.3 <- Meth.10.3$b fit.10.3 <- cbind(b0.10.3,t(beta.10.3)) # all coefficients new.10.3 <- beta.10.3[beta.10.3 >=0.01] fit.10.3.nz <- which(beta.10.3 >= 0.01) # indices X.pred <- cbind(matrix(1,(N+u-1-H+1),1),QQ.Meth.9.3) Y.10.3 <- X.pred %*% t(fit.10.3) # 1-step ahead prediction P.err.Meth.10.3[u] <- Y.est[N+u-H+1]-Y.10.3[(N+u-1-H+1)] nr.fit.10.3 <- length(fit.10.3.nz) if (nr.fit.10.3 >= 1){ S.Meth.10.3[u,1:nr.fit.10.3] <- fit.10.3.nz } ## Method 10.4 f.10.4 <- rq(Y.est.np ~ QQ.Meth.9.4, method="lasso",lambda=30,tau=tt) Meth.10.4 <- qrfit.lasso(QQ.Meth.9.4,Y.est.np,tau=tt,method="admm",beta=f.10.4$coefficients) beta.10.4 <- Meth.10.4$beta b0.10.4 <- Meth.10.4$b fit.10.4 <- cbind(b0.10.4,t(beta.10.4)) # all coefficients new.10.4 <- beta.10.4[beta.10.4 >=0.01] fit.10.4.nz <- which(beta.10.4 >= 0.01) # indices X.pred <- cbind(matrix(1,(N+u-1-H+1),1),QQ.Meth.9.4) Y.10.4 <- X.pred %*% t(fit.10.4) # 1-step ahead prediction P.err.Meth.10.4[u] <- Y.est[N+u-H+1]-Y.10.4[(N+u-1-H+1)] nr.fit.10.4 <- length(fit.10.4.nz) if (nr.fit.10.4 >= 1){ S.Meth.10.4[u,1:nr.fit.10.4] <- fit.10.4.nz } } # END BIG LOOP, u=[1:N.total (624)], for ONE quantile level ###### END SECOND BIG LOOP ################################### ###### SUMMARY MEASURES methods 5 - 10 (Method 9 excluded) ### ####### Method 5 ################# # Model 5.1 new benchmark (best selected models) S.5.1 <- S.Meth.5.1 # convert PE matrix in column vector S.5.1.nz <- S.5.1[S.5.1 !=0] # remove zeros x.5.1 <- S.5.1.nz # frequency table cat("Model 5.1") cat("Series ID ", "Freq.", " Cumulative", " Rel. Freq =100*Cum/total(Cum)") cbind(Freq=table(x.5.1),Cumul=cumsum(table(x.5.1)), Relative=prop.table(table(x.5.1))*100) # Model 5.2 new benchmark (best selected models) S.5.2 <- S.Meth.5.2 # convert PE matrix in column vector S.5.2.nz <- S.5.2[S.5.2 !=0] # remove zeros x.5.2 <- S.5.2.nz # frequency table cat("Model 5.2") cat("Series ID ", "Freq.", " Cumulative", " Rel. Freq =100*Cum/total(Cum)") cbind(Freq=table(x.5.2),Cumul=cumsum(table(x.5.2)), Relative=prop.table(table(x.5.2))*100) # Model 5.3 new benchmark (best selected models) S.5.3 <- S.Meth.5.3 # convert PE matrix in column vector S.5.3.nz <- S.5.3[S.5.3 !=0] # remove zeros x.5.3 <- S.5.3.nz # frequency table cat("Model 5.3") cat("Series ID ", "Freq.", " Cumulative", " Rel. Freq =100*Cum/total(Cum)") cbind(Freq=table(x.5.3),Cumul=cumsum(table(x.5.3)), Relative=prop.table(table(x.5.3))*100) # Model 5.4 new benchmark (best selected models) S.5.4 <- S.Meth.5.4 # convert PE matrix in column vector S.5.4.nz <- S.5.4[S.5.4 !=0] # remove zeros x.5.4 <- S.5.4.nz # frequency table cat("Model 5.4") cat("Series ID ", "Freq.", " Cumulative", " Rel. Freq =100*Cum/total(Cum)") cbind(Freq=table(x.5.4),Cumul=cumsum(table(x.5.4)), Relative=prop.table(table(x.5.4))*100) ############## Method 6 ################### # Model 6 (Results are the same for Model (8) new benchmark (best selected models) S.6 <- S.Meth.6 # convert PE matrix in column vector S.6.nz <- S.6[S.6 !=0] # remove zeros x.6 <- S.6.nz # frequency table cat("Model 6 = Model 8") cat("Series ID ", "Freq.", " Cumulative", " Rel. Freq =100*Cum/total(Cum)") cbind(Freq=table(x.6),Cumul=cumsum(table(x.6)), Relative=prop.table(table(x.6))*100) ############### Method 7 ################## # Model 7 new benchmark (best selected models) S.7 <- S.Meth.7 # convert PE matrix in column vector S.7.nz <- S.7[S.7 !=0] # remove zeros x.7 <- S.7.nz # frequency table cat("Model 7 ") cat("Series ID ", "Freq.", " Cumulative", " Rel. Freq =100*Cum/total(Cum)") cbind(Freq=table(x.7),Cumul=cumsum(table(x.7)), Relative=prop.table(table(x.7))*100) ######### Method 10 ######################### # Model 10.1 new benchmark (best selected models) S.10.1 <- S.Meth.10.1 # convert PE matrix in column vector S.10.1.nz <- S.10.1[S.10.1 !=0] # remove zeros x.10.1 <- S.10.1.nz # frequency table cat("Model 10.1") cat("Series ID ", "Freq.", " Cumulative", " Rel. Freq =100*Cum/total(Cum)") cbind(Freq=table(x.10.1),Cumul=cumsum(table(x.10.1)), Relative=prop.table(table(x.10.1))*100) # Model 10.2 new benchmark (best selected models) S.10.2 <- S.Meth.10.2 # convert PE matrix in column vector S.10.2.nz <- S.10.2[S.10.2 !=0] # remove zeros x.10.2 <- S.10.2.nz # frequency table cat("Model 10.2") cat("Series ID ", "Freq.", " Cumulative", " Rel. Freq =100*Cum/total(Cum)") cbind(Freq=table(x.10.2),Cumul=cumsum(table(x.10.2)), Relative=prop.table(table(x.10.2))*100) # Model 10.3 new benchmark (best selected models) S.10.3 <- S.Meth.10.3 # convert PE matrix in column vector S.10.3.nz <- S.10.3[S.10.3 !=0] # remove zeros x.10.3 <- S.10.3.nz # frequency table cat("Model 10.3") cat("Series ID ", "Freq.", " Cumulative", " Rel. Freq =100*Cum/total(Cum)") cbind(Freq=table(x.10.3),Cumul=cumsum(table(x.10.3)), Relative=prop.table(table(x.10.3))*100) # Model 10.4 new benchmark (best selected models) S.10.4 <- S.Meth.10.4 # convert PE matrix in column vector S.10.4.nz <- S.10.4[S.10.4 !=0] # remove zeros x.10.4 <- S.10.4.nz # frequency table cat("Model 10.4") cat("Series ID ", "Freq.", " Cumulative", " Rel. Freq =100*Cum/total(Cum)") cbind(Freq=table(x.10.4),Cumul=cumsum(table(x.10.4)), Relative=prop.table(table(x.10.4))*100) ########### results given in 1st part of Table 4 ############ ########### Thick loss (TI) functions ##### Ti.5.1 = mean((tt-(P.err.Meth.5.1<0))*P.err.Meth.5.1) Ti.5.2 = mean((tt-(P.err.Meth.5.2<0))*P.err.Meth.5.2) Ti.5.3 = mean((tt-(P.err.Meth.5.3<0))*P.err.Meth.5.3) Ti.5.4 = mean((tt-(P.err.Meth.5.4<0))*P.err.Meth.5.4) ## Ti.9.1 = mean((tt-(P.err.Meth.9.1<0))*P.err.Meth.9.1) ## Ti.9.2 = mean((tt-(P.err.Meth.9.2<0))*P.err.Meth.9.2) ## Ti.9.3 = mean((tt-(P.err.Meth.9.3<0))*P.err.Meth.9.3) ## Ti.9.4 = mean((tt-(P.err.Meth.9.4<0))*P.err.Meth.9.4) Ti.M6 = mean((tt-(P.err.Meth.6<0))*P.err.Meth.6) Ti.M7 = mean((tt-(P.err.Meth.7<0))*P.err.Meth.7) P.err.Meth.8.new = P.err.Meth.8[P.err.Meth.8 != 1111] #remove outliers/non convergence Ti.M8.new = mean((tt-(P.err.Meth.8.new<0))*P.err.Meth.8.new) Ti.10.1 = mean((tt-(P.err.Meth.10.1<0))*P.err.Meth.10.1) Ti.10.2 = mean((tt-(P.err.Meth.10.2<0))*P.err.Meth.10.2) Ti.10.3 = mean((tt-(P.err.Meth.10.3<0))*P.err.Meth.10.3) Ti.10.4 = mean((tt-(P.err.Meth.10.4<0))*P.err.Meth.10.4) ############################### Err.5.1 = ((tt-(P.err.Meth.5.1<0))*P.err.Meth.5.1) Err.5.2 = ((tt-(P.err.Meth.5.2<0))*P.err.Meth.5.2) Err.5.3 = ((tt-(P.err.Meth.5.3<0))*P.err.Meth.5.3) Err.5.4 = ((tt-(P.err.Meth.5.4<0))*P.err.Meth.5.4) ## Err.9.1 = ((tt-(P.err.Meth.9.1<0))*P.err.Meth.9.1) ## Err.9.2 = ((tt-(P.err.Meth.9.2<0))*P.err.Meth.9.2) ## Err.9.3 = ((tt-(P.err.Meth.9.3<0))*P.err.Meth.9.3) ## Err.9.4 = ((tt-(P.err.Meth.9.4<0))*P.err.Meth.9.4) Err.M6 = ((tt-(P.err.Meth.6<0))*P.err.Meth.6) Err.M7 = ((tt-(P.err.Meth.7<0))*P.err.Meth.7) P.err.Meth.8.new = P.err.Meth.8[P.err.Meth.8 != 1111] # remove outliers/non convergence Err.M8.new = ((tt-(P.err.Meth.8.new<0))*P.err.Meth.8.new) Err.10.1 = ((tt-(P.err.Meth.10.1<0))*P.err.Meth.10.1) Err.10.2 = ((tt-(P.err.Meth.10.2<0))*P.err.Meth.10.2) Err.10.3 = ((tt-(P.err.Meth.10.3<0))*P.err.Meth.10.3) Err.10.4 = ((tt-(P.err.Meth.10.4<0))*P.err.Meth.10.4) Err.EW = ((tt-(PE.EW<0))*PE.EW) ## Diebold-Mariano (DM), results given in 2nd part of Table 4 ## EW: Alternative=greater, i.e. method 2 is more accurate ## than method 1 ## p-value <= 0.05 reject H0 method 1 = method 2 ############################################################ T <- dm.test(Err.EW,Err.5.1,alternative=c("greater"),h=1,power=2) pv <- T$p.value pv.DM.EW.5.1 = pv T <- dm.test(Err.EW,Err.5.2,alternative=c("greater"),h=1,power=2) pv <- T$p.value pv.DM.EW.5.2 = pv T <- dm.test(Err.EW,Err.5.3,alternative=c("greater"),h=1,power=2) pv <- T$p.value pv.DM.EW.5.3 = pv T <- dm.test(Err.EW,Err.5.4,alternative=c("greater"),h=1,power=2) pv <- T$p.value pv.DM.EW.5.4 = pv T <- dm.test(Err.EW,Err.M6,alternative=c("greater"),h=1,power=2) pv <- T$p.value pv.DM.EW.M6 = pv T <- dm.test(Err.EW,Err.M7,alternative=c("greater"),h=1,power=2) pv <- T$p.value pv.DM.EW.M7 = pv T <- dm.test(Err.EW,Err.M8.new,alternative=c("greater"),h=1,power=2) pv <- T$p.value pv.DM.EW.M8.new = pv T <- dm.test(Err.EW,Err.10.1,alternative=c("greater"),h=1,power=2) pv <- T$p.value pv.DM.EW.10.1 = pv T <- dm.test(Err.EW,Err.10.2,alternative=c("greater"),h=1,power=2) pv <- T$p.value pv.DM.EW.10.2 = pv T <- dm.test(Err.EW,Err.10.3,alternative=c("greater"),h=1,power=2) pv <- T$p.value pv.DM.EW.10.3 = pv T <- dm.test(Err.EW,Err.10.4,alternative=c("greater"),h=1,power=2) pv <- T$p.value pv.DM.EW.10.4 = pv ############### Method 5 ####################################### #### CHANGING ORDER OF RATIOS AS COMPARED WITH ABOVE RESULTS #### (IN TABLE 4) #### Ratios of Tick loss functions; NEW benchmark models #### vs METHODS 5 -- 10 (excluding Method 9) #### Method 5 ########### # M1 R.5.1.n1 = Ti.5.1/Ti.1 # >1 = Ti.naive is worse than Ti.1 R.5.1.n2 = Ti.5.1/Ti.2 R.5.1.n3 = Ti.5.1/Ti.3 R.5.1.n4 = Ti.5.1/Ti.4 Ratio.Ti.5.1.new = cbind(R.5.1.n1,R.5.1.n2,R.5.1.n3,R.5.1.n4) # 1st index = new, 2nd "5" # print(Ratio.Ti.5.1.new) # M2 R.5.2.n1 = Ti.5.2/Ti.1 # >1 = Ti.naive is worse than Ti.1 R.5.2.n2 = Ti.5.2/Ti.2 R.5.2.n3 = Ti.5.2/Ti.3 R.5.2.n4 = Ti.5.2/Ti.4 Ratio.Ti.5.2.new = cbind(R.5.2.n1,R.5.2.n2,R.5.2.n3,R.5.2.n4) # 1st index = new, 2nd "5" # print(Ratio.Ti.5.2.new) # M3 R.5.3.n1 = Ti.5.3/Ti.1 # >1 = Ti.naive is worse than Ti.1 R.5.3.n2 = Ti.5.3/Ti.2 R.5.3.n3 = Ti.5.3/Ti.3 R.5.3.n4 = Ti.5.3/Ti.4 Ratio.Ti.5.3.new = cbind(R.5.3.n1,R.5.3.n2,R.5.3.n3,R.5.3.n4) # 1st index = new, 2nd "5" # print(Ratio.Ti.5.3.new) # M4 R.5.4.n1 = Ti.5.4/Ti.1 # >1 = Ti.naive is worse than Ti.1 R.5.4.n2 = Ti.5.4/Ti.2 R.5.4.n3 = Ti.5.4/Ti.3 R.5.4.n4 = Ti.5.4/Ti.4 Ratio.Ti.5.4.new = cbind(R.5.4.n1,R.5.4.n2,R.5.4.n3,R.5.4.n4) # 1st index = new, 2nd "5" # print(Ratio.Ti.5.4.new) #### (IN TABLE 4) #### Method 6 ########### R.M6.n1 = Ti.M6/Ti.1 R.M6.n2 = Ti.M6/Ti.2 R.M6.n3 = Ti.M6/Ti.3 R.M6.n4 = Ti.M6/Ti.4 Ratio.Ti.M6.new = cbind(R.M6.n1,R.M6.n2,R.M6.n3,R.M6.n4) # print(Ratio.Ti.M6.new) #### (IN TABLE) #### Method 7 ########### R.M7.n1 = Ti.M7/Ti.1 R.M7.n2 = Ti.M7/Ti.2 R.M7.n3 = Ti.M7/Ti.3 R.M7.n4 = Ti.M7/Ti.4 Ratio.Ti.M7.new = cbind(R.M7.n1,R.M7.n2,R.M7.n3,R.M7.n4) # print(Ratio.Ti.M7.new) #### (IN TABLE 4) #### Method 8 ########### R.M8.n1 = Ti.M8.new/Ti.1 R.M8.n2 = Ti.M8.new/Ti.2 R.M8.n3 = Ti.M8.new/Ti.3 R.M8.n4 = Ti.M8.new/Ti.4 Ratio.Ti.M8.new = cbind(R.M8.n1,R.M8.n2,R.M8.n3,R.M8.n4) # print(Ratio.Ti.M8.new) #### Method 10 ######### #### (IN TABLE 4) # M1 R.10.1.n1 = Ti.10.1/Ti.1 # If >1 = Ti.new.1 is worse than Ti.10 R.10.1.n2 = Ti.10.1/Ti.2 R.10.1.n3 = Ti.10.1/Ti.3 R.10.1.n4 = Ti.10.1/Ti.4 Ratio.Ti.10.new1 = cbind(R.10.1.n1,R.10.1.n2,R.10.1.n3,R.10.1.n4) # M2 R.10.2.n1 = Ti.10.2/Ti.1 # If >1 = Ti.new.1 is worse than Ti.10 R.10.2.n2 = Ti.10.2/Ti.2 R.10.2.n3 = Ti.10.2/Ti.3 R.10.2.n4 = Ti.10.2/Ti.4 Ratio.Ti.10.new2 = cbind(R.10.2.n1,R.10.2.n2,R.10.2.n3,R.10.2.n4) # M3 R.10.3.n1 = Ti.10.3/Ti.1 # If >1 = Ti.new.1 is worse than Ti.10 R.10.3.n2 = Ti.10.3/Ti.2 R.10.3.n3 = Ti.10.3/Ti.3 R.10.3.n4 = Ti.10.3/Ti.4 Ratio.Ti.10.new3 = cbind(R.10.3.n1,R.10.3.n2,R.10.3.n3,R.10.3.n4) # M4 R.10.4.n1 = Ti.10.4/Ti.1 # If >1 = Ti.new.1 is worse than Ti.10 R.10.4.n2 = Ti.10.4/Ti.2 R.10.4.n3 = Ti.10.4/Ti.3 R.10.4.n4 = Ti.10.4/Ti.4 Ratio.Ti.10.new4 = cbind(R.10.4.n1,R.10.4.n2,R.10.4.n3,R.10.4.n4) ######################################################## ######################################################## e.BM.1 = (tt-(P.err.min.all[ ,1]<0))*P.err.min.all[ ,1] e.BM.2 = (tt-(P.err.min.all[ ,2]<0))*P.err.min.all[ ,2] e.BM.3 = (tt-(P.err.min.all[ ,3]<0))*P.err.min.all[ ,3] e.BM.4 = (tt-(P.err.min.all[ ,4]<0))*P.err.min.all[ ,4] e.5.1 = (tt-(P.err.Meth.5.1<0))*P.err.Meth.5.1 e.5.2 = (tt-(P.err.Meth.5.2<0))*P.err.Meth.5.2 e.5.3 = (tt-(P.err.Meth.5.3<0))*P.err.Meth.5.3 e.5.4 = (tt-(P.err.Meth.5.4<0))*P.err.Meth.5.4 e.M6 = (tt-(P.err.Meth.6<0))*P.err.Meth.6 e.M7 = (tt-(P.err.Meth.7<0))*P.err.Meth.7 e.M8 = (tt-(P.err.Meth.8.new<0))*P.err.Meth.8.new e.10.1 = (tt-(P.err.Meth.10.1<0))*P.err.Meth.10.1 e.10.2 = (tt-(P.err.Meth.10.2<0))*P.err.Meth.10.2 e.10.3 = (tt-(P.err.Meth.10.3<0))*P.err.Meth.10.3 e.10.4 = (tt-(P.err.Meth.10.4<0))*P.err.Meth.10.4 ###### Method BM ########## DM.H.1.BM.1 <- dm.test(e.10.1, e.BM.1,alternative = c("less"), h = 1, power = 2) pv <- DM.H.1.BM.1$p.value pv.DM.10.1.BM.1 = pv DM.H.2.BM.2 <- dm.test(e.10.2, e.BM.2,alternative = c("less"), h = 1, power = 2) pv <- DM.H.2.BM.2$p.value pv.DM.10.2.BM.2 = pv DM.H.3.BM.3 <- dm.test(e.10.3, e.BM.3,alternative = c("less"), h = 1, power = 2) pv <- DM.H.3.BM.3$p.value pv.DM.10.3.BM.3 = pv DM.H.4.BM.4 <- dm.test(e.10.4, e.BM.4,alternative = c("less"), h = 1, power = 2) pv <- DM.H.24BM.4$p.value pv.DM.10.4.BM.4 = pv ##################################################################### ################ DM p-values included in Table 4 ################### DM.H.EW.BM.4 <- dm.test(Err.EW,e.BM.4,alternative=c("greater"),h=1,power=2) pv <- DM.H.EW.BM.4$p.value pv.DM.EW.BM.4= pv DM.H.5.4.BM.4 <- dm.test(e.5.4, e.BM.4,alternative = c("less"), h = 1, power = 2) pv <- DM.H.5.4.BM.4$p.value pv.DM.5.4.BM.4 = pv DM.H.M6.BM.4 <- dm.test(e.M6, e.BM.4,alternative = c("less"), h = 1, power = 2) pv <- DM.H.M6.BM.4$p.value pv.DM.M6.BM.4 = pv DM.H.M7.BM.4 <- dm.test(e.M7, e.BM.4,alternative = c("less"), h = 1, power = 2) pv <- DM.H.M7.BM.4$p.value pv.DM.M7.BM.4 = pv DM.H.M8.BM.4 <- dm.test(e.M8, e.BM.4,alternative = c("less"), h = 1, power = 2) pv <- DM.H.M8.BM.4$p.value pv.DM.M8.BM.4 = pv DM.H.10.4.BM.4 <- dm.test(e.10.4, e.BM.4,alternative = c("less"), h = 1, power = 2) pv <- DM.H.10.4.BM.4$p.value pv.DM.10.4.BM.4 = pv ########################################################################## ###### Method 5 ########## DM.H.1.5.1 <- dm.test(e.10.1, e.5.1,alternative = c("less"), h = 1, power = 2) pv <- DM.H.1.5.1$p.value pv.DM.10.1.5.1 = pv DM.H.2.5.2 <- dm.test(e.10.2, e.5.2,alternative = c("less"), h = 1, power = 2) pv <- DM.H.2.5.2$p.value pv.DM.10.2.5.2 = pv DM.H.3.5.3 <- dm.test(e.10.3, e.5.3,alternative = c("less"), h = 1, power = 2) pv <- DM.H.3.5.3$p.value pv.DM.10.3.5.3 = pv DM.H.4.5.4 <- dm.test(e.10.4, e.5.4,alternative = c("less"), h = 1, power = 2) pv <- DM.H.4.5.4$p.value pv.DM.10.4.5.4 = pv ###### Method 6 ########## DM.H.1.M6 <- dm.test(e.10.1, e.M6,alternative = c("less"), h = 1, power = 2) pv <- DM.H.1.M6$p.value pv.DM.10.1.M6 = pv DM.H.2.M6 <- dm.test(e.10.2, e.M6,alternative = c("less"), h = 1, power = 2) pv <- DM.H.2.M6$p.value pv.DM.10.2.M6 = pv DM.H.3.M6 <- dm.test(e.10.3, e.M6,alternative = c("less"), h = 1, power = 2) pv <- DM.H.3.M6$p.value pv.DM.10.3.M6 = pv DM.H.4.M6 <- dm.test(e.10.4, e.M6,alternative = c("less"), h = 1, power = 2) pv <- DM.H.4.M6$p.value pv.DM.10.4.M6 = pv ###### Method 7 ########## DM.H.1.M7 <- dm.test(e.10.1, e.M7,alternative = c("less"), h = 1, power = 2) pv <- DM.H.1.M7$p.value pv.DM.10.1.M7 = pv DM.H.2.M7 <- dm.test(e.10.2, e.M7,alternative = c("less"), h = 1, power = 2) pv <- DM.H.2.M7$p.value pv.DM.10.2.M7 = pv DM.H.3.M7 <- dm.test(e.10.3, e.M7,alternative = c("less"), h = 1, power = 2) pv <- DM.H.3.M7$p.value pv.DM.10.3.M7 = pv DM.H.4.M7 <- dm.test(e.10.4, e.M7,alternative = c("less"), h = 1, power = 2) pv <- DM.H.4.M7$p.value pv.DM.10.4.M7 = pv ######### Method 8 ############ # NOTE; e.M8 has smaller length than 624 DM.H.1.M8 <- dm.test(e.10.1[1:length(e.M8)], e.M8,alternative = c("less"), h = 1, power = 2) pv <- DM.H.1.M8$p.value pv.DM.10.1.M8 = pv DM.H.2.M8 <- dm.test(e.10.2[1:598], e.M8,alternative = c("less"), h = 1, power = 2) pv <- DM.H.2.M8$p.value pv.DM.10.2.M8 = pv DM.H.3.M8 <- dm.test(e.10.3[1:598], e.M8,alternative = c("less"), h = 1, power = 2) pv <- DM.H.3.M8$p.value pv.DM.10.3.M8 = pv DM.H.4.M8 <- dm.test(e.10.4[1:598], e.M8,alternative = c("less"), h = 1, power = 2) pv <- DM.H.4.M8$p.value pv.DM.10.4.M8 = pv # For alternative="less", the alternative hypothesis is that # method 2 is less accurate than method 1. For alternative="greater", # the alternative hypothesis is that method 2 is more # accurate than method 1. ############################################################### # cat("Diebold-Mariano:") # print(pv.DM.EW.BM.4) # print(pv.DM.5.4.BM.4) # print(pv.DM.M6.BM.4) # print(pv.DM.M7.BM.4) # print(pv.DM.M8.BM.4) # print(pv.DM.10.4.BM.4) ################################### # cat("DM 10.1 - 10.4 vs. BM.1 - BM.4","\n\n") # print(pv.DM.10.1.BM.1) # print(pv.DM.10.2.BM.2) # print(pv.DM.10.3.BM.3) # print(pv.DM.10.4.BM.4) # cat("DM 10.1 - 10.4 vs. 5.1 - 5.4","\n\n") # print(pv.DM.10.1.5.1) # print(pv.DM.10.2.5.2) # rint(pv.DM.10.3.5.3) # print(pv.DM.10.4.5.4) # cat("DM 10.1 - 10.4 vs. M6","\n\n") # print(pv.DM.10.1.M6) # print(pv.DM.10.2.M6) # print(pv.DM.10.3.M6) # print(pv.DM.10.4.M6) # cat("DM 10.1 - 10.4 vs. M7","\n\n") # print(pv.DM.10.1.M7) # print(pv.DM.10.2.M7) # print(pv.DM.10.3.M7) # print(pv.DM.10.4.M7) # cat("DM 10.1 - 10.4 vs. M8","\n\n") # print(pv.DM.10.1.M8) # print(pv.DM.10.2.M8) # print(pv.DM.10.3.M8) # print(pv.DM.10.4.M8) ######################################################### ######################################################### ############ PRINT OUTPUT ###################### cat("H-step ahead forecasts, H= ", H, "\n\n") cat("tau = ", tt, "\n\n") # quantile level ######### Diebold-Mariano ########### cat("DM EW vs. 5.1 - 5.4","\n\n") print(pv.DM.EW.5.1) print(pv.DM.EW.5.2) print(pv.DM.EW.5.3) print(pv.DM.EW.5.4) cat("DM EW vs. M6","\n\n") print(pv.DM.EW.M6) cat("DM EW vs. M7","\n\n") print(pv.DM.EW.M7) cat("DM vs. M8.new","\n\n") print(pv.DM.EW.M8.new) cat("DM EW vs. 10.1 - 10.4","\n\n") print(pv.DM.EW.10.1) print(pv.DM.EW.10.2) print(pv.DM.EW.10.3) print(pv.DM.EW.10.4) ###################################### cat("Ratios of Ti:") cat("EW 1-4 1st col vs. BM 1-4 2nd col:","\n\n") print(Ratio.naive1.new) # naive = EM, Ti.1 - Ti.4 = BM print(Ratio.naive2.new) print(Ratio.naive3.new) print(Ratio.naive4.new) cat("PA-L 1-4 1st col vs. BM 1-4 2nd col:") print(Ratio.Ti.5.1.new) print(Ratio.Ti.5.2.new) print(Ratio.Ti.5.3.new) print(Ratio.Ti.5.4.new) print(Ratio.Ti.M6.new) print(Ratio.Ti.M7.new) print(Ratio.Ti.M8.new) print(Ratio.Ti.10.new1) print(Ratio.Ti.10.new2) print(Ratio.Ti.10.new3) print(Ratio.Ti.10.new4) ######################################################### ################## Codes of selected indices ########### # (Parametric Model no.,var) Macroeconomic variables # (1,1) (2,29) (3,57) (4,85) DP # (1,2) (2,30) (3,58) (4,86) DY # (1,3) (2,31) (3,59) (4,87) EP # (1,4) (2,32) (3,60) (4,88) DE # (1,5) (2,33) (3,61) (4,89) RVOL # (1,6) (2,34) (3,62) (4,90) BM # (1,7) (2,35) (3,63) (4,91) NTIS # (1,8) (2,36) (3,64) (4,92) TBL # (1,9) (2,37) (3,65) (4,93) LTY # (1,10) (2,38) (3,66) (4,94) LTR # (1,11) (2,39) (3,67) (4,95) TMS # (1,12) (2,40) (3,68) (4,96) DFY # (1,13) (2,41) (3,69) (4,97) DFR # (1,14) (2,42) (3,70) (4,98) INFL # # Technical variables # (1,15) (2,43) (3,71) (4,99) MA(1,9) # (1,16) (2,44) (3,72) (4,100) MA(1,12) # (1,17) (2,45) (3,73) (4,101) MA(2,9) # (1,18) (2,46) (3,74) (4,102) MA(2,12) # (1,19) (2,47) (3,75) (4,103) MA(3,9) # (1,20) (2,48) (3,76) (4,104) MA(3,12) # (1,21) (2,49) (3,77) (4,105) MOM(9) # (1,22) (2,50) (3,78) (4,106) MOM(12) # (1,23) (2,51) (3,79) (4,107) VOL(1,9) # (1,24) (2,52) (3,80) (4,108) VOL(1,12) # (1,25) (2,53) (3,81) (4,109) VOL(2,9) # (1,26) (2,54) (3,82) (4,110) VOL(2,12) # (1,27) (2,55) (3,83) (4,111) VOL(3,9) # (1,28) (2,56) (3,84) (4,112) VOL(3,12) # Linear fit: # 113-DP, 114-DY, 115-EP, 116-DE, 117-RVOL, 118-BM, 119-NTIS, 120-TBL , # 121-LTY, 122-LTR, 123-TMS, 124-DFY, 125-DFR, 126-INFL, # # 127-MA(1,9), 128-MA(1,12), 129-MA(2,9), 130-MA(2,12), 131-MA(3,9), # 132-MA(3,12),133-MOM(9), 134-MOM(19), 135-MOM(12), 136-VOL(1,9), # 137-VOL(2,9),138-VOL(2,12),139-VOL(3,9),140-VOL(3,12) # Nonparametric, BOOSTING: # 113-DP, 114-DY, 115-EP, 116-DE, 117-RVOL, 118-BM, 119-NTIS 120-TBL, # 121-LTY, 122-LTR, 123-TMS, 124-DFY, 125-DFR, 126-INFL # Direct # 1-DP, 2-DY, 3-EP, 4-DE, 5-RVOL, 6-BM, 7-NTIS, 8-TBL, 9-LTY, 10-LTR, # 11-TMS, 12-DFY, 13-DFR, 14-INFL # # 15-MA(1,9), 16-MA(1,12), 17-MA(2,9), 18-MA(2,12), 19-MA(3,9), # 20-MA(3,12), 21-MOM(9), 22-MOM(19), 23-MOM(12), 24-VOL(1,9), # 25-VOL(2,9), 26-VOL(2,12), 27-VOL(3,9), 28-VOL(3,12) # ################### END PRINTING OUTPUT #################################