import os, sys, time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# Force reload: drop cached module first so Python reads the file on disk
if 'mnprobit_gibbs' in sys.modules:
    del sys.modules['mnprobit_gibbs']

import mnprobit_gibbs
from mnprobit_gibbs import mnprobit_gibbs as gibbs_sampler, make_prior, make_init

# Verify n_sweeps is in the module source (it's a Mcmc dict key, not a function arg)
import inspect
src = inspect.getsource(mnprobit_gibbs)
assert "n_sweeps" in src, "Module not updated — check mnprobit_gibbs.py was saved"
print("Module loaded OK — n_sweeps confirmed in source")

np.random.seed(2025)
print("Imports OK")

Module loaded OK — n_sweeps confirmed in source
Imports OK

y_synth = pd.read_csv("mnp_synth_y.csv")["y"].values.astype(int)
X_synth = pd.read_csv("mnp_synth_X.csv").values.astype(float)

N   = len(y_synth)
p   = int(y_synth.max())
pm1 = p - 1
k   = X_synth.shape[1]

# True parameters (for recovery check)
BETA_TRUE  = -2.0
SIGMA_TRUE = np.array([[1.00, 0.50, 0.30],
                        [0.50, 1.00, 0.40],
                        [0.30, 0.40, 1.00]])

print(f"Synthetic data: N={N}, p={p}, k={k}")
print(f"Choice frequencies: {dict(zip(*np.unique(y_synth, return_counts=True)))}")

Synthetic data: N=2000, p=4, k=1
Choice frequencies: {np.int64(1): np.int64(493), np.int64(2): np.int64(517), np.int64(3): np.int64(471), np.int64(4): np.int64(519)}

NDRAWS_S = 20000
BURN_S   =  5000

prior_s = make_prior(pm1, k)
init_s  = make_init(p, k)

Data_s  = {'p': p, 'y': y_synth, 'X': X_synth}
Mcmc_s  = {'R': NDRAWS_S, 'keep': 1, 'nprint': 3000, **init_s}

np.random.seed(2024)   # fix seed for reproducibility
print(f"Running synthetic Gibbs (N={N}, {NDRAWS_S} draws) ...")
t0 = time.time()
out_s = gibbs_sampler(Data_s, prior_s, Mcmc_s)
print(f"Done in {time.time()-t0:.1f} s")

Running synthetic Gibbs (N=2000, 20000 draws) ...
MNP Gibbs: n=2000, p=4, k=1, R=20000, keep=1, n_sweeps=1, brand_price=False
  Iteration   3000 / 20000
  Iteration   6000 / 20000
  Iteration   9000 / 20000
  Iteration  12000 / 20000
  Iteration  15000 / 20000
  Iteration  18000 / 20000
Done in 53.2 s

beta_draws_s  = out_s['betadraw'][BURN_S:]
sigma_draws_s = out_s['sigmadraw'][BURN_S:]

beta_mean_s  = float(beta_draws_s.mean())
sigma_mat_s  = sigma_draws_s.mean(axis=0).reshape(pm1, pm1)

print(f"sigma_11 (should = 1): {sigma_mat_s[0, 0]:.6f}")
print(f"\nbeta* — true: {BETA_TRUE:.3f}  |  posterior mean: {beta_mean_s:.3f}  |  diff: {beta_mean_s - BETA_TRUE:.4f}")

print("\nSigma* — posterior mean vs true:")
df_sig = pd.DataFrame({
    'element'  : [f'Sigma[{r+1},{c+1}]' for r in range(pm1) for c in range(pm1)],
    'true_val' : SIGMA_TRUE.ravel().round(4),
    'post_mean': sigma_mat_s.ravel().round(4),
    'diff'     : (sigma_mat_s - SIGMA_TRUE).ravel().round(4),
})
print(df_sig.to_string(index=False))
print(f"Max |diff| Sigma: {np.abs(sigma_mat_s - SIGMA_TRUE).max():.4f}")

sigma_11 (should = 1): 1.000000

beta* — true: -2.000  |  posterior mean: -1.851  |  diff: 0.1485

Sigma* — posterior mean vs true:
   element  true_val  post_mean    diff
Sigma[1,1]       1.0     1.0000  0.0000
Sigma[1,2]       0.5     0.4454 -0.0546
Sigma[1,3]       0.3     0.3378  0.0378
Sigma[2,1]       0.5     0.4454 -0.0546
Sigma[2,2]       1.0     0.7416 -0.2584
Sigma[2,3]       0.4     0.2298 -0.1702
Sigma[3,1]       0.3     0.3378  0.0378
Sigma[3,2]       0.4     0.2298 -0.1702
Sigma[3,3]       1.0     0.8190 -0.1810
Max |diff| Sigma: 0.2584

np.savez("synth_gibbs.npz",
         beta_mean  = beta_mean_s,
         beta_std   = float(beta_draws_s.std()),
         beta_ci    = np.percentile(beta_draws_s, [2.5, 97.5]),
         sigma_mean = sigma_mat_s,
         n_draws    = len(beta_draws_s))
print("Saved synth_gibbs.npz")

# CSV versions for R (MNP_bayesm.ipynb §7 load_stats)
_ci_s = np.percentile(beta_draws_s, [2.5, 97.5])
pd.DataFrame({
    "beta_mean": [beta_mean_s],
    "beta_std":  [float(beta_draws_s.std())],
    "ci_lo":     [_ci_s[0]],
    "ci_hi":     [_ci_s[1]],
    "n_draws":   [len(beta_draws_s)]
}).to_csv("synth_gibbs_stats.csv", index=False)
pd.DataFrame(sigma_mat_s).to_csv("synth_gibbs_sigma.csv")
print("Saved synth_gibbs_stats.csv  and  synth_gibbs_sigma.csv")

Saved synth_gibbs.npz
Saved synth_gibbs_stats.csv  and  synth_gibbs_sigma.csv

# Figure — identification: choice data pin down beta but not Sigma
sig3d_s = sigma_draws_s.reshape(len(sigma_draws_s), pm1, pm1)
fig, ax = plt.subplots(1, 3, figsize=(12, 3.6))
panels = [(beta_draws_s.ravel(), -2.0, r'price coefficient  $\beta^*$',      'true $-2.0$'),
          (sig3d_s[:, 1, 1],      1.0,  r'$\Sigma^*_{22}$  (a free diagonal)', 'true $1.0$'),
          (sig3d_s[:, 1, 2],      0.4,  r'$\Sigma^*_{23}$  (a correlation)',   'true $0.4$')]
for a, (draws, truth, title, tl) in zip(ax, panels):
    a.hist(draws, bins=45, density=True, color='steelblue', edgecolor='white', alpha=0.85)
    a.axvline(truth, color='firebrick', lw=2, ls='--', label=tl)
    a.axvline(draws.mean(), color='black', lw=1.2, label=f'post mean {draws.mean():.2f}')
    a.set_title(title, fontsize=10); a.set_yticks([]); a.legend(fontsize=8)
fig.suptitle('Choice data identify $\\beta$ sharply, $\\Sigma$ weakly\n'
             'the $\\beta$ posterior brackets the truth; the $\\Sigma$ posteriors are wide and biased, their true values out in the tails',
             fontsize=11)
plt.tight_layout(); plt.show()

y_marg = pd.read_csv("mnp_marg4_y.csv")["y"].values.astype(int)
X_marg = pd.read_csv("mnp_marg4_X.csv").values.astype(float)

BRANDS4 = ["Parkay Stk", "Blue Bonnet Stk", "House Stk", "Shedd Tub"]
N4  = len(y_marg)
p4  = 4
pm4 = p4 - 1

print(f"Margarine (top 4): N={N4}, p={p4}")
shares = {b: (y_marg == i+1).mean() for i, b in enumerate(BRANDS4)}
print("Observed shares:", {k: round(v, 3) for k, v in shares.items()})

Margarine (top 4): N=3377, p=4
Observed shares: {'Parkay Stk': np.float64(0.523), 'Blue Bonnet Stk': np.float64(0.207), 'House Stk': np.float64(0.176), 'Shedd Tub': np.float64(0.094)}

NDRAWS_M = 20000
BURN_M   =  5000

prior_m = make_prior(pm4, 1)
init_m  = make_init(p4, 1)

Data_m  = {'p': p4, 'y': y_marg, 'X': X_marg}
Mcmc_m  = {'R': NDRAWS_M, 'keep': 1, 'nprint': 5000, 'n_sweeps': 5, **init_m}

