Elections and economy: Hibbs bread-and-peace model

Source: ElectionsEconomy/hibbs.Rmd

This ports the core linear regression and predictive simulation to Python.

Load data

Code 
from pathlib import Path
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.formula.api as smf
from scipy.stats import norm

root = Path("../../ROS-Examples")
hibbs = pd.read_table(root / "ElectionsEconomy/data/hibbs.dat", sep=r"\s+")
hibbs.head()
year growth vote inc_party_candidate other_candidate
0 1952 2.40 44.60 Stevenson Eisenhower
1 1956 2.89 57.76 Eisenhower Stevenson
2 1960 0.85 49.91 Nixon Kennedy
3 1964 4.21 61.34 Johnson Goldwater
4 1968 3.02 49.60 Humphrey Nixon

Scatter plot

Code 
fig, ax = plt.subplots()
ax.scatter(hibbs["growth"], hibbs["vote"], color="black")
for _, r in hibbs.iterrows():
    ax.text(r["growth"], r["vote"], str(int(r["year"])), fontsize=8)
ax.axhline(50, color="gray", linewidth=0.8)
ax.set_xlabel("Average recent growth in personal income")
ax.set_ylabel("Incumbent party's vote share")
Text(0, 0.5, "Incumbent party's vote share")

Linear regression

R original:

M1 <- stan_glm(vote ~ growth, data = hibbs)

Python classical fit:

Code 
M1 = smf.ols("vote ~ growth", data=hibbs).fit()
M1.summary()
OLS Regression Results
Dep. Variable: vote R-squared: 0.580
Model: OLS Adj. R-squared: 0.550
Method: Least Squares F-statistic: 19.32
Date: Sun, 31 May 2026 Prob (F-statistic): 0.000610
Time: 23:01:57 Log-Likelihood: -42.839
No. Observations: 16 AIC: 89.68
Df Residuals: 14 BIC: 91.22
Df Model: 1
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept 46.2476 1.622 28.514 0.000 42.769 49.726
growth 3.0605 0.696 4.396 0.001 1.567 4.554
Omnibus: 5.392 Durbin-Watson: 2.379
Prob(Omnibus): 0.067 Jarque-Bera (JB): 2.828
Skew: -0.961 Prob(JB): 0.243
Kurtosis: 3.738 Cond. No. 4.54


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
Code 
fig, ax = plt.subplots()
ax.scatter(hibbs["growth"], hibbs["vote"], color="black")
xs = np.linspace(-0.7, 4.5, 200)
ax.plot(xs, M1.params["Intercept"] + M1.params["growth"] * xs, color="black")
ax.axhline(50, color="gray", linewidth=0.8)
ax.set_xlabel("Average recent growth in personal income")
ax.set_ylabel("Incumbent party's vote share")
Text(0, 0.5, "Incumbent party's vote share")

Posterior-ish uncertainty from asymptotic normal approximation

This reproduces the book’s informal Bayesian regression simulation logic.

Code 
rng = np.random.default_rng(123)
beta_draws = rng.multivariate_normal(M1.params.values, M1.cov_params().values, size=1000)
Code 
fig, ax = plt.subplots()
ax.scatter(hibbs["growth"], hibbs["vote"], color="black", zorder=2)
for b0, b1 in beta_draws[:50]:
    ax.plot(xs, b0 + b1 * xs, color="gray", alpha=0.25, linewidth=0.8)
ax.plot(xs, M1.params["Intercept"] + M1.params["growth"] * xs, color="black")
ax.set_xlabel("growth")
ax.set_ylabel("vote")
Text(0, 0.5, 'vote')

Prediction at 2% growth

Code 
x_new = 2.0
pred_mean = M1.params["Intercept"] + M1.params["growth"] * x_new
pred_sd = np.sqrt(M1.mse_resid)
print(pred_mean, pred_sd)
52.36870412557468 3.763287642229796
Code 
grid = np.linspace(35, 70, 500)
plt.plot(grid, norm.pdf(grid, pred_mean, pred_sd))
plt.fill_between(grid[grid >= 50], norm.pdf(grid[grid >= 50], pred_mean, pred_sd), alpha=0.3)
plt.xlabel("Incumbent party two-party vote share")
plt.title("Predictive distribution at 2% growth")
Text(0.5, 1.0, 'Predictive distribution at 2% growth')

CmdStanPy and BlackJAX notes

  • CmdStanPy gives the closest replacement for rstanarm::stan_glm if exact posterior draws are needed.
  • BlackJAX can sample the same Gaussian regression posterior from an explicit log density, but that adds plumbing without much pedagogical gain for this simple example.