Chapter 14 14 Bayesian Structural Time Series Models
14.1 Introduction
Time series data pervade scientific and applied domains, from economics and epidemiology to environmental monitoring and finance. Extracting meaningful patterns from such data requires statistical models that can handle temporal dependencies, seasonality, trends, and structural breaks while providing interpretable decompositions. Bayesian Structural Time Series (BSTS) models offer a principled framework for addressing these challenges by decomposing time series into latent components while incorporating prior information and quantifying uncertainty.
Initially popularized by Google for causal impact analysis, BSTS models extend classical structural time series approaches by embedding them within a fully Bayesian inference framework. However, many real-world time series exhibit periods of varying volatility—economic indicators may experience stable growth followed by high volatility due to macroeconomic shocks, while financial markets alternate between calm and turbulent periods. Standard BSTS models assume constant volatility in trend components, which can be overly restrictive and obscure critical temporal dynamics.
In this comprehensive post, we explore both canonical BSTS models and their extension to time-varying trend volatility. We provide complete implementations in Stan, demonstrate fitting procedures in R using RStan, and illustrate applications across different scenarios. This unified treatment enables practitioners to understand the full spectrum of BSTS modeling capabilities and choose appropriate specifications for their specific time series challenges.
14.1.0.1 Model Specification: Core BSTS Framework
The canonical BSTS model represents an observed time series as the sum of interpretable latent components:
\[y_t = \mu_t + \tau_t + \gamma_t + \varepsilon_t, \quad \varepsilon_t \sim \mathcal{N}(0, \sigma^2)\]
where: - \(\mu_t\) is the local level component - \(\tau_t\) represents the stochastic trend (slope or drift) - \(\gamma_t\) captures seasonal effects - \(\varepsilon_t\) is the observation error term
The local level and trend components follow a local linear trend specification:
\[\mu_t = \mu_{t-1} + \tau_{t-1} + \eta_t, \quad \eta_t \sim \mathcal{N}(0, \sigma_\mu^2)\]
\[\tau_t = \tau_{t-1} + \zeta_t, \quad \zeta_t \sim \mathcal{N}(0, \sigma_\tau^2)\]
This formulation allows the trend to evolve smoothly over time while maintaining flexibility in both level and slope dynamics. The seasonal component \(\gamma_t\) is typically modeled using periodic structures that ensure the seasonal effects sum to zero over each complete cycle.
14.1.0.2 Extension to Time-Varying Trend Volatility
In many applications, the assumption of constant trend volatility (\(\sigma_\tau^2\)) proves inadequate. Economic time series, for instance, may exhibit stable growth periods punctuated by high-volatility episodes during recessions or policy changes. To accommodate such dynamics, we extend the standard framework by allowing the volatility of the trend innovation to evolve over time.
We introduce a stochastic volatility process for the trend innovation variance:
\[\log \sigma_{\tau, t}^2 = h_t, \quad h_t = h_{t-1} + \nu_t, \quad \nu_t \sim \mathcal{N}(0, \sigma_h^2)\]
\[\zeta_t \sim \mathcal{N}(0, \exp(h_t))\]
This specification implements a log-random walk for the trend volatility, ensuring positive variance while allowing for gradual or abrupt changes in trend smoothness. The parameter \(\sigma_h^2\) controls how quickly the trend volatility can change: larger values permit more rapid volatility shifts, while smaller values enforce smoother volatility evolution.
14.1.0.3 Data Simulation: Building Synthetic Time Series
To demonstrate both modeling approaches, we begin by simulating realistic time series data that incorporates the key components of interest.
- Basic BSTS Simulation
set.seed(42)
n <- 120 # number of time points
m <- 12 # seasonal period
level <- numeric(n)
trend <- numeric(n)
season <- numeric(n)
y <- numeric(n)
# Initialize components
level[1] <- 10
trend[1] <- 0.5
seasonal_pattern <- sin(2 * pi * (1:m) / m)
season[1:m] <- seasonal_pattern
# Generate time series
for (t in 2:n) {
level[t] <- level[t-1] + trend[t-1] + rnorm(1, 0, 0.3)
trend[t] <- trend[t-1] + rnorm(1, 0, 0.05)
season[t] <- season[t %% m + 1] # repeat seasonal pattern
y[t] <- level[t] + season[t] + rnorm(1, 0, 1)
}
ts.plot(y, main = "Simulated BSTS Time Series")- Time-Varying Volatility Simulation
For the extended model with stochastic volatility, we modify the simulation to include time-varying trend volatility:
set.seed(123)
n <- 120
m <- 12
level <- numeric(n)
trend <- numeric(n)
h <- numeric(n)
sigma_trend_t <- numeric(n)
season <- numeric(n)
y <- numeric(n)
# Initialize components
level[1] <- 10
trend[1] <- 0.5
h[1] <- log(0.05^2) # initial log-variance
seasonal_pattern <- sin(2 * pi * (1:m) / m)
season[1:m] <- seasonal_pattern
# Generate time series with time-varying trend volatility
for (t in 2:n) {
# Update log-volatility
h[t] <- h[t - 1] + rnorm(1, 0, 0.1)
sigma_trend_t[t] <- sqrt(exp(h[t]))
# Update components
trend[t] <- trend[t - 1] + rnorm(1, 0, sigma_trend_t[t])
level[t] <- level[t - 1] + trend[t - 1] + rnorm(1, 0, 0.3)
season[t] <- season[t %% m + 1]
y[t] <- level[t] + season[t] + rnorm(1, 0, 1)
}
ts.plot(y, main = "Simulated Time Series with Time-Varying Trend Volatility")This simulation generates periods of varying trend smoothness, providing a realistic testbed for the extended model. ## Stan Implementation
14.1.1 Stan Implementation: Basic BSTS Model
We implement the core BSTS model in Stan, emphasizing clarity and computational efficiency:
functions {
vector seasonal_component(vector s, int N, int m) {
vector[N] seasonal;
for (t in 1:N) {
if (t > m)
seasonal[t] = s[t - m];
else
seasonal[t] = s[t];
}
return seasonal;
}
}
data {
int<lower=1> N;
int<lower=1> m; // seasonal period
vector[N] y;
}
parameters {
real<lower=0> sigma_obs;
real<lower=0> sigma_level;
real<lower=0> sigma_trend;
real<lower=0> sigma_seasonal;
vector[N] level;
vector[N] trend;
vector[N] seasonal;
}
model {
// State evolution equations
level[2:N] ~ normal(level[1:(N-1)] + trend[1:(N-1)], sigma_level);
trend[2:N] ~ normal(trend[1:(N-1)], sigma_trend);
seasonal[(m+1):N] ~ normal(seasonal[1:(N-m)], sigma_seasonal);
// Observation equation
y ~ normal(level + seasonal_component(seasonal, N, m), sigma_obs);
// Priors
sigma_obs ~ normal(0, 1);
sigma_level ~ normal(0, 0.5);
sigma_trend ~ normal(0, 0.1);
sigma_seasonal ~ normal(0, 0.5);
}14.1.2 Stan Implementation: Time-Varying Trend Volatility
The extended model incorporates stochastic volatility into the trend component:
functions {
vector seasonal_component(vector s, int N, int m) {
vector[N] seasonal;
for (t in 1:N) {
if (t > m)
seasonal[t] = s[t - m];
else
seasonal[t] = s[t];
}
return seasonal;
}
}
data {
int<lower=1> N;
int<lower=1> m;
vector[N] y;
}
parameters {
real<lower=0> sigma_obs;
real<lower=0> sigma_level;
real<lower=0> sigma_h;
vector[N] level;
vector[N] trend;
vector[N] seasonal;
vector[N] h; // log variance of trend innovation
}
model {
// Priors
sigma_obs ~ normal(0, 1);
sigma_level ~ normal(0, 0.5);
sigma_h ~ normal(0, 0.2);
h[1] ~ normal(log(0.05^2), 0.1);
// Seasonal component evolution
seasonal[(m+1):N] ~ normal(seasonal[1:(N-m)], 0.5);
// Latent state processes
for (t in 2:N) {
h[t] ~ normal(h[t-1], sigma_h);
trend[t] ~ normal(trend[t-1], exp(0.5 * h[t]));
level[t] ~ normal(level[t-1] + trend[t-1], sigma_level);
}
// Observation equation
y ~ normal(level + seasonal_component(seasonal, N, m), sigma_obs);
}The key innovation lies in the stochastic evolution of h[t], which governs the time-varying volatility of trend innovations through the exponential transformation exp(0.5 * h[t]).
- Model Fitting in R with RStan
Both models can be fitted using RStan with similar workflows, differing primarily in their complexity and computational requirements.
- Basic BSTS Model Fitting
library(rstan)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
# Prepare data
data_list_basic <- list(N = length(y), m = m, y = y)
# Fit basic model
fit_basic <- stan(model_code = stan_model_bsts_basic,
data = data_list_basic,
chains = 4, iter = 2000, warmup = 1000, seed = 42)
print(fit_basic, pars = c("sigma_obs", "sigma_level", "sigma_trend", "sigma_seasonal"))- Time-Varying Volatility Model Fitting
# Fit extended model
fit_extended <- stan(model_code = stan_model_bsts_extended,
data = data_list_basic,
chains = 4, iter = 2000, warmup = 1000, seed = 123)
print(fit_extended, pars = c("sigma_obs", "sigma_level", "sigma_h"))The posterior for sigma_h provides insights into volatility dynamics: larger values indicate more variable trend volatility, while smaller values suggest smoother volatility evolution.
- Visualization and Interpretation
Effective visualization of BSTS results requires displaying both point estimates and uncertainty quantification for the latent components.
- Basic Model Visualization
# Extract posterior samples
post_basic <- extract(fit_basic)
# Compute posterior summaries
level_mean <- apply(post_basic$level, 2, mean)
level_lower <- apply(post_basic$level, 2, quantile, probs = 0.025)
level_upper <- apply(post_basic$level, 2, quantile, probs = 0.975)
trend_mean <- apply(post_basic$trend, 2, mean)
library(ggplot2)
df_basic <- data.frame(
time = 1:n,
observed = y,
level = level_mean,
level_lower = level_lower,
level_upper = level_upper,
trend = trend_mean
)
# Plot level component with uncertainty
ggplot(df_basic, aes(x = time)) +
geom_line(aes(y = observed), color = "black", linetype = "dashed", alpha = 0.7) +
geom_line(aes(y = level), color = "blue", size = 1) +
geom_ribbon(aes(ymin = level_lower, ymax = level_upper),
fill = "blue", alpha = 0.2) +
labs(title = "BSTS Level Component with 95% Credible Intervals",
x = "Time", y = "Value") +
theme_minimal()
# Plot trend evolution
ggplot(df_basic, aes(x = time)) +
geom_line(aes(y = trend), color = "red", size = 1) +
labs(title = "Estimated Trend Component",
x = "Time", y = "Trend") +
theme_minimal()- Time-Varying Volatility Visualization
# Extract posterior samples for extended model
post_extended <- extract(fit_extended)
# Compute posterior summaries
h_mean <- apply(post_extended$h, 2, mean)
h_lower <- apply(post_extended$h, 2, quantile, probs = 0.025)
h_upper <- apply(post_extended$h, 2, quantile, probs = 0.975)
trend_mean_ext <- apply(post_extended$trend, 2, mean)
df_extended <- data.frame(
time = 1:n,
log_vol = h_mean,
log_vol_lower = h_lower,
log_vol_upper = h_upper,
trend = trend_mean_ext
)
# Plot time-varying log-volatility
ggplot(df_extended, aes(x = time)) +
geom_line(aes(y = log_vol), color = "red", size = 1) +
geom_ribbon(aes(ymin = log_vol_lower, ymax = log_vol_upper),
fill = "red", alpha = 0.2) +
labs(title = "Time-Varying Trend Log-Volatility with 95% Credible Intervals",
y = "log(σ²_τ)", x = "Time") +
theme_minimal()
# Plot corresponding trend with varying smoothness
ggplot(df_extended, aes(x = time)) +
geom_line(aes(y = trend), color = "blue", size = 1) +
labs(title = "Trend Component with Time-Varying Volatility",
y = "Trend", x = "Time") +
theme_minimal()These visualizations reveal how uncertainty in trend evolution changes over time, enabling analysts to identify periods of instability or structural transitions.
14.1.2.1 Applications and Extensions
BSTS models find applications across diverse domains where interpretable time series decomposition is valuable. The basic framework handles missing data naturally and provides a foundation for causal impact analysis through counterfactual prediction. The time-varying volatility extension proves particularly valuable in financial econometrics, where volatility clustering and regime changes are common phenomena. In financial time series, periods of market calm alternate with high-volatility episodes during crises or major economic announcements. The time-varying trend volatility model can capture these dynamics while maintaining the interpretable decomposition that BSTS models provide. Economic indicators often exhibit structural breaks and changing volatility regimes due to policy interventions, technological changes, or external shocks. The flexible volatility specification enables more adaptive forecasting in such environments. Climate time series frequently display non-stationary variance patterns due to changing atmospheric or oceanic conditions. The stochastic volatility framework can accommodate such dynamics while preserving seasonal and trend interpretability.
Both modeling frameworks can be extended to hierarchical structures where multiple related time series share common components, or to multivariate settings where cross-series dependencies are of interest. The Stan implementations provide a flexible foundation for such extensions. BSTS models naturally accommodate regression components with spike-and-slab priors for variable selection. This capability enables automatic identification of relevant predictors while maintaining the structural decomposition that makes BSTS models interpretable. The computational complexity of BSTS models scales with the number of time points and latent states. The basic model typically converges within standard MCMC runs, while the time-varying volatility extension requires longer chains due to the additional stochastic volatility component.
Stan’s efficient gradient-based sampling (HMC/NUTS) handles the high-dimensional latent state spaces more effectively than traditional Gibbs sampling approaches. However, practitioners should monitor convergence diagnostics carefully, particularly for the volatility parameters in the extended model.
For very long time series, computational efficiency can be improved through state space formulations that marginalize over latent states, though this sacrifices direct access to the latent component estimates that make BSTS models interpretable.
14.2 Model Comparison and Selection
Comparing basic and extended BSTS models requires balancing goodness-of-fit against complexity. The Widely Applicable Information Criterion (WAIC) provides a practical tool for model comparison within the Bayesian framework:
library(loo)
# Compute WAIC for model comparison
waic_basic <- waic(extract_log_lik(fit_basic))
waic_extended <- waic(extract_log_lik(fit_extended))
# Compare models
loo_compare(waic_basic, waic_extended)Additionally, posterior predictive checks help assess whether the added complexity of time-varying volatility improves model adequacy for the specific application.
14.3 Conclusion
Bayesian Structural Time Series models provide a powerful and interpretable framework for analyzing complex temporal data. The basic BSTS specification offers transparent decomposition into level, trend, and seasonal components while maintaining full uncertainty quantification through the Bayesian approach. The extension to time-varying trend volatility addresses limitations of constant volatility assumptions, enabling more adaptive modeling of time series with changing uncertainty regimes.
The Stan implementations presented here provide complete, working examples that practitioners can adapt to their specific applications. The state-space framework naturally handles missing data and irregular observation patterns, while the Bayesian approach enables principled incorporation of prior information and coherent uncertainty propagation.
As computational tools continue to improve and the demand for interpretable time series models grows, BSTS models are positioned to play an increasingly important role across scientific and applied domains. The flexibility demonstrated through the time-varying volatility extension illustrates how the basic framework can be adapted to address specific modeling challenges while preserving the interpretability that makes BSTS models valuable for practical time series analysis.
The unified treatment presented here equips practitioners with both foundational understanding and practical tools for implementing BSTS models across a spectrum of applications, from basic seasonal decomposition to advanced volatility modeling in dynamic environments.