Bayesian Exponential Smoothing and Holt-Winters Models

In the previous post, we worked through Bayesian AR, ARMA, and ARIMA models. Those models describe a time series through its relationship with past observations and past innovations. Exponential smoothing takes a different but closely related perspective: instead of modeling lags directly, it recursively updates unobserved components such as a local level, a trend, and a seasonal pattern.

This is often a natural way to think about forecasting problems. A sales series may have a changing baseline, a gradual upward trend, and a recurring monthly pattern. Exponential smoothing separates these features into interpretable components and gives more weight to recent observations than older ones. In a Bayesian formulation, we retain that familiar recursive structure while placing probability distributions on the unknown smoothing parameters and observation noise.

We will implement three increasingly flexible models using Stan and RStan:

1. Bayesian Simple Exponential Smoothing (SES) for a changing level.
2. Bayesian Holt’s linear trend model for level and trend.
3. Bayesian additive Holt-Winters model for level, trend, and seasonality.

The implementations below use a conditional likelihood: we estimate the initial state and then recursively calculate the remaining states from observed data and model parameters. This approach is computationally efficient in Stan because the recursive states are deterministic functions of the parameters rather than a large collection of latent variables that must each be sampled.

Bayesian Simple Exponential Smoothing

Simple Exponential Smoothing is appropriate when a series fluctuates around a level that may change gradually over time, but has no persistent trend or seasonal cycle. Let $\ell_t$ denote the local level at time $t$. The one-step-ahead forecast is simply the current level:

After observing $y_t$, we update the level according to

where $0 < \alpha < 1$ is the smoothing parameter. When $\alpha$ is close to one, the level responds sharply to the newest observation. When $\alpha$ is close to zero, the model smooths aggressively and retains more memory of earlier values.

The update equation can be expanded recursively:

This makes the name “exponential smoothing” literal: observations further in the past receive exponentially decreasing weight.

We first simulate a series from this model. The latent level is updated after each noisy observation, which creates a sequence that is locally stable but still able to adapt.

set.seed(1)

n <- 100
alpha_true <- 0.30
sigma_true <- 1.00

level_true <- numeric(n)
y <- numeric(n)

level_true[1] <- 10
y[1] <- rnorm(1, mean = level_true[1], sd = sigma_true)

for (t in 2:n) {
  level_true[t] <- alpha_true * y[t - 1] +
    (1 - alpha_true) * level_true[t - 1]

  y[t] <- rnorm(1, mean = level_true[t], sd = sigma_true)
}

ts.plot(y, main = "Simulated Simple Exponential Smoothing Series")

The Stan model estimates the initial level $\ell_1$, the smoothing parameter $\alpha$, and the observation scale $\sigma$. The remaining level values are calculated in a transformed parameters block. The likelihood begins at $t=2$, because the first observation is used to condition the recursive model.

data {
  int<lower=2> N;
  vector[N] y;
}

parameters {
  real level_1;
  real<lower=0, upper=1> alpha;
  real<lower=0> sigma;
}

transformed parameters {
  vector[N] level;

  level[1] = level_1;

  for (t in 2:N) {
    level[t] = alpha * y[t - 1] +
      (1 - alpha) * level[t - 1];
  }
}

model {
  level_1 ~ normal(y[1], 5 * sd(y));
  alpha ~ beta(2, 2);
  sigma ~ exponential(1);

  y[2:N] ~ normal(level[2:N], sigma);
}

generated quantities {
  vector[N] y_rep;
  vector[N] one_step_mean;

  one_step_mean[1] = level_1;
  y_rep[1] = normal_rng(one_step_mean[1], sigma);

  for (t in 2:N) {
    one_step_mean[t] = level[t];
    y_rep[t] = normal_rng(one_step_mean[t], sigma);
  }
}

The beta(2, 2) prior is weakly regularizing: it permits values throughout $[0,1]$ but mildly discourages the two extreme cases where the level either ignores nearly all new information or follows each new observation almost exactly. The posterior may still concentrate near either boundary if the data support it.

