1"""F5-conf-04: VFI vs closed form for optimal growth with log utility.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py"""3import json, os4import numpy as np5import pandas as pd67df = pd.read_csv("data/model.csv")8p = dict(zip(df["parameter"], df["value"]))9alpha = float(p["alpha"]); beta = float(p["beta"])10y_min = float(p["y_min"]); y_max = float(p["y_max"]); n = int(float(p["grid_points"]))11tol = float(p["vfi_tolerance"])12evals = [float(p["eval_y_1"]), float(p["eval_y_2"]), float(p["eval_y_3"])]13assert str(p["grid_spacing"]).strip().lower() == "log"1415# ---- closed form (guess V = A + B ln y, verified in report) ----16ab = alpha * beta17B = 1.0 / (1.0 - ab)18A = (np.log(1.0 - ab) + ab / (1.0 - ab) * np.log(ab)) / (1.0 - beta)19c_closed = lambda y: (1.0 - ab) * np.asarray(y)20V_closed = lambda y: A + B * np.log(y)2122# ---- VFI on log-spaced grid, choice = next-period output y' on the same grid ----23y = np.exp(np.linspace(np.log(y_min), np.log(y_max), n))24k_needed = y ** (1.0 / alpha) # capital needed to reach y' : k = y'^(1/alpha)25C = y[:, None] - k_needed[None, :] # consumption for (y, y')26feasible = C > 027U = np.full_like(C, -np.inf)28U[feasible] = np.log(C[feasible])2930V = np.zeros(n)31it = 0; sup = np.inf32while sup >= tol:33 V_new = np.max(U + beta * V[None, :], axis=1)34 sup = float(np.max(np.abs(V_new - V)))35 V = V_new; it += 136 if it > 100000:37 raise RuntimeError("VFI did not converge")38idx = np.argmax(U + beta * V[None, :], axis=1)39c_grid = y - k_needed[idx]40at_upper = int(np.sum(idx == n - 1)); at_lower = int(np.sum(idx == 0))4142# interpolate (linear in ln y) to evaluation points43ly = np.log(y)44c_num = [float(np.interp(np.log(e), ly, c_grid)) for e in evals]45c_cf = [float(c_closed(e)) for e in evals]46V_num_1 = float(np.interp(0.0, ly, V))47on_grid = [bool(np.any(np.isclose(y, e, rtol=1e-12))) for e in evals]4849# diagnostics over the whole grid50max_err_grid = float(np.max(np.abs(c_grid - c_closed(y))))51max_V_err_grid = float(np.max(np.abs(V - V_closed(y))))5253issues = []54if not all(on_grid):55 issues.append({"issue": "Evaluation points " + ", ".join(str(e) for e, g in zip(evals, on_grid) if not g)56 + " are not grid nodes of the prescribed log grid.",57 "handling": "VFI policy and value linearly interpolated in ln y between adjacent grid nodes, as TASK.md allows."})58issues.append({"issue": "TASK does not prescribe the choice set for VFI.",59 "handling": "Discrete choice of next-period output y' on the same 400-point grid (c = y - y'^(1/alpha)); no off-grid optimisation. Initial guess V0 = 0."})6061res = {62 "task_id": "F5-conf-04",63 "c_closed_y05": c_cf[0], "c_closed_y1": c_cf[1], "c_closed_y2": c_cf[2],64 "c_num_y05": c_num[0], "c_num_y1": c_num[1], "c_num_y2": c_num[2],65 "max_abs_error_policy": float(max(abs(a - b) for a, b in zip(c_num, c_cf))),66 "value_closed_y1": float(V_closed(1.0)),67 "value_num_y1": V_num_1,68 "savings_rate_closed_form": float(ab),69 "contraction_modulus": beta,70 "grid_points_used": n,71 "n_iterations": it,72 "final_sup_change": sup,73 "feasibility": "complete",74 "feasibility_reason": None,75 "data_issues": issues,76 "n_input_rows": {"model.csv": int(len(df))},77 "n_analysed": n,78 "random_seed": None,79 "diagnostics": {80 "max_abs_policy_error_on_grid": max_err_grid,81 "max_abs_value_error_on_grid": max_V_err_grid,82 "value_error_y1": V_num_1 - float(V_closed(1.0)),83 "error_bound_from_stopping_rule": beta * sup / (1 - beta),84 "grid_log_step": float(ly[1] - ly[0]),85 "n_grid_points_choosing_upper_bound": at_upper,86 "n_grid_points_choosing_lower_bound": at_lower,87 "closed_form_A": float(A), "closed_form_B": float(B),88 },89}90os.makedirs("out", exist_ok=True)91with open("out/results.json", "w") as f:92 json.dump(res, f, indent=2)93print(json.dumps(res, indent=2))94