Bayesian Hierarchical SUR

Python  ·  Download Gibbs sampler

Model

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