Bayesian Hierarchical SUR

Python  ·  Download Gibbs sampler

Model

A chain selling 11 brands of orange juice across 83 stores faces two questions at once, and neither can be answered properly without the other. Within a store the 11 brands compete, so whatever is missing from the model — a display at the end of the aisle, a hot week — hits all 11 at the same moment. Across stores, the same brand behaves differently in different neighbourhoods: a commodity filler in one, the preferred regional brand in another. This page is what it takes to hold both facts at once.

The hSUR model stacks two sources of structure. Across equations within a store: SUR cross-equation error correlation Ω\boldsymbol{\Omega}. Across stores within an equation: a hierarchical prior that shrinks store-level coefficients β(i)\boldsymbol{\beta}^{(i)} toward a population mean Δ⊤zi\Delta^\top z_i.

FeatureSURHLMhSUR
Units1 market, M eqsN stores, 1 eqN stores, M eqs
Cross-unit link—βi∼N(Δ⊤zi,Vβ)\beta_i \sim \mathcal{N}(\Delta^\top z_i, V_\beta)Both
Cross-eq linkΩ\boldsymbol{\Omega}—Ω\boldsymbol{\Omega}
Gibbs blocks23–44

Neither SUR nor HLM alone captures both dimensions. hSUR requires all four Gibbs blocks because no closed-form joint posterior exists.

Observation model

yijt=xijt⊤βj(i)+εijt,εt(i)∼NM(0, Ω)y_{ijt} = \mathbf{x}_{ijt}^\top \boldsymbol{\beta}_j^{(i)} + \varepsilon_{ijt}, \qquad \boldsymbol{\varepsilon}_t^{(i)} \sim \mathcal{N}_M(\mathbf{0},\, \boldsymbol{\Omega})

Hierarchy

β(i)∼N ⁣(Δ⊤zi,  Vβ),i=1,…,N\boldsymbol{\beta}^{(i)} \sim \mathcal{N}\!\left(\Delta^\top z_i,\; V_\beta\right), \qquad i = 1,\ldots,N

Priors

Ω∼IW(ν, V),Vβ∼IW(νβ, Vβ0)\boldsymbol{\Omega} \sim \mathcal{IW}(\nu,\, \mathbf{V}), \qquad V_\beta \sim \mathcal{IW}(\nu_\beta,\, V_\beta^0)

4-Block Gibbs

The 4-block Gibbs sampler cycles through exact full conditionals:

  1. β(i)∣Ω,Δ,Vβ\boldsymbol{\beta}^{(i)} \mid \boldsymbol{\Omega}, \Delta, V_\beta — Store-specific coefficient vectors. Each store's K×1K \times 1 coefficient vector has a Gaussian posterior shrunk toward the population mean Δ⊤zi\Delta^\top z_i with precision Vβ−1V_\beta^{-1} plus the within-store SUR precision.
  2. Ω∣{β(i)}\boldsymbol{\Omega} \mid \{\boldsymbol{\beta}^{(i)}\} — Cross-equation error covariance. Inverse-Wishart on the N×MN \times M residual matrix after removing store-level fits.
  3. Δ∣{β(i)},Vβ\Delta \mid \{\boldsymbol{\beta}^{(i)}\}, V_\beta — Population regression. Matrix-Normal update: how unit-level covariates ziz_i map to average coefficient levels.
  4. Vβ∣{β(i)},ΔV_\beta \mid \{\boldsymbol{\beta}^{(i)}\}, \Delta — Store-level covariance. Inverse-Wishart on deviations β(i)−Δ⊤zi\boldsymbol{\beta}^{(i)} - \Delta^\top z_i — controls how much stores differ from the population.

Notebooks

