1"""F5-conf-04: VFI vs closed form for log-utility Cobb-Douglas growth (full depreciation).2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3"""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 = float(p["alpha"]); beta = float(p["beta"])12y_min = float(p["y_min"]); y_max = float(p["y_max"]); N = 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"])]1617# ---------- closed form (guess V = A + B ln y, verified in report) ----------18ab = alpha * beta19B = 1.0 / (1.0 - ab)20A = (np.log(1.0 - ab) + (ab / (1.0 - ab)) * np.log(ab)) / (1.0 - beta)21c_closed = lambda y: (1.0 - ab) * np.asarray(y)22V_closed = lambda y: A + B * np.log(y)2324# ---------- VFI: discrete choice of next-period output y' on the grid ----------25lny = np.linspace(np.log(y_min), np.log(y_max), N)26y = np.exp(lny)27k_needed = y ** (1.0 / alpha) # savings needed to reach y' = grid[j]28C = y[:, None] - k_needed[None, :] # consumption for (state i, choice j)29feasible = C > 030U = np.full_like(C, -np.inf)31U[feasible] = np.log(C[feasible])3233V = np.zeros(N)34n_it = 035while True:36 RHS = U + beta * V[None, :]37 Vn = RHS.max(axis=1)38 sup = float(np.max(np.abs(Vn - V)))39 V = Vn40 n_it += 141 if sup < tol:42 break43pol_j = (U + beta * V[None, :]).argmax(axis=1)44c_grid = C[np.arange(N), pol_j]4546# check whether the choice ever hits grid bounds47hits_lower = bool(np.any(pol_j == 0)); hits_upper = bool(np.any(pol_j == N - 1))4849# interpolate (linear in ln y) at evaluation points50interp = lambda f, x: float(np.interp(np.log(x), lny, f))51c_num = [interp(c_grid, e) for e in evals]52c_cl = [float(c_closed(e)) for e in evals]53max_err = float(max(abs(a - b) for a, b in zip(c_num, c_cl)))54v_num_1 = interp(V, 1.0)55v_cl_1 = float(V_closed(1.0))5657# whole-grid diagnostics (not required keys)58sup_policy_err_grid = float(np.max(np.abs(c_grid - c_closed(y))))59sup_value_err_grid = float(np.max(np.abs(V - V_closed(y))))6061issues = []62for e in evals:63 if not np.any(np.isclose(y, e, rtol=0, atol=1e-12)):64 issues.append({"issue": f"evaluation point y={e} is not a grid node",65 "handling": "VFI policy/value linearly interpolated in ln y between adjacent grid nodes"})66issues.append({"issue": "VFI choice variable not specified in model.csv",67 "handling": "next-period output y' chosen from the same 400-point log grid (discrete choice, no off-grid interpolation of V); V initialised at 0"})6869res = {70 "task_id": "F5-conf-04",71 "c_closed_y05": c_cl[0], "c_closed_y1": c_cl[1], "c_closed_y2": c_cl[2],72 "c_num_y05": c_num[0], "c_num_y1": c_num[1], "c_num_y2": c_num[2],73 "max_abs_error_policy": max_err,74 "value_closed_y1": v_cl_1, "value_num_y1": v_num_1,75 "savings_rate_closed_form": float(ab),76 "contraction_modulus": float(beta),77 "grid_points_used": int(N),78 "n_iterations": int(n_it),79 "final_sup_change": sup,80 "feasibility": "complete", "feasibility_reason": None,81 "data_issues": issues,82 "n_input_rows": {"model.csv": int(len(df))},83 "n_analysed": int(N),84 "random_seed": None,85 "diagnostics": {86 "grid_log_step": float(lny[1] - lny[0]),87 "sup_abs_policy_error_on_grid": sup_policy_err_grid,88 "sup_abs_value_error_on_grid": sup_value_err_grid,89 "value_error_bound_from_tolerance": float(beta * tol / (1 - beta)),90 "policy_hits_lower_grid_bound": hits_lower,91 "policy_hits_upper_grid_bound": hits_upper,92 "closed_form_A": float(A), "closed_form_B": float(B),93 "steady_state_y": float(ab ** (alpha / (1 - alpha))),94 },95}96with open(os.path.join(ROOT, "out", "results.json"), "w") as f:97 json.dump(res, f, indent=2)98print(json.dumps(res, indent=2))99