# Bayesian order identification of ARMA models with projection predictive
inference

Yann McLatchie  Department of Computer Science, Aalto University, Espoo,
Finland  Asael Alonzo Matamoros  Department of Computer Science, Aalto
University, Espoo, Finland  David Kohns  Department of Computer Science, Aalto
University, Espoo, Finland  Aki Vehtari  Department of Computer Science, Aalto
University, Espoo, Finland

###### Abstract

Auto-regressive moving-average (ARMA) models are ubiquitous forecasting tools.
Parsimony in such models is highly valued for their interpretability and
computational tractability, and as such the identification of model orders
remains a fundamental task. We propose a novel method of ARMA order
identification through projection predictive inference, which benefits from
improved stability through the use of a reference model. The procedure
consists of two steps: in the first, the practitioner incorporates their
understanding of underlying data-generating process into a reference model,
which we latterly project onto possibly parsimonious submodels. These
submodels are optimally inferred to best replicate the predictive performance
of the reference model. We further propose a search heuristic amenable to the
ARMA framework. We show that the submodels selected by our procedure exhibit
predictive performance at least as good as those chosen by AIC over simulated
and real-data experiments, and in some cases out-perform the latter. Finally
we show that our procedure is robust to noise, and scales well to larger data.

###### keywords:

ARMA order identification; Bayesian model comparison; projection predictive
inference.

##  1 Introduction

Since their introduction by Box and Jenkins, (1970), auto-regressive moving-
average models (ARMA) have become ubiquitous forecasting tools. This is thanks
to their high predictive power, ease of implementation, and intuitive
interpretation. Their use is often argued from a theoretical perspective in
time series analysis in that any stationary time series process can be
represented by an infinite order moving-average (MA) model, according to the
Wold decomposition theorem (Wold,, 1938). Such MA models can then be
arbitrarily well approximated by appropriate finite order ARMA models. From an
empirical perspective, they play a pivotal role in the social sciences (among
many others), and particularly economics. Indeed, macroeconomic time series
are often sums of underlying sub-time series (e.g. disaggregated inflation
items feeding into headline inflation) and thus ARMA models naturally arise
due to the sum of auto-regressive time series theorem (Granger and Morris,,
1976). Further, many structural economic models have a moving averaging
representation (Giacomini,, 2013), and ARMA models have been shown to have
competitive forecast performance for aggregate economic time series (Stock and
Watson,, 2007; Koop,, 2013; Chan,, 2013; Zhang et al.,, 2020).

While we don’t claim that ARMA models are better than any other in terms of
their predictive performance, they can provide the statistician with
information on correlation structures and baseline predictive performance.
When building more complex models, such baseline models can operate as sanity
checks and are commonly employed in statistical analyses due to their
effectiveness given their relative simplicity.

Indeed, parsimony in these models is highly valued, not least due to the fact
that latent MA components often result in difficult likelihoods with
increasing order (Chib and Greenberg,, 1994; Ives et al.,, 2010; Chan,, 2013).
Thus one primary goal when implementing ARMA models is to identify small
orders capable of good predictive performance. One popular approach consists
of fitting many such models with different orders and selecting the best one
according to some criterion, most famously the Akaike information criterion
(AIC; Hyndman and Khandakar,, 2008). This approach, however, has been known to
not always select parsimonious models, and may select models with unexpectedly
ragged temporal structures.

In this paper we propose a new method of selecting ARMA and seasonal ARMA
orders motivated from fully Bayesian decision theory. The proposed methodology
contributes to projection predictive inference, originally defined by Goutis
and Robert, (1998) and later developed by Dupuis and Robert, (2003), Piironen
and Vehtari, 2016a , and Piironen et al., (2018) by making it amenable to
detecting relevant ARMA and seasonal ARMA orders from a predictive
perspective. The method follows a two step procedure in which the modeller
first specifies a possibly large model, which incorporates all relevant
knowledge about the underlying data generating process of the data, and passes
posterior predictive checks (Pavone et al.,, 2022; Gelman et al.,, 2020; Gabry
et al.,, 2019).111With provision of a reference model, the model space is
formally “completed”. Crucially, we make no assumption on the existence of a
true model, solely that the reference model passes posterior predictive
checks. This stands in contrast to Bayesian model averaging for which one
assumes an open model space (Vehtari and Ojanen,, 2012). Next, the posterior
predictive information is projected optimally onto possibly parsimonious
submodels. Previous implementations of projective predictive inference have
not been considered for ARMA models and time-series models more generally,
where the projection step is challenging due to the latency of the MA
component and submodel search is complicated by a preference for exploring
increasingly complex ARMA orders, however not at random. This is motivated by
typical stationary time-series having rather smooth than ragged serial
correlation structures, especially when any seasonal time-series correlations
are appropriately modelled. To avoid selection of such ragged ARMA orders, we
propose a forward search heuristic which iteratively increases the order until
reaching the predictive performance of the reference model. This provides a
safeguard against over-fitting, as the projected models typically won’t
exhibit better fit than that of the reference model.

The projection predictive paradigm differs from conventional selection
approaches based on information criteria (Watanabe,, 2013; Spiegelhalter et
al.,, 2002), cross validation (Geisser and Eddy,, 1979), or subset selection
from sparsity inducing priors (Barbieri and Berger,, 2004; Piironen and
Vehtari,, 2017) in that selection is done with respect to point predictions
from a reference model as opposed to the target data directly. Piironen and
Vehtari, 2016a  and Piironen et al., (2018) show that this improves stability
of selection and is coherent with Bayesian decision theory rather than reliant
on asymptotic approximations.222The WAIC is asymptotically equal to Bayesian
LOO-CV, although both induce a small bias in their utility estimates. This can
lead to high variance in the estimates of predictive utility, sub-optimal
model selection, and over-fitting when the model space is large Piironen and
Vehtari, 2016a . Although Piironen and Vehtari, 2016a  show that integrating
over model uncertainty may provide the best predicting model, one may still
use such a model in the proposed methodology as the reference for projection
on more parsimonious submodels.

We aim to convince the reader of the benefits of using reference models for
ARMA order selection by comparing the stability and predictive performance of
reference model-based selection to that achieved by auto.arima (Hyndman and
Khandakar,, 2008). We find through various simulation and real word data
exercises that our procedure identifies models with predictive performance no
worse than those found by auto.arima, and in some cases out-performs the
latter. We also find that our procedure is robust to instances of noisy data,
and avoids over-fitting in those situations where auto.arima liable to do so.
By motivating a robust search heuristic, we show how our procedure is able to
always produce well-performing models in cases where auto.arima fails.

This paper is organised as follows: we begin by discussing the theory behind
ARMA models in Section 2 before defining projection predictive model selection
in Section 3. Having done so, we will be equipped to define our novel order
identification procedure in Section 4 before discussing its position in the
relevant literature in Section 5. We then justify the utility of our procedure
over the most widely used alternative through experiments in Section 6 before
summarising our contributions in Section 7.

##  2 Auto-regressive moving-average processes

In this section, we briefly present the Bayesian ARMA model along with its
decomposition, a central aspect of our proposed order selection approach later
presented in section 4.

###  2.1 The ARMA model

Let the variable of interest $y_{t}$ be observed at integer time points
$t\in\\{1,\dotsc,T\\}$ and assume its time series dynamics are described by an
ARMA model of order $p,q$, denoted as $\operatorname{ARMA}(p,q)$:

$\phi(L)y_{t}=c+\theta(L)\varepsilon_{t},\;\varepsilon_{t}\sim\operatorname{N}(0,\sigma^{2}),$
(2.1)

where $\phi(L)=1-\phi_{1}L-\dotsc-\phi_{p}L^{p}$ and
$\theta(L)=1+\theta_{1}L+\dotsc+\theta_{q}L^{q}$ are lag polynomials for the
AR and MA component respectively, and $L$ is defined as the lag operator
($Lx_{t}=x_{t-1}$). We assume for simplicity throughout the paper that $c=0$
and that the initial conditions
$\\{y_{0},\dotsc,y_{-p},\varepsilon_{0},\dotsc,\varepsilon_{-q}\\}$ are zero
and given.333Note that one can treat these alternatively as unknowns to be
estimated from the data if needed. However, when $T\gg\max(p,q)$, they have
typically little impact on inference on the ARMA dynamics, even when ignored
in the model. Stationarity of Equation 2.1 is not required to achieve proper
posteriors via Bayesian updating with appropriate priors (Sims,, 1988; Sims
and Uhlig,, 1991; Schotman and Van Dijk,, 1991).444These authors in fact show
that the Bayesian approach is valid even when the AR and MA polynomials
display unit roots, and are less prone to incur bias asymptotically for near
unit-root data-generating processes. In order to be able to identify the
optimal ARMA orders for our projection algorithm, we will assume that the
roots of $\phi(L)$ and $\theta(L)$ lie outside the unit circle and thus we
have weak stationarity (Chib and Greenberg,, 1994).555Note that the methods in
this paper are amendable to non-stationary components, $\mu_{t}$ in the data-
generating process, if one proceeds with $y^{*}_{t}=y_{t}-\mu_{t}$.

To define the likelihood, we follow Chan, (2013) and Zhang et al., (2020) by
first stacking all observations according to

$H_{\phi}\bm{y}=H_{\theta}\bm{\varepsilon},\;\bm{\varepsilon}\sim\operatorname{N}(\bm{0},\Omega),$
(2.2)

where $\bm{\varepsilon}=(\varepsilon_{1},\dotsc,\varepsilon_{T})^{T}$,
$\bm{y}=(y_{1},\dotsc,y_{T})^{T}$, $\Omega=\sigma^{2}I_{T}$, and $H_{\phi}$
and $H_{\theta}$ appropriately defined as $T\times T$ difference matrices
(Chan,, 2013). Since these are by definition lower triangular and invertible
for any $\phi=(\phi_{1},\dotsc,\phi_{p})$, we can write
$\bm{y}=H^{-1}_{\phi}H_{\theta}\bm{\varepsilon}$. Then the log-likelihood
function is compactly written as

$\log(\bm{y}|\phi,\theta,\Omega)=-\frac{T}{2}\log(2\pi)-\frac{1}{2}|\Omega_{y}|-\frac{1}{2}(\bm{y}-\bm{c})^{T}\Omega_{y}(\bm{y}-\bm{c}),$
(2.3)

