Bayesian Estimation & Estimation Risk

Risk and Asset Allocation — the from-scratch engine

"The efficient frontier you compute from historical data is a mirage. The portfolios on it look wonderful — until you actually hold them."

Self-contained; no prior reading of Meucci's Risk and Asset Allocation (Springer 2005, chapters 7–9) required.

Where we are

Project 1 showed that the sample mean and covariance are noisy, and that shrinkage — pulling them toward a target — helps. This notebook turns that idea into a full probability model and then confronts the question a portfolio manager should lose sleep over:

How much does estimation error actually cost me, and can I do something about it?

We answer in four movements:

1. Bayesian estimation. Put a Normal-Inverse-Wishart (NIW) prior on $(\mu,\Sigma)$ and update it with data. The posterior is closed-form, and its mean is exactly a shrinkage estimator — the Bayesian foundation under Project 1.
2. Validate with MCMC. The closed form is a lucky special case. We build a from-scratch random-walk Metropolis sampler and confirm it reproduces the analytic posterior — the tool we will need when conjugacy is lost.
3. The predictive distribution. Integrating out the unknown $(\mu,\Sigma)$ gives the distribution of a future return. Its covariance is inflated by parameter uncertainty — and optimising against it is what "Bayesian allocation" means.
4. Estimation risk. We measure the deception of the sample efficient frontier (it promises more than it delivers) and show that Bayesian allocation cuts the opportunity cost of estimation error, dramatically so when data is scarce.

Throughout we use both a synthetic market (known truth → we can measure error) and the real sector-ETF panel.

import numpy as np, pandas as pd, matplotlib.pyplot as plt
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, PURP = "#2b6cb0","#dd6b20","#2f855a","#c53030","#718096","#6b46c1"
np.set_printoptions(precision=3, suppress=True)

R = pd.read_csv("sector_etf_weekly.csv", index_col=0, parse_dates=True)
SECTORS = [c for c in R.columns if c != "SPY"]
X = R[SECTORS].values
T, N = X.shape
print("Real panel:", N, "sector ETFs over", T, "weeks (weekly log-returns, %)")

Real panel: 10 sector ETFs over 310 weeks (weekly log-returns, %)

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 (via yfinance). Each ETF holds the S&P 500 members of one GICS sector:

Ticker	Sector	Ticker	Sector
XLB	Materials	XLP	Consumer Staples
XLE	Energy	XLU	Utilities
XLF	Financials	XLV	Health Care
XLI	Industrials	XLY	Consumer Discretionary
XLK	Technology	XLC	Communication Services
		SPY	Market (S&P 500 benchmark)

The sectors are strongly, positively correlated (they all ride the market), so their sample covariance is noisy and their sample means unreliable — exactly the setting where Bayesian estimation earns its keep. We validate every claim first on a synthetic market whose true $(\mu,\Sigma)$ we know.

1. Bayesian estimation with the Normal-Inverse-Wishart prior

The data model. Weekly returns are treated as independent draws from a normal distribution: $$x_t \mid \mu,\Sigma \;\sim\; \mathcal N(\mu,\Sigma), \qquad t=1,\dots,T.$$ The unknowns are the mean vector $\mu$ and covariance $\Sigma$.

The prior. Before seeing the data we hold beliefs about $(\mu,\Sigma)$. The conjugate choice — the one that makes the update clean — is the Normal-Inverse-Wishart: $$\Sigma \sim \mathcal{IW}(\nu_0,\ \nu_0\Sigma_0), \qquad \mu \mid \Sigma \sim \mathcal N\!\big(\mu_0,\ \Sigma/T_0\big).$$ Here $\mu_0,\Sigma_0$ are our best guesses and the two scalars are confidence dials:

• $T_0$ = how many observations' worth of conviction we have in the prior mean $\mu_0$,
• $\nu_0$ = the same for the prior covariance $\Sigma_0$.

Large $T_0,\nu_0$ = a strong prior that the data must fight hard to move; small values = a weak prior that lets the data speak.

