C07 · Stochastic volatility and VaR — Bayesian Modelling in Practice

The brief

You have joined the market-risk desk of a small asset manager. Every evening the desk publishes a one-day 99% Value-at-Risk (VaR) for its S&P 500 book: the loss that tomorrow's return should exceed on only 1 day in 100. Capital is held against that number, and the regulator back-tests it the simple way - by counting the days on which the loss was worse than the VaR.

The desk's current model is one line: daily returns are Normal with a constant mean and a constant volatility, estimated from history. The head of risk suspects it is wrong in two ways at once - "too relaxed the week everything falls apart, and far too scared for years afterwards" - and asks you for a replacement, plus the evidence that it is better.

You will build the model on data up to 30 June 2011 and be judged on the 2.5 years that follow, which the model never sees during fitting.

How this notebook works

Each task states what to deliver, not how. Write your code in the YOUR CODE HERE cells.
Stuck? h.hint("task3") reveals hints one level at a time: nudge → approach → code skeleton. Try to get by on nudges.
h.check("task3", vol_of_vol=...) compares your numbers with the reference solution.
A full worked solution lives in notebooks/solutions/. Open it only when you are done (or truly stuck).

import logging

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm
import pytensor.tensor as pt
from scipy import stats

from pymc_challenges import Hints, data

# PyTensor 3.3 logs a long but harmless "Rewrite failure" traceback when it compiles
# pm.GaussianRandomWalk with more than 127 steps. The model is unaffected; keep the log quiet.
logging.getLogger("pytensor.graph.rewriting.basic").setLevel(logging.CRITICAL)

RANDOM_SEED = 2008
rng = np.random.default_rng(RANDOM_SEED)
az.style.use("arviz-variat")

h = Hints("C07")
h.tasks()

C07 - Stochastic volatility: a Value-at-Risk number you can defend

task0 - What do returns look like? (3 hints, 2 check(s))
task1 - The desk's model and checks that expose it (3 hints, 2 check(s))
task2 - Fat tails (3 hints, 1 check(s))
task3 - Latent log-volatility and a sampler in trouble (3 hints, 2 check(s))
task4 - Does the volatility path make sense? (3 hints, 2 check(s))
task5 - Tomorrow's VaR (3 hints, 2 check(s))
task6 - Back-test and recommend (3 hints, 2 check(s))

The data

change is the daily log-return of the index. Returns are easier to read (and to put priors on) in percent, so that is the unit used everywhere below, including in the checks: ret = 100 * change.

The file runs to 2019. We use only the first 5.5 years: a latent volatility per trading day makes the model grow with the data, and 800 training days keep each fit to about a minute. train is all you may fit on; test is for Tasks 5-6 only.

data.describe("sp500")
sp500 = data.load("sp500")
sp500["ret"] = 100 * sp500["change"]

train = sp500.loc[:"2011-06-30"]
test = sp500.loc["2011-07-01":"2013-12-31"]
print(f"train: {train.index[0].date()} to {train.index[-1].date()}, {len(train)} days")
print(f"test:  {test.index[0].date()} to {test.index[-1].date()}, {len(test)} days")
train.head()

sp500: Daily close and log-return ('change').
  source: S&P 500 index daily closes (2008-2019), via pymc-examples.
  url:    https://raw.githubusercontent.com/pymc-devs/pymc-examples/main/examples/data/SP500.csv
train: 2008-05-02 to 2011-06-30, 798 days
test:  2011-07-01 to 2013-12-31, 629 days

	Close	change	ret
Date
2008-05-02	1413.900024	0.003230	0.323038
2008-05-05	1407.489990	-0.004544	-0.454389
2008-05-06	1418.260010	0.007623	0.762281
2008-05-07	1392.569946	-0.018280	-1.827985
2008-05-08	1397.680054	0.003663	0.366284

Task 0 · What do returns look like?

Deliver (on train) the evidence for or against each half of the desk's assumption:

Normal: compare the distribution of returns with a Normal of the same mean and standard deviation, in a way that shows the tails. Report the excess kurtosis.
Constant volatility: plot the autocorrelation function of the returns and of their absolute values, side by side. What does each say?

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

A histogram with a Normal curve on top looks fine for almost any return series, because the eye is drawn to the middle. Which plot puts the tails on the diagonal? And for the second half: if volatility were constant, would knowing that yesterday's move was large tell you anything about the size of today's?

