Introduction to arimasel: Cartesian Product ARIMA Model Selection

1. Motivation

The Box-Jenkins methodology remains the dominant framework for univariate time series modelling. Its classical identification step relies on the Autocorrelation Function (ACF) and Partial Autocorrelation Function (PACF) plots. However, these plots often present conflicting or ambiguous pictures, particularly for mixed ARIMA(p,d,q) models, making reliable identification difficult in practice.

The arimasel package proposes an alternative algorithm grounded in the mathematical theory of Cartesian products of sets. Rather than examining ACF/PACF plots or using stepwise search (as in auto.arima()), the algorithm:

Accepts user-defined index sets \(P\), \(D\), and \(Q\).
Computes the Cartesian product \(P \times D \times Q\).
Fits every ARIMA\((p,d,q)\) triple exhaustively.
Ranks all converged models by four information criteria simultaneously: AIC, AICc, BIC, and HQIC.
Quantifies model uncertainty via Akaike weights.
Supports best-model and ensemble forecasting.

2. Exploratory Data Analysis with `ts_eda()`

Before searching over any Cartesian product, it is good practice to look at the series itself. ts_eda() produces a time plot, distribution, ACF/PACF, and (for seasonal series) a seasonal subseries plot and seasonal lag plot, alongside a console summary of scale-free time series characteristics computed by ts_features().

eda <- ts_eda(gdp_ng)
#> Registered S3 method overwritten by 'quantmod':
#>   method            from
#>   as.zoo.data.frame zoo
#> Warning in tseries::kpss.test(x_vec): p-value greater than printed p-value

print(eda)
#> 
#> Exploratory Data Analysis
#> ==========================
#> n = 34  |  frequency = 1 (non-seasonal)
#> mean = 4.2176  |  sd = 2.8260  |  skewness = -0.0340  |  kurtosis = -0.4341
#> spectral entropy = 0.7273  |  ACF(1) = 0.5921
#> ADF statistic = -1.8808  (p ~ 0.1795)
#> Suggested d = 2

ts_features() can also be used on its own – for example, to compare several series’ characteristics side by side, or as an input to automated model recommendation, as in smart_arima() (Section 15).

ts_features(inflation_ng)
#>       n    mean     sd skewness kurtosis trend_strength seasonal_strength
#> xt1 288 14.1957 3.0576   0.0508  -0.6994         0.7901            0.0381
#>     entropy   acf1 lumpiness stability adf_stat adf_pvalue
#> xt1  0.6273 0.8691   28.5642    3.0594  -3.9536     0.0128

3. Enumerating Candidates with `cp_sets()`

cp_sets(p_set = 0:2, d_set = 0:1, q_set = 0:2)
#> Cartesian Product  P x D x Q
#>   P = {0, 1, 2}
#>   D = {0, 1}
#>   Q = {0, 1, 2}
#>   Total candidate model(s): 18
#> 
#> Candidate ARIMA models:
#> ARIMA(0,0,0)  ARIMA(0,0,1)  ARIMA(0,0,2)  ARIMA(1,0,0)  ARIMA(1,0,1)  ARIMA(1,0,2)  ARIMA(2,0,0)  ARIMA(2,0,1)  ARIMA(2,0,2)  ARIMA(0,1,0)  ARIMA(0,1,1)  ARIMA(0,1,2)  ARIMA(1,1,0)  ARIMA(1,1,1)  ARIMA(1,1,2)  ARIMA(2,1,0)  ARIMA(2,1,1)  ARIMA(2,1,2)

The function confirms \(|P| \times |D| \times |Q| = 3 \times 2 \times 3 = 18\) candidate models.

4. Data

The package ships with three Nigerian macroeconomic time series.

data(gdp_ng)
data(exchange_ng)
data(inflation_ng)

oldpar <- par(mfrow = c(1, 3), mar = c(4, 4, 3, 1))
plot(gdp_ng,      main = "GDP Growth (%)",    ylab = "%")
plot(exchange_ng, main = "USD/NGN Rate",      ylab = "NGN")
plot(inflation_ng,main = "CPI Inflation (%)", ylab = "%")

par(oldpar)

5. Pre-screening: Stationarity and Differencing

