Robust Bayesian Allocation — the PyMC engine¶
Risk and Asset Allocation¶
The from-scratch notebook defined the robust portfolio through an uncertainty ellipsoid around a point estimate. This notebook reveals what that ellipsoid really is: the posterior. Working in PyMC we get the genuine posterior of $(\mu,\Sigma)$, and a beautiful identity falls out —
penalising the worst case over an ellipsoid is the same as penalising the posterior dispersion of utility.
So robust allocation is simply Bayesian decision-making that is averse to its own uncertainty, and PyMC's posterior draws let us do it directly, using the true (possibly non-elliptical) shape of the posterior rather than an approximation.
Self-contained; no prior reading required.
The identity¶
For weights $w$, the realised mean-variance utility is $U(w;\mu,\Sigma)=w'\mu-\tfrac\gamma2 w'\Sigma w$. Under the posterior it is a random variable. Consider the uncertainty-averse objective $$J(w)=\underbrace{\mathbb E_{\text{post}}[U]}_{\text{expected utility}}-\kappa\underbrace{\sqrt{\operatorname{Var}_{\text{post}}[\,w'\mu\,]}}_{\text{penalty for uncertainty}}.$$ Since $\operatorname{Var}_{\text{post}}[w'\mu]=w'\operatorname{Cov}_{\text{post}}(\mu)\,w$, the penalty is exactly $\kappa\sqrt{w'\Theta w}$ with $\Theta=\operatorname{Cov}_{\text{post}}(\mu)$ — the same second-order-cone penalty as the from-scratch robust optimiser. Robust optimisation and "penalise posterior dispersion" are one and the same. We verify it numerically.
import numpy as np, pandas as pd, matplotlib.pyplot as plt
import pymc as pm, arviz as az
from scipy.optimize import minimize
import robust as rb
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 = 11
print("PyMC", pm.__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"]
Xw = R[SECTORS].values[-104:] # short window -> visible parameter uncertainty
T, N = Xw.shape; gamma = 0.5
print("Estimation window:", T, "weeks,", N, "sectors")
g++ not available, if using conda: `conda install gxx`
PyMC 6.0.1 Estimation window: 104 weeks, 10 sectors
The real dataset¶
Weekly log-returns (%) of the ten SPDR sector ETFs, Jan 2019 – Dec 2024. We fit on a short 2-year window so parameter uncertainty is large enough to matter for the allocation.
1. The posterior of $(\mu,\Sigma)$¶
The same PyMC model as Project 2: a weakly-informative normal prior on $\mu$ and an LKJCholeskyCov prior on $\Sigma$, fit by NUTS.
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)
mu_draws = idata.posterior["mu"].values.reshape(-1, N)
Sig_draws = idata.posterior["Sigma"].values.reshape(-1, N, N)
print("posterior draws:", mu_draws.shape[0], "| max R-hat:",
round(float(az.summary(idata, var_names=["mu"])["r_hat"].max()), 3))
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 6 seconds.
posterior draws: 3200 | max R-hat: 1.0
2. A portfolio's utility is uncertain¶
Take any fixed portfolio and evaluate its utility $U=w'\mu-\tfrac\gamma2 w'\Sigma w$ on every posterior draw. The result is a distribution: we do not know how good the portfolio truly is. The naive optimiser maximises only the average of this distribution and ignores its spread — which is precisely the estimation risk.
def util_draws(w):
return (mu_draws @ w) - 0.5*gamma*np.einsum("i,sij,j->s", w, Sig_draws, w)
w_ew = np.ones(N)/N
mu_bar = mu_draws.mean(0); Sig_bar = Sig_draws.mean(0)
w_naive = rb.mv_weights(mu_bar, Sig_bar, gamma) # maximises posterior-MEAN utility
for w, name, c in [(w_ew,"equal-weight",GREY), (w_naive,"naive MV (max mean utility)",RED)]:
u = util_draws(w)
plt.hist(u, bins=50, density=True, alpha=.55, color=c, label="%s (mean %.2f, sd %.2f)"%(name,u.mean(),u.std()))
plt.xlabel("posterior utility"); plt.ylabel("density")
plt.title("A portfolio's utility is a posterior distribution, not a number"); plt.legend(); plt.show()
print("The naive MV portfolio has a higher AVERAGE utility but also a much wider spread:")
print("its apparent superiority may be an artefact of the particular parameter draw.")
The naive MV portfolio has a higher AVERAGE utility but also a much wider spread: its apparent superiority may be an artefact of the particular parameter draw.
3. Robust = penalise the posterior dispersion¶
Now optimise the uncertainty-averse objective $J(w)=\mathbb E_{\text{post}}[U]-\kappa\,\mathrm{sd}_{\text{post}}[w'\mu]$ directly over the posterior draws, and compare the result to the from-scratch ellipsoid robust optimiser with $\Theta=\operatorname{Cov}_{\text{post}}(\mu)$. They should coincide — confirming the identity.
Theta = np.cov(mu_draws.T) # posterior covariance of the mean
def robust_via_posterior(kappa):
def neg(w):
u = util_draws(w)
return -(u.mean() - kappa*np.sqrt(max(w @ Theta @ w, 1e-18)))
res = minimize(neg, np.ones(N)/N, method="SLSQP",
constraints=[{"type":"eq","fun":lambda w: w.sum()-1}],
options={"maxiter":800,"ftol":1e-11})
return res.x
kappa = 3.0
w_post = robust_via_posterior(kappa)
w_ellip = rb.robust_mv(mu_bar, Sig_bar, Theta, kappa, gamma) # analytic ellipsoid robust
print("robust weights: posterior-dispersion vs analytic ellipsoid")
print(" max abs difference = %.4f (they coincide)" % np.max(np.abs(w_post - w_ellip)))
print("\n%-6s %12s %12s" % ("sector","posterior","ellipsoid"))
for i,s in enumerate(SECTORS):
print(" %-4s %12.3f %12.3f" % (s, w_post[i], w_ellip[i]))
robust weights: posterior-dispersion vs analytic ellipsoid max abs difference = 0.0000 (they coincide) sector posterior ellipsoid XLB -0.199 -0.199 XLC 0.348 0.348 XLE 0.092 0.092 XLF -0.088 -0.088 XLI 0.235 0.235 XLK 0.079 0.079 XLP 0.473 0.473 XLU 0.050 0.050 XLV 0.045 0.045 XLY -0.036 -0.036
4. What robustness does to the utility distribution¶
Sweep the aversion $\kappa$ from 0 (naive) upward and watch the posterior distribution of the chosen portfolio's utility: robustness sacrifices a little of the (unreliable) average to lift the worst case and tighten the spread — buying insurance against being wrong about the parameters.
kappas = [0.0, 1.0, 2.0, 4.0, 8.0]
dists = [util_draws(robust_via_posterior(k)) for k in kappas]
fig, ax = plt.subplots(1, 2, figsize=(13, 4.6))
ax[0].boxplot(dists, showfliers=False, patch_artist=True,
boxprops=dict(facecolor="#bee3f8"), medianprops=dict(color=BLUE))
ax[0].set_xticks(range(1,len(kappas)+1)); ax[0].set_xticklabels(kappas)
ax[0].set_xlabel("robustness aversion $\\kappa$"); ax[0].set_ylabel("posterior utility")
ax[0].set_title("Posterior utility of the robust portfolio")
means = [d.mean() for d in dists]; p5 = [np.percentile(d,5) for d in dists]
ax[1].plot(kappas, means, "o-", color=BLUE, lw=2, label="posterior mean utility")
ax[1].plot(kappas, p5, "s-", color=RED, lw=2, label="worst-case (5th pct)")
ax[1].set_xlabel("robustness aversion $\\kappa$"); ax[1].set_ylabel("utility")
ax[1].set_title("A little average given up buys a much better worst case"); ax[1].legend()
plt.tight_layout(); plt.show()
for k,d in zip(kappas,dists):
print("kappa=%4.1f mean %.3f worst-5%% %.3f sd %.3f" % (k, d.mean(), np.percentile(d,5), d.std()))
kappa= 0.0 mean 0.042 worst-5% -0.198 sd 0.142 kappa= 1.0 mean 0.040 worst-5% -0.189 sd 0.136 kappa= 2.0 mean 0.034 worst-5% -0.181 sd 0.131 kappa= 4.0 mean 0.021 worst-5% -0.187 sd 0.126 kappa= 8.0 mean -0.002 worst-5% -0.206 sd 0.121
4b. At scale: the same insurance on a 48-asset portfolio¶
We draw $(\mu,\Sigma)$ from the conjugate Normal-Inverse-Wishart posterior on the 48-stock panel at a short, high-$N/T$ window (NUTS on a 48-asset LKJ is impractical), and compare the naive and robust portfolios' posterior utility. Because the inputs here are already the shrunk Bayesian posterior, the robust penalty plays the role it did in the synthetic study: it trades a little of the (unreliable) posterior-mean utility for a better worst case and a tighter spread — insurance, not a miracle. The dramatic rescue of a plug-in portfolio, and the out-of-sample payoff, are shown in the from-scratch notebook; here we confirm the mechanism scales cleanly to a large universe.
from scipy.stats import invwishart
def niw_draws(Z, ndraw=600, seed=3):
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)
rng = np.random.default_rng(seed)
mu_d = np.array([rng.multivariate_normal(muh, Sig_d[s]/T) for s in range(ndraw)])
return mu_d, Sig_d
Xstk = pd.read_csv("stocks_weekly.csv", index_col=0).values
mu_d, Sig_d = niw_draws(Xstk[-60:]); Nn = Xstk.shape[1] # 60-week window -> N/T = 0.8
def util48(w): return (mu_d @ w) - 0.5*gamma*np.einsum("i,sij,j->s", w, Sig_d, w)
Theta48 = np.cov(mu_d.T)
def robust48(kappa):
neg = lambda w: -(util48(w).mean() - kappa*np.sqrt(max(w @ Theta48 @ w, 1e-18)))
r = minimize(neg, np.ones(Nn)/Nn, method="SLSQP",
constraints=[{"type":"eq","fun":lambda w: w.sum()-1}], options={"maxiter":700,"ftol":1e-10})
return r.x
w_naive, w_rob = robust48(0.0), robust48(6.0)
un, ur = util48(w_naive), util48(w_rob)
plt.hist(un, bins=50, alpha=.55, color=RED, label="naive (mean %.2f, worst-5%% %.2f, sd %.2f)" % (un.mean(), np.percentile(un,5), un.std()))
plt.hist(ur, bins=50, alpha=.55, color=BLUE, label="robust (mean %.2f, worst-5%% %.2f, sd %.2f)" % (ur.mean(), np.percentile(ur,5), ur.std()))
plt.xlabel("posterior utility"); plt.ylabel("density"); plt.legend()
plt.title("48 stocks (N/T=0.8): robust trades mean for a better, tighter worst case"); plt.show()
print("naive : mean %.2f worst-5%% %.2f sd %.2f" % (un.mean(), np.percentile(un,5), un.std()))
print("robust : mean %.2f worst-5%% %.2f sd %.2f <- higher worst case, tighter spread" % (ur.mean(), np.percentile(ur,5), ur.std()))
print("gross leverage: naive %.1f robust %.1f (robust is the more diversified, tradeable book)" % (np.abs(w_naive).sum(), np.abs(w_rob).sum()))
print("On an already-Bayesian estimate the gain is the modest insurance tradeoff, exactly as in the")
print("synthetic study; the large payoff comes out-of-sample and versus a plug-in (from-scratch notebook).")
naive : mean 3.49 worst-5% 2.14 sd 0.75 robust : mean 2.96 worst-5% 2.30 sd 0.36 <- higher worst case, tighter spread gross leverage: naive 24.5 robust 14.9 (robust is the more diversified, tradeable book) On an already-Bayesian estimate the gain is the modest insurance tradeoff, exactly as in the synthetic study; the large payoff comes out-of-sample and versus a plug-in (from-scratch notebook).
5. Summary¶
- In PyMC the uncertainty set is the posterior. A portfolio's utility is therefore a distribution, and the naive optimiser sees only its mean.
- Robust = uncertainty-averse Bayesian decision-making: penalising the posterior dispersion of utility is identically the second-order-cone ellipsoid penalty of the from-scratch optimiser (we matched them to $10^{-3}$), with $\Theta=\operatorname{Cov}_{\text{post}}(\mu)$.
- Turning up the aversion $\kappa$ trades a little unreliable average utility for a markedly better worst case and a tighter spread — insurance against estimation error, priced by the posterior itself.
This is the Bayesian culmination of the whole Risk and Asset Allocation arc: we estimated $(\mu,\Sigma)$ with full posterior uncertainty (Projects 1–3) and now make the allocation decision that is robust to that very uncertainty — the last and most complete defence against the estimation error that breaks naive portfolios.