Bayesian Structural Time Series Models in Stan
In the previous post, we fitted exponential smoothing models in Stan using a conditional likelihood, where the level, trend, and seasonal states were deterministic functions of the parameters and the data. Here we take the related but more general state-space route: the level, slope, and seasonal effects become stochastic latent states with their own innovation variances, and the posterior is computed jointly over states and parameters.
Bayesian structural time series (BSTS) models decompose a series into interpretable components and quantify uncertainty in each. The approach became widely known through the work of Scott and Varian at Google, notably in causal-impact analysis. The models we build below are simpler than that full framework (there is no regression component or spike-and-slab selection), but they share its core idea.
One limitation of the standard formulation is that the innovation variances are constant. Economic and financial series often alternate between calm and turbulent periods, and a fixed slope variance may then be too rigid. We therefore fit two models: a basic BSTS with a local linear trend and a seasonal component, and an extension in which the variance of the slope innovations follows a stochastic-volatility process.
Model Specification
The observed series is the sum of a level, a seasonal effect, and observation noise:
\[y_t = \mu_t + \gamma_t + \varepsilon_t, \qquad \varepsilon_t \sim \mathcal{N}(0, \sigma_\varepsilon^2).\]The level follows a local linear trend, in which the slope $\tau_t$ enters through the level equation rather than directly into $y_t$:
\[\mu_t = \mu_{t-1} + \tau_{t-1} + \eta_t, \qquad \eta_t \sim \mathcal{N}(0, \sigma_\mu^2),\] \[\tau_t = \tau_{t-1} + \zeta_t, \qquad \zeta_t \sim \mathcal{N}(0, \sigma_\tau^2).\]For a seasonal period $m$, we use the dummy-variable form, in which the seasonal effects over any $m$ consecutive periods sum to a small random innovation:
\[\gamma_t = -\sum_{j=1}^{m-1} \gamma_{t-j} + \omega_t, \qquad \omega_t \sim \mathcal{N}(0, \sigma_\gamma^2).\]When $\sigma_\gamma = 0$ the pattern repeats exactly and sums to zero over each cycle. Positive values allow it to drift slowly. This constraint also separates the seasonal pattern from the level. Without it, a constant could be moved between the two without changing the fit.
Time-Varying Trend Volatility
To let the smoothness of the slope change over time, we replace the constant $\sigma_\tau^2$ with a time-varying variance whose logarithm follows a random walk:
\[\zeta_t \sim \mathcal{N}\!\left(0, e^{h_t}\right), \qquad h_t = h_{t-1} + \nu_t, \qquad \nu_t \sim \mathcal{N}(0, \sigma_h^2).\]The exponential guarantees a positive variance. The parameter $\sigma_h$ controls how fast the log-variance can move, with larger values allowing quicker shifts between stable and unstable slope regimes.
This extension is demanding for the data. The slope is itself a latent state, and its innovation variance is a latent state of that state, so a series of about 100 observations carries little direct information about $h_t$. We should expect wide credible intervals and some sensitivity to the prior on $\sigma_h$.
Simulating the Data
The first simulation follows the basic model. The seasonal effects use the same sum-to-zero recursion as the Stan program, with a small innovation.
set.seed(42)
n <- 120
m <- 12
level <- numeric(n)
trend <- numeric(n)
season <- numeric(n)
y <- numeric(n)
level <- 10
trend <- 0.5
season[1:m] <- sin(2 * pi * (1:m) / m) # sums to zero over one cycle
y <- level + season + rnorm(1, 0, 1)
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)
if (t > m) {
season[t] <- -sum(season[(t - m + 1):(t - 1)]) + rnorm(1, 0, 0.05)
}
y[t] <- level[t] + season[t] + rnorm(1, 0, 1)
}
y_basic <- y
truth_basic <- data.frame(level = level, trend = trend, season = season)
ts.plot(y_basic, main = "Simulated BSTS Series")
The second simulation lets the trend innovation variance follow a log random walk. We use $\sigma_h = 0.2$, which is large enough to produce visible changes in smoothness over 120 periods. A smaller value would give a trend that is almost indistinguishable from the constant-variance case.
set.seed(123)
n <- 120
m <- 12
sigma_h_true <- 0.2
level <- numeric(n)
trend <- numeric(n)
season <- numeric(n)
h <- numeric(n)
y <- numeric(n)
level <- 10
trend <- 0.5
h <- log(0.05^2)
season[1:m] <- sin(2 * pi * (1:m) / m)
y <- level + season + rnorm(1, 0, 1)
for (t in 2:n) {
h[t] <- h[t - 1] + rnorm(1, 0, sigma_h_true)
trend[t] <- trend[t - 1] + rnorm(1, 0, exp(0.5 * h[t]))
level[t] <- level[t - 1] + trend[t - 1] + rnorm(1, 0, 0.3)
if (t > m) {
season[t] <- -sum(season[(t - m + 1):(t - 1)]) + rnorm(1, 0, 0.05)
}
y[t] <- level[t] + season[t] + rnorm(1, 0, 1)
}
y_sv <- y
truth_sv <- data.frame(level = level, trend = trend, season = season, h = h)
ts.plot(y_sv, main = "Simulated Series with Time-Varying Trend Volatility")
Stan Implementation: Basic BSTS
Two implementation choices matter here. First, the program standardises $y$ internally, so the priors are expressed in standard-deviation units regardless of the scale of the data. Posterior quantities are transformed back to the original scale in generated quantities.
Second, the latent states use a non-centred parameterisation: Stan samples standard-normal innovations, and the states are built from them. Centred random walks with unknown innovation scales tend to produce funnel-shaped posteriors, particularly when a scale is close to zero, and non-centring usually reduces divergences.
Save this as bsts_basic.stan.
data {
int<lower=2> N;
int<lower=2, upper=N> m;
vector[N] y;
}
transformed data {
real y_mean = mean(y);
real y_sd = sd(y);
vector[N] y_std = (y - y_mean) / y_sd;
}
parameters {
real<lower=0> sigma_obs;
real<lower=0> sigma_level;
real<lower=0> sigma_trend;
real<lower=0> sigma_seasonal;
real level_1;
real trend_1;
vector[m - 1] season_init;
vector[N - 1] z_level;
vector[N - 1] z_trend;
vector[N - m + 1] z_season;
}
transformed parameters {
vector[N] level;
vector[N] trend;
vector[N] seasonal;
level = level_1;
trend = trend_1;
seasonal[1:(m - 1)] = season_init;
for (t in 2:N) {
trend[t] = trend[t - 1] + sigma_trend * z_trend[t - 1];
level[t] = level[t - 1] + trend[t - 1] + sigma_level * z_level[t - 1];
}
for (t in m:N) {
seasonal[t] = -sum(seasonal[(t - m + 1):(t - 1)])
+ sigma_seasonal * z_season[t - m + 1];
}
}
model {
// Scale priors (standardised units)
sigma_obs ~ normal(0, 1);
sigma_level ~ normal(0, 0.1);
sigma_trend ~ normal(0, 0.05);
sigma_seasonal ~ normal(0, 0.1);
// Initial states
level_1 ~ normal(y_std, 1);
trend_1 ~ normal(0, 0.1);
season_init ~ normal(0, 0.5);
// Non-centred innovations
z_level ~ std_normal();
z_trend ~ std_normal();
z_season ~ std_normal();
y_std ~ normal(level + seasonal, sigma_obs);
}
generated quantities {
vector[N] level_orig = y_mean + y_sd * level;
vector[N] trend_orig = y_sd * trend;
vector[N] seasonal_orig = y_sd * seasonal;
real sigma_obs_orig = y_sd * sigma_obs;
vector[N] y_rep;
for (t in 1:N) {
y_rep[t] = normal_rng(level_orig[t] + seasonal_orig[t], sigma_obs_orig);
}
}
The scale priors are weakly informative on the standardised scale, and the choices are judgement calls. For a series with a very different signal-to-noise ratio, they should be checked with prior predictive simulations.
Stan Implementation: Time-Varying Trend Volatility
The extended model replaces sigma_trend with a log-variance path h. The seasonal scale is now a parameter rather than a hard-coded constant. Save this as bsts_sv.stan.
data {
int<lower=2> N;
int<lower=2, upper=N> m;
vector[N] y;
}
transformed data {
real y_mean = mean(y);
real y_sd = sd(y);
vector[N] y_std = (y - y_mean) / y_sd;
}
parameters {
real<lower=0> sigma_obs;
real<lower=0> sigma_level;
real<lower=0> sigma_seasonal;
real<lower=0> sigma_h;
real level_1;
real trend_1;
real h_1;
vector[m - 1] season_init;
vector[N - 1] z_level;
vector[N - 1] z_trend;
vector[N - 1] z_h;
vector[N - m + 1] z_season;
}
transformed parameters {
vector[N] level;
vector[N] trend;
vector[N] seasonal;
vector[N] h;
level[1] = level_1;
trend[1] = trend_1;
h[1] = h_1;
seasonal[1:(m - 1)] = season_init;
for (t in 2:N) {
h[t] = h[t - 1] + sigma_h * z_h[t - 1];
trend[t] = trend[t - 1] + exp(0.5 * h[t]) * z_trend[t - 1];
level[t] = level[t - 1] + trend[t - 1] + sigma_level * z_level[t - 1];
}
for (t in m:N) {
seasonal[t] = -sum(seasonal[(t - m + 1):(t - 1)])
+ sigma_seasonal * z_season[t - m + 1];
}
}
model {
sigma_obs ~ normal(0, 1);
sigma_level ~ normal(0, 0.1);
sigma_seasonal ~ normal(0, 0.1);
sigma_h ~ normal(0, 0.5);
level_1 ~ normal(y_std[1], 1);
trend_1 ~ normal(0, 0.1);
season_init ~ normal(0, 0.5);
h_1 ~ normal(2 * log(0.01), 2); // log-variance of slope innovations
z_level ~ std_normal();
z_trend ~ std_normal();
z_h ~ std_normal();
z_season ~ std_normal();
y_std ~ normal(level + seasonal, sigma_obs);
}
generated quantities {
vector[N] level_orig = y_mean + y_sd * level;
vector[N] trend_orig = y_sd * trend;
vector[N] seasonal_orig = y_sd * seasonal;
vector[N] h_orig = h + 2 * log(y_sd); // log-variance on the original scale
real sigma_obs_orig = y_sd * sigma_obs;
vector[N] y_rep = to_vector(normal_rng(level_orig + seasonal_orig, sigma_obs_orig));
}
The logarithm of the slope variance is shifted by $2\log(\text{sd}(y))$ when moving back to the original scale. h_orig is therefore directly comparable with the h used in the simulation.
Fitting the Models in R
Each model is fitted to its own series. The volatility model gets a higher adapt_delta, since its posterior geometry is harder.
library(rstan)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
fit_basic <- stan(
file = "bsts_basic.stan",
data = list(N = length(y_basic), m = 12, y = y_basic),
chains = 4, iter = 2000, warmup = 1000, seed = 42,
control = list(adapt_delta = 0.95)
)
fit_sv <- stan(
file = "bsts_sv.stan",
data = list(N = length(y_sv), m = 12, y = y_sv),
chains = 4, iter = 3000, warmup = 1500, seed = 123,
control = list(adapt_delta = 0.97)
)
print(fit_basic, pars = c("sigma_obs_orig", "sigma_level", "sigma_trend", "sigma_seasonal"))
print(fit_sv, pars = c("sigma_obs_orig", "sigma_level", "sigma_seasonal", "sigma_h"))
check_hmc_diagnostics(fit_basic)
check_hmc_diagnostics(fit_sv)
Check $\hat{R}$, effective sample sizes, and divergences before interpreting any component. State-space models with several variance parameters can be slow to mix, and a clean-looking plot from a poorly converged chain would be misleading.
Note that sigma_level and sigma_trend are on the standardised scale. To compare them with the simulation values (0.3 and 0.05), multiply them by sd(y).
Visualising the Posterior
A small helper keeps the summaries short. In the plots below, the true simulated components are drawn alongside the posterior estimates, which is possible only because the data are synthetic.
library(ggplot2)
post_basic <- rstan::extract(fit_basic)
post_sv <- rstan::extract(fit_sv)
summ <- function(x) data.frame(
mean = apply(x, 2, mean),
lower = apply(x, 2, quantile, probs = 0.025),
upper = apply(x, 2, quantile, probs = 0.975)
)
plot_component <- function(draws, truth, title, ylab) {
s <- summ(draws)
s$time <- seq_len(nrow(s))
s$truth <- truth
ggplot(s, aes(time)) +
geom_ribbon(aes(ymin = lower, ymax = upper), fill = "steelblue", alpha = 0.25) +
geom_line(aes(y = mean), color = "firebrick", linewidth = 0.9) +
geom_line(aes(y = truth), color = "black", linetype = "dashed") +
labs(title = title, subtitle = "Posterior mean, 95% interval, and true value (dashed)",
x = "Time", y = ylab) +
theme_minimal()
}
plot_component(post_basic$level_orig, truth_basic$level,
"Level: basic BSTS", "Level")
plot_component(post_basic$trend_orig, truth_basic$trend,
"Slope: basic BSTS", "Slope")
plot_component(post_basic$seasonal_orig, truth_basic$season,
"Seasonal effect: basic BSTS", "Seasonal effect")
For the stochastic-volatility model, the quantity of interest is the log-variance path and its effect on the slope:
plot_component(post_sv$h_orig, truth_sv$h,
"Log-variance of slope innovations", "h[t]")
plot_component(post_sv$trend_orig, truth_sv$trend,
"Slope: stochastic-volatility BSTS", "Slope")
The slope and level are usually well recovered, while the h[t] path will probably have wide intervals. That outcome reflects how little information a short series carries about second-order dynamics, and it shouldn’t be read as a failure of the sampler. If the interval for h[t] is wide enough to contain a flat line, the data do not clearly favour time-varying volatility over the basic model.
As a further check, compare the two models on the volatility series with y_rep and a predictive criterion such as PSIS-LOO (this needs a log_lik block, which we omitted for brevity).
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.