Bayesian Estimation & Estimation Risk — the PyMC engine¶

Risk and Asset Allocation¶

The from-scratch notebook (niw_python.ipynb) built the Normal-Inverse-Wishart posterior by hand and showed that estimation error deceives the sample efficient frontier. This notebook does the same estimation in PyMC and then uses it for the one thing point-estimators simply cannot do:

propagate parameter uncertainty all the way into the portfolio — producing a distribution over the optimal weights and over the efficient frontier itself, not a single fragile answer.

That distribution is the honest picture a manager should see, and it is the direct motivation for robust allocation (Project 4).

Self-contained; no prior reading required.

Roadmap¶

  1. Fit $(\mu,\Sigma)$ with a PyMC model and cross-check the posterior mean against the closed-form NIW.
  2. The posterior predictive distribution of next week's returns.
  3. Allocation under uncertainty — turn the posterior over $\Sigma$ into a posterior over minimum-variance weights.
  4. The efficient frontier with error bars.
In [1]:
import numpy as np, pandas as pd, matplotlib.pyplot as plt
import pymc as pm, arviz as az
import niw

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 = "#2b6cb0","#dd6b20","#2f855a","#c53030","#718096"
RNG = 42
print("PyMC", pm.__version__, "| ArviZ", az.__version__)

R = pd.read_csv("sector_etf_weekly.csv", index_col=0, parse_dates=True)
SECTORS = [c for c in R.columns if c != "SPY"]
# use a 2-year window: short enough that parameter uncertainty is visible
Xw = R[SECTORS].values[-104:]
T, N = Xw.shape
print("Estimation window:", T, "weeks,", N, "sectors")

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

PyMC 6.0.1 | ArviZ 1.2.0
Estimation window: 104 weeks, 10 sectors

The real dataset in detail: US sector ETFs¶

Weekly log-returns (%) of the ten SPDR sector ETFs plus SPY, Jan 2019 – Dec 2024 (312 weeks), dividend/split-adjusted. Each ETF holds the S&P 500 members of one GICS sector — Materials XLB, Energy XLE, Financials XLF, Industrials XLI, Technology XLK, Staples XLP, Utilities XLU, Health-care XLV, Discretionary XLY, Communications XLC — with SPY as the market benchmark. Here we deliberately estimate on a short 2-year (104-week) window so that parameter uncertainty is large enough to see its effect on the portfolio.

1. The model in PyMC¶

We separate the covariance into volatilities and a correlation matrix (an LKJCholeskyCov prior, numerically friendlier than a raw Inverse-Wishart) and give the mean a weakly-informative normal prior: $$ \mu \sim \mathcal N(0,\ 3^2\mathbf I), \qquad \Sigma = LL' \ \text{with}\ L\sim \text{LKJCholeskyCov}(\eta=2), \qquad x_t \sim \mathcal N(\mu,\Sigma). $$ The LKJ shape $\eta=2$ mildly favours lower correlations (a gentle prior regularisation). NUTS samples the joint posterior of $(\mu,\Sigma)$; we then read off whatever function of them we need — including portfolio weights.