The posterior. Because the prior is conjugate, the posterior is again Normal-Inverse-Wishart, with hyper-parameters that blend prior and data: $$ \begin{aligned} T_1 &= T_0 + T, & \mu_1 &= \frac{T_0\mu_0 + T\,\hat\mu}{T_0+T},\\ \nu_1 &= \nu_0 + T, & \nu_1\Sigma_1 &= \nu_0\Sigma_0 + T\hat\Sigma + \frac{T\,T_0}{T_0+T}(\mu_0-\hat\mu)(\mu_0-\hat\mu)'. \end{aligned} $$ Look at $\mu_1$: it is a weighted average of the prior mean and the sample mean — a shrinkage estimator, with weight $T_0/(T_0+T)$ on the prior. This is the precise Bayesian statement of Project 1's idea.

# A weak, honest prior for the sector ETFs: no directional view on returns (mu0=0),
# covariance centred on the average sample variance, confidence of a few observations.
mu_hat, S_hat = niw.sample_moments(X)
mu0    = np.zeros(N)
Sigma0 = np.mean(np.diag(S_hat)) * np.eye(N)          # sphere with the average variance
T0, nu0 = 5.0, N + 4                                  # weak prior (~5 obs of conviction)

post = niw.niw_posterior(X, mu0, T0, Sigma0, nu0)
mu1, Sig1_mean = niw.post_mean_cov(post)

print("prior mean weight on mu  = T0/(T0+T) = %.3f" % (T0/(T0+T)))
print("\nsector   sample mean   posterior mean")
for i,s in enumerate(SECTORS):
    print("  %-4s   %8.3f      %8.3f" % (s, mu_hat[i], mu1[i]))

prior mean weight on mu  = T0/(T0+T) = 0.016

sector   sample mean   posterior mean
  XLB       0.205         0.202
  XLC       0.291         0.286
  XLE       0.215         0.212
  XLF       0.257         0.253
  XLI       0.262         0.257
  XLK       0.469         0.461
  XLP       0.184         0.181
  XLU       0.177         0.174
  XLV       0.183         0.180
  XLY       0.274         0.269

# Prior vs posterior, seen as distributions: draw from each NIW and look at marginals
rng = np.random.default_rng(1)
prior_hp = dict(mu1=mu0, T1=T0, Sigma1=Sigma0, nu1=nu0)           # prior as a NIW 'posterior'
mus_pri, Sig_pri = niw.niw_rand(prior_hp, 8000, rng)
mus_pos, Sig_pos = niw.niw_rand(post,     8000, rng)

j = SECTORS.index("XLK")                                          # technology
fig, ax = plt.subplots(1, 2, figsize=(12, 4.3))
ax[0].hist(mus_pri[:,j], bins=60, density=True, color=GREY,  alpha=.6, label="prior")
ax[0].hist(mus_pos[:,j], bins=60, density=True, color=BLUE,  alpha=.7, label="posterior")
ax[0].axvline(mu_hat[j], color=RED, ls="--", label="sample mean")
ax[0].set_title("Belief about mean return of %s" % SECTORS[j]); ax[0].set_xlabel("weekly %"); ax[0].legend()
ax[1].hist(Sig_pri[:,j,j], bins=60, density=True, color=GREY, alpha=.6, label="prior")
ax[1].hist(Sig_pos[:,j,j], bins=60, density=True, color=BLUE, alpha=.7, label="posterior")
ax[1].axvline(S_hat[j,j], color=RED, ls="--", label="sample variance")
ax[1].set_title("Belief about variance of %s" % SECTORS[j]); ax[1].set_xlabel("variance"); ax[1].legend()
plt.tight_layout(); plt.show()
print("The data (weak prior, T=312) dominates: the posterior concentrates near the sample values,")
print("but as a full distribution -- we now carry uncertainty about mu and Sigma, not point guesses.")

The data (weak prior, T=312) dominates: the posterior concentrates near the sample values,
but as a full distribution -- we now carry uncertainty about mu and Sigma, not point guesses.