stationarity_test(gdp_ng)
#> Warning in tseries::kpss.test(x_vec): p-value greater than printed p-value
#> 
#> Stationarity Test Battery  (alpha = 0.05)
#> n = 34 observations
#> 
#>                   Test Statistic p_value     Conclusion
#>    ADF (H0: unit root)   -1.8808  0.1795 Non-stationary
#>     PP (H0: unit root)  -15.0002  0.1639 Non-stationary
#>  KPSS (H0: stationary)    0.1995  0.1000     Stationary
#> 
#> Consensus : Non-stationary  |  Suggested d : 1
d_rec <- suggest_d(gdp_ng)
cat("Recommended d:", d_rec, "\n")
#> Recommended d: 2

6. Cartesian Product ARIMA Selection

result <- cart_arima(gdp_ng,
                     p_set     = 0:3,
                     d_set     = 0:1,
                     q_set     = 0:3,
                     criterion = "AIC",
                     top_n     = 10L)
print(result)
#> 
#> Cartesian Product ARIMA Model Selection
#> =========================================
#> Call      : cart_arima(x = gdp_ng, p_set = 0:3, d_set = 0:1, q_set = 0:3, 
#>  Call      :     criterion = "AIC", top_n = 10L)
#> Obs (n)   : 34
#> Candidates: 32  |  Converged: 32  |  Failed: 0
#> Criterion : AIC
#> Best model: ARIMA(0,1,0)
#> 
#> Top 10 models (ranked by AIC):
#>  Rank        Model AIC AICc BIC HQIC Delta Weight
#>     1 ARIMA(0,1,0) 155  155 156  155 0.000 0.1381
#>     2 ARIMA(1,1,1) 155  156 159  156 0.053 0.1345
#>     3 ARIMA(0,1,1) 155  156 158  156 0.576 0.1036
#>     4 ARIMA(0,1,2) 156  156 160  157 1.012 0.0833
#>     5 ARIMA(1,1,0) 156  156 159  157 1.267 0.0733
#>     6 ARIMA(2,1,0) 156  157 161  158 1.381 0.0692
#>     7 ARIMA(1,0,0) 157  157 161  158 1.899 0.0534
#>     8 ARIMA(2,1,1) 157  158 163  159 2.052 0.0495
#>     9 ARIMA(1,1,2) 157  158 163  159 2.052 0.0495
#>    10 ARIMA(0,1,3) 158  159 164  160 2.998 0.0309

summary(result)
#> 
#> === arimasel: Cartesian Product ARIMA Summary ===
#> 
#> Call: cart_arima(x = gdp_ng, p_set = 0:3, d_set = 0:1, q_set = 0:3, 
#> 
#>  Call:     criterion = "AIC", top_n = 10L)
#> 
#> Search space:
#>   P = {0,1,2,3}, D = {0,1}, Q = {0,1,2,3}
#>   Total candidates: 32  |  Converged: 32
#> 
#> Best model by primary criterion:
#>   ARIMA(0,1,0)  (criterion: AIC)
#> 
#> Criterion vote table:
#>         Model Votes          Criteria
#>  ARIMA(0,1,0)     4 AIC/AICc/BIC/HQIC
#> 
#> Top 10 model comparison:
#>         Model p d q  P  D  Q LogLik AIC AICc BIC HQIC Rank Delta Weight
#>  ARIMA(0,1,0) 0 1 0 NA NA NA  -76.3 155  155 156  155    1 0.000 0.1381
#>  ARIMA(1,1,1) 1 1 1 NA NA NA  -74.4 155  156 159  156    2 0.053 0.1345
#>  ARIMA(0,1,1) 0 1 1 NA NA NA  -75.6 155  156 158  156    3 0.576 0.1036
#>  ARIMA(0,1,2) 0 1 2 NA NA NA  -74.8 156  156 160  157    4 1.012 0.0833
#>  ARIMA(1,1,0) 1 1 0 NA NA NA  -76.0 156  156 159  157    5 1.267 0.0733
#>  ARIMA(2,1,0) 2 1 0 NA NA NA  -75.0 156  157 161  158    6 1.381 0.0692
#>  ARIMA(1,0,0) 1 0 0 NA NA NA  -75.3 157  157 161  158    7 1.899 0.0534
#>  ARIMA(2,1,1) 2 1 1 NA NA NA  -74.4 157  158 163  159    8 2.052 0.0495
#>  ARIMA(1,1,2) 1 1 2 NA NA NA  -74.4 157  158 163  159    9 2.052 0.0495
#>  ARIMA(0,1,3) 0 1 3 NA NA NA  -74.8 158  159 164  160   10 2.998 0.0309
#> 
#> Best model coefficients:
#> numeric(0)
#> 
#> sigma^2: 5.9821  |  Log-likelihood: -76.340

