I have a code which tests each possible order of ARIMA and selects the best model by choosing the one with the absolute minimum sum of lags from the PACF graph. The code then proceeds to add weight to recent errors and runs an optimization on the parameters to get the minimum mean absolute error.

The code runs fine and gives excellent results (e.g 0.2% MAPE etc) however once the parameters have been optimized the ACF and PACf graphs show lags outside the threshold.

I would like to add into my code a loop which does the following:

if any of the first 4 lags of the ACF or PACF graphs of the residuals found from the optimized ARIMA model are outside the threshold (2/sqrt(n)) then the optimization is re-run but doesn't allow those parameters to be selected/those parameters are skipped in the optimization process.

Here is my code:

Data.col<-c(5403.676,6773.505, 7231.117, 7835.552, 5236.710, 5526.619, 6555.782,11464.727, 7210.069, 7501.610, 8670.903,10872.935, 8209.023, 8153.393,10196.448,13244.502, 8356.733,10188.442,10601.322,12617.821, 11786.526,10044.987,11006.005,15101.946,10992.273,11421.189,10731.312)
# Turns the Data.col into a time-series

Data.col.ts <- ts(Data.col, deltat=(1/4), start = c(8,1))

# Starts the testing to see if the data should be logged

trans<- BoxCox.lambda(Data.col, method = "loglik")
categ<-as.character( c(cut(trans,c(0,0.25,0.75,Inf),right=FALSE)) )