2. Validating the posterior with a from-scratch MCMC sampler

The NIW posterior above is available in closed form only because the prior is conjugate. Change the model even slightly — fat-tailed returns, a prior on correlations only, parameter constraints — and no formula exists. Then we must sample the posterior numerically with Markov-chain Monte Carlo (MCMC).

To show the machinery (and to check our formulas), we build a random-walk Metropolis sampler from scratch on a small 2-asset market. It targets the log-posterior $$\log p(\mu,\Sigma \mid \text{data}) = \underbrace{\sum_t \log \mathcal N(x_t\mid\mu,\Sigma)}_{\text{likelihood}} + \underbrace{\log\mathcal N(\mu\mid\mu_0,\Sigma/T_0) + \log\mathcal{IW}(\Sigma\mid\nu_0,\nu_0\Sigma_0)}_{\text{prior}},$$ parameterising $\Sigma = LL'$ through its Cholesky factor $L$ so every proposal stays a valid covariance. At each step we perturb $\mu$ and $L$ and accept with the Metropolis probability $\min(1, e^{\Delta\log p})$. If our derivation is right, the sampler's output must match the analytic NIW posterior.

# A 2-asset market so we can visualise the joint posterior
np.random.seed(0)
mu_true2 = np.array([1.0, 0.7]); Sig_true2 = np.array([[4.0, 1.2],[1.2, 2.0]])
rng2 = np.random.default_rng(0)
X2 = rng2.multivariate_normal(mu_true2, Sig_true2, size=52)
prior2 = dict(mu0=np.zeros(2), T0=10.0, Sigma0=3*np.eye(2), nu0=10.0)
post2 = niw.niw_posterior(X2, **prior2)

mc = niw.metropolis_niw(X2, **prior2, n_draws=20000, burn=4000, rng=rng2, thin=2)
mus_an, Sig_an = niw.niw_rand(post2, 20000, rng2)          # analytic draws
print("Metropolis acceptance rate: %.2f  (a healthy 20-60%% band)" % mc["accept"])

Metropolis acceptance rate: 0.49  (a healthy 20-60% band)

fig, ax = plt.subplots(1, 3, figsize=(14, 4))
ax[0].plot(mc["mu"][:1500,0], color=BLUE, lw=.6); ax[0].set_title("MCMC trace: $\\mu_1$"); ax[0].set_xlabel("iteration")
# posterior of mu_1: analytic vs MCMC
ax[1].hist(mus_an[:,0], bins=60, density=True, color=GREY, alpha=.6, label="analytic NIW")
ax[1].hist(mc["mu"][:,0], bins=60, density=True, color=BLUE, alpha=.6, label="MCMC")
ax[1].set_title("Posterior of $\\mu_1$"); ax[1].legend()
# posterior of covariance term Sigma_12
ax[2].hist(Sig_an[:,0,1], bins=60, density=True, color=GREY, alpha=.6, label="analytic NIW")
ax[2].hist(mc["Sigma"][:,0,1], bins=60, density=True, color=GREEN, alpha=.6, label="MCMC")
ax[2].set_title("Posterior of covariance $\\Sigma_{12}$"); ax[2].legend()
plt.tight_layout(); plt.show()
print("analytic vs MCMC posterior means:")
print("  mu    : %s  vs  %s" % (mus_an.mean(0).round(3), mc["mu"].mean(0).round(3)))
print("  Sig12 : %.3f       vs  %.3f" % (Sig_an[:,0,1].mean(), mc["Sigma"][:,0,1].mean()))
print("The from-scratch sampler reproduces the closed-form posterior -> our derivation checks out.")

analytic vs MCMC posterior means:
  mu    : [0.776 0.742]  vs  [0.783 0.743]
  Sig12 : 0.690       vs  0.663
The from-scratch sampler reproduces the closed-form posterior -> our derivation checks out.

3. The posterior predictive: what a future return looks like