A single, long notebook that builds the 4-block sampler (hsur_gibbs.py) and takes it through five stages. Synthetic validation first (50 stores × 3 equations × 30 periods) — all four blocks recover the generating Δ\Delta and Ω\Omega inside their 95% intervals. Then the orange-juice panel: Dominick's Finer Foods weekly scanner data (83 stores × 11 brands × ≈117 weeks, TiT_i ranging 87–121), regressing log unit-sales on own log-price and a deal flag, with cross-brand correlated shocks. The three output blocks each answer a different question. Δ\Delta — population elasticities: all 11 credible intervals sit entirely below zero with SD ≈ 0.09–0.12, precision bought by pooling across 83 stores. A clear size-format pattern emerges — within each manufacturer the 96 oz jug is materially less elastic than the 64 oz (Tropicana −2.00 vs −2.80; Minute Maid −2.47 vs −3.18) — and a private-label paradox: Dominick's own 96 oz is the least elastic brand in the sample (−1.50), below even premium Tropicana. Ω\Omega — residual structure: the same 96 oz formats have residual variances under 0.10, roughly 2.5–7× below their 64 oz siblings — habitual bulk buyers are highly predictable once price and deal are known, while 64 oz purchases are impulsive and full of unmodelled promotion. β(i)\beta^{(i)} — store-level heterogeneity: shrinkage is unbiased (store means land on the population mean), but the spread is large and it reverses conclusions: FloGold ranges from −5.91 to −1.54 across stores — a commodity filler in some, the preferred regional brand in others — so a single chain-wide price set from the −4.06 average would badly misprice both ends. Two extensions close it out. Putting store demographics in ZZ (396 parameters from 83 stores) turns Δ\Delta into a regression of each coefficient on income, education and competition — and the honest finding is that it explains little: only 5 effects are distinguishable from zero, and the store-level spread barely moves. The unrestricted cross-price model (K=143K=143 per store, a ~90-minute run) then puts all 11 log-prices in every equation, and overturns the naive diagnosis: the truly confounded brand is Trop64b (−3.37 → −4.42, jumping from 4th to 2nd most elastic), while Dom96 is essentially unchanged (+0.010) and robustly the least elastic — a pooled-OLS scoping check had predicted the opposite, and the hierarchy resolves the paradox.

Python From-scratch 4-block Gibbs — synthetic validation, the orange-juice panel (population elasticities, Ω residual structure, store-level heterogeneity), store demographics in Z, and the unrestricted cross-price model (Python) Open →

Downloads

hsur_gibbs.py Gibbs sampler module oj_yx.csv Orange juice log-sales, prices, deal indicators oj_storedemo.csv Store-level demographic covariates

Gibbs Sampler — Source Code

"""
Bayesian Hierarchical SUR (hSUR) -- Gibbs sampler.

Model
-----
Store i (i = 1,...,nreg), M equations:
  y_{ij} = X_{ij} @ beta_j^(i) + eps_{ij},   j = 1,...,M
  eps_{it} ~ N_M(0, Omega)      cross-equation correlated errors per period

Each store has its own coefficient vector:
  beta^(i) = (beta_1^(i)', ..., beta_M^(i)')  of length K = sum(k_j)

Hierarchy
---------
  beta^(i) = Delta.T @ z_i + u_i,   u_i ~ N(0, Vbeta)

  z_i : (nz,)  store-level covariates  (include an intercept column)
  Delta: (nz, K)  population regression of beta on z

Priors
------
  vec(Delta) | Vbeta ~ N(vec(Deltabar), Vbeta (x) A^{-1})   (matrix normal)
  Vbeta  ~ IW(nu_beta, V_beta)
  Omega  ~ IW(nu,      V)

4-block Gibbs
-------------
  (1) beta^(i) | Omega, Delta, Vbeta, y_i, X_i  -- N  (one draw per store)
      Q_i = Vbeta^{-1} + [[omega^{pq} X_{ip}' X_{iq}]]
      b_i = Vbeta^{-1} (Delta.T z_i) + [[sum_q omega^{pq} X_{ip}' y_{iq}]]

  (2) Omega | {beta^(i)}, Y  -- IW(nu + N_total, V + sum_i E_i'E_i)
      N_total = sum_i T_i,  E_i[:,j] = y_{ij} - X_{ij} beta_j^(i)

  (3) Delta | {beta^(i)}, Vbeta  -- MN (same as rhierLinearModel)
      A_tilde = A + Z'Z,  Delta_tilde = A_tilde^{-1}(A Deltabar + Z' B)
      B = [beta^(1)', ..., beta^(nreg)']   (nreg x K)

  (4) Vbeta | {beta^(i)}, Delta  -- IW(nu_beta + nreg, V_beta + D'D)
      D_i = beta^(i) - Delta.T z_i
"""

