4
$\begingroup$

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 
$\endgroup$
2
$\begingroup$

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.
$\endgroup$
  • $\begingroup$ Thank you for your answer. This is what I suspected, you specify the models up front and then use the rolling forecast to compare their accuracy. I think what I have in mind is to use a loop to create models for all combinations of p = 10 to 1, d=2 to 0, q = 10 to 1, then use a rolling forecast to evaluate those models, and pick out the one that say has the lowest mape. Do you know of any examples like that, or can you suggest how to do it? $\endgroup$ – ndderwerdo Apr 24 '16 at 19:01
  • $\begingroup$ I don't know concrete examples, but the idea seems very natural. With regards to implementation, function arfimacv from "rugarch" package in R could be useful. $\endgroup$ – Richard Hardy Apr 24 '16 at 19:04
  • $\begingroup$ Thank you, I'll check those packages out. I've also added some pseudo code to make it clearer what I haven in mind. I'm not the greatest with loops so if you have any tips I'm all ears. $\endgroup$ – ndderwerdo Apr 24 '16 at 19:31
0
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ ACF/PACF work fine in simple cases like AR(p) or MA(q), but not for ARMA(p,q). Even so the typical model identification based on ACF/PACF would yield a too complex model, e.g. more complex than AIC-based selection which should be preferred when forecasting due to the efficiency property of the AIC. Rolling windows could be superior when the error distribution is quite different from the assumed distribution in which case AIC can be distorted, and also when the data generating process evolves over time (which will not be captured by fitting a model on the whole sample). $\endgroup$ – Richard Hardy Feb 13 '17 at 17:51
  • $\begingroup$ Also, what do you mean by your last paragraph? Clearly, there is a difference between in-sample and out-of-sample model validation (just like in cross-sectional setting). $\endgroup$ – Richard Hardy Feb 13 '17 at 17:52
  • $\begingroup$ @RichardHardy there is a difference, but we blur it when using cross validation the way you describe to a point where it become almost indistinguishable in terms of power for model selection. the ustility of out of sample is in dealing with overfitting, and in your procedure it doesn't help much in this regard $\endgroup$ – Aksakal Feb 13 '17 at 18:15
  • $\begingroup$ OK, in the sense that leave-one-out cross validation and AIC are asymptotically equivalent in terms of model selection, unless (if I am not mistaken) we misspecify the likelihood which causes AIC to misbehave. But I maintain that there is a clear difference in model evaluation when doing out-of-sample vs. in-sample forecast comparison. In-sample comparison can be way too optimistic while out-of-sample (using rolling windows) will be much more realistic. $\endgroup$ – Richard Hardy Feb 13 '17 at 18:19
  • $\begingroup$ If you deal with the business data, e.g. retail revenues, then it's easy to see how time separated holdout samples would be an issue because there could be changes in the business between sub samples. $\endgroup$ – Aksakal Feb 13 '17 at 18:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.