Bayesian Stochastic Volatility

Model

Where GARCH makes the conditional variance a deterministic function of past returns, a stochastic-volatility (SV) model gives the log-variance its own latent AR(1) process, driven by an independent shock. Returns are noisy readings of this hidden volatility state. In the Cox (1981) taxonomy, GARCH is observation-driven (vol is known once the data are), SV is parameter-driven (vol has its own innovation and is never observed) — so the fitted volatility path carries a genuine posterior credible band, which a GARCH filter cannot produce. That is the headline difference.

r_t = e^{h_t/2}\,\varepsilon_t, \qquad h_t = \mu + \phi\,(h_{t-1}-\mu) + \sigma_\eta\,\eta_t, \qquad \varepsilon_t,\eta_t \sim N(0,1)

$h_t$ — latent log-variance (the state); inferred, not observed
$\phi$ — persistence of log-volatility (the GARCH $\alpha+\beta$ analogue); $\sigma_\eta$ — vol-of-vol
$\mu$ — mean log-variance; extensions add $\nu$ (fat tails) and $\rho$ (leverage)

From scratch — the state-space view (Kim–Shephard–Chib)

SV is a state-space model, and the state-space toolkit fits it. The from-scratch notebook builds the Kim–Shephard–Chib (1998) sampler by hand: squaring and log-transforming the returns linearises the measurement equation at the cost of a $\log\chi^2_1$ error, which KSC approximate by a 7-component Gaussian mixture. Conditional on the mixture indicators the system is linear-Gaussian, so the entire log-variance path $h_{1:T}$ is drawn in one block by the Kalman filter + backward sampling (Carter–Kohn / Frühwirth-Schnatter FFBS). A Gibbs sweep then cycles $s\mid h$ → $h\mid s,\theta$ (FFBS) → $\theta\mid h$ , with conjugate/Metropolis updates for $(\mu, \phi, \sigma_\eta)$ .

Block	Draw	Method
(1)	$s_t \mid y^*_t, h_t$	Mixture indicators — multinomial over the 7 KSC components
(2)	$h_{1:T} \mid s, \theta$	Whole path in one block — Kalman filter + backward sample (FFBS)
(3)	$\theta \mid h$	$\mu$ Gaussian, $\sigma_\eta^2$ inverse-Gamma, $\phi$ independence Metropolis

Notebooks

On daily S&P 500 returns, three independent engines agree on the base SV posterior: the from-scratch KSC sampler ( $\phi=0.988$ , $\sigma_\eta=0.147$ ), NumPyro/NUTS (non-centered, to defuse the SV funnel: $0.987$ / $0.155$ ), and R's stochvol (Kastner's ASIS interweaving sampler: $0.987$ / $0.152$ ). Volatility is highly persistent ( $\phi\approx 0.99$ ) but genuinely stochastic — $\sigma_\eta\approx 0.15$ is a real vol-of-vol, and the inferred path comes with a credible band that widens through the 1987 and 2008 crises. Extensions (a 2×2 family). Two ingredients on top of the plain SV, the SV analogues of the GARCH extensions: Student-t innovations (SV-t — fat tails survive even under stochastic vol, $\nu\approx 9$ ) and leverage (asymmetric SV — correlate the return and vol shocks, $\rho\approx -0.5\text{ to }-0.7$ , so a negative return pushes next-period volatility up). Ranked by WAIC, SV-t-leverage wins, with leverage the larger gain. This closes the loop on the volatility arc: GARCH, GJR-GARCH and SV all recover the same three stylized facts — persistence, fat tails, and the $\rho<0$ / $\gamma>0$ leverage asymmetry — written three ways (an observation-driven recursion, a latent-state correlation, and the option-pricing skew). Both extension notebooks are cross-validated between NumPyro and stochvol (svtsample / svlsample / svtlsample). SV vs GARCH — the out-of-sample horse race. A final notebook settles the practical question: does the latent-state flexibility actually forecast better? Parameters are estimated once on a training window (1987–2001), then each model's volatility filter is rolled forward across the 2001–2009 test set (including 2008) to produce genuine 1-step-ahead predictive densities — closed-form for GARCH, but for SV via a bootstrap particle filter (a non-linear non-Gaussian state-space model has no closed-form predictive; the particle filter yields a fat-tailed scale-mixture for free). Matched feature-for-feature the two paradigms are a near dead heat on log predictive score (GJR-t best by a whisker, SV-t-leverage 0.004 nats behind), but on VaR calibration SV-t-leverage is the only model clean at both 1% and 5% — its latent-vol dynamics plus a milder $\nu$ get the whole left tail right, not just the extreme.