In [2]:
sd_guess = Xw.std(0).mean()
with pm.Model() as model:
    chol, corr, sds = pm.LKJCholeskyCov("L", n=N, eta=2.0,
                                        sd_dist=pm.HalfNormal.dist(sd_guess), compute_corr=True)
    mu = pm.Normal("mu", 0.0, 3.0, shape=N)
    Sigma = pm.Deterministic("Sigma", chol @ chol.T)
    pm.MvNormal("obs", mu=mu, chol=chol, observed=Xw)
    idata = pm.sample(800, tune=1000, 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"])["r_hat"].max()))

Initializing NUTS using jitter+adapt_diag...

Sequential sampling (4 chains in 1 job)

NUTS: [L, mu]

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

max R-hat: 1.0
In [3]:
# Cross-check the PyMC posterior mean against the closed-form NIW posterior
mu_post = idata.posterior["mu"].mean(("chain","draw")).values
Sig_post = idata.posterior["Sigma"].mean(("chain","draw")).values
post_niw = niw.niw_posterior(Xw, np.zeros(N), 5.0, np.mean(np.diag(niw.sample_moments(Xw)[1]))*np.eye(N), N+4)
mu_niw, Sig_niw = niw.post_mean_cov(post_niw)
mu_hat, S_hat = niw.sample_moments(Xw)

df = pd.DataFrame({"sample": mu_hat, "PyMC posterior": mu_post, "conjugate NIW": mu_niw}, index=SECTORS).round(3)
print(df)
print("\nPyMC and the closed-form NIW agree on the posterior mean (both shrink the sample mean).")

     sample  PyMC posterior  conjugate NIW
XLB   0.143           0.141          0.137
XLC   0.764           0.759          0.729
XLE   0.137           0.143          0.131
XLF   0.385           0.389          0.368
XLI   0.340           0.338          0.324
XLK   0.682           0.669          0.650
XLP   0.096           0.099          0.092
XLU   0.141           0.147          0.134
XLV   0.057           0.060          0.054
XLY   0.548           0.542          0.522

PyMC and the closed-form NIW agree on the posterior mean (both shrink the sample mean).

2. Posterior predictive: next week's returns¶

The posterior predictive draws parameters from the posterior and then a return from the model, so it folds in both market randomness and our uncertainty about $(\mu,\Sigma)$. Its spread is what an honest risk model should use.

In [4]:
with model:
    ppc = pm.sample_posterior_predictive(idata, var_names=["obs"], progressbar=False, random_seed=RNG)
xnew = ppc.posterior_predictive["obs"].values.reshape(-1, N)
pred_cov = np.cov(xnew.T)
infl = np.diag(pred_cov) / np.diag(S_hat)
print("Predictive vs plug-in variance (inflation from parameter + sampling uncertainty):")
for i,s in enumerate(SECTORS):
    print("  %-4s  plug-in %6.3f   predictive %6.3f   x%.3f" % (s, S_hat[i,i], pred_cov[i,i], infl[i]))

j = SECTORS.index("XLK")
plt.hist(Xw[:,j], bins=25, density=True, color=GREY, alpha=.6, label="historical %s" % SECTORS[j])
plt.hist(xnew[:,j], bins=60, density=True, color=BLUE, alpha=.5, label="posterior predictive")
plt.title("Posterior predictive vs historical returns: %s" % SECTORS[j]); plt.xlabel("weekly %"); plt.legend(); plt.show()

Sampling: [obs]

Predictive vs plug-in variance (inflation from parameter + sampling uncertainty):
  XLB   plug-in  4.647   predictive  4.069   x0.876
  XLC   plug-in  4.450   predictive  4.255   x0.956
  XLE   plug-in  7.272   predictive  7.423   x1.021
  XLF   plug-in  5.652   predictive  4.845   x0.857
  XLI   plug-in  3.854   predictive  3.213   x0.834
  XLK   plug-in  7.611   predictive  7.183   x0.944
  XLP   plug-in  1.738   predictive  1.719   x0.989
  XLU   plug-in  5.646   predictive  5.706   x1.011
  XLV   plug-in  2.719   predictive  2.622   x0.964
  XLY   plug-in  6.927   predictive  6.216   x0.897
No description has been provided for this image

3. Allocation under uncertainty — a distribution over weights¶

Here is the payoff of doing this in a Bayesian framework. The minimum-variance weights depend only on the covariance, $$w(\Sigma) = \frac{\Sigma^{-1}\mathbf 1}{\mathbf 1'\Sigma^{-1}\mathbf 1}.$$ A point estimator gives one weight vector. But we have a posterior over $\Sigma$, so we can push every posterior draw through the formula and obtain a posterior over the weights. The width of that posterior tells the manager how much to trust the allocation — information the plug-in approach throws away.

In [5]:
Sig_draws = idata.posterior["Sigma"].values.reshape(-1, N, N)
Sig_draws = Sig_draws[::3]                                  # thin for speed
W = np.array([niw.min_var_weights(S) for S in Sig_draws])   # posterior of min-var weights
w_plugin = niw.min_var_weights(S_hat)                       # single plug-in answer
w_postmean = niw.min_var_weights(Sig_post)

order = np.argsort(W.mean(0))
plt.figure(figsize=(11,5))
plt.boxplot([W[:,i] for i in order], vert=True, showfliers=False,
            patch_artist=True, boxprops=dict(facecolor="#bee3f8"), medianprops=dict(color=BLUE))
plt.plot(range(1,N+1), w_plugin[order], "D", color=RED, ms=7, label="plug-in (sample) weight")
plt.xticks(range(1,N+1), np.array(SECTORS)[order]); plt.axhline(0, color="k", lw=.6)
plt.ylabel("portfolio weight"); plt.title("Posterior distribution of minimum-variance weights (each sector)")
plt.legend(); plt.show()
print("Every sector's weight is a DISTRIBUTION, often straddling zero -> the data cannot actually")
print("tell us whether to be long or short. The plug-in (red) reports a single confident number the")
print("data does not support. This uncertainty is what robust allocation (Project 4) will exploit.")
No description has been provided for this image 
Every sector's weight is a DISTRIBUTION, often straddling zero -> the data cannot actually
tell us whether to be long or short. The plug-in (red) reports a single confident number the
data does not support. This uncertainty is what robust allocation (Project 4) will exploit.
In [6]:
# How different are the plug-in and posterior-mean allocations, and how uncertain?
tbl = pd.DataFrame({
    "plug-in": w_plugin,
    "posterior mean": W.mean(0),
    "posterior sd": W.std(0)}, index=SECTORS).round(3)
print(tbl)
print("\nTypical weight uncertainty (posterior sd): %.3f -- comparable to the weights themselves." % W.std(0).mean())

     plug-in  posterior mean  posterior sd
XLB   -0.136          -0.073         0.077
XLC    0.186           0.161         0.059
XLE    0.136           0.115         0.042
XLF   -0.239          -0.154         0.076
XLI    0.174           0.151         0.109
XLK    0.039           0.044         0.053
XLP    0.696           0.576         0.082
XLU   -0.001           0.020         0.051
XLV    0.163           0.158         0.078
XLY   -0.017           0.003         0.063

Typical weight uncertainty (posterior sd): 0.069 -- comparable to the weights themselves.

4. The efficient frontier with error bars¶

The same logic applies to the whole efficient frontier. For each posterior draw of $(\mu,\Sigma)$ we trace a frontier; the collection gives a band of plausible frontiers. Compare that band with the single plug-in frontier — and remember from the from-scratch notebook that even the plug-in frontier is optimistic about true performance.

In [7]:
mu_draws = idata.posterior["mu"].values.reshape(-1, N)[::3]
grid = np.linspace(-0.1, 0.9, 30)
curves = []
for md_, Sd_ in zip(mu_draws, Sig_draws):
    try:
        vols = [np.sqrt(w @ Sd_ @ w) for w in (niw.frontier_weights(md_, Sd_, m) for m in grid)]
        curves.append(vols)
    except np.linalg.LinAlgError:
        continue
curves = np.array(curves)
lo, med, hi = np.percentile(curves, [5,50,95], axis=0)
plug_vols = [np.sqrt(w @ S_hat @ w) for w in (niw.frontier_weights(mu_hat, S_hat, m) for m in grid)]

plt.figure(figsize=(8.5,5.5))
plt.fill_betweenx(grid, lo, hi, color="#bee3f8", alpha=.7, label="posterior 5–95% band")
plt.plot(med, grid, color=BLUE, lw=2, label="posterior median frontier")
plt.plot(plug_vols, grid, "--", color=RED, lw=2, label="plug-in (sample) frontier")
plt.xlabel("volatility (weekly %)"); plt.ylabel("expected return (weekly %)")
plt.title("The efficient frontier is uncertain — a band, not a line")
plt.legend(); plt.show()
print("The plug-in frontier is a single optimistic line; the posterior reveals a wide band of")
print("frontiers consistent with the data. Committing to the plug-in ignores that uncertainty.")
No description has been provided for this image 
The plug-in frontier is a single optimistic line; the posterior reveals a wide band of
frontiers consistent with the data. Committing to the plug-in ignores that uncertainty.

4b. At scale: allocation uncertainty explodes with dimension¶

The distribution over weights is the whole point of the Bayesian approach, and it matters far more at high dimension, where the data must pin down 48 assets' worth of covariance from a short history. NUTS on a 48-asset LKJ is impractical, so we draw $(\mu,\Sigma)$ from the equivalent conjugate Normal-Inverse-Wishart posterior and push each draw through the minimum-variance formula. The weight posteriors are so wide that most positions cannot even be signed.

In [8]:
from scipy.stats import invwishart
def niw_weight_draws(Z, ndraw=800, seed=1):
    T, Nn = Z.shape; muh = Z.mean(0); Sc = (Z-muh).T @ (Z-muh)
    nu0 = Nn + 2; Psi0 = np.mean(np.diag(np.cov(Z.T))) * np.eye(Nn)
    Sig_d = invwishart.rvs(df=nu0+T, scale=Psi0+Sc, size=ndraw, random_state=seed)
    W = np.empty((ndraw, Nn))
    for s in range(ndraw):
        inv1 = np.linalg.solve(Sig_d[s], np.ones(Nn)); W[s] = inv1/inv1.sum()
    return W

Xstk = pd.read_csv("stocks_weekly.csv", index_col=0).values
W48 = niw_weight_draws(Xstk[-104:])
lo, hi = np.percentile(W48, 5, axis=0), np.percentile(W48, 95, axis=0)
straddle = np.mean((lo < 0) & (hi > 0))
order = np.argsort(W48.mean(0))
plt.figure(figsize=(12,5))
plt.boxplot([W48[:,i]*100 for i in order], showfliers=False,
            patch_artist=True, boxprops=dict(facecolor="#bee3f8"), medianprops=dict(color=BLUE))
plt.axhline(0, color="k", lw=.6); plt.xticks([], [])
plt.xlabel("48 stocks (sorted by mean weight)"); plt.ylabel("posterior weight (%)")
plt.title("Posterior of minimum-variance weights on 48 stocks: almost none can be signed")
plt.show()
print("Fraction of the 48 assets whose 90%% posterior weight interval straddles zero: %.0f%%" % (100*straddle))
print("At high dimension the data barely determine the sign of each position, let alone its size.")
print("A point optimiser hides this entirely; the Bayesian posterior lays it bare -- and it is the")
print("raw material for the robust allocation of Project 4.")
No description has been provided for this image 
Fraction of the 48 assets whose 90% posterior weight interval straddles zero: 60%
At high dimension the data barely determine the sign of each position, let alone its size.
A point optimiser hides this entirely; the Bayesian posterior lays it bare -- and it is the
raw material for the robust allocation of Project 4.

5. Summary¶

  • PyMC fits the full posterior of $(\mu,\Sigma)$; its posterior mean agrees with the closed-form NIW and shrinks the sample estimates.
  • The posterior predictive distribution inflates risk to account for parameter uncertainty.
  • The decisive Bayesian contribution: uncertainty flows into the allocation. Minimum-variance weights become a posterior distribution — frequently straddling zero — and the efficient frontier becomes a band, not a line. The plug-in method reports single confident numbers the data cannot support.

This uncertainty is not a nuisance to be hidden; it is the raw material of robust portfolio construction. Project 4 will optimise against this posterior — choosing allocations that are good across the whole band rather than optimal for one fragile point. First, though, Project 3 shows where a good prior comes from in the first place: the market equilibrium of Black–Litterman.