1"""F5-conf-04: VFI vs closed form, optimal growth with log utility.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Reads only data/model.csv; writes out/results.json. Deterministic (no randomness)."""4import json, os5import numpy as np6import pandas as pd78ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))9df = pd.read_csv(os.path.join(ROOT, "data", "model.csv"))10p = dict(zip(df["parameter"], df["value"]))11alpha, beta = float(p["alpha"]), float(p["beta"])12y_min, y_max, N = float(p["y_min"]), float(p["y_max"]), int(float(p["grid_points"]))13spacing, tol = str(p["grid_spacing"]).strip(), float(p["vfi_tolerance"])14eval_y = [float(p["eval_y_1"]), float(p["eval_y_2"]), float(p["eval_y_3"])]15assert spacing == "log" and 0 < alpha < 1 and 0 < beta < 116data_issues = []1718# ---------- closed form (guess V = A + B ln y, verified in report) ----------19ab = alpha * beta20B = 1.0 / (1.0 - ab)21A = (np.log(1.0 - ab) + ab * B * np.log(ab)) / (1.0 - beta)22c_closed = lambda y: (1.0 - ab) * np.asarray(y)23V_closed = lambda y: A + B * np.log(np.asarray(y))24savings_rate = ab # (y - c)/y2526# ---------- grid ----------27grid = np.exp(np.linspace(np.log(y_min), np.log(y_max), N))28lng = np.log(grid)2930# ---------- VFI: choose next-period output y' on the grid ----------31# k = y'^(1/alpha) is the investment needed; c = y - k must be > 0.32k_needed = grid ** (1.0 / alpha) # indexed by j (y')33C = grid[:, None] - k_needed[None, :] # rows i: current y, cols j: next y'34U = np.where(C > 0, np.log(np.where(C > 0, C, 1.0)), -np.inf)35V = np.zeros(N)36n_iter, sup = 0, np.inf37while sup >= tol:38 Vn = np.max(U + beta * V[None, :], axis=1)39 sup = float(np.max(np.abs(Vn - V)))40 V = Vn41 n_iter += 142 if n_iter > 100000:43 raise RuntimeError("VFI did not converge")44jstar = np.argmax(U + beta * V[None, :], axis=1)45c_pol = grid - k_needed[jstar]46binding_low = int(np.sum(jstar == 0)); binding_high = int(np.sum(jstar == N - 1))47if binding_low or binding_high:48 data_issues.append({"issue": f"policy hits grid bound at {binding_low} (low) / {binding_high} (high) nodes",49 "handling": "reported; grid kept as prescribed"})5051def interp(vals, y): # linear interpolation in ln y52 return float(np.interp(np.log(y), lng, vals))5354on_node = [bool(np.min(np.abs(grid - y)) < 1e-12) for y in eval_y]55c_num = [interp(c_pol, y) for y in eval_y]56c_cl = [float(c_closed(y)) for y in eval_y]57errs = [abs(a - b) for a, b in zip(c_num, c_cl)]5859# diagnostics (not required keys)60yprime_cf = (ab * grid) ** alpha61diag = {62 "eval_points_are_grid_nodes": on_node,63 "abs_error_policy_each_eval_point": errs,64 "max_abs_error_policy_all_grid_nodes": float(np.max(np.abs(c_pol - c_closed(grid)))),65 "max_abs_error_value_all_grid_nodes": float(np.max(np.abs(V - V_closed(grid)))),66 "vfi_truncation_bound_beta_over_1mbeta_times_last_change": beta / (1 - beta) * sup,67 "log_grid_step": float(lng[1] - lng[0]),68 "closed_form_yprime_range_on_grid": [float(yprime_cf.min()), float(yprime_cf.max())],69 "policy_nodes_at_lower_bound": binding_low, "policy_nodes_at_upper_bound": binding_high,70 "closed_form_A": float(A), "closed_form_B": float(B),71 "steady_state_y_closed_form": float(ab ** (alpha / (1 - alpha))),72 "interpolation": "linear in ln y between grid nodes",73 "vfi_choice_variable": "next-period output y' restricted to grid nodes (discrete VFI), initial V=0",74}7576# refinement diagnostic (sensitivity only; headline results use the prescribed N)77def vfi_on(Nr):78 g = np.exp(np.linspace(np.log(y_min), np.log(y_max), Nr)); kn = g ** (1.0 / alpha)79 Cr = g[:, None] - kn[None, :]80 Ur = np.where(Cr > 0, np.log(np.where(Cr > 0, Cr, 1.0)), -np.inf)81 Vr, s = np.zeros(Nr), np.inf82 while s >= tol:83 Vq = np.max(Ur + beta * Vr[None, :], axis=1); s = np.max(np.abs(Vq - Vr)); Vr = Vq84 cr = g - kn[np.argmax(Ur + beta * Vr[None, :], axis=1)]85 return max(abs(float(np.interp(np.log(y), np.log(g), cr)) - float(c_closed(y))) for y in eval_y)86diag["sensitivity_max_abs_error_policy_by_grid_points"] = {str(n): vfi_on(n) for n in (200, 800)}8788res = {89 "task_id": "F5-conf-04",90 "c_closed_y05": c_cl[0], "c_closed_y1": c_cl[1], "c_closed_y2": c_cl[2],91 "c_num_y05": c_num[0], "c_num_y1": c_num[1], "c_num_y2": c_num[2],92 "max_abs_error_policy": float(max(errs)),93 "value_closed_y1": float(V_closed(1.0)),94 "value_num_y1": interp(V, 1.0),95 "savings_rate_closed_form": float(savings_rate),96 "baseline_savings_rate_closed_form": float(savings_rate),97 "contraction_modulus": float(beta),98 "grid_points_used": int(N),99 "n_iterations": int(n_iter),100 "final_sup_change": float(sup),101 "feasibility": "complete",102 "feasibility_reason": None,103 "data_issues": data_issues,104 "n_input_rows": {"model.csv": int(len(df))},105 "n_analysed": int(N),106 "random_seed": None,107 "diagnostics": diag,108}109with open(os.path.join(ROOT, "out", "results.json"), "w") as f:110 json.dump(res, f, indent=2)111print(json.dumps({k: v for k, v in res.items() if k != "diagnostics"}, indent=2))112