Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
762,340 characters · 11 sections · 72 citation commands
The ARR2 prior: flexible predictive prior definition for Bayesian auto-regressions
\thispagestyle{empty}
Independent priors over model components have the undesirable tendency to inflate explained variance with increasing model complexity. As a remedy to this, the recent literature has advocated for the use of joint priors over model components. Most notably, the R2D2 prior zhang_bayesian_2022 allows the modeller to encode some prior belief on the model's performance as defined by the coefficient of determination, $R^2$, which then trickles down to the parameter level through a mapping between the two. This idea has since been extended to multi-level models aguilar_intuitive_2023, generalised linear models yanchenko_r2d2_2021 and spatial models yanchenko_r2d2_2023. While these priors have found rich application, it is not clear whether and how they can be applied to time-series settings. The auto-covariance structures defined by lagged observations and latent dynamics create dependencies between the model parameters that govern auto-correlations and the $R^2$. This necessitates adaptation of the previously presented $R^2$ prior frameworks. The kernel of this work is an extension of the joint shrinkage framework over $R^2$ to the auto-regressive (AR) components of Bayesian time-series models, including latent state-space models. Concretely, in this paper we:
In Section (ref), we begin by discussing relevant background and previous work. In Section (ref) we motivate and define the $R^2$ for $\operatorname{AR}$ models, and establish connections to the previous literature. Section (ref) in particular summarises the properites of the ARR2 prior compared to popular alternative shrinkage priors for $\operatorname{AR}$ models. In Section (ref) we adapt the ARR2 prior to ARX and state-space models. We then compare performance in simulated data experiments in Section (ref), and a real-world data forecasting exercise to US inflation in Section (ref). Finally, we conclude in Section (ref) with some recommendations and discussion of promising future research directions.
For a model with many regression coefficients, independent priors, which do not encode dependence between coefficients, can imply a high degree of explained variance gelmanRegressionOtherStories2020. In terms of $R^2$, this is exhibited as a prior concentrated near the upper bound of $1$, indicating that the model is expected to almost perfectly fit the data. If the likelihood does not dominate the prior, this can lead to posterior estimates of explained variance that overestimate the actual explanatory power of the model at the population level\footnote{This can be seen by in-sample $R^2$ being much higher than out-of-sample $R^2$.}, negatively impacting both predictive performance and inference.
Imposing sparsity with discrete mixture priors ishwaran2005spike or shrinkage with continuous priors carvalho_handling_2009,bhattacharya2015dirichlet,piironen_sparsity_2017 may limit these tendencies by regularising the variance of the predictor term. Notably, the R2D2 prior zhang_bayesian_2022 takes the approach of directly encoding prior beliefs about the $R^2$ as a function of total variance and residual variance. The prior total variance explained is then decomposed and assigned to the additive components of the model through a simplex prior. By choosing appropriate priors over the simplex, one may also induce sparsity-favouring patterns in coefficients, which is separate from the information encoded by the prior on $R^2$.
A similar prior framework is proposed by fuglstad_intuitive_2020 in which a variance parameter is shared by all coefficients across model components that is either divided according to a simplex prior or regularised via a penalised complexity prior simpson2017penalising. Penalised complexity priors have been extended to AR models by sorbye_penalised_2017. This approach specifies priors on partial autocorrelations of the AR process, and does not directly set a prior on $R^2$. Such penalised complexity priors are not straightforwardly extended to the different time-series considered in this paper. As such, we will leave any comparisons to such priors for future investigation.
The previous time-series literature has remained wary of the usefulness of $R^2$ as a measure of model fit. This is due to the non-standard frequency distribution of $R^2$ under non-stationarity phillips1986understanding, and the often condemned property of $R^2$, conditional on parameters and design, to increase monotonically with the number of covariates included. Here though, we will assume stationarity and are interested in finite sample properties under different priors heaps_enforcing_2022. Indeed, we will show that popular priors used in the time-series literature can be extremely informative on the $R^2$ space, despite having no explicit interest in the quantity\footnote{kovalBayesianReconciliationReturn2024 show however how the Bayesian $R^2$ is a useful metric to describe stock return explainability.}. As such, even if one does not believe that they should be reasoning about $R^2$, they may find themselves actually doing so unknowingly. Our proposed prior, on the other hand, makes prior knowledge on the $R^2$ of an auto-regressive model explicit, while benefiting from the same shrinkage properties usually sought out in previous work.
The priors suggested in this paper also differ substantially from the default or Jeffreys priors which have been considered for simple time-series model in the past zellner1996introduction,berger1994noninformative. liseo2013objective recommend not using such priors for AR processes of order larger than 4, so we do not entertain comparison to these priors any further.
In the following, we will denote the population variance operator $\var{x}$ by $\sigma^2_{x}$, its conditional variant $\var{x\mid\theta}$ by $\sigma^2_{x\mid \theta}$, and their unbiased sample estimators by $\hat\sigma^2_{x}$ and $ \hat{\sigma}^2_{x\mid \theta}$ respectively. Further, we denote by $\theta$ the vector of all parameters entering the observation model. Say, $y_t = \phi y_{t-1} + \epsilon_t$ and $\epsilon_t \sim \normal(0,\sigma^2)$, assuming weak stationarity, then $\theta=(\phi,\sigma^2)$, and $\sigma^2_{y|\theta} = \frac{\sigma^2}{1-\phi^2}$. Further, we suppress the conditioning set of parameters and data when using the “$\sim$” symbol for readability.
Consider a pure auto-regressive time-series model of order $p$, denoted by $\operatorname{AR}(p)$, as follows
where $y_t$ is the observed value at time $t$, $\phi_i$ are the regression coefficients, and $\epsilon_t$ is white noise with variance $\sigma^2>0$. In the following, denote by $\mu_t$ the term of the linear model independent of the contemporaneous white noise term, $\mu_t=\sum_{i=1}^p \phi_i y_{t-i}$.
Assume that $\phi$ is in the stationary region $\phi \in M^{p \times 1} \subseteq \mathbb{R}^{p\times1}$ such that the roots of the characteristic polynomial $\phi(u) = 1 -\phi_1u-\dotsc-\phi_pu^{p}, u \in \mathbb{C} $ lie outside the unit circle hamilton2020time. This implies a stable, weakly stationary time-series with $\vert\E{y_t}\vert < \infty$ and $\sigma_{y_t}^2 < \infty$. Without loss of generality, we will assume throughout that the time-series process is centred on zero. The fraction of variance that cannot be explained by our model right before time-point $t$ is $\sigma^2/\var{y_t}$. Similar to nelson_interpretation_1976, we define $R^2$ for auto-regressive models as
where the final equality stems from the white noise definition, implying that current white noise and lags are independent. Hence, setting a prior on the model parameters coherent with a prior on $R^2$ necessitates a joint prior on the variance of the predictor term and of the noise variance.
We now derive the hierarchies and distributional assumptions needed on the model parameters in order to induce a $\betadist(\mu_{R^2},\varphi_{R^2})$ for $R^2$. Throughout, we parameterise the beta distribution in terms of a location $\mu_{R^2}$ and precision $\varphi_{R^2}$ for convenience of interpretation aguilar_intuitive_2023. Consider a multivariate normal prior centred on the zero vector with covariance $\Lambda = \text{diag}(\lambda_1^2,\dots,\lambda_p^2)$
Denote by $y_{-p} = (y_{t-1} \cdots y_{t-p})^{\prime}$ the vector of $p$ lagged observations and by $\theta$ a vector of all parameters of the observation model $\theta = (\phi,\sigma^2)$. Then, the variance expression of the predictor term is
where $\Sigma_y$ is the lag covariance matrix with diagonal $(\E{y_{t-1}y_{t-1}\vert\theta},\dotsc,\E{y_{t-p}y_{t-p}\vert\theta})$. Equation (ref) follows from the law of total variance and the second term vanishes due to the prior mean of zero. Equation (ref) follows from the fact that the variance term is a scalar. Due to the assumption of weak stationarity of the AR process, the conditional expectation functions are the same for all lags $i \in \{1,\dotsc,p \}$ and given by
The auto-covariance function is $\gamma(k) = \text{cov}(y_t, y_{t-k}\vert \theta) = \E{y_{t}y_{t-k} \mid \theta}$, $\{\gamma(k) \in \mathbb{R} : k \in \mathbb{Z}\}$, where $k$ stands for the order. Again, due to the weak stationarity assumption, these are only functions of $k$ (the time distance between lags) and not $t$ (the time index itself).\footnote{The set of all auto-covariances are described by a set of homogeneous difference equations, whose system is named after yule_vii_1927 and walker_periodicity_1931.} The total variance of the predictor term in Equation (ref) is therefore the scaled sum of prior variances
The $R^2$ in Equation (ref) reduces to
By a change of variables, a beta prior on $R^2$ then implies a generalised beta-prime distribution johnson1995continuous for the sum of scaled prior variances of the $\operatorname{AR}$ coefficients
The GBP distribution can be obtained through a transformation of a beta-prime distributed variable, so if $x \sim \text{BP}(a,b)$, then $y = d x^{1/c} \sim \text{GBP}(a,b,c,d)$ for $a,b,c,d > 0$ johnson1995continuous. Therefore, an equivalent prior can be found by scaling the prior variance of $\phi_i$ by $\sigma^2 / \sigma^2_{y_t\mid \theta}$, implying $R^2 = \tau^2/(\tau^2+1)$ and the prior for $\tau^2$ reduces to the simpler beta-prime distribution johnson1995continuous, $\tau^2 \sim \text{BP}(\mu_{R^2},\varphi_{R^2})$. Hence, prior (ref) needs to be scaled by $\sigma^2/\sigma^2_{y_t}$ in order for the implied prior in $R^2$ to remain non-constant with respect to changes in $\sigma^2$ or $\sigma^2_{y_t \mid \theta}$.
In a final step, we allow for adaptivity of shrinkage at the coefficient level, as per the global-local prior framework polson2010shrink, and decompose the sum of prior scales as
with $\psi_i\geq0$ lying on the probabilistic simplex so that $\sum_{i=1}^p\psi_i = 1$. We posit the natural Dirichlet prior over the vector $\psi$ following zhang_bayesian_2022. Each $\psi_i$ determines the fraction of total variance allocated to the $i\text{th}$ lag. We summarise the model hierarchy in Definition (ref).
To make inference on the total variance term $\tau^2$, we formulate the equivalent prior in $R^2$ space and transform back appropriately via Equation (ref).
To avoid $\sigma^2_{y_t \mid \theta}$ being dependent on parameters since this can complicate the posterior geometry, we propose the simplification of replacing $\sigma^2_{y_t \mid \theta}$ with the single data-based estimate $\hat{\sigma}^2_{y}$. This is common practice within the previous literature for $R^2$ priors aguilar_intuitive_2023,aguilar2024generalized. In fact, $\hat{\sigma}^2_{y}$ can be shown to be a conservative estimate of the conditional variance, since $ \sigma^2_{y} = \mathbb{E}[\sigma^2_{y_t\vert \theta}]+ \var{\mathbb{E}[y_t\vert \theta]}$ $\implies \sigma^2_{y}\geq \mathbb{E}[\sigma^2_{y_t\vert \theta}]$. From Equation (ref) one can see that this will tend to bias $R^2$ slightly upwards. Yet, shrinkage properties of the prior on $\phi$, that is, the behaviour of the prior near the origin and the tails polson2010shrink, remain unaltered by this modification. These are controlled by the GBP parameters $ac$ and $bc$ respectively yanchenko_r2d2_2021. Additionally, we expect the data-based estimate to be reasonable, as the term $\E{y_t\vert \theta}$ as defined by Equation (ref) is 0.
For the ARR2 prior, the hyperparameters controlling the beta prior on the $R^2$ and the Dirichlet prior on $\psi$ need to be specified by the modeller. The prior on the $R^2$ encodes beliefs about the variance explained by the model, while the prior on $\psi$ encodes beliefs about the contribution of each lag to the explained variance, and the degree of sparsity. Their hyperparameters alter the prior's properties for $\phi$ at the origin as well as the tails zhang_bayesian_2022. zhang_bayesian_2022 use a uniform prior on $R^2$, induced by $(\mu_{R^2}, \varphi_{R^2})$ = (0.5, 2), while aguilar_intuitive_2023 instead recommend $ (\mu_{R^2},\varphi_{R^2})=(0.5,1)$. They show that setting $(\mu_{R^2},\varphi_{R^2})$ such that $(1-\mu_{R^2})\varphi_{R^2} \in [0,0.5]$ induces fat enough tails so that the marginal prior for $\phi$ is of bounded influence, meaning not shrinking sufficiently large signals. This additionally causes the $R^2$ prior to exhibit a bathtub-like shape with relatively high mass at 0 and 1. Large mass at 1, however, is undesirable for stationary AR models, since the stationary region for $\phi$ implies $R^2 < 1$\footnote{Suppose $\phi \in M^{p\times 1}$ and $\sigma^2 > 0$. Then, by the properties of weak stationarity, $\text{var}(y_{-p}^{\prime}\phi) < \infty$. This renders $\text{var}(y_{-p}^{\prime}\phi)/ (\text{var}(y_{-p}^{\prime}\phi) + \sigma^2) < 1$.}. It is important to recognise that non-stationary processes may occur even when \(R^2 < 1\), particularly in finite samples, and the decomposition defined by the Dirichlet prior will also impact the prior probability of stationarity. Supplementary Material Section A shows that for the ARR2 prior, the majority of the prior mass lies in the stationary region, however this can be influenced by the hyperparameters of the \(R^2\) prior and its simplex decomposition. Decompositions that heavily regularise higher order lag polynomials to zero reduce the probability of prior non-stationarity since lower dimensional $\operatorname{AR}$ models have less complex restrictions on the stationary parameter space huerta1999priors,heaps_enforcing_2022.
As a default, we prefer $R^2 \sim \betadist(1/3,3)$ which exerts shrinkage toward lower $R^2$ values, yet has a small, nearly constant gradient. This implies a long tail with reasonable mass on larger $R^2$ values. Based on our experiments, this is a sensible default, but potential prior-likelihood conflict may be checked, for example, with power-scaling sensitivity analysis kallioinenDetectingDiagnosingPrior2023.Supplementary Material Section B shows a plot of our recommended $R^2$ prior along with two alternatives.
As discussed by zhang_bayesian_2022, the larger the prior concentrations $\xi_i$, the more the Dirichlet distribution will resemble a uniform $\psi \approx (1/p,\dots,1/p)$, whereas small concentration values lead to more mass in the edges of the simplex, inducing stronger sparsity patterns in $\phi$. aguilar_intuitive_2023 recommends setting $\xi_i=0.5$ to encourage sparse posteriors\footnote{In practice, we have found $\xi_i =1$ to yield clear sparsity patterns while remaining computationally feasible.}. By setting these hyperparameters appropriately, the variance decompositions can resemble priors commonly used in time-series modelling (Table (ref)). Importantly, unlike these other priors, the ARR2 allows independently changing the prior on the $R^2$ without affecting this decomposition (see Figure (ref)). We consider two types of decompositions motivated by priors popular in time-series analysis, discussed below.
The Minnesota prior, often used for economic data, was originally proposed by doan1984forecasting and is motivated by the finding that stationary time-series exhibit auto-correlations which are often well described by exponential decay. We consider the hierarchical version presented by carriero_bayesian_2015 which may be viewed as a non-hierarchical normal-gamma prior brown2010inference. Increasing shrinkage with each lag polynomial favours initial lag components to contribute most to total variance.
The regularised horseshoe prior piironen_hyperprior_2017 on the other hand is a more general sparsity-inducing prior which has been recently adopted for many time-series models in economics chan_minnesota-type_2021,huber2021inducing,kohns2024horseshoe. It belongs to the class of global-local priors polson_half-cauchy_2012 in which fat tails favour either shrinking a coefficient strongly toward zero, or only very little. Following piironen_hyperprior_2017, we scale the global scale $\tau_{\text{RHS}}$ by the hyperparameter $\tau_0$, defined by a prior on number of active coefficients (set here to half the number of lags included in the model). Additionally, a default choice of independent Gaussians with unit variance is commonly used as a weakly informative prior in linear regression models, and we compare to this as a baseline.
Figure (ref) shows in the upper panel the induced prior for $R^2$, and in the lower, relative $R^2$ contributions\footnote{Conditional on the parameters, the conditional $R^2$ for the $i\text{th}$ lag can also be shown to reduce to the squared $i\text{th}$-degree partial auto-correlation (see Supplementary Material Section F.2). In the frequentist treatment this is used to choose the appropriate AR order box_time_1994.} from each lag polynomial which we denote by
As shown, with appropriate hyperparameter settings, the ARR2 prior can resemble the other priors, or can encode different beliefs entirely. In the simulated data experiments (Section (ref)) and case study (Section (ref)), we will set the prior concentration values to generate comparable shapes to the Minnesota ($\xi_i = (p^2/10 \cdot 1/i^2)$) and RHS ($\xi = (0.1,\dots,0.1)$).\footnote{Scaling the concentration values for the Minnesota ARR2 by $p^2/10$ enforces a lower bound of 0.1, which we have found empirically to be a threshold under which computational issues arise with sampling algorithms.}
The previous discussion has focused on the prior properties in $R^2$-space. It is also possible to analyse the induced priors for the AR coefficients, the partial autocorrelations of the AR process, and the roots of the characteristic polynomial. This is shown in Supplementary Material Section A. Also in these spaces, the ARR2 (Minn.) mostly resembles the Minnesota prior and ARR2 (flat) the RHS, respectively.
For simplicity, such analyses consider the observation noise variance to be known and fixed. With respect to the AR coefficients, we conjecture that the marginal properties of our prior are similar to those discussed in aguilar_intuitive_2023 under the additional assumption of fixing $\sigma^2_{y_t \mid \theta}$ to a known scalar. We refer the reader to aguilar_intuitive_2023's work for further comparison to other popular shrinkage priors.
\FloatBarrier
Auto-regressive models with exogenous covariates (ARX models) are an extension to pure AR models which allow the modeller to incorporate exogenous information (independent of the lagged covariates). They are used widely in economic analysis as building blocks for multivariate extensions to AR models litterman1986forecasting, and for monitoring in various fields including building and structure engineering saitoBayesianModelSelection2010,matsuokaBayesianEstimationInstantaneous2021,matsuokaBayesianTimeFrequency2020,kimLongtermBridgeHealth2018,barraza-barrazaAdaptiveARXModel2017, medicine fangBayesianInferenceFederated2021,nunesARXModelingDrug2013, and environmental sciences zanottiChoosingLinearNonlinear2019,zhangEstimatingDynamicSolar2022. ARX models are defined formally as
where $X \in \mathbb{R}^{T\times m}$ is the exogenous design matrix, $\beta$ is an $m\times 1$ vector of parameters and $\epsilon_t$ is some white-noise term with variance $\sigma^2$.
Here, the predictor term $\mu_t = \sum_{i=1}^p\phi_iy_{t-i} + x_t^{\prime}\beta$, now defines the total variance $\sigma^2_{\mu_t} = \var{\phi^{\prime}y_{-p}} + \var{x_t^{\prime} \beta}$ which we can now decompose as in Section (ref). Assume that the covariates $X$ are scaled to 0 mean and unit variance. Since $\E{x_t^{\prime}\beta y_{t-k}} = 0$ for all lag polynomials, it is easy to verify that $\var{\phi^{\prime} y_{-p}}$ induces the same set of conditional variance functions as a pure AR model (Supplementary Material Section C). From this, we can repeat the probabilistic arguments of Section (ref) to form an $R^2$ prior on ARX models (which we will refer to as the ARR2 prior for ARX models).
A more complete derivation of this prior is provided in Supplementary Material Section E. This prior makes explicit that both the $\operatorname{AR}$ and exogenous regression components contribute to the total variance of the model $\tau^2$, which is then decomposed probabilistically via $\psi$. In addition to the $p$ components for each lag, $\psi$ now includes $m$ additional components attributed to each exogenous covariate. These two sets of components compete for the relative contribution to the total variance, and the hyperparameters $\xi$ can be used to nudge the model toward particular decompositions.
The hyperparameters $\left(\xi_1,\dotsc,\xi_p\right)$ may be set according to the proposed time-series lag structure, as in Figure (ref), and $\left(\xi_{p+1},\dotsc,\xi_m\right)$ according to prior knowledge of sparsity and correlation structure. For example, one may set $(\xi_1,\dotsc,\xi_p)$ according to Minnesota decay as for the AR models, while imposing sparsity-favouring decomponsition for $x_t$ by setting $(\xi_{p+1},\dotsc, \xi_{p+m})$ to a relatively low value of 0.1.
Similar to the derived $R^2$ prior for observable time-series dynamics, we present in this section an extension to the state-space framework which is a flexible family of models for joint estimation of observable and latent time-series. These have a long tradition in Bayesian estimation kitagawa1984smoothness,west2006bayesian and probabilistic filtering sarkka2023bayesian and are routinely used in the social and hard sciences chan2023bayesian.
We consider state-space models of the following form:
where $y_t \in \mathbb{R}$ is the scalar valued target and $x_t \in \mathbb{R}^{m\times 1}$ exogenous covariates. $s_t \in \mathbb{R}^{Q\times 1}$ are unknown states, $G$ is a state coefficient vector, $\Phi \in \mathbb{R}^{Q\times Q}$ is some invertible state transition matrix, and $\Sigma_s$ the state error covariance matrix with diagonal elements $(\sigma_{s_1} ,\dots,\sigma_{s_Q})$. Assume, further, that there are $Q$ initial conditions to the states, $s_0 \sim N(0,\Sigma_s)$. The equation for $y_t$ is called the observation equation in the state-space literature and that for $s_t$ the transition equation harvey_forecasting_1990. The generality of the state-space model in Equations (ref)-(ref) can be seen from the fact that AR, MA, ARMA, ARX, Bayesian structural time-series models brodersen2015inferring,scott2014predicting, and many more, are special cases harvey_forecasting_1990.\footnote{While Equation (ref) is formulated as an $\operatorname{AR}(1)$ process, any scalar or vector valued auto-regressive time-series of order $p$ can be written as an $\operatorname{AR}(1)$ process via its companion form hamilton2020time.}
Define the predictor term of the observation equation as $\mu_t = x_t^{\prime}\beta + s_t^{\prime}G$. Under the assumption of independence between $x_t$ and states, the variance of the observation equation can be factored as $\sigma^2_{\mu_t} = \var{x_t^{\prime} \beta} + \var{s_t^{\prime}G}$. As before, the variance contribution of the exogenous covariate component simplifies to $\var{x_t^{\prime} \beta} = \sum_{j=1}^m\lambda_j^2\sigma_{x_j}^2$. The marginal variance contribution of the states is complicated by the law of motion in Equation (ref), however, we make three simplifying assumptions for ease of exposition of the following defintion:
where $\boldsymbol{I_{Q}}$ denotes the $Q$-dimensional identity matrix. While assumptions (ref) and (ref) may seem limiting for the prior on a state-space model's $R^2$, a large class of models such as dynamic regressions chan2020reducing, unobserved component models in economics fruhwirth2010stochastic and forecasting models kohns_nowcasting_2022,potjagailo2023flexible similarly assume independent $\operatorname{AR}(1)$ processes for the states. In Supplementary Material Section F, we derive the $R^2$ relaxing some of these assumptions. Assumption (ref) is needed in order for the $R^2$ to be well defined, as $R^2 \rightarrow 1$ with non-stationarity of the states.
A more complete derivation of this prior is provided in Supplementary Material Section F. The hierarchy of the ARR2 prior makes three types of dependencies explicit. First, unlike previous priors for state-space models, the states for all $t$ depend on the noise variance, $\sigma^2$. This is needed in order for the relative shares of explained variance between the regression and state component to stay constant to changes in the scale of the observation noise. One may think of this as a state component analog to the prior scaling of $\beta$ by $\sigma^2$ in Definition (ref). Second, a state's prior variance decreases with its AR parameter $\phi_q$. This makes intuitive sense, since not much can be learned from states which are random noise with low serial correlation. Lastly, the ARR2 prior makes joint shrinkage of the state and regression component explicit via decomposition of the total variance $\tau^2$. Previous approaches on the other hand, assume independence between the state and other model components in the observation equation cadonna2020triple. This neglects the fact that states and other components in the observation equation compete to explain shares in the total variance of $y_t$.
The dimensionality of the estimation problem in Definition (ref) can easily get very high. Consider the addition of $p$ lags per covariate. In these situations, one can reduce computational complexity by decomposing the prior variance at the group-level: $\beta_{i,j} \sim \normal(0,\tilde{\psi}_{i,j}\tau^2) $, where $\tilde{\psi}_{i,j} = w_j \psi_i $, $\psi \sim \Dirichlet(\xi_1,\dots,\xi_m)$ and $\sum_{j=1}^p w_j = 1$ are some deterministic weights for $j \in \{1,\dots,p\}$ and $i \in \{1,\dots,m\}$. In its simplest form, set $w_j=1/p \forall j$ such that the simplex dimensionality reduces to $m+Q$ instead of $mp + Q$. Inspired by a Minnesota-like decomposition for a group of lags, one may then set $w_j$ to $(1/j^2) / (\sum_{s =1 }^p 1/s^2)$. \paragraph{Local-linear trend model with covariates} To illustrate the above, consider the following simplified Bayesian structural equation model, referred to as the local-trend model (LTX):
where $\delta_t$ is an unobserved time-trend and $\phi \in (-1,1)$. It can be shown that Equation $\ref{eq:llt_def}$ implies an ARMA-X model for $y_t$ harvey_forecasting_1990. The $R^2$ for this prior is then defined as:
Hence, $R^2 \sim \betadist(\mu_{R^2},\varphi_{R^2})$ again implies $\tau^2 \sim \beta^\prime(\mu_{R^2},\varphi_{R^2})$. Then
In Supplementary Material Section F.1, we show how the ARR2 prior can be applied to the popular dynamic regression model.
In this section, we primarily aim to understand the quality of parameter estimation and predictions when estimating models of increasing complexity for different data-generating processes (DGPs) simulated from AR, ARX and LTX models. We investigate the behaviour when the complexity of the DGP stays fixed as the complexity of the estimated model grows. This is a common problem in the Bayesian workflow gelman2020bayesian in which model building will often involve sequentially adding more complexity. For all simulations, we fit the proposed ARR2 prior with two hyperparameter settings, Minnesota ARR2 and flat ARR2 (see Section (ref)). We compare these models to two types of commonly used prior for time-series analysis, the Minnesota giannone_prior_2012 and regularised horseshoe piironen_hyperprior_2017, as well as a relatively wide non-hierarchical standard Gaussian prior. For all simulation experiments, we simulate 25 sets of data from each DGP.
We measure goodness of parameter recovery by average root-mean-squared error ($\mathrm{RMSE}$) between the posterior mean and the true coefficients. Let $\hat{\theta} = 1/S \sum^S_{s = 1 }\theta^{(s)}$, where $S$ are the number of retained posterior draws. We define $\mathrm{RMSE}$ as:
where $K$ is the dimensionality of $\theta$ and $||\theta||^2_2$ calculates the squared Euclidean norm.
Predictions are evaluated using the leave-future-out (LFO) expected log predictive density vehtari_practical_2017,burkner_approximate_2020, for $M$--step-ahead future observations, which we compute as
which is estimated by
In all simulations below, we focus for simplicity on 1-step-ahead LFO predictions $(M=1)$. For a set of integers $A\subset\mathbb{Z}$, we denote $y_{-A} = \{y_j: j = 1,\dotsc,n,\,j\notin A\}$.
In order to compare predictions between simulation exercises with differing time-series lengths, we report mean log predictive density (MLPD) which is \(\mathrm{elpd}_\mathrm{LFO} / (T-L)\). Predictions are always evaluated on a hold-out set where $L$ is equal to half the number of observations ($T/2$).
All experiments were performed using \proglang{Stan} stan using the cmdstanr interface in \proglang{R}, and the source code is freely available at: \url{https://github.com/n-kall/arr2}. Supplementary Material Section C has \proglang{Stan} stan code for the ARR2 prior for AR models.
We generate data from three different AR processes whose AR coefficients mimic 3 types of processes typically found in stationary time-series models: Minnesota, delayed relevance, and dampened oscillations. Minnesota type and hump shape (delayed relevance) AR parameter profiles are commonly found in economic time-series doan1984forecasting which we adapt from the article by mogliani2021bayesian. These DGPs favour priors such as the Minnesota and RHS prior. Dampened oscillations, on the other hand, are more common to physical applications west2006bayesian and are included as an example of complicated auto-covariances with many large coefficients. We expect this DGP to disadvantage priors which assume exchangeability in the AR parameters. The parameters are defined in Table (ref) and are chosen to induce true $R^2$ of around 0.7. For all DGPs, we set $\sigma^2 = 1$. The true lag order is always fixed to eight. Under each DGP, we simulate $T = 120$ observations, and fit AR models of increasing order (up to $p = T / 2 = 60$).
The results for the parameter recovery are shown in Figure (ref). As expected, priors whose shrinkage across lag polynomials matches the DGP's coefficient profile perform best. The Minnesota ARR2 and Minnesota prior perform best for the DGPs with decaying AR coefficients, while the flat ARR2 and RHS perform best for the delayed DGP.
Interestingly, all priors bar the Gaussian exhibit estimation error which is nearly independent of the lag order. Independent Gaussians, result in much higher RMSE than the any of the other priors. This is particularly noticeable as the model size increases (more irrelevant lags are included).