Code for 3 Part Series -- The Summer Slide Economy
Executable Code for series
This is the Python executable code for all simulations in the series starting today. The actual note follows in a few minutes, but this allows me to create a reference link.
Appendix: Python Code
Can be executed in Google AiStudio: https://aistudio.google.com/prompts/new_chat?model=gemini-3-pro-preview. Make sure you have “Code Execution” selected on the lower right portion of your screen.
Phase 1: The “Talent” Extinction Event
import numpy as np
import matplotlib.pyplot as plt
def run_phase_1_final_three_panels_v2():
# Parameters
n_agents = 5000
T = 40
dt = 0.1
steps = int(T / dt)
W0 = 1000.0
W_barrier = 200.0 # Ruin Barrier
# Groups: (Name, Mu, Sigma, Color, Style)
# Using distinct colors/styles to ensure all 4 are visible
groups = [
("Normie / High SES", 0.05, 0.10, "tab:blue", "-"), # Safe / Safe
("Normie / Low SES", 0.05, 0.20, "tab:cyan", "--"), # Safe / Risky
("Gifted / High SES", 0.08, 0.20, "tab:green", "-"), # Risky / Safe
("Gifted / Low SES", 0.08, 0.40, "tab:red", "--") # Risky / Risky
]
# Store data
results_avg = {}
results_median = {}
results_survival = {}
np.random.seed(42)
for name, mu, sigma, color, style in groups:
wealth = np.full(n_agents, W0)
alive = np.ones(n_agents, dtype=bool)
avg_hist = [W0]
median_hist = [W0]
survival_hist = [1.0]
for t in range(steps):
if np.any(alive):
z = np.random.normal(0, 1, n_agents)
growth = (mu - 0.5 * sigma**2) * dt + sigma * np.sqrt(dt) * z
# Update wealth for alive agents
wealth[alive] *= np.exp(growth[alive])
# Check Ruin
ruined = (wealth < W_barrier) & alive
alive[ruined] = False
wealth[ruined] = 0.0
avg_hist.append(np.mean(wealth))
median_hist.append(np.median(wealth))
survival_hist.append(np.mean(alive))
results_avg[name] = avg_hist
results_median[name] = median_hist
results_survival[name] = survival_hist
# Plotting
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(18, 6))
time = np.linspace(0, T, steps + 1)
# Common function to plot with distinct styles
def plot_lines(ax, data_dict):
for name, data in data_dict.items():
# Retrieve style info
grp = next(g for g in groups if g[0] == name)
color = grp[3]
style = grp[4]
# Plot
ax.plot(time, data, label=name, color=color, linestyle=style, linewidth=2.5, alpha=0.9)
# Panel 1: Average
plot_lines(ax1, results_avg)
ax1.set_title('Average Wealth (The Illusion)', fontweight='bold')
ax1.set_ylabel('Wealth ($)')
ax1.grid(True, alpha=0.3)
ax1.legend(loc='upper left', fontsize='small')
# Panel 2: Median
plot_lines(ax2, results_median)
ax2.set_title('Median Wealth (The Reality)', fontweight='bold')
ax2.grid(True, alpha=0.3)
ax2.legend(loc='upper left', fontsize='small')
# Panel 3: Survival
# Multiply by 100 for percentage
results_survival_pct = {k: np.array(v)*100 for k, v in results_survival.items()}
plot_lines(ax3, results_survival_pct)
ax3.set_title('Survival Rate (The Mechanism)', fontweight='bold')
ax3.set_ylabel('% Surviving')
ax3.set_ylim(0, 105)
ax3.grid(True, alpha=0.3)
ax3.legend(loc='lower left', fontsize='small')
plt.tight_layout()
plt.show()
run_phase_1_final_three_panels_v2()Phase 2: The Summer Slide Economy
import numpy as np
import matplotlib.pyplot as plt
def run_phase_2_random_shocks():
# Parameters
n_agents = 5000
T = 20 # Years
dt = 1/12 # Monthly steps
steps = int(T / dt)
W0 = 1000.0
W_barrier = 200.0
# Common Baseline (Stable Times)
mu_stable = 0.08
sigma_stable = 0.25
# Shock Parameters
shock_prob = 1/12 # Probability per month
# High SES (Buffered) Shock
mu_shock_H = 0.00
sigma_shock_H = 0.25
# Low SES (Unbuffered) Shock
mu_shock_L = -0.20
sigma_shock_L = 0.45
# Setup Groups
groups = [
("High SES (Buffered)", "green", "-"),
("Low SES (Unbuffered)", "red", "--")
]
# Data Storage
results_avg = {g[0]: [W0] for g in groups}
results_median = {g[0]: [W0] for g in groups}
results_survival = {g[0]: [1.0] for g in groups}
wealths = {
groups[0][0]: np.full(n_agents, W0),
groups[1][0]: np.full(n_agents, W0)
}
alives = {
groups[0][0]: np.ones(n_agents, dtype=bool),
groups[1][0]: np.ones(n_agents, dtype=bool)
}
np.random.seed(42)
for t in range(steps):
for name, color, style in groups:
w = wealths[name]
alive = alives[name]
if np.any(alive):
# Determine who gets shocked (Individual Shocks)
is_shocked = np.random.rand(n_agents) < shock_prob
# Parameters for each agent
# Default to stable
mus = np.full(n_agents, mu_stable)
sigmas = np.full(n_agents, sigma_stable)
# Apply shock params
if name == "High SES (Buffered)":
mus[is_shocked] = mu_shock_H
sigmas[is_shocked] = sigma_shock_H
else:
mus[is_shocked] = mu_shock_L
sigmas[is_shocked] = sigma_shock_L
# Geometric Brownian Motion Step
z = np.random.normal(0, 1, n_agents)
growth = (mus - 0.5 * sigmas**2) * dt + sigmas * np.sqrt(dt) * z
# Update Wealth
w[alive] *= np.exp(growth[alive])
# Check Ruin
ruined = (w < W_barrier) & alive
alive[ruined] = False
w[ruined] = 0.0
results_avg[name].append(np.mean(w))
results_median[name].append(np.median(w))
results_survival[name].append(np.mean(alive))
# Plotting
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(18, 6))
time = np.linspace(0, T, steps + 1)
def plot_panel(ax, data_dict, title, ylabel=None, ylim=None):
for name, color, style in groups:
y_vals = np.array(data_dict[name]) * 100 if "Survival" in title else data_dict[name]
ax.plot(time, y_vals, label=name, color=color, linestyle=style, linewidth=2.5)
ax.set_title(title, fontweight='bold')
if ylabel: ax.set_ylabel(ylabel)
if ylim: ax.set_ylim(ylim)
ax.grid(True, alpha=0.3)
plot_panel(ax1, results_avg, 'Average Wealth', ylabel='Wealth ($)')
plot_panel(ax2, results_median, 'Median Wealth')
plot_panel(ax3, results_survival, 'Survival Rate', ylabel='% Surviving', ylim=(0, 105))
# Legend
lines, labels = ax1.get_legend_handles_labels()
fig.legend(lines, labels, loc='lower center', ncol=2, bbox_to_anchor=(0.5, -0.05))
plt.tight_layout()
plt.subplots_adjust(bottom=0.15)
plt.show()
run_phase_2_random_shocks()Phase 3: Eating
import numpy as np
import matplotlib.pyplot as plt
def run_phase_3_4pct_unemployment():
# Parameters
n_agents = 5000
T = 100 # 100 Years
dt = 1/12 # Monthly steps
steps = int(T / dt)
# Calibration for 4% Unemployment Rate:
# If duration (D) = 3 months, and prob (p) is monthly:
# Steady State Rate (u) = p*D / (1 + p*D - p)
# 0.04 = 3p / (1 + 2p) -> 0.04 + 0.08p = 3p -> 0.04 = 2.92p -> p ≈ 0.0137
shock_prob = 0.0137
unemployment_duration = 3 # months
# Endowments (Latent Heat)
w_high = np.full(n_agents, 500.0)
w_low = np.full(n_agents, 50.0)
# Economic Physics
salary_yr = 55.0
col_adult_yr = 50.0
job_penalty = 10.0 # Restart penalty
mu_cap = 0.08 # 8% annual return on savings
sigma_cap = 0.25 # 25% annual volatility
results = {
"High SES": {"avg": [], "median": [], "survival": []},
"Low SES": {"avg": [], "median": [], "survival": []}
}
np.random.seed(42)
for name, wealth in [("High SES", w_high), ("Low SES", w_low)]:
w = wealth.copy()
alive = np.ones(n_agents, dtype=bool)
unemp_timer = np.zeros(n_agents)
for t in range(steps):
if np.any(alive):
# 1. Labor Market Dynamics
# Check for new shocks
new_shock = (np.random.rand(n_agents) < shock_prob) & (unemp_timer == 0)
unemp_timer[new_shock] = unemployment_duration
# Monthly Salary
monthly_salary = np.full(n_agents, salary_yr * dt)
monthly_salary[unemp_timer > 0] = 0.0
# Apply Job Penalty at re-employment (end of month 3)
just_hired = (unemp_timer == 1)
w[just_hired & alive] -= job_penalty
# Decrement timer
unemp_timer[unemp_timer > 0] -= 1
# 2. Capital Growth (Multiplicative)
z = np.random.normal(0, 1, n_agents)
growth_mask = (w > 0) & alive
growth = (mu_cap - 0.5 * sigma_cap**2) * dt + sigma_cap * np.sqrt(dt) * z
w[growth_mask] *= np.exp(growth[growth_mask])
# 3. Fixed Costs (Eating)
current_col = col_adult_yr * dt
w[alive] += monthly_salary[alive] - current_col
# 4. Survival Check
ruined = (w <= 0) & alive
alive[ruined] = False
w[ruined] = 0.0
results[name]["avg"].append(np.mean(w))
results[name]["median"].append(np.median(w))
results[name]["survival"].append(np.mean(alive) * 100)
# Visualization
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(20, 6))
time = np.linspace(0, T, steps)
for name, color, style in [("High SES", "green", "-"), ("Low SES", "red", "--")]:
ax1.plot(time, results[name]["avg"], label=name, color=color, linestyle=style, linewidth=2)
ax2.plot(time, results[name]["median"], label=name, color=color, linestyle=style, linewidth=2)
ax3.plot(time, results[name]["survival"], label=name, color=color, linestyle=style, linewidth=2)
ax1.set_title('Average Wealth (The Illusion)', fontweight='bold')
ax2.set_title('Median Wealth (The Struggle)', fontweight='bold')
ax3.set_title('Survival Rate (The Attrition)', fontweight='bold')
for ax in [ax1, ax2, ax3]:
ax.grid(True, alpha=0.3)
ax.legend()
ax3.set_ylim(0, 105)
plt.tight_layout()
plt.show()
run_phase_3_4pct_unemployment()Phase 4: The Sandwich Generation
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
def run_phase_4_silver_aristocracy():
# ---------------------------------------------------------
# 1. Global Parameters (Phase 4 Logic)
# ---------------------------------------------------------
n_agents = 5000
T = 100 # 100-year Post-New Deal Era
dt = 1/12 # Monthly steps
steps = int(T / dt)
np.random.seed(42)
# Endowments (The Inheritance Gap)
# High SES: $500 (Latent Heat / Buffer)
# Low SES: $0 (The Sub-Critical Floor)
w = np.concatenate([np.full(n_agents//2, 500.0), np.full(n_agents//2, 0.0)])
is_high_ses = np.concatenate([np.ones(n_agents//2, dtype=bool), np.zeros(n_agents//2, dtype=bool)])
# Demographics
age = np.random.randint(0, 85, n_agents).astype(float)
unemp_timer = np.zeros(n_agents)
# Child Tracking (Max 3 Children per household)
n_children = np.zeros(n_agents, dtype=int)
child_timers = [np.zeros(n_agents) for _ in range(3)]
# Economic Constants
salary_yr = 55.0 # Labor Income
col_adult_yr = 50.0 # Baseline "Eating" Cost
child_cost_yr = 25.0 # Fixed 50% "Child Tax"
job_penalty = 10.0 # Restart Haircut (delta)
parent_threshold = 375.0 # Required wealth (3yrs CoL)
# Market Risk/Return
mu_cap = 0.08 # 8% "Gifted" Potential
sigma_cap = 0.25 # 25% "Ferrari" Volatility
shock_prob = 0.0137 # Calibrated for 4% steady-state unemployment
unemployment_duration = 3
# Storage for time-series visualization
results_median = []
results_avg = []
# ---------------------------------------------------------
# 2. Main Simulation Loop
# ---------------------------------------------------------
for t in range(steps):
age += dt
# A. Death & Dynastic Replacement (Inheritance)
dead = (age >= 85)
if np.any(dead):
# The next generation begins - High SES inherit, Low SES reset to 0
w[dead & is_high_ses] = 500.0
w[dead & ~is_high_ses] = 0.0
age[dead] = 0
n_children[dead] = 0
for i in range(3): child_timers[i][dead] = 0
unemp_timer[dead] = 0
# B. Labor Market (The Effort Bucket)
is_worker = (age >= 20) & (age < 70)
is_retiree = (age >= 70) & (age < 85)
# Unemployment shocks
new_shock = (np.random.rand(n_agents) < shock_prob) & (unemp_timer == 0) & is_worker
unemp_timer[new_shock] = unemployment_duration
salary_mo = np.zeros(n_agents)
salary_mo[is_worker & (unemp_timer == 0)] = salary_yr * dt
# Apply Job Penalty (Restart Cost) at the moment of re-employment
just_hired = (unemp_timer == 1)
w[just_hired] -= job_penalty
unemp_timer[unemp_timer > 0] -= 1
# C. The Sandwich Transfer (10% Wage tax from workers to retirees)
total_transfer_pool = np.sum(salary_mo) * 0.10
salary_mo[is_worker] *= 0.90
num_retirees = np.sum(is_retiree)
if num_retirees > 0:
w[is_retiree] += (total_transfer_pool / num_retirees)
# D. Capital Growth (The Luck Bucket)
z = np.random.normal(0, 1, n_agents)
growth = (mu_cap - 0.5 * sigma_cap**2) * dt + sigma_cap * np.sqrt(dt) * z
w[w > 0] *= np.exp(growth[w > 0])
# E. The Parent Trap (Reproduction Mechanic)
can_have_kid = (age >= 25) & (age < 45) & (w >= parent_threshold) & (n_children < 3)
# Random birth chance if threshold is met
new_kid = can_have_kid & (np.random.rand(n_agents) < (1/48))
for i in np.where(new_kid)[0]:
idx = n_children[i]
child_timers[idx][i] = 18.0
n_children[i] += 1
# F. Fixed Drains (Eating + Children)
active_child_costs = np.zeros(n_agents)
for i in range(3):
is_active = child_timers[i] > 0
active_child_costs[is_active] += (child_cost_yr * dt)
child_timers[i][is_active] -= dt
w += salary_mo - (col_adult_yr * dt) - active_child_costs
# G. The Ruin Barrier (Volatility Trap)
w[w < 0.01] = 0.01
results_median.append(np.median(w))
results_avg.append(np.mean(w))
# ---------------------------------------------------------
# 3. Visualization
# ---------------------------------------------------------
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 7))
time = np.linspace(0, T, steps)
# Figure 1: Average vs Median
ax1.plot(time, results_avg, label='Average Wealth (Ensemble)', color='green', linewidth=2)
ax1.plot(time, results_median, label='Median Wealth (Typical Agent)', color='red', linestyle='--', linewidth=2)
ax1.set_title('Figure 1: 100-Year Wealth Divergence (Post-New Deal Era)', fontweight='bold')
ax1.set_xlabel('Years')
ax1.set_ylabel('Wealth ($)')
ax1.legend()
ax1.grid(True, alpha=0.3)
# Figure 2: Log-Quintile Distribution
retiree_w = w[(age >= 70) & (age < 85)]
parent_w = w[(age >= 25) & (age < 45)]
def get_quintile_means(data):
if len(data) == 0: return [0.01]*5
data_sorted = np.sort(data)
chunks = np.array_split(data_sorted, 5)
return [np.mean(c) for c in chunks]
p_quint = get_quintile_means(parent_w)
r_quint = get_quintile_means(retiree_w)
x = np.arange(5)
width = 0.35
ax2.bar(x - width/2, p_quint, width, label='Parents (Ages 25-45)', color='navy', alpha=0.8)
ax2.bar(x + width/2, r_quint, width, label='Retirees (Ages 70-85)', color='darkorange', alpha=0.8)
ax2.set_title('Figure 2: Wealth Concentration by Quintile (Log Scale)', fontweight='bold')
ax2.set_xticks(x)
ax2.set_xticklabels(['Q1', 'Q2', 'Q3', 'Q4', 'Q5'])
ax2.set_yscale('log')
ax2.set_ylabel('Wealth ($)')
ax2.legend()
ax2.grid(True, which="both", ls="-", alpha=0.2)
plt.tight_layout()
plt.show()
# Execution
run_phase_4_silver_aristocracy()Phase 5: Tax Policy
import numpy as np
import matplotlib.pyplot as plt
def run_comparative_tax_simulation():
# ---------------------------------------------------------
# 1. Global Parameters
# ---------------------------------------------------------
n_agents = 5000
T = 100 # 100 Years
dt = 1/12 # Monthly time steps
steps = int(T / dt)
# Common Economic Parameters
salary_yr = 55.0
col_adult_yr = 50.0
child_cost_yr = 25.0
job_penalty = 10.0
parent_threshold = 375.0
mu_cap = 0.08 # 8% Return (The "Ferrari")
sigma_cap = 0.25 # 25% Volatility
shock_prob = 0.0137 # 4% Unemployment Rate
unemployment_duration = 3
# The Two Scenarios
scenarios = [
{
"name": "Scenario A: Annual Tax (30% Cap / 0% Estate)",
"labor_tax": 0.10,
"annual_cap_tax": 0.30,
"estate_tax": 0.0
},
{
"name": "Scenario B: Deferred Tax (0% Cap / 50% Estate)",
"labor_tax": 0.10,
"annual_cap_tax": 0.0,
"estate_tax": 0.50
}
]
results = {}
print(f"Running Simulation for {T} Years with N={n_agents} Agents...")
for sc in scenarios:
# Reset Agents
# High SES start at $500; Low SES start at $0
w = np.concatenate([np.full(n_agents//2, 500.0), np.full(n_agents//2, 0.0)])
is_high_ses = np.concatenate([np.ones(n_agents//2, dtype=bool), np.zeros(n_agents//2, dtype=bool)])
# Consistent demographics reset
np.random.seed(42)
age = np.random.randint(0, 85, n_agents).astype(float)
unemp_timer = np.zeros(n_agents)
n_children = np.zeros(n_agents, dtype=int)
child_timers = [np.zeros(n_agents) for _ in range(3)]
median_history = []
aggregate_history = []
# Re-seed for loop dynamics
np.random.seed(42)
for t in range(steps):
age += dt
# --- DEATH & ESTATE TAX (The Harvest) ---
dead = (age >= 85)
estate_tax_revenue = 0.0
if np.any(dead):
dead_wealth = w[dead].copy()
# Apply Estate Tax
tax_collected = dead_wealth * sc["estate_tax"]
estate_tax_revenue = np.sum(tax_collected)
# Inheritance (Heir gets the remainder)
w[dead] -= tax_collected
# Reset Demographics for Heir
age[dead] = 0
n_children[dead] = 0
for i in range(3): child_timers[i][dead] = 0
unemp_timer[dead] = 0
# --- LABOR MARKET ---
is_worker = (age >= 20) & (age < 70)
new_shock = (np.random.rand(n_agents) < shock_prob) & (unemp_timer == 0) & is_worker
unemp_timer[new_shock] = unemployment_duration
salary_mo = np.zeros(n_agents)
salary_mo[is_worker & (unemp_timer == 0)] = salary_yr * dt
just_hired = (unemp_timer == 1)
w[just_hired] -= job_penalty
unemp_timer[unemp_timer > 0] -= 1
# --- LABOR TAX ---
tax_labor_rev = np.sum(salary_mo * sc["labor_tax"])
salary_mo *= (1 - sc["labor_tax"])
# --- CAPITAL GROWTH & ANNUAL TAX ---
current_w = w.copy()
z = np.random.normal(0, 1, n_agents)
growth = (mu_cap - 0.5 * sigma_cap**2) * dt + sigma_cap * np.sqrt(dt) * z
w[w > 0] *= np.exp(growth[w > 0])
# Annual Capital Tax (Only for Scenario A)
annual_tax_rev = 0.0
if sc["annual_cap_tax"] > 0:
gains = w - current_w
taxable_gains = np.maximum(gains, 0)
tax_collected = taxable_gains * sc["annual_cap_tax"]
annual_tax_rev = np.sum(tax_collected)
w -= tax_collected
# --- DIVIDEND (The Mechanism) ---
# Pool all three sources: Labor Tax + Annual Cap Tax + Estate Tax
total_revenue = tax_labor_rev + annual_tax_rev + estate_tax_revenue
dividend = total_revenue / n_agents
w += dividend
# --- COSTS ---
active_child_costs = np.zeros(n_agents)
for i in range(3):
is_active = child_timers[i] > 0
active_child_costs[is_active] += (child_cost_yr * dt)
child_timers[i][is_active] -= dt
w += salary_mo - (col_adult_yr * dt) - active_child_costs
# Floor (Absorbing Barrier)
w[w < 0.01] = 0.01
# Reproduction
can_have_kid = (age >= 25) & (age < 45) & (w >= parent_threshold) & (n_children < 3)
new_kid = can_have_kid & (np.random.rand(n_agents) < (1/48))
for i in np.where(new_kid)[0]:
idx = n_children[i]
child_timers[idx][i] = 18.0
n_children[i] += 1
median_history.append(np.median(w))
aggregate_history.append(np.sum(w))
results[sc["name"]] = {"median": median_history, "aggregate": aggregate_history}
# ---------------------------------------------------------
# 2. Visualization
# ---------------------------------------------------------
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6))
time = np.linspace(0, T, steps)
# Plot Aggregate Wealth
for name, data in results.items():
color = 'red' if "Annual" in name else 'green'
style = '--' if "Annual" in name else '-'
ax1.plot(time, data["aggregate"], label=name, color=color, linestyle=style, linewidth=2)
ax1.set_title('Figure 1: Aggregate Wealth (The Pie)', fontweight='bold')
ax1.set_ylabel('Total Wealth ($) - Log Scale')
ax1.set_xlabel('Years')
ax1.set_yscale('log') # Log scale needed to see the magnitude difference
ax1.legend()
ax1.grid(True, which="both", ls="-", alpha=0.2)
# Plot Median Wealth
for name, data in results.items():
color = 'red' if "Annual" in name else 'green'
style = '--' if "Annual" in name else '-'
ax2.plot(time, data["median"], label=name, color=color, linestyle=style, linewidth=2)
ax2.set_title('Figure 2: Median Wealth (The Household)', fontweight='bold')
ax2.set_ylabel('Median Wealth ($)')
ax2.set_xlabel('Years')
ax2.legend()
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Print Final Stats
print("\nFinal Statistics (Year 100):")
for name, data in results.items():
print(f"{name}:")
print(f" Aggregate Wealth: ${data['aggregate'][-1]:,.0f}")
print(f" Median Wealth: ${data['median'][-1]:,.0f}")
print("-" * 30)
run_comparative_tax_simulation Keep reading with a 7-day free trial
Subscribe to Yes, I give a fig... thoughts on markets from Michael Green to keep reading this post and get 7 days of free access to the full post archives.