2 Approach

Tails: a Normal Q-Q plot (scipy.stats.probplot(x, dist="norm", plot=ax)) and the excess kurtosis, scipy.stats.kurtosis(x) (Fisher definition: a Normal scores 0, not 3).

Clustering: the sample autocorrelation of ret and of np.abs(ret) for lags 1-30. For white noise, about 95% of the bars fall within ±2 / sqrt(T).

3 Skeleton

def acf(x, max_lag):
    x = x - x.mean()
    return np.array([np.dot(x[:-k], x[k:]) for k in range(1, max_lag + 1)]) / ___

r = train["ret"].values
stats.probplot(r, dist="norm", plot=ax1)
ax2.bar(lags, acf(r, 30)); ax2.bar(lags, acf(___, 30))
ax2.axhspan(-2 / np.sqrt(len(r)), 2 / np.sqrt(len(r)), alpha=0.1)
stats.kurtosis(r)

Self-checks in the notebook (h.check("task0", …)):

excess_kurtosis — excess (Fisher) kurtosis of the training returns, about 6.7 (a Normal has 0 - did you subtract 3?)
acf_abs_lag1 — lag-1 autocorrelation of the absolute training returns, about 0.25

Task 1 · The desk's model, and checks that expose it

Fit the desk's model: $r_t \sim \text{Normal}(\mu, \sigma)$, both constant.

Deliver

The fit, with priors that make sense for daily returns in percent.
A posterior predictive check built on test statistics you choose yourself - at least three, each aimed at something a VaR number depends on. For each one: the distribution of the statistic over replicated datasets, the observed value, and the posterior predictive p-value $P(T(y^\text{rep}) \ge T(y))$.
One sentence per statistic: what does the model get wrong, and would a density overlay (az.plot_ppc_dist) have shown it?

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

A posterior predictive check compares a statistic of the real data with the same statistic computed on each replicated dataset. The model fits mean and variance by construction, so those statistics can never fail. What does a VaR number depend on that the Normal does not get for free? Think tails, think worst day, think time.

2 Approach

After pm.sample_posterior_predictive(idata, extend_inferencedata=True), az.extract(idata, group="posterior_predictive", var_names="ret") gives an array of replicated histories. Apply a function along the date axis to get one value of T per replication; the posterior predictive p-value is (T_rep >= T_obs).mean(). Values near 0 or 1 mean the model cannot reproduce that feature.

Good candidates: scipy.stats.kurtosis, np.abs(x).max(), and the autocorrelation of np.abs(x) (one lag, or the mean over the first ten). az.plot_ppc_tstat(idata, t_stat=...) accepts "max", a quantile or your own function if you prefer ArviZ's plot.

3 Skeleton

with pm.Model(coords={"date": train.index}) as normal_model:
    mu = pm.Normal("mu", 0, ___)
    sigma = pm.HalfNormal("sigma", ___)              # percent per day
    pm.Normal("ret", mu, sigma, observed=r, dims="date")
    idata = pm.sample(random_seed=RANDOM_SEED)
    pm.sample_posterior_predictive(idata, extend_inferencedata=True, random_seed=RANDOM_SEED)

rep = az.extract(idata, group="posterior_predictive", var_names="ret")
rep = rep.transpose("sample", "date").values          # (replications, days)
t_rep = stats.kurtosis(rep, axis=1)
ppp = (t_rep >= ___).mean()

Self-checks in the notebook (h.check("task1", …)):

sigma — posterior mean of the constant daily volatility in PERCENT, about 1.83
ppp_kurtosis — posterior predictive p-value P(T_rep >= T_obs) for excess kurtosis under the Normal model - essentially 0

Task 2 · Fat tails

The textbook fix for fat tails is a Student-t likelihood. Keep $\mu$ and the scale constant.

Deliver the fit, the same predictive checks as in Task 1, and a verdict: which failures are fixed, which are not, and is any statistic now wrong in the other direction? Explain what the posterior of the degrees of freedom $\nu$ is telling you.

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

Change one thing - the likelihood - and rerun exactly the same checks. Watch for a statistic whose p-value flips from 0 to nearly 1: that is also a failure. Then ask what a single scale parameter has to do to cover both October 2008 and a quiet week in 2010.

2 Approach

