Bayesian Student-t Linear Model

Model

A handful of years in any long US macroeconomic series — the Depression, the war years — simply do not behave like the rest. A procedure that treats every year as equally informative has no way of saying so: it splits the difference, and a few extraordinary years quietly bend the answer for all the ordinary ones. The question here is how to let the model notice for itself which observations are exceptional, and discount each one in proportion to how exceptional it is.

The standard linear model assumes Gaussian errors. When data contain occasional large shocks — outliers from economic crises, measurement error, fat-tailed processes — the Normal assumption overstates the influence of extreme observations on $\beta$ . Replacing $\mathcal{N}(0, \sigma^2)$ errors with a Student- $t$ distribution provides automatic robustness: the heavier tails absorb outliers without requiring any manual flagging.

y_i = x_i'\beta + \varepsilon_i, \qquad \varepsilon_i \stackrel{\text{i.i.d.}}{\sim} t(\nu,\, 0,\, \sigma^2)

Scale-Mixture Representation

The key to tractable Bayesian inference is the scale-mixture-of-normals representation of the Student- $t$ . Introducing a latent precision weight $\lambda_i$ per observation makes all four full conditionals conjugate. Outlier observations receive low $\lambda_i$ (high local variance $\sigma^2/\lambda_i$ ), limiting their leverage on $\beta$ automatically.

