import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from hsur_gibbs import rhsur_gibbs, simulate_hsur, posterior_summary

np.set_printoptions(precision=4, suppress=True)
pd.set_option('display.float_format', '{:.4f}'.format)
print(f'NumPy {np.__version__}')

NumPy 2.4.6

nreg_s = 50
M_s    = 3
k_s    = 2   # intercept + slope
T_s    = 30
K_s    = M_s * k_s   # 6

# True parameters
true_Delta_s = np.array([[3.0, -2.0, 2.0, -1.5, 1.0, -1.0]])  # (nz=1, K=6)
true_Vbeta_s = 0.25 * np.eye(K_s)
true_Omega_s = np.array([[1.0, 0.5, 0.2],
                          [0.5, 1.0, 0.3],
                          [0.2, 0.3, 0.8]])

print('True Delta (population means per coefficient):')
eq_names  = [f'eq{j+1}' for j in range(M_s)]
par_names = ['intercept', 'slope']
col_names = [f'{e}.{p}' for e in eq_names for p in par_names]
print(pd.DataFrame(true_Delta_s, columns=col_names, index=['intercept-only z']))
print('\nTrue Omega:')
print(pd.DataFrame(true_Omega_s, index=eq_names, columns=eq_names))

True Delta (population means per coefficient):
                  eq1.intercept  eq1.slope  eq2.intercept  eq2.slope  \
intercept-only z         3.0000    -2.0000         2.0000    -1.5000   

                  eq3.intercept  eq3.slope  
intercept-only z         1.0000    -1.0000  

True Omega:
       eq1    eq2    eq3
eq1 1.0000 0.5000 0.2000
eq2 0.5000 1.0000 0.3000
eq3 0.2000 0.3000 0.8000

regdata_s, betas_s, Z_s = simulate_hsur(
    nreg_s, M_s, k_s, T_s,
    true_Delta_s, true_Vbeta_s, true_Omega_s,
    seed=0
)

print(f'Simulated: {nreg_s} stores × {M_s} equations × {T_s} periods')
print(f'Total obs: {nreg_s * T_s}   K = {K_s} coefficients per store')
print(f'\nStore-level beta spread (std across stores):')
print(pd.Series(betas_s.std(0).round(3), index=col_names, name='std'))

Simulated: 50 stores × 3 equations × 30 periods
Total obs: 1500   K = 6 coefficients per store

Store-level beta spread (std across stores):
eq1.intercept   0.5310
eq1.slope       0.4790
eq2.intercept   0.4920
eq2.slope       0.5520
eq3.intercept   0.4760
eq3.slope       0.4460
Name: std, dtype: float64

prior_s = dict(
    nu      = M_s + 3,
    V       = (M_s + 2.0) * np.eye(M_s),
    nu_beta = K_s + 3,
    V_beta  = (K_s + 2.0) * np.eye(K_s)
)
mcmc_s = dict(R=10000, burn=2000, keep=1, nprint=2000)

t0     = time.perf_counter()
out_s  = rhsur_gibbs(regdata_s, Z=Z_s, Prior=prior_s, Mcmc=mcmc_s, seed=42)
print(f'\nDone in {time.perf_counter()-t0:.1f}s  |  kept draws: {out_s["betadraw"].shape[2]}')

  Iter   2000/10000

  Iter   4000/10000

  Iter   6000/10000

  Iter   8000/10000

  Iter  10000/10000

Done in 31.6s  |  kept draws: 8000

Ddraw_s = out_s['Deltadraw']   # (R_kept, 6)  -- nz=1 so vec(Delta) = Delta[0,:]

true_delta_flat = true_Delta_s.ravel(order='F')  # (6,)
summ_D = posterior_summary(Ddraw_s, col_names, true_vals=true_delta_flat)

print('Delta (population means) — recovery:')
print(summ_D.round(4))
n_cov = summ_D['covered'].sum() if 'covered' in summ_D else None
if n_cov is not None:
    print(f'\nCoverage: {n_cov}/{len(summ_D)} parameters covered at 95% CrI')

Delta (population means) — recovery:
                 mean     sd    q025    q975    true  covered
param                                                        
eq1.intercept  2.9668 0.1033  2.7636  3.1697  3.0000     True
eq1.slope     -1.8726 0.0937 -2.0550 -1.6883 -2.0000     True
eq2.intercept  1.9103 0.0961  1.7222  2.1009  2.0000     True
eq2.slope     -1.6201 0.0981 -1.8138 -1.4280 -1.5000     True
eq3.intercept  0.8671 0.0956  0.6844  1.0528  1.0000     True
eq3.slope     -0.9539 0.0928 -1.1331 -0.7691 -1.0000     True

Coverage: 6/6 parameters covered at 95% CrI

Odraw_s  = out_s['Omegadraw']   # (R_kept, 9)
Omega_pm_s = Odraw_s.mean(0).reshape(M_s, M_s, order='F')

print('Omega recovery (posterior mean vs true):')
tbl = []
for i in range(M_s):
    for j in range(i, M_s):
        d   = Odraw_s[:, j * M_s + i]   # col-major vec: Omega[i,j] at index j*M+i
        lo, hi = np.percentile(d, [2.5, 97.5])
        tr  = true_Omega_s[i, j]
        tbl.append({'entry': f'Omega[{i+1},{j+1}]',
                    'true': tr, 'mean': d.mean(), 'q025': lo, 'q975': hi,
                    'covered': lo <= tr <= hi})
print(pd.DataFrame(tbl).set_index('entry').round(4))

print('\nVbeta diagonal (posterior mean std of each coefficient across stores):')
Vbdraw_s = out_s['Vbetadraw']  # (R_kept, 36)
Vb_pm_s  = Vbdraw_s.mean(0).reshape(K_s, K_s, order='F')
print(pd.Series(np.diag(Vb_pm_s).round(4), index=col_names, name='Vbeta_diag'))
print(f'True diagonal: {np.diag(true_Vbeta_s)}')

Omega recovery (posterior mean vs true):
             true   mean   q025   q975  covered
entry                                          
Omega[1,1] 1.0000 1.0413 0.9668 1.1225     True
Omega[1,2] 0.5000 0.5044 0.4472 0.5658     True
Omega[1,3] 0.2000 0.2048 0.1553 0.2534     True
Omega[2,2] 1.0000 0.9638 0.8957 1.0374     True
Omega[2,3] 0.3000 0.2705 0.2228 0.3199     True
Omega[3,3] 0.8000 0.7940 0.7386 0.8554     True