Python SV from scratch — Kim–Shephard–Chib 7-component mixture + FFBS state-space sampler (Python) Open → Python NumPyro / NUTS — non-centered base SV, the SV funnel defused (Python) Open → R R stochvol (Kastner ASIS interweaving) — two-engine validation of the SV posterior (R) Open → Python Extensions — SV-t & leverage (asymmetric SV), the 2×2 family ranked by WAIC, NumPyro (Python) Open → R Extensions in R stochvol — svtsample / svlsample / svtlsample, independent check (R) Open → Python SV vs GARCH — out-of-sample density forecasts (LPS) + VaR backtest, particle-filter predictive (Python) Open →

Downloads

sv_scratch.py From-scratch KSC (1998) mixture / FFBS stochastic-volatility Gibbs sampler (NumPy) svgarch_forecast.py SV-vs-GARCH out-of-sample engine — rolling filter, particle-filter predictive, VaR backtests (needs sv_scratch.py + garch_scratch.py) sp500ret.csv Daily S&P 500 returns, 1987–2009

References

Taylor, S. J. (1986). Modelling Financial Time Series. Wiley. — the original stochastic-volatility model
Jacquier, E., Polson, N. G. & Rossi, P. E. (1994). Bayesian analysis of stochastic volatility models. Journal of Business & Economic Statistics 12(4), 371–389. — Bayesian inference for the latent-volatility model
Kim, S., Shephard, N. & Chib, S. (1998). Stochastic volatility: likelihood inference and comparison with ARCH models. Review of Economic Studies 65(3), 361–393. — the 7-component mixture / FFBS sampler built from scratch here
Carter, C. K. & Kohn, R. (1994). On Gibbs sampling for state space models. Biometrika 81(3), 541–553. — the forward-filter backward-sampling (FFBS) draw of the whole log-variance path
Kastner, G. & Frühwirth-Schnatter, S. (2014). Ancillarity-sufficiency interweaving strategy (ASIS) for boosting MCMC estimation of stochastic volatility models. Computational Statistics & Data Analysis 76, 408–423. — the interweaving sampler behind R's stochvol
Kastner, G. (2016). Dealing with stochastic volatility in time series using the R package stochvol. Journal of Statistical Software 69(5), 1–30. — the stochvol implementation and its SV-t / leverage extensions

SV Sampler — Source Code

"""
Stochastic volatility from scratch — Kim, Shephard & Chib (1998) mixture sampler.

State-space form of the SV model (returns already de-meaned):
    measurement:  r_t = exp(h_t / 2) * eps_t,      eps_t ~ N(0, 1)
    state (AR1):  h_t = mu + phi (h_{t-1} - mu) + sigma_eta * eta_t,   eta_t ~ N(0, 1)

h_t is the latent LOG-VARIANCE. The model is a NON-LINEAR state-space model, but squaring and
taking logs linearises the measurement equation at the cost of a non-Gaussian error:
    y*_t = log(r_t^2) = h_t + z_t,        z_t = log(eps_t^2) ~ log chi^2_1   (mean -1.2704, var pi^2/2).

Kim-Shephard-Chib approximate z_t by a 7-component Gaussian mixture (indicator s_t in 1..7). GIVEN the
indicators the system is a LINEAR GAUSSIAN state-space model, so the whole path h_{1:T} is drawn in one
block by the Kalman filter + backward sampling (Carter-Kohn / Fruhwirth-Schnatter forward-filter
backward-sample, FFBS). A Gibbs sweep then cycles:  s | h  ->  h | s, theta (FFBS)  ->  theta | h.

References
----------
Taylor (1982, 1986) original SV model; Cox (1981) observation- vs parameter-driven taxonomy;
Jacquier, Polson & Rossi (1994, JBES) Bayesian SV (single-move); Carter & Kohn (1994, Biometrika) and
Fruhwirth-Schnatter (1994, JTSA) FFBS; Kim, Shephard & Chib (1998, RES) the mixture/multi-move sampler
implemented here; Durbin & Koopman (2012) state-space textbook; Omori et al. (2007) 10-component +
leverage; Kastner & Fruhwirth-Schnatter (2014, CSDA) ASIS (the stochvol package).
"""
import numpy as np

