Selecting ARIMA Order using Rolling Forecast

Question

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 

Answer 1

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.

Answer 2

You can do this, but be aware that it's not that much different from inferring the lags from in-sample fit procedures such as ACF/PACF analysis. Why? Because you're using the forecast (supposedly out of sample) in the same way as in-sample for model selection. You're running multiple specifications through the same sample to pick the best in terms of some fit measure.

You could as well skip random points or chunks of sample in the middle of the full sample, then use state space estimation of arima to infer those value and compare them to the actuals. All these techniques make out-of-sample (forecast) testing very similar to in-sample, in that they weaken the power of these tests.

I claim that it's better to be upfront about what we're doing and simply extract the lag structure using in-sample analysis techniques. We won't pretend that we're doing out-of-sample, and it makes the results cleaner.

Otherwise, rolling forecasts especially one step ahead ones are linked to Chow test of parameter constancy. These are legitimate techniques, when used appropriately. For instance, you got your sample estimated the model in 2016, then in 2017 you got new data, which was not available at the time of modeling. So, you do rolling forecast, keeping the estimates from original model, and compare one-step ahead forecasts with new data. Chow test will provide you with a statistical measure of parameter constancy, e.g. it can detect intercept change.

You could try to apply this to the rolling forecasts, and it will look fine on the surface, but the truth is that in rolling forecasts, you have the forecast period's actuals available to you. It's just you pretend that you don't know them for the purpose of cross validation.