--- title: "Introduction to arimasel: Cartesian Product ARIMA Model Selection" author: "Olushina Olawale Awe (PStat, CStat, FRSS)" date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Introduction to arimasel} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 5, fig.align = "center" ) library(arimasel) set.seed(2026) ``` ## 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: 1. Accepts user-defined index sets $P$, $D$, and $Q$. 2. Computes the Cartesian product $P \times D \times Q$. 3. Fits every ARIMA$(p,d,q)$ triple exhaustively. 4. Ranks all converged models by four information criteria simultaneously: AIC, AICc, BIC, and HQIC. 5. Quantifies model uncertainty via Akaike weights. 6. 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()`. ```{r ts-eda, fig.height=6} eda <- ts_eda(gdp_ng) print(eda) ``` `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). ```{r ts-features} ts_features(inflation_ng) ``` ## 3. Enumerating Candidates with `cp_sets()` ```{r cp-sets} cp_sets(p_set = 0:2, d_set = 0:1, q_set = 0: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. ```{r data-plots, fig.height=4} 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 ```{r stationarity} stationarity_test(gdp_ng) d_rec <- suggest_d(gdp_ng) cat("Recommended d:", d_rec, "\n") ``` ## 6. Cartesian Product ARIMA Selection ```{r cart-arima} result <- cart_arima(gdp_ng, p_set = 0:3, d_set = 0:1, q_set = 0:3, criterion = "AIC", top_n = 10L) print(result) ``` ```{r summary} summary(result) ``` ## 7. Re-ranking by Different Criteria ```{r rerank} arima_table(result, criterion = "BIC", top_n = 5) ``` ## 8. Akaike Weights ```{r weights} ic_vals <- setNames(result$full_table$AIC, result$full_table$Model) arima_weights(ic_vals[1:6]) ``` ## 9. Criterion Vote Table ```{r vote} result$vote_table ``` ## 10. Visualising the Results ```{r plot-criteria, fig.height=5} plot(result, type = "criteria", top_n = 8) ``` ```{r plot-weights, fig.height=5} plot(result, type = "weights", top_n = 8) ``` ```{r plot-surface, fig.height=5} plot(result, type = "surface", criterion = "AIC") ``` ```{r plot-fitted, fig.height=5} plot(result, type = "fitted") ``` ## 11. Residual Diagnostics ```{r diagnose} arima_diagnose(result, lags = 15) ``` ## 12. Forecasting ### 11.1 Best-model forecast ```{r forecast-best, fig.height=5} fc <- arima_forecast(result, h = 8, level = c(80, 95), plot = TRUE) cat("Point forecasts:\n") round(fc$mean, 2) ``` ### 11.2 Ensemble forecast ```{r forecast-ensemble, fig.height=5} fc_ens <- arima_forecast(result, h = 8, ensemble = TRUE, top_k = 4L, plot = TRUE) ``` ## 13. Monthly Series: Exchange Rate ```{r exchange} res_ex <- cart_arima(exchange_ng, p_set = 0:2, d_set = 1L, q_set = 0:2, criterion = "BIC") print(res_ex) arima_forecast(res_ex, h = 12, plot = TRUE) ``` ## 14. Comparing with auto.arima() ```{r compare, eval=requireNamespace("forecast", quietly=TRUE)} compare_arima(gdp_ng, p_set = 0:3, d_set = 0:1, q_set = 0:3, holdout = 5L) ``` ## 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. ```{r seasonal-strength} seasonal_strength(inflation_ng)$seasonal_strength suggest_D(inflation_ng) ``` `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$: ```{r seasonal-search} 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) ``` ```{r plot-seasonal, fig.height=5} 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 = )`: ```{r xreg-demo} 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) ``` ## 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. ```{r cv-demo} cv <- arima_cv(result, h = 3, initial = floor(0.7 * result$n_obs)) print(cv) ``` ```{r plot-cv, fig.height=5} 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. ```{r smart-arima} res_smart <- smart_arima(inflation_ng, p_set = 0:1, q_set = 0:1) print(res_smart) ``` ## 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.