y_i \mid \lambda_i \sim N\!\left(x_i'\beta,\; \frac{\sigma^2}{\lambda_i}\right), \qquad \lambda_i \mid \nu \sim \text{Gamma}\!\left(\tfrac{\nu}{2},\, \tfrac{\nu}{2}\right)

$\beta$ — regression coefficients; flat or $N(b_0, B_0)$ prior
$\sigma^2$ — error scale; $\text{IG}(a_0/2, s_0/2)$ prior (flat if $a_0 = s_0 = 0$ )
$\lambda_i$ — latent precision weight; low values flag outliers automatically
$\nu$ — degrees of freedom; discrete uniform on $\{2, \ldots, 40\}$

The degrees-of-freedom parameter $\nu$ controls tail heaviness: $\nu \to \infty$ recovers the Normal; $\nu = 4$ has substantially heavier tails. Rather than reparametrising $\nu$ for continuous HMC, Geweke places a discrete uniform prior on a grid $\{2, 3, \ldots, 40\}$ and draws $\nu$ as a categorical variable each iteration — fast and numerically stable near the $\nu = 2$ boundary.

4-Block Gibbs — Geweke (1993)

Block	Draw	Distribution
(1)	$\beta \mid \lambda, \sigma^2, y$	Normal conjugate — WLS posterior $N(\bar{b}, \bar{B})$ with $\bar{B}^{-1} = X'\Lambda X / \sigma^2$
(2)	$\sigma^2 \mid \beta, \lambda, y$	Inverse-Gamma conjugate $\text{IG}(a_n/2,\; s_n/2)$ where $s_n = s_0 + \textstyle\sum \lambda_i e_i^2$
(3)	$\lambda_i \mid \beta, \sigma^2, \nu, y$	Gamma conjugate $\text{Gamma}\!\left(\tfrac{\nu+1}{2},\; \tfrac{\nu + e_i^2/\sigma^2}{2}\right)$ independently
(4)	$\nu \mid \lambda$	Discrete draw on $\nu$ -grid via Gamma log-likelihood

Notebooks

Applied to two datasets. Section 1: synthetic contaminated data ( $n=120$ , $\nu^\star=4$ ), where the sampler recovers the true degrees of freedom and correctly identifies downweighted outliers. Section 2: Nelson-Plosser (1982) 14 annual U.S. macroeconomic series (1909–1970), where low posterior $\nu$ for series like Real GNP, Industrial Production, and Employment corresponds to Depression- and WWII-era shocks — observations the model down-weights rather than attributes to a stochastic trend. Section 3: PyMC comparison with continuous $\nu$ via NUTS on both datasets. Two companion notebooks cross-check the model against general-purpose tools. The PyMC notebook treats $\nu$ as continuous with a weakly-informative Gamma prior and samples $(\beta, \sigma, \nu)$ jointly by NUTS — the $\beta$ posteriors match the discrete-ν Gibbs to within a posterior SD, with the trade-off that NUTS marginalises out the latent $\lambda_i$ (so it loses the per-observation outlier weights the scale-mixture Gibbs exposes directly). The R notebook adds the frequentist counterparts — hett::tlm (Student-t by maximum likelihood, recovering $\nu\approx 3.9$ on the synthetic $t(4)$ data), MASS::rlm (Huber M-estimation), and lm (the non-robust OLS baseline that gets dragged toward the outliers). Gibbs, NUTS, and ML all agree the macro series are heavy-tailed.

Python Student-t Gibbs sampler (discrete ν), Nelson-Plosser application (Python) Open → Python PyMC cross-check — continuous ν via pm.StudentT / NUTS, vs the discrete-ν Gibbs (Python) Open → R R cross-check — hett::tlm (Student-t ML), MASS::rlm (robust), lm (OLS baseline) (R) Open →

Downloads

student_t_lm_gibbs.py 4-block Gibbs sampler module nelplo.csv Nelson-Plosser (1982) macroeconomic series — 14 annual U.S. series, 1909–1970

Gibbs Sampler — Source Code

"""
student_t_lm_gibbs.py

Bayesian inference for the linear model with i.i.d. Student-t errors.
Implements Geweke (1993) "Bayesian treatment of the independent Student-t
linear model."

Model
-----
    y_i = x_i' β + ε_i,    ε_i ~ i.i.d. t(ν, 0, σ²)

Scale-mixture representation (Student-t = Normal scale-mixture)
---------------------------------------------------------------
    y_i | λ_i  ~  N(x_i' β,  σ²/λ_i)
    λ_i | ν    ~  Gamma(ν/2, ν/2)          [shape / rate]
    ε_i  marginally  ~  t(ν, 0, σ²)

Four-block Gibbs sampler
------------------------
    (1)  β    | λ, σ², y  ~  N(b̄, B̄)          WLS posterior
    (2)  σ²   | β, λ, y   ~  IG(aₙ/2, sₙ/2)
    (3)  λᵢ   | β, σ², ν, y  ~  Gamma((ν+1)/2,  (ν + eᵢ²/σ²)/2)
    (4)  ν    | λ          ~  discrete on ν_grid  (via Gamma log-likelihood)

Default priors
--------------
    β   ~  flat  (or N(b0, B0) if supplied)
    σ²  ~  IG(a0/2, s0/2)       flat when a0=s0=0
    ν   ~  discrete uniform on {2, 3, …, 40}

Reference
---------
    Geweke, J. (1993). "Bayesian treatment of the independent Student-t
    linear model." Journal of Applied Econometrics, 8(S1), S19-S40.
"""

import time
import numpy as np
from scipy.special import gammaln


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

def student_t_gibbs(
    y, X,
    nu_grid=None,
    b0=None, B0=None,
    a0=0.0, s0=0.0,
    R=10000, burn=2000,
    seed=42, nprint=1000,
):
    """
    Four-block Gibbs sampler for the Student-t linear model.

    Parameters
    ----------
    y       : (n,)        observations
    X       : (n, k)      design matrix (include intercept column if desired)
    nu_grid : (G,) array  grid of candidate ν values  [default: 2, 3, …, 40]
    b0      : (k,)        prior mean for β             [flat if omitted]
    B0      : (k, k)      prior covariance for β       [flat if omitted]
    a0, s0  : float       IG(a0/2, s0/2) prior for σ²  [flat if both 0]
    R, burn : int         total and burn-in iterations
    seed    : int
    nprint  : int         print progress every nprint iterations (0 = silent)

    Returns
    -------
    dict with keys
      'beta'   : (R_kept, k)  posterior draws for β
      'sigma2' : (R_kept,)    posterior draws for σ²
      'lam'    : (R_kept, n)  posterior draws for the mixing weights λᵢ
      'nu'     : (R_kept,)    posterior draws for ν  (discrete)
    """
    rng   = np.random.default_rng(seed)
    y     = np.asarray(y, float).ravel()
    X     = np.asarray(X, float)
    n, k  = X.shape

    if nu_grid is None:
        nu_grid = np.arange(2.0, 41.0)
    nu_grid = np.asarray(nu_grid, float)
    G = len(nu_grid)

    # ── prior ─────────────────────────────────────────────────────────────────
    if b0 is not None and B0 is not None:
        B0_inv = np.linalg.inv(np.asarray(B0, float))
        b0     = np.asarray(b0, float)
    else:
        B0_inv = None

    # ── initial values ────────────────────────────────────────────────────────
    beta   = np.linalg.lstsq(X, y, rcond=None)[0]
    e      = y - X @ beta
    sigma2 = max(float(e @ e) / max(n - k, 1), 1e-8)
    lam    = np.ones(n)
    nu     = float(nu_grid[G // 2])

    # ── pre-compute ν-dependent log-normaliser (constant across iterations) ──
    # log p(λ | ν) = n × [(ν/2) log(ν/2) − logΓ(ν/2)]
    #              + (ν/2 − 1) Σ log λᵢ  −  ν/2 Σ λᵢ
    half_nu    = nu_grid / 2.0
    log_norm   = n * (half_nu * np.log(half_nu) - gammaln(half_nu))  # (G,)

    # ── storage ───────────────────────────────────────────────────────────────
    R_kept       = R - burn
    beta_draws   = np.empty((R_kept, k))
    sigma2_draws = np.empty(R_kept)
    lam_draws    = np.empty((R_kept, n))
    nu_draws     = np.empty(R_kept)
    kept = 0
    t0   = time.time()

    for r in range(1, R + 1):

        # ── Block 1: β | λ, σ², y  (WLS posterior) ───────────────────────────
        XtLX = X.T @ (lam[:, None] * X)   # (k, k)  X'ΛX
        XtLy = X.T @ (lam * y)            # (k,)    X'Λy

        if B0_inv is not None:
            Prec = XtLX / sigma2 + B0_inv
            rhs  = XtLy / sigma2 + B0_inv @ b0
        else:
            Prec = XtLX / sigma2
            rhs  = XtLy / sigma2

        Bbar   = np.linalg.inv(Prec)
        bbar   = Bbar @ rhs
        L_chol = np.linalg.cholesky(Bbar)
        beta   = bbar + L_chol @ rng.standard_normal(k)

        # ── Block 2: σ² | β, λ, y  (IG conjugate) ───────────────────────────
        e      = y - X @ beta
        SS_w   = float(lam @ (e * e))     # Σ λᵢ eᵢ²
        a_n    = a0 + n
        s_n    = s0 + SS_w
        sigma2 = s_n / rng.chisquare(a_n)

        # ── Block 3: λᵢ | β, σ², ν, y  ~ Gamma((ν+1)/2, (ν+eᵢ²/σ²)/2) ─────
        shape_lam = (nu + 1.0) / 2.0
        rate_lam  = (nu + e * e / sigma2) / 2.0   # (n,)
        lam       = rng.gamma(shape_lam, 1.0 / rate_lam)

        # ── Block 4: ν | λ  (discrete draw) ──────────────────────────────────
        sum_log_lam = float(np.log(lam).sum())
        sum_lam     = float(lam.sum())
        log_p = log_norm + (half_nu - 1.0) * sum_log_lam - half_nu * sum_lam
        log_p -= log_p.max()       # shift for numerical stability
        p = np.exp(log_p)
        p /= p.sum()
        nu = float(rng.choice(nu_grid, p=p))

        # ── store ─────────────────────────────────────────────────────────────
        if r > burn:
            beta_draws[kept]   = beta
            sigma2_draws[kept] = sigma2
            lam_draws[kept]    = lam
            nu_draws[kept]     = nu
            kept += 1

        if nprint > 0 and r % nprint == 0:
            elapsed = time.time() - t0
            print(f'  Iter {r:5d}/{R}  ({elapsed:.1f}s elapsed)')

    elapsed = time.time() - t0
    print(f'Done in {elapsed / 60:.1f} min  |  kept draws: {R_kept}')
    return {
        'beta':   beta_draws,
        'sigma2': sigma2_draws,
        'lam':    lam_draws,
        'nu':     nu_draws,
    }


def normal_gibbs(y, X, b0=None, B0=None, a0=0.0, s0=0.0,
                 R=10000, burn=2000, seed=42, nprint=0):
    """
    Gibbs sampler for the Normal linear model (ν → ∞ limit of Student-t).

    Identical to student_t_gibbs but with λᵢ ≡ 1 (no mixing weights).
    Used as a baseline for model comparison.

    Returns
    -------
    dict with keys 'beta' and 'sigma2'.
    """
    rng   = np.random.default_rng(seed)
    y     = np.asarray(y, float).ravel()
    X     = np.asarray(X, float)
    n, k  = X.shape

    if b0 is not None and B0 is not None:
        B0_inv = np.linalg.inv(np.asarray(B0, float))
        b0     = np.asarray(b0, float)
    else:
        B0_inv = None

    beta   = np.linalg.lstsq(X, y, rcond=None)[0]
    e      = y - X @ beta
    sigma2 = max(float(e @ e) / max(n - k, 1), 1e-8)

    R_kept       = R - burn
    beta_draws   = np.empty((R_kept, k))
    sigma2_draws = np.empty(R_kept)
    kept = 0
    t0   = time.time()

    XtX = X.T @ X
    Xty = X.T @ y

    for r in range(1, R + 1):
        if B0_inv is not None:
            Prec = XtX / sigma2 + B0_inv
            rhs  = Xty / sigma2 + B0_inv @ b0
        else:
            Prec = XtX / sigma2
            rhs  = Xty / sigma2

        Bbar   = np.linalg.inv(Prec)
        bbar   = Bbar @ rhs
        L_chol = np.linalg.cholesky(Bbar)
        beta   = bbar + L_chol @ rng.standard_normal(k)

        e      = y - X @ beta
        a_n    = a0 + n
        s_n    = s0 + float(e @ e)
        sigma2 = s_n / rng.chisquare(a_n)

        if r > burn:
            beta_draws[kept]   = beta
            sigma2_draws[kept] = sigma2
            kept += 1

        if nprint > 0 and r % nprint == 0:
            print(f'  Iter {r:5d}/{R}  ({time.time()-t0:.1f}s elapsed)')

    if nprint > 0:
        print(f'Done in {(time.time()-t0)/60:.1f} min  |  kept draws: {R_kept}')
    return {'beta': beta_draws, 'sigma2': sigma2_draws}


def posterior_summary(draws, param_names=None, true_vals=None):
    """
    Posterior mean, SD, and 95% credible interval for each column.

    Parameters
    ----------
    draws       : (R, p) array of posterior draws
    param_names : list of p strings       [optional]
    true_vals   : (p,) true values; adds 'true' and 'covered' columns

    Returns
    -------
    pandas.DataFrame
    """
    import pandas as pd
    p = draws.shape[1]
    if param_names is None:
        param_names = [f'b{j}' for j in range(p)]
    lo = np.quantile(draws, 0.025, axis=0)
    hi = np.quantile(draws, 0.975, axis=0)
    df = pd.DataFrame({
        'mean': draws.mean(0),
        'sd':   draws.std(0),
        'q025': lo,
        'q975': hi,
    }, index=param_names)
    if true_vals is not None:
        tv = np.asarray(true_vals)
        df['true']    = tv
        df['covered'] = (tv >= lo) & (tv <= hi)
    return df.round(4)


def log_lik_t(y, X, beta_draws, sigma2_draws, nu_draws):
    """
    Posterior mean log-likelihood under the Student-t model.

    Evaluates log p(y | βᵣ, σ²ᵣ, νᵣ) for each draw r and returns the mean.
    Used for informal model comparison (not a proper Bayes factor).
    """
    from scipy.stats import t as t_dist
    R = len(nu_draws)
    lls = np.empty(R)
    for r in range(R):
        e = y - X @ beta_draws[r]
        lls[r] = t_dist.logpdf(
            e, df=nu_draws[r], loc=0, scale=np.sqrt(sigma2_draws[r])
        ).sum()
    return float(lls.mean()), float(lls.std())


def log_lik_normal(y, X, beta_draws, sigma2_draws):
    """Posterior mean log-likelihood under the Normal model."""
    from scipy.stats import norm
    R = len(sigma2_draws)
    lls = np.empty(R)
    for r in range(R):
        e = y - X @ beta_draws[r]
        lls[r] = norm.logpdf(e, loc=0, scale=np.sqrt(sigma2_draws[r])).sum()
    return float(lls.mean()), float(lls.std())

References

Geweke, J. (1993). Bayesian treatment of the independent Student-t linear model. Journal of Applied Econometrics 8(S1), S19–S40. — the model this example implements
Nelson, C. R. & Plosser, C. I. (1982). Trends and random walks in macroeconomic time series: some evidence and implications. Journal of Monetary Economics 10(2), 139–162. — the 14 annual U.S. macroeconomic series used in the application
Perron, P. (1989). The great crash, the oil price shock, and the unit root hypothesis. Econometrica 57(6), 1361–1401. — the structural-break re-reading of the same Nelson–Plosser series