where
$\Omega_{y}=H^{-1}_{\phi}H_{\theta}\Omega(H^{-1}_{\phi}H_{\theta})^{T},\,\bm{c}\in\mathbb{R}^{T}$.

Importantly, the observational model of $\bm{y}$ belongs to the exponential
family of distributions, namely the Gaussian, such that much of the theory for
projection predictive inference in Section 3 follows through immediately. The
assumption of Gaussianity is made for simplicity and analytic tractability,
although our proposed procedure makes no assumption on the observation family
and thus could be naturally extended to other observational models, which we
leave for future research.

To complete the model specification, we follow the recommendation of Matamoros
and Torres, (2021) who assume independent priors: rl ϕ &∼N(0,Λ_Φ)
θ ∼N(0,Λ_Θ)
c ∼Student-t(0,σ_c,ν)
σ ∼Student-t_+(0,σ_σ,ν_σ), where
$\operatorname{Student-}t_{+}(\cdot,\cdot,\cdot)$ refers to the half-
Student-$t$ distribution with positive support, and
$\Lambda_{\Phi},\,\Lambda_{\Theta}$ denote the diagonal prior covariances
motivated by Matamoros and Torres, (2021). While a plethora of variable
selection and shrinkage priors have been proposed for highly parameterised
time series forecasting models such as vector auto-regressive (VAR) models
(Bańbura et al.,, 2010; Koop,, 2013; Carriero et al.,, 2015; Giannone et al.,,
2015; Chan et al.,, 2016), this is less of a concern for relatively more
parsimonious ARMA models. In this paper, we purposefully consider a
statistician who uses default priors to compare our procedure with alternative
(potentially non-Bayesian) ARMA selection techniques such as the popular
auto.arima.

Independently to the choice of priors, it is well known that ARMA likelihoods
can be multi-modal (Chan and Chen,, 2011) and additionally pose computational
challenges due to the latent MA components.666Traditionally used Kalman
filters for maximum likelihood formulations of ARMA models (Harvey,, 1985) can
become computationally demanding, particularly with large $T$ (Kim et al.,,
1999). To aid posterior inference, all models in this paper are estimated via
Hamiltonian Monte Carlo (Neal et al.,, 2011) using the so-called no U-turn
sampling (NUTS) algorithm as implemented in Stan (version 2.26.1; Hoffman and
Gelman,, 2011; Carpenter et al.,, 2017).

###  2.2 Decomposing the ARMA model

ARMA models pose the computational problem that the MA component
$\theta(L)\varepsilon_{t}$ are unobserved. To deal with this, we follow the
logic presented in least-square literature on ARMA estimation (Kapetanios,,
2003) by splitting the $\operatorname{ARMA}(p,q)$ model into its auto-
regressive and moving average components and then estimating these in two
separate steps as proposed by Hannan and Rissanen, (1982).

To illustrate this, assume oracle knowledge of the order from the respective
AR and MA components, and begin by fitting $\phi(L)y_{t}=\delta_{t}$, where
$\delta_{t}$ is some white noise process, without observing $\varepsilon_{t}$.
Then, we fit a linear model to the residuals, $\hat{\delta}_{t}$, to
approximate the MA component. Formally, rl ϕ(L)y_t &= δ_t
θ(L)^δ_t = ξ_t, noting that $\hat{\delta}_{t}\approx\varepsilon_{t}$ from
Equation 2.1, and $\xi_{t}$ are from a white noise. Ng and Perron, (1995)
detail the sufficient conditions required to estimate $\hat{\delta}_{t}$
consistently, which in turn allows for consistent estimation of the ARMA order
from a frequentist perspective. In performing this sequential model fitting,
we also simplify our search heuristic when traversing the model space as we
shall see later in Section 3.2, and forms the basis of our two-step procedure
later defined in Section 4.

###  2.3 Seasonal ARMA models

One innovation on the base ARMA model previously discussed is to model
recurring seasonal trends. We can do so by defining structural relationships
occurring at regular lag intervals, for example each week, quarter, year, and
so on. The multiplicative seasonal auto-regressive moving-average (SARMA)
model with seasonal patter repeating every $s$ lags, non-seasonal ARMA
components $p,q$ and seasonal components $P,Q$ is denoted
$\operatorname{ARMA}(p,q)\times(P,Q)_{s}$ and written in terms of polynomials
in $L$, as

$\Phi(L^{s})\phi(L)y_{t}=\Theta(L^{s})\theta(L)\varepsilon_{t}.$

These SARMA models extend the theory of ARMA models to fit and forecast more
complex, macro trends in time series data and thus lend themselves to more
valuable application.

##  3 Projection predictive model selection

In this section we outline the underlying theory of projection predictive
inference, as well as our contributions to this approach in making it amenable
to ARMA models.

###  3.1 The projection

The idea of projection predictive inference is to separate prediction and
model selection into two stages: we first identify a model that produces the
best predictive performance given the information set of the statistician,
which we call the reference model (Vehtari and Ojanen,, 2012); we then
construct smaller models capable of replicating the predictive performance of
this reference model. These smaller models are essentially fit to the fit of
the reference model (a procedure we call the projection), and are then used
for their improved interpretability, or to decrease data collection cost
(Piironen et al.,, 2018). We usually identify these smaller models via
sparsity-favouring search heuristics (Hahn and Carvalho,, 2015; Ray and
Bhattacharya,, 2018; Kohns and Potjagailo,, 2022). Using a reference model in
the model selection process has been shown previously by Piironen and Vehtari,
2016a , and more recently by Pavone et al., (2022), to lead to more stable
submodel discovery in which submodels are less prone to over-fit the data than
those found through other procedures that conduct selection directly on
$y_{t}$. A salient reason for the stability of the projection is the fact that
a well-defined reference model is able to separate signals of the data-
generating process from noise (Piironen et al.,, 2018).

We start, then, with the model containing all available data which we fit
using reasonable priors, and which has passed posterior checks (Gelman et
al.,, 1996, 2020; Gabry et al.,, 2019). This model is written in terms of the
full parameter space $\bm{\theta}^{\ast}\in\bm{\Theta}^{\ast}$, and is fit to
the observed data $\mathcal{D}=\\{X,y\\}$. From this reference model, we wish
to achieve some more parsimonious model in the restricted parameter space
$\bm{\theta}^{\perp}\in\bm{\Theta}^{\perp}$ that will often incorporate
sparsity.777In general, we do not require this restricted parameter space to
be a subset of the reference parameter space, although it is usually chosen to
induce sparsity in the submodels. Concretely, we replace the posterior
distribution over the reference model parameters
$p(\bm{\theta}^{\ast}\mid\mathcal{D})$ with some simpler distribution
$q(\bm{\theta}^{\perp})$ such that its induced covariate-conditional
predictive distribution, $q(\tilde{y}\mid\bm{\theta}^{\perp})$, is not
significantly different to that of the reference model. We quantify this with
some distance measure,

$d\left\\{p(\tilde{y}\mid\mathcal{D},\bm{\theta}^{\ast}),\,q(\tilde{y}\mid\bm{\theta}^{\perp})\right\\}\leq\varepsilon,$
(3.1)

where $\tilde{y}$ represents predictions of the variate, $d$ is a distance
measure, and $\varepsilon>0$ is small. Vehtari and Ojanen, (2012) consider
projected posteriors from a decision-theoretic standpoint and reason that
given a logarithmic utility function, the optimal values of the restricted
model parameters $\bm{\theta}^{\perp}$ are achieved by minimising the
Kullback-Leibler (KL) distance between the posterior predictive distributions
of the reference and restricted model. This divergence benefits from its
analytical tractability and efficient computation within the exponential
family of distributions, as well as the guarantee of a unique optimum and
natural position within a Bayesian workflow (Goutis and Robert,, 1998;
Piironen et al.,, 2018).

Assume we have collected $S$ posterior draws from the reference model,
$\\{\bm{\theta}_{(s)}^{*}\\}_{s=1}^{S}$. The projection is then simply the
solution to the optimisation problem,

$\bm{\theta}^{\perp}_{(s)}=\operatorname*{arg\,min}_{\bm{\theta}\in\bm{\Theta}^{\perp}}\mathbb{KL}\left\\{\,p(\tilde{y}\mid\mathcal{D},\bm{\theta}^{\ast}_{(s)})\,||\,q(\tilde{y}\mid\bm{\theta}_{(s)})\,\right\\},$

for which one can show analytical solutions exit when the models are contained
within the exponential family of distributions. Piironen et al., (2018)
propose to solve this problem sample-wise for computational reasons, leaving
us with some set of projected samples we can consider as samples from the
submodel’s posterior (and which we achieved at a negligible cost since no
Monte Carlo methods were used in projection step). We may then evaluate the
new submodel’s predictive performance with such metrics as expected log-
predictive density (elpd; Vehtari et al.,, 2016) using leave-one-out cross-
validation (LOO-CV) to compare it to the reference model.

Since LOO-CV is computationally demanding, particularly when $T$ is large, we
will use the recently proposed approximate LOO-CV using Pareto smoothed
importance sampling (PSIS; Vehtari et al.,, 2016). This represents a fully
probabilistic way of doing LOO-CV that avoids repeatedly fitting the reference
model by re-weighting the posterior draws
$\\{\bm{\theta}_{(s)}^{*}\\}_{s=1}^{S}$ with importance weights. In
particular, the weight for draw $\bm{\theta}^{*}_{(s)}$, leaving the $t^{th}$
out, denoted $\omega_{(s)}^{-t}$ where $t$ indexes the left-out observation,
is given by

$\omega_{(s)}^{-t}\propto\frac{1}{p(y_{t}|\bm{\theta}^{*}_{(s)})}.$

These weights are then stabilised with Pareto smoothing for instances in which
the importance weight distribution has a thick tail (Vehtari et al.,, 2016).

###  3.2 Search strategies for traversing the ARMA model space

$\small\\{\text{Intercept}\\}$$\\{\text{Intercept},\phi_{1}\\}$$\cdots$$\cdots$$\\{\text{Intercept},\phi_{1},\phi_{4}\\}$$\cdots$$\cdots$$\cdots$$\cdots$$\cdots$

(a) A forward search heuristic in the auto-regressive space, in which each
step considers all previously unconsidered parameters to be appended to the
path regardless of their lag or the feasibility of the resulting model.

