Invariance & Horizon Projection — the PyMC engine¶

The from-scratch notebook projected the invariant to the horizon as if its distribution were known exactly. It is not — we estimate it from finite data, and that estimation error propagates all the way to the horizon risk number. This notebook does the projection Bayesianly: it fits the invariant's distribution in PyMC, then carries the full parameter uncertainty into a posterior distribution of the horizon Value-at-Risk. The punchline connects to Project 2 — parameter uncertainty adds to risk, and the honest horizon VaR is fatter than the plug-in.

Self-contained; no prior reading required.

The plan¶

  1. Fit a Student-$t$ distribution to the daily invariant (SPY log-returns) — the posterior over location, scale, and the tail-fatness $\nu$.
  2. A Bayesian invariance check: is the lag-1 autocorrelation credibly zero?
  3. Bayesian horizon projection: propagate the posterior into the horizon distribution, giving the horizon VaR with error bars, and show the Bayesian-predictive VaR is more conservative than the plug-in.
In [1]:
import numpy as np, pandas as pd, matplotlib.pyplot as plt
import pymc as pm, arviz as az
import horizon as hz

plt.rcParams.update({"figure.figsize": (9,4.5), "axes.grid": True, "grid.alpha": .25,
                     "axes.spines.top": False, "axes.spines.right": False, "font.size": 11})
BLUE, ORANGE, GREEN, RED, GREY, PURP = "#2b6cb0","#dd6b20","#2f855a","#c53030","#718096","#6b46c1"
RNG = 3

D = pd.read_csv("market_levels_daily.csv", index_col=0, parse_dates=True)
ret = np.diff(np.log(D["SPY"].values)) * 100        # the daily invariant (log-return, %)
print("Daily invariant: %d obs, mean %.3f%%, std %.3f%%, excess-kurtosis %.1f"
      % (len(ret), ret.mean(), ret.std(), hz.sample_moments(ret)["exkurt"]))

g++ not available, if using conda: `conda install gxx`

Daily invariant: 3770 obs, mean 0.051%, std 1.078%, excess-kurtosis 11.5

The dataset¶

Daily SPY log-returns (%), 2010–2024 (3,770 observations) — the equity invariant identified in the from-scratch notebook (autocorrelation $\approx 0$, distribution stable across sub-samples). It is fat-tailed, which is why we model it with a Student-$t$ rather than a normal.

1. Bayesian fit of the invariant¶

We model the daily invariant as Student-$t$: $r_t\sim t(\nu,\mu,\sigma)$. The degrees of freedom $\nu$ control the tail fatness — and, crucially, $\nu$ is poorly identified from data (the tail is where observations are scarce), so its posterior is wide. That uncertainty is the whole point: it is what makes the horizon tail risk uncertain.

In [2]:
with pm.Model() as tfit:
    mu    = pm.Normal("mu", 0.0, 0.5)
    sigma = pm.HalfNormal("sigma", 3.0)
    nu    = pm.Gamma("nu", alpha=2.0, beta=0.2)          # dof > 0, mean ~10
    pm.StudentT("r", nu=nu, mu=mu, sigma=sigma, observed=ret)
    idata = pm.sample(1500, tune=1500, chains=4, cores=1, target_accept=0.9,
                      progressbar=False, random_seed=RNG)
print("max R-hat:", float(az.summary(idata, var_names=["mu","sigma","nu"])["r_hat"].max()))
print(az.summary(idata, var_names=["mu","sigma","nu"]).iloc[:, :5].to_string())

Initializing NUTS using jitter+adapt_diag...

Sequential sampling (4 chains in 1 job)

NUTS: [mu, sigma, nu]

Sampling 4 chains for 1_500 tune and 1_500 draw iterations (6_000 + 6_000 draws total) took 4 seconds.

max R-hat: 1.0
         mean      sd eti89_lb eti89_ub  ess_bulk
mu     0.0935  0.0128    0.073     0.11      3840
sigma  0.6378  0.0149     0.61     0.66      3572
nu      2.761   0.153      2.5        3      3498
In [3]:
fig, ax = plt.subplots(1, 2, figsize=(12, 4.2))
nud = idata.posterior["nu"].values.flatten()
ax[0].hist(nud, bins=50, density=True, color=BLUE, alpha=.8)
ax[0].axvline(nud.mean(), color=RED, lw=2, label="mean %.1f" % nud.mean())
ax[0].set_xlabel("degrees of freedom $\\nu$ (tail fatness)"); ax[0].set_title("Posterior of the tail parameter"); ax[0].legend()
# posterior predictive daily density vs data
mud = idata.posterior["mu"].values.flatten(); sigd = idata.posterior["sigma"].values.flatten()
from scipy.stats import t as tdist
xx = np.linspace(-8, 8, 400)
ax[1].hist(ret, bins=120, density=True, color=GREY, alpha=.5, label="data")
pp = np.mean([tdist.pdf(xx, nud[k], mud[k], sigd[k]) for k in range(0, len(nud), 20)], axis=0)
ax[1].plot(xx, pp, color=BLUE, lw=2, label="Student-t posterior fit")
ax[1].set_xlim(-8,8); ax[1].set_yscale("log"); ax[1].set_xlabel("daily return (%)"); ax[1].set_title("Fit captures the fat tails (log scale)"); ax[1].legend()
plt.tight_layout(); plt.show()
print("The posterior of nu is wide (%.1f-%.1f): the data cannot pin down the tail thickness, so any" % (np.percentile(nud,3), np.percentile(nud,97)))
print("risk number that depends on the tail -- like VaR -- must inherit that uncertainty.")
No description has been provided for this image 
The posterior of nu is wide (2.5-3.1): the data cannot pin down the tail thickness, so any
risk number that depends on the tail -- like VaR -- must inherit that uncertainty.

