# Forecasting accuracy of bi and multivariate timeseries

First of all I'd like to point out that I was unsure whether to post this question here or on stack overflow, considering it is rather code heavy; but as I think it plays more on the logic and application of cross validation, a topic this forum has more expertise in than the other, I chose to post it here. Please move it to stack-overflow if deemed inappropriate.

My question is as follows:

When attempting to run the below code, adapted from this tutorial by Rob Hyndman, I run into trouble when populating the mae matrix in the below for loop. This I gather is from the handling of dx and its length, which due to its bivariate nature is double the actual time frame (more of course if multivariate).

Now, while I have an understanding of the principle behind using training and testing sets, I am not very familiar with how to apply it, and as a result am thoroughly stuck at the above problem.

Any and all help is greatly appreciated as I'm at my wits end, this being the last major hurdle in what to me is an important project.

library("forecast")
library("vars")

d <- rnorm(70)
x <- rnorm(70)

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

# Forecast Accuracy
k <- 58 # data length less forecast horison (as minimum)
n <- length(dx)
mae <- matrix(NA, n-k, 12)
st <- tsp(dx)[1]+(k-2)/12

for (i in 1:(n-k)) {
dxshort <- window(dx, end=st+i/12)
dxnext <- window(dx, start=st + (i+1)/12, end=st+(i+12)/12)
fit <- VAR(dxshort, p = 2)
fcast <- forecast(fit, h = 12)
fcastmean <- do.call('cbind', fcast[['mean']])
mae[i,1:length(dxnext)] <- (fcastmean - dxnext)
}


The above code uses random numbers with a normal distribution for illustration and should run as is.