$\\{\text{Intercept}\\}$$\\{\text{Intercept},\phi_{1}\\}$$\cdots$$\\{\text{Intercept},\phi_{1},\phi_{2},\phi_{3},\phi_{4}\\}$$\cdots$

(b) The proposed search path through the parameter space, where each next
parameter is that one situated just one lag further from the observed point,
thus considering only valid ARMA models.

Figure 1: A comparison of parameter search heuristics discussed in Section
3.2. Each node in the tree represents a submodel’s parameters, and the path
from the intercept-only model to the final model (from left to right in the
graphs) represents the search path.

In order to determine which parameter subsets to project our reference model
onto, we need some search heuristic to propose a collection of submodels.
Piironen et al., (2018) propose a forward search heuristic in the case of
generalised linear models, wherein candidate submodels are found by
iteratively appending the variable minimising the distance in Equation 3.1
starting with the intercept-only model. Such search heuristics do not adapt
well to the ARMA model.

For instance, suppose we have an $\operatorname{AR}(1)$ model and to it we add
the $\phi_{4}$ parameter in our search. In doing so, one may be prone to
skipping intermediate lags which creates a seasonal behaviour in the AR
dynamics. This, however, should not appear, particularly when modelling
seasonal components as we do in section 4.2 and making sure that outliers and
departures from the stationary component of the data-generating process are
adequately modelled in the reference model. This example is seen in Figure
1(a).

Different from order identification schemes such as Hyndman and Khandakar,
(2008), we implement a search heuristic that iteratively appends the next
lagged variable to the set, since this is directly equivalent to increasing
the order of our AR or MA model by one with each added parameter. For example
in Figure 1(b), the submodel containing $\phi_{4}$ is indeed an
$\operatorname{AR}(4)$ model.

Alternative order selection procedured proposed by Nardi and Rinaldo, (2011)
and Chan and Chen, (2011) instead conduct variable selection directly on the
ARMA polynomials via frequentist lasso style regularisation. They prove that
under some regularity conditions, an adaptive Lasso regression of the time
series on its lags enjoys oracle properties asymptotically. We do not consider
this any further since we remain primarily interested in the finite data
regime.

###  3.3 Submodel acceptance heuristics

Having projected our reference model onto different parameter subsets, we then
move on to identify the smallest submodel whose predictive performance is
comparable to that of the reference model. We propose the submodel selection
heuristic outlined by Piironen and Vehtari, 2016a  based on the sensible
cross-validation of submodels following a forward search through the parameter
space. Denote the elpd of the reference model as
$\operatorname{elpd}(\mathcal{M}^{\ast}\mid y)$ and that of the submodel on
$p$ parameters as $\operatorname{elpd}(\mathcal{M}_{p}\mid y)$. Having
projected our reference model onto a set of submodels, we choose the submodel
with the smallest $p$ for which the upper bound of its normal-approximation
$68\%$ elpd confidence interval (the one standard deviation range) is at least
as good as the elpd point estimate of the reference model:

$\mathbb{P}\left(\operatorname{elpd}(\mathcal{M}^{\ast}\mid
y)-\operatorname{elpd}(\mathcal{M}_{p}\mid y)\leq 0\right)\geq 0.68.$ (3.2)

##  4 Projection predictive ARMA order identification

Having now addressed both projection predictive model selection and the ARMA
model, we combine the former with the decomposition of ARMA models outlined in
Section 2.2, the selection heuristic proposed in Equation 3.2 and the search
heuristic from Section 3.2 to build our order identification procedure.

###  4.1 A fully Bayesian ARMA order identification procedure

Suppose we have some data we wish to model with an ARMA model

$\phi(L)y_{t}=c+\theta(L)\varepsilon_{t},$

where $\\{\varepsilon_{t}\\}$ are some noise in our data and assume that we
have fit an ARMA reference model that passes predictive checks.888More
generally, one might consult the empirical auto-correlation functions to
determine the maximum $\operatorname{ARMA}(p,q)$ to start the analysis..
Denote the orders of the reference model as $(p^{*},q^{*})$. Taking
$(\\{y_{t}\\},p^{\ast},q^{\ast})$ as inputs, we propose the application of
projection predictive inference as outlined in Algorithm 1, and which we
briefly describe below.999An implementation of Algorithm 1 in R is available
at https://github.com/yannmclatchie/projpred-arma.

Algorithm 1 Projection predictive order identification of ARMA models

1:procedure ProjpredARMA($\\{y_{t}\\},p^{\ast},q^{\ast}$)

2: Fit an auto-regressive model with HMC to data observations $\\{y_{t}\\}$
with the lag $p^{\ast}$ from the reference model,
$\operatorname{AR}(p^{\ast})$, using MCMC and reasonable priors.

3: Treat this $\operatorname{AR}(p^{\ast})$ model as the auto-regressive
reference model, and apply projection predictive model selection to it using
the temporal search method defined in Section 3.2 and extract the restricted
lag value $p^{\perp}$ with the selection heuristic defined in Equation 3.2.

4: Fit a linear model to the residuals of the $\operatorname{AR}(p^{\perp})$
model $\\{\varepsilon_{t}^{\perp}\\}$ (through posterior predictive mean point
estimates) with HMC, and using the reference model MA lag, $q^{\ast}$.

5: Apply the same modified projection predictive model selection as in step 3
to this linear reference model of order $q^{\ast}$ from step 4 and retrieve
$q^{\perp}$.

6: return $p^{\perp},q^{\perp}$.

7:end procedure

First, refit an $\operatorname{AR}(p^{*})$ model to $y_{t}$. For instance,
should we find that a reference $\operatorname{ARMA}(5,5)$ fits the data well,
then we would fit a AR reference model of structure $\operatorname{AR}(5)$.

Second, we perform projection predictive model inference on this
$\operatorname{AR}(5)$ model to find a possibly more parsimonious restricted
model of order $\operatorname{AR}(p^{\perp})$, say an $\operatorname{AR}(1)$
for illustration, using the temporal search heuristic in 3.2. We save the
residuals of the projected AR model, denoted as
$\\{\varepsilon_{t}^{\perp}\\}$ since these will serve as the observed errors
to project the MA component next. Note that this step produces an implied
distribution over the residuals upon prediction. For simplicity, we define
$\\{\varepsilon_{t}^{\perp}\\}$ here as the posterior mean of the residuals.

Third, we fit an linear model of the order defined by the reference model (in
this case $q^{\ast}=5$) to the residuals $\\{\varepsilon_{t}^{\perp}\\}$.

Finally, we perform projection predictive model inference on this linear model
fit to the residuals to retrieve a restricted MA order $q^{\perp}$ (reducing
an $\operatorname{MA}(5)$ to, say, an $\operatorname{MA}(3)$). The combined
orders of the restricted submodels identify the order of the projected
$\operatorname{ARMA}(p^{\perp},q^{\perp})$ that most closely replicates the
predictive performance of the much larger $\operatorname{ARMA}(p,q)$ reference
model as measured by Kullback-Leibler distance.

This methodology allows us to search the parameter space in an intuitive
manner in terms of two linear models sequentially and with minimal information
loss while benefiting fully from the stability afforded by a reference model.

While these $p^{\ast}$ and $q^{\ast}$ can in theory take any value, we will
limit them in our experiments to $p^{\ast}=q^{\ast}=5$ in line with the
default values proposed by Hyndman and Khandakar, (2008). Indeed, if we were
to set these reference lags to larger values, we would have to choose our
priors appropriately to communicate the fact that more distant lags are less
likely to have an effect on the present and in order to enforce model
stationarity.

###  4.2 Extension to SARMA models

We presently show how we might naturally extend the procedure presented in
Algorithm 1 to SARMA models. In this scenario, we define our reference model
as a function of the seasonality being modelled. Specifically, given a
reference model
$\operatorname{ARMA}(p^{\ast},q^{\ast})\times(P^{\ast},Q^{\ast})_{s}$, with
both seasonal and non-seasonal components and seasonality $s$, we produce two
datasets: one seasonal and one non-seasonal to which we apply our procedure
independently. The outputs of these two runs provide us the non-seasonal and
seasonal restricted parameters respectively. Again, in line with Hyndman and
Khandakar, (2008), we choose as default values $P^{\ast}=Q^{\ast}=3$ in our
experiments.

###  4.3 Cross-validation with time series data

In previously seen experiments with projection predictive model selection,
such as those carried out by Piironen and Vehtari, 2016a , Piironen et al.,
(2018), Catalina et al., (2021), and Catalina et al., (2022), evaluation of
each of the submodels was performed with LOO-CV for efficiency. Approximate
leave-future-out cross-validation (LFO-CV; Bürkner et al.,, 2020) provides an
alternative to this for time series data cross-validation.

In order to convince the reader that the computationally cheaper LOO-CV is
justified in our case, consider a time series model $p(y_{i}\mid
f_{i},\bm{\theta})$ with latent process values $\bm{f}$ and prior $p(\bm{f})$
such that $(y_{i}-f_{i})$ are independent and identically distributed. Now, if
we are interested in predicting unseen future observations, then LFO-CV is the
natural choice. However, should we instead be interested in reasoning on the
structure of our time series, that is the conditional observation model
$p(y_{i}\mid f_{i})$, then Bürkner et al., (2020) show that it is reasonable
to use LOO-CV instead. The authors reason that LOO-CV can be thought of as a
biased approximation to LFO-CV, and further that the biases in LOO-CV and LFO-
CV are likely in the same direction (Bürkner et al.,, 2020). This means that
since the difference between model comparison through LOO-CV and LFO-CV is
small, and since it has been empirically noted that covariate ordering is
similar between the two, we are able to rely on the computationally cheaper
LOO-CV for variable selection in time series models.

The use of the AIC and the Bayesian information criterion (BIC; Schwarz,,
1978), both commonly implemented in analyses, are themselves asymptotically
equivalent to LOO-CV and $K$-fold-CV respectively (Stone,, 1977; Shao,, 1997;
Arlot and Celisse,, 2010; Vehtari and Ojanen,, 2012). As such, using LOO-CV
can be considered the finite-sample analogue of the asymptotically equivalent
AIC.

###  4.4 Alternative projections

In general, projection predictive inference accommodates the possibility of
projecting an arbitrary reference model structure onto a set of models with
different structures. We briefly discuss how our procedure in Algorithm 1
could be adapted to achieve differently nuanced model selection results.

