1"""F5-conf-04: VFI vs closed form for Brock-Mirman growth with log utility.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py"""3import json, os4import numpy as np, pandas as pd56df = pd.read_csv("data/model.csv")7p = dict(zip(df["parameter"], df["value"]))8alpha = float(p["alpha"]); beta = float(p["beta"])9y_min = float(p["y_min"]); y_max = float(p["y_max"]); N = int(float(p["grid_points"]))10tol = float(p["vfi_tolerance"]); spacing = str(p["grid_spacing"]).strip()11ev = [float(p["eval_y_1"]), float(p["eval_y_2"]), float(p["eval_y_3"])]12assert spacing == "log"1314# ---- closed form (guess V = A + B ln y, verify) ----15ab = alpha * beta16B = 1.0 / (1.0 - ab)17A = (np.log(1 - ab) + (ab / (1 - ab)) * np.log(ab)) / (1 - beta)18c_closed = lambda y: (1 - ab) * np.asarray(y)19V_closed = lambda y: A + B * np.log(y)2021# ---- VFI: discrete choice of next-period output y' on the grid ----22grid = np.exp(np.linspace(np.log(y_min), np.log(y_max), N))23k_needed = grid ** (1.0 / alpha) # savings needed to reach y'24C = grid[:, None] - k_needed[None, :] # consumption for (y, y')25U = np.where(C > 0, np.log(np.where(C > 0, C, 1.0)), -np.inf)26V = np.zeros(N)27it = 0; sup = np.inf28while sup >= tol:29 RHS = U + beta * V[None, :]30 Vn = RHS.max(axis=1)31 sup = float(np.max(np.abs(Vn - V)))32 V = Vn; it += 133 if it > 100000:34 raise RuntimeError("VFI did not converge")35pol_idx = (U + beta * V[None, :]).argmax(axis=1)36c_grid = grid - k_needed[pol_idx]3738c_num = np.interp(ev, grid, c_grid) # linear interpolation in y39v_num_1 = float(np.interp(1.0, grid, V))40c_cl = c_closed(ev)41err = np.abs(c_num - c_cl)4243# diagnostics (not required keys)44sup_V_err = float(np.max(np.abs(V - V_closed(grid))))45sup_c_err_grid = float(np.max(np.abs(c_grid - c_closed(grid))))4647res = {48 "task_id": "F5-conf-04",49 "c_closed_y05": float(c_cl[0]), "c_closed_y1": float(c_cl[1]), "c_closed_y2": float(c_cl[2]),50 "c_num_y05": float(c_num[0]), "c_num_y1": float(c_num[1]), "c_num_y2": float(c_num[2]),51 "max_abs_error_policy": float(err.max()),52 "value_closed_y1": float(V_closed(1.0)),53 "value_num_y1": v_num_1,54 "savings_rate_closed_form": float(ab),55 "contraction_modulus": beta,56 "grid_points_used": N,57 "n_iterations": it,58 "final_sup_change": sup,59 "feasibility": "complete",60 "feasibility_reason": None,61 "data_issues": [62 {"issue": "Evaluation points 0.5, 1.0, 2.0 are not nodes of the 400-point log grid.",63 "handling": "VFI policy c(y) and value V(y) linearly interpolated in y between adjacent grid nodes, as allowed by the task."},64 {"issue": "Choice variable discretised: next-period output y' restricted to the same grid (discrete-choice VFI); log utility unbounded so -inf assigned to infeasible c<=0.",65 "handling": "Standard grid-search VFI; every grid y has feasible choices (y'=y_min always feasible), so no node is dropped. Discretisation error is reported."},66 ],67 "n_input_rows": {"model.csv": int(len(df))},68 "n_analysed": N,69 "random_seed": None,70 "diagnostics": {71 "sup_abs_error_value_on_grid": sup_V_err,72 "sup_abs_error_policy_on_grid": sup_c_err_grid,73 "abs_error_policy_each_eval": [float(x) for x in err],74 "iteration_error_bound_beta_tol_over_1mbeta": beta * tol / (1 - beta),75 "grid_log_step": float(np.log(grid[1] / grid[0])),76 "steady_state_y": float(ab ** (alpha / (1 - alpha))),77 "A": float(A), "B": float(B),78 },79}80os.makedirs("out", exist_ok=True)81with open("out/results.json", "w") as f:82 json.dump(res, f, indent=2)83print(json.dumps(res, indent=2))84