Detecting patterns in white noise

Back in 2015, when I was working on my paper on Complex Exponential Smoothing, I conducted a simple simulation experiment to check how ARIMA and ETS select components/orders in time series. And I found something interesting…

One of the important steps in forecasting with statistical models is identifying the existing structure. In the case of ETS, it comes to selecting trend/seasonal components, while for ARIMA, it’s about order selection. In R, several functions automatically handle this based on information criteria (Hyndman & Khandakar, 2006; Svetunkov & Boylan (2017); Chapter 15 of ADAM). I decided to investigate how this mechanism works.

I generated data from the Normal distribution with a fixed mean of 5000 and a standard deviation of 50. Then, I asked ETS and ARIMA (from the forecast package in R) to automatically select the appropriate model for each of 1000 time series. Here is the R code for this simple experiment:

Some R code
# Set random seed for reproducibility
set.seed(41, kind="L'Ecuyer-CMRG")
# Number of iterations
nsim <- 1000
# Number of observations
obsAll <- 120
# Generate data from N(5000, 50)
rnorm(nsim*obsAll, 5000, 50) |>
  matrix(obsAll, nsim) |>
  ts(frequency=12) -> x

# Load forecast package
library(forecast)
# Load doMC for parallel calculations
# doMC is only available on Linux and Max
# Use library(doParallel) on Windows
library(doMC)
registerDoMC(detectCores())

# A loop for ARIMA, recording the orders
matArima <- foreach(i=1:nsim, .combine=cbind, .packages=c("forecast")) %dopar% {
    testModel <- auto.arima(x[,i])
    # The element number 5 is just m, period of seasonality
    return(c(testModel$arma[-5],(!is.na(testModel$coef["drift"]))*1))
}
rownames(matArima) <- c("AR","MA","SAR","SMA","I","SI","Drift")

# A loop for ETS, recording the model types
matEts <- foreach(i=1:nsim, .combine=cbind, .packages=c("forecast")) %dopar% {
    testModel <- ets(x[,i], allow.multiplicative.trend=TRUE)
    return(testModel[13]$method)
}

The findings of this experiment are summarised using the following chunk of the R code:

R code for the analysis of the results
#### Auto ARIMA ####
# Non-seasonal ARIMA elements
mean(apply(matArima[c("AR","MA","I","Drift"),]!=0, 2, any))
# Seasonal ARIMA elements
mean(apply(matArima[c("SAR","SMA","SI"),]!=0, 2, any))

#### ETS ####
# Trend in ETS
mean(substr(matEts,7,7)!="N")
# Seasonality in ETS
mean(substr(matEts,nchar(matEts)-1,nchar(matEts)-1)!="N")

I summarised them in the following table:

ARIMA ETS
Non-seasonal elements 24.8% 2.3%
Seasonal elements 18.0% 0.2%
Any type of structure 37.9% 2.4%

So, ARIMA detected some structure (had non-zero orders) in almost 40% of all time series, even though the data was designed to have no structure (just white noise). It also captured non-seasonal orders in a quarter of the series and identified seasonality in 18% of them. ETS performed better (only 0.2% of seasonal models identified on the white noise), but still captured trends in 2.3% of cases.

Does this simple experiment suggest that ARIMA is a bad model and ETS is a good one? No, it does not. It simply demonstrates that ARIMA tends to overfit the data if allowed to select whatever it wants. How can we fix that?

My solution: restrict the pool of ARIMA models to check, preventing it from going crazy. My personal pool includes ARIMA(0,1,1), (1,1,2), (0,2,2), along with the seasonal orders of (0,1,1), (1,1,2), and (0,2,2), and combinations between them. This approach is motivated by the connection between ARIMA and ETS. Additionally, we can check whether the addition of AR/MA orders detected by ACF/PACF analysis of the best model reduces the AICc. If not, they shouldn't be included.

This algorithm can be written in the following simple function that uses msarima() function from the smooth package in R (note that the reason why this function is used is because all ARIMA models implemented in the function are directly comparable via information criteria):

