#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Luck vs talent vs effort over a working life ============================================ An agent-based Monte Carlo model in which each simulated person's wealth at the end of a 40-year working life depends on * TALENT - an innate, unchosen ability (standard-normal score z_T), * EFFORT - (a) working, which yields a steady stream of earnings, part of which is saved, and (b) a persistent effort score z_E that raises earnings and the number of opportunities a person goes after, * LUCK - random positive and negative events: earnings shocks, job displacement, health shocks, graduating into a recession, business/investment opportunities whose payoff is a lottery, windfalls, and large wealth losses, * BACKGROUND - unchosen starting conditions (family, schooling, networks; z_B). Evidence source: ``research_report.pdf`` in this repository ("Luck, Skill and Effort in Business and Career Success: What the Evidence Shows"). Every parameter below is tagged with one of: [EVIDENCE] value taken directly from a figure quoted in the report [CALIBRATED] value solved for inside this script so that the model reproduces a figure quoted in the report [ASSUMPTION] the report has no usable evidence; the value is our choice, with the justification written next to it, and it is varied in the sensitivity analysis Three model versions are run (all with common random numbers, so the same simulated people experience the same luck in A and B): A "with work income" - the full model. B "no work income" - identical people and identical random events, but nobody earns or saves a wage; wealth changes only through returns and random events. C "TvL reference" - a re-implementation of the rules of the Talent-vs-Luck toy model described in the report (lucky event -> capital doubles with probability = talent; unlucky event -> capital halves). Used as a benchmark only. Usage ----- python simulation/model.py # full run: 1,000 runs x 1,000 agents python simulation/model.py --quick # smoke test with 40 runs (minutes -> seconds) Outputs go to ``results/`` (results.json and charts/*.png). The script prints progress and timing. Requires numpy and matplotlib only. """ from __future__ import annotations import argparse import json import math import os import sys import time from copy import deepcopy import numpy as np # ============================================================================= # 1. PARAMETERS (every model parameter lives in this section) # ============================================================================= # Size of the experiment ------------------------------------------------------ N_AGENTS = 1000 # agents per simulated population (TvL also used 1,000) N_RUNS = 1000 # independent Monte Carlo runs per model version / variant MASTER_SEED = 20261112 # all randomness flows from this seed (fully reproducible) # Each tuple: (name, baseline value, kind, source / justification) PARAMETER_TABLE = [ # ---------------------------------------------------------------- horizon ("n_years", 40, "ASSUMPTION", "Working life of 40 annual steps, ages 23-62. Start age 23 is the age at which " "Huggett, Ventura & Yaron (2011) measure initial conditions (report s.3.4); 40 years " "is also the Talent-vs-Luck horizon (report s.3.1)."), # ----------------------------------------------------------------- traits ("talent_effort_corr", 0.0, "ASSUMPTION", "Correlation between talent and effort scores. The report gives no estimate; 0 keeps " "the two channels separable. 0.3 is tested."), # ------------------------------------------- persistent earnings component ("hvy_target_share", 0.615, "EVIDENCE", "Share of the variance of (log) lifetime earnings due to conditions fixed at age 23: " "61.5% in the Huggett, Ventura & Yaron (2011) benchmark (report s.3.4 and table s.4; " "range 53.8-66%; Keane & Wolpin give ~90%)."), ("sigma_alpha", None, "CALIBRATED", "SD of the persistent person-specific log-earnings level. Solved by simulation so that " "the model's initial-conditions share of lifetime-earnings variance equals hvy_target_share."), ("w_talent", 1 / 3, "ASSUMPTION", "Share of the persistent earnings component's variance coming from talent. The report " "cannot split 'initial conditions' into talent, effort and background, so the three get " "equal shares. Todd & Zhang (report s.3.4: family background <10%, plus cognitive " "ability 14%) suggest talent+background could be smaller; 0.10-0.60 is tested."), ("w_effort", 1 / 3, "ASSUMPTION", "Share from effort (see w_talent)."), ("w_background", 1 / 3, "ASSUMPTION", "Share from unchosen background (see w_talent)."), ("talent_growth", 0.004, "EVIDENCE", "Extra annual log-earnings growth per SD of talent (learning ability). HVY (report s.3.4): " "+1 SD learning ability raises expected lifetime wealth by ~8%; 0.004/yr x ~20 yrs average " "exposure = ~8%. Tested 0 and 0.008."), ("age_slope", 0.05, "ASSUMPTION", "Hump-shaped age-earnings profile (linear term). Not in the report; a conventional shape " "that roughly doubles median earnings between 23 and the early 50s."), ("age_curvature", -0.0009, "ASSUMPTION", "Quadratic term of the age profile (peak about age 51)."), # ------------------------------------------------ year-to-year earnings luck ("rho_persistent", 0.97, "ASSUMPTION", "Annual persistence of the persistent earnings shock. The report says many shocks are " "persistent (Guvenen et al. 2021) but gives no number; 0.97 is a near-unit-root value."), ("shock_large_prob", 0.10, "CALIBRATED", "Normal-mixture earnings shock (report s.4 guidance: 'mostly small changes, with rare large " "ones (a normal mixture fits this)'). Probability of a draw from the large component. The " "five mixture/transitory values are set jointly so that annual log-earnings changes have " "kurtosis inside the report's 10-30 band and strong negative skew (Guvenen et al. 2021)."), ("shock_small_sd", 0.04, "CALIBRATED", "SD of the small (common) mixture component."), ("shock_large_mean", -0.15, "CALIBRATED", "Mean of the large mixture component (negative -> negative skew). The small component's " "mean is set so the overall shock has mean zero."), ("shock_large_sd", 0.40, "CALIBRATED", "SD of the large mixture component."), ("sigma_transitory", 0.04, "CALIBRATED", "SD of a purely transitory (one-year) log-earnings shock."), ("guvenen_asymmetry", True, "EVIDENCE", "Guvenen et al. 2021 (report s.3.3): large positive shocks to high earners are transitory " "while negative ones persist; the reverse for low earners. Implemented for large shocks " "using earnings terciles (the tercile cut-offs are our choice). Switched off in sensitivity."), ("earnings_sd_target", 0.60, "ASSUMPTION", "Cross-sectional SD of log annual earnings, averaged over ages 23-62. The report pins the SHAPE " "of earnings shocks (kurtosis, skew) and the variance SHARE due to initial conditions (HVY), but " "not the overall spread of earnings. 0.6 is our approximate reading of typical values for " "working-age U.S. earners (published estimates span roughly 0.5-0.9 depending on sample); it is " "NOT from the report. Tested 0.45 and 0.80."), ("earnings_shock_scale", None, "CALIBRATED", "Multiplier on the five earnings-shock values above (mixture SDs/means and transitory SD). Solved " "jointly with sigma_alpha so that the model hits BOTH earnings_sd_target and hvy_target_share. " "Discrete events (displacement, health, recession entry) are not scaled: their sizes are evidence-based."), ("shock_distribution", "mixture", "ASSUMPTION", "'mixture' (baseline) or 'gaussian' (same variance, no fat tails) - sensitivity only."), # ------------------------------------------------- discrete negative events ("displacement_rate", 0.02, "ASSUMPTION", "Annual probability of a mass-layoff job loss. The report flags this as a gap ('no good " "direct estimate'). 2%/yr means ~55% of people are displaced at least once in 40 years. " "Tested 1% and 4%."), ("displacement_pv_loss_low", 1.4, "EVIDENCE", "Present-value earnings loss from displacement when unemployment <6%: 1.4 years of " "pre-displacement earnings (20-yr horizon, 5% discount; Davis & von Wachter 2011, report s.3.5)."), ("displacement_pv_loss_high", 2.8, "EVIDENCE", "Same when unemployment >8%: 2.8 years (Davis & von Wachter 2011)."), ("displacement_high_unemp_share", 0.25, "ASSUMPTION", "Share of displacements that happen in high-unemployment years. Not in the report. Tested 0 and 1."), ("displacement_decay", 0.90, "ASSUMPTION", "Annual persistence of the displacement earnings scar. The year-0 loss is then solved so the " "20-year PV loss matches 1.4 / 2.8 years exactly (gives ~20% / ~39% initial losses)."), ("health_rate_start", 0.003, "ASSUMPTION", "Annual probability of a major health shock (hospital admission) at age 23; rises linearly " "to health_rate_end at 62. The report has no frequency estimate (gap). About 35% of people " "are hit at least once. Scaled x0.5 and x2 in sensitivity."), ("health_rate_end", 0.020, "ASSUMPTION", "Annual probability at age 62 (see health_rate_start)."), ("health_rate_scale", 1.0, "ASSUMPTION", "Multiplier on both health rates (sensitivity knob)."), ("health_earnings_loss", 0.20, "EVIDENCE", "A hospital admission cuts annual earnings by ~20% (Dobkin et al. 2018, report s.3.5)."), ("health_decay", 1.0, "EVIDENCE", "Persistence of the health earnings loss. The report says the decline 'is persistent', so " "baseline = no recovery. A recovering version (0.90/yr) is tested."), ("entry_recession_prob", 0.20, "ASSUMPTION", "Chance of entering the labour market in a recession. Not in the report; roughly one " "entry year in five."), ("entry_recession_loss", 0.10, "ASSUMPTION", "Initial earnings loss for recession graduates. The report gives the duration (~10 years) " "but not the size. Tested 0 and 0.20."), ("entry_recession_years", 10, "EVIDENCE", "The loss fades linearly to zero over 10 years: 'persistent earnings declines lasting ten " "years' (Oreopoulos, von Wachter & Heisz 2012, report s.3.5)."), ("entry_recession_disadvantage", 0.30, "ASSUMPTION", "Losses are larger for the less advantaged (direction: report s.3.5). Size is our choice: " "loss x exp(-0.3 z_B)."), # ------------------------------------------------------- saving and wealth ("work_income", True, "ASSUMPTION", "Version switch: True = people work and save part of their earnings (version A); False = no " "earnings enter wealth (version B)."), ("savings_rate", 0.10, "ASSUMPTION", "Share of each year's earnings added to wealth. Not in the report. Tested 0.05 and 0.20."), ("savings_income_elasticity", 0.0, "ASSUMPTION", "0 = everyone saves the same share of earnings. >0 = higher earners save a larger share " "(rate proportional to (earnings/typical earnings)^elasticity, rescaled so the aggregate saving rate " "is unchanged, capped at 50%). Not in the report; 1.0 is tested."), ("initial_wealth", 1.0, "ASSUMPTION", "Starting wealth at 23, in units of one year of median entry-level earnings (all money in " "the model uses this unit). Equal for everyone at baseline so 'worse off than at the start' " "is well defined. Tested 0.25 and 4."), ("initial_wealth_sigma", 0.0, "ASSUMPTION", "Log-SD of unequal starting wealth, tied to background z_B (HVY include initial wealth in " "'initial conditions'). 0 at baseline; 1.0 tested."), ("return_rate", 0.03, "ASSUMPTION", "Real annual return on wealth. Not in the report. Tested 0 and 0.05."), ("le_discount", 0.03, "ASSUMPTION", "Discount rate for the present value of lifetime earnings (used only for the HVY variance " "decomposition)."), # --------------------------------------- random positive events: opportunities ("opp_rate", 0.10, "ASSUMPTION", "Annual probability that an average-effort person meets a business/investment opportunity " "(~4 per career). The report: 'no good direct evidence' on the frequency of major positive " "events (gap); Elvidge (2020) finds the lucky-event probability is the most influential TvL " "input. Tested 0.05 and 0.20."), ("effort_opp_elasticity", 0.30, "ASSUMPTION", "Effort raises the number of 'shots' (report s.5: 'maximise the number of independent " "shots'): rate x exp(0.3 z_E - 0.045). No magnitude in the report. Tested 0 and 0.6."), ("stake_fraction", 0.25, "ASSUMPTION", "Share of current wealth put into each opportunity. Not in the report. Tested 0.10 and 0.50."), ("venture_base_success", 0.20, "EVIDENCE", "Average success probability of a venture: 18% (NBER) / 21% (JFE) for first-time founders " "(Gompers et al.); ~25% of VC-backed founders get something at exit (Hall & Woodward). " "Tested 0.18, 0.21, 0.25."), ("venture_skill_ratio", 1.6, "EVIDENCE", "Success probability of a +1 SD talent person relative to an average one (logit model). " "Report s.4 guidance: skill shifts success odds by ~1.5-1.7x for strong vs average " "('strong' = +1 SD is our reading). Tested 1.0 (no skill), 1.5, 1.7 and the ~2.1 implied " "by Gompers' raw persistence."), ("venture_repeat_success", 0.30, "EVIDENCE", "Success rate of founders whose previous venture succeeded (Gompers et al., report s.3.2). " "Used for validation and for the 'persistence-implied' skill variant."), ("venture_share_over_5x", 0.06, "EVIDENCE", "~6% of VC-backed ventures return more than 5x (Kerr, Nanda & Rhodes-Kropf 2014, report s.3.3)."), ("venture_top_return_share", 0.50, "EVIDENCE", "Those >5x ventures supply roughly half of gross returns (Kerr et al. 2014). Together with " "the two figures above this pins down a Pareto payoff for successes. Tested 0.4 and 0.6."), # failure -> the stake is lost (Hall & Woodward: ~3/4 of founders get nothing at exit) # ---------------------------------------- random positive events: windfalls ("windfall_rate", 0.01, "ASSUMPTION", "Annual probability of a pure-luck windfall (inheritance, gift, prize). Not in the report. " "Tested 0.005 and 0.02."), ("windfall_median", 1.0, "ASSUMPTION", "Median windfall size (years of median entry earnings)."), ("windfall_sigma", 1.0, "ASSUMPTION", "Log-SD of windfall size (lognormal)."), ("windfall_labour_response", 0.01, "EVIDENCE", "Windfalls reduce annual labour earnings by ~1% of the prize, persistently (Cesarini et al. " "2017, report s.3.5)."), # ---------------------------------------------- random negative wealth events ("loss_rate", 0.02, "ASSUMPTION", "Annual probability of a large wealth loss (failed asset, divorce, uninsured loss, fraud). " "Not in the report. Tested 0.01 and 0.05."), ("loss_fraction", 0.50, "EVIDENCE", "Share of wealth lost: the TvL unlucky event halves capital (report s.3.1). TvL is " "uncalibrated, so 0.25 and 0.75 are tested."), # ---------------------------------------------- C: Talent-vs-Luck reference ("tvl_talent_mean", 0.6, "EVIDENCE", "TvL talent ~ Normal(0.6, 0.1) truncated to [0,1]. The report says talent is normal but " "does not quote the moments; these are the original paper's values (Pluchino et al. 2018)."), ("tvl_talent_sd", 0.1, "EVIDENCE", "See tvl_talent_mean."), ("tvl_initial_capital", 10.0, "EVIDENCE", "Everyone starts with capital 10 (original TvL paper)."), ("tvl_steps", 80, "EVIDENCE", "40 years in six-month steps (report s.3.1)."), ("tvl_event_prob", 0.05, "ASSUMPTION", "Per-step probability of meeting a lucky event, and separately an unlucky one. TvL " "generates encounters from moving points on a grid; the report gives no rate. 0.05 gives " "~4 of each per career. Tested 0.02 and 0.10."), ] BASE = {name: value for name, value, _kind, _why in PARAMETER_TABLE} PARAM_DOCS = {name: {"baseline": value, "kind": kind, "source": why} for name, value, kind, why in PARAMETER_TABLE} # Diagnostics only: person-years with earnings below this (in units of median entry earnings) are # excluded from the earnings-change statistics, mirroring the minimum-earnings threshold used in # administrative-data studies such as Guvenen et al. (2021). Does not affect the simulation. MIN_EARNINGS_DIAG = 0.05 # Discount rate and horizon used by Davis & von Wachter for their PV figures. DVW_DISCOUNT, DVW_HORIZON = 0.05, 20 # ============================================================================= # 2. CALIBRATION HELPERS # ============================================================================= _GH_X, _GH_W = np.polynomial.hermite_e.hermegauss(120) # quadrature for E[f(z)], z~N(0,1) _GH_W = _GH_W / _GH_W.sum() def _sigmoid(v): return 1.0 / (1.0 + np.exp(-v)) def _solve_logit_intercept(b, p_bar): """Intercept a such that E[sigmoid(a + b z)] = p_bar for z ~ N(0,1).""" lo, hi = -30.0, 30.0 for _ in range(200): a = 0.5 * (lo + hi) if (_GH_W * _sigmoid(a + b * _GH_X)).sum() < p_bar: lo = a else: hi = a return 0.5 * (lo + hi) def venture_logit_from_ratio(p_bar, ratio): """(a, b) with mean success p_bar and sigmoid(a+b)/sigmoid(a) = ratio.""" if ratio <= 1.0 + 1e-12: return _solve_logit_intercept(0.0, p_bar), 0.0 lo, hi = 0.0, 6.0 for _ in range(100): b = 0.5 * (lo + hi) a = _solve_logit_intercept(b, p_bar) if _sigmoid(a + b) / _sigmoid(a) < ratio: lo = b else: hi = b b = 0.5 * (lo + hi) return _solve_logit_intercept(b, p_bar), b def venture_logit_from_persistence(p_bar, p_repeat): """(a, b) such that a pure latent-skill model reproduces P(success | previous success).""" lo, hi = 0.0, 6.0 for _ in range(100): b = 0.5 * (lo + hi) a = _solve_logit_intercept(b, p_bar) p = _sigmoid(a + b * _GH_X) if (_GH_W * p * p).sum() / (_GH_W * p).sum() < p_repeat: lo = b else: hi = b b = 0.5 * (lo + hi) return _solve_logit_intercept(b, p_bar), b def venture_logit_diagnostics(a, b): """Persistence statistics implied by a logit skill model (exact, by quadrature).""" p = _sigmoid(a + b * _GH_X) ep, ep2 = (_GH_W * p).sum(), (_GH_W * p * p).sum() return { "mean_success": ep, "success_after_success": ep2 / ep, "success_after_failure": (ep - ep2) / (1 - ep), "consecutive_outcome_correlation": (ep2 - ep * ep) / (ep * (1 - ep)), "ratio_plus1sd_vs_average": _sigmoid(a + b) / _sigmoid(a), } def pareto_from_kerr(p_bar, share_over_5x, top_return_share): """Pareto (x_min, alpha) for the payoff multiple of successful ventures. Failures pay 0. Successes (probability p_bar) pay M ~ Pareto(x_min, alpha). Two Kerr et al. facts pin both values: P(M > 5) * p_bar = share_over_5x -> (x_min/5)^alpha = q share of all returns from M > 5 = top_share -> (5/x_min)^(1-alpha) = h Dividing gives 5/x_min = h/q. """ q = share_over_5x / p_bar h = top_return_share if not (0 < q < h < 1): raise ValueError("Kerr-calibrated Pareto needs share_over_5x/p_bar < top_return_share < 1") k = h / q x_min = 5.0 / k alpha = math.log(1.0 / q) / math.log(k) return x_min, alpha def displacement_log_scar(pv_target, decay, discount=DVW_DISCOUNT, horizon=DVW_HORIZON): """Year-0 log earnings scar d (<0) such that the PV of losses equals pv_target years. The scar decays as d * decay**k, so the loss in year k is 1 - exp(d * decay**k). """ k = np.arange(horizon) lo, hi = 1e-9, 1 - 1e-12 for _ in range(200): L = 0.5 * (lo + hi) pv = ((1 - np.exp(np.log1p(-L) * decay ** k)) / (1 + discount) ** k).sum() if pv < pv_target: lo = L else: hi = L L = 0.5 * (lo + hi) return math.log1p(-L), L def derive(P, fixed_earnings=None, calib_seed=None, verbose=False): """Compute derived / calibrated quantities for parameter set P. fixed_earnings = (earnings_shock_scale, sigma_alpha) reuses earlier calibrated values; otherwise both are solved by simulation (see calibrate_earnings). """ D = {} a, b = venture_logit_from_ratio(P["venture_base_success"], P["venture_skill_ratio"]) D["venture_a"], D["venture_b"] = a, b D["venture_diag"] = venture_logit_diagnostics(a, b) D["pareto_xmin"], D["pareto_alpha"] = pareto_from_kerr( P["venture_base_success"], P["venture_share_over_5x"], P["venture_top_return_share"]) D["disp_logscar_low"], D["disp_initial_loss_low"] = displacement_log_scar( P["displacement_pv_loss_low"], P["displacement_decay"]) D["disp_logscar_high"], D["disp_initial_loss_high"] = displacement_log_scar( P["displacement_pv_loss_high"], P["displacement_decay"]) # normal mixture: small-component mean chosen so that the shock has mean zero pi = P["shock_large_prob"] D["shock_small_mean"] = -pi * P["shock_large_mean"] / (1 - pi) if fixed_earnings is None and P["sigma_alpha"] is not None and P["earnings_shock_scale"] is not None: fixed_earnings = (P["earnings_shock_scale"], P["sigma_alpha"]) if fixed_earnings is None: k, sa, info = calibrate_earnings(P, D, seed=calib_seed, verbose=verbose) D["earnings_calibration_check"] = info else: k, sa = fixed_earnings D["earnings_shock_scale"], D["sigma_alpha"] = k, sa return D # ============================================================================= # 3. THE SIMULATION # ============================================================================= def simulate(P, D, n_runs, n_agents, trait_seed, shock_seed, diagnostics=False, track_age=False, earnings_only=False, track_earnings_sd=False): """Simulate n_runs independent populations of n_agents people over P['n_years'] years. All arrays have shape (n_runs, n_agents). Random numbers are always drawn in the same order whatever the switches, so two calls with the same seeds give the same people the same luck ("common random numbers"); that makes version comparisons paired. Traits come from ``trait_seed``, life events from ``shock_seed``: re-running with a new ``shock_seed`` gives the same people a new set of luck ("twin lives"). """ R, N, T = int(n_runs), int(n_agents), int(P["n_years"]) shp = (R, N) f32 = np.float32 # ---------------- traits: fixed at age 23 ------------------------------- rt = np.random.default_rng(trait_seed) zT = rt.standard_normal(shp) # talent c = P["talent_effort_corr"] zE = c * zT + math.sqrt(1 - c * c) * rt.standard_normal(shp) # effort zB = rt.standard_normal(shp) # unchosen background rs = np.random.default_rng(shock_seed) entry_recession = rs.random(shp) < P["entry_recession_prob"] # persistent individual earnings level (talent + effort + background) w = np.array([P["w_talent"], P["w_effort"], P["w_background"]], float) w = w / w.sum() alpha = D["sigma_alpha"] * (np.sqrt(w[0]) * zT + np.sqrt(w[1]) * zE + np.sqrt(w[2]) * zB) entry_loss = np.clip(P["entry_recession_loss"] * np.exp(-P["entry_recession_disadvantage"] * zB), 0.0, 0.9) s0 = P["initial_wealth_sigma"] W0 = P["initial_wealth"] * np.exp(s0 * zB - 0.5 * s0 * s0) W = W0.copy() g = P["effort_opp_elasticity"] opp_rate_i = np.clip(P["opp_rate"] * np.exp(g * zE - 0.5 * g * g), 0.0, 1.0) p_success_i = _sigmoid(D["venture_a"] + D["venture_b"] * zT) x_min, a_par = D["pareto_xmin"], D["pareto_alpha"] # state variables p = np.zeros(shp) # persistent earnings shock d_disp = np.zeros(shp) # displacement log-scar d_health = np.zeros(shp) # health log-scar cum_windfall = np.zeros(shp) LE = np.zeros(shp) # PV of lifetime earnings cnt = {k: np.zeros(shp, np.int16) for k in ("opportunities", "successes", "windfalls", "losses", "displacements", "health_shocks")} k_sc = D["earnings_shock_scale"] earn_sd = [] if track_earnings_sd else None pi, mu_L, sd_L = P["shock_large_prob"], k_sc * P["shock_large_mean"], k_sc * P["shock_large_sd"] mu_S, sd_S = k_sc * D["shock_small_mean"], k_sc * P["shock_small_sd"] sd_trans = k_sc * P["sigma_transitory"] gaussian = P["shock_distribution"] == "gaussian" if gaussian: # same mean (0) and variance as the mixture sd_gauss = math.sqrt((1 - pi) * (sd_S ** 2 + mu_S ** 2) + pi * (sd_L ** 2 + mu_L ** 2)) rho, rho_d, rho_h = P["rho_persistent"], P["displacement_decay"], P["health_decay"] log_health = math.log1p(-P["health_earnings_loss"]) diag_n = min(R, 20) # runs kept for earnings-change diagnostics logY_hist = [] if diagnostics else None venture_draws = [] if diagnostics else None age_track = {"mean": [], "p10": [], "p50": [], "p90": [], "p99": []} if track_age else None def _record_age(Wnow): q = np.percentile(Wnow, [10, 50, 90, 99]) age_track["mean"].append(float(Wnow.mean())) for key, v in zip(("p10", "p50", "p90", "p99"), q): age_track[key].append(float(v)) if track_age: _record_age(W) for t in range(T): # ---------- random draws (fixed order: common random numbers) ---------- u_mix = rs.random(shp, dtype=f32) n_shock = rs.standard_normal(shp, dtype=f32) n_trans = rs.standard_normal(shp, dtype=f32) u_disp = rs.random(shp, dtype=f32) u_sev = rs.random(shp, dtype=f32) u_health = rs.random(shp, dtype=f32) if not earnings_only: u_opp = rs.random(shp, dtype=f32) u_succ = rs.random(shp, dtype=f32) u_par = rs.random(shp) # float64: keeps the Pareto tail exact u_wf = rs.random(shp, dtype=f32) n_wf = rs.standard_normal(shp, dtype=f32) u_loss = rs.random(shp, dtype=f32) # ---------- earnings ---------------------------------------------------- if gaussian: shock = sd_gauss * n_shock.astype(float) large = np.zeros(shp, bool) else: large = u_mix < pi shock = np.where(large, mu_L + sd_L * n_shock, mu_S + sd_S * n_shock).astype(float) if P["guvenen_asymmetry"] and not gaussian: level = alpha + p q1, q2 = np.quantile(level[: min(R, 100)], [1 / 3, 2 / 3]) transitory = large & (((shock > 0) & (level > q2)) | ((shock < 0) & (level < q1))) p = rho * p + np.where(transitory, 0.0, shock) t_shock = np.where(transitory, shock, 0.0) else: p = rho * p + shock t_shock = 0.0 eps = sd_trans * n_trans displaced = u_disp < P["displacement_rate"] scar = np.where(u_sev < P["displacement_high_unemp_share"], D["disp_logscar_high"], D["disp_logscar_low"]) d_disp = rho_d * d_disp + np.where(displaced, scar, 0.0) h_rate = P["health_rate_scale"] * (P["health_rate_start"] + (P["health_rate_end"] - P["health_rate_start"]) * t / max(T - 1, 1)) sick = u_health < h_rate d_health = rho_h * d_health + np.where(sick, log_health, 0.0) fade = max(0.0, 1.0 - t / P["entry_recession_years"]) entry_scar = np.where(entry_recession, np.log1p(-entry_loss * fade), 0.0) if fade > 0 else 0.0 logY = (alpha + P["age_slope"] * t + P["age_curvature"] * t * t + P["talent_growth"] * zT * t + p + t_shock + eps + d_disp + d_health + entry_scar) Y = np.maximum(np.exp(logY) - P["windfall_labour_response"] * cum_windfall, 0.0) LE += Y / (1 + P["le_discount"]) ** t if track_earnings_sd: earn_sd.append(float(logY.std())) if diagnostics: # like Guvenen et al., drop person-years below a minimum earnings threshold Yd = Y[:diag_n] logY_hist.append(np.where(Yd >= MIN_EARNINGS_DIAG, np.log(np.maximum(Yd, 1e-300)), np.nan) .astype(np.float32)) cnt["displacements"] += displaced cnt["health_shocks"] += sick if earnings_only: continue # ---------- wealth ------------------------------------------------------ W = W * (1 + P["return_rate"]) # opportunity: stake a fraction of wealth in a lottery-like venture opp = u_opp < opp_rate_i succ = u_succ < p_success_i M = np.where(succ, x_min * (1.0 - u_par) ** (-1.0 / a_par), 0.0) f = P["stake_fraction"] W = np.where(opp, W * (1 - f + f * M), W) if diagnostics: sel = opp[:diag_n] venture_draws.append(np.stack([zT[:diag_n][sel], M[:diag_n][sel]])) # windfall: pure luck, independent of talent wf = u_wf < P["windfall_rate"] amount = np.where(wf, P["windfall_median"] * np.exp(P["windfall_sigma"] * n_wf), 0.0) W = W + amount cum_windfall += amount # large wealth loss: pure luck lost = u_loss < P["loss_rate"] W = np.where(lost, W * (1 - P["loss_fraction"]), W) # steady work income: save a share of this year's earnings if P["work_income"]: eta = P["savings_income_elasticity"] if eta == 0: W = W + P["savings_rate"] * Y else: x = Y / math.exp(P["age_slope"] * t + P["age_curvature"] * t * t) raw = x ** eta scale = P["savings_rate"] * Y.sum(1, keepdims=True) / np.maximum((Y * raw).sum(1, keepdims=True), 1e-300) W = W + np.minimum(scale * raw, 0.5) * Y cnt["opportunities"] += opp cnt["successes"] += opp & succ cnt["windfalls"] += wf cnt["losses"] += lost if track_age: _record_age(W) out = {"zT": zT, "zE": zE, "zB": zB, "W0": W0, "W": W, "LE": LE, "entry_recession": entry_recession, "counts": cnt} if diagnostics: out["logY_hist"] = np.stack(logY_hist, axis=-1) # (diag_n, N, T) out["venture_draws"] = np.concatenate(venture_draws, axis=1) if venture_draws else None if track_age: out["age_track"] = age_track if track_earnings_sd: out["earnings_sd_by_age"] = earn_sd return out def simulate_tvl(P, n_runs, n_agents, seed): """Version C: the Talent-vs-Luck rules (report s.3.1) with a fixed encounter probability.""" R, N = int(n_runs), int(n_agents) rng = np.random.default_rng(seed) T = np.clip(P["tvl_talent_mean"] + P["tvl_talent_sd"] * rng.standard_normal((R, N)), 0, 1) C = np.full((R, N), P["tvl_initial_capital"]) n_lucky = np.zeros((R, N), np.int16) n_used = np.zeros((R, N), np.int16) n_unlucky = np.zeros((R, N), np.int16) pe = P["tvl_event_prob"] for _ in range(int(P["tvl_steps"])): lucky = rng.random((R, N), dtype=np.float32) < pe used = lucky & (rng.random((R, N), dtype=np.float32) < T) unlucky = rng.random((R, N), dtype=np.float32) < pe C = np.where(used, 2 * C, C) C = np.where(unlucky, 0.5 * C, C) n_lucky += lucky n_used += used n_unlucky += unlucky zT = (T - P["tvl_talent_mean"]) / P["tvl_talent_sd"] W0 = np.full((R, N), P["tvl_initial_capital"]) return {"zT": zT, "talent_raw": T, "W0": W0, "W": C, "counts": {"lucky_events": n_lucky, "lucky_events_exploited": n_used, "unlucky_events": n_unlucky}} def calibrate_earnings(P, D, seed=None, verbose=False, n_runs=40, n_agents=1000): """Jointly solve (earnings_shock_scale, sigma_alpha) so that the model reproduces (1) hvy_target_share: the share of log lifetime-earnings variance due to conditions fixed at age 23 (Huggett, Ventura & Yaron 2011), and (2) earnings_sd_target: the average cross-sectional SD of log annual earnings (assumption). (1) is measured exactly with 'twin lives': the same people (same traits) are simulated twice with independent luck; Cov(log LE_1, log LE_2) / Var(log LE) is the share of variance explained by everything fixed at the start. Nested bisection: inner on sigma_alpha, outer on the scale. """ seed = MASTER_SEED + 7 if seed is None else seed Dl = dict(D) def measure(k, sa): Dl["earnings_shock_scale"], Dl["sigma_alpha"] = k, sa s1 = simulate(P, Dl, n_runs, n_agents, seed, seed + 1, earnings_only=True, track_earnings_sd=True) s2 = simulate(P, Dl, n_runs, n_agents, seed, seed + 2, earnings_only=True) return twin_share(np.log(s1["LE"]), np.log(s2["LE"])), float(np.mean(s1["earnings_sd_by_age"])) def solve_sigma_alpha(k): lo, hi = 0.0, 3.0 for _ in range(14): mid = 0.5 * (lo + hi) if measure(k, mid)[0] < P["hvy_target_share"]: lo = mid else: hi = mid sa = 0.5 * (lo + hi) return (sa,) + measure(k, sa) lo, hi = 0.05, 5.0 for _ in range(12): k = 0.5 * (lo + hi) if solve_sigma_alpha(k)[2] < P["earnings_sd_target"]: lo = k else: hi = k k = 0.5 * (lo + hi) sa, share, sd = solve_sigma_alpha(k) if verbose: print(f" earnings calibrated: shock scale {k:.4f}, sigma_alpha {sa:.4f} " f"(share {share:.3f} vs {P['hvy_target_share']}, SD {sd:.3f} vs {P['earnings_sd_target']})") return k, sa, {"hvy_share_achieved": share, "earnings_sd_achieved": sd} # ============================================================================= # 4. METRICS # ============================================================================= def rowcorr(x, y): """Pearson correlation computed separately for each run (row).""" xm = x - x.mean(axis=1, keepdims=True) ym = y - y.mean(axis=1, keepdims=True) return (xm * ym).sum(1) / np.sqrt((xm * xm).sum(1) * (ym * ym).sum(1)) def rank_rows(x): return np.argsort(np.argsort(x, axis=1, kind="stable"), axis=1, kind="stable") def twin_share(y1, y2): """Share of Var(y) explained by fixed traits: Cov(y1, y2) / mean variance, pooled.""" y1 = (y1 - y1.mean(1, keepdims=True)).ravel() y2 = (y2 - y2.mean(1, keepdims=True)).ravel() return float((y1 * y2).mean() / (0.5 * ((y1 * y1).mean() + (y2 * y2).mean()))) def r2_on_traits(y, traits): """R^2 of a pooled OLS of y on the given trait arrays (within-run demeaned).""" yv = (y - y.mean(1, keepdims=True)).ravel() X = np.column_stack([(t - t.mean(1, keepdims=True)).ravel() for t in traits]) beta, *_ = np.linalg.lstsq(X, yv, rcond=None) resid = yv - X @ beta return float(1 - resid.var() / yv.var()) def run_metrics(sim): """All headline metrics, one value per run (arrays of length n_runs).""" W, W0, zT = sim["W"], sim["W0"], sim["zT"] R, N = W.shape Ws = np.sort(W, axis=1) tot = Ws.sum(1) k1, k10, k20 = max(1, round(0.01 * N)), round(0.10 * N), round(0.20 * N) i = np.arange(1, N + 1) rT = rank_rows(zT) rW = rank_rows(W) order = np.argsort(-W, axis=1) top1, top2_10 = order[:, :k1], order[:, k1:k10] m = { "median": np.median(W, axis=1), "mean": W.mean(1), "p10": np.percentile(W, 10, axis=1), "p90": np.percentile(W, 90, axis=1), "p99": np.percentile(W, 99, axis=1), "max": Ws[:, -1], "top1_share": Ws[:, -k1:].sum(1) / tot, "top10_share": Ws[:, -k10:].sum(1) / tot, "top20_share": Ws[:, -k20:].sum(1) / tot, "bottom50_share": Ws[:, : N // 2].sum(1) / tot, "gini": 2 * (Ws * i).sum(1) / (N * tot) - (N + 1) / N, "share_worse_off": (W < W0).mean(1), "pearson_talent": rowcorr(zT, W), "pearson_talent_logw": rowcorr(zT, np.log(np.maximum(W, 1e-300))), "spearman_talent": rowcorr(rT.astype(float), rW.astype(float)), "richest_talent_pct": (rT[np.arange(R), order[:, 0]] + 0.5) / N, "talent_top1": np.take_along_axis(zT, top1, 1).mean(1), "talent_top2_10": np.take_along_axis(zT, top2_10, 1).mean(1), "top1_from_top_talent_decile": (np.take_along_axis(rT, top1, 1) >= 0.9 * N).mean(1), } m["denrell_liu_gap"] = m["talent_top1"] - m["talent_top2_10"] m["talent_richest"] = zT[np.arange(R), order[:, 0]] m["talent_sd_top1"] = np.take_along_axis(zT, top1, 1).std(1, ddof=1) if k1 > 1 else np.full(R, np.nan) m["talent_sd_top2_10"] = np.take_along_axis(zT, top2_10, 1).std(1, ddof=1) if "zE" in sim: rE, rB = rank_rows(sim["zE"]), rank_rows(sim["zB"]) m["spearman_effort"] = rowcorr(rE.astype(float), rW.astype(float)) m["spearman_background"] = rowcorr(rB.astype(float), rW.astype(float)) return m def summarize(x): """Mean across runs, 95% CI of that mean, and the 2.5-97.5% run-to-run range.""" x = np.asarray(x, float) x = x[np.isfinite(x)] n = len(x) mean = float(x.mean()) se = float(x.std(ddof=1) / math.sqrt(n)) if n > 1 else float("nan") return {"mean": mean, "ci95_low": mean - 1.96 * se, "ci95_high": mean + 1.96 * se, "run_p2_5": float(np.percentile(x, 2.5)), "run_p97_5": float(np.percentile(x, 97.5)), "n_runs": n} def by_decile(sim, key="zT"): """Mean / median final wealth and P(top 1%) by within-run decile of a trait.""" W, z = sim["W"], sim[key] R, N = W.shape srt = np.argsort(z, axis=1) Wz = np.take_along_axis(W, srt, 1).reshape(R, 10, N // 10) thr = np.quantile(W, 0.99, axis=1, keepdims=True) topz = np.take_along_axis(W >= thr, srt, 1).reshape(R, 10, N // 10) W0z = np.take_along_axis(sim["W0"], srt, 1).reshape(R, 10, N // 10) out = [] for d in range(10): out.append({"decile": d + 1, "mean_wealth": summarize(Wz[:, d].mean(1)), "median_wealth": summarize(np.median(Wz[:, d], axis=1)), "prob_top1pct": summarize(topz[:, d].mean(1)), "share_worse_off": summarize((Wz[:, d] < W0z[:, d]).mean(1))}) return out def talent_by_wealth(sim, n_bins=100, top_ranks=50): """Mean talent z-score by wealth percentile (100 bins) and by exact wealth rank (1 = richest).""" W, zT = sim["W"], sim["zT"] R, N = W.shape srt = np.argsort(-W, axis=1) # richest first zs = np.take_along_axis(zT, srt, 1) by_pct = zs.reshape(R, n_bins, N // n_bins).mean(2) # bin 0 = richest 1% rows_pct = [{"wealth_percentile_from_top": i + 1, **summarize(by_pct[:, i])} for i in range(n_bins)] rows_rank = [{"wealth_rank": i + 1, **summarize(zs[:, i])} for i in range(min(top_ranks, N))] return {"by_wealth_percentile": rows_pct, "by_top_rank": rows_rank, "note": "Mean talent z-score (0 = average, 1 = one SD above). Percentile 1 = richest 1% of a run."} def histogram_data(W, W0, lo=-3.0, hi=6.0, width=0.1): """Histogram of log10(final wealth / starting wealth), pooled over all agents.""" x = np.log10(np.maximum(W / W0, 1e-300)).ravel() edges = np.round(np.arange(lo, hi + width / 2, width), 6) counts, _ = np.histogram(np.clip(x, lo, hi - 1e-9), bins=edges) return {"log10_bin_edges": edges.tolist(), "counts": counts.astype(int).tolist(), "share": (counts / x.size).tolist(), "n_agents": int(x.size), "n_below_range": int((x < lo).sum()), "n_above_range": int((x >= hi).sum()), "note": "Values outside the range are clipped into the first/last bin."} def pooled_percentiles(W, W0): qs = [1, 5, 10, 25, 50, 75, 90, 95, 99, 99.9] return {"percentiles": qs, "final_wealth": np.percentile(W, qs).tolist(), "final_over_initial": np.percentile(W / W0, qs).tolist()} def earnings_change_stats(logY_hist): """Moments of annual log-earnings changes (validation vs Guvenen et al. 2021).""" share_below = float(np.isnan(logY_hist).mean()) dy = np.diff(logY_hist, axis=-1).astype(float).ravel() dy = dy[np.isfinite(dy)] m, s = dy.mean(), dy.std() z = (dy - m) / s edges = np.round(np.arange(-2.0, 2.0001, 0.05), 4) counts, _ = np.histogram(np.clip(dy, -2, 2 - 1e-9), bins=edges) return {"n": int(dy.size), "share_person_years_below_min_earnings": share_below, "min_earnings_threshold": MIN_EARNINGS_DIAG, "mean": float(m), "sd": float(s), "skewness": float((z ** 3).mean()), "kurtosis": float((z ** 4).mean()), "share_abs_change_below_0.05": float((np.abs(dy) < 0.05).mean()), "share_abs_change_above_0.5": float((np.abs(dy) > 0.5).mean()), "p1": float(np.percentile(dy, 1)), "p99": float(np.percentile(dy, 99)), "hist_edges": edges.tolist(), "hist_counts": counts.astype(int).tolist()} def venture_validation(venture_draws): zt, M = venture_draws return {"n_ventures": int(M.size), "share_zero_payoff": float((M == 0).mean()), "share_over_5x": float((M > 5).mean()), "share_of_returns_from_over_5x": float(M[M > 5].sum() / M.sum()), "mean_multiple": float(M.mean()), "corr_talent_with_multiple": float(np.corrcoef(zt, M)[0, 1])} # ============================================================================= # 5. EXPERIMENTS # ============================================================================= KEY_METRICS = ["median", "mean", "top1_share", "top10_share", "top20_share", "gini", "share_worse_off", "pearson_talent", "pearson_talent_logw", "spearman_talent", "spearman_effort", "spearman_background", "richest_talent_pct", "talent_top1", "talent_top2_10", "denrell_liu_gap", "top1_from_top_talent_decile", "talent_richest", "talent_sd_top1", "talent_sd_top2_10", "bottom50_share", "p10", "p90", "p99", "max"] def metric_summaries(m, keys=None): keys = keys or [k for k in KEY_METRICS if k in m] return {k: summarize(m[k]) for k in keys if k in m} def twin_analysis(P, D, R, N, trait_seed, shock_seed, shock_seed_2): """Same people, two independent sets of luck.""" s1 = simulate(P, D, R, N, trait_seed, shock_seed) s2 = simulate(P, D, R, N, trait_seed, shock_seed_2) W1, W2 = s1["W"], s2["W"] k1, k10 = max(1, round(0.01 * N)), round(0.10 * N) r1, r2 = rank_rows(W1), rank_rows(W2) top1_1, top1_2 = r1 >= N - k1, r2 >= N - k1 top10_1, top10_2 = r1 >= N - k10, r2 >= N - k10 below_med_2 = r2 < N // 2 res = { "log_wealth_share_from_fixed_traits": twin_share(np.log(W1), np.log(W2)), "wealth_level_share_from_fixed_traits": twin_share(W1, W2), "rank_correlation_between_lives": summarize(rowcorr(r1.astype(float), r2.astype(float))), "p_top1_again_given_top1": float((top1_1 & top1_2).sum() / top1_1.sum()), "p_top10_again_given_top10": float((top10_1 & top10_2).sum() / top10_1.sum()), "p_below_median_given_top1_in_other_life": float((top1_1 & below_med_2).sum() / top1_1.sum()), "log_lifetime_earnings_share_from_fixed_traits": twin_share(np.log(s1["LE"]), np.log(s2["LE"])), "r2_log_wealth_on_traits": r2_on_traits(np.log(W1), [s1["zT"], s1["zE"], s1["zB"]]), "r2_log_wealth_on_talent_only": r2_on_traits(np.log(W1), [s1["zT"]]), "r2_log_wealth_on_effort_only": r2_on_traits(np.log(W1), [s1["zE"]]), "r2_log_wealth_on_background_only": r2_on_traits(np.log(W1), [s1["zB"]]), "sample_pairs": {"log10_wealth_life1": np.log10(W1[0]).round(4).tolist(), "log10_wealth_life2": np.log10(W2[0]).round(4).tolist()}, } return res def sensitivity_variants(P_base, D_base): """One-at-a-time variants: (parameter, label, overrides, applies_to, recalibrate_sigma_alpha).""" a, b = venture_logit_from_persistence(P_base["venture_base_success"], P_base["venture_repeat_success"]) persist_ratio = venture_logit_diagnostics(a, b)["ratio_plus1sd_vs_average"] V = [] def add(param, values, applies="AB", recal=False, label_fmt="{}"): for v in values: ov = v if isinstance(v, dict) else {param: v} label = label_fmt.format(v if not isinstance(v, dict) else list(v.values())[0]) V.append({"parameter": param, "label": label, "overrides": ov, "applies_to": applies, "recalibrate_sigma_alpha": recal}) # evidence-based parameters across their evidence ranges add("hvy_target_share", [0.538, 0.66, 0.90], "A", recal=True) add("venture_skill_ratio", [1.0, 1.5, 1.7, round(persist_ratio, 3)]) add("venture_base_success", [0.18, 0.21, 0.25]) add("venture_top_return_share", [0.40, 0.60]) add("loss_fraction", [0.25, 0.75]) add("talent_growth", [0.0, 0.008], "A") add("displacement_high_unemp_share", [0.0, 1.0], "A") add("health_decay", [0.90], "A") add("guvenen_asymmetry", [False], "A") # assumptions add("w_talent", [{"w_talent": 0.10, "w_effort": 0.45, "w_background": 0.45}, {"w_talent": 0.60, "w_effort": 0.20, "w_background": 0.20}], "A") add("talent_effort_corr", [0.3]) add("opp_rate", [0.05, 0.20]) add("effort_opp_elasticity", [0.0, 0.6]) add("stake_fraction", [0.10, 0.50]) add("loss_rate", [0.01, 0.05]) add("windfall_rate", [0.005, 0.02]) add("return_rate", [0.0, 0.05]) add("initial_wealth", [0.25, 4.0]) add("initial_wealth_sigma", [1.0]) add("savings_rate", [0.05, 0.20], "A") add("displacement_rate", [0.01, 0.04], "A") add("health_rate_scale", [0.5, 2.0], "A") add("entry_recession_loss", [0.0, 0.20], "A") add("shock_distribution", ["gaussian"], "A") add("earnings_sd_target", [0.45, 0.80], "A", recal=True) add("savings_income_elasticity", [1.0], "A") add("rho_persistent", [0.90, 0.995], "A", recal=True) add("n_years", [30, 45]) # structural switches: turn single luck channels off add("switch:no_opportunities", [{"opp_rate": 0.0}], label_fmt="opportunities off") add("switch:no_windfalls", [{"windfall_rate": 0.0}], label_fmt="windfalls off") add("switch:no_wealth_losses", [{"loss_rate": 0.0}], label_fmt="wealth-loss events off") add("switch:no_wealth_events", [{"opp_rate": 0.0, "windfall_rate": 0.0, "loss_rate": 0.0}], "A", label_fmt="all wealth events off") add("switch:no_earnings_luck", [{"shock_large_prob": 0.0, "shock_small_sd": 0.0, "sigma_transitory": 0.0, "displacement_rate": 0.0, "health_rate_scale": 0.0, "entry_recession_prob": 0.0}], "A", label_fmt="all earnings luck off") return V SENS_METRICS = ["median", "mean", "top1_share", "top10_share", "top20_share", "gini", "share_worse_off", "pearson_talent", "pearson_talent_logw", "spearman_talent", "spearman_effort", "spearman_background", "richest_talent_pct", "denrell_liu_gap", "top1_from_top_talent_decile"] def run_sensitivity(P_base, D_base, R, N, seeds, versions=("A", "B"), log=print): """Run every one-at-a-time variant for versions A and B with common random numbers.""" out = {"A_with_income": [], "B_no_income": []} variants = sensitivity_variants(P_base, D_base) recal_cache = {} for vi, var in enumerate(variants): for ver in versions: if ver not in var["applies_to"]: continue P = deepcopy(P_base) P.update(var["overrides"]) P["work_income"] = ver == "A" t0 = time.time() if var["recalibrate_sigma_alpha"]: D = recal_cache.get(var["label"] + var["parameter"]) if D is None: D = derive(P, fixed_earnings=None, calib_seed=seeds["calibration"]) recal_cache[var["label"] + var["parameter"]] = D else: D = derive(P, fixed_earnings=(D_base["earnings_shock_scale"], D_base["sigma_alpha"])) sim = simulate(P, D, R, N, seeds["trait"], seeds["shock"]) m = run_metrics(sim) res = {"parameter": var["parameter"], "label": var["label"], "overrides": {k: (v if not isinstance(v, float) else round(v, 6)) for k, v in var["overrides"].items()}, "metrics": metric_summaries(m, SENS_METRICS)} if ver == "A": res["lifetime_earnings_r2_on_traits"] = r2_on_traits( np.log(np.maximum(sim["LE"], 1e-12)), [sim["zT"], sim["zE"], sim["zB"]]) res["log_wealth_r2_on_traits"] = r2_on_traits(np.log(np.maximum(sim["W"], 1e-300)), [sim["zT"], sim["zE"], sim["zB"]]) res["sigma_alpha_used"] = D["sigma_alpha"] res["earnings_shock_scale_used"] = D["earnings_shock_scale"] out["A_with_income" if ver == "A" else "B_no_income"].append(res) log(f" [{vi + 1}/{len(variants)}] {ver} {var['parameter']}={var['label']}: " f"gini={res['metrics']['gini']['mean']:.3f} " f"spearman={res['metrics']['spearman_talent']['mean']:.3f} " f"({time.time() - t0:.1f}s)") return out def tvl_sensitivity(P_base, R, N, seed, log=print): out = [] for pe in (0.02, 0.10): P = deepcopy(P_base) P["tvl_event_prob"] = pe m = run_metrics(simulate_tvl(P, R, N, seed)) out.append({"parameter": "tvl_event_prob", "label": str(pe), "overrides": {"tvl_event_prob": pe}, "metrics": metric_summaries(m, [k for k in SENS_METRICS if k in m]), "max_capital_run_range_over_mean": float( (np.percentile(m["max"], 97.5) - np.percentile(m["max"], 2.5)) / m["max"].mean())}) log(f" TvL event prob {pe}: gini={out[-1]['metrics']['gini']['mean']:.3f}") return out def robustness_table(sens, baseline): """For each metric: baseline, min/max over all variants, and the most influential parameter.""" table = {} for ver, rows in sens.items(): table[ver] = {} for k in SENS_METRICS: if k not in baseline[ver]: continue base = baseline[ver][k]["mean"] vals = [(r["metrics"][k]["mean"], r["parameter"], r["label"]) for r in rows if k in r["metrics"]] if not vals: continue lo = min(vals, key=lambda x: x[0]) hi = max(vals, key=lambda x: x[0]) swing = {} for v, prm, _ in vals: rng = swing.setdefault(prm, [base, base]) rng[0], rng[1] = min(rng[0], v), max(rng[1], v) ranked = sorted(swing.items(), key=lambda kv: kv[1][1] - kv[1][0], reverse=True) table[ver][k] = { "baseline": base, "min": lo[0], "min_variant": f"{lo[1]}={lo[2]}", "max": hi[0], "max_variant": f"{hi[1]}={hi[2]}", "swing_by_parameter": [{"parameter": prm, "low": r[0], "high": r[1], "swing": r[1] - r[0]} for prm, r in ranked], } return table def paired_checks(sens, base_A, base_B): """For every variant run in both A and B: does each qualitative A-vs-B conclusion still hold?""" B_rows = {(r["parameter"], r["label"]): r for r in sens["B_no_income"]} pairs = [((r["parameter"], r["label"]), r, B_rows[(r["parameter"], r["label"])]) for r in sens["A_with_income"] if (r["parameter"], r["label"]) in B_rows] pairs.append((("baseline", "baseline"), {"metrics": base_A}, {"metrics": base_B})) checks = { "gini_lower_with_income": lambda a, b: a["gini"]["mean"] < b["gini"]["mean"], "top1_share_lower_with_income": lambda a, b: a["top1_share"]["mean"] < b["top1_share"]["mean"], "spearman_talent_higher_with_income": lambda a, b: a["spearman_talent"]["mean"] > b["spearman_talent"]["mean"], "share_worse_off_lower_with_income": lambda a, b: a["share_worse_off"]["mean"] < b["share_worse_off"]["mean"], "median_higher_with_income": lambda a, b: a["median"]["mean"] > b["median"]["mean"], } res = {} for name, fn in checks.items(): ok = [fn(a["metrics"], b["metrics"]) for _, a, b in pairs] fails = [f"{k[0]}={k[1]}" for (k, a, b), o in zip(pairs, ok) if not o] res[name] = {"holds_in": int(sum(ok)), "of": len(ok), "fails_in": fails} # single-version checks for ver, rows, base in (("A_with_income", sens["A_with_income"], base_A), ("B_no_income", sens["B_no_income"], base_B)): allm = [r["metrics"] for r in rows] + [base] res[f"{ver}:spearman_talent_positive"] = { "holds_in": int(sum(m["spearman_talent"]["ci95_low"] > 0 for m in allm)), "of": len(allm)} res[f"{ver}:top1_talent_below_top2_10_talent (Denrell-Liu)"] = { "holds_in": int(sum(m["denrell_liu_gap"]["ci95_high"] < 0 for m in allm)), "of": len(allm)} res[f"{ver}:richest_is_not_top_talent_decile_on_average"] = { "holds_in": int(sum(m["richest_talent_pct"]["mean"] < 0.9 for m in allm)), "of": len(allm)} res[f"{ver}:gini_above_0.5"] = { "holds_in": int(sum(m["gini"]["mean"] > 0.5 for m in allm)), "of": len(allm)} return res # ============================================================================= # 6. OUTPUT HELPERS # ============================================================================= def _clean(o, nd=6): """Make numpy objects JSON-serialisable and round floats to nd significant digits.""" if isinstance(o, dict): return {str(k): _clean(v, nd) for k, v in o.items()} if isinstance(o, (list, tuple)): return [_clean(v, nd) for v in o] if isinstance(o, np.ndarray): return _clean(o.tolist(), nd) if isinstance(o, (np.bool_, bool)): return bool(o) if isinstance(o, (np.integer,)): return int(o) if isinstance(o, (float, np.floating)): f = float(o) if not math.isfinite(f): return None return float(f"{f:.{nd}g}") return o def agent_sample(simA, simB, D, n=1000): """Run 0 of versions A and B: the same 1,000 people with and without work income.""" N = simA["W"].shape[1] n = min(n, N) rT = rank_rows(simA["zT"][:1])[0] pS = _sigmoid(D["venture_a"] + D["venture_b"] * simA["zT"][0]) rows = [] for i in range(n): rows.append({ "id": i, "talent_z": simA["zT"][0, i], "talent_percentile": (rT[i] + 0.5) / N, "effort_z": simA["zE"][0, i], "background_z": simA["zB"][0, i], "venture_success_prob": pS[i], "initial_wealth": simA["W0"][0, i], "final_wealth_with_income": simA["W"][0, i], "final_wealth_no_income": simB["W"][0, i], "lifetime_earnings_pv": simA["LE"][0, i], "entered_in_recession": bool(simA["entry_recession"][0, i]), **{f"n_{k}": int(v[0, i]) for k, v in simA["counts"].items()}, }) return rows def tvl_sample(simC, n=1000): N = simC["W"].shape[1] rT = rank_rows(simC["zT"][:1])[0] return [{"id": i, "talent": simC["talent_raw"][0, i], "talent_percentile": (rT[i] + 0.5) / N, "initial_capital": simC["W0"][0, i], "final_capital": simC["W"][0, i], **{f"n_{k}": int(v[0, i]) for k, v in simC["counts"].items()}} for i in range(min(n, N))] SCHEMA = { "meta": "Run settings: number of runs, agents per run, seeds, runtime, software versions.", "parameters": "Every model parameter: baseline value, kind (EVIDENCE / CALIBRATED / ASSUMPTION) and source or justification.", "derived_parameters": "Values solved inside the model from report evidence (venture logit, Pareto tail, displacement scars, sigma_alpha).", "validation": "Model outputs compared with the report figures they were calibrated to or should reproduce.", "versions": ("Object keyed by version id (A_with_income, B_no_income, C_tvl_reference). Each has: " "description; metrics {name: {mean, ci95_low, ci95_high, run_p2_5, run_p97_5, n_runs}} where " "mean is the average over runs, ci95_* is the 95% confidence interval of that mean, and run_p* " "is the 2.5-97.5% range across individual runs (run-to-run variability of a 1,000-person " "population); histogram {log10_bin_edges, counts, share} of log10(final wealth / initial wealth) " "pooled over all agents in all runs; pooled_percentiles; by_talent_decile and by_effort_decile " "(list of 10 rows: decile 1 = lowest; mean_wealth, median_wealth, prob_top1pct, share_worse_off); " "by_background_decile (same, by background); talent_by_wealth {by_wealth_percentile: 100 rows " "(1 = richest 1%) and by_top_rank: 50 rows (1 = richest person), each the mean talent z-score with CIs}; " "richest_talent_pct_histogram (counts of runs by the talent percentile of the richest agent, 20 bins " "of 5 percentiles from 0 to 100); wealth_by_age (A and B: pooled mean/p10/p50/p90/p99 of wealth at " "ages 23..63); twin_lives (A and B: same people re-simulated with new luck; sample_pairs holds " "log10 wealth in both lives for the 1,000 agents of run 0)."), "metric_definitions": "Plain-language definition of each metric name.", "comparison_A_minus_B": "Paired differences (A minus B, same people and same luck) for each metric with 95% CI.", "sensitivity": ("{A_with_income: [...], B_no_income: [...], C_tvl_reference: [...]}: one row per one-at-a-time " "variant: parameter, label (value tried), overrides (all parameters changed), metrics (same " "format as versions.*.metrics)."), "robustness": ("For each version and metric: baseline value, min and max over all sensitivity variants (with the " "variant that produced them) and swing_by_parameter sorted from most to least influential."), "robustness_checks": "Qualitative conclusions and in how many sensitivity variants each one holds.", "agents_sample": ("1,000 agents from run 0. Fields: id, talent_z (standard-normal talent score), talent_percentile " "(0-1 within the population), effort_z, background_z, venture_success_prob, initial_wealth, " "final_wealth_with_income (version A), final_wealth_no_income (version B: the same person with the " "same luck but no wage), lifetime_earnings_pv, entered_in_recession, and counts of opportunities, " "successes, windfalls, losses, displacements and health_shocks. Wealth is in units of one year of " "median entry-level earnings."), "agents_sample_tvl": ("1,000 agents from run 0 of the TvL reference: talent (0-1 scale), talent_percentile, " "initial_capital (10), final_capital, counts of lucky events, exploited lucky events and unlucky events."), "units": "All money amounts are in units of one year of median earnings at age 23. TvL capital is in TvL units (start = 10).", } METRIC_DEFINITIONS = { "median": "Median final wealth across the 1,000 agents of a run.", "mean": "Mean final wealth.", "p10/p90/p99/max": "Percentiles / maximum of final wealth within a run.", "top1_share / top10_share / top20_share": "Share of total final wealth held by the richest 1% / 10% / 20%.", "bottom50_share": "Share of total wealth held by the poorest 50%.", "gini": "Gini coefficient of final wealth (0 = equal, 1 = one person has everything).", "share_worse_off": "Share of agents whose final wealth is below their own starting wealth.", "pearson_talent": "Pearson correlation between talent score and final wealth (raw levels; sensitive to outliers).", "pearson_talent_logw": "Pearson correlation between talent score and log final wealth.", "spearman_talent / spearman_effort / spearman_background": "Spearman rank correlation of each trait with final wealth.", "richest_talent_pct": "Talent percentile (0-1) of the single richest agent in the run.", "talent_top1 / talent_top2_10": "Mean talent z-score of the richest 1% / of agents ranked in the top 2-10% by wealth.", "denrell_liu_gap": "talent_top1 minus talent_top2_10 (negative = the very top is less talented than the next tier).", "top1_from_top_talent_decile": "Share of the richest 1% who are in the top 10% of talent (10% = no link).", } # ============================================================================= # 7. MAIN # ============================================================================= def main(argv=None): ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) ap.add_argument("--runs", type=int, default=N_RUNS, help="Monte Carlo runs per version/variant") ap.add_argument("--agents", type=int, default=N_AGENTS, help="agents per run") ap.add_argument("--quick", action="store_true", help="smoke test: 40 runs, sensitivity included") ap.add_argument("--no-sensitivity", action="store_true", help="skip the sensitivity analysis") ap.add_argument("--no-charts", action="store_true", help="skip chart generation") ap.add_argument("--out", default=None, help="results directory (default: /results)") args = ap.parse_args(argv) R = 40 if args.quick else args.runs N = args.agents repo = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) outdir = args.out or os.path.join(repo, "results") os.makedirs(os.path.join(outdir, "charts"), exist_ok=True) seeds = {"trait": MASTER_SEED + 1, "shock": MASTER_SEED + 2, "shock_twin": MASTER_SEED + 3, "tvl": MASTER_SEED + 4, "calibration": MASTER_SEED + 7} T0 = time.time() def log(msg): print(f"[{time.time() - T0:7.1f}s] {msg}", flush=True) log(f"Runs per version: {R}; agents per run: {N}") P = deepcopy(BASE) log("Calibrating...") D = derive(P, calib_seed=seeds["calibration"], verbose=True) log(f" venture logit a={D['venture_a']:.3f} b={D['venture_b']:.3f}; Pareto x_min={D['pareto_xmin']:.3f} " f"alpha={D['pareto_alpha']:.3f}; displacement initial loss {D['disp_initial_loss_low']:.3f}/" f"{D['disp_initial_loss_high']:.3f}") # ---------------- main versions ---------------- log("Version A (with work income)...") PA = dict(P, work_income=True) simA = simulate(PA, D, R, N, seeds["trait"], seeds["shock"], diagnostics=True, track_age=True) mA = run_metrics(simA) log("Version B (no work income)...") PB = dict(P, work_income=False) simB = simulate(PB, D, R, N, seeds["trait"], seeds["shock"], diagnostics=False, track_age=True) mB = run_metrics(simB) log("Version C (Talent-vs-Luck reference)...") simC = simulate_tvl(P, R, N, seeds["tvl"]) mC = run_metrics(simC) log("Twin lives (same people, new luck)...") twinA = twin_analysis(PA, D, R, N, seeds["trait"], seeds["shock"], seeds["shock_twin"]) twinB = twin_analysis(PB, D, R, N, seeds["trait"], seeds["shock"], seeds["shock_twin"]) # ---------------- validation against the report ---------------- ec = earnings_change_stats(simA["logY_hist"]) vv = venture_validation(simA["venture_draws"]) LE = simA["LE"] logLE = np.log(LE) validation = { "lifetime_earnings_variance_share_initial_conditions": { "target_HVY_2011": P["hvy_target_share"], "report_range": [0.538, 0.66], "model_twin_lives": twinA["log_lifetime_earnings_share_from_fixed_traits"], "model_r2_on_traits": r2_on_traits(logLE, [simA["zT"], simA["zE"], simA["zB"]]), "r2_talent_only": r2_on_traits(logLE, [simA["zT"]]), "r2_effort_only": r2_on_traits(logLE, [simA["zE"]]), "r2_background_only": r2_on_traits(logLE, [simA["zB"]])}, "annual_log_earnings_changes": { "report_target": "strongly negatively skewed; kurtosis ~10-30 (Gaussian = 3); most changes small, rare large jumps", **ec}, "venture_payoffs": { "report_targets": {"share_zero_or_loss": "~0.73-0.75 founders get nothing (Hall & Woodward); 55% of VC deals at a loss (Kerr)", "share_over_5x": 0.06, "share_of_returns_from_over_5x": 0.5, "expert_score_correlation_with_returns": 0.1}, **vv}, "venture_persistence_implied_by_baseline_skill": { "report_targets": {"first_time": "0.18 (NBER) / 0.21 (JFE)", "after_success": 0.30, "after_failure": "0.20 (NBER) / 0.22 (JFE)", "consecutive_outcome_correlation": "~0.10-0.15 (report's interpretation)"}, **D["venture_diag"]}, "venture_persistence_calibrated_to_gompers": venture_logit_diagnostics( *venture_logit_from_persistence(P["venture_base_success"], P["venture_repeat_success"])), "displacement": {"report_target_pv_years": [P["displacement_pv_loss_low"], P["displacement_pv_loss_high"]], "initial_loss_low_unemployment": D["disp_initial_loss_low"], "initial_loss_high_unemployment": D["disp_initial_loss_high"], "share_displaced_at_least_once": float((simA["counts"]["displacements"] > 0).mean()), "mean_displacements_per_career": float(simA["counts"]["displacements"].mean())}, "health": {"share_with_health_shock": float((simA["counts"]["health_shocks"] > 0).mean())}, "event_counts_per_career": {k: float(v.mean()) for k, v in simA["counts"].items()}, "share_with_no_successful_venture": float((simA["counts"]["successes"] == 0).mean()), "median_lifetime_earnings_pv": float(np.median(LE)), "lifetime_earnings_p90_over_p10": float(np.percentile(LE, 90) / np.percentile(LE, 10)), "sd_log_lifetime_earnings": float(logLE.std()), "sd_log_earnings_age23": float(np.nanstd(simA["logY_hist"][..., 0])), "sd_log_earnings_mean_over_ages": float(np.nanstd(simA["logY_hist"], axis=(0, 1)).mean()), "earnings_sd_target": P["earnings_sd_target"], "sd_log_earnings_age62": float(np.nanstd(simA["logY_hist"][..., -1])), } # ---------------- assemble version blocks ---------------- def version_block(desc, sim, m, twin=None, age=True): blk = {"description": desc, "metrics": metric_summaries(m), "histogram": histogram_data(sim["W"], sim["W0"]), "pooled_percentiles": pooled_percentiles(sim["W"], sim["W0"]), "by_talent_decile": by_decile(sim, "zT"), "talent_by_wealth": talent_by_wealth(sim)} if "zE" in sim: blk["by_effort_decile"] = by_decile(sim, "zE") blk["by_background_decile"] = by_decile(sim, "zB") if age and "age_track" in sim: blk["wealth_by_age"] = {"age": list(range(23, 23 + int(P["n_years"]) + 1)), **sim["age_track"]} if twin is not None: blk["twin_lives"] = twin blk["richest_talent_pct_histogram"] = np.histogram(m["richest_talent_pct"], bins=np.linspace(0, 1, 21))[0].tolist() return blk versions = { "A_with_income": version_block( "Full model: talent, effort, steady work income (10% saved), earnings shocks, displacement, " "health and recession-entry shocks, venture opportunities, windfalls and wealth losses.", simA, mA, twinA), "B_no_income": version_block( "Same people and same random events as A, but no work income enters wealth. Wealth changes only " "through returns, venture opportunities, windfalls and wealth losses.", simB, mB, twinB), "C_tvl_reference": version_block( "Talent-vs-Luck rules (report s.3.1): talent ~ N(0.6, 0.1); per six-month step a lucky event " "(prob 0.05) doubles capital with probability = talent, an unlucky event (prob 0.05) halves it; " "no income, no returns.", simC, mC, age=False), } log("Paired A-B comparison...") comparison = {k: summarize(mA[k] - mB[k]) for k in KEY_METRICS if k in mA and k in mB} results = { "schema": SCHEMA, "metric_definitions": METRIC_DEFINITIONS, "meta": {"n_runs": R, "n_agents_per_run": N, "n_years": P["n_years"], "seeds": seeds, "numpy": np.__version__, "python": sys.version.split()[0], "generated_by": "simulation/model.py", "quick_mode": bool(args.quick)}, "parameters": PARAM_DOCS, "derived_parameters": {k: v for k, v in D.items() if k != "venture_diag"}, "validation": validation, "versions": versions, "comparison_A_minus_B": comparison, } # ---------------- sensitivity ---------------- if not args.no_sensitivity: log("Sensitivity analysis...") sens = run_sensitivity(P, D, R, N, seeds, log=log) log("TvL sensitivity...") sens["C_tvl_reference"] = tvl_sensitivity(P, R, N, seeds["tvl"], log=log) results["sensitivity"] = sens base_metrics = {"A_with_income": versions["A_with_income"]["metrics"], "B_no_income": versions["B_no_income"]["metrics"]} results["robustness"] = robustness_table({k: sens[k] for k in ("A_with_income", "B_no_income")}, base_metrics) results["robustness_checks"] = paired_checks(sens, base_metrics["A_with_income"], base_metrics["B_no_income"]) tvl_max = mC["max"] results["validation"]["tvl_run_to_run_variability"] = { "report_claim": "Elvidge (2020): identical inputs give a range of outputs over 50% of the mean", "max_capital_mean": float(tvl_max.mean()), "max_capital_p2_5_to_p97_5_range_over_mean": float( (np.percentile(tvl_max, 97.5) - np.percentile(tvl_max, 2.5)) / tvl_max.mean())} results["agents_sample"] = agent_sample(simA, simB, D) results["agents_sample_tvl"] = tvl_sample(simC) results["meta"]["runtime_seconds"] = time.time() - T0 path = os.path.join(outdir, "results.json") with open(path, "w") as fh: json.dump(_clean(results), fh, indent=1) log(f"Wrote {path} ({os.path.getsize(path) / 1e6:.1f} MB)") if not args.no_charts: log("Charts...") sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) import charts # local module: simulation/charts.py charts.make_all(_clean(results), simA, simB, simC, twinA, twinB, os.path.join(outdir, "charts")) log("Done.") return results if __name__ == "__main__": main()