We can fit the model from R by saving the Stan program as ses_model.stan.

library(rstan)
library(ggplot2)

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

data_list <- list(
  N = length(y),
  y = y
)

fit_ses <- stan(
  file = "ses_model.stan",
  data = data_list,
  chains = 4,
  iter = 2000,
  warmup = 1000,
  seed = 123
)

print(fit_ses, pars = c("level_1", "alpha", "sigma"))

The posterior draws for level give an uncertainty interval for the smoothed signal, rather than a single fitted curve. We can compare the posterior level estimates with the simulated level used to generate the data.

post_ses <- rstan::extract(fit_ses)

level_mean <- apply(post_ses$level, 2, mean)
level_lower <- apply(post_ses$level, 2, quantile, probs = 0.025)
level_upper <- apply(post_ses$level, 2, quantile, probs = 0.975)

ses_df <- data.frame(
  time = 1:n,
  observed = y,
  true_level = level_true,
  level_mean = level_mean,
  level_lower = level_lower,
  level_upper = level_upper
)

ggplot(ses_df, aes(x = time)) +
  geom_line(aes(y = observed), color = "black", linetype = "dashed") +
  geom_line(aes(y = true_level), color = "steelblue", linewidth = 0.8) +
  geom_ribbon(
    aes(ymin = level_lower, ymax = level_upper),
    fill = "firebrick",
    alpha = 0.20
  ) +
  geom_line(aes(y = level_mean), color = "firebrick", linewidth = 0.9) +
  labs(
    title = "Bayesian Simple Exponential Smoothing",
    subtitle = "Posterior local-level estimate and 95% credible interval",
    x = "Time",
    y = "Value"
  ) +
  theme_minimal()

The dashed line represents the observed data, the blue line is the latent level used in simulation, and the red line is the posterior mean estimate. The shaded region is a 95% pointwise credible interval. Notice that this uncertainty is about the estimated smoothed level conditional on the observed series; posterior predictive draws such as y_rep also include future observation noise.

Bayesian Holt’s Trend Model

SES is deliberately conservative: it interprets sustained increases or decreases as a changing level, but it has no separate representation for trend. Holt’s method extends the model by adding a trend component $b_t$. The local level and trend are updated recursively as

The one-step-ahead observation model is

Here, $\alpha$ controls how quickly the level responds to new information, while $\beta$ determines how quickly the estimated slope changes. A small $\beta$ produces a stable trend estimate; a large $\beta$ lets the slope react rapidly to local movements in the series.

The following code simulates a series with a positive underlying trend.

set.seed(2)

n <- 100
alpha_true <- 0.30
beta_true <- 0.10
sigma_true <- 1.00

level_true <- numeric(n)
trend_true <- numeric(n)
y <- numeric(n)

level_true[1] <- 5
trend_true[1] <- 0.50
y[1] <- rnorm(
  1,
  mean = level_true[1] + trend_true[1],
  sd = sigma_true
)

for (t in 2:n) {
  level_true[t] <- alpha_true * y[t - 1] +
    (1 - alpha_true) * (level_true[t - 1] + trend_true[t - 1])

  trend_true[t] <- beta_true * (level_true[t] - level_true[t - 1]) +
    (1 - beta_true) * trend_true[t - 1]

  y[t] <- rnorm(
    1,
    mean = level_true[t] + trend_true[t],
    sd = sigma_true
  )
}

ts.plot(y, main = "Simulated Holt Linear Trend Series")

In the Stan implementation, the trend and level are deterministic recursive quantities. We estimate only their initial values, the two smoothing parameters, and the observation noise scale.

data {
  int<lower=2> N;
  vector[N] y;
}

parameters {
  real level_1;
  real trend_1;
  real<lower=0, upper=1> alpha;
  real<lower=0, upper=1> beta;
  real<lower=0> sigma;
}

