Question

I'm trying to learn a vector autoregressive model using the vars package in R. This package doesn't have any way to measure the accuracy of the returned model.

Specifically, I want to use MASE as defined in the accuracy function from the forecast package in R to compare forecasting with VAR with forecasting using Arima models on each component time series (I'm using 4 possibly correlated time series). accuracy doesn't recognize the varest object returned by vars. How can I get the MASE for each forecasted component? I want to calculate both in-sample and out-of-sample accuracy

Code example:

library(vars)
library(forecast)
data(Canada)
v<- VAR(window(Canada, end=c(1998,4)), p=2)
accuracy(v$varresult[[1]])

The argument of accuracy is an lm object and returns the training accuracy for series 1 as:

                       ME      RMSE       MAE           MPE      MAPE       MASE
Training set 1.536303e-15 0.3346096 0.2653946 -1.288309e-05 0.0281736 0.03914555

I want to get the out-of-sample test accuracy using something like (not exactly this, since the forecast period needs to be specified):

 accuracy(v$varresult[[1]], window(Canada[,1], start=c(1999,1)))

but this is not supported for lm objects and returns the error

 Error in testaccuracy(f, x, test) : Unknown list structure

And if I use the values directly as follows, I don't get the MASE, which needs information about the training set. This is also prone to off-by-one errors because values are used and not ts objects, for which accuracy will match the stored times directly:

 p<-predict(v, n.ahead=8)
 accuracy(p$fcst[[1]][,"fcst"],window(Canada[,1], start=c(1999,1)))

             ME      RMSE       MAE         MPE       MAPE      ACF1 Theil's U
Test set -0.1058358 0.8585455 0.7385238 -0.01114099 0.07694492 0.5655117  1.359761

Ideally, I should like to forecast it as:

fr<-forecast(v$varresult[[1]], h=8)

but this can't work because it needs the other series for the prediction, and gives:

Error in eval(expr, envir, enclos) : object 'datamat' not found

I could try copying the functionality of forecast.Arima etc and try writing a forecast.varresult package, but is there a simpler way out?

Was it helpful?

Solution

Why don't you try reading the documentation. Here is what it says about the first argument, f:

An object of class "forecast", or a numerical vector containing forecasts. It will also work with Arima, ets and lm objects if x is omitted – in which case in-sample accuracy measures are returned.

VAR does not return an object of class "forecast", but you can compute a numerical vector containing forecasts.

Now read about the second argument, x.

An optional numerical vector containing actual values of the same length as object, or a time series overlapping with the times of f.

OK, that's pretty straightforward. Just give it the actual values in x and the forecast values in f.

But that won't give you the MASE as further down the help page it explains that the "MASE calculation is scaled using MAE of in-sample naive forecasts for non-seasonal time series, in-sample seasonal naive forecasts for seasonal time series and in-sample mean forecasts for non-time series data." So it can't do that calculation without the historical data, and unless you are passing an object of class 'forecast' it won't know about them.

However, it is not hard to trick it into giving what you want. Here is some code that does it:

trainingdata <- window(Canada, end=c(1998,4))
testdata <- window(Canada, start=c(1999,1))
v <- VAR(trainingdata, p=2)
p <- predict(v, n.ahead=8)
res <- residuals(v)
fits <- fitted(v)
for(i in 1:4)
{
  fc <- structure(list(mean=p$fcst[[i]][,"fcst"], x=trainingdata[,i],
     fitted=c(NA,NA,fits[,i])),class="forecast")
  print(accuracy(fc,testdata[,i]))
}
Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top