Bayesian Forecasting

Definition and motivation

Bayesian forecasting is the process of producing probability distributions for future observations of a time series by conditioning on historical data and integrating over posterior uncertainty in both parameters and latent states. Formally, for observed data $y_{1:T}$ and future horizon $h\ge 1$, the posterior predictive distribution is

This identity is the backbone: it separates aleatory uncertainty (conditional observation noise given states and parameters) from epistemic uncertainty (uncertainty about parameters $\theta$ and future latent states $x_{T+h}$).

Why care: uncertainty quantification, decision-making, and model comparison

• •Inventory/demand applications require not just a point forecast but tail probabilities to set stock levels and service-level constraints.
• •Bayesian forecasts naturally produce predictive intervals, predictive densities, and full probabilistic statements useful in expected-cost decisions.
• •Because we explicitly model parameter and structural uncertainty, Bayesian forecasting supports Bayesian model averaging (BMA) to hedge against model misspecification.

Relation to prerequisites

• •In State-Space Models, we learned how the Kalman filter and smoothing produce filtering and smoothing distributions for the latent state when parameters are known; Bayesian forecasting builds on that by treating parameters as uncertain and integrating them out.
• •In Hierarchical Bayesian Models, we learned partial pooling and hyperpriors; in forecasting this allows sharing information across related series (e.g., stores) so that forecasts borrow strength.
• •In Conjugate Priors, we learned closed-form posterior updates (e.g., Gamma-Poisson); those conjugate cases produce closed-form posterior predictives, which are computationally cheap and interpretable.

Concrete simple formula and numeric example

Consider IID Poisson counts $y_t\sim\mathrm{Pois}(\lambda)$ with prior $\lambda\sim\mathrm{Gamma}(\alpha,\beta)$ (shape-rate). The posterior is $\lambda\mid y_{1:T}\sim\mathrm{Gamma}(\alpha+\sum y_t,\ \beta+T)$ and the posterior predictive for a future single count is Negative Binomial (Gamma-Poisson mixture):

where $S=\sum_{t=1}^T y_t$. Numeric example: if $T=10$, $S=50$, prior $\alpha=2,\beta=1$, then posterior shape $=52$, rate $=11$. The predictive mean is $(\alpha+S)/(\beta+T) = 52/11 \approx 4.727$ and the predictive variance is higher than the mean due to parameter uncertainty; the distribution is Negative Binomial with those parameters.

Key conceptual point

The posterior predictive is the unique object you should use for decision-making: it already accounts for parameter and latent-state uncertainty. The different computational strategies (closed-form conjugate, analytic Kalman-filter propagation, Monte Carlo integration, Sequential Monte Carlo) are simply ways to evaluate or approximate this integral.

Two canonical, illuminating cases give intuition and practical recipes: (A) conjugate observation models (e.g., Poisson-Gamma, Normal-Inverse-Gamma) and (B) linear Gaussian state-space models where Kalman filter analytic formulas produce predictive distributions conditional on parameters. Both yield explicit predictive forms that we can extend to parameter uncertainty by analytic marginalization or Monte Carlo.

A. Conjugate observation models (closed-form predictive)

If the observation model and prior are conjugate, integration over parameters yields closed-form posterior predictives. Example: Normal mean with unknown variance (Normal-Inverse-Gamma conjugacy). Suppose

In this case the one-step posterior predictive is a Student-t:

with posterior hyperparameters $\kappa_T=\kappa_0+T$, $\mu_T=(\kappa_0\mu_0+T\bar y)/\kappa_T$, $\alpha_T=\alpha_0+T/2$, $\beta_T=\beta_0+\tfrac12\sum (y_t-\bar y)^2+\tfrac{\kappa_0 T(\bar y-\mu_0)^2}{2\kappa_T}$. Numeric example: prior $\mu_0=0,\kappa_0=1,\alpha_0=2,\beta_0=2$, data $T=4$ with $y=(3,4,2,5)$ so $\bar y=3.5$, sample variance $\sum(y_t-\bar y)^2=5$. Then $\kappa_T=5$, $\mu_T=(1*0+4*3.5)/5=2.8$, $\alpha_T=4$, and $\beta_T=2+2.5+\tfrac{1*4*(3.5-0)^2}{2*5}=4.5+\tfrac{49}{25}=6.46$ (approx). The predictive is approximately $t_8(2.8,\:\text{scale}^2\approx 6.46*(6/5)/(4*5))$ — compute numerically to get predictive mean 2.8 and heavier tails than Gaussian.

Why this matters: the predictive Student-t reflects uncertainty in both mean and variance; credible intervals are wider than plug-in intervals.