transformed parameters {
  vector[N] level;
  vector[N] trend;
  vector[N] forecast_mean;

  level[1] = level_1;
  trend[1] = trend_1;
  forecast_mean[1] = level[1] + trend[1];

  for (t in 2:N) {
    level[t] = alpha * y[t - 1] +
      (1 - alpha) * (level[t - 1] + trend[t - 1]);

    trend[t] = beta * (level[t] - level[t - 1]) +
      (1 - beta) * trend[t - 1];

    forecast_mean[t] = level[t] + trend[t];
  }
}

model {
  level_1 ~ normal(y[1], 5 * sd(y));
  trend_1 ~ normal(0, sd(y));
  alpha ~ beta(2, 2);
  beta ~ beta(2, 2);
  sigma ~ exponential(1);

  y[2:N] ~ normal(forecast_mean[2:N], sigma);
}

generated quantities {
  vector[N] y_rep;

  y_rep[1] = normal_rng(forecast_mean[1], sigma);

  for (t in 2:N) {
    y_rep[t] = normal_rng(forecast_mean[t], sigma);
  }
}

The model uses the filtered state at time $t$ to define forecast_mean[t], which represents the expected next observation after the level and trend updates based on data through time $t-1$. This indexing convention is important. Exponential smoothing models are often written with slightly different time indices across textbooks and software packages, but the central idea remains the same: forecast first, observe the new value, then update the state.

Fit the model in R as follows.

data_list <- list(
  N = length(y),
  y = y
)

fit_holt <- stan(
  file = "holt_model.stan",
  data = data_list,
  chains = 4,
  iter = 2000,
  warmup = 1000,
  seed = 456
)

print(
  fit_holt,
  pars = c("level_1", "trend_1", "alpha", "beta", "sigma")
)

We can plot the posterior mean of the fitted signal, $\ell_t + b_t$, along with its credible interval.

post_holt <- rstan::extract(fit_holt)

forecast_mean <- apply(post_holt$forecast_mean, 2, mean)
forecast_lower <- apply(
  post_holt$forecast_mean,
  2,
  quantile,
  probs = 0.025
)
forecast_upper <- apply(
  post_holt$forecast_mean,
  2,
  quantile,
  probs = 0.975
)

holt_df <- data.frame(
  time = 1:n,
  observed = y,
  true_signal = level_true + trend_true,
  fitted_mean = forecast_mean,
  fitted_lower = forecast_lower,
  fitted_upper = forecast_upper
)

ggplot(holt_df, aes(x = time)) +
  geom_line(aes(y = observed), color = "black", linetype = "dashed") +
  geom_line(aes(y = true_signal), color = "steelblue", linewidth = 0.8) +
  geom_ribbon(
    aes(ymin = fitted_lower, ymax = fitted_upper),
    fill = "firebrick",
    alpha = 0.20
  ) +
  geom_line(aes(y = fitted_mean), color = "firebrick", linewidth = 0.9) +
  labs(
    title = "Bayesian Holt Linear Trend Model",
    subtitle = "Posterior fitted signal and 95% credible interval",
    x = "Time",
    y = "Value"
  ) +
  theme_minimal()

A trend model can produce implausibly large long-horizon forecasts if a short-run increase is mistaken for a permanent slope. This is one reason that damped-trend variants are common in practical forecasting. In a Bayesian model, posterior uncertainty in the trend parameter also becomes especially important as the forecast horizon grows.

Bayesian Seasonal Holt-Winters Model

Many series contain repeating patterns that cannot be represented by a level and linear trend alone. Monthly sales may rise every December, electricity demand may follow an annual cycle, and web traffic may vary systematically by day of the week. The additive Holt-Winters model adds a seasonal state $s_t$ with period $m$.

For an additive seasonal model, the one-step forecast is

After observing $y_t$, the model updates its components according to