We might be tempted, for example, to project our reference ARMA model directly
onto the AR component rather than perform the initial decomposition, and then
continue the procedure from step 4 of Algorithm 1. We would expect this to
then over-select the AR size and under-select the restricted MA component,
since we then project the information communicated by both the AR and MA onto
solely the AR component of the restricted model. There do exist, however,
instances where such projections may afford the statistician a reduced
computational cost. Chan, (2013) discusses the difficulty of inferring MA
parameters often present in economic models. As a remedy to this, and
understanding that the inclusion of an MA component in the ARMA model induces
an infinite AR order, one might project their ARMA model onto only an AR
component purposefully to achieve the minimal order necessary to replicate the
behaviour of an $\operatorname{AR}(\infty)$ model at hopefully less
computational cost.

Further, if instead of finding parsimonious submodels, our interest is in
identifying the true order of the ARMA with our approach, one might consider
application of so-called complete variable selection as presented by Pavone et
al., (2022), where the statistician is interested in identifying all
covariates (theoretically) relevant to predicting the target. Here, the
projective inference step is repeatedly applied until some stoppage criterion
is met.101010Pavone et al., (2022) suggest using local false discovery rates
(Efron,, 2008, 2010), empirical Bayes median (Johnstone and Silverman,, 2004),
and posterior predictive credible intervals. Such an approach can similarly be
applied to update the projected residuals with each projection step.

We show how these slightly different projection techniques can be implemented
in Section 6.3. We remain primarily interested in identifying parsimonious
ARMA orders in this paper, and so do not consider complete variable selection
any further.

##  5 Related work

Other model selection procedures exist in the case of ARMA models. We brielfy
discuss the most competitive alternative to our procedure in auto.arima, and
motivate why other procedures are inadequate to identify and communicate ARMA
subset selection.

###  5.1 Automatic order selection with unit root tests and AIC

The forecast package in R, developed by Hyndman and Khandakar, (2008), has
long been the modus operandi of statistical practitioners for the order
identification of ARMA models. In particular, the auto.arima module automates
an order selection heuristic based on unit root tests and the AIC.

In the original algorithm proposed by Hyndman and Khandakar, (2008), some unit
tests are performed on the data and based on their results, they fit four
models by maximum likelihood. These four initial models are pre-defined before
we see any data and without any prior information. Of these four models, they
select the one with the lowest (best) AIC, and consider a further thirteen
variations of it. The model with the lowest AIC of these resulting models is
chosen as the best submodel. Various constraints are imposed to ensure
convergence or near unit root, and these constraints combined with the finite
search space guarantees that at least one of the models considered will be
valid (Hyndman and Khandakar,, 2008).

Nevertheless, due to its competitive performance and popularity, we will
consider auto.arima for comparison in later experiments. As such, we use the
default algorithm in Hyndman and Khandakar, (2008) to identify the ARMA orders
and then conduct inference via Hamiltonian Monte Carlo, using the same priors
as in Algorithm 1 so that we can compare predictive performance with our
procedure. This is summarised in Algorithm 2.

Contrasting Algorithms 1 and 2, we highlight two main reasons why projective
inference may result in most stable submodel discovery. Firstly, it is well
known that priors can have a regularising impact on the posterior, and that
MCMC methods such as HMC allow us to explore the posterior more efficiently
than maximising the likelihood directly. Secondly, projective inference
conducts submodel selection based on predictions of a reference model, whereas
auto.arima conducts its selection directly on the observations $\\{y_{t}\\}$.
This use of a reference model helps submodel discovery by typically filtering
out noise from the observations, and can also help avoid over-fitting to the
data since the fit of the submodels is bounded by that of the reference
model.111111In fact McQuarrie and Tsai, (1998) note the tendency of AIC-based
model selection to over-fit data in instances of small sample size. This will
be formally investigated in Section 6.1.

Algorithm 2 MCMC fitting of an auto.arima-selected model.

1:procedure MCMC AutoARIMA($\\{y_{t}\\}$)

2: Identify some orders $p_{\operatorname{auto}},q_{\operatorname{auto}}$ from
differenced, stationary data with auto.arima.

3: Fit the
$\operatorname{ARMA}(p_{\operatorname{auto}},q_{\operatorname{auto}})$ model
with previously found orders using MCMC, employing some reasonable priors.

4:end procedure

###  5.2 Cross-validation

Another common approach to model selection is to fit some collection of models
and estimate their respective elpd scores (or any other scoring rule, e.g.
those proposed by Gneiting and Raftery,, 2007) with cross-validation,
whereupon the optimal subset according to the highest score is chosen.

While this procedure has gained popularity (Arlot and Celisse,, 2010), and
while it has been shown to be a robust method when dealing with relatively few
models, if the number of models being compared is relatively large or the
number of covariates is relatively small, various issues arise. Piironen et
al., (2018) showed that when the number of models is large, then cross-
validation without a reference model is liable to over-fit and result in the
selection of a sub-optimal model. Indeed, Piironen and Vehtari, 2016a
compared this approach directly with projection predictive model selection,
finding that the latter is significantly more resilient to these issues.
Pavone et al., (2022) further showed that the use of a reference model in
model selection affords a greatly improved stability in selection which is due
to that fact that the reference model is able to filter out noise before
arriving at the submodel selection stage.

A computational hurdle for Bayesian workflows is additionally that full cross-
validation requires fitting many models for which MCMC is re-conducted for
each left-out observation. This is a highly expensive endeavour. To summarise,
we present in Table 1 algorithmic complexities of the reviewed methods.

Procedure | Worst-case complexity
---|---
ProjpredARMA | $\mathcal{O}(1)$
MCMC AutoARIMA | $\mathcal{O}(1)$
Cross-validation | $\mathcal{O}(p^{\ast}\cdot q^{\ast})$
Table 1: A comparison of the algorithmic complexity of the four proposed model
selection procedures as a function of the number of models needed to be fit by
MCMC.

Since fitting a model with MCMC represents the largest user and computational
cost, we measure the complexity of each procedure as a function of the number
of models needing to be fit by MCMC and present their worst case algorithmic
complexities. In ProjpredARMA the number of models needed to be fit by MCMC is
independent of the size of the reference model, namely, we will always fit two
models (AR and MA components) with MCMC regardless of the maximum lags
$p^{\ast}$ and $q^{\ast}$. Similarly, in our MCMC AutoARIMA procedure, we
leverage the speed of auto.arima to only fit one model by MCMC. However, when
using a pairwise selection criterion such as cross-validated elpd, in the
worst case we need to fit all possible models encompassed by maximum lags
$p^{\ast}$ and $q^{\ast}$. Consequently we find that the number of models
needed to be fit has complexity $\mathcal{O}(p^{\ast}\cdot q^{\ast})$.

###  5.3 Using sparsifiying priors for ARMA models

A commonly used alternative to variable selection for generalised linear
models are sparsifying priors over the regression coefficients such as the
regularised horseshoe prior (Piironen and Vehtari,, 2017), the spike-and-slab
prior (Mitchell and Beauchamp,, 1988) or R2D2 prior (Zhang et al.,, 2022) (see
e.g. Polson and Scott,, 2012; Bhadra et al.,, 2019, for excellent reviews on
further shrinkage priors). Such priors force certain parameter values close to
zero based on their relevance to posterior predictions.

Such priors face issues in our case. Namely, as is noted by Catalina et al.,
(2021), the posterior of a model fitted with a sparsifying prior is not truly
sparse in that parameter posteriors are not generally point masses at zero,
and manual effort is required from the statistician to first conceive a
threshold of posterior relevancy, and then to prune those parameters beneath
it. Then, it is not clear that these sparsified posteriors represent intuitive
or desirable ARMA models in general, as such priors are usually used under the
assumption of exchangeability of the regression weights. When the data are
highly correlated, as we expect with AR and MA lags, aggressive shrinkage
priors may cause the marginal posteriors of the regression weights to overlap
strongly with zero, thus falsely indicating insignificant lags. For instance,
we might identify that the first three lags may be pruned based on the
individual marginal posteriors, yet they are important to include so as not to
create unintended seasonal patterns in the predicted ARMA (much like in Figure
1(a)). And in the particular case of the spike-and-slab prior, the optimality
of selecting the median probability model for prediction assumes orthogonality
of covariates, which is not the case in ARMA models. Indeed Barbieri and
Berger, (2004) put forward a case in which the median probability model with
correlated covariates is clearly sub-optimal. We therefore do to not consider
sparsifying priors any further as an alternative model selection procedure in
our case. For a comprehensive comparison of projection predictive inference
with sparsifying priors and median probability model, and a comparison with
projection predictive inference, see the review by Piironen and Vehtari, 2016a
.

Alternatively to using shrinkage priors, the so-called Minnesota prior
(Giannone et al.,, 2015) and its adaptive variant (Chan,, 2021) can be used to
shrink lag effects as a function of their distance from the current
realisation through time. One might also consider functional restrictions on
AR and MA lag polynomials such as used in mixed-frequency applications (Kohns
and Potjagailo,, 2022; Mogliani and Simoni,, 2021). While useful to some
applications, such restrictions create irreducible bias when the restrictions
are not approximately correct. If, however, sufficient knowledge of such
restrictions exist, they may provide an adequate reference model within the
proposed framework.

##  6 Experiments

We presently demonstrate the practical value of our proposed procedure
compared to auto.arima. We do so first by using it to identify predictive
submodels from multiple simulations of different data-generating processes and
comparing the stability of submodel selection of the two procedures, as well
as their closeness to the true process. We then fit some reference models to a
selection of well-known datasets and compare our procedure to auto.arima in
achieving predictive submodels. Having investigated the procedures’ stability
and the predictive performance of their submodels, we then compare the
performance of the two procedures under different noise regimes. Finally, we
illustrate the behaviour of the search heuristic motivated in Section 3.3 in
an auto-regressive example with many distant, near-zero lags compared to
auto.arima, and conclude with a demonstration of how our procedure scales to
larger data.

The models used in the experiments were fitted with the defaults priors
suggested by Matamoros and Torres, (2021) and all experiments were performed
with a modified version of projpred (Piironen et al.,, 2022).

###  6.1 Stability in model selection

(a) ARMA(0,2) model.

(b) ARMA(1,0) model.

(c) ARMA(1,2) model.

(d) ARMA(2,1) model.