A portfolio is held into the future, so the object we actually care about is the distribution of the next return $x_{\text{new}}$ given the data — with the unknown parameters integrated out: $$p(x_{\text{new}}\mid\text{data}) = \int p(x_{\text{new}}\mid\mu,\Sigma)\,p(\mu,\Sigma\mid\text{data})\,d\mu\,d\Sigma.$$ For the NIW model this is a multivariate Student-$t$ distribution, and its first two moments are $$\mathbb E[x_{\text{new}}] = \mu_1, \qquad \operatorname{Cov}[x_{\text{new}}] = \Big(1+\tfrac1{T_1}\Big)\,\frac{\nu_1\Sigma_1}{\nu_1-N-1}.$$ The predictive covariance is larger than the plug-in covariance: it adds the risk that comes from not knowing the parameters. This is the crucial difference between a naive optimiser (which pretends the estimates are exact) and a Bayesian one (which admits they are not).

mu_pred, S_pred = niw.predictive_moments(post)
_, S_postmean   = niw.post_mean_cov(post)
infl = np.diag(S_pred) / np.diag(S_hat)
print("Predictive vs plug-in (sample) variances -- inflation from parameter uncertainty:")
for i,s in enumerate(SECTORS):
    print("  %-4s  sample var %7.3f   predictive var %7.3f   x%.3f" % (s, S_hat[i,i], S_pred[i,i], infl[i]))
print("\nWith T=312 the inflation is small; with short samples it is large -- and that is exactly")
print("when ignoring it is most dangerous, as the estimation-risk experiment below shows.")

Predictive vs plug-in (sample) variances -- inflation from parameter uncertainty:
  XLB   sample var   9.553   predictive var   9.966   x1.043
  XLC   sample var   8.180   predictive var   8.603   x1.052
  XLE   sample var  23.265   predictive var  23.590   x1.014
  XLF   sample var  11.910   predictive var  12.308   x1.033
  XLI   sample var   9.657   predictive var  10.069   x1.043
  XLK   sample var  11.350   predictive var  11.754   x1.036
  XLP   sample var   4.114   predictive var   4.562   x1.109
  XLU   sample var   9.975   predictive var  10.385   x1.041
  XLV   sample var   6.179   predictive var   6.614   x1.070
  XLY   sample var  11.472   predictive var  11.874   x1.035

With T=312 the inflation is small; with short samples it is large -- and that is exactly
when ignoring it is most dangerous, as the estimation-risk experiment below shows.

4. Estimation risk: the deception of the sample efficient frontier

The efficient frontier is the set of portfolios with the lowest risk for each level of expected return. With budget constraint $w'\mathbf 1 = 1$ and short sales allowed, the minimum-variance weights achieving a target return $m$ have the closed form $$w(m) = \Sigma^{-1}\big[\mathbf 1\,(C - Bm) + \mu\,(Am - B)\big]/D,\quad A=\mathbf 1'\Sigma^{-1}\mathbf 1,\ B=\mathbf 1'\Sigma^{-1}\mu,\ C=\mu'\Sigma^{-1}\mu,\ D=AC-B^2.$$

Now the experiment that exposes estimation risk. Fix a known synthetic market. Repeatedly:

1. draw a sample of $T$ returns;
2. compute the efficient frontier from the sample $(\hat\mu,\hat\Sigma)$ — this is the frontier the manager sees and reports (the claimed frontier);
3. take those same sample-optimal weights and score them at the true $(\mu,\Sigma)$ — this is what they actually deliver.

Averaged over samples, the claimed frontier sits above and to the left of the truth (it looks too good), while the realised performance sits below and to the right of the true frontier. The gap is pure estimation error — the "deception".

# Known 8-asset market; weak prior; short samples make estimation error bite
Na = 8
mu_t, Sig_t, rngd = niw.make_true_market(Na, rho=0.5, seed=3)
prior_d = dict(mu0=np.zeros(Na), T0=20.0, Sigma0=np.mean(np.diag(Sig_t))*np.eye(Na), nu0=20.0)
res = niw.deception_experiment(mu_t, Sig_t, T=24, prior=prior_d, n_rep=600, n_pts=25, rng=rngd)

