My aim is to compare the forecast performance of several time series models. I have a bivariate dataset, and applied three different models to it:
1) A univariate Arima model (applied to the first variable) using the automatic order selection function 'auto.arima()'. The estimated model is Arima(1,1,1)
2) A Verctorautoregression using both variables. The recognized model (using the package 'vars') is VAR(1), and
3) A univariate state space model (aggain applied to the first variable) using the package 'dlm'. I specified a state space form of Arima(1,1,1) model, as it was suggested by 'auto.arima()', namely I constructed a model consisting of a stochastic trend and of an arma model with one parameter ar and one ma (for details see the code).
I then generated forecasts compared the results graphically and was surprised to see, how poorly my state space model performs. The results of Arima are quite similar, only the VAR(2) model performs relatively well.
Is this poor result of state space model realistic or is my model specification wrong?
data <- read.table(...)
library(vars)
library(forecast)
library("dlm", lib.loc="C:/Users/incognito/Documents/R/win-library/3.0")
# subsetting the data:
data.s<-data[1:528,1]
# Estimation of univariate Arima model and generating a forecast:
arima.m<-auto.arima(data.s.g)
arima.f<-forecast(arima.m,h=30)
# Estimation of a state space representation of Arima(1,1,1) model and forecast:
level0 <- data.s.g[1]
slope0 <- mean(diff(data.s.g))
buildGap <- function(u) {
trend <- dlmModPoly(dV = 1e-7, dW = exp(u[1 : 2]),
m0 = c(level0, slope0),
C0 = 2 * diag(2))
gap <- dlmModARMA(ar = ARtransPars(u[4]),ma=u[5], sigma2 = exp(u[3]))
return(trend + gap)}
init <- c(-3, -1, -3, .4, .4)
outMLE <- dlmMLE(data.s.g, init, buildGap)
dlmGap <- buildGap(outMLE$par)
filt<-dlmFilter(data.s.g,dlmGap)
forc<-dlmForecast(filt,nAhead=30)
# A bivariate VAR model and forecast:
var<-VAR(data.s)
var.f<-predict(var,n.ahead=30)
# Plotting the results:
plot(data.s.g,xlim=c(400,560),ylim=c(1.5,4),type="l")
lines(529:558,forc$f)
lines(529:558,var.f$fcst$gas[,1],col=3)
lines(529:558,data$gas[529:558],col=4)
lines(529:558,arima.f$mean,col=2)
legend("topleft",legend=c("state space","arima","var"),lty=1,col=c(1,2,3))