Figure 2: Simulated data. Histograms showing the frequency of model orders
selected by projection predictive inference ($p^{\perp}$ and $q^{\perp}$ in
green) compared to auto.arima ($p_{\operatorname{auto}}$ and
$q_{\operatorname{auto}}$ in orange). The reference model parameters used,
$p^{\ast}$ and $q^{\ast}$, and the true data-generating parameters
$p^{\text{true}}$ and $q^{\text{true}}$ are labelled and shown by the dotted
black lines. We find that across non-seasonal simulated data experiments,
projection predictive inference benefits from far greater stability than
auto.arima. Figure 3: Simulated data. Histograms showing the frequency of
model orders selected by projection predictive inference (in green) compared
to auto.arima (in orange) for data sampled from a
$\text{SARMA}(1,2)\times(1,0)_{12}$ model. The reference model parameters and
the true data-generating parameters are shown by the dotted black lines. We
find that projection predictive inference again benefits from increased
stability compared to auto.arima across seasonal simulated data experiments.

We presently simulate data according to four $\operatorname{ARMA}(p,q)$
processes, each with different orders $p$ and $q$, and aim to recover some
restricted model with predictive performance close to a reference
$\operatorname{ARMA}(5,5)$ model. Namely, the true data generating processes
we sample from are the $\operatorname{ARMA}(1,0),$ $\operatorname{ARMA}(0,2)$,
$\operatorname{ARMA}(2,1)$, and $\operatorname{ARMA}(1,2)$ models. We will
then demonstrate the validity of the extension of the procedure to SARMA
processes by repeating this for an $\operatorname{ARMA}(1,2)\times(1,0)_{12}$
model. For each model, we run $100$ simulations each generating $500$ data
points. For each of these simulated series, we employ both projection
predictive inference and auto.arima to identify parsimonious model orders.

We define stability in model selection concretely in terms of the
concentration of the distribution of model orders selected. In a word, we call
a model selection procedure “stable” if the orders of the ARMA model selected
are similar across repetitions.

In Figure 2 we see that the auto-regressive orders identified by projection
predictive inference are significantly more stable than those selected by
auto.arima. The moving-average component is also broadly more stable, but not
as significantly so. This is likely a symptom of the model decomposition we
propose, where some noise is perhaps leaking into the residuals of the AR
component and contaminating our decision. We further note that these highly-
stable values of $p^{\perp}$ and $q^{\perp}$ are either precisely the true
model values or were close to the true values. When our procedure
“incorrectly” selects lag values (by which we mean not exactly the true
model), most selections still fall within one lag from the actual value,
showing that our procedure has learned some of the underlying structure in the
time series data even in the worst cases. The AR component retrieves exactly
the correctly lag in the absence of an MA component in the ARMA data-
generating process. Indeed, since the inclusion of any MA component leads to
an infinite AR component, it is expected that our restricted AR order will
marginally over-select – as is the case. We include this analysis to
demonstrate that despite the aim of the procedure not explicitly being about
retrieving a true model (since we do not assume one exists), it is able to
identify some underlying structure important to predictive performance.

Figure 3 shows the results of the SARMA experiment. Interestingly, while our
procedure is able to identify seasonal lags that are both highly stable and
close to the true model, auto.arima fails to recognise a seasonal component at
all. This behaviour was seen across other simulated examples not shown here,
and is perhaps an artifact of the AIC-based selection criterion. As well as
this out-performance in the seasonal component, projection predictive
inference is once more closer to the truth and more stable in its non-seasonal
selection when compared to auto.arima.

###  6.2 Predictive performance

Having shown that projection predictive model selection is considerably more
stable than auto.arima, we presently show that this stability does not come at
the expense of submodel predictive performance.

To this end, we use both projection predictive model selection and auto.arima
to identify parsimonious submodels for well-studied datasets curated by
Hyndman and Athanasopoulos, (2021). We then record the orders selected and the
predictive performance of their MCMC fitted models as measured by LOO-CV elpd.
Finally, we measure the difference between the models’ elpds in an effort to
ascertain the magnitude of difference between the two model’s predictive
performances.

We can understand difference in predictive performance between the two
procedures through the elpd difference (diff.) column of tables 2 and 3. When
this difference is positive, the mean elpd of the model selected by our
procedure (Algorithm 1) was higher (better) than that selected by auto.arima
(Algorithm 2), and vice versa. We embolden differences for which zero is not
the approximate 90% normal interval over the mean elpd difference estimate
(the $1.64$ standard deviation range).121212See the work of Sivula et al.,
(2020) for the properties of the sampling distribution of LOO-CV elpd
differences.

| ProjpredARMA | MCMC AutoARIMA |  |  |
---|---|---|---|---|---
Data | $p^{\perp}$ | $q^{\perp}$ |  elpd $\pm$ s.e. | $p_{\operatorname{auto}}$ | $q_{\operatorname{auto}}$ |  elpd $\pm$ s.e. |  elpd diff. $\pm$ s.e.
A | $2$ | $3$ | ${}-100{}$ | ${}\pm{}$ | ${}12{}$ | $0$ | $1$ | ${}-130{}$ | ${}\pm{}$ | ${}13{}$ | ${}\bm{30}{}$ | ${}\bm{\pm}{}$ | ${}\bm{7.2}{}$
B | $1$ | $0$ | ${}7.3{}$ | ${}\pm{}$ | ${}3.3{}$ | $1$ | $0$ | ${}7.4{}$ | ${}\pm{}$ | ${}3.3{}$ | ${}0.0{}$ | ${}\pm{}$ | ${}0.1{}$
C | $2$ | $3$ | ${}-108{}$ | ${}\pm{}$ | ${}7{}$ | $2$ | $1$ | ${}-110{}$ | ${}\pm{}$ | ${}7{}$ | ${}1.4{}$ | ${}\pm{}$ | ${}1.0{}$
D | $1$ | $1$ | ${}-80{}$ | ${}\pm{}$ | ${}5{}$ | $1$ | $1$ | ${}-80{}$ | ${}\pm{}$ | ${}5{}$ | ${}0.0{}$ | ${}\pm{}$ | ${}0.1{}$
E | $3$ | $0$ | ${}-608{}$ | ${}\pm{}$ | ${}25{}$ | $0$ | $2$ | ${}-613{}$ | ${}\pm{}$ | ${}25{}$ | ${}\bm{5.1}{}$ | ${}\bm{\pm}{}$ | ${}\bm{3.1}{}$
F | $2$ | $0$ | ${}-110{}$ | ${}\pm{}$ | ${}7{}$ | $0$ | $1$ | ${}-109{}$ | ${}\pm{}$ | ${}6{}$ | ${}-1.0{}$ | ${}\pm{}$ | ${}1.4{}$
G | $1$ | $3$ | ${}1.9{}$ | ${}\pm{}$ | ${}5.9{}$ | $1$ | $1$ | ${}1.4{}$ | ${}\pm{}$ | ${}6.2{}$ | ${}0.5{}$ | ${}\pm{}$ | ${}0.9{}$

  * •

Dataset lookup: A $=$ airline passengers ($n=47$), B $=$ international
visitors ($n=36$), C $=$ Lake Huron bathymetry ($n=98$), D $=$ insurance
quotes ($n=40$), E $=$ Ansett Airline passengers ($n=269$), F $=$ maximum
annual temperature ($n=46$), G $=$ female murder rate ($n=55$).

Table 2: Non-seasonal real world data. Predictive performance comparison of
ProjpredARMA with MCMC AutoARIMA on non-seasonal data sets curated by Hyndman
and Athanasopoulos, (2021). For both procedures, we report the order of the
suggested model (denoted $p^{\perp}$ and $q^{\perp}$ in the case of the
projected orders, and $p_{\operatorname{auto}}$ and $q_{\operatorname{auto}}$
for auto.arima) as well as the models’ elpd and elpd standard error (s.e.). We
embolden the elpd differences such that zero is not included within the
approximate 90% normal interval over the mean elpd difference estimate (the
$1.64$ standard deviation range). The number of observations in each dataset
is reported by $n$ in the dataset lookup. Our results show that our procedure
always produces models with predictive performance at least as good as
auto.arima, and in three cases identified considerably better models.

| ProjpredARMA | MCMC AutoARIMA |  |  |
---|---|---|---|---|---
Data | $p^{\perp}$ | $q^{\perp}$ | $P^{\perp}$ | $Q^{\perp}$ |  elpd $\pm$ s.e. | $p_{\operatorname{auto}}$ | $q_{\operatorname{auto}}$ | $P_{\operatorname{auto}}$ | $Q_{\operatorname{auto}}$ |  elpd $\pm$ s.e. |  elpd diff. $\pm$ s.e.
H | $3$ | $3$ | $1$ | $0$ | ${}-153{}$ | ${}\pm{}$ | ${}14{}$ | $1$ | $2$ | $0$ | $0$ | ${}-190{}$ | ${}\pm{}$ | ${}15{}$ | ${}\bm{38}{}$ | ${}\bm{\pm}{}$ | ${}\bm{9.3}{}$
I | $3$ | $4$ | $3$ | $2$ | ${}281{}$ | ${}\pm{}$ | ${}12{}$ | $4$ | $1$ | $0$ | $2$ | ${}278{}$ | ${}\pm{}$ | ${}11{}$ | ${}2.8{}$ | ${}\pm{}$ | ${}4.5{}$
J | $3$ | $4$ | $1$ | $0$ | ${}-273{}$ | ${}\pm{}$ | ${}19{}$ | $1$ | $1$ | $0$ | $1$ | ${}-282{}$ | ${}\pm{}$ | ${}22{}$ | ${}8.7{}$ | ${}\pm{}$ | ${}8.4{}$
K | $4$ | $0$ | $2$ | $3$ | ${}-494{}$ | ${}\pm{}$ | ${}11{}$ | $4$ | $1$ | $0$ | $1$ | ${}-492{}$ | ${}\pm{}$ | ${}11{}$ | ${}-1.8{}$ | ${}\pm{}$ | ${}2.4{}$
L | $1$ | $4$ | $2$ | $1$ | ${}-1470{}$ | ${}\pm{}$ | ${}32{}$ | $1$ | $2$ | $2$ | $2$ | ${}-1465{}$ | ${}\pm{}$ | ${}32{}$ | ${}-3.0{}$ | ${}\pm{}$ | ${}4.1{}$

  * •

