Seemingly Unrelated Regressions — bayesm (rsurGibbs)¶
Reference: Zellner, A. (1962). An efficient method of estimating seemingly unrelated regressions. Journal of the American Statistical Association, 57, 348–368.
Sampler: Rossi, P.E., Allenby, G.M. & McCulloch, R. (2005). Bayesian Statistics and Marketing. Wiley. Implemented as bayesm::rsurGibbs.
Model¶
M equations, N observations each: $$y_j = X_j \boldsymbol{\beta}_j + \boldsymbol{\varepsilon}_j, \qquad j = 1,\ldots,M$$
Cross-equation correlated errors: $\boldsymbol{\varepsilon}_t = (\varepsilon_{1t},\ldots,\varepsilon_{Mt})^\top \sim N_M(\mathbf{0}, \boldsymbol{\Omega})$
Priors (independent): $$\boldsymbol{\beta} \sim N(\bar{\boldsymbol{\beta}},\, \mathbf{A}^{-1}), \qquad \boldsymbol{\Omega} \sim \mathcal{IW}(\nu,\, \mathbf{V})$$
Full conditionals — 2-block Gibbs:
(i) $\boldsymbol{\beta} \mid \boldsymbol{\Omega}, Y$: Gaussian with precision $\mathbf{Q} = \mathbf{A} + \sum_{i,j} \omega^{ij} X_i^\top X_j$ and right-hand side $\mathbf{b} = \mathbf{A}\bar{\boldsymbol{\beta}} + \sum_{i,j} \omega^{ij} X_i^\top y_j$
(ii) $\boldsymbol{\Omega} \mid \boldsymbol{\beta}, Y$: $\mathcal{IW}(\nu+N,\; \mathbf{V}+\mathbf{E}^\top\mathbf{E})$ where $E_{tj} = y_{jt} - x_{jt}^\top\boldsymbol{\beta}_j$
Notebooks¶
| File | Language | Content |
|---|---|---|
| sur_bayesm.ipynb (this) | R | rsurGibbs — basic SUR, synthetic + tuna data |
| sur_gibbs.py | Python | Gibbs sampler module |
| sur_python.ipynb | Python | Python replication |
1. Setup¶
In [1]:
userLib <- file.path(Sys.getenv("USERPROFILE"), "R", "win-library", "4.6")
.libPaths(c(userLib, .libPaths()))
library(bayesm)
library(MASS)
cat("bayesm version:", as.character(packageVersion("bayesm")), "\n")
cat("Working dir:", getwd(), "\n")
bayesm version: 3.1.7 Working dir: c:/Users/user/project
2. Synthetic Data¶
Simulate $M=3$ equations, $k=2$ regressors each (intercept + continuous $x$), $N=200$ observations. True $\boldsymbol{\beta}$ parameters and true $\boldsymbol{\Omega}$ ($3\times3$ positive-definite) are known, so we can check recovery.
Save $Y$ and $X$ CSVs for the Python replication notebook.
In [2]:
set.seed(2025L)
M_s <- 3L; k_s <- 2L; N_s <- 200L; K_s <- M_s * k_s
# True parameters: k x M (column j = beta_j)
B_TRUE <- matrix(c(5.0, -1.0,
4.0, -1.5,
3.0, -2.0), k_s, M_s)
rownames(B_TRUE) <- c("intercept", "slope")
colnames(B_TRUE) <- paste0("eq", seq_len(M_s))
# True Omega: 3x3 positive definite
OMEGA_TRUE <- matrix(c(1.0, 0.6, 0.3,
0.6, 1.2, 0.4,
0.3, 0.4, 0.8), M_s, M_s)
# Simulate
X_sim <- cbind(1, rnorm(N_s))
E_sim <- mvrnorm(N_s, rep(0, M_s), OMEGA_TRUE) # N x M correlated errors
regdata_sim <- lapply(seq_len(M_s), function(j)
list(y = drop(X_sim %*% B_TRUE[, j]) + E_sim[, j],
X = X_sim))
cat(sprintf("Synthetic SUR: M=%d k=%d N=%d\n", M_s, k_s, N_s))
cat("True B (rows=params, cols=equations):\n")
print(B_TRUE)
cat("\nTrue Omega:\n")
print(round(OMEGA_TRUE, 3))
# Save for Python notebook
Y_sim_mat <- do.call(cbind, lapply(regdata_sim, `[[`, "y")) # N x M
write.csv(Y_sim_mat, "sur_synth_Y.csv", row.names = FALSE)
write.csv(as.data.frame(X_sim), "sur_synth_X.csv", row.names = FALSE)
cat("\nSaved: sur_synth_Y.csv sur_synth_X.csv\n")
Synthetic SUR: M=3 k=2 N=200
True B (rows=params, cols=equations):
eq1 eq2 eq3
intercept 5 4.0 3
slope -1 -1.5 -2
True Omega:
[,1] [,2] [,3]
[1,] 1.0 0.6 0.3
[2,] 0.6 1.2 0.4
[3,] 0.3 0.4 0.8
Saved: sur_synth_Y.csv sur_synth_X.csv
3. rsurGibbs — Synthetic Data¶
Default diffuse prior: $\bar{\boldsymbol{\beta}} = \mathbf{0}$, $\mathbf{A} = 0.01\,\mathbf{I}_K$ (variance 100), $\nu = M+3$, $\mathbf{V} = \nu\,\mathbf{I}_M$ (bayesm default).
rsurGibbs returns $betadraw (R × K) and $Sigmadraw (R × M²), where M² entries are vec(Sigma) in column-major order.
Note (bayesm 3.1.7 bug): when
Prior$Vis supplied, bayesm assigns it to an internal variablessqinstead ofV, leavingVundefined and crashing atprint(V). The workaround is to omitVfrom the Prior list and rely on the default.
In [3]:
Prior_s <- list(betabar = rep(0, K_s),
A = 0.01 * diag(K_s),
nu = M_s + 3L)
# V omitted: bayesm default is V = nu * diag(M), which is diffuse and avoids a 3.1.7 bug
Mcmc_s <- list(R = 20000L, keep = 1L, nprint = 5000L)
set.seed(42L)
cat("Running rsurGibbs on synthetic data ...\n")
t0 <- proc.time()["elapsed"]
out_s <- rsurGibbs(Data = list(regdata = regdata_sim),
Prior = Prior_s,
Mcmc = Mcmc_s)
cat(sprintf("Done in %.1f s\n", proc.time()["elapsed"] - t0))
cat("Output names:", paste(names(out_s), collapse = ", "), "\n")
cat(sprintf("betadraw: %d x %d Sigmadraw: %d x %d\n",
nrow(out_s$betadraw), ncol(out_s$betadraw),
nrow(out_s$Sigmadraw), ncol(out_s$Sigmadraw)))
Running rsurGibbs on synthetic data ...
Starting Gibbs Sampler for SUR Regression Model
with 3 regressions
and 200 observations for each regression
Prior Parms:
betabar
[1] 0 0 0 0 0 0
A
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0.01 0.00 0.00 0.00 0.00 0.00
[2,] 0.00 0.01 0.00 0.00 0.00 0.00
[3,] 0.00 0.00 0.01 0.00 0.00 0.00
[4,] 0.00 0.00 0.00 0.01 0.00 0.00
[5,] 0.00 0.00 0.00 0.00 0.01 0.00
[6,] 0.00 0.00 0.00 0.00 0.00 0.01
nu = 6
V =
[,1] [,2] [,3]
[1,] 6 0 0
[2,] 0 6 0
[3,] 0 0 6
MCMC parms:
R= 20000 keep= 1 nprint= 5000
MCMC Iteration (est time to end - min)
5000 (0.0)
10000 (0.0)
15000 (0.0)
20000 (0.0)
Total Time Elapsed: 0.00
Done in 0.1 s
Output names: betadraw, Sigmadraw
betadraw: 20000 x 6 Sigmadraw: 20000 x 9
In [4]:
BURN_S <- 5000L
keep_s <- seq(BURN_S + 1L, 20000L)
bdraw_s <- out_s$betadraw[keep_s, ] # (15000, 6)
odraw_s <- out_s$Sigmadraw[keep_s, ] # (15000, 9)
B_pm_s <- matrix(colMeans(bdraw_s), k_s, M_s) # k x M
Omega_pm_s <- matrix(colMeans(odraw_s), M_s, M_s) # M x M (column-major reconstruction)
rownames(B_pm_s) <- c("intercept", "slope")
colnames(B_pm_s) <- paste0("eq", seq_len(M_s))
cat("Posterior mean B:\n")
print(round(B_pm_s, 4))
cat("\nSlope recovery (true vs posterior):\n")
cat(sprintf(" %-6s %7s %8s %8s %8s %8s %s\n",
"Eq", "True", "Mean", "Diff", "CI_lo", "CI_hi", "Covered"))
cat(" ", paste(rep("-", 65), collapse = ""), "\n", sep = "")
for (j in seq_len(M_s)) {
col_j <- (j - 1L) * k_s + 2L # slope is 2nd param in each k-block
d <- bdraw_s[, col_j]
lo <- quantile(d, 0.025); hi <- quantile(d, 0.975)
tr <- B_TRUE[2, j]
covered <- if (tr >= lo && tr <= hi) "YES" else "NO"
cat(sprintf(" eq%-4d %7.3f %8.4f %+8.4f %8.4f %8.4f %s\n",
j, tr, mean(d), mean(d) - tr, lo, hi, covered))
}
cat("\nOmega recovery (upper triangle):\n")
cat(sprintf(" %-12s %7s %8s %8s\n", "Element", "True", "PostMean", "Diff"))
cat(" ", paste(rep("-", 45), collapse = ""), "\n", sep = "")
for (i in seq_len(M_s)) {
for (j in i:M_s) {
tr <- OMEGA_TRUE[i, j]; pm <- Omega_pm_s[i, j]
cat(sprintf(" Omega[%d,%d] %7.3f %8.4f %+8.4f\n",
i, j, tr, pm, pm - tr))
}
}
# Save for Python comparison
write.csv(as.data.frame(t(B_pm_s)), "sur_synth_bayesm_B.csv", row.names = FALSE)
write.csv(as.data.frame(Omega_pm_s), "sur_synth_bayesm_Sigma.csv", row.names = FALSE)
cat("\nSaved: sur_synth_bayesm_B.csv sur_synth_bayesm_Sigma.csv\n")
Posterior mean B:
eq1 eq2 eq3
intercept 4.9506 3.9874 3.0108
slope -1.0240 -1.7255 -2.0732
Slope recovery (true vs posterior):
Eq True Mean Diff CI_lo CI_hi Covered
-----------------------------------------------------------------
eq1 -1.000 -1.0240 -0.0240 -1.1631 -0.8876 YES
eq2 -1.500 -1.7255 -0.2255 -1.8773 -1.5767 NO
eq3 -2.000 -2.0732 -0.0732 -2.2114 -1.9328 YES
Omega recovery (upper triangle):
Element True PostMean Diff
---------------------------------------------
Omega[1,1] 1.000 0.9681 -0.0319
Omega[1,2] 0.600 0.5047 -0.0953
Omega[1,3] 0.300 0.3260 +0.0260
Omega[2,2] 1.200 1.1247 -0.0753
Omega[2,3] 0.400 0.4217 +0.0217
Omega[3,3] 0.800 0.9805 +0.1805
Saved: sur_synth_bayesm_B.csv sur_synth_bayesm_Sigma.csv
Synthetic results (validated)¶
Posterior mean B:
eq1 eq2 eq3
intercept 4.9506 3.9874 3.0108
slope -1.0240 -1.7255 -2.0732
Slope recovery:
Eq True Mean Diff CI_lo CI_hi Covered
eq1 -1.000 -1.0240 -0.0240 -1.1631 -0.8876 YES
eq2 -1.500 -1.7255 -0.2255 -1.8773 -1.5767 NO
eq3 -2.000 -2.0732 -0.0732 -2.2114 -1.9328 YES
Omega recovery (upper triangle):
Omega[1,1]: true= 1.000 post= 0.9681 diff=-0.0319
Omega[1,2]: true= 0.600 post= 0.5047 diff=-0.0953
Omega[1,3]: true= 0.300 post= 0.3260 diff=+0.0260
Omega[2,2]: true= 1.200 post= 1.1247 diff=-0.0753
Omega[2,3]: true= 0.400 post= 0.4217 diff=+0.0217
Omega[3,3]: true= 0.800 post= 0.9805 diff=+0.1805
Coverage 2/3 for slopes (eq2 just misses: CI=[-1.877, -1.577] vs true=-1.5). P(≥1 miss out of 3) ≈ 14% under the null — expected finite-sample noise, not a bias. Omega diagonals within ±0.18, off-diagonals within ±0.10.
4. Tuna Data — Restricted Own-Brand Model¶
Source: Dominick's Finer Foods weekly scanner data, aggregated to chain level. Chevalier, J., Kashyap, A. & Rossi, P.E. (2003). Why don't prices rise during periods of peak demand? American Economic Review, 93(1), 15–37. Referenced in Rossi, Allenby & McCulloch (2005), Ch. 7.
Data: N=338 weeks × M=7 brands of canned tuna. Variables per brand:
MOVE: unit salesLPRICE: log retail priceNSALE: display/feature activity indicator (0–1)LWHPRIC: log wholesale price (not used here)
Brands:
| # | Label | Full name |
|---|---|---|
| 1 | StarKist | Star Kist 6 oz. |
| 2 | ChickSea | Chicken of the Sea 6 oz. |
| 3 | BBeeSolid | Bumble Bee Solid 6.12 oz. |
| 4 | BBeeChunk | Bumble Bee Chunk 6.12 oz. |
| 5 | Geisha | Geisha 6 oz. |
| 6 | BBeeLarge | Bumble Bee Large Cans |
| 7 | HHChunk | HH Chunk Lite 6.5 oz. |
SUR specification: restricted own-brand model with k=3 regressors per equation (intercept, own log-price, own feature). Each brand's log-sales is modelled as a separate equation; cross-brand contemporaneous correlation is captured by Ω.
In [5]:
Y_tuna <- as.matrix(read.csv("tuna_lmove.csv")) # 338 x 7
X_full <- as.matrix(read.csv("tuna_X.csv")) # 338 x 15
M_t <- 7L; k_t <- 3L; N_t <- nrow(Y_tuna); K_t <- M_t * k_t
# Brand names from bayesm tuna documentation (Dominick's Finer Foods)
BRANDS <- c("StarKist", "ChickSea", "BBeeSolid", "BBeeChunk",
"Geisha", "BBeeLarge", "HHChunk")
# Build per-brand design matrices: intercept, own-price, own-feature
regdata_tuna <- lapply(seq_len(M_t), function(j)
list(y = Y_tuna[, j],
X = cbind(X_full[, 1L], # intercept
X_full[, j + 1L], # LPRICE_j
X_full[, j + 8L]))) # NSALE_j
cat(sprintf("Tuna SUR: M=%d k=%d N=%d\n", M_t, k_t, N_t))
cat("\nLog-sales summary (mean and SD per brand):\n")
summ <- apply(Y_tuna, 2L, function(x) c(mean = mean(x), sd = sd(x)))
colnames(summ) <- BRANDS
print(round(summ, 3))
Tuna SUR: M=7 k=3 N=338
Log-sales summary (mean and SD per brand):
StarKist ChickSea BBeeSolid BBeeChunk Geisha BBeeLarge HHChunk
mean 9.521 9.000 7.713 8.952 7.908 6.864 8.757
sd 0.746 0.839 0.830 0.871 0.345 0.562 0.689
5. rsurGibbs — Tuna¶
In [6]:
Prior_t <- list(betabar = rep(0, K_t),
A = 0.01 * diag(K_t),
nu = M_t + 3L)
# V omitted: bayesm default V = nu * diag(M) is diffuse
Mcmc_t <- list(R = 25000L, keep = 1L, nprint = 5000L)
set.seed(42L)
cat("Running rsurGibbs on tuna data ...\n")
t0 <- proc.time()["elapsed"]
out_t <- rsurGibbs(Data = list(regdata = regdata_tuna),
Prior = Prior_t,
Mcmc = Mcmc_t)
cat(sprintf("Done in %.1f s\n", proc.time()["elapsed"] - t0))
Running rsurGibbs on tuna data ...
Starting Gibbs Sampler for SUR Regression Model
with 7 regressions
and 338 observations for each regression
Prior Parms:
betabar
[1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
A
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
[1,] 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[2,] 0.00 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[3,] 0.00 0.00 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[4,] 0.00 0.00 0.00 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[5,] 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[6,] 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[7,] 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.00 0.00 0.00 0.00
[8,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.00 0.00 0.00
[9,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.00 0.00
[10,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.00
[11,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.00
[12,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.00
[13,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01
[14,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[15,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[16,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[17,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[18,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[19,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[20,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[21,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21]
[1,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[2,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[3,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[4,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[5,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[6,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[7,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[8,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[9,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[10,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[11,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[12,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[13,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[14,] 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[15,] 0.00 0.01 0.00 0.00 0.00 0.00 0.00 0.00
[16,] 0.00 0.00 0.01 0.00 0.00 0.00 0.00 0.00
[17,] 0.00 0.00 0.00 0.01 0.00 0.00 0.00 0.00
[18,] 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.00
[19,] 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.00
[20,] 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.00
[21,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01
nu = 10
V =
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 10 0 0 0 0 0 0
[2,] 0 10 0 0 0 0 0
[3,] 0 0 10 0 0 0 0
[4,] 0 0 0 10 0 0 0
[5,] 0 0 0 0 10 0 0
[6,] 0 0 0 0 0 10 0
[7,] 0 0 0 0 0 0 10
MCMC parms:
R= 25000 keep= 1 nprint= 5000
MCMC Iteration (est time to end - min)
5000 (0.0)
10000 (0.0)
15000 (0.0)
20000 (0.0)
25000 (0.0)
Total Time Elapsed: 0.02
Done in 0.6 s
In [7]:
BURN_T <- 5000L
keep_t <- seq(BURN_T + 1L, 25000L)
bdraw_t <- out_t$betadraw[keep_t, ] # (20000, 21)
odraw_t <- out_t$Sigmadraw[keep_t, ] # (20000, 49)
B_pm_t <- matrix(colMeans(bdraw_t), k_t, M_t) # 3 x 7
Omega_pm_t <- matrix(colMeans(odraw_t), M_t, M_t) # 7 x 7
rownames(B_pm_t) <- c("intercept", "own_price", "own_feature")
colnames(B_pm_t) <- BRANDS
cat("Own-price elasticities (posterior mean and 95% CrI):\n")
cat(sprintf(" %-10s %8s %7s %8s %8s\n",
"Brand", "Mean", "Std", "CI_lo", "CI_hi"))
cat(" ", paste(rep("-", 50), collapse = ""), "\n", sep = "")
for (j in seq_len(M_t)) {
col_j <- (j - 1L) * k_t + 2L # own_price is 2nd param in each 3-block
d <- bdraw_t[, col_j]
cat(sprintf(" %-10s %8.4f %7.4f %8.4f %8.4f\n",
BRANDS[j], mean(d), sd(d),
quantile(d, 0.025), quantile(d, 0.975)))
}
cat("\nPosterior mean Omega (7x7):\n")
print(round(Omega_pm_t, 4))
# Save for Python comparison
write.csv(as.data.frame(t(B_pm_t)), "tuna_sur_bayesm_B.csv", row.names = FALSE)
write.csv(as.data.frame(Omega_pm_t), "tuna_sur_bayesm_Sigma.csv", row.names = FALSE)
cat("\nSaved: tuna_sur_bayesm_B.csv tuna_sur_bayesm_Sigma.csv\n")
Own-price elasticities (posterior mean and 95% CrI):
Brand Mean Std CI_lo CI_hi
--------------------------------------------------
StarKist -3.6654 0.3062 -4.2657 -3.0684
ChickSea -4.2708 0.2926 -4.8410 -3.6933
BBeeSolid -2.9077 1.0418 -4.9651 -0.8661
BBeeChunk -4.5927 0.2593 -5.1030 -4.0780
Geisha -4.7180 0.4675 -5.6311 -3.7995
BBeeLarge 0.6087 1.0808 -1.5047 2.7263
HHChunk -2.2899 0.4019 -3.0838 -1.5011
Posterior mean Omega (7x7):
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 0.3081 0.0342 -0.0640 0.0357 0.0178 -0.0588 0.0045
[2,] 0.0342 0.3411 -0.0106 0.0700 0.0316 -0.0213 0.0406
[3,] -0.0640 -0.0106 0.6348 0.0305 -0.0454 0.2781 -0.0238
[4,] 0.0357 0.0700 0.0305 0.3367 0.0183 -0.0238 0.0351
[5,] 0.0178 0.0316 -0.0454 0.0183 0.0846 -0.0257 0.0125
[6,] -0.0588 -0.0213 0.2781 -0.0238 -0.0257 0.3030 -0.0402
[7,] 0.0045 0.0406 -0.0238 0.0351 0.0125 -0.0402 0.3860
Saved: tuna_sur_bayesm_B.csv tuna_sur_bayesm_Sigma.csv
Tuna results (validated)¶
Own-price elasticities (posterior mean and 95% CrI):
Brand Mean Std CI_lo CI_hi
--------------------------------------------------
StarKist -3.6654 0.3062 -4.2657 -3.0684
ChickSea -4.2708 0.2926 -4.8410 -3.6933
BBeeSolid -2.9077 1.0418 -4.9651 -0.8661
BBeeChunk -4.5927 0.2593 -5.1030 -4.0780
Geisha -4.7180 0.4675 -5.6311 -3.7995
BBeeLarge 0.6087 1.0808 -1.5047 2.7263
HHChunk -2.2899 0.4019 -3.0838 -1.5011
Posterior mean Omega (7x7):
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 0.3081 0.0342 -0.0640 0.0357 0.0178 -0.0588 0.0045
[2,] 0.0342 0.3411 -0.0106 0.0700 0.0316 -0.0213 0.0406
[3,] -0.0640 -0.0106 0.6348 0.0305 -0.0454 0.2781 -0.0238
[4,] 0.0357 0.0700 0.0305 0.3367 0.0183 -0.0238 0.0351
[5,] 0.0178 0.0316 -0.0454 0.0183 0.0846 -0.0257 0.0125
[6,] -0.0588 -0.0213 0.2781 -0.0238 -0.0257 0.3030 -0.0402
[7,] 0.0045 0.0406 -0.0238 0.0351 0.0125 -0.0402 0.3860
Elasticity notes:
- Six of seven brands show clearly negative own-price elasticities with CIs well below zero.
- BBeeLarge (brand 6): point estimate +0.61, CI [−1.50, +2.73] straddles zero — wide uncertainty, likely reflecting sparse or noisy price variation for the large-can SKU.
- BBeeSolid (brand 3): SD=1.04, CI [−4.97, −0.87] — negative overall but imprecisely estimated.
- Most elastic: Geisha (−4.72) and BBeeChunk (−4.59); least elastic (among clearly identified): HHChunk (−2.29) and StarKist (−3.67).
Omega notes:
- Largest off-diagonal: Ω[3,6] = 0.278 (BBeeSolid ↔ BBeeLarge) — the two Bumble Bee sizes share correlated demand shocks, reflecting common promotions or shelf placement.
- Geisha has the smallest variance (Ω[5,5] = 0.085), indicating highly predictable log-sales.
- BBeeSolid has the largest variance (Ω[3,3] = 0.635), consistent with its wide elasticity CI.
6. Compare with Python¶
Loads CSVs saved by sur_python.ipynb. Agreement within Monte Carlo noise confirms the Python sur_gibbs.py module is a faithful replication of rsurGibbs.
In [8]:
py_synth_ok <- file.exists("sur_synth_python_B.csv") && file.exists("sur_synth_python_Sigma.csv")
py_tuna_ok <- file.exists("sur_tuna_python_B.csv") && file.exists("sur_tuna_python_Sigma.csv")
if (!py_synth_ok && !py_tuna_ok) {
cat("Python CSVs not found -- run sur_python.ipynb first.\n")
}
if (py_synth_ok) {
py_B_s <- as.matrix(read.csv("sur_synth_python_B.csv"))
py_Omega_s <- as.matrix(read.csv("sur_synth_python_Sigma.csv"))
cat("=== SYNTHETIC (M=3, k=2, N=200) ===\n")
cat(sprintf(" %-8s %7s %8s %8s %10s\n",
"Eq", "True", "R", "Python", "R-Python"))
cat(" ", paste(rep("-", 50), collapse = ""), "\n", sep = "")
for (j in seq_len(M_s)) {
tr <- B_TRUE[2, j]
r_v <- B_pm_s[2, j]
py_v <- py_B_s[j, 2]
cat(sprintf(" eq%-5d %7.3f %8.4f %8.4f %+10.5f\n",
j, tr, r_v, py_v, r_v - py_v))
}
cat(sprintf(" max|R-Python|: slope=%.5f Omega=%.5f\n",
max(abs(B_pm_s[2, ] - py_B_s[, 2])),
max(abs(Omega_pm_s - py_Omega_s))))
}
if (py_tuna_ok) {
py_B_t <- as.matrix(read.csv("sur_tuna_python_B.csv"))
py_Omega_t <- as.matrix(read.csv("sur_tuna_python_Sigma.csv"))
cat("\n=== TUNA (M=7, k=3, N=338) ===\n")
cat(sprintf(" %-12s %8s %8s %10s\n",
"Brand", "R", "Python", "R-Python"))
cat(" ", paste(rep("-", 45), collapse = ""), "\n", sep = "")
for (j in seq_len(M_t)) {
r_v <- B_pm_t[2, j]
py_v <- py_B_t[j, 2]
cat(sprintf(" %-12s %8.4f %8.4f %+10.5f\n",
BRANDS[j], r_v, py_v, r_v - py_v))
}
cat(sprintf(" max|R-Python|: price=%.5f Omega=%.5f\n",
max(abs(B_pm_t[2, ] - py_B_t[, 2])),
max(abs(Omega_pm_t - py_Omega_t))))
}
=== SYNTHETIC (M=3, k=2, N=200) === Eq True R Python R-Python -------------------------------------------------- eq1 -1.000 -1.0240 -1.0234 -0.00062 eq2 -1.500 -1.7255 -1.7240 -0.00150 eq3 -2.000 -2.0732 -2.0720 -0.00128 max|R-Python|: slope=0.00150 Omega=0.00630 === TUNA (M=7, k=3, N=338) === Brand R Python R-Python --------------------------------------------- StarKist -3.6654 -3.6614 -0.00395 ChickSea -4.2708 -4.2666 -0.00417 BBeeSolid -2.9077 -2.9142 +0.00647 BBeeChunk -4.5927 -4.5947 +0.00205 Geisha -4.7180 -4.7249 +0.00688 BBeeLarge 0.6087 0.5978 +0.01088 HHChunk -2.2899 -2.2908 +0.00093 max|R-Python|: price=0.01088 Omega=0.00303