Stochastic Volatility — R (stochvol, Kastner)¶
Two-engine validation of the NumPyro SV fit¶
Companion to sv_numpyro.ipynb. Kastner's stochvol is the gold-standard SV sampler: it uses the ASIS / ancillarity–sufficiency interweaving strategy (Kastner & Frühwirth-Schnatter 2014), which mixes far better on the SV funnel than a naïve Gibbs sampler. Same model, same S&P 500 data — the point is that a completely different algorithm lands on the same posterior for $\phi$ and $\sigma_\eta$.
$$r_t = e^{h_t/2}\,\varepsilon_t,\qquad h_t = \mu + \phi\,(h_{t-1}-\mu) + \sigma_\eta\,\eta_t.$$
suppressMessages(library(stochvol))
d <- read.csv("sp500ret.csv"); y <- d$ret - mean(d$ret); dts <- as.Date(d$date)
set.seed(1)
res <- svsample(y, draws = 8000, burnin = 2000, quiet = TRUE) # ASIS sampler; ~30s
sm <- summary(res, showlatent = FALSE)$para
cat("stochvol posterior (mu, phi, sigma):\n"); print(round(sm[, c("mean","5%","95%")], 4))
cat(sprintf("\nphi %.4f sigma %.4f mu %.4f\n", sm["phi","mean"], sm["sigma","mean"], sm["mu","mean"]))
cat("NumPyro (sv_numpyro.ipynb): phi 0.987 sigma 0.155 mu -0.233 -> two engines agree\n")
stochvol posterior (mu, phi, sigma):
mean 5% 95% mu -0.2314 -0.5016 0.0447 phi 0.9869 0.9816 0.9917 sigma 0.1524 0.1305 0.1749 exp(mu/2) 0.8938 0.7782 1.0226 sigma^2 0.0234 0.0170 0.0306
phi 0.9869 sigma 0.1524 mu -0.2314
NumPyro (sv_numpyro.ipynb): phi 0.987 sigma 0.155 mu -0.233 -> two engines agree
# SV volatility path with a 90% credible band (the stochvol posterior of exp(h/2))
options(repr.plot.width = 11, repr.plot.height = 4.8)
h <- latent(res) # draws x T (log-variance)
vq <- apply(exp(h/2), 2, quantile, probs = c(.05, .5, .95)) # 3 x T (daily % vol)
plot(dts, pmin(abs(y), 12), type = "h", col = "grey86", xlab = "", ylab = "daily volatility (%)",
ylim = c(0, 12), main = "stochvol (Kastner ASIS sampler): SV volatility with a 90% credible band")
polygon(c(dts, rev(dts)), c(vq[1, ], rev(vq[3, ])), col = rgb(.27, .51, .71, .35), border = NA)
lines(dts, vq[2, ], col = "steelblue", lwd = 1)
legend("topright", c("|return|", "SV vol (median)", "90% band"),
col = c("grey70", "steelblue", rgb(.27, .51, .71, .6)), lty = c(1, 1, NA), lwd = c(1, 1, NA),
pch = c(NA, NA, 15), bty = "n", cex = .8)
# Parameter posteriors with the NumPyro estimates overlaid (dashed) = the validation
options(repr.plot.width = 11, repr.plot.height = 4.2); par(mfrow = c(1, 2))
pd <- as.matrix(para(res))
hist(pd[, "phi"], breaks = 50, col = rgb(.27, .51, .71, .7), border = "white", freq = FALSE,
xlab = expression(phi), main = sprintf("persistence phi = %.3f", mean(pd[, "phi"])))
abline(v = 0.9865, col = "firebrick", lty = 2, lwd = 2); legend("topleft", "NumPyro 0.987", text.col = "firebrick", bty = "n", cex = .8)
hist(pd[, "sigma"], breaks = 50, col = rgb(.18, .55, .34, .7), border = "white", freq = FALSE,
xlab = expression(sigma[eta]), main = sprintf("vol-of-vol sigma = %.3f", mean(pd[, "sigma"])))
abline(v = 0.1546, col = "firebrick", lty = 2, lwd = 2); legend("topright", "NumPyro 0.155", text.col = "firebrick", bty = "n", cex = .8)
par(mfrow = c(1, 1))
Results — two engines agree¶
| parameter | NumPyro (NUTS) | stochvol (ASIS) |
|---|---|---|
| persistence $\phi$ | 0.987 | 0.987 |
| vol-of-vol $\sigma_\eta$ | 0.155 | 0.152 |
| mean log-var $\mu$ | −0.233 | −0.231 |
A Hamiltonian sampler (NumPyro) and a bespoke interweaving Gibbs sampler (stochvol) — different algorithms, different languages — recover the same SV posterior, so the fit is trustworthy. stochvol is markedly faster (~20–30 s) because its sampler is purpose-built for the SV funnel; NumPyro is the more flexible engine for the extensions to come. The volatility band tells the same story as the NumPyro notebook: latent vol carrying genuine uncertainty.
Hooks for the next notebooks: para(res) already exposes nu and rho columns — stochvol::svtsample (SV-t, fat-tailed $\varepsilon$) and svlsample (leverage SV, $\rho<0$) drop straight in, mirroring the NumPyro extensions.