np.random.seed(2025)   # fix seed for reproducibility
print(f"Running M1 ({NDRAWS_M} draws, n_sweeps=5) ...")
t0 = time.time()
out_m = gibbs_sampler(Data_m, prior_m, Mcmc_m)
print(f"Done in {time.time()-t0:.1f} s")

Running M1 (20000 draws, n_sweeps=5) ...
MNP Gibbs: n=3377, p=4, k=1, R=20000, keep=1, n_sweeps=5, brand_price=False
  Iteration   5000 / 20000
  Iteration  10000 / 20000
  Iteration  15000 / 20000
  Iteration  20000 / 20000
Done in 260.7 s

# Extract posterior draws (drop burn-in)
beta_draws_m  = out_m['betadraw'][BURN_M:]     # (n_post, 1)
sigma_draws_m = out_m['sigmadraw'][BURN_M:]    # (n_post, 9)

n_post       = len(beta_draws_m)
beta_chain   = beta_draws_m.ravel()            # (n_post,)
S3d          = sigma_draws_m.reshape(n_post, pm4, pm4)
sigma_mat_m  = S3d.mean(axis=0)

ALT_LABELS = ['BB Stk', 'House Stk', 'Shedd Tub']

print(f"Post-burn-in draws : {n_post:,}")
print(f"sigma_11 (should=1): {sigma_mat_m[0,0]:.6f}")
print(f"\nPrice coefficient beta*  mean = {beta_chain.mean():.4f}   SD = {beta_chain.std():.4f}")
print(f"  95% CI: [{np.percentile(beta_chain,2.5):.4f}, {np.percentile(beta_chain,97.5):.4f}]")
print(f"\nPosterior mean Sigma*:")
print(pd.DataFrame(sigma_mat_m, index=ALT_LABELS, columns=ALT_LABELS).round(4))

Post-burn-in draws : 15,000
sigma_11 (should=1): 1.000000

Price coefficient beta*  mean = -1.5179   SD = 0.0956
  95% CI: [-1.7014, -1.3291]

Posterior mean Sigma*:
           BB Stk  House Stk  Shedd Tub
BB Stk     1.0000     0.6608     1.2198
House Stk  0.6608     0.5045     0.8557
Shedd Tub  1.2198     0.8557     1.7033

from statsmodels.tsa.stattools import acf as compute_acf