# ---- KSC (1998) Table 4: 7-component mixture for log chi^2_1 -------------------------------------
# means are deviations (weighted sum ~ 0); the constant offset -1.2704 is added so the mixture matches
# E[log chi^2_1] = -1.2704 and Var = pi^2/2 = 4.9348 (both reproduced to 4 dp -> internal check).
_Q  = np.array([0.00730, 0.10556, 0.00002, 0.04395, 0.34001, 0.24566, 0.25750])
_M  = np.array([-10.12999, -3.97281, -8.56686, 2.77786, 0.61942, 1.79518, -1.08819])
_V2 = np.array([5.79596, 2.61369, 5.17950, 0.16735, 0.64009, 0.34023, 1.26261])
_OFF = -1.2704


def mixture():
    """Return (probabilities, component means incl. offset, component variances)."""
    return _Q.copy(), _M + _OFF, _V2.copy()


# ---- data generating process (for validation) ---------------------------------------------------
def simulate_sv(T, mu, phi, sigma, seed=0):
    rng = np.random.default_rng(seed)
    h = np.empty(T)
    h[0] = mu + rng.standard_normal() * sigma / np.sqrt(1 - phi**2)     # stationary start
    for t in range(1, T):
        h[t] = mu + phi * (h[t-1] - mu) + sigma * rng.standard_normal()
    r = np.exp(h / 2) * rng.standard_normal(T)
    return r, h


# ---- Gibbs building blocks ----------------------------------------------------------------------
def _sample_indicators(ystar, h, rng):
    """s_t | y*_t, h_t : multinomial over the 7 components (vectorised)."""
    resid = (ystar - h)[:, None] - (_M[None, :] + _OFF)
    logp = np.log(_Q)[None, :] - 0.5 * np.log(2 * np.pi * _V2)[None, :] - 0.5 * resid**2 / _V2[None, :]
    logp -= logp.max(1, keepdims=True)
    p = np.exp(logp); p /= p.sum(1, keepdims=True)
    cdf = np.cumsum(p, 1)
    return (rng.random(len(ystar))[:, None] > cdf).sum(1)               # chosen component index 0..6


def _ffbs(ystar, a_obs, R_obs, mu, phi, sig2, rng):
    """Forward-filter backward-sample: draw h_{1:T} from a linear Gaussian state-space model.

    obs:   y*_t = h_t + a_obs_t + N(0, R_obs_t)
    state: h_t  = (1-phi) mu + phi h_{t-1} + N(0, sig2),   h_1 ~ N(mu, sig2/(1-phi^2)) stationary.
    """
    T = len(ystar); d = (1 - phi) * mu
    af = np.empty(T); Pf = np.empty(T)
    ap, Pp = mu, sig2 / (1 - phi**2)                                   # t=0 prior = stationary
    for t in range(T):
        if t > 0:
            ap = d + phi * af[t-1]; Pp = phi**2 * Pf[t-1] + sig2       # predict
        F = Pp + R_obs[t]; K = Pp / F                                  # update
        af[t] = ap + K * (ystar[t] - a_obs[t] - ap)
        Pf[t] = Pp * (1 - K)
    h = np.empty(T)
    h[T-1] = af[T-1] + np.sqrt(Pf[T-1]) * rng.standard_normal()
    for t in range(T-2, -1, -1):                                       # backward sample
        Pp_next = phi**2 * Pf[t] + sig2
        J = phi * Pf[t] / Pp_next
        mean = af[t] + J * (h[t+1] - d - phi * af[t])
        var = Pf[t] * (1 - J * phi)
        h[t] = mean + np.sqrt(max(var, 1e-12)) * rng.standard_normal()
    return h