B. Linear-Gaussian state-space: Kalman predictive and marginalization over parameters

Consider the Dynamic Linear Model (DLM):

State evolution: $x_{t+1}=F_t x_t + w_t$, $w_t\sim\mathcal N(0,W_t)$

Observation: $y_t=H_t x_t + v_t$, $v_t\sim\mathcal N(0,V_t)$

If parameters $(F_t,H_t,W_t,V_t)$ are known, the Kalman filter gives the predictive distribution for $y_{T+1}$:

Concrete numeric example: local level model with scalar state $x_{t+1}=x_t+\eta_t$ ($W=0.5$), observations $y_t=x_t+\epsilon_t$ ($V=1$), prior $x_0\sim\mathcal N(0,1)$. After observing $T=3$ data $y=(1.2,0.9,1.5)$, the Kalman filter produces $m_{3\mid3}\approx1.2$, $C_{3\mid3}\approx0.4$ and the one-step ahead predictive mean for $y_4$ is $1.2$ with variance $C_{3\mid3}+V\approx0.4+1=1.4$.

Marginalizing over parameters: combine Kalman conditionals with posterior over parameters. If unknown variances have conjugate priors (Inverse-Gamma), marginal predictive can be Student-t-like; otherwise sample from $p(\theta\mid y)$ (via MCMC) and average Kalman predictive densities:

This Monte Carlo marginalization is simple and effective: run MCMC for parameters, and for each draw run a Kalman update to compute the conditional predictive, then average or simulate draws.

Tips for practice

• •Always distinguish predictive variance due to observation noise vs parameter uncertainty; report both when you can.
• •Use conjugacy for fast benchmarks and sanity checks; use MCMC/SMC for nonconjugate structural models.

Structural time series models and their expressiveness

A structural time series specifies interpretable latent components: level, trend, seasonality, regression effects, interventions, and time-varying coefficients. For example, a commonly used local linear trend plus seasonality model for scalar observations is

State vector: $x_t=[\ell_t,\, b_t,\, s_t^1,\dots,s_t^{S-1}]^\top$ with evolution

Observation: $y_t=\ell_t + s_t + Z_t \beta + v_t$ where $Z_t\beta$ are regression covariates with coefficients $\beta$. This structure is richer than ARIMA: components have direct interpretation, which is invaluable in business contexts (trend vs seasonal vs promotion lift).

Parameter uncertainty and hierarchical priors

In practice you often have multiple related time series (SKU/store pairs). In Hierarchical Bayesian Models we learned partial pooling — use a hierarchical prior on coefficients $\beta_j$ or innovation variances to borrow strength. For example, for store-specific promotion effect $\beta_i\sim\mathcal N(\mu_\beta,\tau_\beta^2)$ with hyperpriors on $(\mu_\beta,\tau_\beta)$. Posterior sharing reduces overfitting for low-data series while allowing high-data series to dominate.

Computational strategies: MCMC, SMC, and particle MCMC

Structural models rarely conjugate fully. Two common approaches:

1) MCMC with forward-filter backward-sample (FFBS): sample parameters $\theta$ (e.g., variances, regression coefficients) using Metropolis or Gibbs; conditionally, sample the latent states via the Kalman smoother (or FFBS for linear-Gaussian) to produce draws from $p(x_{1:T}\mid y_{1:T},\theta)$. Each sample yields a draw from the joint posterior; simulate future $y_{T+h}$ by forward-propagating the state and adding observation noise.

2) Sequential Monte Carlo / Particle Filters: useful for online updating and non-linear/non-Gaussian models. Particle MCMC (e.g., Particle Gibbs) provides exact (up to Monte Carlo error) draws from the joint posterior.

Numeric sketch: FFBS for DLM

Given parameter draw $\theta^{(i)}$, run Kalman filter forward to obtain filtering distributions, then run backward sampling to draw $x_{1:T}^{(i)}$. Forward-simulate $x_{T+1}^{(i)}$ and $y_{T+1}^{(i)}$. Repeat across $i=1\dots M$ samples to obtain Monte Carlo predictive draws. Example numbers: 1000 posterior draws produce 1000 predictive draws, from which empirical predictive mean, 90% interval, and predictive quantiles are computed.

Bayesian model averaging (BMA) and model selection

BMA weights models by posterior model probability:

with

The marginal likelihood $p(y\mid M_k)$ is often computed by Laplace approximation, bridge sampling, or thermodynamic integration. Numeric example: two models with log marginal likelihoods $\log p(y\mid M_1)=-120$, $\log p(y\mid M_2)=-125$, equal priors $p(M_1)=p(M_2)=0.5$. Then posterior odds favor $M_1$ by factor $e^{5}\approx148.4$, so $p(M_1\mid y)\approx 0.993$ and $p(M_2\mid y)\approx0.007$.