2. A Bayesian invariance check¶

Is the invariant really serially independent? We fit an AR(1) coefficient $\phi$ to it; if the series is i.i.d., the posterior of $\phi$ should sit tightly around zero. (A wandering level would give $\phi\approx1$.)

In [4]:
with pm.Model() as ar1:
    phi = pm.Normal("phi", 0.0, 0.3)
    c   = pm.Normal("c", 0.0, 0.5)
    s   = pm.HalfNormal("s", 3.0)
    pm.Normal("r", c + phi*ret[:-1], s, observed=ret[1:])
    idata_ar = pm.sample(1000, tune=1000, chains=4, cores=1, progressbar=False, random_seed=RNG)
phid = idata_ar.posterior["phi"].values.flatten()
print("Posterior of AR(1) coefficient phi: mean %.3f, 94%% HDI [%.3f, %.3f]"
      % (phid.mean(), np.percentile(phid,3), np.percentile(phid,97)))
print("phi sits close to zero -> the invariant is (approximately) serially independent, as required.")
print("Contrast: the same model on the log-PRICE level would return phi ~ 1 (a persistent random walk).")

Initializing NUTS using jitter+adapt_diag...

Sequential sampling (4 chains in 1 job)

NUTS: [phi, c, s]

C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\pymc\step_methods\hmc\quadpotential.py:321: RuntimeWarning: overflow encountered in dot
  return 0.5 * np.dot(x, v_out)

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 1 seconds.

Posterior of AR(1) coefficient phi: mean -0.098, 94% HDI [-0.128, -0.068]
phi sits close to zero -> the invariant is (approximately) serially independent, as required.
Contrast: the same model on the log-PRICE level would return phi ~ 1 (a persistent random walk).

3. Bayesian horizon projection: VaR with error bars¶

Now the projection. For each posterior draw $(\mu,\sigma,\nu)$ the invariant is a specific Student-$t$; because it is i.i.d., we simulate the $T$-day horizon return as the sum of $T$ draws and read off its 1% Value-at-Risk. Doing this for every posterior draw gives the posterior distribution of the horizon VaR. Pooling all the draws together yields the Bayesian-predictive horizon distribution, which integrates the parameter uncertainty — and is therefore fatter-tailed than the plug-in projection that fixes the parameters at their posterior mean.

In [5]:
rng = np.random.default_rng(RNG)
K = 300; nsim = 20000
idx = rng.integers(0, len(nud), K)
mu_k, sig_k, nu_k = mud[idx], sigd[idx], nud[idx]

def horizon(T, alpha=0.01):
    var_draws = np.empty(K); pooled = []
    for k in range(K):
        s = (mu_k[k] + sig_k[k]*rng.standard_t(nu_k[k], size=(nsim, T))).sum(1)
        var_draws[k] = np.percentile(s, 100*alpha); pooled.append(s)
    pooled = np.concatenate(pooled)
    # plug-in: posterior-mean parameters
    sp = (mud.mean() + sigd.mean()*rng.standard_t(nud.mean(), size=(K*nsim//4, T))).sum(1)
    return var_draws, np.percentile(pooled, 100*alpha), np.percentile(sp, 100*alpha)

T = 21
var_post, var_pred, var_plug = horizon(T)
plt.hist(var_post, bins=40, density=True, color=BLUE, alpha=.75, label="posterior of 1-month VaR")
plt.axvline(var_plug, color=GREY, lw=2, ls="--", label="plug-in VaR %.2f%%" % var_plug)
plt.axvline(var_pred, color=RED, lw=2, label="Bayesian-predictive VaR %.2f%%" % var_pred)
plt.xlabel("1-month 1% VaR (cumulative %)"); plt.ylabel("density")
plt.title("The 1-month horizon VaR is itself uncertain"); plt.legend(); plt.show()
print("Posterior of the 1-month VaR: mean %.2f%%, 94%% HDI [%.2f%%, %.2f%%] -- a full percentage point"
      % (var_post.mean(), np.percentile(var_post,3), np.percentile(var_post,97)))
print("of uncertainty in the risk number itself. The Bayesian-PREDICTIVE VaR (%.2f%%) is more extreme" % var_pred)
print("than the plug-in (%.2f%%): integrating parameter uncertainty fattens the tail -- estimation risk," % var_plug)
print("now at the horizon. Reporting a single plug-in VaR quietly overstates how much we actually know.")
No description has been provided for this image 
Posterior of the 1-month VaR: mean -11.47%, 94% HDI [-13.14%, -10.05%] -- a full percentage point
of uncertainty in the risk number itself. The Bayesian-PREDICTIVE VaR (-11.50%) is more extreme
than the plug-in (-11.43%): integrating parameter uncertainty fattens the tail -- estimation risk,
now at the horizon. Reporting a single plug-in VaR quietly overstates how much we actually know.

4. Summary¶

  • A Bayesian Student-$t$ fit of the invariant returns a wide posterior for the tail parameter $\nu$ — the data genuinely cannot pin down how fat the tail is.
  • A one-line AR(1) check confirms the invariant is serially independent ($\phi\approx0$), the defining property from step one.
  • Propagating the posterior through the projection makes the horizon VaR a distribution, not a number — here a full percentage point of uncertainty at the one-month horizon — and the Bayesian-predictive VaR is more conservative than the plug-in, because integrating parameter uncertainty fattens the tail.

This is the estimation-risk lesson of Project 2, now attached to the horizon machinery of the from-scratch notebook: the honest horizon risk carries the uncertainty of the invariant that generates it.