pm.StudentT("ret", nu=nu, mu=mu, sigma=sigma, observed=...) with a weak prior such as pm.Gamma("nu", 2, 0.1). For nu <= 2 the variance is infinite and for nu <= 4 the kurtosis is, so replicated kurtosis and maxima are wildly skewed - plot them on a log axis. An independent heavy-tailed draw per day still has no memory.

3 Skeleton

with pm.Model(coords={"date": train.index}) as student_model:
    mu = pm.Normal("mu", 0, 0.5)
    sigma = pm.HalfNormal("sigma", 3)
    nu = pm.Gamma("nu", ___, ___)
    pm.StudentT("ret", nu=nu, mu=mu, sigma=sigma, observed=r, dims="date")
    ...
# same T_STATS loop as in Task 1; use ax.set_xscale("log") for kurtosis and max

Self-checks in the notebook (h.check("task2", …)):

nu — posterior mean of the degrees of freedom with constant scale, about 2.1 (extremely heavy tails)

Task 3 · Let volatility move

Model the log-volatility as a latent random walk: each day has its own volatility, tomorrow's is today's plus a small shock, and returns are Student-t around $\mu$ with that day's scale. That is one latent variable per trading day - about 800 parameters.

Deliver

The model, with a prior on the random walk's innovation scale (the "vol of vol") that you can defend: how much can volatility plausibly change in one day?
Your first fit and an honest reading of its diagnostics. Do not stop at "no divergences": check r_hat and ESS of every parameter, and find which one is in trouble.
A diagnosis (why is this posterior hard to sample?), a fix, and a final fit with clean diagnostics. If a fix trades one problem for another, show that too, and show that your final answer did not change with the sampler settings.

Budget: a fit takes about a minute. If you need more than three or four, think rather than wait.

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

Write the model as in the maths first and read the whole diagnostics table. A model can have zero divergences and still be unusable: which single parameter has an ESS in the tens? It is the one that ties 800 latent steps together. You have seen this geometry before, in hierarchical models - what was the cure there, and what is the time-series equivalent of "offset times scale"?

2 Approach

Centred: pm.GaussianRandomWalk("log_vol", sigma=vol_of_vol, init_dist=pm.Normal.dist(0, 1), dims="date"). Given the path, vol_of_vol is pinned to the size of its steps, so the sampler can only move it together with all the steps - slowly.

Non-centred: z ~ Normal(0, 1) per day and log_vol = log_vol_0 + pt.cumsum(vol_of_vol * z). This fixes the mixing but may produce divergences, because through the cumulative sum every z now influences every later day. Find out where they occur (compare vol_of_vol for divergent and other draws). If they sit in a region of higher curvature rather than a funnel neck, a smaller step size (target_accept=0.95) removes them - verify that they are gone and that the posterior did not move. nuts_sampler="nutpie" is fast for models of this size.

The prior on vol_of_vol is on the log scale, i.e. a relative daily change: 0.1 means "about 10% a day". A scale of several units would let volatility explain every return on its own.

3 Skeleton

with pm.Model(coords={"date": train.index}) as sv_model:
    vol_of_vol = pm.HalfNormal("vol_of_vol", ___)
    nu = pm.Gamma("nu", 2, 0.1)
    mu = pm.Normal("mu", 0, 0.5)
    log_vol_0 = pm.Normal("log_vol_0", 0, 1)
    z = pm.Normal("z", 0, 1, dims="date")
    log_vol = pm.Deterministic("log_vol", log_vol_0 + pt.cumsum(___), dims="date")
    pm.StudentT("ret", nu=nu, mu=mu, sigma=pm.math.exp(___), observed=r, dims="date")
    idata = pm.sample(target_accept=___, random_seed=RANDOM_SEED)

div = idata.sample_stats["diverging"].values.ravel()
vov = idata.posterior["vol_of_vol"].values.ravel()
np.median(vov[div]), np.median(vov)
az.summary(idata, var_names=["vol_of_vol", "nu", "mu"])

Self-checks in the notebook (h.check("task3", …)):

vol_of_vol — posterior mean of the daily innovation sd of log-volatility in the FINAL model, about 0.085
nu — posterior mean of the Student-t degrees of freedom once volatility moves, in the teens (not 2!)

Task 4 · Does the volatility path make sense?

Deliver