cl, ts, tf = res["claimed_sample"], res["true_sample"], res["true_frontier"]
plt.figure(figsize=(8.5,5.5))
plt.plot(tf[:,0], tf[:,1], "-",  color=GREEN, lw=2.5, label="TRUE frontier (unattainable ideal)")
plt.plot(cl[:,0], cl[:,1], "--", color=BLUE,  lw=2,   label="CLAIMED frontier (what the sample shows)")
plt.plot(ts[:,0], ts[:,1], "-o", color=RED,   lw=2, ms=3, label="TRUE performance of sample-optimal portfolios")
plt.xlabel("volatility (weekly %)"); plt.ylabel("expected return (weekly %)")
plt.title("The sample efficient frontier is a mirage (N=8, T=24, 600 samples)")
plt.legend(); plt.show()
print("Blue (claimed) is optimistic; red (reality) is worse on BOTH risk and return.")
print("The manager who trusts the blue curve is systematically disappointed.")

Blue (claimed) is optimistic; red (reality) is worse on BOTH risk and return.
The manager who trusts the blue curve is systematically disappointed.

5. Bayesian allocation reduces the cost of estimation error

The deception is diagnosed; now the cure. Instead of plugging in the sample moments, plug in the posterior-predictive moments $(\mu_1,\ \text{predictive }\Sigma)$. Because the predictive mean is shrunk and the predictive covariance is inflated, the resulting portfolio is more conservative — and, crucially, it performs better out of sample.