The seasonal parameter $\gamma$ controls how quickly the seasonal pattern is allowed to evolve. As with the other smoothing parameters, values close to zero imply a stable, slowly changing component, while values close to one make the component highly responsive to recent observations.

We simulate monthly data with a twelve-period seasonal cycle. The seasonal pattern is additive, so its magnitude remains roughly constant as the level of the series changes.

set.seed(3)

n <- 120
m <- 12

alpha_true <- 0.30
beta_true <- 0.10
gamma_true <- 0.20
sigma_true <- 1.00

level_true <- numeric(n)
trend_true <- numeric(n)
season_true <- numeric(n)
y <- numeric(n)

initial_season <- 2 * sin(2 * pi * (1:m) / m)

level_true[1] <- 10
trend_true[1] <- 0.20
season_true[1:m] <- initial_season

y[1] <- rnorm(
  1,
  mean = level_true[1] + trend_true[1] + season_true[1],
  sd = sigma_true
)

for (t in 2:n) {
  seasonal_lag <- if (t > m) season_true[t - m] else season_true[t]

  level_true[t] <- alpha_true * (y[t - 1] - seasonal_lag) +
    (1 - alpha_true) * (level_true[t - 1] + trend_true[t - 1])

  trend_true[t] <- beta_true * (level_true[t] - level_true[t - 1]) +
    (1 - beta_true) * trend_true[t - 1]

  if (t > m) {
    season_true[t] <- gamma_true *
      (y[t - 1] - level_true[t - 1] - trend_true[t - 1]) +
      (1 - gamma_true) * season_true[t - m]
  }

  current_season <- if (t <= m) season_true[t] else season_true[t - m]

  y[t] <- rnorm(
    1,
    mean = level_true[t] + trend_true[t] + current_season,
    sd = sigma_true
  )
}

ts.plot(y, main = "Simulated Additive Seasonal Holt-Winters Series")

Seasonal states and the level are not separately identifiable without a constraint: adding a constant to every seasonal effect and subtracting the same constant from the level gives the same fitted values. We address this by constraining the initial seasonal factors to sum to zero. Stan represents the first $m-1$ initial factors as parameters and constructs the final factor so that the complete initial seasonal cycle has mean zero.

data {
  int<lower=2> N;
  int<lower=2> m;
  vector[N] y;
}

parameters {
  real level_1;
  real trend_1;
  vector[m - 1] season_init_free;

  real<lower=0, upper=1> alpha;
  real<lower=0, upper=1> beta;
  real<lower=0, upper=1> gamma;
  real<lower=0> sigma;
}

transformed parameters {
  vector[m] season_init;
  vector[N] level;
  vector[N] trend;
  vector[N] season;
  vector[N] forecast_mean;

  for (j in 1:(m - 1)) {
    season_init[j] = season_init_free[j];
  }
  season_init[m] = -sum(season_init_free);

  level[1] = level_1;
  trend[1] = trend_1;
  season[1] = season_init[1];
  forecast_mean[1] = level[1] + trend[1] + season[1];

  for (t in 2:N) {
    real seasonal_lag;

    if (t <= m) {
      season[t] = season_init[t];
      seasonal_lag = season[t];
    } else {
      seasonal_lag = season[t - m];

      season[t] = gamma *
        (y[t - 1] - level[t - 1] - trend[t - 1]) +
        (1 - gamma) * season[t - m];
    }

    level[t] = alpha * (y[t - 1] - seasonal_lag) +
      (1 - alpha) * (level[t - 1] + trend[t - 1]);

    trend[t] = beta * (level[t] - level[t - 1]) +
      (1 - beta) * trend[t - 1];

    forecast_mean[t] = level[t] + trend[t] + seasonal_lag;
  }
}

model {
  level_1 ~ normal(y[1], 5 * sd(y));
  trend_1 ~ normal(0, sd(y));
  season_init_free ~ normal(0, 2 * sd(y));

  alpha ~ beta(2, 2);
  beta ~ beta(2, 2);
  gamma ~ beta(2, 2);
  sigma ~ exponential(1);

  y[2:N] ~ normal(forecast_mean[2:N], sigma);
}