R code for the compact ARIMA function
arimaCompact <- function(y, lags=c(1,frequency(y)), ic=c("AICc","AIC","BIC","BICc"), ...){

    # Start measuring the time of calculations
    startTime <- Sys.time();

    # If there are no lags for the basic components, correct this.
    if(sum(lags==1)==0){
        lags <- c(1,lags);
    }

    orderLength <- length(lags);
    ic <- match.arg(ic);
    IC <- switch(ic,
                 "AIC"=AIC,
                 "AICc"=AICc,
                 "BIC"=BIC,
                 "BICc"=BICc);

    # We consider the following list of models:
    # ARIMA(0,1,1), (1,1,2), (0,2,2),
    # ARIMA(0,0,0)+c, ARIMA(0,1,1)+c,
    # seasonal orders (0,1,1), (1,1,2), (0,2,2)
    # And all combinations between seasonal and non-seasonal parts
    # 
    # Encode all non-seasonal parts
    nNonSeasonal <- 5
    arimaNonSeasonal <- matrix(c(0,1,1,0, 1,1,2,0, 0,2,2,0, 0,0,0,1, 0,1,1,1), nNonSeasonal,4,
                               dimnames=list(NULL, c("ar","i","ma","const")), byrow=TRUE)
    # Encode all seasonal parts ()
    nSeasonal <- 4
    arimaSeasonal <- matrix(c(0,0,0, 0,1,1, 1,1,2, 0,2,2), nSeasonal,3,
                               dimnames=list(NULL, c("sar","si","sma")), byrow=TRUE)

    # Check all the models in the pool
    testModels <- vector("list", nSeasonal*nNonSeasonal);
    m <- 1;
    for(i in 1:nSeasonal){
        for(j in 1:nNonSeasonal){
            testModels[[m]] <- msarima(y, orders=list(ar=c(arimaNonSeasonal[j,1],arimaSeasonal[i,1]),
                                                      i=c(arimaNonSeasonal[j,2],arimaSeasonal[i,2]),
                                                      ma=c(arimaNonSeasonal[j,3],arimaSeasonal[i,3])),
                                       constant=arimaNonSeasonal[j,4]==1, lags=lags, ...);
            m[] <- m+1;
        }
    }

    # Find the best one
    m <- which.min(sapply(testModels, IC));
    # Amend computational time
    testModels[[m]]$timeElapsed <- Sys.time()-startTime;

    return(testModels[[m]]);
}

Additionally, we can check whether the addition of AR/MA orders detected by ACF/PACF analysis of the best model reduces the AICc. If not, they shouldn't be included. I have not added that part in the code above. Still, this algorithm brings some improvements:

R code for the application of compact ARIMA to the data
#### Load the smooth package
library(smooth)

# A loop for the compact ARIMA, recording the orders
matArimaCompact <- foreach(i=1:nsim, .packages=c("smooth")) %dopar% {
    testModel <- arimaCompact(x[,i])
    return(orders(testModel))
}

#### Auto MSARIMA from smooth ####
# Non-seasonal ARIMA elements
mean(sapply(sapply(matArimaCompact, "[[", "ar"), function(x){x[1]!=0}) |
  sapply(sapply(matArimaCompact, "[[", "i"), function(x){x[1]!=0}) |
  sapply(sapply(matArimaCompact, "[[", "ma"), function(x){x[1]!=0}))

# Seasonal ARIMA elements
mean(sapply(sapply(matArimaSmooth, "[[", "ar"), function(x){length(x)==2 && (x[2]!=0)}) |
  sapply(sapply(matArimaSmooth, "[[", "i"), function(x){length(x)==2 && (x[2]!=0)}) |
  sapply(sapply(matArimaSmooth, "[[", "ma"), function(x){length(x)==2 && (x[2]!=0)}))

In my case, it resulted in the following:

ARIMA ETS Compact ARIMA
Non-seasonal elements 24.8% 2.3% 2.4%
Seasonal elements 18.0% 0.2% 0.0%
Any type of structure 37.9% 2.4% 2.4%

As we see, when we impose restrictions on order selection in ARIMA, it avoids fitting seasonal models to non-seasonal data. While it still makes minor mistakes in terms of non-seasonal structure, it's nothing compared to the conventional approach. What about accuracy? I don't know. I'll have to write another post on this :).

Note that the models were applied to samples of 120 observations, which is considered "small" in statistics, while in real life is sometimes a luxury to have...

Leave a Reply