Model averaging vs selection in forecasting

• •BMA accounts for model uncertainty and often yields better calibrated predictive distributions than selecting a single model.
• •If a single model is strongly dominant (as above), BMA reduces to that model. If models have complementary strengths (e.g., one captures promotions, another long-run trend), BMA blends predictive densities to reduce worst-case risk.

Model stacking (an alternative) optimizes weights to minimize a predictive loss (e.g., log-score) on holdout data rather than using marginal likelihoods; it often improves predictive performance in practice when models are misspecified.

Practical warnings and diagnostics

• •Check calibration (PIT histograms, predictive coverage). A well-calibrated model yields uniform PIT and correct empirical coverage of posterior predictive intervals.
• •Inspect predictive tails. For inventory, tail underestimation is costly.
• •Use predictive scoring rules (log score, CRPS) on holdout sets to compare models and weighting schemes.

Demand forecasting is the canonical application where Bayesian forecasting shines: decisions (ordering, pricing, promotions) depend on tails and asymmetric loss. Below I show concrete modeling patterns, loss-based decision rules, and numeric illustrations.

Count data and overdispersion: Poisson-Gamma & Negative Binomial predictive

For counts, start with Poisson likelihood and add hierarchical priors or latent overdispersion. The Gamma-Poisson conjugacy yields a Negative Binomial predictive that captures extra-Poisson variance due to parameter uncertainty.

Numeric worked mini-case: single SKU daily demand

Observations: 30 days with total sales S=150 (mean 5/day). Prior $\lambda\sim\mathrm{Gamma}(2,1)$ (shape-rate). Posterior: $\lambda\mid y\sim\mathrm{Gamma}(152,31)$, posterior mean $152/31\approx4.903$. One-day predictive is Negative Binomial with mean $152/31\approx4.903$ and variance $4.903+4.903^2/(152)\approx 4.903+0.158\approx5.061$. Compare to plug-in Poisson variance 4.903: the predictive variance is slightly larger.

In-store hierarchy: partial pooling for items/stores

Suppose we have 50 stores and want per-store daily rate $\lambda_i$. A hierarchical model

with hyperpriors on $(\alpha,\beta)$ implements partial pooling. In Hierarchical Bayesian Models we learned how posteriors shrink store estimates toward the chain mean proportionally to information; this reduces forecast variance for low-volume stores.

Promotion and regression effects

Add covariates $x_{it}$ (price, promotion flag) in a log-linear rate:

Posterior uncertainty in $\beta_i$ propagates into predictive uncertainty. To compute predictive uplift under a promotion, simulate posterior draws $\beta_i^{(s)}$ and produce predictive draws for future periods under different $Z_{it}$ — this yields posterior distributions over promotional lift.

Inventory decision example: Newsvendor loss

If stock level $q$ incurs underage cost $c_u$ per unit sold but not stocked and overage cost $c_o$ per leftover unit, the Bayes-optimal quantity minimizes posterior expected loss and corresponds to the posterior predictive quantile:

Numeric example: if stockout cost $c_u=10$, leftover cost $c_o=2$, optimal service quantile is $10/(10+2)=0.8333$. If the posterior predictive negative binomial CDF evaluated at $q=9$ is 0.81 and at $q=10$ is 0.87, choose $q^*=10$.

Assessing and communicating uncertainty

• •Report predictive mean, median, and multiple quantiles (e.g., 5%, 50%, 95%). For demand planning, also report expected stockouts given a replenishment policy and the tail probability beyond a large threshold.
• •Use probabilistic forecasts in downstream stochastic optimization (e.g., stochastic inventory, dynamic pricing) rather than plugging in point forecasts.

Scaling considerations

• •For many similar series, exploit hierarchical structure and amortized inference (e.g., variational Bayes, amortized MCMC) to scale.
• •Use model selection/stacking to combine cheap parametric models with expensive simulation-based models; stack weights are fit to maximize predictive performance on held-out recent windows.

Summary of applied recipe

1) Specify structural components you believe matter (trend, seasonality, covariates).

2) Choose priors that express domain knowledge; use hierarchical priors for pools of series.

3) Fit via MCMC/SMC/variational methods; validate predictive calibration on holdout data (coverage, PIT).

4) Use posterior predictive draws for decision tasks (newsvendor quantiles, expected costs, promotion uplift distributions).

This pipeline turns time-series forecasting into a decision-ready probabilistic system that explicitly balances bias-variance and model uncertainty.