7. Re-ranking by Different Criteria

arima_table(result, criterion = "BIC", top_n = 5)
#>          Model p d q  P  D  Q  LogLik     AIC    AICc     BIC    HQIC Rank
#> 1 ARIMA(0,1,0) 0 1 0 NA NA NA -76.340 154.680 154.805 156.206 155.200    1
#> 2 ARIMA(0,1,1) 0 1 1 NA NA NA -75.628 155.256 155.643 158.309 156.297    2
#> 3 ARIMA(1,1,0) 1 1 0 NA NA NA -75.974 155.947 156.334 159.000 156.988    3
#> 4 ARIMA(1,1,1) 1 1 1 NA NA NA -74.366 154.733 155.533 159.312 156.294    4
#> 5 ARIMA(0,1,2) 0 1 2 NA NA NA -74.846 155.692 156.492 160.271 157.253    5
#>   Delta Weight
#> 1 0.000 0.4281
#> 2 2.103 0.1496
#> 3 2.794 0.1059
#> 4 3.106 0.0906
#> 5 4.065 0.0561

8. Akaike Weights

ic_vals <- setNames(result$full_table$AIC,
                    result$full_table$Model)
arima_weights(ic_vals[1:6])
#>          Model IC_value Delta Weight EvidenceRatio
#> 1 ARIMA(0,1,0)  154.680 0.000 0.2294          1.00
#> 2 ARIMA(1,1,1)  154.733 0.053 0.2234          1.03
#> 3 ARIMA(0,1,1)  155.256 0.576 0.1720          1.33
#> 4 ARIMA(0,1,2)  155.692 1.012 0.1383          1.66
#> 5 ARIMA(1,1,0)  155.947 1.267 0.1218          1.88
#> 6 ARIMA(2,1,0)  156.061 1.381 0.1150          1.99

9. Criterion Vote Table

result$vote_table
#>          Model Votes          Criteria
#> 1 ARIMA(0,1,0)     4 AIC/AICc/BIC/HQIC

10. Visualising the Results

plot(result, type = "criteria", top_n = 8)

plot(result, type = "weights", top_n = 8)

plot(result, type = "surface", criterion = "AIC")

plot(result, type = "fitted")

11. Residual Diagnostics

arima_diagnose(result, lags = 15)

#> 
#> Residual diagnostics for ARIMA(0,1,0)
#> n (residuals) = 34
#> sigma^2 = 5.9821
#> Shapiro-Wilk: W = 0.9839, p-value = 0.8859  [Cannot reject normality]
#> Ljung-Box: No significant autocorrelation detected (alpha = 0.05).

12. Forecasting

11.1 Best-model forecast

fc <- arima_forecast(result, h = 8, level = c(80, 95), plot = TRUE)

cat("Point forecasts:\n")
#> Point forecasts:
round(fc$mean, 2)
#> [1] 2.9 2.9 2.9 2.9 2.9 2.9 2.9 2.9

11.2 Ensemble forecast

fc_ens <- arima_forecast(result, h = 8, ensemble = TRUE,
                         top_k = 4L, plot = TRUE)

13. Monthly Series: Exchange Rate

res_ex <- cart_arima(exchange_ng,
                     p_set = 0:2, d_set = 1L, q_set = 0:2,
                     criterion = "BIC")