Vbeta diagonal (posterior mean std of each coefficient across stores):
eq1.intercept   0.4754
eq1.slope       0.3956
eq2.intercept   0.4246
eq2.slope       0.4497
eq3.intercept   0.4233
eq3.slope       0.3926
Name: Vbeta_diag, dtype: float64
True diagonal: [0.25 0.25 0.25 0.25 0.25 0.25]

# Explanatory: does hSUR recover the known synthetic truth?
fig, ax = plt.subplots(1, 3, figsize=(15, 4.2))
dm = Ddraw_s.mean(0); dlo, dhi = np.percentile(Ddraw_s, [2.5, 97.5], axis=0)
ax[0].errorbar(true_delta_flat, dm, yerr=[dm - dlo, dhi - dm], fmt='o', color='steelblue', capsize=3)
lim = [min(true_delta_flat.min(), dm.min()) - .5, max(true_delta_flat.max(), dm.max()) + .5]
ax[0].plot(lim, lim, 'k--', lw=.8); ax[0].set_xlim(lim); ax[0].set_ylim(lim)
ax[0].set_xlabel('true $\\Delta$'); ax[0].set_ylabel('posterior mean (95% CrI)'); ax[0].set_title('Population means $\\Delta$')
tt, em, elo, ehi = [], [], [], []
for i in range(M_s):
    for j in range(i, M_s):
        d = Odraw_s[:, j * M_s + i]; tt.append(true_Omega_s[i, j]); em.append(d.mean())
        lo, hi = np.percentile(d, [2.5, 97.5]); elo.append(lo); ehi.append(hi)
tt = np.array(tt); em = np.array(em)
ax[1].errorbar(tt, em, yerr=[em - np.array(elo), np.array(ehi) - em], fmt='o', color='firebrick', capsize=3)
lim = [min(tt.min(), em.min()) - .2, max(tt.max(), em.max()) + .2]; ax[1].plot(lim, lim, 'k--', lw=.8)
ax[1].set_xlabel('true $\\Omega$ entry'); ax[1].set_ylabel('posterior mean'); ax[1].set_title('Residual covariance $\\Omega$')
bts = out_s['betadraw'].mean(2)
ax[2].scatter(betas_s.ravel(), bts.ravel(), s=10, alpha=.4, color='seagreen')
lim = [min(betas_s.min(), bts.min()) - .5, max(betas_s.max(), bts.max()) + .5]; ax[2].plot(lim, lim, 'k--', lw=.8)
ax[2].set_xlabel('true store $\\beta$'); ax[2].set_ylabel('posterior mean store $\\beta$'); ax[2].set_title('Store-level $\\beta$ (50×6)')
plt.suptitle('Synthetic validation — hSUR recovers known $\\Delta$, $\\Omega$, and store-level $\\beta$', y=1.02)
plt.tight_layout(); plt.savefig('hsur_synth_recovery.png', dpi=120, bbox_inches='tight'); plt.show()

yx   = pd.read_csv('oj_yx.csv')
demo = pd.read_csv('oj_storedemo.csv')

M_oj = 11
k_oj = 3   # intercept, log-own-price, deal
K_oj = M_oj * k_oj   # 33

BRANDS_OJ = [
    'Trop64', 'Trop96', 'FloNat', 'Trop64b', 'MM64',
    'MM96',   'CitHill', 'TreeFr', 'FloGold', 'Dom64', 'Dom96'
]

stores = sorted(yx['store'].unique())
nreg   = len(stores)
Ts     = [yx[yx['store'] == s]['week'].nunique() for s in stores]

print(f'Stores: {nreg}   Brands: {M_oj}   k: {k_oj}   K: {K_oj}')
print(f'T per store: min={min(Ts)}  median={int(np.median(Ts))}  max={max(Ts)}')
print(f'Total obs (store-brand-week): {sum(T * M_oj for T in Ts):,}')
print(f'\nLog-sales summary (across all store-brand-weeks):')
print(yx['logmove'].describe().round(3))

Stores: 83   Brands: 11   k: 3   K: 33
T per store: min=87  median=117  max=121
Total obs (store-brand-week): 106,139

Log-sales summary (across all store-brand-weeks):
count   106139.0000
mean         8.4170
std          1.1410
min          4.1590
25%          7.7000
50%          8.3930
75%          9.0490
max         13.4820
Name: logmove, dtype: float64

regdata_oj = []

for store_id in stores:
    store_df = yx[yx['store'] == store_id]
    y_list = []
    X_list = []
    for j in range(M_oj):
        brand_id  = j + 1
        brand_df  = store_df[store_df['brand'] == brand_id].sort_values('week')
        y_j       = brand_df['logmove'].values.astype(float)
        log_p_j   = np.log(brand_df[f'price{brand_id}'].values.astype(float))
        deal_j    = brand_df['deal'].values.astype(float)
        X_j       = np.column_stack([np.ones(len(y_j)), log_p_j, deal_j])
        y_list.append(y_j)
        X_list.append(X_j)
    regdata_oj.append({'y_list': y_list, 'X_list': X_list})

print(f'regdata built: {len(regdata_oj)} stores')
print(f'Store 0: T={len(regdata_oj[0]["y_list"][0])}  X shape={regdata_oj[0]["X_list"][0].shape}')

regdata built: 83 stores
Store 0: T=110  X shape=(110, 3)

# Explanatory: the OJ demand relationship (log quantity vs log price), coloured by deal
fig, axes = plt.subplots(1, 2, figsize=(13, 4.6))
for ax, j in zip(axes, [0, 10]):
    lp, lm, dl = [], [], []
    for st in regdata_oj:
        X = st['X_list'][j]; y = st['y_list'][j]
        lp.append(X[:, 1]); lm.append(y); dl.append(X[:, 2])
    lp = np.concatenate(lp); lm = np.concatenate(lm); dl = np.concatenate(dl)
    sc = ax.scatter(lp, lm, c=dl, cmap='coolwarm', s=6, alpha=.35)
    b = np.polyfit(lp, lm, 1); xs = np.array([lp.min(), lp.max()])
    ax.plot(xs, b[0] * xs + b[1], 'k-', lw=2)
    ax.set_title(f'{BRANDS_OJ[j]}: demand curve (OLS slope = {b[0]:.2f})')
    ax.set_xlabel('log price'); ax.set_ylabel('log quantity (logmove)')