Dataset lookup: H $=$ Mona Loa $CO_{2}$ ($n=468,\,s=12$), I $=$ corticosteroid
subsidy ($n=204,\,s=12$), J $=$ anti-diabetic drug subsidy ($n=204,\,s=12$), K
$=$ equipment manufacturing ($n=195,\,s=12$), L $=$ daily electricity demand
($n=365,\,s=7$).

Table 3: Seasonal real world data. Predictive performance comparison of
ProjpredARMA with MCMC AutoARIMA on seasonal data sets curated by Hyndman and
Athanasopoulos, (2021). For both procedures, we report the order of the
suggested model as well as the models’ elpd and elpd standard error (s.e.). We
embolden the differences between elpds similarly to as in Table 2. As before,
we list the number of data points in each series under $n$ in the dataset
lookup, as well as the identified seasonality $s$. These results show again
that our procedure produces models at least as good as those chosen by
auto.arima with regards to predictive performance, and twice out-performed the
latter.

We perform this experiment with data from the fpp2 package (Hyndman and
Athanasopoulos,, 2021) since they are well-studied, clean, and we can easily
fit a reference model to them with the default priors previously discussed. If
the data are deemed to be seasonal by auto.arima, we use the seasonal lag $s$
provided therein for our reference models. Similarly we use the same order of
differencing (seasonal and non-seasonal) identified by auto.arima in both
procedures for consistency. We disallow the possibility of drift in the models
identified by auto.arima and concern ourselves only with inference in the
stationary case.

We only consider data where the suggested differencing and seasonalities are
inline with our prior knowledge, and appear to be reasonable upon closer
inspection of the ACF and PACF plots. As is mentioned by Hyndman and
Khandakar, (2008), we favour data where minimal differencing is required for
improved predictive performance. Further, we use data such that the ACF and
PACFs plots lead naturally to reference models through the Box-Jenkins
approach (Box and Jenkins,, 1970) that pass posterior checks when our default
priors are used. It is worth mentioning that seasonality can also be directly
ascertained from the ACF and PACF plots, but we prefer to use them only as
sanity checks. Plots of the data used can be found in Appendix A.

It is important to remember that the results we achieve are entirely dependent
on the reference model used in the projection procedure. Indeed, this allows
the statistician to encode prior beliefs into their model selection procedure
through the reference model, a luxury unavailable in auto.arima.

In Table 2 we present the predictive performance of the models selected by
both procedures across a set of non-seasonal fpp2 datasets. We find that the
difference between the mean elpds is almost always positive (meaning that the
mean elpd of the ProjpredARMA-achieved submodel is almost always higher than
that of the submodel suggested by MCMC AutoARIMA), apart from in datasets A
and E (airline passenger data and Ansett Airline data) where we find that the
difference between elpds is skewed more than 1.64 standard deviations in our
procedure’s favour. We thus find that our procedure is able to achieve
increased stability without sacrificing predictive performance, and in some
cases even out-performs auto.arima.

Moving to seasonal experiments, we tabulate the results in Table 3. In one of
the five examples (Mona Loa $CO_{2}$), projection predictive model selection
was able to find a submodel out-performing that selected by auto.arima. In the
remaining instances, the mean difference between elpds remains mostly positive
(in our favour). Again, we find that our improved stability does not come at
the cost of predictive performance, and indeed some improvement in the latter
is also felt.

###  6.3 Alternative projections

Figure 4: Simulated data. A frequency plot of the selected model orders under
the procedure defined in Algorithm 1 in green and the modification to it
described in Section 6.3 in red. The $100$ series were generated from an
$\operatorname{ARMA}(2,4)$ process. The reference model parameters used,
$p^{\ast}$ and $q^{\ast}$, and the true data-generating parameters
$p^{\text{true}}$ and $q^{\text{true}}$ shown by the dotted black lines. We
find that our procedure as defined in Algorithm 1 correctly handles auto-
regressive and moving-average information, and that the modification described
in Section 6.3 risks over-selecting auto-regressive components and under-
selecting moving-average components. Figure 5: Simulated data. The normal-
approximation $68\%$ elpd confidence intervals of $\operatorname{AR}(p)$
models produced by projection predictive inference are shown in the black
point ranges. The predictive performance of the reference model is seen in the
horizontal dashed yellow line, and model size selected by our search heuristic
from Section 3.3 is shown in the vertical dotted green line. The predictive
performance of the $\operatorname{ARMA}(3,3)$ model found by auto.arima is
seen in the horizontal dashed orange line.

As was previously stated in Section 4.4, it is necessary to separate the
reference model into its AR and MA components before performing projection
predictive inference. Indeed, in Figure 4 we show that when we do not seperate
these components in the reference model, we are liable to vastly over-select
the size of the restricted AR model and slightly under-select the size of the
restricted MA model, since we project too much information onto this
restricted AR component. This results in submodels requiring larger AR
components, which in turn produce noisier residuals and smaller MA components.
We show this by performing the modification to our procedure discussed in
Section 4.4 to $100$ simulated series from an $\operatorname{ARMA}(2,4)$
process. For comparison, we also show the selection frequencies of our
original procedure in Algorithm 1 and the true process values, confirming our
suspicions that such modifications incorrectly incorporate information from
the reference model in submodels.

There are instances, as we have previously discussed, where a projection
directly from an ARMA model to an AR model can reduce computational burden in
inference. What is then of interest to us is to understand how large the AR
order must be to replicate the predictive performance of the
$\operatorname{AR}(\infty)$ model implied by the non-zero MA component of the
ARMA model. To demonstrate this, we presently sample a series from an
$\operatorname{ARMA}(2,3)$ model, and fit to it a reference
$\operatorname{ARMA}(5,5)$ model. We then perform the projection from this
$\operatorname{ARMA}(5,5)$ onto $\operatorname{AR}(p),\,p=1,\dotsc,20$, cross-
validating the models’ performances, and show the results in Figure 5. We
conclude that our $\operatorname{ARMA}(5,5)$ model’s predictive performance
can be replicated by as small a model as an $\operatorname{AR}(4)$. Not only
this, but when we use auto.arima to identify a more parsimonious model for
these data, we are returned an $\operatorname{ARMA}(3,3)$. This model is
neither easy to infer, nor does it match our reference model – which we are
able to do by construction. Thus we find that this $\operatorname{ARMA}$ to
$\operatorname{AR}$ projection can not only out-perform auto.arima, but can do
so at much reduce cost.

###  6.4 Robustness to noise

Figure 6: Simulated data. In the bottom axes are shown the elpds of models
selected by projection predictive model selection in green and by auto.arima
in orange from a simulated $\operatorname{AR}(3)$ model, corrupted with
different levels of noise, $\sigma^{2}$, seen along the $x$-axis. The one
standard error interval of observed mean elpds is also shown, and we note the
scale of the axis. The top axes show the size of AR component identified by
the the submodels at each value of $\sigma^{2}$ (in the top left panel), and
the MA order identified likewise (in the top right). We find that across all
levels of noise, auto.arima fails to recognise the true temporal structure
while projection predictive inference chooses the same (and true) model orders
independently of noise. Further, we see that auto.arima produces models with
significantly worse poor predictive performance than our procedure.

A major advantage to the Bayesian workflow is the ability to manage
uncertainty in data. It is then interesting for us to understand how our two
competing procedures behave when different levels of noise are injected into
the same underlying process.

We produce data from an $\operatorname{AR}(3)$ model with different levels of
noise $\sigma$ injected. Formally, we sample from rcl y_t &= ϕ_1y_t-1 +
ϕ_2y_t-2 + ϕ_3y_t-3 + ε_t
ε_t ∼ N(0,σ^2). We then apply both projection predictive model selection and
auto.arima procedures to identify submodels for each of the series, compute
their predictive performance, and show the results in Figure 6.

We note three important results from this experiment convincing us of the
robustness to noise afforded by our proposed procedure. First we find that
projection predictive inference produces models with elpd at least as good as
those identified by auto.arima under all noise levels.

Second we find that auto.arima is liable to produce ill-performing models,
seen in the instances of very low and uncertain elpd compared to projection
predictive inference. Uncertainty in the selected orders is additionally much
larger compared to our approach, in line with the results in Section 6.1. It
is important to note the scale of the elpd plot; as a heuristic, we might
consider that two elpd values with a difference of less that $4$ are
essentially equivalent, whereas in this plot we show differences of several
orders of magnitude indicating an infeasible model. These results suggest that
some of the models identified by auto.arima may be misspecified given the
data.

Finally, projection predictive model selection is exceptionally stable in
model size, which does not vary greatly among the different levels of
$\sigma^{2}$, whereas auto.arima exhibits large variability when selecting the
moving-average component. Further, we find that projection predictive
inference identifies the true model under all noise regimes.

Such behaviour is expected from projection predictive inference. Indeed,
Piironen and Vehtari, 2016a  discuss how the use of a reference model in model
selection (and more so in projection predictive model selcetion) is able to
filter out noise from the data. We are thus confident in projection predictive
inference’s robustness to noisy data.

###  6.5 The effect of distant lags on submodel selection

Figure 7: Simulated data. The top axes show the values of the true data
generating $\operatorname{AR}(6)$ process parameters $\theta_{p}$ at each lag
$p$, below which are the normal-approximation $68\%$ elpd confidence intervals
of models selected by projection predictive inference. The predictive
performance of the reference model is seen in the dashed yellow line. The
model size selected by our search heuristic from Section 3.3 is shown in the
dotted green line, while that of auto.arima is indicated by the dashed orange
line. This plot shows our model selection heuristic behaviour’s in a simple
case where it finds distant lags to have predictive effect close to zero,
while auto.arima has a tendency to select slightly larger models for
negligible performance gain.

As has been previously discussed, the primary aim of projection predictive
model selection is to find the smallest submodel such that its posterior
predictive performance is not significantly different to that of a reference
model. In the case of additive models, the addition of an additive component
can only ever improve the posterior predictive performance. Thus, if we have
many distant lags with low influence in some true data generating process,
then our procedure will select only the first few necessary to encapsulate the
expressiveness of the reference model and allow it to generalise well to new
data.

