Bayesian Estimation & Estimation Risk — the R engine¶
Risk and Asset Allocation¶
The independent R implementation of niw_python.ipynb. It re-derives the Normal-Inverse-Wishart posterior from scratch in base R, validates it against a from-scratch Gibbs sampler (a different algorithm — if both agree, the derivation and code are right), and reproduces the two estimation-risk results: the deception of the sample efficient frontier and the opportunity cost that Bayesian allocation reduces.
Self-contained; no prior reading required. Uses only base R + MASS (for the multivariate normal and rWishart), so there are no fragile package dependencies.
The model. Returns $x_t\sim\mathcal N(\mu,\Sigma)$; conjugate prior $\Sigma\sim\mathcal{IW}(\nu_0,\nu_0\Sigma_0)$, $\mu\mid\Sigma\sim\mathcal N(\mu_0,\Sigma/T_0)$. The posterior is again NIW with $$T_1=T_0+T,\quad \mu_1=\tfrac{T_0\mu_0+T\hat\mu}{T_0+T},\quad \nu_1=\nu_0+T,\quad \nu_1\Sigma_1=\nu_0\Sigma_0+T\hat\Sigma+\tfrac{T T_0}{T_0+T}(\mu_0-\hat\mu)(\mu_0-\hat\mu)'.$$
suppressMessages(library(MASS))
options(repr.plot.width=9, repr.plot.height=4.6)
BLUE<-"#2b6cb0"; ORANGE<-"#dd6b20"; GREEN<-"#2f855a"; RED<-"#c53030"; GREY<-"#718096"
R <- read.csv("sector_etf_weekly.csv", row.names=1)
SECTORS <- setdiff(colnames(R), "SPY"); X <- as.matrix(R[, SECTORS])
T <- nrow(X); N <- ncol(X)
cat(sprintf("Real panel: %d sector ETFs over %d weeks\n", N, T))
Real panel: 10 sector ETFs over 310 weeks
The real dataset: 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; SPY is the market benchmark. Estimation-risk claims are validated on a synthetic market whose true $(\mu,\Sigma)$ we set ourselves.
# --- core NIW machinery in base R ---
riw <- function(nu, Psi) { W <- rWishart(1, nu, solve(Psi))[,,1]; solve(W) } # InvWishart draw
sample_mom <- function(X){ mu<-colMeans(X); Xc<-sweep(X,2,mu); list(mu=mu, S=crossprod(Xc)/nrow(X)) }
niw_post <- function(X, mu0, T0, Sigma0, nu0){
Tn<-nrow(X); sm<-sample_mom(X)
T1<-T0+Tn; mu1<-(T0*mu0+Tn*sm$mu)/T1; nu1<-nu0+Tn
d<-matrix(mu0-sm$mu, ncol=1)
Sigma1<-(nu0*Sigma0 + Tn*sm$S + (Tn*T0/T1)*(d%*%t(d)))/nu1
list(mu1=mu1, T1=T1, Sigma1=Sigma1, nu1=nu1)
}
niw_draw <- function(p, n){ N<-length(p$mu1); mus<-matrix(0,n,N); Sg<-array(0,c(n,N,N))
for(i in 1:n){ S<-riw(p$nu1, p$nu1*p$Sigma1); Sg[i,,]<-S; mus[i,]<-mvrnorm(1,p$mu1,S/p$T1) }
list(mu=mus, Sigma=Sg) }
post_mean_cov <- function(p){ N<-length(p$mu1); list(mu=p$mu1, Sigma=p$nu1*p$Sigma1/(p$nu1-N-1)) }
predictive <- function(p){ N<-length(p$mu1); ES<-p$nu1*p$Sigma1/(p$nu1-N-1); list(mu=p$mu1, cov=(1+1/p$T1)*ES) }
# weak prior for the sector ETFs
mu0<-rep(0,N); Sigma0<-mean(diag(sample_mom(X)$S))*diag(N); T0<-5; nu0<-N+4
p <- niw_post(X, mu0, T0, Sigma0, nu0)
sm <- sample_mom(X)
cat(sprintf("prior weight on the mean = T0/(T0+T) = %.3f\n\n", T0/(T0+T)))
print(round(data.frame(sample=sm$mu, posterior=p$mu1), 3))
prior weight on the mean = T0/(T0+T) = 0.016
sample posterior 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
1. Validation: analytic NIW vs a from-scratch Gibbs sampler¶
The closed-form posterior exists only because the prior is conjugate. As an independent check (and a template for models where no formula exists) we build a Gibbs sampler that alternates the two exact full-conditionals $$\mu\mid\Sigma,X \sim \mathcal N\!\Big(\tfrac{T_0\mu_0+T\bar x}{T_0+T},\ \tfrac{\Sigma}{T_0+T}\Big),\qquad \Sigma\mid\mu,X \sim \mathcal{IW}\!\Big(\nu_0+T+1,\ \nu_0\Sigma_0+\textstyle\sum_t(x_t-\mu)(x_t-\mu)'+T_0(\mu-\mu_0)(\mu-\mu_0)'\Big).$$ If the sampler's output matches the analytic NIW draws, our derivation is confirmed. We use a small 2-asset market so the posterior is easy to visualise.
gibbs_niw <- function(X, mu0, T0, Sigma0, nu0, ndraw=8000, burn=1500){
Tn<-nrow(X); N<-ncol(X); xbar<-colMeans(X); mu<-xbar; Sig<-cov(X)
MU<-matrix(0,ndraw,N); SG<-array(0,c(ndraw,N,N)); T1<-T0+Tn
for(it in 1:(ndraw+burn)){
mu <- mvrnorm(1, (T0*mu0+Tn*xbar)/T1, Sig/T1)
Xc<-sweep(X,2,mu); dm<-matrix(mu-mu0,ncol=1)
Sig <- riw(nu0+Tn+1, nu0*Sigma0 + crossprod(Xc) + T0*(dm%*%t(dm)))
if(it>burn){ MU[it-burn,]<-mu; SG[it-burn,,]<-Sig }
}
list(mu=MU, Sigma=SG)
}
set.seed(1)
mu_t2<-c(1,0.7); S_t2<-matrix(c(4,1.2,1.2,2),2); X2<-mvrnorm(52, mu_t2, S_t2)
pri2<-list(mu0=c(0,0), T0=10, Sigma0=3*diag(2), nu0=10)
p2 <- niw_post(X2, pri2$mu0, pri2$T0, pri2$Sigma0, pri2$nu0)
da <- niw_draw(p2, 8000)
dg <- gibbs_niw(X2, pri2$mu0, pri2$T0, pri2$Sigma0, pri2$nu0, 8000, 1500)
par(mfrow=c(1,2))
hist(da$mu[,1], breaks=50, freq=FALSE, col=rgb(.44,.42,.46,.5), border=NA,
main=expression(paste("Posterior of ", mu[1])), xlab="")
hist(dg$mu[,1], breaks=50, freq=FALSE, col=rgb(.17,.42,.69,.5), border=NA, add=TRUE)
legend("topright", c("analytic NIW","Gibbs"), fill=c(GREY,BLUE), bty="n")
hist(da$Sigma[,1,2], breaks=50, freq=FALSE, col=rgb(.44,.42,.46,.5), border=NA,
main=expression(paste("Posterior of ", Sigma[12])), xlab="")
hist(dg$Sigma[,1,2], breaks=50, freq=FALSE, col=rgb(.18,.52,.36,.5), border=NA, add=TRUE)
legend("topright", c("analytic NIW","Gibbs"), fill=c(GREY,GREEN), bty="n")
par(mfrow=c(1,1))
cat(sprintf("mu mean analytic %s gibbs %s\n", paste(round(colMeans(da$mu),3),collapse=","), paste(round(colMeans(dg$mu),3),collapse=",")))
cat(sprintf("S12 mean analytic %.3f gibbs %.3f -> the two agree\n", mean(da$Sigma[,1,2]), mean(dg$Sigma[,1,2])))
mu mean analytic 0.729,0.434 gibbs 0.73,0.434
S12 mean analytic 0.826 gibbs 0.827 -> the two agree
2. The posterior predictive and its variance inflation¶
Integrating out $(\mu,\Sigma)$ gives the distribution of a future return; its covariance $\big(1+\tfrac1{T_1}\big)\tfrac{\nu_1\Sigma_1}{\nu_1-N-1}$ is inflated relative to the plug-in, by the amount of parameter uncertainty. Optimising against it is Bayesian allocation.
pred <- predictive(p)
cat("Predictive vs plug-in variance (a few sectors):\n")
for(i in c(1,5,10)) cat(sprintf(" %-4s plug-in %6.3f predictive %6.3f (x%.3f)\n",
SECTORS[i], sm$S[i,i], pred$cov[i,i], pred$cov[i,i]/sm$S[i,i]))
Predictive vs plug-in variance (a few sectors):
XLB plug-in 9.553 predictive 9.966 (x1.043) XLI plug-in 9.657 predictive 10.069 (x1.043) XLY plug-in 11.472 predictive 11.874 (x1.035)
3. Estimation risk: the deception of the sample frontier¶
The efficient frontier weights (budget-constrained, shorts allowed) are $w(m)=\Sigma^{-1}[\mathbf 1(C-Bm)+\mu(Am-B)]/D$ with $A=\mathbf 1'\Sigma^{-1}\mathbf 1,\ B=\mathbf 1'\Sigma^{-1}\mu,\ C=\mu'\Sigma^{-1}\mu,\ D=AC-B^2$. On a known synthetic market we compute the sample frontier and compare what it claims with what those portfolios truly deliver.
make_true <- function(N, rho=0.5, vlo=1, vhi=4, mu_scale=0.6, seed=3){
set.seed(seed); v<-seq(vlo,vhi,length.out=N); C<-(1-rho)*diag(N)+rho*matrix(1,N,N)
Sig<-outer(v,v)*C; list(mu=as.numeric(mu_scale*(Sig%*%rep(1,N))/N), Sigma=Sig) }
fw <- function(mu,Sig,m){ Si<-solve(Sig); one<-rep(1,length(mu))
A<-c(one%*%Si%*%one); B<-c(one%*%Si%*%mu); C<-c(mu%*%Si%*%mu); D<-A*C-B*B
as.numeric(Si%*%(one*(C-B*m)+mu*(A*m-B))/D) }
Na<-8; tm<-make_true(Na); pri<-list(mu0=rep(0,Na), T0=20, Sigma0=mean(diag(tm$Sigma))*diag(Na), nu0=20)
grid<-seq(min(tm$mu)-0.2, max(tm$mu)*1.6, length.out=25)
Tn<-24; nrep<-500
claimed<-matrix(0,25,2); trueS<-matrix(0,25,2)
set.seed(7)
for(r in 1:nrep){ Xi<-mvrnorm(Tn, tm$mu, tm$Sigma); sm2<-sample_mom(Xi)
for(i in seq_along(grid)){ w<-fw(sm2$mu, sm2$S, grid[i])
claimed[i,]<-claimed[i,]+c(sqrt(w%*%sm2$S%*%w), w%*%sm2$mu)
trueS[i,]<-trueS[i,]+c(sqrt(w%*%tm$Sigma%*%w), w%*%tm$mu) } }
claimed<-claimed/nrep; trueS<-trueS/nrep
tv<-sapply(grid, function(m){ w<-fw(tm$mu,tm$Sigma,m); sqrt(w%*%tm$Sigma%*%w) })
plot(tv, grid, type="l", col=GREEN, lwd=3, xlab="volatility (weekly %)", ylab="expected return (weekly %)",
main="The sample frontier is a mirage (N=8, T=24, 500 samples)", xlim=range(c(tv,claimed[,1],trueS[,1])))
lines(claimed[,1], claimed[,2], col=BLUE, lwd=2, lty=2)
points(trueS[,1], trueS[,2], col=RED, pch=19, cex=.6); lines(trueS[,1], trueS[,2], col=RED, lwd=2)
legend("bottomright", c("TRUE frontier (ideal)","CLAIMED (sample shows)","TRUE performance of sample portfolios"),
col=c(GREEN,BLUE,RED), lwd=2, lty=c(1,2,1), bty="n")
cat("Blue (claimed) sits above-left of green (truth); red (reality) sits below -- the deception.\n")
Blue (claimed) sits above-left of green (truth); red (reality) sits below -- the deception.
4. Bayesian allocation reduces the opportunity cost¶
For risk aversion $\gamma$ the ideal investor holds $w^\star=\arg\max_{w'\mathbf 1=1} w'\mu-\tfrac{\gamma}{2}w'\Sigma w$. A real investor estimates the inputs; the opportunity cost $\text{CE}_{\text{true}}-\text{CE}_{\text{realised}}$ is the price of estimation error. We compare sample vs Bayesian-predictive inputs across sample sizes.
opt_mv <- function(mu,Sig,gamma){ Si<-solve(Sig); one<-rep(1,length(mu))
A<-c(one%*%Si%*%one); B<-c(one%*%Si%*%mu); lam<-(B-gamma)/A; as.numeric(Si%*%(mu-lam*one)/gamma) }
ce <- function(w,mu,Sig,gamma) as.numeric(w%*%mu - 0.5*gamma*(w%*%Sig%*%w))
gamma<-0.5; w_star<-opt_mv(tm$mu, tm$Sigma, gamma); ce_true<-ce(w_star, tm$mu, tm$Sigma, gamma)
Ts<-c(16,24,40,80,160,320); oc_s<-oc_b<-numeric(length(Ts)); set.seed(11)
for(k in seq_along(Ts)){ Tn<-Ts[k]; s<-b<-0; nrep<-800
for(r in 1:nrep){ Xi<-mvrnorm(Tn, tm$mu, tm$Sigma); smi<-sample_mom(Xi)
pp<-niw_post(Xi, pri$mu0, pri$T0, pri$Sigma0, pri$nu0); pr<-predictive(pp)
ws<-opt_mv(smi$mu, smi$S, gamma); wb<-opt_mv(pr$mu, pr$cov, gamma)
s<-s+(ce_true-ce(ws,tm$mu,tm$Sigma,gamma)); b<-b+(ce_true-ce(wb,tm$mu,tm$Sigma,gamma)) }
oc_s[k]<-s/nrep; oc_b[k]<-b/nrep }
par(mfrow=c(1,2))
plot(Ts, oc_s, type="b", pch=15, col=RED, lwd=2, log="y", xlab="sample size T",
ylab="opportunity cost (log)", main="Cost of estimation error")
lines(Ts, oc_b, type="b", pch=19, col=BLUE, lwd=2); abline(h=ce_true, col=GREY, lty=3)
legend("topright", c("sample (plug-in)","Bayesian (predictive)"), col=c(RED,BLUE), pch=c(15,19), lwd=2, bty="n")
plot(Ts, 100*(1-oc_b/oc_s), type="b", pch=19, col=GREEN, lwd=2, ylim=c(0,100),
xlab="sample size T", ylab="reduction (%)", main="Bayesian reduction in opportunity cost")
par(mfrow=c(1,1))
for(k in seq_along(Ts)) cat(sprintf("T=%3d sample OC %8.3f bayes OC %6.3f (%.0f%% lower)\n",
Ts[k], oc_s[k], oc_b[k], 100*(1-oc_b[k]/oc_s[k])))
T= 16 sample OC 14.628 bayes OC 0.173 (99% lower) T= 24 sample OC 3.166 bayes OC 0.166 (95% lower) T= 40 sample OC 1.059 bayes OC 0.144 (86% lower) T= 80 sample OC 0.299 bayes OC 0.101 (66% lower) T=160 sample OC 0.125 bayes OC 0.067 (47% lower) T=320 sample OC 0.057 bayes OC 0.040 (29% lower)
4b. Dimensionality: the estimation-risk tax scales with N/T¶
The opportunity-cost above used a small market with a short sample. The same severity arrives naturally through dimensionality. We confirm on 48 stocks (weekly) and the canonical Fama–French 48 industries (monthly, DeMiguel–Garlappi–Uppal): a rolling out-of-sample backtest of the sample vs Bayesian mean-variance portfolio against the estimation-error-free 1/N.
oosR <- function(Xd, window, which, ppy, g=0.5){
N <- ncol(Xd); rets <- c()
for (t in (window+1):nrow(Xd)) {
Xw <- Xd[(t-window):(t-1), , drop=FALSE]
if (which == "1/N") { w <- rep(1/N, N) }
else {
if (which == "sample") { sm <- sample_mom(Xw); mu <- sm$mu; S <- sm$S }
else { p <- niw_post(Xw, rep(0,N), 5, mean(diag(sample_mom(Xw)$S))*diag(N), N+4)
pr <- predictive(p); mu <- pr$mu; S <- pr$cov }
w <- opt_mv(mu, S, g)
}
rets <- c(rets, sum(Xd[t,]*w))
}
c(mean(rets)/sd(rets)*sqrt(ppy), sd(rets)*sqrt(ppy))
}
Xstk <- as.matrix(read.csv("stocks_weekly.csv", row.names=1))
FF <- as.matrix(read.csv("ff_industries_monthly.csv", row.names=1))
cat("OOS annualised Sharpe (volatility %):\n\n")
cat(sprintf("%-13s %18s %18s\n", "strategy", "48 stocks T=60wk", "FF 48 ind T=60mo"))
for (which in c("sample","bayes","1/N")) {
a <- oosR(Xstk, 60, which, 52); b <- oosR(FF, 60, which, 12)
cat(sprintf("%-13s %11.2f (%4.0f%%) %11.2f (%4.0f%%)\n", which, a[1], a[2], b[1], b[2]))
}
cat("\nAt high N/T the sample optimiser deteriorates while the Bayesian portfolio overtakes 1/N on\n")
cat("both universes -- most emphatically on the canonical Fama-French benchmark (Sharpe ~1.0 vs ~0.83).\n")
OOS annualised Sharpe (volatility %):
strategy 48 stocks T=60wk FF 48 ind T=60mo
sample 0.43 ( 199%) 0.30 ( 88%) bayes 0.94 ( 16%) 1.01 ( 12%) 1/N 0.77 ( 20%) 0.83 ( 17%)
At high N/T the sample optimiser deteriorates while the Bayesian portfolio overtakes 1/N on
both universes -- most emphatically on the canonical Fama-French benchmark (Sharpe ~1.0 vs ~0.83).
5. Summary¶
- The from-scratch NIW posterior matches a from-scratch Gibbs sampler on both the mean and the covariance — two independent routes, same answer.
- The posterior predictive inflates risk by parameter uncertainty.
- The sample efficient frontier is a mirage: it claims risk/return it cannot deliver.
- Bayesian (predictive) inputs cut the opportunity cost of estimation error — dramatically when data is scarce, converging to the sample method only as $T\to\infty$.
This matches the Python and PyMC engines exactly. Project 3 (Black–Litterman) answers the natural next question: where should the prior come from? Answer — the market equilibrium, updated by the investor's views.