print(res_ex)
#> 
#> Cartesian Product ARIMA Model Selection
#> =========================================
#> Call      : cart_arima(x = exchange_ng, p_set = 0:2, d_set = 1L, q_set = 0:2, 
#>  Call      :     criterion = "BIC")
#> Obs (n)   : 168
#> Candidates: 9  |  Converged: 9  |  Failed: 0
#> Criterion : BIC
#> Best model: ARIMA(0,1,0)
#> 
#> Top 9 models (ranked by BIC):
#>  Rank        Model  AIC AICc  BIC HQIC Delta Weight
#>     1 ARIMA(0,1,0) 1624 1624 1627 1625  0.00 0.8478
#>     2 ARIMA(0,1,1) 1626 1626 1632 1628  5.10 0.0661
#>     3 ARIMA(1,1,0) 1626 1626 1632 1628  5.11 0.0660
#>     4 ARIMA(0,1,2) 1627 1627 1637 1631  9.67 0.0068
#>     5 ARIMA(2,1,0) 1627 1627 1637 1631  9.68 0.0067
#>     6 ARIMA(1,1,1) 1628 1628 1637 1631 10.03 0.0056
#>     7 ARIMA(1,1,2) 1629 1629 1642 1634 14.79 0.0005
#>     8 ARIMA(2,1,1) 1629 1629 1642 1634 14.80 0.0005
#>     9 ARIMA(2,1,2) 1631 1632 1647 1638 19.91 0.0000
arima_forecast(res_ex, h = 12, plot = TRUE)

14. Comparing with auto.arima()

compare_arima(gdp_ng, p_set = 0:3, d_set = 0:1, q_set = 0:3,
              holdout = 5L)
#> Fitting cart_arima...
#> Fitting auto.arima()...
#> 
#> --- Model Comparison ---
#>      Method        Model AIC AICc BIC HQIC OOS_RMSE OOS_MAE
#>  cart_arima ARIMA(0,1,0) 129  129 131  130     2.07    1.84
#>  auto.arima ARIMA(1,0,0) 133  134 137  134     2.24    1.53
#> 
#> Winner (most metric wins): cart_arima

15. Seasonal ARIMA Search

Monthly series such as inflation_ng often display genuine seasonal structure. seasonal_strength() quantifies this (via STL decomposition) before committing to a seasonal search, and suggest_D() turns that into a recommended seasonal differencing order.

seasonal_strength(inflation_ng)$seasonal_strength
#> [1] 0.03810152
suggest_D(inflation_ng)
#> [1] 0

cart_arima() accepts a seasonal specification – a list with index sets P, D, Q and the seasonal period – and exhaustively searches every combination of \((p,d,q)\) and \((P,D,Q)_m\):

res_seas <- cart_arima(inflation_ng,
                       p_set = 0:1, d_set = 0:1, q_set = 0:1,
                       seasonal = list(P = 0:1, D = 0:1, Q = 0:1, period = 12),
                       criterion = "AIC", top_n = 8L)
print(res_seas)
#> 
#> Cartesian Product ARIMA Model Selection
#> =========================================
#> Call      : cart_arima(x = inflation_ng, p_set = 0:1, d_set = 0:1, q_set = 0:1, 
#>  Call      :     seasonal = list(P = 0:1, D = 0:1, Q = 0:1, period = 12), 
#>  Call      :     criterion = "AIC", top_n = 8L)
#> Obs (n)   : 288
#> Seasonal  : period = 12, P = {0,1}, D = {0,1}, Q = {0,1}
#> Candidates: 64  |  Converged: 64  |  Failed: 0
#> Criterion : AIC
#> Best model: ARIMA(1,0,0)(0,1,1)[12]
#> 
#> Top 8 models (ranked by AIC):
#>  Rank                   Model  AIC AICc  BIC HQIC Delta Weight
#>     1 ARIMA(1,0,0)(0,1,1)[12] 1055 1055 1066 1059 0.000 0.2279
#>     2 ARIMA(1,0,0)(1,1,1)[12] 1055 1055 1070 1061 0.371 0.1893
#>     3 ARIMA(1,0,1)(1,1,1)[12] 1056 1056 1075 1064 1.661 0.0993
#>     4 ARIMA(1,0,1)(0,1,1)[12] 1056 1056 1071 1062 1.705 0.0972
#>     5 ARIMA(1,1,1)(0,1,1)[12] 1056 1057 1071 1062 1.888 0.0887
#>     6 ARIMA(1,1,1)(1,1,1)[12] 1057 1057 1075 1064 2.116 0.0791
#>     7 ARIMA(1,1,1)(0,0,0)[12] 1058 1058 1069 1062 3.315 0.0434
#>     8 ARIMA(1,1,1)(0,0,1)[12] 1058 1058 1073 1064 3.676 0.0363

plot(res_seas, type = "seasonal")