We quantify with the certainty-equivalent opportunity cost. An investor with risk aversion $\gamma$ ideally holds $$w^\star = \arg\max_{w'\mathbf 1=1}\ w'\mu - \tfrac{\gamma}{2}w'\Sigma w,$$ earning certainty-equivalent $\text{CE}_{\text{true}}$ when $(\mu,\Sigma)$ are the truth. A real investor must estimate the inputs; scoring their chosen weights at the truth gives a realised CE. The opportunity cost $\text{CE}_{\text{true}} - \text{CE}_{\text{realised}} \ge 0$ is the price of estimation error. We compare the sample and Bayesian input estimators across sample sizes $T$.

gamma = 0.5
Ts = [16, 24, 40, 80, 160, 320]
oc_sample, oc_bayes = [], []
for Tn in Ts:
    s, b, ce_true = niw.opportunity_cost_experiment(mu_t, Sig_t, Tn, prior_d, gamma, n_rep=1500, rng=rngd)
    oc_sample.append(s.mean()); oc_bayes.append(b.mean())

fig, ax = plt.subplots(1, 2, figsize=(13, 4.6))
ax[0].plot(Ts, oc_sample, "s-", color=RED,  lw=2, label="sample (plug-in)")
ax[0].plot(Ts, oc_bayes,  "o-", color=BLUE, lw=2, label="Bayesian (predictive)")
ax[0].set_yscale("log"); ax[0].set_xlabel("sample size T"); ax[0].set_ylabel("opportunity cost (log)")
ax[0].axhline(ce_true, color=GREY, ls=":", label="CE of the true optimum")
ax[0].set_title("Cost of estimation error vs sample size"); ax[0].legend()
red = [100*(1-b/s) for s,b in zip(oc_sample, oc_bayes)]
ax[1].plot(Ts, red, "o-", color=GREEN, lw=2)
ax[1].set_xlabel("sample size T"); ax[1].set_ylabel("Bayesian reduction in opportunity cost (%)")
ax[1].set_title("Bayesian advantage is largest when data is scarce"); ax[1].set_ylim(0, 100)
plt.tight_layout(); plt.show()
for Tn,s,b in zip(Ts, oc_sample, oc_bayes):
    print("T=%3d   sample OC %8.3f   bayes OC %6.3f   (%.0f%% lower)" % (Tn, s, b, 100*(1-b/s)))

T= 16   sample OC   13.201   bayes OC  0.170   (99% lower)
T= 24   sample OC    3.155   bayes OC  0.166   (95% lower)
T= 40   sample OC    0.927   bayes OC  0.145   (84% lower)
T= 80   sample OC    0.310   bayes OC  0.104   (66% lower)
T=160   sample OC    0.126   bayes OC  0.066   (48% lower)
T=320   sample OC    0.057   bayes OC  0.039   (31% lower)

The same curve on real data — without a known truth

On real data we cannot compute the opportunity cost, because there is no known $(\mu,\Sigma)$ to define the ideal $\text{CE}_{\text{true}}$. But we can draw the honest analog: the realized out-of-sample certainty-equivalent of each strategy as a function of the estimation-window length $W$ (the real-data stand-in for sample size $T$). Each week we estimate the inputs from the trailing $W$ weeks, form the optimal portfolio, hold it one week; the strategy's realized-return series gives $$\text{CE}^{\text{OOS}} = \operatorname{mean}\big(r^{\text{port}}\big) - \tfrac{\gamma}{2}\operatorname{var}\big(r^{\text{port}}\big),\qquad \text{higher is better.}$$ We expect the same shape as the synthetic curve: the sample strategy collapses when the window is short (unstable weights), the Bayesian one stays stable, and the two converge as $W$ grows.

def strat_returns(X, W, gamma, bayes):
    rets = []
    for t in range(W, len(X)):
        Xw = X[t-W:t]
        if bayes:
            p = niw.niw_posterior(Xw, np.zeros(N), 5.0,
                                  np.mean(np.diag(niw.sample_moments(Xw)[1]))*np.eye(N), N+4)
            mu_e, S_e = niw.predictive_moments(p)
        else:
            mu_e, S_e = niw.sample_moments(Xw)
        rets.append(X[t] @ niw.optimal_mv(mu_e, S_e, gamma))
    return np.array(rets)

gamma_r = 1.0
Ws = [16, 20, 26, 40, 52, 78, 104, 156]
ce_s, ce_b = [], []
for W in Ws:
    rs = strat_returns(X, W, gamma_r, False); rb = strat_returns(X, W, gamma_r, True)
    ce_s.append(rs.mean() - 0.5*gamma_r*rs.var())
    ce_b.append(rb.mean() - 0.5*gamma_r*rb.var())
ce_s, ce_b = np.array(ce_s), np.array(ce_b)

fig, ax = plt.subplots(1, 2, figsize=(13, 4.6))
ax[0].plot(Ws, ce_s, "s-", color=RED,  lw=2, label="sample (plug-in)")
ax[0].plot(Ws, ce_b, "o-", color=BLUE, lw=2, label="Bayesian (predictive)")
ax[0].set_ylim(-14, 0.5); ax[0].axhline(0, color="k", lw=.5)
ax[0].set_xlabel("estimation window W (weeks)"); ax[0].set_ylabel("realized OOS certainty-equivalent")
ax[0].set_title("Real sector ETFs: OOS CE vs window (higher = better)"); ax[0].legend()
for W, c in zip(Ws, ce_s):                                  # annotate off-scale sample points
    if c < -14:
        ax[0].annotate("%.0f" % c, (W, -13.6), color=RED, ha="center", fontsize=8, fontweight="bold")
ax[1].plot(Ws, ce_b - ce_s, "o-", color=GREEN, lw=2)
ax[1].set_yscale("symlog"); ax[1].set_xlabel("estimation window W (weeks)")
ax[1].set_ylabel("Bayesian advantage  CE$_{bayes}$ − CE$_{sample}$")
ax[1].set_title("Bayesian advantage — largest when the window is short")
plt.tight_layout(); plt.show()
for W, cs, cb in zip(Ws, ce_s, ce_b):
    print("  W=%3d   sample CE %10.3f   bayes CE %8.3f" % (W, cs, cb))
print("\nSame shape as the synthetic opportunity-cost curve: the SAMPLE strategy's CE collapses at")
print("short windows (off-scale values annotated in red), while the Bayesian strategy stays stable.")
print("The advantage shrinks as the window lengthens -- the estimation-risk story, on live data.")

  W= 16   sample CE    -17.654   bayes CE   -2.220
  W= 20   sample CE     -7.514   bayes CE   -2.167
  W= 26   sample CE     -3.741   bayes CE   -2.102
  W= 40   sample CE     -2.690   bayes CE   -2.223
  W= 52   sample CE     -2.561   bayes CE   -2.447
  W= 78   sample CE     -1.555   bayes CE   -1.404
  W=104   sample CE     -1.696   bayes CE   -1.534
  W=156   sample CE     -1.634   bayes CE   -1.581

Same shape as the synthetic opportunity-cost curve: the SAMPLE strategy's CE collapses at
short windows (off-scale values annotated in red), while the Bayesian strategy stays stable.
The advantage shrinks as the window lengthens -- the estimation-risk story, on live data.

6. Dimensionality: the estimation-risk tax scales with N/T

The opportunity-cost experiment above used a small market where estimation error was made severe by a short sample. In practice the same severity arrives naturally through dimensionality: the more assets $N$ relative to the window $T$, the noisier the inputs and the heavier the "tax" a naive optimiser pays. Ten sectors ($N/T\approx0.1$) is an easy problem; a realistic universe of dozens of names is not.

We confirm it on 48 individual large-cap US stocks (same weekly period). A rolling out-of-sample backtest compares the sample and Bayesian mean-variance portfolios against the estimation-error-free 1/N benchmark, at a long ($T=104$, $N/T\approx0.46$) and a short ($T=60$, $N/T\approx0.80$) window.

Sstk = pd.read_csv("stocks_weekly.csv", index_col=0, parse_dates=True).values
def oos(Xd, window, which, gamma=0.5, ppy=52):
    Nn = Xd.shape[1]; rets = []
    for t in range(window, len(Xd)):
        Xw = Xd[t-window:t]
        if which == "1/N":
            w = np.ones(Nn)/Nn
        else:
            if which == "sample":
                mu, S = niw.sample_moments(Xw)
            else:  # bayes
                p = niw.niw_posterior(Xw, np.zeros(Nn), 5.0,
                                      np.mean(np.diag(niw.sample_moments(Xw)[1]))*np.eye(Nn), Nn+4)
                mu, S = niw.predictive_moments(p)
            w = niw.optimal_mv(mu, S, gamma)
        rets.append(Xd[t] @ w)
    r = np.array(rets); return r.mean()/r.std()*np.sqrt(ppy), r.std()*np.sqrt(ppy)

print("OOS annualised Sharpe (volatility %) on 48 stocks:\n")
print("%-14s %20s %20s" % ("strategy", "T=104 (N/T=0.46)", "T=60 (N/T=0.80)"))
res = {}
for which in ["sample", "bayes", "1/N"]:
    a = oos(Sstk, 104, which); b = oos(Sstk, 60, which); res[which] = (a, b)
    print("%-14s %14.2f (%4.0f%%) %14.2f (%4.0f%%)" % (which, a[0], a[1], b[0], b[1]))

OOS annualised Sharpe (volatility %) on 48 stocks:

strategy           T=104 (N/T=0.46)      T=60 (N/T=0.80)
sample                   0.58 (  29%)           0.43 ( 199%)
bayes                    0.77 (  14%)           0.94 (  16%)
1/N                      0.86 (  16%)           0.77 (  20%)

labels = ["sample MV", "Bayesian MV", "1/N"]
keys = ["sample", "bayes", "1/N"]
xx = np.arange(3)
plt.bar(xx-0.2, [res[k][0][0] for k in keys], width=0.4, color=BLUE,   label="T=104 (N/T=0.46)")
plt.bar(xx+0.2, [res[k][1][0] for k in keys], width=0.4, color=ORANGE, label="T=60 (N/T=0.80)")
plt.xticks(xx, labels); plt.ylabel("OOS annualised Sharpe"); plt.axhline(0, color="k", lw=.6)
plt.title("48 stocks: the Bayesian advantage grows as N/T rises"); plt.legend(); plt.show()
print("At the shorter window (higher N/T) the sample optimiser deteriorates sharply, while the")
print("Bayesian portfolio improves and overtakes even the formidable 1/N benchmark -- the clean win")
print("that the low-dimensional sector universe could not deliver. Estimation risk, and the value of")
print("addressing it, both scale with dimensionality.")

At the shorter window (higher N/T) the sample optimiser deteriorates sharply, while the
Bayesian portfolio improves and overtakes even the formidable 1/N benchmark -- the clean win
that the low-dimensional sector universe could not deliver. Estimation risk, and the value of
addressing it, both scale with dimensionality.

The canonical benchmark: Fama–French 48 industries

We confirm the same on the dataset the academic estimation-risk literature uses (DeMiguel–Garlappi–Uppal 2009): Ken French's 48 industry portfolios, value-weighted monthly, 1970–2024. Rolling windows of 60 and 120 months give $N/T=0.8$ and $0.4$.

FF = pd.read_csv("ff_industries_monthly.csv", index_col=0).values
print("Fama-French 48 industries, OOS annualised Sharpe (volatility %):\n")
print("%-14s %20s %20s" % ("strategy", "T=120mo (N/T=0.40)", "T=60mo (N/T=0.80)"))
ffres = {}
for which in ["sample", "bayes", "1/N"]:
    a = oos(FF, 120, which, ppy=12); b = oos(FF, 60, which, ppy=12); ffres[which] = (a, b)
    print("%-14s %14.2f (%4.0f%%) %14.2f (%4.0f%%)" % (which, a[0], a[1], b[0], b[1]))
print("\nThe canonical result: at N/T=0.8 the Bayesian portfolio reaches Sharpe %.2f, decisively above" % ffres["bayes"][1][0])
print("1/N (%.2f) and far above the sample optimiser (%.2f) -- with a fraction of the volatility."
      % (ffres["1/N"][1][0], ffres["sample"][1][0]))

Fama-French 48 industries, OOS annualised Sharpe (volatility %):

strategy         T=120mo (N/T=0.40)    T=60mo (N/T=0.80)
sample                   0.67 (  18%)           0.30 (  88%)

bayes                    1.00 (  12%)           1.01 (  12%)
1/N                      0.79 (  17%)           0.83 (  17%)

The canonical result: at N/T=0.8 the Bayesian portfolio reaches Sharpe 1.01, decisively above
1/N (0.83) and far above the sample optimiser (0.30) -- with a fraction of the volatility.

7. Summary and the bridge to Black–Litterman

• Bayesian estimation places a Normal-Inverse-Wishart prior on $(\mu,\Sigma)$; the posterior is closed-form and its mean is a shrinkage estimator — the rigorous version of Project 1.
• A from-scratch Metropolis sampler reproduces the analytic posterior, giving us a tool that survives when conjugacy is lost.
• The posterior predictive distribution inflates covariance by parameter uncertainty; optimising against it is Bayesian allocation.
• Estimation risk is real and large. The sample efficient frontier is a mirage — it promises risk/return it cannot deliver. Bayesian allocation cuts the opportunity cost of this error by up to ~99% when data is scarce, converging to the sample method only as $T\to\infty$.

The bridge. So far our prior was deliberately vague ($\mu_0=0$). But where should the prior come from? A powerful answer: let the market itself supply it. If we believe markets are roughly in equilibrium, the capitalisation-weighted portfolio implies a set of expected returns (via reverse optimisation). That equilibrium becomes the prior, and an investor's subjective views become the "data" that updates it. That is precisely the Black–Litterman model — Project 3 — the most famous Bayesian recipe in asset management, and a direct application of everything here.