import numpy as np
from scipy.stats import invwishart
from scipy.linalg import cholesky, solve_triangular


# ── private helpers ───────────────────────────────────────────────────────────

def _chol_draw(Q, b, rng):
    """Draw x ~ N(Q^{-1} b, Q^{-1}) via Cholesky of the precision matrix."""
    L  = cholesky(Q, lower=True)
    v  = solve_triangular(L,   b, lower=True)
    mu = solve_triangular(L.T, v, lower=False)
    return mu + solve_triangular(L.T, rng.standard_normal(len(b)), lower=False)


def _draw_Delta(betas, Z, Vbeta, Deltabar, A, rng):
    """Delta | Vbeta, betas ~ MN(Delta_tilde, Vbeta, A_tilde^{-1}).
    betas : (nreg, K);  Z : (nreg, nz);  Delta : (nz, K)
    """
    A_tilde     = A + Z.T @ Z
    A_tilde_inv = np.linalg.inv(A_tilde)
    Delta_tilde = A_tilde_inv @ (A @ Deltabar + Z.T @ betas)
    L_A = np.linalg.cholesky(A_tilde_inv)
    L_V = np.linalg.cholesky(Vbeta)
    return Delta_tilde + L_A @ rng.standard_normal(Delta_tilde.shape) @ L_V.T


def _draw_Vbeta(betas, Z, Delta, nu, V, rng):
    """Vbeta | betas, Delta ~ IW(nu + nreg, V + D'D)."""
    resid = betas - Z @ Delta         # (nreg, K)
    return invwishart.rvs(df=nu + len(betas),
                          scale=V + resid.T @ resid,
                          random_state=rng)


def _build_store_stats(y_list, X_list, M):
    """Precompute S[p,q] = X_p'X_q and XtY[p,q] = X_p'y_q for one store."""
    S   = [[X_list[p].T @ X_list[q] for q in range(M)] for p in range(M)]
    XtY = [[X_list[p].T @ y_list[q] for q in range(M)] for p in range(M)]
    return S, XtY


def _draw_beta_i(S_i, XtY_i, starts, Omega_inv, Vbeta_inv, mu_i, rng):
    """Draw beta^(i) ~ N(Q_i^{-1} b_i, Q_i^{-1}).
    S_i    : M x M list of (k_p, k_q) matrices  X_p' X_q
    XtY_i  : M x M list of (k_p,)    vectors    X_p' y_q
    mu_i   : (K,) prior mean  Delta.T @ z_i
    """
    K = starts[-1]
    M = Omega_inv.shape[0]
    Q = Vbeta_inv.copy()
    b = Vbeta_inv @ mu_i
    for p in range(M):
        sp = starts[p]; ep = starts[p + 1]
        for q in range(M):
            sq = starts[q]; eq = starts[q + 1]
            Q[sp:ep, sq:eq] += Omega_inv[p, q] * S_i[p][q]
            b[sp:ep]        += Omega_inv[p, q] * XtY_i[p][q]
    return _chol_draw(Q, b, rng)


# ── public API ────────────────────────────────────────────────────────────────