16. Exogenous Regressors (Regression with ARIMA Errors)

cart_arima() also accepts an xreg matrix of external predictors, fitted alongside every candidate model, with matching support in arima_forecast(..., newxreg = ):

set.seed(42)
trend <- as.numeric(time(gdp_ng)) - 1990
res_xreg <- cart_arima(gdp_ng, p_set = 0:2, d_set = 0:1, q_set = 0:2,
                       xreg = trend, criterion = "AIC")
fc_xreg <- arima_forecast(res_xreg, h = 5,
                          newxreg = matrix(max(trend) + 1:5, ncol = 1),
                          plot = FALSE)
round(fc_xreg$mean, 2)
#> [1] 2.85 2.79 2.72 2.65 2.57

17. Rolling-Origin Cross-Validation

Information criteria are in-sample measures. arima_cv() complements them with rolling-origin (expanding-window) cross-validation – refitting the selected model order at successive origins and recording genuine out-of-sample forecast errors, a standard “backtesting” practice in applied and industrial forecasting.

cv <- arima_cv(result, h = 3, initial = floor(0.7 * result$n_obs))
print(cv)
#> 
#> Rolling-Origin Cross-Validation
#> ================================
#> Model order : ARIMA(0,1,0)
#> Origins     : 9  (initial = 23, step = 1, h = 3)
#> 
#> Out-of-sample accuracy by horizon:
#>  Horizon   RMSE    MAE     MAPE N
#>        1 3.0428 2.5222 128.4525 9
#>        2 3.7477 3.1000 156.0349 9
#>        3 3.6220 2.7889 171.0843 9

plot(cv, type = "rmse")

18. Feature-Guided Automatic Search with `smart_arima()`

smart_arima() combines the exploratory feature analysis of Sections 2 and 15 with cart_arima(): it computes time series features, uses them to narrow the differencing order (and, for seasonal series, to decide whether a seasonal search is worthwhile and what seasonal differencing order to use), and then runs the usual exhaustive search within that narrowed space. It is a convenient default entry point when you do not want to specify d_set (or seasonal) by hand, while remaining fully transparent about the choices it makes.

res_smart <- smart_arima(inflation_ng, p_set = 0:1, q_set = 0:1)
#> smart_arima: feature-guided search space
#>   seasonal strength = 0.038 | trend strength = 0.790 | entropy = 0.627
#>   d_set = {0,1}  |  seasonal D = {0}, period = 12
print(res_smart)
#> 
#> Cartesian Product ARIMA Model Selection
#> =========================================
#> Call      : cart_arima(x = x_ts, p_set = p_set, d_set = d_set, q_set = q_set, 
#>  Call      :     seasonal = seasonal_arg, criterion = criterion, top_n = top_n, 
#>  Call      :     parallel = parallel, n_cores = n_cores)
#> Obs (n)   : 288
#> Seasonal  : period = 12, P = {0,1}, D = {0}, Q = {0,1}
#> Candidates: 32  |  Converged: 32  |  Failed: 0
#> Criterion : AIC
#> Best model: ARIMA(1,1,1)(0,0,0)[12]
#> 
#> Top 10 models (ranked by AIC):
#>  Rank                   Model  AIC AICc  BIC HQIC Delta Weight
#>     1 ARIMA(1,1,1)(0,0,0)[12] 1058 1058 1069 1062 0.000 0.2047
#>     2 ARIMA(1,1,1)(0,0,1)[12] 1058 1058 1073 1064 0.361 0.1709
#>     3 ARIMA(1,1,1)(1,0,0)[12] 1058 1058 1073 1064 0.475 0.1614
#>     4 ARIMA(1,0,0)(0,0,0)[12] 1059 1059 1070 1064 1.490 0.0972
#>     5 ARIMA(1,0,0)(0,0,1)[12] 1060 1060 1074 1066 1.938 0.0777
#>     6 ARIMA(1,0,0)(1,0,0)[12] 1060 1060 1075 1066 2.042 0.0737
#>     7 ARIMA(1,1,1)(1,0,1)[12] 1060 1060 1078 1067 2.239 0.0668
#>     8 ARIMA(1,0,1)(0,0,0)[12] 1061 1061 1076 1067 3.425 0.0369
#>     9 ARIMA(1,0,1)(0,0,1)[12] 1062 1062 1080 1069 3.693 0.0323
#>    10 ARIMA(1,0,1)(1,0,0)[12] 1062 1062 1080 1069 3.802 0.0306

