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?
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]))
}