fig.colorbar(sc, ax=axes[1], label='deal (0 = no, 1 = yes)')
plt.suptitle('What hSUR fits — a price/quantity demand curve per brand, lifted by promotion', y=1.01)
plt.tight_layout(); plt.savefig('hsur_demand_eda.png', dpi=120, bbox_inches='tight'); plt.show()

Z_oj = np.ones((nreg, 1))   # intercept-only: common population mean

prior_oj = dict(
    nu      = M_oj + 3,
    V       = (M_oj + 2.0) * np.eye(M_oj),
    nu_beta = K_oj + 3,
    V_beta  = (K_oj + 2.0) * np.eye(K_oj)
)
mcmc_oj = dict(R=5000, burn=1000, keep=1, nprint=1000)

t0     = time.perf_counter()
out_oj = rhsur_gibbs(regdata_oj, Z=Z_oj, Prior=prior_oj, Mcmc=mcmc_oj, seed=42)
print(f'\nDone in {time.perf_counter()-t0:.1f}s  |  kept draws: {out_oj["betadraw"].shape[2]}')

  Iter   1000/5000

  Iter   2000/5000

  Iter   3000/5000

  Iter   4000/5000

  Iter   5000/5000

Done in 166.7s  |  kept draws: 4000

price_idx = [j * k_oj + 1 for j in range(M_oj)]   # 1, 4, 7, ..., 31
deal_idx  = [j * k_oj + 2 for j in range(M_oj)]   # 2, 5, 8, ..., 32

Ddraw_oj = out_oj['Deltadraw']   # (R_kept, 33)

# Parameter index layout: [int_1, price_1, deal_1, int_2, price_2, deal_2, ...]
price_idx = [j * k_oj + 1 for j in range(M_oj)]   # 1,4,7,...,31
deal_idx  = [j * k_oj + 2 for j in range(M_oj)]   # 2,5,8,...,32

tbl = []
for j in range(M_oj):
    p_draws = Ddraw_oj[:, price_idx[j]]
    d_draws = Ddraw_oj[:, deal_idx[j]]
    p_lo, p_hi = np.percentile(p_draws, [2.5, 97.5])
    d_lo, d_hi = np.percentile(d_draws, [2.5, 97.5])
    tbl.append({
        'brand'    : BRANDS_OJ[j],
        'price_mean': p_draws.mean(),
        'price_sd'  : p_draws.std(),
        'price_lo'  : p_lo,
        'price_hi'  : p_hi,
        'deal_mean' : d_draws.mean(),
        'deal_sd'   : d_draws.std(),
    })

df_pop = pd.DataFrame(tbl).set_index('brand')
print('Population-level own-price elasticities and deal effects:')
print(df_pop.round(4))

Population-level own-price elasticities and deal effects:
         price_mean  price_sd  price_lo  price_hi  deal_mean  deal_sd
brand                                                                
Trop64      -2.7992    0.0894   -2.9735   -2.6253     0.0660   0.0739
Trop96      -1.9969    0.0876   -2.1656   -1.8241     0.1894   0.0727
FloNat      -3.6726    0.1096   -3.8834   -3.4556     0.2456   0.0735
Trop64b     -3.3736    0.1002   -3.5709   -3.1762     0.6546   0.0774
MM64        -3.1769    0.0946   -3.3602   -2.9923     0.4197   0.0748
MM96        -2.4653    0.0921   -2.6472   -2.2830     0.1159   0.0711
CitHill     -4.3581    0.1078   -4.5691   -4.1458     0.1940   0.0742
TreeFr      -2.8157    0.0996   -3.0147   -2.6232     0.0665   0.0738
FloGold     -4.0570    0.1232   -4.2983   -3.8144     0.4026   0.0771
Dom64       -2.9569    0.0907   -3.1364   -2.7768     0.5203   0.0783
Dom96       -1.5009    0.0978   -1.6862   -1.3105     0.3652   0.0721

fig, ax = plt.subplots(figsize=(10, 5))
y_pos = np.arange(M_oj)
means = df_pop['price_mean'].values
lo    = df_pop['price_lo'].values
hi    = df_pop['price_hi'].values

ax.barh(y_pos, means, xerr=[means - lo, hi - means],
        align='center', color='steelblue', alpha=0.7,
        error_kw={'linewidth': 1.5, 'capsize': 4})
ax.axvline(0, color='black', linewidth=0.8, linestyle='--')
ax.set_yticks(y_pos)
ax.set_yticklabels(BRANDS_OJ)
ax.set_xlabel('Posterior mean (95% CrI)')
ax.set_title('Population-level own-price elasticities\n(hSUR, intercept-only hierarchy)')
plt.tight_layout()
plt.show()

Odraw_oj = out_oj['Omegadraw']   # (R_kept, 121)
Omega_pm = Odraw_oj.mean(0).reshape(M_oj, M_oj, order='F')

D_inv = np.diag(1.0 / np.sqrt(np.diag(Omega_pm)))
Corr  = D_inv @ Omega_pm @ D_inv

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

im0 = axes[0].imshow(Omega_pm, cmap='RdBu_r')
axes[0].set_title('Posterior mean $\\Omega$ (residual covariance)', fontsize=11)
axes[0].set_xticks(range(M_oj)); axes[0].set_xticklabels(BRANDS_OJ, rotation=45, ha='right', fontsize=8)
axes[0].set_yticks(range(M_oj)); axes[0].set_yticklabels(BRANDS_OJ, fontsize=8)
plt.colorbar(im0, ax=axes[0])
for i in range(M_oj):
    for j in range(M_oj):
        axes[0].text(j, i, f'{Omega_pm[i,j]:.2f}', ha='center', va='center', fontsize=6)

im1 = axes[1].imshow(Corr, vmin=-1, vmax=1, cmap='RdBu_r')
axes[1].set_title('Implied cross-brand error correlation', fontsize=11)
axes[1].set_xticks(range(M_oj)); axes[1].set_xticklabels(BRANDS_OJ, rotation=45, ha='right', fontsize=8)
axes[1].set_yticks(range(M_oj)); axes[1].set_yticklabels(BRANDS_OJ, fontsize=8)
plt.colorbar(im1, ax=axes[1])
for i in range(M_oj):
    for j in range(M_oj):
        axes[1].text(j, i, f'{Corr[i,j]:.2f}', ha='center', va='center', fontsize=6)

plt.suptitle('OJ hSUR — cross-brand residual structure (after controlling own price + deal)', fontsize=11)
plt.tight_layout()
plt.show()