19. Summary

The arimasel package provides a transparent, exhaustive, multi-criteria approach to ARIMA model identification. Key advantages over stepwise methods:

Exhaustive: every \((p,d,q)\) in the Cartesian product is evaluated.
Multi-criteria: simultaneous ranking by AIC, AICc, BIC, and HQIC.
Uncertainty quantification: Akaike weights and evidence ratios.
Ensemble forecasting: weighted averaging across top-\(k\) models.
Seasonal-aware: exhaustive search over \((p,d,q)(P,D,Q)_m\), with STL-based seasonal-strength diagnostics.
Regression-capable: exogenous regressors via xreg/newxreg.
Out-of-sample validated: rolling-origin cross-validation via arima_cv().
Feature-driven: exploratory ts_eda()/ts_features() and a feature-guided smart_arima() search.
Reproducible: fully deterministic given the input sets.

References

Awe, O. O. (2026). arimasel: Cartesian Product-Based ARIMA Model Identification and Selection. R package version 0.2.0.

Akaike, H. (1974). A new look at the statistical model identification. IEEE Transactions on Automatic Control, 19(6), 716–723.

Box, G. E. P., Jenkins, G. M., Reinsel, G. C., and Ljung, G. M. (2015). Time Series Analysis: Forecasting and Control (5th ed.). John Wiley & Sons.

Burnham, K. P. and Anderson, D. R. (2002). Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach (2nd ed.). Springer.

Hannan, E. J. and Quinn, B. G. (1979). The determination of the order of an autoregression. Journal of the Royal Statistical Society, Series B, 41(2), 190–195.

Hyndman, R. J. and Athanasopoulos, G. (2021). Forecasting: Principles and Practice (3rd ed.). OTexts.

Hyndman, R. J., Wang, E. and Laptev, N. (2015). Large-scale unusual time series detection. 2015 IEEE International Conference on Data Mining Workshop, 1616–1623.

Wang, X., Smith, K. A. and Hyndman, R. J. (2006). Characteristic-based clustering for time series data. Data Mining and Knowledge Discovery, 13(3), 335–364.

Introduction to arimasel: Cartesian Product ARIMA Model Selection

Olushina Olawale Awe (PStat, CStat, FRSS)

2026-09-09

1. Motivation

2. Exploratory Data Analysis with `ts_eda()`

3. Enumerating Candidates with `cp_sets()`

4. Data

5. Pre-screening: Stationarity and Differencing

6. Cartesian Product ARIMA Selection

7. Re-ranking by Different Criteria

8. Akaike Weights

9. Criterion Vote Table

10. Visualising the Results

11. Residual Diagnostics

12. Forecasting

11.1 Best-model forecast

11.2 Ensemble forecast

13. Monthly Series: Exchange Rate

14. Comparing with auto.arima()

15. Seasonal ARIMA Search

16. Exogenous Regressors (Regression with ARIMA Errors)

17. Rolling-Origin Cross-Validation

18. Feature-Guided Automatic Search with `smart_arima()`

19. Summary

References

Introduction to arimasel: Cartesian Product ARIMA Model Selection

Olushina Olawale Awe (PStat, CStat, FRSS)

2026-09-09

1. Motivation

2. Exploratory Data Analysis with ts_eda()

3. Enumerating Candidates with cp_sets()

4. Data

5. Pre-screening: Stationarity and Differencing

6. Cartesian Product ARIMA Selection

7. Re-ranking by Different Criteria

8. Akaike Weights

9. Criterion Vote Table

10. Visualising the Results

11. Residual Diagnostics

12. Forecasting

11.1 Best-model forecast

11.2 Ensemble forecast

13. Monthly Series: Exchange Rate

14. Comparing with auto.arima()

15. Seasonal ARIMA Search

16. Exogenous Regressors (Regression with ARIMA Errors)

17. Rolling-Origin Cross-Validation

18. Feature-Guided Automatic Search with smart_arima()

19. Summary

References

2. Exploratory Data Analysis with `ts_eda()`

3. Enumerating Candidates with `cp_sets()`

18. Feature-Guided Automatic Search with `smart_arima()`