[Paper Review] varstan: An R package for Bayesian analysis of structured time series models with Stan

Izhar Asael Alonzo Matamoros, Cristian Andrés Cruz Torres|arXiv (Cornell University)|May 20, 2020

Statistical Methods and Inference19 references4 citations

TL;DR

This paper introduces varstan, an R package that enables Bayesian analysis of structured time series models using Stan’s efficient Hamiltonian Monte Carlo sampling. It supports a wide range of models—including SARIMA, GARCH, dynamic regression, and VARMA—offers flexible prior specification, and provides advanced model comparison via WAIC, LOOIC, and Bayes factors, demonstrating superior performance on U.S. monthly birth data with a selected SARIMA(1,1,1)(1,1,1)[12] model.

ABSTRACT

varstan is an \proglang{R} package for Bayesian analysis of time series models using \proglang{Stan}. The package offers a dynamic way to choose a model, define priors in a wide range of distributions, check model's fit, and forecast with the m-steps ahead predictive distribution. The users can widely choose between implemented models such as extit{multiplicative seasonal ARIMA, dynamic regression, random walks, GARCH, dynamic harmonic regressions,VARMA, stochastic Volatility Models, and generalized t-student with unknown degree freedom GARCH models}. Every model constructor in \pkg{varstan} defines weakly informative priors, but prior specifications can be changed in a dynamic and flexible way, so the prior distributions reflect the parameter's initial beliefs. For model selection, the package offers the classical information criteria: AIC, AICc, BIC, DIC, Bayes factor. And more recent criteria such as Widely-applicable information criteria ( extit{WAIC}), and the Bayesian leave one out cross-validation ( extit{loo}). In addition, a Bayesian version for automatic order selection in seasonal ARIMA and dynamic regression models can be used as an initial step for the time series analysis.

Motivation & Objective

• To develop a flexible, user-friendly R package for Bayesian time series analysis using Stan’s efficient sampling algorithms.
• To support a broad range of structured time series models, including SARIMA, GARCH, dynamic regression, and VARMA.
• To enable dynamic and flexible prior specification that reflects user beliefs while maintaining weakly informative defaults.
• To provide comprehensive model selection tools, including WAIC, LOOIC, DIC, AICc, BIC, and Bayes factors.
• To offer a practical, accessible interface inspired by classic R time series packages like forecast and astsa, lowering the barrier to Bayesian time series modeling.

Proposed method

• The package uses Stan’s No-U-Turn Sampler (NUTS) for efficient Bayesian inference via Hamiltonian Monte Carlo, enabling fast convergence in high-dimensional models.
• It implements a modular interface where each model type (e.g., SARIMA, GARCH) is constructed via a dedicated constructor function with default weakly informative priors.
• Users can customize priors dynamically, allowing for subjective belief incorporation while maintaining computational stability.
• Model fitting is performed via MCMC sampling using Stan’s R interface, with diagnostics such as R-hat and effective sample size (ESS) reported for convergence assessment.
• Model comparison is supported through multiple criteria: WAIC, LOOIC, DIC, AICc, BIC, and Bayes factors computed via marginal likelihoods.
• The package includes automated order selection for SARIMA and dynamic regression models, serving as a starting point for model exploration.

Experimental results

Research questions

• RQ1 How can Bayesian time series modeling be made more accessible and flexible for practitioners using R?
• RQ2 Can Stan’s NUTS sampler effectively handle complex structured time series models like SARIMA and GARCH with non-conjugate priors?
• RQ3 How do model comparison criteria like WAIC and LOOIC perform in selecting the best-fitting time series model compared to classical criteria?
• RQ4 To what extent does the automated order selection function in varstan improve model identification compared to manual specification?
• RQ5 How does the performance of a manually selected SARIMA model compare to an automatically selected one using the auto.sarima function in terms of predictive accuracy and model fit?

Key findings

• The SARIMA(1,1,1)(1,1,1)[12] model (model1) was selected as the best-fitting model for U.S. monthly birth data based on a log Bayes factor of 199.14, providing strong evidence in its favor.
• The model comparison via LOOIC showed that the auto.sarima-selected model (SARIMA(0,1,2)(1,1,1)[12]) had a lower LOOIC (2470) than model1 (2472), indicating better predictive performance.
• The LOOIC difference was 2 units, with a standard error of 5.8, suggesting that the auto.sarima model is slightly preferred, though the difference is small.
• The package successfully estimated posterior distributions for all parameters with R-hat values close to 1.0, indicating good convergence across chains.
• The harmonic regression model, while fitting more parameters, was outperformed by the seasonal ARIMA model in model selection criteria, confirming the latter’s better fit for the sinusoidal birth pattern.
• The package demonstrated robust performance in handling complex models, including GARCH and stochastic volatility variants, with flexible prior specification and efficient sampling.