print('Per-brand residual variance (diagonal of Omega):')
print(pd.Series(np.diag(Omega_pm).round(4), index=BRANDS_OJ, name='Var'))

Per-brand residual variance (diagonal of Omega):
Trop64    0.1914
Trop96    0.0751
FloNat    0.2798
Trop64b   0.7113
MM64      0.4651
MM96      0.0816
CitHill   0.3994
TreeFr    0.2525
FloGold   0.7822
Dom64     0.6608
Dom96     0.0936
Name: Var, dtype: float64

# Explanatory: cross-equation residual correlation = the SUR efficiency lever
beta_pm_oj = out_oj['betadraw'].mean(2)                          # (nreg, K)
def _resids(j):
    out = []
    for si, st in enumerate(regdata_oj):
        X = st['X_list'][j]; y = st['y_list'][j]
        bj = beta_pm_oj[si, j * k_oj:(j + 1) * k_oj]
        out.append(y - X @ bj)
    return np.concatenate(out)
C2 = Corr.copy(); np.fill_diagonal(C2, 0.0)
ia, ja = np.unravel_index(np.argmax(C2), C2.shape)              # most positively correlated pair
absC = np.abs(C2) + np.eye(M_oj) * 9                            # mask diagonal
il, jl = np.unravel_index(np.argmin(absC), absC.shape)         # near-independent pair
fig, ax = plt.subplots(1, 2, figsize=(13, 5))
for axi, (p, q, lab) in zip(ax, [(ia, ja, 'most correlated'), (il, jl, 'near-independent')]):
    ra, rb = _resids(p), _resids(q); r = np.corrcoef(ra, rb)[0, 1]
    axi.scatter(ra, rb, s=5, alpha=.25, color='slateblue')
    axi.axhline(0, color='k', lw=.6); axi.axvline(0, color='k', lw=.6)
    axi.set_xlabel(f'{BRANDS_OJ[p]} residual'); axi.set_ylabel(f'{BRANDS_OJ[q]} residual')
    axi.set_title(f'{lab}: corr = {r:.2f}')
plt.suptitle('Why SUR? Cross-equation residual correlation that joint estimation exploits', y=1.0)
plt.tight_layout(); plt.savefig('hsur_sur_resid.png', dpi=120, bbox_inches='tight'); plt.show()

Vbdraw_oj = out_oj['Vbetadraw']   # (R_kept, 33*33)
Vb_pm     = Vbdraw_oj.mean(0).reshape(K_oj, K_oj, order='F')

# Diagonal: variance of each coefficient across stores
vb_diag = np.diag(Vb_pm)
param_labels = [f'{b}.{p}'
                for b in BRANDS_OJ
                for p in ['intercept', 'log_price', 'deal']]

# Focus on own-price variance
price_var = vb_diag[price_idx]
price_sd  = np.sqrt(price_var)

print('Cross-store std of own-price elasticity  (sqrt of Vbeta diagonal):')
print(pd.Series(price_sd.round(4), index=BRANDS_OJ, name='std_price'))

print('\nCross-store std of deal effect:')
deal_var = vb_diag[deal_idx]
print(pd.Series(np.sqrt(deal_var).round(4), index=BRANDS_OJ, name='std_deal'))

Cross-store std of own-price elasticity  (sqrt of Vbeta diagonal):
Trop64    0.7797
Trop96    0.7394
FloNat    0.8763
Trop64b   0.7338
MM64      0.7524
MM96      0.7856
CitHill   0.8840
TreeFr    0.8383
FloGold   0.9829
Dom64     0.7285
Dom96     0.8527
Name: std_price, dtype: float64

Cross-store std of deal effect:
Trop64    0.6537
Trop96    0.6513
FloNat    0.6608
Trop64b   0.6761
MM64      0.6756
MM96      0.6495
CitHill   0.6644
TreeFr    0.6582
FloGold   0.6828
Dom64     0.6737
Dom96     0.6561
Name: std_deal, dtype: float64

# Store-level posterior mean own-price elasticities
bdraw_oj = out_oj['betadraw']              # (nreg, K, R_kept)
beta_pm  = bdraw_oj.mean(axis=2)           # (nreg, K)
price_by_store = beta_pm[:, price_idx]     # (nreg, M_oj)

# Distribution across stores for each brand
print('Store-level own-price elasticity distribution (across 83 stores):')
tbl_store = []
for j, brand in enumerate(BRANDS_OJ):
    e = price_by_store[:, j]
    tbl_store.append({'brand': brand,
                      'pop_mean': df_pop.loc[brand, 'price_mean'],
                      'store_mean': e.mean(),
                      'store_std': e.std(),
                      'store_min': e.min(),
                      'store_max': e.max()})
print(pd.DataFrame(tbl_store).set_index('brand').round(3))

Store-level own-price elasticity distribution (across 83 stores):
         pop_mean  store_mean  store_std  store_min  store_max
brand                                                         
Trop64    -2.7990     -2.8000     0.3970    -3.8280    -1.9940
Trop96    -1.9970     -1.9980     0.3010    -3.0520    -1.3120
FloNat    -3.6730     -3.6730     0.5120    -5.1640    -2.6000
Trop64b   -3.3740     -3.3740     0.2310    -4.0920    -2.8520
MM64      -3.1770     -3.1760     0.3050    -3.9800    -2.5280
MM96      -2.4650     -2.4660     0.4110    -3.2350    -1.3340
CitHill   -4.3580     -4.3590     0.5420    -5.8770    -3.3090
TreeFr    -2.8160     -2.8180     0.4780    -4.3090    -1.7950
FloGold   -4.0570     -4.0570     0.6670    -5.9130    -1.5430
Dom64     -2.9570     -2.9580     0.2490    -3.7320    -2.3470
Dom96     -1.5010     -1.5010     0.5340    -3.2120    -0.5830

fig, ax = plt.subplots(figsize=(12, 5))

ax.boxplot(price_by_store, labels=BRANDS_OJ, vert=True,
           medianprops={'color': 'red', 'linewidth': 1.5},
           whiskerprops={'linewidth': 1},
           flierprops={'markersize': 3, 'alpha': 0.5})
ax.axhline(0, color='black', linewidth=0.8, linestyle='--')

pop_means = [df_pop.loc[b, 'price_mean'] for b in BRANDS_OJ]
ax.scatter(range(1, M_oj + 1), pop_means, color='steelblue', zorder=5,
           s=60, label='Population mean $\\Delta$')