def ess(chain, max_lag=500):
    n = len(chain)
    acf_v = compute_acf(chain, nlags=min(max_lag, n // 2), fft=True)
    ci = 1.96 / np.sqrt(n)
    rho = 0.0
    for k in range(1, len(acf_v)):
        if abs(acf_v[k]) < ci:
            break
        rho += acf_v[k]
    return n / max(1.0, 1 + 2 * rho)

# All monitored parameters
monitor = {
    'beta* (price)'    : beta_chain,
    'Sigma[BB,House]'  : S3d[:, 0, 1],
    'Sigma[BB,Shedd]'  : S3d[:, 0, 2],
    'Sigma[House,House]': S3d[:, 1, 1],
    'Sigma[House,Shedd]': S3d[:, 1, 2],
    'Sigma[Shedd,Shedd]': S3d[:, 2, 2],
}

print(f"{'Parameter':<22}  {'Mean':>8}  {'SD':>8}  {'ESS':>7}  {'ESS%':>7}")
print("-" * 60)
ess_vals = {}
for name, chain in monitor.items():
    e = ess(chain)
    ess_vals[name] = e
    print(f"{name:<22}  {chain.mean():>8.4f}  {chain.std():>8.4f}  {e:>7.0f}  {e/n_post:>7.2%}")

Parameter                   Mean        SD      ESS     ESS%
------------------------------------------------------------
beta* (price)            -1.5179    0.0956      522    3.48%
Sigma[BB,House]           0.6608    0.0247      189    1.26%
Sigma[BB,Shedd]           1.2198    0.0307      279    1.86%
Sigma[House,House]        0.5045    0.0349      190    1.27%
Sigma[House,Shedd]        0.8557    0.0412      226    1.51%
Sigma[Shedd,Shedd]        1.7033    0.0851      384    2.56%

iters  = np.arange(1, n_post + 1)
names  = list(monitor.keys())
chains = list(monitor.values())
n_par  = len(names)
NLAGS  = 60
ci_acf = 1.96 / np.sqrt(n_post)

fig, axes = plt.subplots(n_par, 3, figsize=(15, 2.1 * n_par))

for i, (name, chain) in enumerate(zip(names, chains)):
    # Trace
    axes[i, 0].plot(iters, chain, lw=0.3, color='steelblue', alpha=0.8)
    axes[i, 0].axhline(chain.mean(), color='red', lw=1, ls='--')
    axes[i, 0].set_ylabel(name, fontsize=8, rotation=0, ha='right', labelpad=80)
    if i == 0:
        axes[i, 0].set_title('Trace', fontsize=10)

    # ACF
    acf_v = compute_acf(chain, nlags=NLAGS, fft=True)[1:]
    axes[i, 1].bar(np.arange(1, NLAGS + 1), acf_v, color='steelblue', width=0.8, alpha=0.7)
    axes[i, 1].axhline( ci_acf, color='red', lw=1, ls='--')
    axes[i, 1].axhline(-ci_acf, color='red', lw=1, ls='--')
    axes[i, 1].axhline(0, color='k', lw=0.5)
    axes[i, 1].set_ylim(-0.3, 1.0)
    if i == 0:
        axes[i, 1].set_title('ACF', fontsize=10)

    # Running mean
    rm = np.cumsum(chain) / iters
    axes[i, 2].plot(iters, rm, lw=1.2, color='steelblue')
    axes[i, 2].axhline(chain.mean(), color='red', lw=1, ls='--')
    if i == 0:
        axes[i, 2].set_title('Running mean', fontsize=10)

for ax in axes[-1]:
    ax.set_xlabel('Post-burn-in draw' if ax != axes[-1, 1] else 'Lag')

plt.suptitle('MCMC diagnostics — margarine MNP', fontsize=12, y=1.01)
plt.tight_layout()
plt.show()

ci95 = np.percentile(beta_chain, [2.5, 97.5])

fig, ax = plt.subplots(figsize=(7, 4))
ax.hist(beta_chain, bins=60, density=True, color='steelblue', edgecolor='white')
ax.axvline(beta_chain.mean(), color='red',  ls='--', lw=1.5,
           label=f'Mean = {beta_chain.mean():.3f}')
ax.axvline(ci95[0], color='orange', ls=':', lw=1.5)
ax.axvline(ci95[1], color='orange', ls=':', lw=1.5,
           label=f'95% CI [{ci95[0]:.3f}, {ci95[1]:.3f}]')
ax.axvline(0, color='black', lw=1, label='zero')
ax.set_xlabel('Price coefficient β*')
ax.set_title('M1 posterior: price sensitivity')
ax.legend()
plt.tight_layout(); plt.show()

print(f"β* posterior mean : {beta_chain.mean():.4f}")
print(f"β* posterior SD   : {beta_chain.std():.4f}")
print(f"95% CI            : [{ci95[0]:.4f}, {ci95[1]:.4f}]")
print(f"P(β* < 0)         : {(beta_chain < 0).mean():.4f}")
print()
print("Price-change implications (Δprice vs Parkay):")
for dp in [0.05, 0.10, 0.20]:
    dmu = beta_chain.mean() * dp
    print(f"  Δprice = +${dp:.2f}  →  Δutility diff = {dmu:+.4f}  "
          f"(≈ {abs(dmu)/beta_chain.std():.1f} posterior SDs)")

β* posterior mean : -1.5179
β* posterior SD   : 0.0956
95% CI            : [-1.7014, -1.3291]
P(β* < 0)         : 1.0000

Price-change implications (Δprice vs Parkay):
  Δprice = +$0.05  →  Δutility diff = -0.0759  (≈ 0.8 posterior SDs)
  Δprice = +$0.10  →  Δutility diff = -0.1518  (≈ 1.6 posterior SDs)
  Δprice = +$0.20  →  Δutility diff = -0.3036  (≈ 3.2 posterior SDs)

### 7b. Covariance / correlation structure (Σ*)
import seaborn as sns

D_inv  = np.diag(1.0 / np.sqrt(np.diag(sigma_mat_m)))
corr_m = D_inv @ sigma_mat_m @ D_inv

print("Posterior mean Σ*  (covariance of utility diffs vs Parkay):")
print(pd.DataFrame(sigma_mat_m, index=ALT_LABELS, columns=ALT_LABELS).round(4))
print("\nImplied correlation matrix  [ρ_jk = Σ_jk / √(Σ_jj Σ_kk)]:")
print(pd.DataFrame(corr_m, index=ALT_LABELS, columns=ALT_LABELS).round(4))

rho12_draws = S3d[:,0,1] / np.sqrt(S3d[:,1,1])
rho13_draws = S3d[:,0,2] / np.sqrt(S3d[:,2,2])
rho23_draws = S3d[:,1,2] / np.sqrt(S3d[:,1,1]*S3d[:,2,2])

fig, axes = plt.subplots(1, 3, figsize=(13, 4))
for ax, rho, lbl in zip(axes,
                         [rho12_draws, rho13_draws, rho23_draws],
                         ['ρ(BB, House)', 'ρ(BB, Shedd)', 'ρ(House, Shedd)']):
    ax.hist(rho, bins=50, color='steelblue', edgecolor='white', density=True)
    ax.axvline(rho.mean(), color='red', ls='--', lw=1.5, label=f'Mean={rho.mean():.3f}')
    ci = np.percentile(rho, [2.5, 97.5])
    ax.axvline(ci[0], color='orange', ls=':', lw=1)
    ax.axvline(ci[1], color='orange', ls=':', lw=1,
               label=f'95%CI [{ci[0]:.2f},{ci[1]:.2f}]')
    ax.axvline(0, color='k', lw=0.8)
    ax.set_title(lbl); ax.legend(fontsize=8)
plt.suptitle('M1 — posterior distributions of implied correlations', fontsize=11)
plt.tight_layout(); plt.show()

fig, axes = plt.subplots(1, 2, figsize=(11, 4))
for ax, mat, title in zip(axes, [sigma_mat_m, corr_m],
                           ['Σ* (covariance)', 'Implied correlations ρ']):
    sns.heatmap(mat, annot=True, fmt='.3f', center=0, cmap='RdBu_r',
                xticklabels=ALT_LABELS, yticklabels=ALT_LABELS, ax=ax,
                annot_kws={'size': 13}, linewidths=0.5)
    ax.set_title(title, fontsize=11)
plt.suptitle('Margarine M1 — posterior mean Σ*', y=1.02)
plt.tight_layout(); plt.show()

Posterior mean Σ*  (covariance of utility diffs vs Parkay):
           BB Stk  House Stk  Shedd Tub
BB Stk     1.0000     0.6608     1.2198
House Stk  0.6608     0.5045     0.8557
Shedd Tub  1.2198     0.8557     1.7033

Implied correlation matrix  [ρ_jk = Σ_jk / √(Σ_jj Σ_kk)]:
           BB Stk  House Stk  Shedd Tub
BB Stk     1.0000     0.9304     0.9346
House Stk  0.9304     1.0000     0.9231
Shedd Tub  0.9346     0.9231     1.0000

### 7c. Choice probabilities — posterior predictive

obs_shares = np.array([(y_marg == j + 1).mean() for j in range(4)])
X_cube_m   = X_marg.reshape(N4, pm4, 1)
mean_diffs  = X_cube_m[:, :, 0].mean(axis=0)

print(f"Mean price diffs vs Parkay:  BB={mean_diffs[0]:+.3f}  "
      f"House={mean_diffs[1]:+.3f}  Shedd={mean_diffs[2]:+.3f}")

def choice_probs_mc(price_diffs, beta_sample, sigma_sample, pm1, n_mc=1500, seed=42):
    rng = np.random.default_rng(seed)
    p = pm1 + 1;  counts = np.zeros(p)
    for beta_row, sig_flat in zip(beta_sample, sigma_sample):
        S  = sig_flat.reshape(pm1, pm1)
        mu = np.asarray(price_diffs) * float(beta_row.ravel()[0])
        w  = rng.multivariate_normal(mu, S, size=n_mc)
        u  = np.column_stack([np.zeros(n_mc), w])
        counts += np.bincount(u.argmax(axis=1), minlength=p)
    return counts / counts.sum()

rng_sub = np.random.default_rng(0)
idx    = rng_sub.choice(n_post, 400, replace=False)
b_sub  = beta_draws_m[idx];  s_sub = sigma_draws_m[idx]

print("\nComputing choice probabilities (≈30 s) ...", flush=True)
p_equal = choice_probs_mc(np.zeros(pm4),                         b_sub, s_sub, pm4)
p_mean  = choice_probs_mc(mean_diffs,                             b_sub, s_sub, pm4)
p_bb10  = choice_probs_mc(mean_diffs + np.array([-0.10, 0, 0]),  b_sub, s_sub, pm4)
p_hse10 = choice_probs_mc(mean_diffs + np.array([0, -0.10, 0]),  b_sub, s_sub, pm4)

hdr = f"{'Brand':<18}  {'Observed':>8}  {'Equal px':>9}  {'Mean px':>8}  {'BB-$0.10':>9}  {'House-$0.10':>12}"
print("\nPosterior-predictive choice probabilities:")
print(hdr); print("-" * len(hdr))
for j, b in enumerate(BRANDS4):
    print(f"{b:<18}  {obs_shares[j]:>8.3f}  {p_equal[j]:>9.3f}  "
          f"{p_mean[j]:>8.3f}  {p_bb10[j]:>9.3f}  {p_hse10[j]:>12.3f}")

fig, ax = plt.subplots(figsize=(11, 4))
x = np.arange(4);  w = 0.16
ax.bar(x-2*w, obs_shares, w, label='Observed',    color='dimgray')
ax.bar(x-  w, p_equal,   w, label='Equal prices', color='steelblue')
ax.bar(x,     p_mean,    w, label='Mean prices',  color='darkorange')
ax.bar(x+  w, p_bb10,   w, label='BB −$0.10',    color='seagreen')
ax.bar(x+2*w, p_hse10,  w, label='House −$0.10', color='orchid')
ax.set_xticks(x); ax.set_xticklabels(BRANDS4, rotation=10, fontsize=9)
ax.set_ylabel('Choice probability')
ax.set_title('M1 — posterior-predictive market shares under price scenarios')
ax.legend(fontsize=9)
plt.tight_layout(); plt.show()

Mean price diffs vs Parkay:  BB=+0.032  House=-0.076  Shedd=+0.326

Computing choice probabilities (≈30 s) ...

Posterior-predictive choice probabilities:
Brand               Observed   Equal px   Mean px   BB-$0.10   House-$0.10
--------------------------------------------------------------------------
Parkay Stk             0.523      0.411     0.407      0.387         0.344
Blue Bonnet Stk        0.207      0.153     0.225      0.343         0.159
House Stk              0.176      0.081     0.250      0.184         0.399
Shedd Tub              0.094      0.356     0.118      0.087         0.098

brand_dummies = np.tile(np.eye(pm4), (N4, 1))
X_marg_m2    = np.hstack([brand_dummies, X_marg])

k_m2 = X_marg_m2.shape[1]
print(f"M2 design matrix: {X_marg_m2.shape}  (N×pm1={N4*pm4}, k={k_m2})")
print(f"beta: [α_BB, α_House, α_Shedd, γ_price]")

prior_m2 = make_prior(pm4, k_m2)

NDRAWS_M2 = 20000
BURN_M2   =  5000

Data_m2  = {'p': p4, 'y': y_marg, 'X': X_marg_m2}
Mcmc_m2  = {'R': NDRAWS_M2, 'keep': 1, 'nprint': 5000,
             'n_sweeps': 5, 'brand_price': True, **make_init(p4, k_m2)}

np.random.seed(2026)   # fix seed for reproducibility
print(f"\nRunning M2 ({NDRAWS_M2} draws, n_sweeps=5, brand_price=True) ...")
t0 = time.time()
out_m2 = gibbs_sampler(Data_m2, prior_m2, Mcmc_m2)
print(f"Done in {time.time()-t0:.1f} s")

M2 design matrix: (10131, 4)  (N×pm1=10131, k=4)
beta: [α_BB, α_House, α_Shedd, γ_price]

Running M2 (20000 draws, n_sweeps=5, brand_price=True) ...
MNP Gibbs: n=3377, p=4, k=4, R=20000, keep=1, n_sweeps=5, brand_price=True
  Iteration   5000 / 20000
  Iteration  10000 / 20000
  Iteration  15000 / 20000
  Iteration  20000 / 20000
Done in 261.1 s

### M2 results

n2       = len(out_m2['betadraw']) - BURN_M2
b2_draws = out_m2['betadraw'][BURN_M2:]     # (n2, 4)
s2_draws = out_m2['sigmadraw'][BURN_M2:]    # (n2, 9)
S2_3d    = s2_draws.reshape(n2, pm4, pm4)

b2_mean  = b2_draws.mean(axis=0)
s2_mean  = S2_3d.mean(axis=0)

PARAM_NAMES_M2 = ['α_BB', 'α_House', 'α_Shedd', 'γ_price']
print("M2 posterior means:")
for name, m, sd in zip(PARAM_NAMES_M2, b2_mean,
                        b2_draws.std(axis=0)):
    ci = np.percentile(b2_draws[:, list(PARAM_NAMES_M2).index(name)], [2.5, 97.5])
    print(f"  {name:<10}  mean={m:+.4f}  SD={sd:.4f}  95%CI=[{ci[0]:+.4f},{ci[1]:+.4f}]")

print(f"\nM1 γ_price: {beta_chain.mean():.4f}   M2 γ_price: {b2_mean[-1]:.4f}")

# Sigma comparison
D2   = np.diag(1/np.sqrt(np.diag(s2_mean)))
c2   = D2 @ s2_mean @ D2

print("\nM2 Sigma* (posterior mean):")
print(pd.DataFrame(s2_mean, index=ALT_LABELS, columns=ALT_LABELS).round(4))
print("\nM2 implied correlations:")
print(pd.DataFrame(c2, index=ALT_LABELS, columns=ALT_LABELS).round(4))

print("\nCorrelation change M1 → M2:")
D1   = np.diag(1/np.sqrt(np.diag(sigma_mat_m)))
c1   = D1 @ sigma_mat_m @ D1
diff = c2 - c1
print(pd.DataFrame(diff, index=ALT_LABELS, columns=ALT_LABELS).round(4))

# ESS for M2 key params
monitor_m2 = {'γ_price': b2_draws[:, -1]}
monitor_m2.update({f'Σ[{ALT_LABELS[r]},{ALT_LABELS[c]}]': S2_3d[:,r,c]
                   for r,c in [(0,1),(0,2),(1,1),(1,2),(2,2)]})
print("\nM2 ESS:")
print(f"{'Parameter':<22}  {'Mean':>8}  {'ESS':>7}  {'ESS%':>7}")
print("-" * 50)
for name, chain in monitor_m2.items():
    e = ess(chain)
    print(f"{name:<22}  {chain.mean():>8.4f}  {e:>7.0f}  {e/n2:>7.2%}")

M2 posterior means:
  α_BB        mean=-0.5570  SD=0.0417  95%CI=[-0.6396,-0.4774]
  α_House     mean=-1.1695  SD=0.0909  95%CI=[-1.3636,-1.0053]
  α_Shedd     mean=-0.0842  SD=0.1540  95%CI=[-0.4050,+0.1877]
  γ_price     mean=-4.4832  SD=0.1941  95%CI=[-4.8734,-4.1086]

M1 γ_price: -1.5179   M2 γ_price: -4.4832

M2 Sigma* (posterior mean):
           BB Stk  House Stk  Shedd Tub
BB Stk     1.0000     0.3810     0.6709
House Stk  0.3810     1.1420     0.2347
Shedd Tub  0.6709     0.2347     2.0093

M2 implied correlations:
           BB Stk  House Stk  Shedd Tub
BB Stk     1.0000     0.3565     0.4733
House Stk  0.3565     1.0000     0.1549
Shedd Tub  0.4733     0.1549     1.0000

Correlation change M1 → M2:
           BB Stk  House Stk  Shedd Tub
BB Stk     0.0000    -0.5738    -0.4613
House Stk -0.5738     0.0000    -0.7681
Shedd Tub -0.4613    -0.7681     0.0000

M2 ESS:
Parameter                   Mean      ESS     ESS%
--------------------------------------------------
γ_price                  -4.4832      564    3.76%
Σ[BB Stk,House Stk]       0.3810      109    0.73%
Σ[BB Stk,Shedd Tub]       0.6709       98    0.65%
Σ[House Stk,House Stk]    1.1420      308    2.06%
Σ[House Stk,Shedd Tub]    0.2347       48    0.32%
Σ[Shedd Tub,Shedd Tub]    2.0093      157    1.05%

# Save M2 margarine results for cross-notebook comparison (MNP_bayesm.ipynb §6)
pd.DataFrame({
    "coef": ["alpha_BB", "alpha_House", "alpha_Shedd", "gamma_price"],
    "beta": b2_mean
}).to_csv("mnp_marg4_Py_beta.csv", index=False)

pd.DataFrame(s2_mean, index=ALT_LABELS, columns=ALT_LABELS
             ).to_csv("mnp_marg4_Py_sigma.csv")

print("Saved mnp_marg4_Py_beta.csv  (4 M2 coefficients)")
print("Saved mnp_marg4_Py_sigma.csv  (3x3 Sigma*)")
print(f"  alpha_BB={b2_mean[0]:.4f}  alpha_House={b2_mean[1]:.4f}  "
      f"alpha_Shedd={b2_mean[2]:.4f}  gamma={b2_mean[3]:.4f}")

Saved mnp_marg4_Py_beta.csv  (4 M2 coefficients)
Saved mnp_marg4_Py_sigma.csv  (3x3 Sigma*)
  alpha_BB=-0.5570  alpha_House=-1.1695  alpha_Shedd=-0.0842  gamma=-4.4832

# Figure — omitted-variable attenuation: the price coefficient M1 vs M2
g_m1 = beta_draws_m.ravel()      # M1 price coefficient draws (price-only model)
g_m2 = b2_draws[:, 3]            # M2 price coefficient draws (with brand intercepts)
fig, ax = plt.subplots(figsize=(8.5, 4.2))
ax.hist(g_m1, bins=50, density=True, color='#f5a623', edgecolor='white', alpha=0.8,
        label=f'M1 (price only):  $\gamma$ = {g_m1.mean():.2f}')
ax.hist(g_m2, bins=50, density=True, color='firebrick', edgecolor='white', alpha=0.8,
        label=f'M2 (+ brand intercepts):  $\gamma$ = {g_m2.mean():.2f}')
ax.axvline(g_m1.mean(), color='#c77f00', lw=1.5); ax.axvline(g_m2.mean(), color='darkred', lw=1.5)
ax.annotate('', xy=(g_m2.mean(), 1.4), xytext=(g_m1.mean(), 1.4),
            arrowprops=dict(arrowstyle='<->', color='gray', lw=1.3))
ax.text((g_m1.mean() + g_m2.mean()) / 2, 1.55,
        f'M2 is {g_m2.mean()/g_m1.mean():.1f}x M1', ha='center', fontsize=9, color='gray')
ax.set_xlabel('price coefficient  $\gamma$'); ax.set_ylabel('posterior density')
ax.set_title('Omitted-variable attenuation, visualised\n'
             'without brand intercepts $\gamma$ must absorb brand preference too, pulling it toward zero')
ax.legend(fontsize=9, loc='upper left'); plt.tight_layout(); plt.show()

obs_shares2 = np.array([(y_marg == j + 1).mean() for j in range(p4)])

X_cube_m2_view = X_marg_m2.reshape(N4, pm4, k_m2)
mean_diffs2    = X_cube_m2_view[:, :, -1].mean(axis=0)   # same as M1 price diffs

print(f"Mean price diffs vs Parkay:  BB={mean_diffs2[0]:+.3f}  "
      f"House={mean_diffs2[1]:+.3f}  Shedd={mean_diffs2[2]:+.3f}")

def choice_probs_mc_m2(price_diffs, beta_sample, sigma_sample, pm1, n_mc=2000, seed=42):
    """beta_sample rows: [alpha_1, ..., alpha_pm1, gamma_price]."""
    rng = np.random.default_rng(seed)
    p = pm1 + 1; counts = np.zeros(p)
    for beta_row, sig_flat in zip(beta_sample, sigma_sample):
        alpha = beta_row[:pm1]
        gamma = beta_row[-1]
        S     = sig_flat.reshape(pm1, pm1)
        mu    = alpha + gamma * np.asarray(price_diffs)
        w     = rng.multivariate_normal(mu, S, size=n_mc)
        u     = np.column_stack([np.zeros(n_mc), w])
        counts += np.bincount(u.argmax(axis=1), minlength=p)
    return counts / counts.sum()

rng_sub2 = np.random.default_rng(1)
idx2   = rng_sub2.choice(n2, 400, replace=False)
b2_sub = b2_draws[idx2]; s2_sub = s2_draws[idx2]

print("\nComputing M2 choice probabilities ...", flush=True)
p2_equal = choice_probs_mc_m2(np.zeros(pm4),                            b2_sub, s2_sub, pm4)
p2_mean  = choice_probs_mc_m2(mean_diffs2,                              b2_sub, s2_sub, pm4)
p2_bb10  = choice_probs_mc_m2(mean_diffs2 + np.array([-0.10, 0, 0]),   b2_sub, s2_sub, pm4)
p2_hse10 = choice_probs_mc_m2(mean_diffs2 + np.array([0, -0.10, 0]),   b2_sub, s2_sub, pm4)
p2_shd10 = choice_probs_mc_m2(mean_diffs2 + np.array([0, 0, -0.10]),   b2_sub, s2_sub, pm4)

hdr = (f"{'Brand':<18}  {'Observed':>8}  {'Equal px':>9}  {'Mean px':>8}  "
       f"{'BB-$0.10':>9}  {'House-$0.10':>12}  {'Shedd-$0.10':>12}")
print("\nM2 — Posterior-predictive choice probabilities:")
print(hdr); print("-" * len(hdr))
for j, b in enumerate(BRANDS4):
    print(f"{b:<18}  {obs_shares2[j]:>8.3f}  {p2_equal[j]:>9.3f}  "
          f"{p2_mean[j]:>8.3f}  {p2_bb10[j]:>9.3f}  {p2_hse10[j]:>12.3f}  "
          f"{p2_shd10[j]:>12.3f}")

fig, ax = plt.subplots(figsize=(13, 4))
x = np.arange(4); w = 0.13
ax.bar(x - 2.5*w, obs_shares2, w, label='Observed',     color='dimgray')
ax.bar(x - 1.5*w, p2_equal,   w, label='Equal prices',  color='steelblue')
ax.bar(x - 0.5*w, p2_mean,    w, label='Mean prices',   color='darkorange')
ax.bar(x + 0.5*w, p2_bb10,   w, label='BB −$0.10',     color='seagreen')
ax.bar(x + 1.5*w, p2_hse10,  w, label='House −$0.10',  color='orchid')
ax.bar(x + 2.5*w, p2_shd10,  w, label='Shedd −$0.10',  color='tomato')
ax.set_xticks(x); ax.set_xticklabels(BRANDS4, rotation=10, fontsize=9)
ax.set_ylabel('Choice probability')
ax.set_title('M2 — posterior-predictive market shares under price scenarios')
ax.legend(fontsize=9)
plt.tight_layout(); plt.show()

Mean price diffs vs Parkay:  BB=+0.032  House=-0.076  Shedd=+0.326

Computing M2 choice probabilities ...

M2 — Posterior-predictive choice probabilities:
Brand               Observed   Equal px   Mean px   BB-$0.10   House-$0.10   Shedd-$0.10
----------------------------------------------------------------------------------------
Parkay Stk             0.523      0.403     0.579      0.480         0.487         0.539
Blue Bonnet Stk        0.207      0.123     0.166      0.313         0.138         0.146
House Stk              0.176      0.064     0.163      0.134         0.292         0.154
Shedd Tub              0.094      0.410     0.092      0.072         0.083         0.161

# Figure — cross-price substitution: where share flows when a brand's price rises
def _cp(price_diffs, bs_, ss_, n_mc=2500, seed=7):
    rng = np.random.default_rng(seed); counts = np.zeros(p4)
    for br, sf in zip(bs_, ss_):
        mu = br[:pm4] + br[pm4] * np.asarray(price_diffs)
        w  = rng.multivariate_normal(mu, sf.reshape(pm4, pm4), size=n_mc)
        counts += np.bincount(np.column_stack([np.zeros(n_mc), w]).argmax(1), minlength=p4)
    return counts / counts.sum()

_idx = np.random.default_rng(0).choice(len(b2_draws), 300, replace=False)
_bs, _ss = b2_draws[_idx], s2_draws[_idx]
BR = ['Parkay', 'Blue Bonnet', 'House', 'Shedd']; DELTA = 0.20
_base = _cp(mean_diffs2, _bs, _ss)
Msub = np.zeros((p4, p4))                                  # row: brand whose price +$0.20 ; col: delta-share
for r in range(p4):
    dd = mean_diffs2.astype(float).copy()
    if r == 0: dd = dd - DELTA                             # Parkay is the base -> every differential shrinks
    else:      dd[r - 1] = dd[r - 1] + DELTA
    Msub[r] = _cp(dd, _bs, _ss) - _base
fig, ax = plt.subplots(figsize=(6.6, 5.4)); vmax = np.abs(Msub).max()
im = ax.imshow(Msub * 100, cmap='RdBu_r', vmin=-vmax*100, vmax=vmax*100)
ax.set_xticks(range(p4)); ax.set_xticklabels(BR, fontsize=9)
ax.set_yticks(range(p4)); ax.set_yticklabels(BR, fontsize=9)
ax.set_xlabel('change in market share of...'); ax.set_ylabel("...when this brand's price rises \$0.20")
for r in range(p4):
    for c in range(p4):
        ax.text(c, r, f'{Msub[r,c]*100:+.1f}', ha='center', va='center', fontsize=8.5,
                color='white' if abs(Msub[r,c]) > vmax*0.55 else 'black')
ax.set_title('Cross-price substitution (M2)\n'
             'a price rise sends buyers to specific rivals, not proportionally (no IIA)\n'
             'rows sum to ~0; the off-diagonal pattern is what logit cannot produce', fontsize=9.5)
fig.colorbar(im, ax=ax, label='delta share (percentage points)', shrink=0.8)
plt.tight_layout(); plt.show()

### Load full margarine data from R CSVs (generated by MNP_bayesm.ipynb)
# If not yet exported, we build it here from scratch using the same logic.

import subprocess, os

# Try loading from R exports first; otherwise build inline
y_all_file = 'mnp_all10_y.csv'
X_all_file = 'mnp_all10_X.csv'

if os.path.exists(y_all_file):
    y_all = pd.read_csv(y_all_file)['y'].values.astype(int)
    X_all = pd.read_csv(X_all_file).values.astype(float)
    print(f"Loaded from CSV: N={len(y_all)}, shape X={X_all.shape}")
else:
    # Build from bayesm margarine data exported as CSV
    # (requires mnp_marg4_y.csv to already exist — re-uses the price columns)
    print("CSV not found — please run the R data-export cell in MNP_bayesm.ipynb first.")
    print("Add a cell there that exports all 10 brands similarly to cell-10.")

Loaded from CSV: N=4470, shape X=(40230, 10)

if 'y_all' in dir() and y_all is not None:
    p_all   = 10
    pm1_all = p_all - 1
    N_all   = len(y_all)
    k_all   = X_all.shape[1]

    BRANDS_ALL = ["Parkay Stk","Blue Bonnet Stk","Fleischmann Stk","House Stk",
                  "Generic Stk","Imperial Stk","Shedd Tub","Parkay Tub",
                  "Fleischmann Tub","House Tub"]

    print(f"M3: N={N_all}, p={p_all}, k={k_all}, pm1={pm1_all}")
    obs_all = np.array([(y_all == j+1).mean() for j in range(p_all)])
    for b, s in zip(BRANDS_ALL, obs_all):
        print(f"  {b:<22}: {s:.3f}")

    prior_m3 = make_prior(pm1_all, k_all)
    NDRAWS_M3 = 50000
    BURN_M3   =  5000

    Data_m3  = {'p': p_all, 'y': y_all, 'X': X_all}
    Mcmc_m3  = {'R': NDRAWS_M3, 'keep': 1, 'nprint': 10000,
                 'n_sweeps': 5, 'brand_price': True,
                 **make_init(p_all, k_all)}

    np.random.seed(2027)   # fix seed for reproducibility
    print(f"\nRunning M3 ({NDRAWS_M3} draws, n_sweeps=5, brand_price=True) ...")
    t0 = time.time()
    out_m3 = gibbs_sampler(Data_m3, prior_m3, Mcmc_m3)
    elapsed_m3 = time.time() - t0
    print(f"Done in {elapsed_m3:.1f} s  ({elapsed_m3/60:.1f} min)")
else:
    print("Load data first (cell above).")

M3: N=4470, p=10, k=10, pm1=9
  Parkay Stk            : 0.395
  Blue Bonnet Stk       : 0.156
  Fleischmann Stk       : 0.054
  House Stk             : 0.133
  Generic Stk           : 0.070
  Imperial Stk          : 0.017
  Shedd Tub             : 0.071
  Parkay Tub            : 0.045
  Fleischmann Tub       : 0.050
  House Tub             : 0.007

Running M3 (50000 draws, n_sweeps=5, brand_price=True) ...
MNP Gibbs: n=4470, p=10, k=10, R=50000, keep=1, n_sweeps=5, brand_price=True
  Iteration  10000 / 50000
  Iteration  20000 / 50000
  Iteration  30000 / 50000
  Iteration  40000 / 50000
  Iteration  50000 / 50000
Done in 2444.3 s  (40.7 min)

### M3 results

b3_draws = out_m3['betadraw'][BURN_M3:]        # (n3, 10)
s3_draws = out_m3['sigmadraw'][BURN_M3:]       # (n3, 81)
n3       = len(b3_draws)
S3_3d    = s3_draws.reshape(n3, pm1_all, pm1_all)

b3_mean  = b3_draws.mean(axis=0)
s3_mean  = S3_3d.mean(axis=0)

ALT_LABELS_M3  = BRANDS_ALL[1:]               # 9 non-base brand names
PARAM_NAMES_M3 = [f'α_{b}' for b in ALT_LABELS_M3] + ['γ_price']

print(f"Post-burn-in draws : {n3:,}")
print(f"sigma_11 (should=1): {s3_mean[0,0]:.6f}")
print(f"\nM3 posterior means:")
print(f"  {'Parameter':<28}  {'Mean':>8}  {'SD':>8}  {'95% CI'}")
print("-" * 70)
for i, name in enumerate(PARAM_NAMES_M3):
    m  = b3_mean[i]
    sd = b3_draws[:, i].std()
    ci = np.percentile(b3_draws[:, i], [2.5, 97.5])
    print(f"  {name:<28}  {m:>+8.4f}  {sd:>8.4f}  [{ci[0]:+.4f}, {ci[1]:+.4f}]")

# ESS — gamma + all 9 Sigma diagonals
monitor_m3 = {'γ_price': b3_draws[:, -1]}
for j in range(pm1_all):
    monitor_m3[f'Σ[{ALT_LABELS_M3[j]},{ALT_LABELS_M3[j]}]'] = S3_3d[:, j, j]

print(f"\nM3 ESS (γ and Σ diagonals):")
print(f"  {'Parameter':<38}  {'Mean':>8}  {'ESS':>7}  {'ESS%':>7}")
print("  " + "-" * 65)
for name, chain in monitor_m3.items():
    e = ess(chain)
    print(f"  {name:<38}  {chain.mean():>8.4f}  {e:>7.0f}  {e/n3:>7.2%}")

Post-burn-in draws : 45,000
sigma_11 (should=1): 1.000000

M3 posterior means:
  Parameter                         Mean        SD  95% CI
----------------------------------------------------------------------
  α_Blue Bonnet Stk              -0.6188    0.0492  [-0.7149, -0.5192]
  α_Fleischmann Stk              -2.0644    0.5468  [-3.2802, -1.1291]
  α_House Stk                    -1.2029    0.1067  [-1.4158, -0.9988]
  α_Generic Stk                  -2.0944    0.1976  [-2.5360, -1.7474]
  α_Imperial Stk                 -2.6688    0.6792  [-4.1765, -1.5752]
  α_Shedd Tub                    -0.3069    0.1657  [-0.6775, -0.0207]
  α_Parkay Tub                   -0.3621    0.4161  [-1.3274, +0.3201]
  α_Fleischmann Tub              -1.9039    0.7891  [-3.5850, -0.5457]
  α_House Tub                    -3.0866    0.5028  [-4.2509, -2.2625]
  γ_price                        -4.1222    0.2375  [-4.5872, -3.6696]

M3 ESS (γ and Σ diagonals):
  Parameter                                   Mean      ESS     ESS%
  -----------------------------------------------------------------
  γ_price                                  -4.1222      209    0.46%
  Σ[Blue Bonnet Stk,Blue Bonnet Stk]        1.0000    43575   96.83%
  Σ[Fleischmann Stk,Fleischmann Stk]        9.3652       78    0.17%
  Σ[House Stk,House Stk]                    1.2489      284    0.63%
  Σ[Generic Stk,Generic Stk]                1.8546      160    0.36%
  Σ[Imperial Stk,Imperial Stk]              5.3153       62    0.14%
  Σ[Shedd Tub,Shedd Tub]                    2.2944      149    0.33%
  Σ[Parkay Tub,Parkay Tub]                  3.7109       70    0.16%
  Σ[Fleischmann Tub,Fleischmann Tub]       11.6513       64    0.14%
  Σ[House Tub,House Tub]                    2.9766       69    0.15%

# Figure - M3 brand preference forest plot (9 intercepts vs the Parkay Stick base)
amean = b3_draws[:, :9].mean(axis=0)
alo   = np.percentile(b3_draws[:, :9],  2.5, axis=0)
ahi   = np.percentile(b3_draws[:, :9], 97.5, axis=0)
BR9 = ['Blue Bonnet Stk', 'Fleischmann Stk', 'House Stk', 'Generic Stk', 'Imperial Stk',
       'Shedd Tub', 'Parkay Tub', 'Fleischmann Tub', 'House Tub']
labels = ['Parkay Stk (base)'] + BR9
means  = np.concatenate([[0.0], amean])
los    = np.concatenate([[0.0], alo])
his    = np.concatenate([[0.0], ahi])
order  = np.argsort(means)                                    # least preferred at the bottom
fig, ax = plt.subplots(figsize=(8, 5.2))
for yi, k in enumerate(order):
    incl0 = (los[k] < 0 < his[k]) or (k == 0)                 # CI covers Parkay (or is the base)
    col = '#c0392b' if incl0 else 'steelblue'
    ax.plot([los[k], his[k]], [yi, yi], color=col, lw=2.2, alpha=0.85, solid_capstyle='round')
    ax.plot(means[k], yi, 'o', color=col, ms=6, zorder=3)
ax.axvline(0, color='black', lw=1, ls='--')
ax.set_yticks(range(len(order))); ax.set_yticklabels([labels[k] for k in order], fontsize=8)
ax.set_xlabel(r'brand intercept  $\alpha_j$   (utility vs Parkay Stick at equal prices)')
ax.set_title('M3 - brand preference ranking, all 10 brands\n'
             'red = 95% CI covers Parkay (statistically indistinguishable); bars widen as brand share shrinks')
plt.tight_layout(); plt.show()

### M3 — Sigma* heatmaps (covariance and correlation)

D3     = np.diag(1.0 / np.sqrt(np.diag(s3_mean)))
corr3  = D3 @ s3_mean @ D3

SHORT = ["BB Stk","Fl Stk","Hse Stk","Gen Stk","Imp Stk",
         "Shdd Tub","Pk Tub","Fl Tub","Hse Tub"]

print("M3 posterior mean Σ*:")
print(pd.DataFrame(s3_mean, index=SHORT, columns=SHORT).round(3).to_string())
print("\nM3 implied correlations:")
print(pd.DataFrame(corr3, index=SHORT, columns=SHORT).round(3).to_string())

fig, axes = plt.subplots(1, 2, figsize=(16, 6))
for ax, mat, title in zip(axes, [s3_mean, corr3],
                           ['Σ* (covariance)', 'Implied correlations ρ']):
    im = ax.imshow(mat, cmap='RdBu_r', vmin=-mat.max(), vmax=mat.max(), aspect='auto')
    ax.set_xticks(range(9)); ax.set_xticklabels(SHORT, rotation=45, ha='right', fontsize=8)
    ax.set_yticks(range(9)); ax.set_yticklabels(SHORT, fontsize=8)
    for r in range(9):
        for c in range(9):
            ax.text(c, r, f'{mat[r,c]:.2f}', ha='center', va='center',
                    fontsize=6.5, color='black')
    ax.set_title(title, fontsize=11)
    plt.colorbar(im, ax=ax, shrink=0.8)

plt.suptitle('M3 — posterior mean Σ* (all 10 brands)', fontsize=12)
plt.tight_layout(); plt.show()

M3 posterior mean Σ*:
          BB Stk  Fl Stk  Hse Stk  Gen Stk  Imp Stk  Shdd Tub  Pk Tub  Fl Tub  Hse Tub
BB Stk     1.000  -1.290    0.425    0.529    0.613     0.741   0.352  -1.215    0.602
Fl Stk    -1.290   9.365   -1.467   -2.247   -1.728    -3.138   0.369   8.700   -1.839
Hse Stk    0.425  -1.467    1.249    0.825    0.496     0.625  -0.277  -1.411    0.355
Gen Stk    0.529  -2.247    0.825    1.855    0.873     1.242  -0.431  -2.325    0.647
Imp Stk    0.613  -1.728    0.496    0.873    5.315     1.153   0.611  -1.629    1.913
Shdd Tub   0.741  -3.138    0.625    1.242    1.153     2.294  -0.017  -3.260    1.173
Pk Tub     0.352   0.369   -0.277   -0.431    0.611    -0.017   3.711   0.812    0.301
Fl Tub    -1.215   8.700   -1.411   -2.325   -1.629    -3.260   0.812  11.651   -1.841
Hse Tub    0.602  -1.839    0.355    0.647    1.913     1.173   0.301  -1.841    2.977

M3 implied correlations:
          BB Stk  Fl Stk  Hse Stk  Gen Stk  Imp Stk  Shdd Tub  Pk Tub  Fl Tub  Hse Tub
BB Stk     1.000  -0.422    0.380    0.389    0.266     0.489   0.183  -0.356    0.349
Fl Stk    -0.422   1.000   -0.429   -0.539   -0.245    -0.677   0.063   0.833   -0.348
Hse Stk    0.380  -0.429    1.000    0.542    0.193     0.369  -0.129  -0.370    0.184
Gen Stk    0.389  -0.539    0.542    1.000    0.278     0.602  -0.164  -0.500    0.275
Imp Stk    0.266  -0.245    0.193    0.278    1.000     0.330   0.138  -0.207    0.481
Shdd Tub   0.489  -0.677    0.369    0.602    0.330     1.000  -0.006  -0.630    0.449
Pk Tub     0.183   0.063   -0.129   -0.164    0.138    -0.006   1.000   0.123    0.090
Fl Tub    -0.356   0.833   -0.370   -0.500   -0.207    -0.630   0.123   1.000   -0.313
Hse Tub    0.349  -0.348    0.184    0.275    0.481     0.449   0.090  -0.313    1.000

X_cube_m3  = X_all.reshape(N_all, pm1_all, k_all)
mean_diffs3 = X_cube_m3[:, :, -1].mean(axis=0)   # (9,) mean price diffs vs Parkay Stk

print("Mean price diffs vs Parkay Stk:")
for b, d in zip(ALT_LABELS_M3, mean_diffs3):
    print(f"  {b:<22}: {d:+.3f}")

obs_shares3 = np.array([(y_all == j + 1).mean() for j in range(p_all)])

# Reuse choice_probs_mc_m2 — structure is identical (alpha + gamma * price_diffs)
rng_sub3 = np.random.default_rng(2)
idx3     = rng_sub3.choice(n3, 400, replace=False)
b3_sub   = b3_draws[idx3]; s3_sub = s3_draws[idx3]

# Price-cut deltas (index into 9-element price_diffs vector)
# 5: Shedd Tub,  6: Parkay Tub,  2: House Stk
shdd_cut = np.zeros(pm1_all); shdd_cut[5] = -0.10
pktub_cut = np.zeros(pm1_all); pktub_cut[6] = -0.10
hse_cut  = np.zeros(pm1_all); hse_cut[2]  = -0.10

print("\nComputing M3 choice probabilities ...", flush=True)
p3_equal  = choice_probs_mc_m2(np.zeros(pm1_all),           b3_sub, s3_sub, pm1_all, n_mc=1000)
p3_mean   = choice_probs_mc_m2(mean_diffs3,                  b3_sub, s3_sub, pm1_all, n_mc=1000)
p3_shdd10 = choice_probs_mc_m2(mean_diffs3 + shdd_cut,       b3_sub, s3_sub, pm1_all, n_mc=1000)
p3_pktub  = choice_probs_mc_m2(mean_diffs3 + pktub_cut,      b3_sub, s3_sub, pm1_all, n_mc=1000)
p3_hse10  = choice_probs_mc_m2(mean_diffs3 + hse_cut,        b3_sub, s3_sub, pm1_all, n_mc=1000)

hdr = (f"{'Brand':<22}  {'Obs':>6}  {'Equal':>7}  {'Mean':>7}  "
       f"{'Shdd-$.10':>10}  {'PkTub-$.10':>11}  {'Hse-$.10':>9}")
print("\nM3 — Posterior-predictive choice probabilities:")
print(hdr); print("-" * len(hdr))
for j, b in enumerate(BRANDS_ALL):
    print(f"{b:<22}  {obs_shares3[j]:>6.3f}  {p3_equal[j]:>7.3f}  {p3_mean[j]:>7.3f}  "
          f"{p3_shdd10[j]:>10.3f}  {p3_pktub[j]:>11.3f}  {p3_hse10[j]:>9.3f}")

fig, ax = plt.subplots(figsize=(15, 5))
x = np.arange(p_all); w = 0.13
ax.bar(x - 2.5*w, obs_shares3,  w, label='Observed',       color='dimgray')
ax.bar(x - 1.5*w, p3_equal,     w, label='Equal prices',   color='steelblue')
ax.bar(x - 0.5*w, p3_mean,      w, label='Mean prices',    color='darkorange')
ax.bar(x + 0.5*w, p3_shdd10,   w, label='Shedd −$0.10',   color='tomato')
ax.bar(x + 1.5*w, p3_pktub,    w, label='PkTub −$0.10',   color='mediumpurple')
ax.bar(x + 2.5*w, p3_hse10,    w, label='House −$0.10',   color='orchid')
ax.set_xticks(x)
ax.set_xticklabels([b.replace(' ', '\n') for b in BRANDS_ALL], fontsize=8)
ax.set_ylabel('Choice probability')
ax.set_title('M3 — posterior-predictive market shares under price scenarios')
ax.legend(fontsize=8, loc='upper right')
plt.tight_layout(); plt.show()

Mean price diffs vs Parkay Stk:
  Blue Bonnet Stk       : +0.025
  Fleischmann Stk       : +0.497
  House Stk             : -0.081
  Generic Stk           : -0.173
  Imperial Stk          : +0.262
  Shedd Tub             : +0.307
  Parkay Tub            : +0.559
  Fleischmann Tub       : +0.671
  House Tub             : +0.050

Computing M3 choice probabilities ...

M3 — Posterior-predictive choice probabilities:
Brand                      Obs    Equal     Mean   Shdd-$.10   PkTub-$.10   Hse-$.10
------------------------------------------------------------------------------------
Parkay Stk               0.395    0.120    0.428       0.404        0.414      0.370
Blue Bonnet Stk          0.156    0.053    0.127       0.115        0.122      0.110
Fleischmann Stk          0.054    0.101    0.055       0.055        0.054      0.054
House Stk                0.133    0.033    0.126       0.118        0.124      0.228
Generic Stk              0.070    0.005    0.069       0.059        0.068      0.056
Imperial Stk             0.017    0.030    0.017       0.015        0.016      0.015
Shedd Tub                0.071    0.251    0.075       0.134        0.073      0.067
Parkay Tub               0.045    0.234    0.046       0.044        0.073      0.044
Fleischmann Tub          0.050    0.171    0.052       0.052        0.050      0.051
House Tub                0.007    0.002    0.006       0.005        0.006      0.005

Field	Detail
Unit	Individual purchase occasion
Households	516
Purchase occasions	4,470 (≈ 8.7 per household)
Brands / SKUs	10 margarine products
Covariates	Transaction price for each brand at each purchase occasion
Demographics	Household income, family size, education, white-collar, retired (in `demos`)

Code	Brand	Format
PPk_Stk	Parkay	Stick
PBB_Stk	Blue Bonnet	Stick
PFl_Stk	Fleischmann's	Stick
PHse_Stk	House Brand	Stick
PGen_Stk	Generic	Stick
PImp_Stk	Imperial	Stick
PSS_Tub	Shedd's Spread	Tub
PPk_Tub	Parkay	Tub
PFl_Tub	Fleischmann's	Tub
PHse_Tub	House Brand	Tub

Model	Brands (p)	Covariates (k)
M1 — price only	4 top brands	1
M2 — intercepts + price	4 top brands	4
M3 — intercepts + price	All 10 brands	10

Symbol	Dimension	Description
$\gamma$	scalar	Price sensitivity (expected $< 0$)
$\boldsymbol{\Sigma}^*$	$3 \times 3$	Covariance of utility differences; identified by fixing $\sigma_{11} = 1$

	BB Stk	House Stk	Shedd Tub
BB Stk	1.000 (fixed)	0.661	1.220
House Stk	0.661	0.505	0.856
Shedd Tub	1.220	0.856	1.703

Price increase	$\Delta\mu$	Fraction of $\sigma_{\gamma}$
+$0.05	$-0.076$	$0.8\,\sigma_\gamma$
+$0.10	$-0.152$	$1.6\,\sigma_\gamma$
+$0.20	$-0.304$	$3.2\,\sigma_\gamma$

Brand	$\sigma_{jj}^*$	Relative to BB
BB Stick	1.000	reference (fixed)
House Stick	0.505	lowest — most predictable
Shedd Tub	1.703	highest — most heterogeneous

Brand loses	pp lost	Share of diverted volume
House Stk	−6.6	56%
Shedd Tub	−3.1	26%
Parkay	−2.0	17%

Brand	$\alpha_j$	95% CI	Reading
BB Stick	−0.557	[−0.640, −0.477]	Clearly less preferred; CI well away from zero
House Brand	−1.170	[−1.364, −1.005]	Most penalised; strong private-label discount
Shedd Tub	−0.084	[−0.405, +0.188]	CI spans zero — not significantly different from Parkay

Pair	M1 $\rho$	M2 $\rho$	Change
BB – House	0.930	0.357	−0.574
BB – Shedd	0.935	0.473	−0.461
House – Shedd	0.923	0.155	−0.768

Scenario	Description
Equal prices	$\Delta\mathbf{p} = \mathbf{0}$ — pure brand preference, no price effect
Mean prices	Actual average price differentials in the data
BB −$0.10	Blue Bonnet cuts price by $0.10 vs Parkay
House −$0.10	House Brand cuts price by $0.10 vs Parkay
Shedd −$0.10	Shedd cuts price by $0.10 vs Parkay — interesting given $\alpha_\text{Shedd} \approx 0$

Brand	Observed	M1 mean-px	M2 mean-px
Parkay	52.3%	40.7%	57.9%
BB Stick	20.7%	22.5%	16.6%
House Stick	17.6%	25.0%	16.3%
Shedd Tub	9.4%	11.8%	9.2%

Brand loses	pp lost	Share of diverted volume
Parkay	−9.9	67%
House Stk	−2.9	20%
Shedd Tub	−2.0	14%

Brand loses	pp lost	Share of diverted volume
Parkay	−4.0	58%
BB Stk	−2.0	29%
House Stk	−0.9	13%

Brand	$\alpha_j$	95% CI	Status
Shedd Tub	−0.307	[−0.678, −0.021]	CI barely excludes zero
Parkay Tub	−0.362	[−1.327, +0.320]	CI spans zero — format-neutral
Blue Bonnet Stk	−0.619	[−0.715, −0.519]	Clearly negative; tightest CI
House Stk	−1.203	[−1.416, −0.999]	Strong private-label penalty
Fleischmann Stk	−2.064	[−3.280, −1.129]	Clearly below Parkay; very wide CI
Fleischmann Tub	−1.904	[−3.585, −0.546]	Very wide — only 5% share
Generic Stk	−2.094	[−2.536, −1.747]	Deep discount brand; well-identified
Imperial Stk	−2.669	[−4.177, −1.575]	Low share (1.7%); wide CI
House Tub	−3.087	[−4.251, −2.263]	Near-zero share (0.7%)

Price cut	M1: Parkay share of diversion	M2: Parkay share of diversion
BB $-\$0.10$	17%	67%
House $-\$0.10$	42%	71%

Brand	Short run $\alpha$	Long run $\alpha$	Shift
Fleischmann Stk	−1.645	−2.064	−0.42
Imperial Stk	−3.589	−2.669	+0.92

Parameter	Short ESS	Long ESS	ESS %
$\gamma^*_{\text{price}}$	34	209	0.46
$\sigma^*_{jj}$ Hse Stk	27	284	0.63
$\sigma^*_{jj}$ Gen Stk	19	160	0.36
$\sigma^*_{jj}$ Shdd Tub	22	149	0.33
$\sigma^*_{jj}$ Fl Stk	23	78	0.17
$\sigma^*_{jj}$ Pk Tub	18	70	0.16
$\sigma^*_{jj}$ Hse Tub	12	69	0.15
$\sigma^*_{jj}$ Imp Stk	13	62	0.14
$\sigma^*_{jj}$ Fl Tub	11	64	0.14

Brand	Observed	M3 mean-px	Diff
Parkay Stk	39.5%	42.8%	+3.3
Blue Bonnet Stk	15.6%	12.7%	−2.9
Fleischmann Stk	5.4%	5.5%	+0.1
House Stk	13.3%	12.6%	−0.7
Generic Stk	7.0%	6.9%	−0.1
Imperial Stk	1.7%	1.7%	0.0
Shedd Tub	7.1%	7.5%	+0.4
Parkay Tub	4.5%	4.6%	+0.1
Fleischmann Tub	5.0%	5.2%	+0.2
House Tub	0.7%	0.6%	−0.1

Pair	Short-run $\rho$	Long-run $\rho$
BB – Shdd Tub	+0.498	+0.489
Gen – Shdd Tub	+0.645	+0.602
Gen – Hse Stk	+0.577	+0.542
BB – Gen Stk	+0.445	+0.389
BB – Hse Stk	+0.390	+0.380

Pair	Short-run $\rho$	Long-run $\rho$
Imp Stk – Fl Stk	+0.838	−0.245
Imp Stk – Fl Tub	+0.800	−0.207
Imp Stk – BB Stk	−0.265	+0.266
Imp Stk – Gen Stk	−0.422	+0.278
Imp Stk – Shdd Tub	−0.547	+0.330

Pair	Short-run $\rho$	Long-run $\rho$
Fl Stk – Shdd Tub	−0.673	−0.677
Fl Tub – Shdd Tub	−0.608	−0.630
Fl Stk – Gen Stk	−0.576	−0.539
Fl Tub – Gen Stk	−0.532	−0.500

Brand loses	pp lost	Share of diverted
Parkay Stk	−2.4	41%
Blue Bonnet Stk	−1.2	20%
Generic Stk	−1.0	17%
House Stk	−0.8	14%
Others	−0.5	8%

Brand loses	pp lost	Share of diverted
Parkay Stk	−1.4	52%
Blue Bonnet Stk	−0.5	19%
Shedd Tub	−0.2	7%
Others	−0.6	22%

Brand loses	pp lost	Share of diverted
Parkay Stk	−5.8	57%
Blue Bonnet Stk	−1.7	17%
Generic Stk	−1.3	13%
Shedd Tub	−0.8	8%
Others	−0.6	6%

Scenario	Parkay pp lost	Parkay share of diversion
House Stk $-\$0.10$	−5.8	57%
Shedd Tub $-\$0.10$	−2.4	41%
Parkay Tub $-\$0.10$ (own brand)	−1.4	52%

Multinomial Probit — Python Gibbs Sampler¶

Model¶

Gibbs Sampler Algorithm¶

Block 1 — Latent utility differences $\mathbf{W}$ (sequential truncated normal)¶

Block 2 — Regression coefficients $\boldsymbol{\beta}$¶

Block 3 — Covariance matrix $\boldsymbol{\Sigma}$¶

1. Imports¶

2. Synthetic Data (generated by MNP_bayesm.ipynb)¶

3. Run Python Gibbs — Synthetic Data¶

Synthetic recovery — read it the right way¶

4. Margarine Data — Description and Setup¶

Data source¶

Structure¶

The 10 brands¶

Treatment¶

Models estimated¶

5. M1 — Price-Only Model (Top 4 Brands)¶

Model specification¶

Parameters¶

Identification¶

Limitation¶

M1 — Initial reading¶

6. MCMC Diagnostics¶

M1 — MCMC mixing assessment¶

7a. Price coefficient¶

What does $\gamma^*$ measure?¶

Interpretation¶

Interpretation¶

Interpretation¶

8. M2 — Brand Intercepts + Price (Top 4 Brands)¶

Interpretation¶

The attenuation, in one picture¶

8b. Choice probabilities — posterior predictive (M2)¶

What logit cannot do — flexible substitution¶

Interpretation¶

9. M3 — All 10 Brands + Intercepts + Price¶

Interpretation¶

The ranking at a glance¶

Interpretation¶

9b. Choice probabilities — posterior predictive (M3)¶

Interpretation¶