1"""Post-solution analysis: valuation ratios, Campbell–Shiller decomposition along deterministic paths,2wealth accounting, welfare (consumption-equivalent variations) and summary statistics by year."""3from __future__ import annotations45import numpy as np67KINDS = ("N", "Q", "E", "H")8910def valuation_ratios(P):11 """Per-period valuation objects for each capital type and for land/housing."""12 out = {}13 T = len(P["Y"]) - 114 for k in KINDS:15 V, D, R, K = P[f"V_{k}"], P[f"D_{k}"], P[f"R_{k}"], P[f"K_{k}"]16 pk = P["pQ"] if k == "Q" else np.ones(T + 1)17 out[f"q_over_p_{k}"] = P[f"q_{k}"] / pk18 out[f"PE_{k}"] = V / np.where(R * K[:T + 1] > 0, R * K[:T + 1], np.nan) # price-earnings19 out[f"PD_{k}"] = V / np.where(D > 0, D, np.nan) # price-payout (if positive)20 out[f"V_over_GDP_{k}"] = V / P["GDP"]21 out["land_PR"] = P["PT"] / P["RTp"] # production land price / rent22 house_value = P["q_H"] * P["K_H"][1:T + 2] + P["PTH"] * 1.023 house_rent = P["RH"] * P["HS"]24 out["house_PR"] = house_value / house_rent25 out["house_land_share"] = P["PTH"] / house_value26 out["W_over_GDP"] = P["W"] / P["GDP"]27 out["wealth_comp"] = {k: P[f"V_{k}"] / P["W"] for k in KINDS}28 out["wealth_comp"]["land"] = (P["PT"] + P["PTH"]) / P["W"]29 return out303132def campbell_shiller_decomposition(P1, P0, claim="land", horizon=None):33 """Exact log price-dividend decomposition along deterministic paths.3435 For a claim with price V_t and cash flow D_{t+1}: V_t = sum_{j>=1} D_{t+j} / prod_{i=1..j}(1+r_{t+i}).36 We compare scenario 1 to scenario 0 at t = 0 and attribute log(V1/V0) to the cash-flow path and the37 discount-rate path by the counterfactual: cash-flow effect = log V(D1, r0) - log V(D0, r0);38 discount-rate effect = log V(D1, r1) - log V(D1, r0). The two add up exactly to log(V1/V0) when both39 valuations use the same terminal convention (Gordon value at the horizon)."""40 def pv(D, r, g_term):41 T = min(len(D), len(r)) - 142 disc = np.cumprod(1.0 / (1.0 + r[1:T + 1]))43 rT = r[T]44 term = D[T] * (1 + g_term) / max(rT - g_term, 1e-4) * disc[-1]45 return float(np.sum(D[1:T + 1] * disc) + term)4647 def cash_flow(P):48 if claim == "land":49 return P["RTp"]50 if claim == "housing_land":51 return P["RTHp"]52 if claim == "house":53 return P["RH"] * P["HS"]54 return P[f"D_{claim}"] # capital type: net payout5556 D1, D0 = cash_flow(P1), cash_flow(P0)57 r1, r0 = P1["r"], P0["r"]58 g1 = D1[-1] / D1[-2] - 159 g0 = D0[-1] / D0[-2] - 160 V11, V00 = pv(D1, r1, g1), pv(D0, r0, g0)61 V10 = pv(D1, r0, g1)62 return dict(total=np.log(V11 / V00), cash_flow=np.log(V10 / V00), discount_rate=np.log(V11 / V10),63 V1=V11, V0=V00)646566def welfare_cev(x_path1, P_path1, x_path0, P_path0, beta, gamma, g=0.0, horizon=40):67 """Consumption-equivalent variation of path 1 relative to path 0 for a household with CRRA over x/P over68 the first `horizon` years (default 40: the remaining working life of a 2026 adult; with hand-to-mouth69 workers whose consumption goes to zero, infinite-horizon CEVs are dominated by the far future).70 Returns lambda such that scaling path 0 by (1+lambda) gives the utility of path 1."""71 T = min(len(x_path1), len(x_path0), horizon)72 b = beta ** np.arange(T)73 def U(x, Pp):74 c = x[:T] / Pp[:T]75 if abs(gamma - 1) < 1e-10:76 return float(np.sum(b * np.log(c)))77 return float(np.sum(b * c ** (1 - gamma) / (1 - gamma)))78 U1, U0 = U(x_path1, P_path1), U(x_path0, P_path0)79 if abs(gamma - 1) < 1e-10:80 return np.exp((U1 - U0) / np.sum(b)) - 181 return (U1 / U0) ** (1 / (1 - gamma)) - 1828384def summary_by_year(P, years=(2026, 2030, 2040, 2050, 2075, 2100)):85 """Dictionary year -> key statistics (used for tables)."""86 yrs = P["years"]87 T = len(P["Y"]) - 188 vr = valuation_ratios(P)89 rows = {}90 gf = P["gf"]91 for y in years:92 t = int(np.where(yrs == y)[0][0])93 comp_share = P["R_Q"][t] * P["K_Q"][t] / P["GDP"][t]94 rows[y] = dict(95 GDP=P["GDP"][t] / P["GDP"][0], GDP_growth=(P["GDP"][min(t + 1, T)] / P["GDP"][t] - 1),96 labor_share=P["labor_share"][t], compute_income_share=comp_share,97 structures_share=P["R_N"][t] * P["K_N"][t] / P["GDP"][t],98 land_share=(P["RTp"][t] + P["RTHp"][t]) / P["GDP"][t],99 energy_share=(P["PE"][t] * P["E"][t]) / P["GDP"][t],100 wL=P["wL"][t] * gf[t] / (P["wL"][0] * gf[0]), wH=P["wH"][t] * gf[t] / (P["wH"][0] * gf[0]),101 r=P["r"][t + 1], A=P["A"][t], KQ=P["K_Q"][t], VQ_over_GDP=vr["V_over_GDP_Q"][t],102 W_over_GDP=vr["W_over_GDP"][t], land_PR=vr["land_PR"][t], house_PR=vr["house_PR"][t],103 PE_Q=vr["PE_Q"][t], PE_N=vr["PE_N"][t], q_N=vr["q_over_p_N"][t], q_Q=vr["q_over_p_Q"][t],104 util=P["util"][t], PE=P["PE"][t] / P["PE"][0], energy_to_compute=P["EQ"][t] / P["E"][t],105 inv_share=P["inv_goods"][t] / P["GDP"][t], CW_over_GDP=P["C_W"][t] / P["GDP"][t],106 CK_over_GDP=P["C_K"][t] / P["GDP"][t], workers_cons_level=P["C_W"][t] / P["C_W"][0],107 tasks_automated=int(np.sum(P["Qi"][t] > 1e-6 * max(P["Qtot"][t], 1e-12))),108 )109 return rows110111112def arrival_lottery_premia(P_ai, P_no, pi, beta, gamma, theta_H):113 """One-period lottery illustration of AGI-arrival risk (labelled as an approximation in the paper).114115 At t=0 the economy learns with probability pi that the AI path (P_ai) is realized and with 1-pi the frozen116 path (P_no). Both paths share the 2026 stocks. The capitalists' expenditure x_0 is taken from the frozen path117 (status quo), the branch SDFs are M_b = beta (x_1^b/x_0)^{-gamma} (P_1^b/P_0)^{gamma-1}, and each claim's118 2026 price is sum_b pi_b M_b (payout_1^b + V_1^b). Returns risk-free rate, expected returns and premia by claim.119 The exercise ignores the feedback of uncertainty on 2026 investment and consumption."""120 x0 = P_no["X_K"][0]121 P0 = P_no["P_idx"][0]122 branches = [(pi, P_ai), (1 - pi, P_no)]123 M = [b * beta * (P["X_K"][1] / x0) ** (-gamma) * (P["P_idx"][1] / P0) ** (gamma - 1) for b, P in branches]124 Rf = 1.0 / sum(M)125 out = dict(Rf=Rf - 1.0)126 claims = {k: (lambda P, k=k: (P[f"D_{k}"][1] + P[f"V_{k}"][1])) for k in KINDS}127 claims["land"] = lambda P: P["RTp"][1] + P["PT"][1]128 claims["housing_land"] = lambda P: P["RTHp"][1] + P["PTH"][1]129 for name, f in claims.items():130 pay = [f(P) for _, P in branches]131 V0 = sum(m * p for m, p in zip(M, pay))132 ER = sum(b * p for (b, _), p in zip(branches, pay)) / V0133 out[name] = dict(V0=V0, expected_return=ER - 1.0, premium=ER - Rf,134 payoff_ai=pay[0], payoff_no=pay[1])135 return out136