A plot of the posterior volatility path $e^{v_t}$ (mean and a credible band) over the training period, on top of the absolute returns, with the dates you would expect to matter marked - at least the Lehman bankruptcy (15 Sep 2008) and the "flash crash" (6 May 2010). Report the posterior mean of $e^{v_t}$ on its peak day and on the last training day.
The predictive checks of Task 1 for this model. Then think about what that check does and does not prove here, given that the replications reuse the fitted path, and add a check that does not have this weakness.

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

log_vol has a date dimension, so the path and its band are one reduction over ("chain", "draw") away. For the second part: if the replications are generated from the fitted path, which features of the data are they almost guaranteed to reproduce? What should be true of the returns after dividing by the fitted volatility?

2 Approach

vol = np.exp(idata.posterior["log_vol"]), then .mean(("chain", "draw")) and .quantile([0.03, 0.97], dim=("chain", "draw")). Plot over np.abs(ret) bars and add ax.axvline(pd.Timestamp("2008-09-15")) for the events.

The stronger check: standardised residuals (r - mu) / exp(log_vol) should be iid Student-t with nu degrees of freedom. Look at the ACF of their absolute values, and at the histogram of scipy.stats.t.cdf(resid, nu) (the PIT), which should be flat.

3 Skeleton

vol = np.exp(sv_idata.posterior["log_vol"])
vol_mean = vol.mean(("chain", "draw")).to_series()
lo, hi = vol.quantile([0.03, 0.97], dim=("chain", "draw")).values
peak_vol, last_vol = vol_mean.max(), vol_mean.iloc[___]

post = az.extract(sv_idata, var_names=["log_vol", "mu", "nu"], num_samples=500)
resid = (r[:, None] - post["mu"].values) / np.exp(post["log_vol"].transpose("date", "sample").values)
acf(np.abs(resid.mean(axis=1)), 30)
pit = stats.t.cdf(resid, ___).mean(axis=1)

Self-checks in the notebook (h.check("task4", …)):

peak_vol — posterior mean of exp(log_vol) on its highest day (October 2008), in percent - about 4.7
last_vol — posterior mean of exp(log_vol) on the last training day (30 June 2011), in percent - about 0.95

Task 5 · Tomorrow's number

It is the evening of 30 June 2011, the last training day.

Deliver the one-day-ahead 99% VaR for 1 July 2011 - the 1% quantile of the posterior predictive distribution of tomorrow's return, in percent (a negative number) - under the desk's model and under yours. Tomorrow's volatility is not in your posterior: get it there. Propagate all the uncertainty (parameters, today's volatility, tomorrow's shock).

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

VaR is a quantile of the posterior predictive distribution of a return that has not happened yet - not a quantile of a parameter, and not mu - 2.33 * sigma at the posterior mean. For your model, tomorrow's volatility is one more step of the same random walk. What do you need, per posterior draw, to take that step?

2 Approach

Per posterior draw: log_vol_next = log_vol[last day] + vol_of_vol * Normal(0, 1), then ret_next = mu + exp(log_vol_next) * StudentT(nu). Use NumPy on the extracted draws (rng.standard_normal, rng.standard_t). A 1% quantile from 4000 numbers is noisy, so simulate several futures per draw. The VaR is np.quantile(ret_next, 0.01) over everything. Do the same for the Normal model with a new rng.normal(mu, sigma) per draw.

3 Skeleton

p = az.extract(sv_idata, var_names=["log_vol", "vol_of_vol", "nu", "mu"])
log_vol_T = p["log_vol"].isel(date=-1).values                 # (draws,)
size = (50, log_vol_T.size)
log_vol_next = log_vol_T + p["vol_of_vol"].values * rng.standard_normal(size)
ret_next = p["mu"].values + np.exp(___) * rng.standard_t(p["nu"].values, size=size)
var_sv = np.quantile(ret_next, ___)

Self-checks in the notebook (h.check("task5", …)):

