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

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)

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

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

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

• Currently your problem is formulated in terms of proof-reading code, which is off-topic here. If you described the issue at hand and pointed out what exactly is troubling you, it could become on-topic. Mar 8, 2016 at 18:35
• @RichardHardy was afraid of that. My thought process was that as well as code review, it also covers specific things related to proper handling of k fold cross validation and therefore more suitable here. I'll clarify the question as soon as I have a chance! Mar 8, 2016 at 18:47
• @RichardHardy hopefully the question now satisfies the criteria for this forum. Mar 8, 2016 at 20:42