#----- Weighting ---------------------------------------------------------------
fweight <- function(x){
  PatX <- 0.5+x 

#Split the integral to several intervals, and pick the weights accordingly

integvals <- rep(0, length.out = length(Data.new))
for (i in 1:length(Data.new)){
  integi <- integrate(fweight, lower = (i-1)/length(Data.new), upper= i/length(Data.new))
  integvals[i] <- 2*integi$value

#----- Find best ARIMA model ---------------------------------------------------

a <- permutations(n = 3, r = 6, v = c(0:2), repeats.allowed = TRUE)
a <- a[ifelse((a[, 1] + a[, 4] > 2 | a[, 2] + a[, 5] > 2 | a[, 3] + a[, 6] > 2),
              FALSE, TRUE), ]

Arimafit <- matrix(0,
                   ncol  = length(Data.new),
                   nrow  = length(a[, 1]),
                   byrow = TRUE)

totb <- matrix(0, ncol = 1, nrow = length(a[, 1]))
arimaerror <- matrix(0, ncol = length(Data.new), nrow = 1)

for (i in 1:length(a[, 1])){
  ArimaData.new <- try(Arima(Data.new,
                             order    = a[i, c(1:3)],
                             seasonal = list(order = a[i, c(4:6)]),
                             method   = "ML"),
                       silent = TRUE)

  if (is(ArimaData.new, "try-error")){
    ArimaData.new <- arimaerror
  } else {
    ArimaData.new <- ArimaData.new

  arimafitted <- try(fitted(ArimaData.new), silent = TRUE)

  if (is(arimafitted, "try-error")){
    fitarima <- arimaerror
  } else {
    fitarima <- arimafitted

  if (categ=="1"){
    Arimafit[i, ] <- c(exp(fitarima))
    Datanew <- c(exp(Data.new))
  } else {
    if (categ=="2"){
      Arimafit[i, ] <- c((fitarima)^2)
      Datanew <- c((Data.new)^2)
    } else {
      Arimafit[i, ] <- c(fitarima)
      Datanew <- c(Data.new)

  data <- c(Datanew)

  arima.fits <- c(Arimafit[i, ])

  fullres <- data - arima.fits

  v <- acf(fullres, plot = FALSE)

  w <- pacf(fullres, plot = FALSE)

  if (v$acf[2]>(2/sqrt(length(Data.col)))|v$acf[2]<(-(2/sqrt(length(Data.col))))|v$acf[3]>(2/sqrt(length(Data.col)))|v$acf[3]<(-(2/sqrt(length(Data.col))))|v$acf[4]>(2/sqrt(length(Data.col)))|v$acf[4]<(-(2/sqrt(length(Data.col))))|v$acf[5]>(2/sqrt(length(Data.col)))|v$acf[5]<(-(2/sqrt(length(Data.col))))|v$acf[6]>(2/sqrt(length(Data.col)))|v$acf[6]<(-(2/sqrt(length(Data.col))))|v$acf[7]>(2/sqrt(length(Data.col)))|v$acf[7]<(-(2/sqrt(length(Data.col))))|w$acf[1]>(2/sqrt(length(Data.col)))|w$acf[1]<(-(2/sqrt(length(Data.col))))|w$acf[2]>(2/sqrt(length(Data.col)))|w$acf[2]<(-(2/sqrt(length(Data.col))))|w$acf[3]>(2/sqrt(length(Data.col)))|w$acf[3]<(-(2/sqrt(length(Data.col))))|w$acf[4]>(2/sqrt(length(Data.col)))|w$acf[4]<(-(2/sqrt(length(Data.col))))|w$acf[5]>(2/sqrt(length(Data.col)))|w$acf[5]<(-(2/sqrt(length(Data.col))))|w$acf[6]>(2/sqrt(length(Data.col)))|w$acf[6]<(-(2/sqrt(length(Data.col))))){
    totb[i] <- "n"
  } else {
    totb[i] <- sum(abs(w$acf[1:4]))

  j <- match(min(totb), totb)

  order.arima <- a[j, c(1:3)]

  order.seasonal.arima <- a[j, c(4:6)]

#----- ARIMA -------------------------------------------------------------------
# Fits an ARIMA model with the orders set
stAW <- Arima(Data.new, order= order.arima, seasonal=list(order=order.seasonal.arima), method="ML")
parSW <- stAW$coef
WMAEOPT <- function(parSW)
  ArimaW <- Arima(Data.new, order = order.arima, seasonal=list(order=order.seasonal.arima), 
                  include.drift=FALSE, method = "ML", fixed = c(parSW))
  errAR <- c(abs(resid(ArimaW)))
  WMAE <- t(errAR) %*% integvals 
OPTWMAE <- optim(parSW, WMAEOPT, method="SANN", set.seed(2), control = list(fnscale= 1, maxit = 5000))
# Alternatively, set  method="Nelder-Mead" or method="L-BFGS-B" 
parS3 <- OPTWMAE$par
Arima.Data.new <- Arima(Data.new, order = order.arima, seasonal=list(order=order.seasonal.arima), 
                        include.drift=FALSE, method = "ML", fixed = c(parS3))

Before the parameters are optimized it gives a graph like this: enter image description here

After the parameters are optimized it gives a graph like this: enter image description here

I want to stop this happening in the second picture. Is this possible to do using optim?

  • $\begingroup$ Not an answer to you question but... Have you considered using AIC (or AICc) to select the best ARIMA model from the pool of all candidate models? This approach has been advocated by Rob J. Hyndman (e.g. this blog post, last paragraph; and related posts), and this may be a good enough reason to seriously consider it. $\endgroup$ – Richard Hardy Nov 17 '14 at 15:08
  • $\begingroup$ hi yes i have many forecasting models. I found that despite AIC being a good and approved way of choosing the ARIMA order, it still gives correlation. My code finds the order with the least amount of correlation and i have found that to give 'better' results than using auto.arima $\endgroup$ – Summer-Jade Gleek'away Nov 17 '14 at 15:11
  • $\begingroup$ Depending on your loss function, you might prioritize the least amount of correlation over what AIC (or AICc) has to offer; but I suggest you think carefully what your loss function actually is. For practical purposes, it might be more plausible that one cares about accurate point forecasts rather than about minimizing the amount of correlations. Of course, you know better than me what you are aiming at :) $\endgroup$ – Richard Hardy Nov 17 '14 at 15:42
  • $\begingroup$ i am writing lots of different codes which all have a slightly different approach and then creating one big code that combines all of them and examines the different statistical outputs from each and 'scores' them and then will produce a document with the outputs and forecasts of the top 3/4. the codes for the 'normal' forecasting methods are already completed and now im just experimenting with new adapted versions $\endgroup$ – Summer-Jade Gleek'away Nov 17 '14 at 15:46

Your Answer

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

Browse other questions tagged or ask your own question.