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.