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()
ax[1].hist(Tn_r, bins=40, density=True, alpha=.5, color='steelblue', label=f'Normal y_rep  (p={pn:.3f})')
ax[1].hist(Tt_r, bins=40, density=True, alpha=.5, color='firebrick', label=f'Student-t y_rep (p={pt:.2f})')
ax[1].axvline(Tn_o.mean(), color='k', ls='--', label='observed')
ax[1].set_xlabel('max |residual| / scale'); 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))

PPC Bayesian p-value (max|resid|/scale):  Normal=0.043 (model rejected)   Student-t=0.41 (adequate)

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 = r'c:\Users\user\project\nelplo.csv'
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: 17.7s

# ── 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.