var_desk — 1% quantile of the constant-volatility Normal posterior predictive, in percent - about -4.3
var_sv — 1% predictive quantile for 1 July 2011 under the stochastic-volatility model, in percent - about -2.7 (did you add tomorrow's volatility shock and the Student-t noise?)

Now the verdict. For every day in test you need the VaR that would have been published the evening before, using only returns up to that evening.

Deliver

The daily VaR series of your model and of the desk's model over the test period, in one plot with the realised returns, exceedances marked.
The number of exceedances of each against the expected 1%, when they happened, and the average VaR (a proxy for the capital the desk must hold).
A recommendation to the head of risk in a short paragraph - including what your model still gets wrong.

Refitting an 800-parameter model 629 times is not an option, and fitting once on train + test would let each day's volatility estimate peek at the future. Find a way to update your belief about today's volatility day by day that is cheap and does not cheat, and state clearly which shortcut you are taking.

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

The expensive part of your model is learning the parameters; the part that must change every day is the volatility state. If you freeze the parameters at their training posterior, updating the state for one new return is a one-dimensional application of Bayes' rule: prior = yesterday's belief pushed through one random-walk step, likelihood = Student-t density of the new return. How can you do that with samples?

2 Approach

A bootstrap particle filter. Carry a cloud of log-volatility "particles" (several per retained posterior draw, each row with its own frozen vol_of_vol, nu, mu, started at that draw's last log_vol). Every day: (1) add a random-walk step to every particle and simulate a return from each - their 1% quantile is the VaR published the evening before; (2) when the real return arrives, weight particles by scipy.stats.t.logpdf(r_t, nu, mu, exp(v)) and resample them in proportion to the weights. Nothing from the future is used.

Then count ret < VaR days for each model, look at when they happen (clustered exceedances are a failure even if the count is fine), and compare the average VaR. With k exceedances in n days, Beta(1 + k, 1 + n - k) is the posterior of the exceedance rate under a flat prior.

3 Skeleton

v = np.repeat(log_vol_T[keep][:, None], M, axis=1)               # (draws, particles)
for t, r_t in enumerate(r_test):
    v = v + vol_of_vol[:, None] * rng.standard_normal(v.shape)    # predict
    r_rep = mu[:, None] + np.exp(v) * rng.standard_t(nu[:, None], size=v.shape)
    var[t] = np.quantile(r_rep, 0.01)
    logw = stats.t.logpdf(r_t, nu[:, None], mu[:, None], ___)     # update
    w = np.exp(logw - logw.max(axis=1, keepdims=True))
    w /= w.sum(axis=1, keepdims=True)
    idx = np.stack([rng.choice(M, size=M, p=w_row) for w_row in w])   # resample per row
    v = np.take_along_axis(v, idx, axis=1)

exceed = r_test < var
stats.beta.ppf([0.03, 0.97], 1 + exceed.sum(), 1 + len(r_test) - exceed.sum())

Self-checks in the notebook (h.check("task6", …)):

exceed_desk — number of test days on which the return was below the desk's constant VaR - a handful, all within two weeks
exceed_sv — number of test days below the stochastic-volatility VaR (expected at 1% - about 6; a filter that works lands between 6 and 14)

Going further

Leverage. Let the shock to tomorrow's volatility be negatively correlated with today's return, or give the returns a skewed distribution. Does the lower-tail coverage in the back-test improve, and what happens to the upper tail?
Mean reversion. Replace the random walk by an AR(1) log-volatility, $v_t = \bar v + \phi\,(v_{t-1} - \bar v) + \sigma_v \epsilon_t$. How well is $\phi$ identified, and does it change the 10-day-ahead VaR more than the 1-day-ahead one?
Check the shortcut. Refit the model on all data up to the end of 2012 and compare its volatility on the last day, and its parameters, with what the particle filter carried forward from June 2011.
Expected shortfall. Regulators have moved from VaR to the expected loss beyond the VaR. You have the predictive draws: compute it, and think about how you would back-test it.

h.progress()

C07: 0/21 hints used, 13/13 checks passed.

Reference solution

A complete, executed solution with commentary on every task — including the ways the obvious model fails. Try the hint ladder first.

Open the worked solution →

#The brief

#How this notebook works

#The data

#Task 0 · What do returns look like?

#Task 1 · The desk's model, and checks that expose it

#Task 2 · Fat tails

#Task 3 · Let volatility move

#Task 4 · Does the volatility path make sense?

#Task 5 · Tomorrow's number

#Task 6 · Back-test, and recommend

#Going further

Reference solution

The brief

How this notebook works

The data

Task 0 · What do returns look like?

Task 1 · The desk's model, and checks that expose it

Task 2 · Fat tails

Task 3 · Let volatility move

Task 4 · Does the volatility path make sense?

Task 5 · Tomorrow's number

Task 6 · Back-test, and recommend

Going further