Solved – Selecting ARIMA Order using Rolling Forecast

arimaforecastingmodel selectionmoving windowtime series

I'm wondering if a rolling forecast technique like the ones mentioned in Rob Hyndman's blogs, and the example below, could be used to select the order for an ARIMA model?

In the examples I've looked at, like the ones below, it seems like the order of the ARIMA model is already specified, or is determined once by auto.arima and then the single model is evaluated using the forloop in the rolling forecast.

I'm wondering how you could use the rolling forecast technique to select the order of the ARIMA model. If anyone has a suggestion or example, that would be great.

Examples:
http://robjhyndman.com/hyndsight/tscvexample/
http://robjhyndman.com/hyndsight/rolling-forecasts/

Code:

library("fpp")

h <- 5
train <- window(hsales,end=1989.99)
test <- window(hsales,start=1990)
n <- length(test) - h + 1
fit <- auto.arima(train)
fc <- ts(numeric(n), start=1990+(h-1)/12, freq=12)
for(i in 1:n)
{  
  x <- window(hsales, end=1989.99 + (i-1)/12)
  refit <- Arima(x, model=fit)
  fc[i] <- forecast(refit, h=h)$mean[h]
}

Update:

Pseudo code:

library("fpp")

h <- 5
train <- window(hsales,end=1989.99)
test <- window(hsales,start=1990)
n <- length(test) - h + 1

##Create models for all combinations of p 10 to 0, d 2 to 0, q 10 to 0

