Forecasting methodology and k-fold cross validation for a vector autoregression

Question

This is a follow up question the question that can be found here, and is a result of me having implemented (after as careful evaluation as I'm capable of) the alterations and changes suggested.

Below is my method and should be replicable.

My question relates to the implementation of k-fold cross validation and whether the code produces a mean average error value that is reliable and whether there are some aspects of k-fold cross validation I may have neglected, thus skewing any results.

Otherwise any comments, both as to the method as it stands or the logic behind their inclusion (see above link) is welcome.

library(plyr)
library(forecast)
library(vars)

#Read Data
da=read.table("VARdata.txt", header=T)
dac <- c(2,3) # Select variables
x=da[,dac]

plot.ts(x)
summary(x)

#Run Augmented Dickey-Fuller tests to determine stationarity and
#differences to achieve stationarity.
adf1 <- ur.df(x[,"VAR1"], type = "drift", lags = 10, selectlags = "AIC")
adf2 <- ur.df(x[,"VAR2"], type = "drift", lags = 10, selectlags = "AIC")

summary(adf1)
summary(adf2)

#Difference to achieve stationarity
d.x1 = diff(x[, "VAR1"], differences = 1)
d.x2 = diff(x[, "VAR2"], differences = 1)


#Check if differenced variables are stationary
adf1b <- ur.df(d.x1, type = "drift", lags = 10, selectlags = "AIC")
adf2b <- ur.df(d.x2, type = "drift", lags = 10, selectlags = "AIC")

summary(adf1b)
summary(adf2b)

#If variable is stationary I(0), do not difference
#Shorten undifferenced variable by n, so as to make all variables same length
# d.x2 = (x[, "VAR2"])
# d.x2 = d.x2[-c(1:1)]

#Bind variables in time series
dx = cbind(d.x1, d.x2)
dx = as.ts(dx)
plot.ts(dx)

summary(dx)

#Lag optimisation
VARselect(dx, lag.max = 10, type = "both")

#Run VAR 
var = VAR(dx, p=2)

#Test for serial autocorrelation using the Portmanteau test
#Rerun var model with other suggested lags if H0 can be rejected at 0.05
serial.test(var, lags.pt = 10, type = "PT.asymptotic")

#ARCH test (Autoregressive conditional heteroscedasdicity)
arch.test(var, lags.multi = 10)

summary(var)

#Forecasting
prd <- forecast(var, h = 12)

print(prd)
plot(prd)

# Forecast Accuracy
data <- as.data.frame(dx)

k = 10 #Folds

# sample from 1 to k, nrow times (the number of observations in the data)
data$id <- sample(1:k, nrow(data), replace = TRUE)
list <- 1:k

# prediction and testset data frames that we add to with each iteration over
# the folds

prediction <- data.frame()
testsetCopy <- data.frame()

#Creating a progress bar to know the status of CV
progress.bar <- create_progress_bar("text")
progress.bar$init(k)

for (i in 1:k){
  # remove rows with id i from dataframe to create training set
  # select rows with id i to create test set
  trainingset <- subset(data, id %in% list[-i])
  trainingset <- as.ts(trainingset)
  testset <- subset(data, id %in% c(i))

  # run a VAR model
  mymodel <- VAR(trainingset, p = 2)

  # remove response column 1
  temp <- forecast(mymodel, h = nrow(testset))
  temp <- do.call('cbind', temp[['mean']])
  temp <- as.data.frame(temp)

  # append this iteration's predictions to the end of the prediction data frame
  prediction <- rbind(prediction, temp)

  # append this iteration's test set to the test set copy data frame
  # keep only the desired Column
  testsetCopy <- rbind(testsetCopy, as.data.frame(testset[,1]))

  progress.bar$step()
}

# add predictions and actual values
result <- cbind(prediction, testsetCopy[, 1])
names(result) <- c("Predicted", "Actual")
result$Difference <- abs(result$Actual - result$Predicted)

# As an example use Mean Absolute Error as Evalution 
summary(result$Difference)
result

Edit based on answer below:

As per the answer below I have changed the code for the cross validation to this (full test code included for ease):

library(forecast)
library(vars)
library(plyr)

x <- rnorm(70)
y <- rnorm(70)

dx <- cbind(x,y)
dx <- as.ts(dx)

j = 12  #Forecast horizon
k = nrow(dx)-j #length of minimum training set

prediction <- data.frame()
actual <- data.frame()

for (i in j) { 
  trainingset <- window(dx, end = k+i-1)
  testset <- window(dx, start = k-j+i+1, end = k+j)
  fit <- VAR(trainingset, p = 2)                       
  fcast <- forecast(fit, h = j)
  fcastmean <- do.call('cbind', fcast[['mean']])
  fcastmean <- as.data.frame(fcastmean)

  prediction <- rbind(prediction, fcastmean)
  actual <- rbind(actual, as.data.frame(testset[,1]))
}

# add predictions and actual values
result <- cbind(prediction, actual[, 1])
names(result) <- c("Predicted", "Actual")
result$Difference <- abs(result$Actual - result$Predicted)

# Use Mean Absolute Error as Evalution 
summary(result$Difference)

Would this be a better application of cross validation? I realize that it is no longer k-fold, but is based on the link provided in the answer.

Answer 1

I did not read your code too carefully, but so far I think it might not be an appropriate implementation of cross validation in a time series setting. You seem to be first omitting randomly some rows from the data set that only includes the original variables but not their lags and then implicitly forming the design matrix consisting of lagged values of the variables and estimating the VAR model. This way you mess up the time series ordering of series when forming the lagged variables. Time series cross validation does not work that way.

You can only remove rows from the full-sample design matrix (that already includes all the lags to be used in the model), and even that requires some conditions to be fulfilled, see Bergmeir et al. "A Note on the Validity of Cross-Validation for Evaluating Time Series Prediction" (2015, working paper).

More generally, you could be doing cross validation using rolling windows as described in Hyndman and Athanasopoulos "Forecasting: principles and practice" section 2.5 (scroll all the way down). See also my recent answer in this thread.