How tsorder compares with stats::ar and auto.arima
Source:vignettes/comparisons.Rmd
comparisons.RmdThe natural question about a new order selector is how it behaves
next to the tools everyone already uses: stats::ar() (AIC
over Yule-Walker fits) and forecast::auto.arima() (stepwise
AICc search). This vignette runs a small, honest comparison on simulated
data where the true order is known. It is intentionally light — a
handful of replications per scenario so the vignette builds quickly; the
repository’s benchmark (benchmark/ in the sources) runs the
same kind of comparison with hundreds of replications per design and is
where the package’s defaults come from.
Two things are worth keeping in mind while reading:
- On a clean, well-separated AR process every reasonable method agrees; differences show up at short samples, near-cancellation, seasonality and contamination.
- tsorder’s threshold estimators answer a slightly different question than a criterion: not “which model predicts best” but “which lags are distinguishable from noise”. Under-and-over-selection are not symmetric for the two goals.
Setup
library(tsorder)
simulate_ar <- function(n, phi, seed) {
set.seed(seed)
as.numeric(arima.sim(list(ar = phi), n))
}
# One replication: the selected AR order for each method.
select_all <- function(x) {
c(
tsorder_pacf = order_select(x, "pacf")$order,
tsorder_jirak = order_select(x, "jirak_refit")$order,
tsorder_bic = ic_order(x, "bic")$order,
stats_ar = stats::ar(x, method = "yule-walker")$order
)
}A clear AR(2): everyone agrees
x <- simulate_ar(500, c(1.2, -0.6), seed = 1)
select_all(x)
#> tsorder_pacf tsorder_jirak tsorder_bic stats_ar
#> 2 10 2 2Where methods part ways: rate across replications
An AR(4) whose last coefficient is moderate, at a short sample.
stats::ar maximises fit with AIC, which is not consistent
and tends to overshoot; BIC and the threshold estimators pay a growing
price per lag.
reps <- 20
hits <- sapply(1:reps, function(r) {
x <- simulate_ar(300, c(0.5, 0, 0, 0.3), seed = 100 + r)
select_all(x) == 4L
})
rowMeans(hits)
#> tsorder_pacf tsorder_jirak tsorder_bic stats_ar
#> 1.00 1.00 0.95 0.65The exact rates move with the seed set — 20 replications is a smoke test, not an estimate — but the ordering is stable in the full benchmark: AIC-type rules over-select on short samples, the PACF-clip and BIC land on the truth or just below it.
Against auto.arima on ARMA
fit_arima() is tsorder’s end-to-end path: estimate
d, select the orders, hand everything to
stats::arima(). method = "hr" selects
p and q jointly by Hannan-Rissanen, which is
the closest analogue of what auto.arima does, at a fraction
of the fits.
set.seed(7)
y <- as.numeric(arima.sim(list(ar = 0.7, ma = 0.4), 400))
fit_ts <- fit_arima(y, method = "hr")
fit_aa <- forecast::auto.arima(y, seasonal = FALSE)
rbind(
tsorder = fit_ts$arma[c(1, 6, 2)],
auto.arima = fit_aa$arma[c(1, 6, 2)]
)
#> [,1] [,2] [,3]
#> tsorder 1 0 1
#> auto.arima 1 0 1
# Both are plain stats::arima fits, so forecast() applies to either.
h <- 12
f_ts <- forecast::forecast(fit_ts, h = h)
f_aa <- forecast::forecast(fit_aa, h = h)
round(cbind(tsorder = f_ts$mean, auto.arima = f_aa$mean)[1:6, ], 3)
#> tsorder auto.arima
#> [1,] -0.292 -0.326
#> [2,] -0.159 -0.232
#> [3,] -0.064 -0.165
#> [4,] 0.002 -0.117
#> [5,] 0.050 -0.083
#> [6,] 0.083 -0.059In the repository benchmark the two track each other closely on well-specified ARMA — the 1-step prediction error differs by at most about 1% — which is the expected reading: when the model family is right, order selection is not where forecasts are won. The differences that matter are elsewhere:
-
contamination: with a few percent of additive
outliers the classical PACF/ACF paths collapse while
pacf_robustandacf_robustkeep selecting (seevignette("tsorder")); -
observation noise: measurement noise turns an AR(p)
into an ARMA(p, p) and misleads every AR-order selector, but leaves the
McMillan degree —
order_select(x, "svc")— untouched; -
cost: the Hannan-Rissanen proposal is a few
least-squares passes, so it can seed or replace a stepwise search where
thousands of series are fit (the
foresightRust crate uses it exactly that way).
Session info
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.5 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
#> [4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
#> [7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
#> [10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] tsorder_0.1.0
#>
#> loaded via a namespace (and not attached):
#> [1] gtable_0.3.6 jsonlite_2.0.0 dplyr_1.2.1 compiler_4.6.1
#> [5] tidyselect_1.2.1 Rcpp_1.1.2 parallel_4.6.1 jquerylib_0.1.4
#> [9] systemfonts_1.3.2 scales_1.4.0 textshaping_1.0.5 yaml_2.3.12
#> [13] fastmap_1.2.0 lattice_0.22-9 ggplot2_4.0.3 R6_2.6.1
#> [17] generics_0.1.4 knitr_1.52 tibble_3.3.1 desc_1.4.3
#> [21] forecast_9.0.2 timeDate_4052.112 bslib_0.12.0 pillar_1.11.1
#> [25] RColorBrewer_1.1-3 rlang_1.3.0 urca_1.3-4 cachem_1.1.0
#> [29] xfun_0.62 fs_2.1.0 sass_0.4.10 S7_0.2.2
#> [33] otel_0.2.0 cli_3.6.6 pkgdown_2.2.1 magrittr_2.0.5
#> [37] digest_0.6.39 grid_4.6.1 nlme_3.1-169 lifecycle_1.0.5
#> [41] fracdiff_1.5-4 vctrs_0.7.3 evaluate_1.0.5 glue_1.8.1
#> [45] farver_2.1.2 zoo_1.9-1 ragg_1.5.2 colorspace_2.1-4
#> [49] rmarkdown_2.32 pkgconfig_2.0.3 tools_4.6.1 htmltools_0.5.9