generated quantities {
  vector[N] y_rep;

  y_rep[1] = normal_rng(forecast_mean[1], sigma);

  for (t in 2:N) {
    y_rep[t] = normal_rng(forecast_mean[t], sigma);
  }
}

This model conditions on the first seasonal cycle through the estimated initial seasonal factors. In real applications, it is generally helpful to have several complete seasonal cycles. For example, a monthly model with annual seasonality should ideally be fit to multiple years of observations. With only one or two cycles, the model may have difficulty distinguishing a genuine recurring seasonal pattern from a one-time change in level or trend.

We fit the model using m = 12 for monthly data.

data_list <- list(
  N = length(y),
  m = m,
  y = y
)

fit_hw <- stan(
  file = "seasonal_holt_winters.stan",
  data = data_list,
  chains = 4,
  iter = 3000,
  warmup = 1500,
  seed = 789,
  control = list(adapt_delta = 0.95)
)

print(
  fit_hw,
  pars = c(
    "level_1",
    "trend_1",
    "alpha",
    "beta",
    "gamma",
    "sigma"
  )
)

The slightly higher adapt_delta asks Stan to use a more conservative Hamiltonian Monte Carlo step size. Recursive models with several smoothing parameters can have correlated posterior geometries, especially when the series is short or when trend and seasonality are weakly identified.

Finally, we visualize the fitted seasonal signal and its posterior uncertainty.

post_hw <- rstan::extract(fit_hw)

fitted_mean <- apply(post_hw$forecast_mean, 2, mean)
fitted_lower <- apply(
  post_hw$forecast_mean,
  2,
  quantile,
  probs = 0.025
)
fitted_upper <- apply(
  post_hw$forecast_mean,
  2,
  quantile,
  probs = 0.975
)

seasonal_signal_true <- level_true + trend_true + season_true

hw_df <- data.frame(
  time = 1:n,
  observed = y,
  true_signal = seasonal_signal_true,
  fitted_mean = fitted_mean,
  fitted_lower = fitted_lower,
  fitted_upper = fitted_upper
)

ggplot(hw_df, aes(x = time)) +
  geom_line(aes(y = observed), color = "black", linetype = "dashed") +
  geom_line(aes(y = true_signal), color = "steelblue", linewidth = 0.8) +
  geom_ribbon(
    aes(ymin = fitted_lower, ymax = fitted_upper),
    fill = "firebrick",
    alpha = 0.20
  ) +
  geom_line(aes(y = fitted_mean), color = "firebrick", linewidth = 0.9) +
  labs(
    title = "Bayesian Additive Seasonal Holt-Winters Model",
    subtitle = "Posterior fitted signal and 95% credible interval",
    x = "Time",
    y = "Value"
  ) +
  theme_minimal()

Conclusion

Exponential smoothing models provide a useful alternative to ARIMA-style approaches when a series is more naturally described through evolving components than through lag polynomials. Simple Exponential Smoothing estimates a changing local level, Holt’s method adds a trend, and the Holt-Winters extension introduces recurring seasonal effects.

The Bayesian versions preserve the core recursive logic of these models while adding posterior uncertainty for the smoothing parameters, initial states, fitted components, and forecasts. That uncertainty matters in practice: a point forecast may look stable even when the data do not strongly determine whether recent movement is a level shift, a persistent trend, or a temporary seasonal deviation.

These conditional Stan implementations are a useful starting point, but they are not the only Bayesian formulation. A full state-space approach can treat the level, trend, and seasonal terms as stochastic latent states, allowing them to evolve with their own process noise. From there, natural extensions include damped trends, robust Student t observation models, regression effects from external predictors, hierarchical pooling across related time series, and posterior predictive forecasting over future horizons.