Bayesian Student-t Linear Model

Geweke (1993) — Python implementation

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

Key insight — scale-mixture of normals: $$y_i \mid \lambda_i \sim N(x_i'\beta,\; \sigma^2/\lambda_i), \qquad \lambda_i \mid \nu \sim \text{Gamma}(\nu/2,\; \nu/2)$$

Augmenting with the latent weights $\lambda_i$ makes all full conditionals conjugate → four-block Gibbs sampler.

Block	Update
1. $\beta \mid \lambda, \sigma^2, y$	$N(\bar{b}, \bar{B})$ — WLS posterior with precision matrix $\bar{B}^{-1} = X'\Lambda X/\sigma^2$
2. $\sigma^2 \mid \beta, \lambda, y$	$\text{IG}(a_n/2,\; s_n/2)$ where $s_n = s_0 + \sum \lambda_i e_i^2$
3. $\lambda_i \mid \beta, \sigma^2, \nu, y$	$\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\text{-grid}$ via Gamma log-likelihood

Interpretation of $\lambda_i$: outlier observations get low $\lambda_i$ (high variance $\sigma^2/\lambda_i$), downweighting their influence on $\beta$.

Structure:

Section	Topic
Section 1	Synthetic contaminated data: recover $\nu$, identify outliers, compare with Normal
Section 2	Nelson-Plosser (1982) macroeconomic series: posterior for $\nu$ across 14 series
Section 3	PyMC equivalent: continuous $\nu$ via NUTS

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

import os, sys

# Windows DLL fix for conda environments
conda_lib = os.path.join(os.environ.get('CONDA_PREFIX', ''), 'Library', 'bin')
if os.path.isdir(conda_lib) and conda_lib not in os.environ.get('PATH', ''):
    os.environ['PATH'] = conda_lib + ';' + os.environ.get('PATH', '')

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import time
from scipy.stats import t as t_dist, norm as norm_dist

sys.path.insert(0, r'c:\Users\user\project')
from student_t_lm_gibbs import (
    student_t_gibbs, normal_gibbs, posterior_summary,
    log_lik_t, log_lik_normal
)
print('Import OK')

Import OK

---

1. Synthetic contaminated data

Design: $k=3$, $n=120$, $\beta^\star=[2, 0.8, -0.5]$, $\sigma^2=1$, errors drawn from $t(\nu^\star=4)$.

A $t(4)$ distribution has heavy tails — about 6% of draws fall beyond $\pm 3\sigma$ (vs 0.3% for Normal). This mimics real data with occasional large shocks.

We compare:

• Student-t Gibbs (our sampler): correctly models heavy tails, recovers $\nu \approx 4$
• Normal Gibbs (baseline): assumes Gaussian errors, biased by outliers

Theorem 3 (Geweke 1993, improper flat prior on $\beta$, $\sigma^2$): Posterior moments exist for $\beta$ when $\nu > 1$, and for $\sigma^2$ when $\nu > 2$. Since our discrete grid starts at $\nu = 2$, posterior moments for $\sigma^2$ are defined everywhere on the grid; for $\beta$ they are always defined.

# ── Generate data ─────────────────────────────────────────────────────────────
rng1 = np.random.default_rng(0)
n, k      = 120, 3
beta_true = np.array([2.0, 0.8, -0.5])
sigma2_true, nu_true = 1.0, 4.0

X1 = np.column_stack([np.ones(n), rng1.standard_normal((n, k - 1))])
eps = t_dist.rvs(df=nu_true, scale=np.sqrt(sigma2_true), size=n,
                 random_state=rng1)
y1 = X1 @ beta_true + eps

# OLS for reference
beta_ols = np.linalg.lstsq(X1, y1, rcond=None)[0]
resid_ols = y1 - X1 @ beta_ols

print(f'True β:    {beta_true}   σ²={sigma2_true}  ν={nu_true}')
print(f'OLS  β:    {beta_ols.round(3)}')
print(f'|OLS error|: {np.abs(beta_ols - beta_true).round(3)}')
print(f'\nResidual range: [{resid_ols.min():.2f}, {resid_ols.max():.2f}]')
print(f'Obs beyond ±3σ: {(np.abs(resid_ols) > 3).sum()} / {n}')

True β:    [ 2.   0.8 -0.5]   σ²=1.0  ν=4.0
OLS  β:    [ 2.061  0.745 -0.675]
|OLS error|: [0.061 0.055 0.175]

Residual range: [-3.60, 4.88]
Obs beyond ±3σ: 4 / 120

# ── Run both samplers ─────────────────────────────────────────────────────────
out_t = student_t_gibbs(y1, X1, R=15000, burn=3000, seed=1, nprint=0)
out_n = normal_gibbs(   y1, X1, R=15000, burn=3000, seed=2, nprint=0)

names1 = ['β[0]', 'β[1]', 'β[2]']
summ_t = posterior_summary(out_t['beta'], names1, true_vals=beta_true)
summ_n = posterior_summary(out_n['beta'], names1, true_vals=beta_true)

print('Student-t posterior (our Gibbs):')
print(summ_t)
print('\nNormal posterior (baseline):')
print(summ_n)

print(f'\nStudent-t  ν posterior: mean={out_t["nu"].mean():.1f}  '
      f'median={np.median(out_t["nu"]):.0f}  '
      f'mode={float(pd.Series(out_t["nu"]).mode().iloc[0]):.0f}'
      f'   (true ν={nu_true:.0f})')

# Posterior mean log-likelihood comparison
ll_t, se_t = log_lik_t(y1, X1, out_t['beta'], out_t['sigma2'], out_t['nu'])
ll_n, se_n = log_lik_normal(y1, X1, out_n['beta'], out_n['sigma2'])
print(f'\nPosterior mean log-lik:  Student-t = {ll_t:.2f} (±{se_t:.2f})'
      f'   Normal = {ll_n:.2f} (±{se_n:.2f})'
      f'   Δ = {ll_t - ll_n:.2f}')

Done in 0.0 min  |  kept draws: 12000
Student-t posterior (our Gibbs):
        mean      sd    q025    q975  true  covered
β[0]  2.0284  0.1109  1.8118  2.2434   2.0     True
β[1]  0.7939  0.1162  0.5665  1.0268   0.8     True
β[2] -0.6628  0.1082 -0.8753 -0.4507  -0.5     True

Normal posterior (baseline):
        mean      sd    q025    q975  true  covered
β[0]  2.0584  0.1216  1.8178  2.2947   2.0     True
β[1]  0.7449  0.1234  0.5036  0.9894   0.8     True
β[2] -0.6752  0.1174 -0.9081 -0.4445  -0.5     True

Student-t  ν posterior: mean=8.6  median=6  mode=4   (true ν=4)

Posterior mean log-lik:  Student-t = -197.90 (±1.83)   Normal = -202.39 (±1.42)   Δ = 4.49

fig, axes = plt.subplots(2, 3, figsize=(14, 7))

# Row 0: β posteriors — t vs Normal
for j, ax in enumerate(axes[0]):
    all_v = np.concatenate([out_t['beta'][:, j], out_n['beta'][:, j]])
    bins  = np.linspace(np.percentile(all_v, 0.5), np.percentile(all_v, 99.5), 60)
    ax.hist(out_n['beta'][:, j], bins=bins, density=True, alpha=0.5,
            color='steelblue', label='Normal Gibbs')
    ax.hist(out_t['beta'][:, j], bins=bins, density=True, alpha=0.6,
            color='orange', label='Student-t Gibbs')
    ax.axvline(beta_true[j], color='k', ls='--', lw=1.5, label=f'True={beta_true[j]}')
    ax.set_title(names1[j]); ax.legend(fontsize=7)

# Row 1, left: ν posterior
ax = axes[1, 0]
nu_vals, nu_counts = np.unique(out_t['nu'], return_counts=True)
ax.bar(nu_vals, nu_counts / len(out_t['nu']), width=0.8, color='orange', alpha=0.8)
ax.axvline(nu_true, color='k', ls='--', lw=1.5, label=f'True ν={nu_true:.0f}')
ax.set_xlabel('ν'); ax.set_ylabel('Posterior probability')
ax.set_title('ν posterior'); ax.legend()

# Row 1, middle: λᵢ vs outlier status
ax = axes[1, 1]
lam_mean = out_t['lam'].mean(axis=0)
outlier  = np.abs(eps) > 2.5  # "true" outliers (large errors)
ax.scatter(np.arange(n)[~outlier], lam_mean[~outlier], s=20,
           color='steelblue', alpha=0.6, label='Normal obs')
ax.scatter(np.arange(n)[outlier],  lam_mean[outlier],  s=40,
           color='red', marker='x', lw=1.5, label='|εᵢ| > 2.5σ')
ax.axhline(1.0, color='k', ls=':', lw=0.8)
ax.set_xlabel('Observation index'); ax.set_ylabel('E[λᵢ | y]')
ax.set_title('Mixing weights λᵢ (outliers → low λ)')
ax.legend(fontsize=8)

# Row 1, right: residuals vs 1/λᵢ
ax = axes[1, 2]
ax.scatter(np.abs(eps), lam_mean, s=20, alpha=0.6, color='orange')
ax.set_xlabel('|true error εᵢ|'); ax.set_ylabel('E[λᵢ | y]')
ax.set_title('|error| vs posterior λᵢ (inverse relationship)')

fig.suptitle(
    f'§1 Synthetic t({nu_true:.0f}) data — Student-t Gibbs recovers true ν\n'
    f'ν posterior mode = {float(pd.Series(out_t["nu"]).mode().iloc[0]):.0f}'
    f'  (true {nu_true:.0f})',
    fontsize=11
)
plt.tight_layout(); plt.show()

Reading the panels. Top row — the posterior of each coefficient under the Student-t (orange) and Normal (blue) samplers. Both cover the truth (dashed), but the Student-t is slightly tighter because it stops the four large errors from inflating the spread. Bottom-left — the posterior of $\nu$ concentrates at its true value (mode $\approx 4$). Bottom-middle — the estimated mixing weights $\mathbb{E}[\lambda_i\mid y]$: the planted large-error observations (red ×) are pushed to low $\lambda$, i.e. automatically down-weighted. Bottom-right — $\lambda_i$ falls monotonically as $|\varepsilon_i|$ grows: the mechanism in one picture, the bigger the shock the less it is allowed to move $\beta$.

The robustness payoff

# Robustness: OLS/Normal regression line dragged by outliers vs Student-t line that resists them
rngR = np.random.default_rng(7); nR = 60
xR = np.sort(rngR.uniform(-2, 2, nR)); yR = 1.0 + 2.0*xR + 0.6*rngR.standard_normal(nR)
oi = np.array([8, 25, 40, 52]); yR[oi] += np.array([7.5, -8.0, 9.0, -7.0])     # gross outliers
XR = np.column_stack([np.ones(nR), xR])
fT = student_t_gibbs(yR, XR, R=8000, burn=2000, seed=11, nprint=0)
fN = normal_gibbs(yR, XR, R=8000, burn=2000, seed=12, nprint=0)
bT = fT['beta'].mean(0); bN = fN['beta'].mean(0); lamR = fT['lam'].mean(0)
fig, ax = plt.subplots(figsize=(9, 5.5))
sc = ax.scatter(xR, yR, c=lamR, cmap='viridis', s=45, edgecolor='k', lw=.3)
plt.colorbar(sc, label='posterior λ  (low = auto-downweighted outlier)')
gx = np.array([xR.min(), xR.max()])
ax.plot(gx, bN[0]+bN[1]*gx, color='steelblue', lw=2.5, label=f'Normal / OLS  (slope {bN[1]:.2f})')
ax.plot(gx, bT[0]+bT[1]*gx, color='firebrick', lw=2.5, label=f'Student-t   (slope {bT[1]:.2f})')
ax.plot(gx, 1.0+2.0*gx, color='black', ls='--', lw=1.2, label='truth (slope 2.00)')
ax.set_xlabel('x'); ax.set_ylabel('y'); ax.set_title('Robust regression: the Normal fit is pulled by outliers, the Student-t fit holds the bulk')
ax.legend(fontsize=9)
plt.tight_layout(); plt.savefig('student_t_robust.png', dpi=120, bbox_inches='tight'); plt.show()

Done in 0.0 min  |  kept draws: 6000

The robustness payoff. Four gross outliers are added to an otherwise clean line. The Normal / OLS fit (blue) is visibly tilted and shifted toward them, while the Student-t fit (red) stays on the true relationship (slope $\approx 2$, black dashed). The point colour is the posterior weight $\lambda$: the four outliers are the darkest points — the model has identified and down-weighted them on its own, with no manual flagging or robust-loss tuning.

Cost of ignoring tails: σ² inflation and a posterior-predictive tail check

# sigma^2 inflation + posterior-predictive tail check on the section-1 data
fig, ax = plt.subplots(1, 2, figsize=(13, 4.5))
# (A) sigma^2: Normal absorbs the heavy tails into the variance; Student-t keeps a small SCALE
ax[0].hist(out_n['sigma2'], bins=50, density=True, alpha=.5, color='steelblue', label=f'Normal  (mean {out_n["sigma2"].mean():.2f})')
ax[0].hist(out_t['sigma2'], bins=50, density=True, alpha=.6, color='firebrick', label=f'Student-t (mean {out_t["sigma2"].mean():.2f})')
ax[0].axvline(sigma2_true, color='k', ls='--', label=f'true scale σ²={sigma2_true:.0f}')
ax[0].set_xlabel('σ²'); ax[0].set_yticks([]); ax[0].set_title('σ²: Normal inflates (variance≈ν/(ν-2)); t recovers the scale'); ax[0].legend(fontsize=8)
# (B) PPC: statistic = max|residual| / scale; replicate under each model
rng = np.random.default_rng(0); S = 1500
def tstat(res, sc): return np.max(np.abs(res))/sc
idn = rng.choice(len(out_n['sigma2']), S, replace=False); idt = rng.choice(len(out_t['nu']), S, replace=False)
Tn_o, Tn_r, Tt_o, Tt_r = [], [], [], []
for r in idn:
    s = np.sqrt(out_n['sigma2'][r]); mu = X1 @ out_n['beta'][r]
    Tn_o.append(tstat(y1-mu, s)); Tn_r.append(tstat(s*rng.standard_normal(n), s))
for r in idt:
    s = np.sqrt(out_t['sigma2'][r]); nu = out_t['nu'][r]; mu = X1 @ out_t['beta'][r]
    Tt_o.append(tstat(y1-mu, s)); Tt_r.append(tstat(s*t_dist.rvs(df=nu, size=n, random_state=rng), s))
Tn_o, Tn_r, Tt_o, Tt_r = map(np.array, (Tn_o, Tn_r, Tt_o, Tt_r))
pn = (Tn_r >= Tn_o).mean(); pt = (Tt_r >= Tt_o).mean()
# the t replicates reach into the hundreds while both bulks sit below ~10, so a linear axis
# leaves the panel almost entirely empty -- bin geometrically and plot on a log scale instead
lo = max(min(Tn_r.min(), Tt_r.min()), 1e-2); hi = max(Tn_r.max(), Tt_r.max())
bins = np.geomspace(lo, hi*1.05, 45)
ax[1].hist(Tn_r, bins=bins, density=True, alpha=.5, color='steelblue', label=f'Normal y_rep  (p={pn:.3f})')
ax[1].hist(Tt_r, bins=bins, density=True, alpha=.5, color='firebrick', label=f'Student-t y_rep (p={pt:.2f})')
ax[1].axvline(Tn_o.mean(), color='steelblue', ls='--', lw=2, label=f'observed, Normal fit ({Tn_o.mean():.1f})')
ax[1].axvline(Tt_o.mean(), color='firebrick', ls='--', lw=2, label=f'observed, t fit ({Tt_o.mean():.1f})')
ax[1].set_xscale('log')
ax[1].set_xlabel('max |residual| / scale   (log axis)'); ax[1].set_yticks([])
ax[1].set_title('Posterior-predictive tail check'); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.savefig('student_t_ppc.png', dpi=120, bbox_inches='tight'); plt.show()
print('PPC Bayesian p-value (max|resid|/scale):  Normal=%.3f (model rejected)   Student-t=%.2f (adequate)' % (pn, pt))
print(f'The replicate distributions are plotted on a log axis because they span '
      f'{Tt_r.max()/Tt_r.min():.0f}x: the Student-t replicates reach {Tt_r.max():.0f} while both bulks sit '
      f'near {np.median(np.r_[Tn_r, Tt_r]):.1f}. That spread is the diagnostic working rather than a nuisance -- a '
      f'model with genuinely heavy tails must be ABLE to produce the occasional huge standardised '
      f'residual, and the Normal cannot, which is why its observed value of {Tn_o.mean():.1f} sits so far '
      f'out in its own replicate distribution.')

PPC Bayesian p-value (max|resid|/scale):  Normal=0.043 (model rejected)   Student-t=0.41 (adequate)
The replicate distributions are plotted on a log axis because they span 87x: the Student-t replicates reach 197 while both bulks sit near 3.2. That spread is the diagnostic working rather than a nuisance -- a model with genuinely heavy tails must be ABLE to produce the occasional huge standardised residual, and the Normal cannot, which is why its observed value of 3.7 sits so far out in its own replicate distribution.

Two costs of forcing a Normal. Left — the error-variance posteriors: to swallow the outliers a thin-tailed model must widen its variance, so the Normal's $\sigma^2$ sits near 2 ($\approx \nu/(\nu-2)$) while the Student-t recovers the true scale of 1. Right — a posterior-predictive check on the tail statistic $\max|\text{residual}|/\text{scale}$: datasets replicated from the Normal (blue) rarely reach the observed value (dashed), giving a Bayesian $p \approx 0.04$ — the Normal is rejected — whereas the Student-t (red) reproduces the observed tails and is adequate.

Residual diagnostics and pointwise model comparison (WAIC)

# QQ-plot of standardized residuals + pointwise WAIC contribution (t vs Normal)
from scipy import stats as sps
from scipy.special import gammaln
bT1 = out_t['beta'].mean(0); s2T = out_t['sigma2'].mean()
std_res = (y1 - X1 @ bT1) / np.sqrt(s2T)
fig, ax = plt.subplots(1, 2, figsize=(13, 4.5))
sps.probplot(std_res, dist='norm', plot=ax[0]); ax[0].get_lines()[0].set_markersize(4)
ax[0].set_title('Normal QQ-plot of residuals (S-shape ⇒ heavy tails)')

# vectorised pointwise WAIC: elpd_i = lppd_i - p_waic_i  (logsumexp-mean over draws)
def lse_mean(ll):                                   # log( mean_r exp(ll_ri) ) per obs, ll shape (R,n)
    m = ll.max(0); return np.log(np.exp(ll - m).mean(0)) + m
ssub = np.random.default_rng(0).choice(len(out_t['nu']), 3000, replace=False)
Bt = out_t['beta'][ssub]; st = np.sqrt(out_t['sigma2'][ssub]); nt = out_t['nu'][ssub]
Et = y1[None,:] - Bt @ X1.T                          # (3000, n)
z2 = (Et / st[:,None])**2
ll_t = (gammaln((nt+1)/2)-gammaln(nt/2)-0.5*np.log(nt*np.pi))[:,None] - np.log(st)[:,None] - (nt[:,None]+1)/2*np.log1p(z2/nt[:,None])
el_t = lse_mean(ll_t) - ll_t.var(0)
ssn = np.random.default_rng(1).choice(len(out_n['sigma2']), 3000, replace=False)
Bn = out_n['beta'][ssn]; sn = np.sqrt(out_n['sigma2'][ssn]); En = y1[None,:] - Bn @ X1.T
ll_n = -0.5*np.log(2*np.pi) - np.log(sn)[:,None] - 0.5*(En/sn[:,None])**2
el_n = lse_mean(ll_n) - ll_n.var(0)
diff = el_t - el_n; order = np.argsort(diff)
ax[1].bar(range(n), diff[order], color=np.where(diff[order] > 0, 'firebrick', 'steelblue'))
ax[1].axhline(0, color='k', lw=.5); ax[1].set_xlabel('observation (sorted)'); ax[1].set_ylabel('elpd_waic(t) − elpd_waic(Normal)')
ax[1].set_title('Pointwise WAIC gain — t wins on the outlier observations (red)')
plt.tight_layout(); plt.savefig('student_t_diag.png', dpi=120, bbox_inches='tight'); plt.show()
se = np.sqrt(n) * diff.std()
print('WAIC elpd:  Student-t=%.1f   Normal=%.1f   diff=%.1f ± %.1f  (t favored)' % (el_t.sum(), el_n.sum(), el_t.sum()-el_n.sum(), se))

WAIC elpd:  Student-t=-200.9   Normal=-205.1   diff=4.2 ± 2.9  (t favored)

Diagnostics. Left — a normal Q–Q plot of the standardized residuals: the S-shape (points curling above the line at the top and below it at the bottom) is the fingerprint of heavy tails that a Gaussian model cannot accommodate. Right — the pointwise WAIC difference, observation by observation: the Student-t's predictive-density advantage (red bars) is concentrated almost entirely on the outlying points, and summed over the sample the WAIC prefers the Student-t.

---

2. Nelson-Plosser (1982) macroeconomic series

Geweke (1993) fits a linear trend model with Student-t errors to each of the 14 annual series: $$y_t = \alpha + \beta \cdot t + \varepsilon_t, \qquad \varepsilon_t \sim t(\nu,\, 0,\, \sigma^2)$$

where $y_t$ is the log of the series (levels for unemployment rate and interest rate).

Key result (Geweke 1993, Table 2): posterior odds strongly favour the Student-t over the Normal for 13 of 14 series, with posterior $\nu$ estimates typically in the range 3–8, indicating substantial heavy-tailedness in macroeconomic time series.

Interpretation of $\lambda_t$: observations receiving low weight ($\lambda_t \ll 1$) correspond to historically unusual years — Great Depression, WWII, oil shocks. The Student-t model down-weights these years automatically rather than letting them dominate the trend estimate.

# ── Load Nelson-Plosser data ──────────────────────────────────────────────────
data_path = 'nelplo.csv'   # sits beside the notebook; the previous absolute path no longer exists
np_df = pd.read_csv(data_path)

# Series to log-transform (per Geweke/Nelson-Plosser convention)
log_series = ['cpi','ip','gnp.nom','vel','emp','nom.wages',
              'gnp.def','money.stock','gnp.real','stock.prices',
              'gnp.capita','real.wages']
level_series = ['int.rate', 'unemp']   # keep in levels

series_labels = {
    'cpi':         'CPI',
    'ip':          'Industrial production',
    'gnp.nom':     'Nominal GNP',
    'vel':         'Velocity',
    'emp':         'Employment',
    'int.rate':    'Interest rate',
    'nom.wages':   'Nominal wages',
    'gnp.def':     'GNP deflator',
    'money.stock': 'Money stock',
    'gnp.real':    'Real GNP',
    'stock.prices':'Stock prices',
    'gnp.capita':  'GNP per capita',
    'real.wages':  'Real wages',
    'unemp':       'Unemployment rate',
}
all_series = list(series_labels.keys())

# Apply log transform and restrict to 1909-1970 (Geweke's sample window)
np_df = np_df[(np_df['year'] >= 1909) & (np_df['year'] <= 1970)].copy()
for s in log_series:
    np_df[s] = np.log(np_df[s])

print(f'Sample: {np_df["year"].min()}–{np_df["year"].max()}  '
      f'({len(np_df)} obs)')
print(f'\nNon-missing obs per series:')
for s in all_series:
    nn = np_df[s].notna().sum()
    print(f'  {series_labels[s]:25s}: {nn}')

Sample: 1909–1970  (62 obs)

Non-missing obs per series:
  CPI                      : 62
  Industrial production    : 62
  Nominal GNP              : 62
  Velocity                 : 62
  Employment               : 62
  Interest rate            : 62
  Nominal wages            : 62
  GNP deflator             : 62
  Money stock              : 62
  Real GNP                 : 62
  Stock prices             : 62
  GNP per capita           : 62
  Real wages               : 62
  Unemployment rate        : 62

# ── Fit Student-t trend model to each series ──────────────────────────────────
# Model: y_t = α + β·t + ε_t,  ε_t ~ t(ν, 0, σ²)
results_np = {}

t0 = time.time()
for si, s in enumerate(all_series):
    sub = np_df[['year', s]].dropna()
    y_s = sub[s].values
    t_s = (sub['year'].values - sub['year'].mean()) / sub['year'].std()  # standardise time
    X_s = np.column_stack([np.ones(len(y_s)), t_s])

    out = student_t_gibbs(y_s, X_s, R=12000, burn=3000, seed=1000 + si,
                          nprint=0)
    out_norm = normal_gibbs(y_s, X_s, R=12000, burn=3000, seed=5000 + si,
                            nprint=0)

    ll_t_v, _ = log_lik_t(y_s, X_s, out['beta'], out['sigma2'], out['nu'])
    ll_n_v, _ = log_lik_normal(y_s, X_s, out_norm['beta'], out_norm['sigma2'])

    results_np[s] = {
        'nu_mean':   out['nu'].mean(),
        'nu_median': np.median(out['nu']),
        'nu_mode':   float(pd.Series(out['nu']).mode().iloc[0]),
        'P_nu_gt20': (out['nu'] > 20).mean(),   # P(approximately Normal | y)
        'll_t': ll_t_v, 'll_n': ll_n_v,
        'delta_ll': ll_t_v - ll_n_v,
        'n': len(y_s),
        'lam_mean': out['lam'].mean(axis=0),
        'years': sub['year'].values,
    }
    print(f'{series_labels[s]:25s}  n={len(y_s):2d}  '
          f'ν mean={out["nu"].mean():.1f}  mode={results_np[s]["nu_mode"]:.0f}'
          f'  Δlog-lik={ll_t_v - ll_n_v:+.1f}')

print(f'\nTotal time: {time.time()-t0:.1f}s')

Done in 0.0 min  |  kept draws: 9000

CPI                        n=62  ν mean=23.7  mode=36  Δlog-lik=-0.4

Done in 0.0 min  |  kept draws: 9000

Industrial production      n=62  ν mean=2.6  mode=2  Δlog-lik=+9.3

Done in 0.0 min  |  kept draws: 9000

Nominal GNP                n=62  ν mean=3.0  mode=2  Δlog-lik=+6.7

Done in 0.0 min  |  kept draws: 9000

Velocity                   n=62  ν mean=5.5  mode=2  Δlog-lik=+5.9

Done in 0.0 min  |  kept draws: 9000

Employment                 n=62  ν mean=3.1  mode=2  Δlog-lik=+7.0

Done in 0.0 min  |  kept draws: 9000

Interest rate              n=62  ν mean=18.1  mode=7  Δlog-lik=+1.1

Done in 0.0 min  |  kept draws: 9000

Nominal wages              n=62  ν mean=21.3  mode=23  Δlog-lik=-0.0

Done in 0.0 min  |  kept draws: 9000

GNP deflator               n=62  ν mean=21.8  mode=37  Δlog-lik=-0.1

Done in 0.0 min  |  kept draws: 9000

Money stock                n=62  ν mean=26.9  mode=39  Δlog-lik=-0.9

Done in 0.0 min  |  kept draws: 9000

Real GNP                   n=62  ν mean=2.7  mode=2  Δlog-lik=+7.9

Done in 0.0 min  |  kept draws: 9000

Stock prices               n=62  ν mean=26.0  mode=39  Δlog-lik=-0.4

Done in 0.0 min  |  kept draws: 9000

GNP per capita             n=62  ν mean=3.3  mode=2  Δlog-lik=+6.4

Done in 0.0 min  |  kept draws: 9000

Real wages                 n=62  ν mean=15.4  mode=3  Δlog-lik=+0.8

Done in 0.0 min  |  kept draws: 9000

Unemployment rate          n=62  ν mean=20.2  mode=4  Δlog-lik=-0.0

Total time: 18.5s

# ── Summary table ─────────────────────────────────────────────────────────────
rows = []
for s, r in results_np.items():
    rows.append({
        'Series':         series_labels[s],
        'n':              r['n'],
        'ν mode':         r['nu_mode'],
        'ν mean':         round(r['nu_mean'], 1),
        'P(ν>20)':        round(r['P_nu_gt20'], 3),
        'Δ log-lik (t−N)': round(r['delta_ll'], 1),
        't preferred':    '✓' if r['delta_ll'] > 0 else '✗',
    })
tbl = pd.DataFrame(rows).set_index('Series')
print(tbl.to_string())
print(f'\nSeries preferring Student-t: {(tbl["Δ log-lik (t−N)"] > 0).sum()} / {len(tbl)}'
      f'  (Geweke 1993: 13/14)')

                        n  ν mode  ν mean  P(ν>20)  Δ log-lik (t−N) t preferred
Series                                                                         
CPI                    62    36.0    23.7    0.611             -0.4           ✗
Industrial production  62     2.0     2.6    0.004              9.3           ✓
Nominal GNP            62     2.0     3.0    0.019              6.7           ✓
Velocity               62     2.0     5.5    0.054              5.9           ✓
Employment             62     2.0     3.1    0.013              7.0           ✓
Interest rate          62     7.0    18.1    0.379              1.1           ✓
Nominal wages          62    23.0    21.3    0.524             -0.0           ✗
GNP deflator           62    37.0    21.8    0.536             -0.1           ✗
Money stock            62    39.0    26.9    0.737             -0.9           ✗
Real GNP               62     2.0     2.7    0.000              7.9           ✓
Stock prices           62    39.0    26.0    0.702             -0.4           ✗
GNP per capita         62     2.0     3.3    0.005              6.4           ✓
Real wages             62     3.0    15.4    0.319              0.8           ✓
Unemployment rate      62     4.0    20.2    0.484             -0.0           ✗

Series preferring Student-t: 8 / 14  (Geweke 1993: 13/14)

Result: 8 of 14 series clearly prefer Student-t, with 2 more (Nominal wages, Unemployment rate) essentially tied at Δlog-lik ≈ 0 — so the count is 8–9 depending on Monte-Carlo noise (Geweke 1993 reports 13/14 using exact posterior odds). The qualitative story is the same — many macroeconomic series show heavy-tailed residuals — but the borderline series disagree. Two likely causes:

1. Data source. We use tseries::NelPlo (the R package version), which carries post-1982 revisions. Geweke used the original Nelson-Plosser (1982) NBER data. The five series that flip (CPI, GNP deflator, Money stock, Stock prices, Unemployment) all have Δ ≈ 0, so minor data differences are enough to tip the sign.

2. ν grid upper bound. Our grid goes to 40; Geweke likely capped at 30. The four series with posterior mode ≥ 30 (CPI, deflator, money, stocks) are the ones most affected — a tighter cap concentrates their posterior away from the Normal limit and slightly increases Δ log-lik.

What agrees with Geweke: the strongly heavy-tailed series — Industrial production, Real GNP, Employment, Nominal GNP, Velocity, GNP per capita — all have ν mode = 2 and Δ log-lik ≈ 6–9. These are the economically interesting cases corresponding to Depression- and WWII-era shocks.

Connection to stationarity testing (Geweke 1993, Table V). Geweke shows that lower posterior $\nu$ systematically reduces posterior odds in favour of difference stationarity relative to trend stationarity, and shrinks the posterior SD of the trend coefficient $\beta$.

The intuition is closely related to Perron (1989): observations from the Great Depression, WWII, and oil shocks create large deviations from a linear trend that look like a unit root under a Normal error model. The Student-t model instead assigns those observations low mixing weights $\lambda_t \ll 1$, absorbing the shock rather than attributing it to a stochastic trend. With outliers down-weighted, the trend coefficient is estimated more precisely (smaller posterior SD) and the data look more consistent with trend stationarity.

So while the Student-t model is not a formal unit root test, jointly estimating $\nu$ with the trend provides an implicit diagnostic: low posterior $\nu$ signals that apparent non-stationarity may be driven by heavy-tailed shocks rather than a genuine unit root. This is one of the deepest results in the paper. The Nelson-Plosser application is a useful reminder that unit root conclusions are not robust to tail assumptions, a point often neglected in applied work that assumes Gaussian disturbances.

fig, axes = plt.subplots(3, 5, figsize=(18, 10))
axes_flat = axes.flat

# ν posterior bar charts for each series
for ax, s in zip(axes_flat, all_series):
    # Rerun sampler to get full nu draws (stored results only have summaries)
    sub = np_df[['year', s]].dropna()
    y_s = sub[s].values
    t_s = (sub['year'].values - sub['year'].mean()) / sub['year'].std()
    X_s = np.column_stack([np.ones(len(y_s)), t_s])
    out = student_t_gibbs(y_s, X_s, R=12000, burn=3000,
                          seed=hash(s) % (2**31), nprint=0)
    nu_vals, nu_cnt = np.unique(out['nu'], return_counts=True)
    ax.bar(nu_vals, nu_cnt / len(out['nu']), width=0.8, color='orange', alpha=0.8)
    ax.axvline(results_np[s]['nu_mode'], color='k', ls='--', lw=1, alpha=0.7)
    ax.set_title(f'{series_labels[s]}\nmode={results_np[s]["nu_mode"]:.0f}'
                 f'  Δ={results_np[s]["delta_ll"]:+.1f}', fontsize=8)
    ax.set_xlim(1, 42); ax.set_xticks([5, 15, 25, 35])

# Hide unused panels (14 series, 15 subplots)
for ax in list(axes_flat)[14:]:
    ax.set_visible(False)

fig.suptitle('§2 Nelson-Plosser — posterior ν for each series (orange bars)\n'
             'Δ log-lik > 0 means Student-t fits better than Normal', fontsize=12)
plt.tight_layout(); plt.show()

Done in 0.0 min  |  kept draws: 9000

Done in 0.0 min  |  kept draws: 9000

Done in 0.0 min  |  kept draws: 9000

Done in 0.0 min  |  kept draws: 9000

Done in 0.0 min  |  kept draws: 9000

Done in 0.0 min  |  kept draws: 9000

Done in 0.0 min  |  kept draws: 9000

Done in 0.0 min  |  kept draws: 9000

Done in 0.0 min  |  kept draws: 9000

Done in 0.0 min  |  kept draws: 9000

Done in 0.0 min  |  kept draws: 9000

Done in 0.0 min  |  kept draws: 9000

Done in 0.0 min  |  kept draws: 9000

Done in 0.0 min  |  kept draws: 9000

$\nu$ posteriors, series by series. Each panel is the posterior of $\nu$ for one series. The economically turbulent series — Industrial Production, Real GNP, Employment, Nominal GNP, Velocity, GNP per capita — pile their mass at the grid floor ($\nu = 2$), i.e. strongly heavy-tailed, with a log-likelihood gain of $\approx 6\text{–}9$ over the Normal. The smoother nominal series (money stock, GNP deflator, CPI, stock prices) spread toward large $\nu$, where the Student-t is indistinguishable from the Normal and $\Delta$ log-lik $\approx 0$.

# ── λₜ over time for Real GNP — identifying historically unusual years ────────
s_focus = 'gnp.real'
sub_f   = np_df[['year', s_focus]].dropna()
y_f     = sub_f[s_focus].values
t_f     = (sub_f['year'].values - sub_f['year'].mean()) / sub_f['year'].std()
X_f     = np.column_stack([np.ones(len(y_f)), t_f])
out_f   = student_t_gibbs(y_f, X_f, R=15000, burn=5000,
                           seed=777, nprint=0)

lam_f   = out_f['lam'].mean(axis=0)
years_f = sub_f['year'].values

fig, ax = plt.subplots(figsize=(13, 4))
ax.plot(years_f, lam_f, 'o-', ms=4, color='orange', lw=1)
ax.axhline(1.0, color='k', ls=':', lw=0.8)
ax.fill_between(years_f, lam_f, 1.0, where=lam_f < 1.0,
                alpha=0.3, color='red', label='λₜ < 1 (outlier)')
ax.fill_between(years_f, lam_f, 1.0, where=lam_f > 1.0,
                alpha=0.2, color='green', label='λₜ > 1 (inlier)')

# Annotate key events
threshold = lam_f.mean() - 1.5 * lam_f.std()
for yr, lv in zip(years_f, lam_f):
    if lv < threshold:
        ax.annotate(str(yr), (yr, lv), textcoords='offset points',
                    xytext=(0, -12), fontsize=7, ha='center', color='red')

ax.set_xlabel('Year'); ax.set_ylabel('E[λₜ | y]')
ax.set_title('Real GNP — posterior mixing weights λₜ over time\n'
             'Low λₜ = large residual, observation down-weighted by Student-t')
ax.legend(fontsize=8)
plt.tight_layout(); plt.show()

Done in 0.0 min  |  kept draws: 10000

Which years get down-weighted. For Real GNP, the posterior weight $\lambda_t$ dips well below 1 (shaded red) in exactly the historically turbulent years — the early-1930s Depression collapse and the wartime / immediate-postwar swings (annotated). These are the observations the Student-t treats as tail draws rather than as evidence about the trend, which is why the fitted trend is estimated more precisely than under a Normal that lets those same years pull it around.

Macro history in the outlier weights — λ across all series

(dotted red lines mark 1929–33 Great Depression and 1941–45 WWII)

# Nelson-Plosser: posterior outlier weights λ by series x year -- macro history in the weights
years_all = np.arange(int(np_df['year'].min()), int(np_df['year'].max())+1)
M = np.full((len(all_series), len(years_all)), np.nan)
for si, s in enumerate(all_series):
    for yy, lm in zip(results_np[s]['years'], results_np[s]['lam_mean']):
        M[si, int(yy)-years_all[0]] = lm
fig, ax = plt.subplots(figsize=(16, 6))
im = ax.imshow(M, aspect='auto', cmap='viridis', vmin=0, vmax=1.8)
plt.colorbar(im, label='posterior mean λ  (dark = down-weighted = unusual year)')
ax.set_yticks(range(len(all_series))); ax.set_yticklabels([series_labels[s] for s in all_series], fontsize=8)
xt = np.arange(0, len(years_all), 10); ax.set_xticks(xt); ax.set_xticklabels(years_all[xt])
for span,lab in [((1929,1933),'Depression'), ((1941,1945),'WWII')]:
    for yv in span: ax.axvline(yv-years_all[0], color='red', lw=.6, ls=':')
ax.set_xlabel('year'); ax.set_title('Nelson-Plosser Student-t outlier weights λₜ across series (low λ = historically unusual year)')
plt.tight_layout(); plt.savefig('student_t_heatmap.png', dpi=120, bbox_inches='tight'); plt.show()

Macro history in the weights. Each row is a series and each column a year; darker cells are lower posterior weights ($\lambda_t \ll 1$). The dark columns line up across many series on the Great Depression (1929–33) and WWII (1941–45) (dotted red) — in effect the heavy-tail machinery dates the century's major macroeconomic disruptions and quarantines them, so they do not distort the trend estimates. This is the visual counterpart to the stationarity result above: the "unit-root-looking" deviations are concentrated in a handful of identifiable years.

---

Summary

Geweke (1993) key contributions

1. Scale-mixture representation — $\varepsilon_i\sim t(\nu,0,\sigma^2)$ written as $\varepsilon_i\mid\lambda_i\sim N(0,\sigma^2/\lambda_i)$, $\lambda_i\sim\text{Gamma}(\nu/2,\nu/2)$, makes all four full conditionals conjugate → an efficient 4-block Gibbs sampler.
2. Discrete ν — a discrete uniform prior on {2,…,40}, drawn as a categorical each sweep: fast, simple, stable near ν=2.
3. Automatic outlier detection — posterior $\lambda_i$ is an inverse outlier weight: unusual observations get $\lambda_i\ll1$, limiting their leverage on β without manual flagging.

Results

	Section 1 Synthetic	Section 2 Nelson-Plosser
true / expected ν	4	3–8 (Geweke)
Gibbs mode	≈ 4 ✓	varies by series
series favouring t	—	8 clear + 2 tied → 8–9/14 (Geweke: 13/14)
robustness	t line resists outliers; Normal dragged	λ-heatmap flags Depression/WWII
model check	PPC: Normal rejected (p≈0.04), t adequate; WAIC favours t	—

Continuous-ν comparison (PyMC / NUTS): see the companion notebook student_t_lm_pymc.ipynb.