Bayesian Time Series

Bayesian time series modeling puts a prior on the parameters of a forecasting model and produces a posterior over them, rather than a single fitted value. Forecasts then come from integrating over that posterior:

\[ p(y_{T+h} \mid y_{1:T}) = \int p(y_{T+h} \mid \theta, y_{1:T})\, p(\theta \mid y_{1:T})\, d\theta \]

The usual justification is that this quantifies uncertainty properly. That is true, and it is worth seeing how much it actually changes.

What the plug-in forecast leaves out

import numpy as np
from scipy.stats import norm, t as student

rng = np.random.default_rng(0)
PHI, SIGMA, T, RUNS = 0.7, 1.0, 60, 20_000

def series(n):
    """AR(1): y_t = phi * y_{t-1} + noise. The parameter is what we do not know."""
    y = np.zeros((n, T + 1))
    for i in range(1, T + 1):
        y[:, i] = PHI * y[:, i - 1] + rng.normal(scale=SIGMA, size=n)
    return y

y = series(RUNS)
x, z = y[:, :-1], y[:, 1:]
sxx = (x ** 2).sum(axis=1)
phi_hat = (x * z).sum(axis=1) / sxx
resid = z - phi_hat[:, None] * x
s2 = (resid ** 2).sum(axis=1) / (T - 1)
print(f"  {RUNS:,} AR(1) series of length {T}, true phi = {PHI}, true sigma = {SIGMA}")
print(f"    estimated phi: mean {phi_hat.mean():.4f}, sd across series {phi_hat.std():.4f}")
"""Conditional on the fixed initial value y0=0, use the joint prior
   p(phi, sigma^2) proportional to 1/sigma^2. Integrating sigma^2 gives
   a Student-t posterior for phi with T-1 degrees of freedom."""

future = PHI * y[:, -1] + rng.normal(scale=SIGMA, size=RUNS)
print(f"{'interval':>34}{'coverage':>11}{'mean width':>13}")
lo = phi_hat * y[:, -1] - 1.96 * np.sqrt(s2)
hi = phi_hat * y[:, -1] + 1.96 * np.sqrt(s2)
print(f"{'plug-in: phi_hat treated as known':>34}"
      f"{float(((future >= lo) & (future <= hi)).mean()):>11.4f}{(hi - lo).mean():>13.4f}")
se_phi = np.sqrt(s2 / sxx)
scale = np.sqrt(s2 + (se_phi * y[:, -1]) ** 2)
crit = student(T - 1).ppf(0.975)
lo2 = phi_hat * y[:, -1] - crit * scale
hi2 = phi_hat * y[:, -1] + crit * scale
print(f"{'posterior predictive':>34}"
      f"{float(((future >= lo2) & (future <= hi2)).mean()):>11.4f}{(hi2 - lo2).mean():>13.4f}")
print(f"{'the truth, if phi were known':>34}"
      f"{float((np.abs(future - PHI * y[:, -1]) <= 1.96 * SIGMA).mean()):>11.4f}"
      f"{2 * 1.96 * SIGMA:>13.4f}")

i = 0
last, ph, sp, sd = y[i, -1], phi_hat[i], se_phi[i], np.sqrt(s2[i])
print(f"    that series: phi_hat {ph:.4f}, its standard error {sp:.4f}, last value {last:.4f}")
print(f"{'horizon':>9}{'sd, phi fixed':>16}{'sd, phi uncertain':>20}{'ratio':>8}")
PATHS = 200_000
for h in (1, 2, 5, 10, 20):
    r = np.random.default_rng(h)
    fixed = np.full(PATHS, last)
    drawn = np.full(PATHS, last)
    phis = ph + sp * r.standard_t(T - 1, size=PATHS)
    for _ in range(h):
        e1 = r.normal(scale=sd, size=PATHS)
        fixed = ph * fixed + e1
        drawn = phis * drawn + r.normal(scale=sd, size=PATHS)
    print(f"{h:>9}{fixed.std():>16.4f}{drawn.std():>20.4f}{drawn.std() / fixed.std():>8.3f}")
#   20,000 AR(1) series of length 60, true phi = 0.7, true sigma = 1.0
#     estimated phi: mean 0.6783, sd across series 0.0971
#                           interval   coverage   mean width
#  plug-in: phi_hat treated as known     0.9444       3.9054
#               posterior predictive     0.9513       4.0223
#       the truth, if phi were known     0.9521       3.9200
#     that series: phi_hat 0.7278, its standard error 0.0896, last value 0.6364
#   horizon   sd, phi fixed   sd, phi uncertain   ratio
#         1          0.9390              0.9387   1.000
#         2          1.1642              1.1705   1.005
#         5          1.3424              1.3812   1.029
#        10          1.3685              1.4504   1.060
#        20          1.3690              1.4977   1.094

Twenty thousand AR(1) series of length 60, with a true coefficient of 0.7. The estimate averages 0.6783 with a standard deviation of 0.0971 across series — a real spread that a plug-in forecast treats as zero.

The one-step intervals show the consequence, and show it honestly: the plug-in interval covers 0.9444 against a nominal 95%, the posterior predictive covers 0.9513, and an interval built with the true coefficient covers 0.9521. The correction is real and it is small — at one step, parameter uncertainty is minor next to the noise.

The second block is where it stops being minor. For this multi-step comparison, the innovation standard deviation is fixed at its fitted value, while phi is either fixed or drawn from its marginal posterior. This isolates coefficient uncertainty rather than sampling the full joint posterior predictive distribution. The displayed standard deviations agree to two decimals at one step, and diverge with the horizon: a ratio of 1.029 at five steps, 1.060 at ten, 1.094 at twenty.

The displayed widening belongs to this fitted series and horizon range. In a stationary AR(1), the initial-state contribution is \(\phi^h y_T\); its derivative \(h\phi^{h-1}y_T\) can eventually shrink. Parameter uncertainty also changes the distribution of accumulated innovations. Its total contribution depends on the posterior, stability and horizon, and need not grow monotonically forever.

What else the posterior buys

Three things that are awkward to get otherwise, and all of them follow from having a distribution rather than a fit.

• Structure as a prior. Smoothness, a trend that should not reverse, a seasonal pattern with a known period — these are statements about parameters, and a prior is where they go. A structural time series model is exactly this: level, trend and seasonality as components with priors on how fast each may move.
• Missing data and irregular sampling are handled by the same machinery, since an unobserved value is just another quantity to integrate over rather than a gap needing imputation first.
• Online updating comes free, because yesterday's posterior is a valid prior for today's observation.

The forecast output is also a full predictive distribution rather than an interval bolted on, so a quantile the business actually needs — the level demand exceeds 5% of the time, say — can be read directly rather than approximated from a symmetric band.

What it costs

Computation depends on the model and inference method. The example uses conjugate calculations; some state-space models also admit efficient analytic filtering. Models fitted by sampling can be much more expensive and need sampler diagnostics alongside residual checks. Priors and likelihood assumptions remain part of the forecast in every case.

Against that, the measurement above sets the benefit honestly: a coverage correction of under a point at one step, and a multi-step standard-deviation ratio of 1.094 at twenty for the selected series. In this simulation, the one-step coverage difference is small. If you forecast far ahead, or the interval drives an inventory or capacity decision, the correction is the thing you were trying to get right.

References

• SciPy documentation, scipy.stats.t and scipy.stats.norm. docs.scipy.org