One such model is given in Figure 7, in which we have an
$\operatorname{AR}(6)$ where the values of $\theta_{p}$ decrease with $p$
until lag six, after which all lags are equal to exactly zero. We simulate
$N=200$ data points from this process and apply our projection predictive
model selection procedure with an oracle providing an $\operatorname{AR}(6)$
reference model. Consequently, we identify an $\operatorname{AR}(2)$ model as
the smallest submodel with comparable posterior predictive performance, since
the mean reference model elpd lies within one standard error from the mean
submodel elpd (the normal-approximation $68\%$ interval as previously
described). However, using auto.arima to identify a parsimonious submodel,
again with an oracle providing information that the process is a stationary
non-seasonal auto-regressive model, identifies an $\operatorname{AR}(4)$
model. This over-fitting of weak-effect covariates is common in AIC-based
validation, while projection predictive model selection is able to avoid it by
prioritising predictive performance over model complexity and is able to lean
on a reference model in order to do so.

###  6.6 Scalability

We have dealt with examples where the number of observations has been fewer
that $500$. The computation time of Algorithm 1 increases linearly with the
number of data observations. Indeed, even for very large datasets we can still
perform our procedure in reasonable time. It is worth noting that most of this
computational time is spent fitting the AR and MA components with MCMC, and as
a result is dependent on the priors used and machinery chosen. Naturally,
given this Bayesian treatment of order identification our procedure remains
computationally intensive when compared to auto.arima.

##  7 Conclusion

We have motivated an extension of projection predictive model selection to
probabilistic ARMA models and thence developed a novel two-stage Bayesian
order identification procedure. Our procedure was shown through simulated and
real-data experiments to:

  1. 1.

be stable in model selection, and robust to instances of noisy data and
complex data-generating models;

  2. 2.

produce submodels with predictive performance at least as good as auto.arima;

  3. 3.

scale well with increased data size, although it remains computationally
expensive when compared to auto.arima.

Importantly, we have shown how the original idea of Goutis and Robert, (1998)
and Dupuis and Robert, (2003) can be abstracted beyond the realm of
generalised linear models and towards time series applications, wherein
parsimony and model structure is a key concern. In doing so, we have motivated
a robust and efficient ARMA order identification procedure from an information
theory perspective.

##  8 Discussion

Model selection is often motivated as a remedy to over-fitting. However, in a
Bayesian regime the statistician is afforded the luxury of explicating their
prior beliefs. The problem of over-fitting can then be mitigated at least in
part with the use of sensible priors. Instead, we propose our procedure for
use in one of three candidate situations:

  1. 1.

we have a rich model which is good for prediction, but would like to reduce
its size to reduce the computational burden or to improve robustness to
changes in the data-generating distribution;

  2. 2.

we would like to identify predictive submodels to gain a better understanding
of important temporal correlation structures;

  3. 3.

we have an ARMA model with moving-average components but would like to convert
it to a purely auto-regressive model, or more generally we would like to
investigate models with similar predictive performance but different
structures.

Importantly, we do not advocate for model selection in a Bayesian setting for
its own sake. Rather we believe that model selection within a reference model
paradigm can provide the statistician improved interpretation of the
underlying data-generating process, and the opportunity to reduce
computational cost without sacrificing predictive performance. As such, the
creation of the reference model and the use of priors therein is of critical
importance to a good analysis.

The choice of reference model is not unambiguous in general, and the results
of our procedure may vary with the ability of the prior to distinguish noise
from signal in the data (Kohns and Potjagailo,, 2022).

The procedure presented in this paper has dealt with the ARMA model given its
prevalence in literature and empirically proven strength in practice. We
believe that other time series models, notably state space and Gaussian
process time series models could both also benefit from projection predictive
model selection. Indeed Piironen and Vehtari, 2016b  have previously shown
that it is able to deal with Gaussian process model selection, and Catalina et
al., (2022) have likewise shown its extension to generalised additive models.

We have also discussed how there may arise situations in which directly using
standard sparsifying priors such as the horseshoe and spike-and-slab priors
(Piironen and Vehtari,, 2017; Mitchell and Beauchamp,, 1988) for ARMA variable
selection may result in invalid or unintuitive models. Further investigation
into the sparsifying time series priors, possibly similar to the R2D2 prior of
Zhang et al., (2022), the $L_{1}$ ball prior of Xu and Duan, (2020), or fused
lasso by Casella et al., (2010) would surely benefit the field.

### Acknowledgements

We acknowledge the computational resources provided by the Aalto Science-IT
project. This paper was partially funced by the Research Council of Finland
Flagship programme: Finnish Center for Artificial Intelligence, and Research
Council of Finland project “Safe iterative model building” (340721).