def rhsur_gibbs(regdata, Z=None, Prior=None, Mcmc=None, seed=42):
    """
    4-block Gibbs sampler for the Hierarchical SUR model.

    Parameters
    ----------
    regdata : list of nreg dicts, each with keys
                'y_list' : list of M arrays, each (T_i,)
                'X_list' : list of M arrays, each (T_i, k_j)
              Alternatively pass {'Y': (T_i, M), 'X': (T_i, k)} for the
              common-X case (same regressors across all equations in store i).
    Z       : (nreg, nz) store-level covariates including intercept.
              Default: column of ones  (intercept-only hierarchy).
    Prior   : dict with optional keys
                Deltabar : (nz, K)  prior mean of Delta,   default 0
                A        : (nz, nz) prior precision of Delta, default 0.01*I
                nu       : int      IW df for Omega,        default M+3
                V        : (M, M)   IW scale for Omega,     default (M+2)*I
                nu_beta  : int      IW df for Vbeta,        default K+3
                V_beta   : (K, K)   IW scale for Vbeta,     default (K+2)*I
    Mcmc    : dict with optional keys
                R      : int  total draws,          default 10000
                burn   : int  burn-in,              default 0.2*R
                keep   : int  thinning,             default 1
                nprint : int  print frequency,      default R//10

    Returns
    -------
    dict with
      betadraw  : (nreg, K, R_kept)  store-level coefficient draws
      Deltadraw : (R_kept, nz*K)     column-major vec(Delta)
      Vbetadraw : (R_kept, K*K)      column-major vec(Vbeta)
      Omegadraw : (R_kept, M*M)      column-major vec(Omega)
    """
    if Prior is None:
        Prior = {}
    if Mcmc is None:
        Mcmc = {'R': 10000}

    rng   = np.random.default_rng(seed)
    nreg  = len(regdata)

    # ── normalise regdata to y_list / X_list ─────────────────────────────────
    def _normalise(d):
        if 'y_list' in d:
            return d['y_list'], d['X_list']
        Y = d['Y']; X = d['X']
        M_ = Y.shape[1]
        return [Y[:, j] for j in range(M_)], [X] * M_

    y_data = [_normalise(d) for d in regdata]   # list of (y_list, X_list)

    M     = len(y_data[0][0])
    ks    = [y_data[0][1][j].shape[1] for j in range(M)]
    K     = sum(ks)
    starts = np.concatenate([[0], np.cumsum(ks)]).astype(int)

    # ── Z ────────────────────────────────────────────────────────────────────
    if Z is None:
        Z = np.ones((nreg, 1))
    Z = np.asarray(Z, float)
    nz = Z.shape[1]

    # ── priors ───────────────────────────────────────────────────────────────
    Deltabar = np.asarray(Prior.get('Deltabar', np.zeros((nz, K))),       float)
    A        = np.asarray(Prior.get('A',        0.01 * np.eye(nz)),       float)
    nu       = int(Prior.get('nu',      M + 3))
    V        = np.asarray(Prior.get('V',        float(M + 2) * np.eye(M)), float)
    nu_beta  = int(Prior.get('nu_beta', K + 3))
    V_beta   = np.asarray(Prior.get('V_beta',   float(K + 2) * np.eye(K)), float)

    # ── MCMC settings ────────────────────────────────────────────────────────
    R      = int(Mcmc['R'])
    burn   = int(Mcmc.get('burn',   0.2 * R))
    keep   = int(Mcmc.get('keep',   1))
    nprint = int(Mcmc.get('nprint', max(R // 10, 1)))
    R_kept = (R - burn) // keep

    # ── precompute sufficient stats per store ─────────────────────────────────
    Ts    = [len(y_data[i][0][0]) for i in range(nreg)]
    N_tot = sum(Ts)
    stats = [_build_store_stats(y_data[i][0], y_data[i][1], M)
             for i in range(nreg)]           # stats[i] = (S_i, XtY_i)

    # ── storage ───────────────────────────────────────────────────────────────
    betadraw  = np.zeros((nreg, K, R_kept))
    Deltadraw = np.zeros((R_kept, nz * K))
    Vbetadraw = np.zeros((R_kept, K * K))
    Omegadraw = np.zeros((R_kept, M * M))

    # ── initialise ────────────────────────────────────────────────────────────
    betas  = np.zeros((nreg, K))
    Delta  = Deltabar.copy()
    Vbeta  = V_beta  / max(nu_beta - K - 1, 1)
    Omega  = V       / max(nu      - M - 1, 1)

    i_kept = 0
    for r in range(R):
        Omega_inv  = np.linalg.inv(Omega)
        Vbeta_inv  = np.linalg.inv(Vbeta)
        mus        = Z @ Delta            # (nreg, K)  prior means

        # ── block 1: beta^(i) | Omega, Delta, Vbeta ──────────────────────────
        for i in range(nreg):
            S_i, XtY_i = stats[i]
            betas[i] = _draw_beta_i(S_i, XtY_i, starts,
                                    Omega_inv, Vbeta_inv, mus[i], rng)

        # ── block 2: Omega | {beta^(i)} ───────────────────────────────────────
        S_eps = np.zeros((M, M))
        for i in range(nreg):
            ys_i, Xs_i = y_data[i]
            E_i = np.column_stack(
                [ys_i[j] - Xs_i[j] @ betas[i, starts[j]:starts[j + 1]]
                 for j in range(M)])                         # (T_i, M)
            S_eps += E_i.T @ E_i
        Omega = invwishart.rvs(df=nu + N_tot, scale=V + S_eps,
                               random_state=rng)

        # ── block 3: Delta | {beta^(i)}, Vbeta ───────────────────────────────
        Delta = _draw_Delta(betas, Z, Vbeta, Deltabar, A, rng)

        # ── block 4: Vbeta | {beta^(i)}, Delta ───────────────────────────────
        Vbeta = _draw_Vbeta(betas, Z, Delta, nu_beta, V_beta, rng)

        if nprint and (r + 1) % nprint == 0:
            print(f'  Iter {r+1:6d}/{R}')

        if r >= burn and (r - burn) % keep == 0:
            betadraw[:, :, i_kept]  = betas
            Deltadraw[i_kept]       = Delta.ravel(order='F')
            Vbetadraw[i_kept]       = Vbeta.ravel(order='F')
            Omegadraw[i_kept]       = Omega.ravel(order='F')
            i_kept += 1

    return {'betadraw': betadraw, 'Deltadraw': Deltadraw,
            'Vbetadraw': Vbetadraw, 'Omegadraw': Omegadraw}


def simulate_hsur(nreg, M, k, T, true_Delta, true_Vbeta, true_Omega,
                  Z=None, seed=0):
    """
    Simulate data from the hierarchical SUR model.

    Parameters
    ----------
    nreg        : number of stores
    M           : number of equations
    k           : regressors per equation (int, same for all)
    T           : observations per store (int or list of ints)
    true_Delta  : (nz, K) population regression matrix
    true_Vbeta  : (K, K)  store-level coefficient covariance
    true_Omega  : (M, M)  cross-equation error covariance
    Z           : (nreg, nz) store covariates; default ones (intercept only)
    seed        : random seed

    Returns
    -------
    regdata : list of nreg dicts with 'y_list' and 'X_list'
    betas   : (nreg, K) true store-level coefficients
    Z       : (nreg, nz) store covariates used
    """
    rng = np.random.default_rng(seed)
    K   = M * k
    Ts  = [T] * nreg if np.isscalar(T) else list(T)

    if Z is None:
        Z = np.ones((nreg, 1))
    Z = np.asarray(Z, float)

    # draw store-level coefficients from hierarchy
    L_V  = np.linalg.cholesky(true_Vbeta)
    mus  = Z @ true_Delta                              # (nreg, K)
    betas = mus + rng.standard_normal((nreg, K)) @ L_V.T

    L_O  = np.linalg.cholesky(true_Omega)
    starts = np.arange(0, K + 1, k)

    regdata = []
    for i in range(nreg):
        Ti   = Ts[i]
        E_i  = rng.standard_normal((Ti, M)) @ L_O.T    # correlated errors
        y_list = []; X_list = []
        for j in range(M):
            X_ij = np.column_stack([np.ones(Ti),
                                    rng.standard_normal((Ti, k - 1))])
            y_ij = X_ij @ betas[i, starts[j]:starts[j + 1]] + E_i[:, j]
            y_list.append(y_ij)
            X_list.append(X_ij)
        regdata.append({'y_list': y_list, 'X_list': X_list})

    return regdata, betas, Z


def posterior_summary(draws, param_names=None, true_vals=None):
    """Posterior mean, SD, 95% CrI for each column of a (R, p) draws array."""
    import pandas as pd
    p = draws.shape[1]
    if param_names is None:
        param_names = [f'p{j}' for j in range(p)]
    rows = []
    for j, name in enumerate(param_names):
        d      = draws[:, j]
        lo, hi = np.percentile(d, [2.5, 97.5])
        row    = {'param': name, 'mean': d.mean(), 'sd': d.std(),
                  'q025': lo, 'q975': hi}
        if true_vals is not None:
            row['true']    = float(true_vals[j])
            row['covered'] = bool(lo <= true_vals[j] <= hi)
        rows.append(row)
    return pd.DataFrame(rows).set_index('param')

References