fit1 <- Arima(train, order=c(10,2,10)
fit2 <- Arima(train, order=c(9,2,10)
fit3 <- Arima(train, order=c(8,2,10)
.
.
.
fit10 <- Arima(train, order=c(0,2,10)
fc1 <- ts(numeric(n), start=1990+(h-1)/12, freq=12)
fc2 <- ts(numeric(n), start=1990+(h-1)/12, freq=12)
fc3 <- ts(numeric(n), start=1990+(h-1)/12, freq=12)
.
.
.
fc10 <- ts(numeric(n), start=1990+(h-1)/12, freq=12)
for(i in 1:n)
{  
  x <- window(hsales, end=1989.99 + (i-1)/12)
  refit1 <- Arima(x, model=fit1)
  refit2 <- Arima(x, model=fit2)
  refit3 <- Arima(x, model=fit3)
  .
  .
  .
  refit10 <- Arima(x, model=fit10)
  fc1[i] <- forecast(refit1, h=h)$mean[h]
	  fc2[i] <- forecast(refit2, h=h)$mean[h]
  fc3[i] <- forecast(refit3, h=h)$mean[h]
	  .
	  .
	  .
	  fc10[i] <- forecast(refit10, h=h)$mean[h]
}

##Calculating mape for forecasts

Accuracy(fc1$mean,test)[,5]
	Accuracy(fc2$mean,test)[,5]
Accuracy(fc3$mean,test)[,5]
	.
	.
	.
	Accuracy(fc10$mean,test)[,5]

##Return the order of the Arima model that has the lowest mape

Best Answer

I'm wondering if a rolling forecast technique like the ones mentioned in Rob Hyndman's blogs, and the example below, could be used to select the order for an ARIMA model?

Rob J. Hyndman indicates in comments to his blog post "Time series cross-validation: an R example":

You would normally be trying several models and selecting the best based on a cross-validation procedure.

Cross-validation is both a method of measuring accuracy and a method for choosing a model. The model you finally use for forecasting is the one that gives the best cross-validation accuracy.

Also, since cross validation is often used for model selection for cross sectional data*, it is quite natural to do something similar for time series data (where regular cross validation is replaced by rolling-window cross validation).

*From another post called "Why every statistician should know about cross-validation":

Minimizing a CV statistic is a useful way to do model selection such as choosing variables in a regression or choosing the degrees of freedom of a nonparametric smoother.

I'm wondering how you could use the rolling forecast technique to select the order of the ARIMA model. If anyone has a suggestion or example, that would be great.

First, you choose a set of candidate models. For each model in the set, you evaluate forecasting performance based on rolling-window cross validation. Then you choose the model that delivers the best forecasting performance.

Here is an example I ran at some point to compare model selection based on rolling-window cross validation with AIC-based selection. (I wanted to illustrate that model selection based on rolling-window cross validation is asymptotically equivalent to AIC-based choice.)

# Generate a T-long sample of an ARMA(1,1) process
T =10^4
#T =2*10^3 # uncomment for a shorter series (10^3 rolling windows instead of 9*10^3)
T0=1*10^3 # the length of the rolling window
set.seed(1); innov1=rnorm(T); set.seed(2); innov2=rnorm(T)
x1=arima.sim(model=list(ar1=0.5,ma1=0.2),n=T,innov=innov1,n.start=10^3,start.innov=innov2)

# Estimate three candidate models (ARMA(1,1), ARMA(2,1), ARMA(0,1)) on T0-long rolling windows, 
# get 1-step-ahead mean squared forecast errors (MSFEs)
# (The loop below runs for about 15 minutes on a ThinkPad laptop with Sandy Bridge i5 processor produced in 2011.)
err1=err2=err3=rep(NA,T)
print(Sys.time()); for(end in T0:(T-1)){
 if(end%%100==0){
  print(paste("end =",end))
  print(Sys.time())
 }
 range=c((end-T0+1):end)
 model1  =arima(x1[range],order=c(1,0,1),seasonal=list(order=c(0,0,0),period=NA),
                xreg=NULL,include.mean=TRUE,method="CSS-ML",optim.method="BFGS")
 fcst1   =as.numeric(predict(object=model1,n.ahead=1)$pred)
     err1[end]=fcst1-x1[end+1]
     model2=  arima(x1[range],order=c(2,0,1),seasonal=list(order=c(0,0,0),period=NA),
                    xreg=NULL,include.mean=TRUE,method="CSS-ML",optim.method="BFGS")
     fcst2   =as.numeric(predict(object=model2,n.ahead=1)$pred)
 err2[end]=fcst2-x1[end+1]
 model3  =arima(x1[range],order=c(0,0,1),seasonal=list(order=c(0,0,0),period=NA),
                xreg=NULL,include.mean=TRUE,method="CSS-ML",optim.method="BFGS")
 fcst3   =as.numeric(predict(object=model3,n.ahead=1)$pred)
 err3[end]=fcst3-x1[end+1]
}; print(Sys.time())
err1_orig=err1; err1=head(tail(err1,-T0),-1); msfe1=mean(err1^2); print(paste("MSFE1 =",msfe1))
err2_orig=err2; err2=head(tail(err2,-T0),-1); msfe2=mean(err2^2); print(paste("MSFE2 =",msfe2))
err3_orig=err3; err3=head(tail(err3,-T0),-1); msfe3=mean(err3^2); print(paste("MSFE3 =",msfe3))

# Estimate the three candidate models on the full sample, obtain their AICs
 model1  =arima(x1[range],order=c(1,0,1),seasonal=list(order=c(0,0,0),period=NA),
                xreg=NULL,include.mean=TRUE,method="CSS-ML",optim.method="BFGS")
 AIC1=AIC(model1); print(paste("AIC1 =",AIC1))
 model2  =arima(x1[range],order=c(2,0,1),seasonal=list(order=c(0,0,0),period=NA),
                xreg=NULL,include.mean=TRUE,method="CSS-ML",optim.method="BFGS")
 AIC2=AIC(model2); print(paste("AIC2 =",AIC2))
 model3  =arima(x1[range],order=c(0,0,1),seasonal=list(order=c(0,0,0),period=NA),
                xreg=NULL,include.mean=TRUE,method="CSS-ML",optim.method="BFGS")
 AIC3=AIC(model3); print(paste("AIC3 =",AIC3))

# The ranking of the models by their 1-step-ahead mean squared forecast error 
# should ideally match the ranking of the models by their AICs.

# Indeed, model1 gets the lowest AIC and the lowest MSFE; model2 follows; model3 is the last. 
# Both rankings are consistent.
Related Question