ax.set_ylabel('Own-price elasticity')
ax.set_title('Store-level own-price elasticities across 83 stores\n'
             '(boxes = posterior mean per store; blue dot = population mean $\\Delta$)')
ax.legend()
plt.xticks(rotation=30, ha='right')
plt.tight_layout()
plt.show()

C:\Users\user\AppData\Local\Temp\ipykernel_16616\3816208421.py:3: MatplotlibDeprecationWarning: The 'labels' parameter of boxplot() has been renamed 'tick_labels' since Matplotlib 3.9; support for the old name will be dropped in 3.11.
  ax.boxplot(price_by_store, labels=BRANDS_OJ, vert=True,

# Explanatory: hierarchical shrinkage -- per-store OLS vs hSUR partial pooling
ols_price = np.full((nreg, M_oj), np.nan)
for si, st in enumerate(regdata_oj):
    for j in range(M_oj):
        X = st['X_list'][j]; y = st['y_list'][j]
        ols_price[si, j] = np.linalg.lstsq(X, y, rcond=None)[0][1]
bsel = int(np.argmax(ols_price.std(0)))                         # brand with the noisiest OLS
fig, ax = plt.subplots(1, 2, figsize=(14, 5))
o = ols_price[:, bsel]; h = price_by_store[:, bsel]; pm = df_pop['price_mean'].values[bsel]
for s in range(nreg):
    ax[0].plot([0, 1], [o[s], h[s]], color='gray', alpha=.35, lw=.6)
ax[0].scatter(np.zeros(nreg), o, color='firebrick', s=14, label='per-store OLS (no pooling)')
ax[0].scatter(np.ones(nreg), h, color='steelblue', s=14, label='hSUR (partial pooling)')
ax[0].axhline(pm, color='k', ls='--', lw=1, label=f'population mean $\\Delta$ = {pm:.2f}')
ax[0].set_xticks([0, 1]); ax[0].set_xticklabels(['OLS', 'hSUR']); ax[0].set_xlim(-.3, 1.3)
ax[0].set_ylabel('own-price elasticity'); ax[0].legend(fontsize=8)
ax[0].set_title(f'Shrinkage for {BRANDS_OJ[bsel]} (noisiest OLS): estimates pulled toward $\\Delta$')
x = np.arange(M_oj); w = .4
ax[1].bar(x - w / 2, ols_price.std(0), w, label='OLS (no pooling)', color='firebrick', alpha=.85)
ax[1].bar(x + w / 2, price_by_store.std(0), w, label='hSUR (partial pooling)', color='steelblue', alpha=.85)
ax[1].set_xticks(x); ax[1].set_xticklabels(BRANDS_OJ, rotation=45, ha='right', fontsize=8)
ax[1].set_ylabel('cross-store SD of own-price elasticity'); ax[1].legend(fontsize=8)
ax[1].set_title('Partial pooling shrinks cross-store variance for every brand')
plt.suptitle('Hierarchical shrinkage — the core of the hierarchy', y=1.02)
plt.tight_layout(); plt.savefig('hsur_shrinkage.png', dpi=120, bbox_inches='tight'); plt.show()

DEMO_COLS = ['AGE60', 'EDUC', 'ETHNIC', 'INCOME', 'HHLARGE',
             'WORKWOM', 'HVAL150', 'SSTRDIST', 'SSTRVOL', 'CPDIST5', 'CPWVOL5']

# Align demographics to stores order, standardise, add intercept
demo_vals = demo.set_index('STORE').loc[stores][DEMO_COLS].values.astype(float)
Z_std     = (demo_vals - demo_vals.mean(0)) / demo_vals.std(0)
Z_demo    = np.column_stack([np.ones(nreg), Z_std])   # (83, 12)
nz_demo   = Z_demo.shape[1]

print(f'Z shape: {Z_demo.shape}   (intercept + {nz_demo - 1} standardised demographics)')
print(f'Delta parameters: {nz_demo} × {K_oj} = {nz_demo * K_oj}')

Z shape: (83, 12)   (intercept + 11 standardised demographics)
Delta parameters: 12 × 33 = 396

prior_demo = dict(
    nu      = M_oj + 3,
    V       = (M_oj + 2.0) * np.eye(M_oj),
    nu_beta = K_oj + 3,
    V_beta  = (K_oj + 2.0) * np.eye(K_oj)
)
mcmc_demo = dict(R=5000, burn=1000, keep=1, nprint=1000)

t0       = time.perf_counter()
out_demo = rhsur_gibbs(regdata_oj, Z=Z_demo, Prior=prior_demo, Mcmc=mcmc_demo, seed=42)
print(f'\nDone in {time.perf_counter()-t0:.1f}s  |  kept draws: {out_demo["betadraw"].shape[2]}')

  Iter   1000/5000

  Iter   2000/5000

  Iter   3000/5000

  Iter   4000/5000

  Iter   5000/5000

Done in 170.1s  |  kept draws: 4000

Ddraw_demo = out_demo['Deltadraw']   # (R_kept, nz_demo * K_oj)
# col-major vec(Delta): Delta[i, k] at position k*nz_demo + i
# i=0 is the intercept (pop mean at average demographics); i=1..11 are demographic slopes

COEFF_LABELS = ['intercept', 'price', 'deal']

# ── Significant demographic effects (95% CrI excludes zero) ──────────────────
print('Significant demographic effects on coefficients:')
print(f'  {"Brand":<10} {"Coeff":<10} {"Demo":<12} {"Mean":>8} {"95% CrI":>22}')
print('  ' + '-' * 68)
found = False
for j in range(M_oj):
    for p, clabel in enumerate(COEFF_LABELS):
        k_idx = j * k_oj + p
        for i_d, dname in enumerate(DEMO_COLS):
            d      = Ddraw_demo[:, k_idx * nz_demo + (i_d + 1)]
            lo, hi = np.percentile(d, [2.5, 97.5])
            if lo > 0 or hi < 0:
                print(f'  {BRANDS_OJ[j]:<10} {clabel:<10} {dname:<12} {d.mean():>+8.3f} '
                      f' [{lo:>+7.3f}, {hi:>+7.3f}]')
                found = True
if not found:
    print('  (none at 95% level)')

# ── Population means: z=1 vs z=demo ──────────────────────────────────────────
print()
print('Population-mean price elasticity: z=1 vs z=demo')
print(f'  {"Brand":<10}  {"z=1":>8}  {"z=demo":>8}  {"Δ":>8}')
print('  ' + '-' * 40)
for j in range(M_oj):
    k_idx     = j * k_oj + 1
    demo_mean = Ddraw_demo[:, k_idx * nz_demo + 0].mean()
    z1_mean   = df_pop.loc[BRANDS_OJ[j], 'price_mean']
    print(f'  {BRANDS_OJ[j]:<10}  {z1_mean:>8.4f}  {demo_mean:>8.4f}  {demo_mean - z1_mean:>+8.4f}')

# ── Store-level spread: z=1 vs z=demo ────────────────────────────────────────
print()
print('Store-level price elasticity spread (IQR across 83 stores):')
print(f'  {"Brand":<10}  {"z=1 IQR":>9}  {"z=demo IQR":>11}  {"Δ":>8}')
print('  ' + '-' * 46)
for j in range(M_oj):
    k_idx     = j * k_oj + 1
    beta_z1   = out_oj['betadraw'][:, k_idx, :].mean(axis=1)    # (nreg,)
    beta_demo = out_demo['betadraw'][:, k_idx, :].mean(axis=1)
    iqr_z1   = np.percentile(beta_z1,   75) - np.percentile(beta_z1,   25)
    iqr_demo = np.percentile(beta_demo, 75) - np.percentile(beta_demo, 25)
    print(f'  {BRANDS_OJ[j]:<10}  {iqr_z1:>9.3f}  {iqr_demo:>11.3f}  {iqr_demo - iqr_z1:>+8.3f}')

Significant demographic effects on coefficients:
  Brand      Coeff      Demo             Mean                95% CrI
  --------------------------------------------------------------------
  Trop64     intercept  HVAL150        +1.487  [ +0.671,  +2.278]
  Trop64b    intercept  SSTRDIST       -0.908  [ -1.482,  -0.366]
  Trop64b    price      SSTRDIST       -0.253  [ -0.510,  -0.011]
  TreeFr     intercept  EDUC           -1.596  [ -2.905,  -0.293]
  TreeFr     intercept  HVAL150        +1.477  [ +0.319,  +2.634]

Population-mean price elasticity: z=1 vs z=demo
  Brand            z=1    z=demo         Δ
  ----------------------------------------
  Trop64       -2.7992   -2.8102   -0.0110
  Trop96       -1.9969   -1.9960   +0.0009
  FloNat       -3.6726   -3.6696   +0.0030
  Trop64b      -3.3736   -3.3514   +0.0222
  MM64         -3.1769   -3.1806   -0.0037
  MM96         -2.4653   -2.4650   +0.0003
  CitHill      -4.3581   -4.3565   +0.0017
  TreeFr       -2.8157   -2.8266   -0.0109
  FloGold      -4.0570   -4.0636   -0.0066
  Dom64        -2.9569   -2.9465   +0.0104
  Dom96        -1.5009   -1.5046   -0.0037

Store-level price elasticity spread (IQR across 83 stores):
  Brand         z=1 IQR   z=demo IQR         Δ
  ----------------------------------------------
  Trop64          0.563        0.622    +0.059
  Trop96          0.339        0.342    +0.002
  FloNat          0.561        0.586    +0.025
  Trop64b         0.309        0.482    +0.172
  MM64            0.397        0.496    +0.099
  MM96            0.508        0.507    -0.000
  CitHill         0.609        0.805    +0.196
  TreeFr          0.610        0.695    +0.085
  FloGold         0.792        0.783    -0.009
  Dom64           0.355        0.390    +0.035
  Dom96           0.832        0.803    -0.029

k_cross  = 13          # intercept + 11 log-prices + own deal
K_cross  = M_oj * k_cross   # 143

# For brand j (0-indexed): own-price column within X = j+1 (log_p_{j+1})
# flat index in K vector: j*k_cross + (j+1)
own_price_idx_cross  = [j * k_cross + (j + 1) for j in range(M_oj)]
deal_idx_cross       = [j * k_cross + 12       for j in range(M_oj)]

regdata_cross = []
for store_id in stores:
    store_df = yx[yx['store'] == store_id]
    y_list = []; X_list = []
    for j in range(M_oj):
        brand_id = j + 1
        bdf      = store_df[store_df['brand'] == brand_id].sort_values('week')
        y_j      = bdf['logmove'].values.astype(float)
        deal_j   = bdf['deal'].values.astype(float)
        log_p    = np.column_stack([
            np.log(bdf[f'price{q+1}'].values.astype(float))
            for q in range(M_oj)
        ])                                          # (T_i, 11) — all log-prices
        X_j = np.column_stack([np.ones(len(y_j)), log_p, deal_j])   # (T_i, 13)
        y_list.append(y_j); X_list.append(X_j)
    regdata_cross.append({'y_list': y_list, 'X_list': X_list})

print(f'regdata_cross: {len(regdata_cross)} stores')
print(f'k={k_cross}  K={K_cross}  (own-price indices: {own_price_idx_cross[:3]} ...)')
print(f'X shape check: {regdata_cross[0]["X_list"][0].shape}')

regdata_cross: 83 stores
k=13  K=143  (own-price indices: [1, 15, 29] ...)
X shape check: (110, 13)

prior_cross = dict(
    nu      = M_oj + 3,
    V       = (M_oj + 2.0) * np.eye(M_oj),
    nu_beta = K_cross + 3,
    V_beta  = (K_cross + 2.0) * np.eye(K_cross)
)
mcmc_cross = dict(R=5000, burn=1000, keep=1, nprint=1000)

t0        = time.perf_counter()
out_cross = rhsur_gibbs(regdata_cross, Z=np.ones((nreg, 1)),
                        Prior=prior_cross, Mcmc=mcmc_cross, seed=42)
elapsed   = time.perf_counter() - t0
print(f'\nDone in {elapsed/60:.1f} min  |  kept draws: {out_cross["betadraw"].shape[2]}')

  Iter   1000/5000

  Iter   2000/5000

  Iter   3000/5000

  Iter   4000/5000

  Iter   5000/5000

Done in 7.2 min  |  kept draws: 4000

Ddraw_cross = out_cross['Deltadraw']   # (R_kept, 143)  — nz=1 so vec(Delta)=Delta[0,:]

# ── Own-price comparison: restricted vs unrestricted ─────────────────────────
tbl_own = []
for j, brand in enumerate(BRANDS_OJ):
    d_restr = Ddraw_oj[:,   price_idx[j]]
    d_cross = Ddraw_cross[:, own_price_idx_cross[j]]
    lo, hi  = np.percentile(d_cross, [2.5, 97.5])
    tbl_own.append({
        'brand'     : brand,
        'restricted': d_restr.mean(),
        'unrestr'   : d_cross.mean(),
        'shift'     : d_cross.mean() - d_restr.mean(),
        'lo'        : lo,
        'hi'        : hi,
    })
df_own = pd.DataFrame(tbl_own).set_index('brand')
print('Own-price elasticities — restricted vs unrestricted (population means):')
print(df_own.round(4))

# ── Cross-price elasticity matrix (population means) ─────────────────────────
# cprice[j, q] = effect of log(price_q) on log(sales_j)
cprice = np.zeros((M_oj, M_oj))
for j in range(M_oj):
    for q in range(M_oj):
        idx = j * k_cross + (q + 1)
        cprice[j, q] = Ddraw_cross[:, idx].mean()

print('\nCross-price elasticity matrix (rows = equation brand, cols = price brand):')
df_cp = pd.DataFrame(cprice, index=BRANDS_OJ, columns=BRANDS_OJ)
print(df_cp.round(3))

# ── Cross-price heatmap ───────────────────────────────────────────────────────
fig, ax = plt.subplots(figsize=(11, 9))
vmax = max(abs(cprice[cprice != np.diag(cprice)].max()),
           abs(cprice[cprice != np.diag(cprice)].min()))
im = ax.imshow(cprice, cmap='RdBu_r', vmin=-vmax*3, vmax=vmax*3)
ax.set_xticks(range(M_oj)); ax.set_xticklabels(BRANDS_OJ, rotation=45, ha='right', fontsize=9)
ax.set_yticks(range(M_oj)); ax.set_yticklabels(BRANDS_OJ, fontsize=9)
ax.set_xlabel('Price of brand (column)'); ax.set_ylabel('Sales of brand (row)')
ax.set_title('Population-mean cross-price elasticities\n'
             '(diagonal = own-price; off-diagonal > 0 ≈ substitutes)')
for i in range(M_oj):
    for j in range(M_oj):
        ax.text(j, i, f'{cprice[i,j]:.2f}', ha='center', va='center', fontsize=7,
                color='white' if abs(cprice[i,j]) > vmax else 'black')
plt.colorbar(im, ax=ax, shrink=0.8)
plt.tight_layout()
plt.show()

Own-price elasticities — restricted vs unrestricted (population means):
         restricted  unrestr   shift      lo      hi
brand                                               
Trop64      -2.7992  -3.1489 -0.3497 -3.4658 -2.8384
Trop96      -1.9969  -1.9269  0.0700 -2.2518 -1.6135
FloNat      -3.6726  -3.9161 -0.2435 -4.2423 -3.5910
Trop64b     -3.3736  -4.4150 -1.0414 -4.7296 -4.0910
MM64        -3.1769  -3.8829 -0.7060 -4.2096 -3.5644
MM96        -2.4653  -2.6540 -0.1887 -2.9806 -2.3269
CitHill     -4.3581  -4.7483 -0.3901 -5.0821 -4.4247
TreeFr      -2.8157  -2.8912 -0.0754 -3.2190 -2.5538
FloGold     -4.0570  -4.2307 -0.1737 -4.5803 -3.8677
Dom64       -2.9569  -2.9684 -0.0115 -3.2858 -2.6613
Dom96       -1.5009  -1.4912  0.0097 -1.8235 -1.1579

Cross-price elasticity matrix (rows = equation brand, cols = price brand):
         Trop64  Trop96  FloNat  Trop64b    MM64    MM96  CitHill  TreeFr  \
Trop64  -3.1490 -0.0230  0.2800   0.1900  0.1360 -0.0440   0.3030  0.0810   
Trop96   0.4080 -1.9270  0.0770  -0.0040  0.0470  0.1690   0.1320  0.1850   
FloNat   0.8260  0.4370 -3.9160   0.3970  0.1550  0.0190  -0.0040 -0.0790   
Trop64b  1.0140  0.5330 -0.0790  -4.4150  1.2650  2.2800   0.4020  0.3240   
MM64     0.8750  0.2030  0.9630   0.9470 -3.8830 -0.1450   0.5950  0.1410   
MM96     0.0650  0.5780  0.1420   0.0890  0.2270 -2.6540   0.1370  0.2180   
CitHill  0.5060  0.5460  0.1140   0.0730  0.8340  0.6270  -4.7480  0.4010   
TreeFr   0.0930  0.3240  0.1510   0.1430  0.3610  0.0120   0.3820 -2.8910   
FloGold  0.0760 -0.0590  0.0000   0.7820  1.0210  0.1990   0.6790 -0.2200   
Dom64    0.0680 -1.5170  0.4410   1.1620  1.2500  0.2550  -0.2550 -0.0810   
Dom96   -0.1850 -0.0270  0.0420   0.1000 -0.0280  0.3530   0.1330 -0.0920   

         FloGold   Dom64   Dom96  
Trop64   -0.0240  0.1450  0.4510  
Trop96    0.0080  0.0570 -0.1320  
FloNat   -0.2910  0.0920  0.6520  
Trop64b   0.2110  0.0200 -1.2630  
MM64      0.1820  0.6760  0.0630  
MM96      0.0420  0.0290  0.0990  
CitHill   0.6430  0.0890  0.3140  
TreeFr    0.1670  0.4010  0.9040  
FloGold  -4.2310  0.0740  0.4910  
Dom64     0.2670 -2.9680 -0.4780  
Dom96    -0.0290  0.1360 -1.4910

Feature	SUR	HLM	hSUR
Units	1 market, M equations	N stores, 1 equation	N stores, M equations
Cross-unit link	—	Prior $\beta_i \sim N(\Delta' z_i, V_\beta)$	Both
Cross-equation link	$\boldsymbol{\Omega}$	—	$\boldsymbol{\Omega}$
Regressors	Different per equation	Same per unit	Different per equation & unit
Gibbs blocks	2	3–4	4

File	Content
`sur_bayesm.ipynb`	R `rsurGibbs` — basic SUR reference
`sur_python.ipynb`	Python SUR + tuna application
`hsur_gibbs.py`	hSUR Gibbs sampler module
`hsur_python.ipynb` (this)	hSUR validation + OJ application

#	Brand
1	Tropicana Premium 64 oz
2	Tropicana Premium 96 oz
3	Florida's Natural
4	Tropicana 64 oz
5	Minute Maid 64 oz
6	Minute Maid 96 oz
7	Citrus Hill 64 oz
8	Tree Fresh 64 oz
9	Florida Gold 64 oz
10	Dominick's 64 oz
11	Dominick's 96 oz

Tier	Brands	Elasticity	Interpretation
Premium national	Trop Premium 64/96	−2.80, −2.00	Strong loyalty; low price sensitivity
Standard national	MM64, MM96, Trop64b	−3.18, −2.47, −3.37	Moderate loyalty
Regional/value	CitHill, FloGold, FloNat, TreeFr	−4.36, −4.06, −3.67, −2.82	Less loyal; readily substitute
Private label	Dom64, Dom96	−2.96, −1.50	Mixed — see below

Brand	Deal	Brand	Deal
Trop64b	+0.655	Trop Premium 64	+0.066
Dom64	+0.520	TreeFr	+0.067
MM64	+0.420	Trop96	+0.189
FloGold	+0.403	MM96	+0.116
Dom96	+0.365	CitHill	+0.194

Format	Variance	Ratio to 64 oz sibling
Trop Premium 64	0.191	—
Trop Premium 96	0.075	2.5× smaller
Trop 64	0.711	—
MM 64	0.465	—
MM 96	0.082	5.7× smaller
Dom 64	0.661	—
Dom 96	0.094	7.0× smaller

Brand	pop_mean	store_std	Combination
Trop64b	−3.37	0.231	Consistently moderately elastic
Dom64	−2.96	0.249	Consistently moderately elastic
MM64	−3.18	0.305	Consistent moderate
Trop96	−2.00	0.301	Low elastic, moderate spread
Dom96	−1.50	0.534	Least elastic on average, wide spread
FloGold	−4.06	0.667	Highly elastic on average, widest spread

Store type	FloGold elasticity	Optimal price response
"Premium regional" store	≈ −1.54	Has pricing power — modest price increase loses few units
"Discount commodity" store	≈ −5.91	Highly elastic — small price increase collapses volume

Approach	What it misses
83 × 11 separate OLS	No shrinkage; no cross-brand $\Omega$
Pooled SUR (1 β per brand)	No store heterogeneity
83 × separate SUR	No shrinkage; 83 imprecise $\Omega$ matrices
hSUR (this model)	Nothing from the above list

Variable	Description
`AGE60`	Fraction of households whose head is aged 60 or over
`EDUC`	Fraction of households with at least a college education
`ETHNIC`	Fraction of non-white (minority) households
`INCOME`	Median household income in the store's trade area
`HHLARGE`	Fraction of large households (5 or more members)
`WORKWOM`	Fraction of households with a working woman
`HVAL150`	Fraction of households with home value ≥ $150,000
`SSTRDIST`	Distance (miles) to the nearest warehouse/club store
`SSTRVOL`	Volume index of the nearest warehouse/club store
`CPDIST5`	Number of competing grocery stores within 5 miles
`CPWVOL5`	Total volume index of competing grocery stores within 5 miles

Brand	$\Delta$ IQR (z=1 → z=demo)
FloGold	0.898 → 0.894 (−0.5%)
Dom96	0.724 → 0.731 (+1.0%)
CitHill	0.729 → 0.722 (−1.0%)
All others	< 1% change

Brand	Own-price (restricted, pooled OLS)	Own-price (unrestricted, pooled OLS)	Shift	Comment
Dom96	−1.01	−0.26	+0.75	"Least elastic" may be confounding
CitHill	−2.78	−4.66	−1.88	May be more elastic than it appears
Dom64	−3.22	−3.21	−0.01	Clean — private label 64 oz uncorrelated

Rank	Brand	Elasticity	Restricted rank
1	CitHill	−4.748	1
2	Trop64b	−4.415	4 ↑
3	FloGold	−4.231	2
4	FloNat	−3.916	3
5	MM64	−3.883	5
6	Trop64	−3.149	8 ↑
7	Dom64	−2.968	6
8	TreeFr	−2.891	7
9	MM96	−2.654	9
10	Trop96	−1.927	10
11	Dom96	−1.491	11

Hierarchical SUR — Python (hsur_gibbs.py)¶

Model hierarchy¶

Notebooks¶

1. Imports¶

2. Synthetic Data Validation¶

3. rhsur_gibbs — Synthetic Data¶

4. Orange Juice Data¶

5. rhsur_gibbs — Orange Juice¶

6. Population-Level Results¶

Parameter layout — why the own-price elasticity sits at $3j+1$¶

Reading $\Delta$¶

OJ results (validated)¶

Analysis — population-level price elasticities and deal effects¶

Omega diagonal — residual variances (validated)¶

Analysis — residual structure ($\Omega$)¶

96 oz formats are systematically more predictable¶

The high-variance group: Trop64b, FloGold, Dom64¶

The low-variance group: Trop96, MM96, Dom96¶

Link to Omega off-diagonals¶

7. Cross-Store Heterogeneity ($V_\beta$)¶

Store-level heterogeneity — validated¶

Analysis — cross-store heterogeneity in price elasticity¶

The shrinkage check¶

Heterogeneity ranking inverts the loyalty ranking¶

Dom96: least elastic but widest-ranging¶

FloGold: the positioning-ambiguous brand¶

Practical implication: store-segment pricing¶

8. Synthesis¶

Three layers of insight from hSUR¶

What a pooled SUR would miss¶

Natural extension: store demographics in $Z$¶

9. Extension: Store Demographics in $Z$¶

Store demographic variables¶

9. Demographics results (validated)¶

Analysis — what demographics explain (and don't)¶

Five significant effects — all on sales levels, not price elasticities¶

Demographics do NOT explain cross-store elasticity heterogeneity¶

Vbeta inflation: a note on identification¶

10. Unrestricted Model: Cross-Price Elasticities¶

Motivation¶

Unrestricted model specification¶

10. Results (validated)¶

Analysis — cross-price results¶

The OLS paradox for Dom96 is resolved by the hierarchy¶

Trop64b is the most confounded brand (−1.04 shift)¶

MM64 also materially biased (−0.71 shift)¶

The stable brands confirm the hierarchy is doing the right thing¶

Cross-price patterns: substitution follows format and tier¶

Revised brand elasticity story¶

Hierarchical SUR — Python (`hsur_gibbs.py`)¶

3. `rhsur_gibbs` — Synthetic Data¶

5. `rhsur_gibbs` — Orange Juice¶