Bayesian Multivariate Regression — from-scratch conjugate sampler¶
mvr_conjugate.py · the tuna demand system¶
Model. $M$ outcomes share one design matrix: $$Y = XB + E,\qquad E_t \sim N_M(0,\ \Sigma),\qquad B:\ K\times M,\ \ \Sigma:\ M\times M.$$ Every equation uses the same regressors $X$ — that is what defines multivariate regression. The payoff over running $M$ separate regressions is $\Sigma$: the cross-equation error correlation, estimated jointly.
Conjugate posterior (Normal–Inverse-Wishart). With $B\mid\Sigma\sim MN(\bar B,\Sigma\otimes A^{-1})$ and $\Sigma\sim IW(\nu,V)$, the posterior is again Normal–IW, so we draw i.i.d. samples (no MCMC) via rmultireg_np_batch.
Real example — tuna demand: $N=338$ store-weeks, $M=7$ canned-tuna brands. Each brand's log unit-sales on all 7 log-prices + 7 promotion indicators (same $X$). The price coefficients are own- and cross-price elasticities; $\Sigma$ is the cross-brand demand-shock correlation.
import numpy as np, pandas as pd, matplotlib.pyplot as plt
from mvr_conjugate import rmultireg_np_batch, posterior_summary
np.set_printoptions(suppress=True)
1. Synthetic validation¶
Recover a known $(B,\Sigma)$ and check the conjugate draws match OLS (under a diffuse prior) and cover the truth.
Ys = pd.read_csv('mvr_Y.csv').values.astype(float); Xs = pd.read_csv('mvr_X.csv').values.astype(float)
Ns, Ms = Ys.shape; Ks = Xs.shape[1]
trueB = pd.read_csv('mvr_true_B.csv', index_col=0).values
if trueB.shape != (Ks, Ms): trueB = trueB.T
ols = np.linalg.lstsq(Xs, Ys, rcond=None)[0]
Bd, Sd = rmultireg_np_batch(Ys, Xs, np.zeros((Ks,Ms)), 0.01*np.eye(Ks), Ms+3, (Ms+2)*np.eye(Ms), R=10000, seed=1)
summ = posterior_summary(Bd, Ks, Ms, B_true=trueB, ols_B=ols)
print('B: max |posterior mean - OLS| = %.4f (diffuse prior -> should match OLS)' % (summ['bayes_mean'].values - summ['ols'].values).__abs__().max())
print('95%% credible-interval coverage of the true B: %d / %d' % (summ['covers'].sum(), len(summ)))
print('posterior-mean Sigma vs truth:')
print(np.round(Sd.mean(0).reshape(Ms,Ms,order='F'),3)); print('truth:'); print(pd.read_csv('mvr_true_Sigma.csv', index_col=0).values.round(3))
B: max |posterior mean - OLS| = 0.0010 (diffuse prior -> should match OLS) 95% credible-interval coverage of the true B: 12 / 12 posterior-mean Sigma vs truth: [[1.357 0.477 0.067] [0.477 0.809 0.329] [0.067 0.329 0.557]] truth: [[1.45 0.63 0.21] [0.63 0.89 0.35] [0.21 0.35 0.49]]
2. Tuna demand system (real data)¶
The data. tuna_lmove.csv + tuna_X.csv are store-level scanner records from Dominick's Finer Foods (Chicago) for the canned-tuna category — the bayesm tuna dataset (Chevalier, Kashyap & Rossi 2003). N = 338 store-weeks, M = 7 brands:
| # | brand | note |
|---|---|---|
| 1 | StarKist | |
| 2 | Chicken of the Sea (ChickSea) |
|
| 3 | Bumble Bee Solid (BB Solid) |
← Bumble Bee |
| 4 | Bumble Bee Chunk (BB Chunk) |
← Bumble Bee |
| 5 | Geisha | |
| 6 | Bumble Bee Large (BB Large) |
← Bumble Bee |
| 7 | H-H Chunk (HH Chunk) |
store/other |
$Y$ = log unit-sales ("log move"); $X$ = intercept + the 7 log-prices (lp1..lp7) + 7 feature/display promo flags (ns1..ns7) — the same regressors in every equation. The price coefficients are own- and cross-price elasticities; $\Sigma$ is the cross-brand demand-shock correlation. Watch brands 3, 4 and 6 — all Bumble Bee variants. Diffuse Normal–IW prior; 20,000 i.i.d. conjugate draws.
Y = pd.read_csv('tuna_lmove.csv').values.astype(float)
X = pd.read_csv('tuna_X.csv').values.astype(float)
N, M = Y.shape; K = X.shape[1]; brands = ['StarKist','ChickSea','BB Solid','BB Chunk','Geisha','BB Large','HH Chunk']
print('Tuna: N=%d store-weeks, M=%d brands, K=%d regressors' % (N, M, K))
Bd, Sd = rmultireg_np_batch(Y, X, np.zeros((K,M)), 0.01*np.eye(K), M+3, (M+2)*np.eye(M), R=20000, seed=42)
Bmat = Bd.mean(0).reshape((K, M), order='F') # (K, M) posterior-mean coefficients
E_cross = Bmat[1:1+M, :].T # (sales eq j, price brand k): rows lp1..lp7
Sig = np.array([s.reshape(M, M, order='F') for s in Sd]) # (R, M, M)
dd = np.sqrt(np.diagonal(Sig, axis1=1, axis2=2)) # (R, M)
Corr = Sig / (dd[:, :, None] * dd[:, None, :]) # (R, M, M) per-draw correlation
Corr_mean = Corr.mean(0)
print('own-price elasticities (diagonal):', np.round(np.diag(E_cross), 2))
Tuna: N=338 store-weeks, M=7 brands, K=15 regressors
own-price elasticities (diagonal): [-4.29 -4.66 -5.36 -4.87 -4.02 3.11 -2.36]
Graphs 1 & 2 — cross-equation error correlation and the elasticity matrix¶
Left: posterior-mean correlation of the demand shocks across brands — why we model the equations jointly. Right: the price-elasticity matrix — diagonal = own-price (negative), off-diagonal = cross-price (positive ⇒ substitutes).
fig, ax = plt.subplots(1, 2, figsize=(13.5, 5.5))
im0 = ax[0].imshow(Corr_mean, cmap='RdBu_r', vmin=-1, vmax=1)
for i in range(M):
for j in range(M):
ax[0].text(j, i, f'{Corr_mean[i,j]:.2f}', ha='center', va='center', fontsize=7,
color='white' if abs(Corr_mean[i,j])>0.5 else 'black')
ax[0].set_xticks(range(M)); ax[0].set_xticklabels(brands, rotation=45, ha='right'); ax[0].set_yticks(range(M)); ax[0].set_yticklabels(brands)
ax[0].set_title('Σ error correlation across brands\n(off-diagonal ≠ 0 ⇒ joint modeling pays)'); plt.colorbar(im0, ax=ax[0], fraction=.046)
vmax = np.abs(E_cross).max()
im1 = ax[1].imshow(E_cross, cmap='RdBu_r', vmin=-vmax, vmax=vmax)
for i in range(M):
for j in range(M):
ax[1].text(j, i, f'{E_cross[i,j]:.2f}', ha='center', va='center', fontsize=7,
color='white' if abs(E_cross[i,j])>vmax*0.5 else 'black')
ax[1].set_xticks(range(M)); ax[1].set_xticklabels([f'P:{b}' for b in brands], rotation=45, ha='right')
ax[1].set_yticks(range(M)); ax[1].set_yticklabels([f'sales:{b}' for b in brands])
ax[1].set_title('Price-elasticity matrix\n(diagonal = own-price, off-diagonal = cross-price)'); plt.colorbar(im1, ax=ax[1], fraction=.046)
plt.tight_layout(); plt.savefig('mvr_heatmaps.png', dpi=120, bbox_inches='tight'); plt.show()
Graphs 3 & 4 — own-price elasticities and the strongest cross-equation correlation¶
Left: each brand's own-price elasticity with its 95% credible interval. Right: the full posterior of the largest off-diagonal error correlation — is it really $\ne 0$?
fig, ax = plt.subplots(1, 2, figsize=(13.5, 5))
own = np.array([Bd[:, j*K + (1+j)] for j in range(M)]) # own-price draws per brand: B[lp_j, eq_j]
mean = own.mean(1); lo = np.percentile(own, 2.5, 1); hi = np.percentile(own, 97.5, 1)
yy = np.arange(M)
ax[0].errorbar(mean, yy, xerr=[mean-lo, hi-mean], fmt='o', color='firebrick', capsize=4)
ax[0].axvline(0, color='gray', ls=':'); ax[0].axvline(-1, color='steelblue', ls='--', lw=1, label='unit elastic (−1)')
ax[0].set_yticks(yy); ax[0].set_yticklabels(brands); ax[0].invert_yaxis()
ax[0].set_xlabel('own-price elasticity'); ax[0].set_title('Own-price elasticities (95% CrI)'); ax[0].legend(fontsize=8)
# strongest off-diagonal correlation
iu, ju = np.triu_indices(M, 1); k = np.argmax(np.abs(Corr_mean[iu, ju])); i_, j_ = iu[k], ju[k]
cd = Corr[:, i_, j_]
ax[1].hist(cd, bins=50, color='steelblue', edgecolor='white', density=True)
ax[1].axvline(cd.mean(), color='firebrick', lw=2, label='mean %.2f'%cd.mean()); ax[1].axvline(0, color='black', ls=':')
ax[1].set_xlabel('error correlation corr(%s, %s)'%(brands[i_], brands[j_])); ax[1].set_yticks([])
ax[1].set_title('Strongest cross-equation correlation\n95%% CrI [%.2f, %.2f] — excludes 0'%(np.percentile(cd,2.5), np.percentile(cd,97.5))); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.savefig('mvr_elasticities.png', dpi=120, bbox_inches='tight'); plt.show()
Analysis of the findings¶
1. Brand-level tuna demand is highly price-elastic. Own-price elasticities (the diagonal) run about −4 to −5 for six of the seven brands (StarKist −4.3, Chicken of the Sea −4.7, BB Solid −5.4, BB Chunk −4.9, Geisha −4.0, H-H Chunk −2.4). A 1% own-price rise cuts a brand's unit sales by ~4–5% — consumers switch readily among near-substitutable canned-tuna brands. H-H Chunk is the least elastic (−2.4), consistent with a cheaper store/secondary label whose buyers are less price-responsive.
2. The Bumble Bee entanglement is the story. The single largest cross-equation error correlation is corr(BB Solid, BB Large) = 0.61 (95% CrI [0.54, 0.68], far from 0), and that same pair drives the only economically odd coefficient: BB Large's positive own-price elasticity (+3.1). Both are explained by the brands being two Bumble Bee products: their shelf prices are set together (so the price regressors are collinear → BB Large's own-price coefficient is weakly identified and flips sign, while BB Solid loads a large cross-price term on BB Large), and their demand shocks move together (shared promotions, displays, and brand-level demand swings → the 0.61 error correlation). This is the clean cautionary tale of full cross-price systems: for closely-related products the individual elasticities can be unreliable even when the joint fit is fine.
3. Brands are substitutes. Away from the Bumble Bee pair, cross-price elasticities are mostly positive — a rival's price increase lifts your sales — the expected substitution pattern for a commodity-like category.
4. Why model jointly (the point of MVR). With identical regressors the coefficient point estimates equal equation-by-equation OLS, so MVR doesn't sharpen the coefficients. What it adds is Σ — and here Σ is not diagonal: apart from the strong Bumble Bee link, the other off-diagonals are modest (|corr| ≲ 0.3). So for most brands separate regressions would lose little, but the Bumble Bee co-movement is real and only the joint model reveals it — giving coherent joint/predictive inference and honest uncertainty on the whole 7×7 elasticity matrix.
Cross-checks: PyMC (mvr_pymc.ipynb) and R bayesm (mvr_bayesm.ipynb) reproduce all of the above.