def _sample_sigma2(h, mu, phi, a0, b0, rng):
    """sigma_eta^2 | h : inverse-gamma conjugate."""
    T = len(h)
    e = h[1:] - mu - phi * (h[:-1] - mu)
    SS = np.sum(e**2) + (1 - phi**2) * (h[0] - mu)**2
    return 1.0 / rng.gamma(a0 + T / 2.0, 1.0 / (b0 + SS / 2.0))


def _sample_mu(h, phi, sig2, m0, s02, rng):
    """mu | h : Gaussian conjugate (h_1 stationary + AR(1) increments)."""
    T = len(h); c = 1 - phi
    prec = 1.0 / s02 + (1 - phi**2) / sig2 + (T - 1) * c**2 / sig2
    mean_num = m0 / s02 + (1 - phi**2) / sig2 * h[0] + c / sig2 * np.sum(h[1:] - phi * h[:-1])
    var = 1.0 / prec
    return var * mean_num + np.sqrt(var) * rng.standard_normal()


def _sample_phi(h, mu, sig2, phi_cur, a_beta, b_beta, rng):
    """phi | h : independence Metropolis. Proposal = the AR(1) regression Gaussian (which absorbs the
    likelihood exactly), so the accept ratio is just the stationary-h_1 factor x the Beta prior on (phi+1)/2."""
    x = h[:-1] - mu; y = h[1:] - mu
    Sxx = np.sum(x * x)
    phi_hat = np.sum(x * y) / Sxx
    prop = phi_hat + np.sqrt(sig2 / Sxx) * rng.standard_normal()
    if abs(prop) >= 1:
        return phi_cur

    def lp(p):
        return (0.5 * np.log(1 - p**2) - (1 - p**2) * (h[0] - mu)**2 / (2 * sig2)
                + (a_beta - 1) * np.log((1 + p) / 2) + (b_beta - 1) * np.log((1 - p) / 2))

    return prop if np.log(rng.random()) < lp(prop) - lp(phi_cur) else phi_cur


# ---- the sampler --------------------------------------------------------------------------------
def sv_gibbs(r, niter=6000, burn=1000, seed=0, offset=1e-4,
             mu0=0.0, s02=100.0, a_beta=20.0, b_beta=1.5, a_sig=5.0, b_sig=0.05,
             store_h_thin=4):
    """KSC (1998) mixture Gibbs sampler for the SV model. Returns posterior draws + latent-vol summary."""
    rng = np.random.default_rng(seed)
    r = np.asarray(r, float); r = r - r.mean()
    ystar = np.log(r**2 + offset)                                     # log-square transform (offset avoids -inf)
    T = len(r)
    q, mmean, v2 = mixture()

    mu, phi, sig2 = float(np.log(np.var(r))), 0.95, 0.05              # init
    h = np.full(T, mu)
    keep = niter - burn
    D = dict(mu=np.empty(keep), phi=np.empty(keep), sigma=np.empty(keep))
    hsum = np.zeros(T); h2sum = np.zeros(T); hpaths = []
    k = 0
    for it in range(niter):
        s = _sample_indicators(ystar, h, rng)
        h = _ffbs(ystar, mmean[s], v2[s], mu, phi, sig2, rng)
        mu = _sample_mu(h, phi, sig2, mu0, s02, rng)
        phi = _sample_phi(h, mu, sig2, phi, a_beta, b_beta, rng)
        sig2 = _sample_sigma2(h, mu, phi, a_sig, b_sig, rng)
        if it >= burn:
            D['mu'][k] = mu; D['phi'][k] = phi; D['sigma'][k] = np.sqrt(sig2)
            hsum += h; h2sum += h**2
            if k % store_h_thin == 0:
                hpaths.append(h.copy())
            k += 1
    D['h_mean'] = hsum / keep
    D['h_sd'] = np.sqrt(np.maximum(h2sum / keep - (hsum / keep)**2, 0.0))
    D['h_paths'] = np.array(hpaths)                                   # (thinned, T) for credible bands
    D['ystar'] = ystar; D['r'] = r
    return D