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; 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"]))13tol = float(p["vfi_tolerance"])14assert str(p["grid_spacing"]).strip() == "log"15evals = [float(p["eval_y_1"]), float(p["eval_y_2"]), float(p["eval_y_3"])]16data_issues = []1718# ---------- Closed form (guess-and-verify: V = A0 + B ln y) ----------19ab = alpha * beta20B = 1.0 / (1.0 - ab)21A0 = (np.log(1.0 - ab) + ab / (1.0 - ab) * np.log(ab)) / (1.0 - beta)22c_closed = lambda y: (1.0 - ab) * np.asarray(y)23V_closed = lambda y: A0 + B * np.log(np.asarray(y))24savings_rate = ab # (y - c)/y2526# ---------- VFI on the prescribed grid (discrete choice of next-period output on the grid) ----------27grid = np.exp(np.linspace(np.log(y_min), np.log(y_max), N))28k_needed = grid ** (1.0 / alpha) # capital required to reach each y' on the grid29C = grid[:, None] - k_needed[None, :] # consumption for (y today, y' tomorrow)30feas = C > 031U = np.full(C.shape, -np.inf)32U[feas] = np.log(C[feas])33V = np.zeros(N)34n_iter, sup_change = 0, np.inf35while sup_change >= tol:36 RHS = U + beta * V[None, :]37 V_new = RHS.max(axis=1)38 sup_change = float(np.max(np.abs(V_new - V)))39 V = V_new40 n_iter += 141 if n_iter > 100000:42 raise RuntimeError("VFI did not converge")43pol_idx = (U + beta * V[None, :]).argmax(axis=1)44c_pol = grid - k_needed[pol_idx]4546# Diagnostics: is the grid bound binding? (policy at the edge of the y' grid)47edge_hits = int(np.sum((pol_idx == 0) | (pol_idx == N - 1)))48yprime_closed = (ab * grid) ** alpha49if yprime_closed.min() < y_min or yprime_closed.max() > y_max:50 data_issues.append({"issue": "closed-form next-period output leaves the prescribed grid",51 "handling": "reported; truncation may bias VFI near bounds"})5253# Interpolation at evaluation points: c linear in y; V linear in ln y54c_num = [float(np.interp(y, grid, c_pol)) for y in evals]55cc = [float(c_closed(y)) for y in evals]56on_node = [bool(np.any(np.isclose(grid, y, rtol=0, atol=1e-12))) for y in evals]57V_num_1 = float(np.interp(np.log(1.0), np.log(grid), V))5859# Theoretical bound on the remaining VFI error from the stopping rule60vfi_err_bound = beta / (1 - beta) * sup_change6162results = {63 "task_id": "F5-conf-04",64 "c_closed_y05": cc[0], "c_closed_y1": cc[1], "c_closed_y2": cc[2],65 "c_num_y05": c_num[0], "c_num_y1": c_num[1], "c_num_y2": c_num[2],66 "max_abs_error_policy": float(max(abs(a - b) for a, b in zip(c_num, cc))),67 "value_closed_y1": float(V_closed(1.0)),68 "value_num_y1": V_num_1,69 "savings_rate_closed_form": float(savings_rate),70 "contraction_modulus": float(beta),71 "grid_points_used": int(N),72 "n_iterations": int(n_iter),73 "final_sup_change": float(sup_change),74 "feasibility": "complete",75 "feasibility_reason": None,76 "data_issues": data_issues,77 "n_input_rows": {"model.csv": int(len(df))},78 "n_analysed": int(N),79 "random_seed": None,80 # supplementary diagnostics81 "baseline_closed_form_B_coefficient": float(B),82 "diag_eval_points_on_grid_node": on_node,83 "diag_policy_edge_hits": edge_hits,84 "diag_closed_yprime_min": float(yprime_closed.min()),85 "diag_closed_yprime_max": float(yprime_closed.max()),86 "diag_grid_log_step": float(np.log(grid[1] / grid[0])),87 "diag_vfi_tolerance_error_bound": float(vfi_err_bound),88 "diag_value_abs_error_y1": float(abs(V_num_1 - float(V_closed(1.0)))),89 "diag_value_sup_error_grid": float(np.max(np.abs(V - V_closed(grid)))),90 "diag_policy_sup_error_grid": float(np.max(np.abs(c_pol - c_closed(grid)))),91 "diag_policy_max_rel_error_eval": float(max(abs(a - b) / b for a, b in zip(c_num, cc))),92}93os.makedirs(os.path.join(ROOT, "out"), exist_ok=True)94with open(os.path.join(ROOT, "out", "results.json"), "w") as f:95 json.dump(results, f, indent=2)96for k, v in results.items():97 print(f"{k}: {v}")98