## References

  * Arlot and Celisse, (2010) Arlot, S. and Celisse, A. (2010).  A survey of cross-validation procedures for model selection.  Statistics surveys, 4:40–79.
  * Bańbura et al., (2010) Bańbura, M., Giannone, D., and Reichlin, L. (2010).  Large Bayesian vector auto regressions.  Journal of applied Econometrics, 25(1):71–92.
  * Barbieri and Berger, (2004) Barbieri, M. M. and Berger, J. O. (2004).  Optimal Predictive Model Selection.  The Annals of Statistics, 32(3):870–897.
  * Bhadra et al., (2019) Bhadra, A., Datta, J., Polson, N. G., and Willard, B. (2019).  Lasso meets horseshoe: A survey.  Statistical Science, 34(3):405–427.
  * Box and Jenkins, (1970) Box, G. and Jenkins, G. M. (1970).  Time Series Analysis: Forecasting and Control.  Holden-Day.
  * Bürkner et al., (2020) Bürkner, P.-C., Gabry, J., and Vehtari, A. (2020).  Approximate leave-future-out cross-validation for Bayesian time series models.  Journal of Statistical Computation and Simulation, 90(14):2499–2523.
  * Carpenter et al., (2017) Carpenter, B., Gelman, A., Hoffman, M. D., Lee, D., Goodrich, B., Betancourt, M., Brubaker, M., Guo, J., Li, P., and Riddell, A. (2017).  Stan: A probabilistic programming language.  Journal of statistical software, 76(1).
  * Carriero et al., (2015) Carriero, A., Clark, T. E., and Marcellino, M. (2015).  Bayesian VARs: specification choices and forecast accuracy.  Journal of Applied Econometrics, 30(1):46–73.
  * Casella et al., (2010) Casella, G., Ghosh, M., Gill, J., and Kyung, M. (2010).  Penalized regression, standard errors, and bayesian lassos.  Bayesian analysis, 5(2):369–411.
  * Catalina et al., (2022) Catalina, A., Burkner, P.-C., and Vehtari, A. (2022).  Projection predictive inference for generalized linear and additive multilevel models.  In AISTATS.
  * Catalina et al., (2021) Catalina, A., Bürkner, P., and Vehtari, A. (2021).  Latent space projection predictive inference.
  * Chan, (2021) Chan, J. C. (2021).  Minnesota-type adaptive hierarchical priors for large Bayesian VARs.  International Journal of Forecasting, 37(3):1212–1226.
  * Chan et al., (2016) Chan, J. C., Eisenstat, E., and Koop, G. (2016).  Large Bayesian VARMAs.  Journal of Econometrics, 192(2):374–390.
  * Chan, (2013) Chan, J. C. C. (2013).  Moving average stochastic volatility models with application to inflation forecast.  Journal of Econometrics, 176(2):162–172.
  * Chan and Chen, (2011) Chan, K.-S. and Chen, K. (2011).  Subset ARMA selection via the adaptive lasso.  Statistics and Its Interface, 4:197–205.
  * Chib and Greenberg, (1994) Chib, S. and Greenberg, E. (1994).  Bayes inference in regression models with ARMA (p, q) errors.  Journal of Econometrics, 64(1-2):183–206.
  * Dupuis and Robert, (2003) Dupuis, J. and Robert, C. (2003).  Variable selection in qualitative models via an entropic explanatory power.  Journal of Statistical Planning and Inference, 111:77–94.
  * Efron, (2008) Efron, B. (2008).  Microarrays, Empirical Bayes and the Two-Groups Model.  Statistical Science, 23(1).  arXiv:0808.0572 [stat].
  * Efron, (2010) Efron, B. (2010).  Large-Scale Inference: Empirical Bayes Methods for Estimation, Testing, and Prediction.  Cambridge University Press, 1 edition.
  * Gabry et al., (2019) Gabry, J., Simpson, D., Vehtari, A., Betancourt, M., and Gelman, A. (2019).  Visualization in Bayesian workflow.  Journal of the Royal Statistical Society: Series A (Statistics in Society), 182(2):389–402.  tex.ids= gabryVisualizationBayesianWorkflow2019a arXiv: 1709.01449.
  * Geisser and Eddy, (1979) Geisser, S. and Eddy, W. F. (1979).  A Predictive Approach to Model Selection.  Journal of the American Statistical Association, 74(365):153–160.
  * Gelman et al., (1996) Gelman, A., Meng, X.-L., and Stern, H. (1996).  Posterior predictive assessment of model fitness via realized discrepancies.  Statistica sinica, pages 733–760.
  * Gelman et al., (2020) Gelman, A., Vehtari, A., Simpson, D., Margossian, C. C., Carpenter, B., Yao, Y., Kennedy, L., Gabry, J., Bürkner, P.-C., and Modrák, M. (2020).  Bayesian workflow.
  * Giacomini, (2013) Giacomini, R. (2013).  The relationship between dsge and var models.  VAR Models in Macroeconomics–New Developments and Applications: Essays in Honor of Christopher A. Sims.
  * Giannone et al., (2015) Giannone, D., Lenza, M., and Primiceri, G. E. (2015).  Prior Selection for Vector Autoregressions.  The Review of Economics and Statistics, 97(2):436–451.  _eprint: https://direct.mit.edu/rest/article-pdf/97/2/436/1917922/rest_a_00483.pdf.
  * Gneiting and Raftery, (2007) Gneiting, T. and Raftery, A. E. (2007).  Strictly proper scoring rules, prediction, and estimation.  Journal of the American statistical Association, 102(477):359–378.
  * Goutis and Robert, (1998) Goutis, C. and Robert, C. (1998).  Model choice in generalised linear models: A Bayesian approach via Kullback-Leibler projections.  Biometrika, 85:29–37.
  * Granger and Morris, (1976) Granger, C. W. and Morris, M. J. (1976).  Time series modelling and interpretation.  Journal of the Royal Statistical Society: Series A (General), 139(2):246–257.
  * Hahn and Carvalho, (2015) Hahn, P. R. and Carvalho, C. M. (2015).  Decoupling shrinkage and selection in Bayesian linear models: a posterior summary perspective.  Journal of the American Statistical Association, 110(509):435–448.
  * Hannan and Rissanen, (1982) Hannan, E. J. and Rissanen, J. (1982).  Recursive estimation of mixed autoregressive-moving average order.  Biometrika, 69(1):81–94.
  * Harvey, (1985) Harvey, A. C. (1985).  Trends and cycles in macroeconomic time series.  Journal of Business & Economic Statistics, 3(3):216–227.
  * Hoffman and Gelman, (2011) Hoffman, M. D. and Gelman, A. (2011).  The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo.  arXiv:1111.4246 [cs, stat].
  * Hyndman and Athanasopoulos, (2021) Hyndman, R. and Athanasopoulos, G. (2021).  Forecasting: Principles and Practice.  OTexts, Australia, 3rd edition.
  * Hyndman and Khandakar, (2008) Hyndman, R. J. and Khandakar, Y. (2008).  Automatic time series forecasting: The forecast package for R.  Journal of Statistical Software, 27(3):1–22.
  * Ives et al., (2010) Ives, A. R., Abbott, K. C., and Ziebarth, N. L. (2010).  Analysis of ecological time series with arma (p, q) models.  Ecology, 91(3):858–871.
  * Johnstone and Silverman, (2004) Johnstone, I. M. and Silverman, B. W. (2004).  Needles and straw in haystacks: Empirical Bayes estimates of possibly sparse sequences.  The Annals of Statistics, 32(4).  arXiv:math/0410088.
  * Kapetanios, (2003) Kapetanios, G. (2003).  A note on an iterative least-squares estimation method for ARMA and VARMA models.  Economics Letters, 79(3):305–312.
  * Kim et al., (1999) Kim, C.-J., Nelson, C. R., et al. (1999).  State-space models with regime switching: classical and Gibbs-sampling approaches with applications.  MIT Press Books, 1.
  * Kohns and Potjagailo, (2022) Kohns, D. and Potjagailo, G. (2022).  A flexible Bayesian MIDAS approach for interpretable nowcasting and forecasting.
  * Koop, (2013) Koop, G. M. (2013).  Forecasting with medium and large Bayesian VARs.  Journal of Applied Econometrics, 28(2):177–203.
  * Matamoros and Torres, (2021) Matamoros, A. A. and Torres, C. C. (2021).  bayesforecast: Bayesian time series modeling with Stan.  ArxiV Preprint.
  * McQuarrie and Tsai, (1998) McQuarrie, A. D. R. and Tsai, C.-L. (1998).  Regression and Time Series Model Selection.  World Scientific.
  * Mitchell and Beauchamp, (1988) Mitchell, T. J. and Beauchamp, J. J. (1988).  Bayesian variable selection in linear regression.  Journal of the American Statistical Association, 83(404):1023–1032.
  * Mogliani and Simoni, (2021) Mogliani, M. and Simoni, A. (2021).  Bayesian midas penalized regressions: estimation, selection, and prediction.  Journal of Econometrics, 222(1):833–860.
  * Nardi and Rinaldo, (2011) Nardi, Y. and Rinaldo, A. (2011).  Autoregressive process modeling via the Lasso procedure.  Journal of Multivariate Analysis, 102(3):528–549.
  * Neal et al., (2011) Neal, R. M. et al. (2011).  Mcmc using Hamiltonian dynamics.  Handbook of markov chain monte carlo, 2(11):2.
  * Ng and Perron, (1995) Ng, S. and Perron, P. (1995).  Unit root tests in ARMA models with data-dependent methods for the selection of the truncation lag.  Journal of the American Statistical Association, 90(429):268–281.
  * Pavone et al., (2022) Pavone, F., Piironen, J., Bürkner, P.-C., and Vehtari, A. (2022).  Using reference models in variable selection.  Computational Statistics.
  * Piironen et al., (2022) Piironen, J., Paasiniemi, M., Catalina, A., Weber, F., and Vehtari, A. (2022).  projpred: Projection predictive feature selection.  R package version 2.2.1.
  * Piironen et al., (2018) Piironen, J., Paasiniemi, M., and Vehtari, A. (2018).  Projective inference in high-dimensional problems: Prediction and feature selection.  Electronic Journal of Statistics, 14(1).
  * (51) Piironen, J. and Vehtari, A. (2016a).  Comparison of Bayesian predictive methods for model selection.  Statistics and Computing, 27(3):711–735.
  * (52) Piironen, J. and Vehtari, A. (2016b).  Projection predictive model selection for Gaussian processes.  In 2016 IEEE 26th International Workshop on Machine Learning for Signal Processing (MLSP). IEEE.
  * Piironen and Vehtari, (2017) Piironen, J. and Vehtari, A. (2017).  Sparsity information and regularization in the horseshoe and other shrinkage priors.  Electronic Journal of Statistics, 11(2).
  * Polson and Scott, (2012) Polson, N. G. and Scott, J. G. (2012).  On the half-Cauchy prior for a global scale parameter.  Bayesian Analysis, 7(4):887–902.
  * Ray and Bhattacharya, (2018) Ray, P. and Bhattacharya, A. (2018).  Signal adaptive variable selector for the horseshoe prior.  arXiv preprint arXiv:1810.09004.
  * Schotman and Van Dijk, (1991) Schotman, P. C. and Van Dijk, H. K. (1991).  On bayesian routes to unit roots.  Journal of Applied Econometrics, 6(4):387–401.
  * Schwarz, (1978) Schwarz, G. (1978).  Estimating the Dimension of a Model.  The Annals of Statistics, 6(2):461 – 464.  Publisher: Institute of Mathematical Statistics.
  * Shao, (1997) Shao, J. (1997).  An asymptotic theory for linear model selection.  Statistica sinica, pages 221–242.
  * Sims and Uhlig, (1991) Sims, C. and Uhlig, H. (1991).  Understanding unit rooters: A helicopter tour.  Econometrica, 59:1591–99.
  * Sims, (1988) Sims, C. A. (1988).  Bayesian skepticism on unit root econometrics.  Journal of Economic Dynamics and Control, 12(2):463–474.
  * Sivula et al., (2020) Sivula, T., Magnusson, M., and Vehtari, A. (2020).  Uncertainty in bayesian leave-one-out cross-validation based model comparison.  arXiv preprint arXiv:2008.10296.
  * Spiegelhalter et al., (2002) Spiegelhalter, D. J., Best, N. G., Carlin, B. P., and van der Linde, A. (2002).  Bayesian measures of model complexity and fit.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 64(4):583–639.
  * Stock and Watson, (2007) Stock, J. H. and Watson, M. W. (2007).  Why has us inflation become harder to forecast?  Journal of Money, Credit and banking, 39:3–33.
  * Stone, (1977) Stone, M. (1977).  An asymptotic equivalence of choice of model by cross-validation and Akaike’s criterion.  Journal of the Royal Statistical Society: Series B (Methodological), 39(1):44–47.
  * Vehtari et al., (2016) Vehtari, A., Gelman, A., and Gabry, J. (2016).  Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC.  Statistics and Computing, 27(5):1413–1432.
  * Vehtari and Ojanen, (2012) Vehtari, A. and Ojanen, J. (2012).  A survey of Bayesian predictive methods for model assessment, selection and comparison.  Statistics Surveys, 6(none):142 – 228.
  * Watanabe, (2013) Watanabe, S. (2013).  A Widely Applicable Bayesian Information Criterion.  Journal of Machine Learning Research, page 31.
  * Wold, (1938) Wold, H. (1938).  A study in the analysis of stationary time series.  PhD thesis, Almqvist & Wiksell.
  * Xu and Duan, (2020) Xu, M. and Duan, L. L. (2020).  Bayesian inference with the l1-ball prior: Solving combinatorial problems with exact zeros.
  * Zhang et al., (2020) Zhang, B., Chan, J. C., and Cross, J. L. (2020).  Stochastic volatility models with ARMA innovations: An application to G7 inflation forecasts.  International Journal of Forecasting, 36(4):1318–1328.
  * Zhang et al., (2022) Zhang, Y. D., Naughton, B. P., Bondell, H. D., and Reich, B. J. (2022).  Bayesian Regression Using a Prior on the Model Fit: The R2-D2 Shrinkage Prior.  Journal of the American Statistical Association, 117(538):862–874.

##  Appendix A Datasets

We present below the raw data used in the experiments of Section 6.2 in Figure
8, and in Figure 9 we show them after applying the differencing suggested by
auto.arima.

Figure 8: Datasets used in Section 6.2. In the lookup below, we also note the
length of the time series by $n$ and the seasonality identified by auto.arima
by $s$.
Dataset lookup: A $=$ airline passengers ($n=47$), B $=$ international
visitors ($n=36$), C $=$ Lake Huron bathymetry ($n=98$), D $=$ insurance
quotes ($n=40$), E $=$ Ansett Airline passengers ($n=269$), F $=$ maximum
annual temperature ($n=46$), G $=$ female murder rate ($n=55$), H $=$ Mona Loa
$CO_{2}$ ($n=468,\,s=12$), I $=$ corticosteroid subsidy ($n=204,\,s=12$), J
$=$ anti-diabetic drug subsidy ($n=204,\,s=12$), K $=$ equipment manufacturing
($n=195,\,s=12$), L $=$ daily electricity demand ($n=365,\,s=7$). Figure 9:
Differenced datasets (denoted $\Delta y_{t}$) used in Section 6.2. We record
the order of non-seasonal and seasonal differencing denoted $d$ and $D$
respectively.
Dataset lookup: A $=$ airline passengers ($d=2$), B $=$ international visitors
($d=1$), C $=$ Lake Huron bathymetry ($d=1$), D $=$ insurance quotes ($d=0$),
E $=$ Ansett Airline passengers ($d=1$), F $=$ maximum annual temperature
($d=1$), G $=$ female murder rate ($d=2$), H $=$ Mona Loa $CO_{2}$
($d=1,\,D=1$), I $=$ corticosteroid subsidy ($d=1,\,D=1$), J $=$ anti-diabetic
drug subsidy ($d=1,\,D=1$), K $=$ equipment manufacturing ($d=0,\,D=1$), L $=$
daily electricity demand ($d=0,\,D=1$).