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.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 = 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"]))13spacing = str(p["grid_spacing"]).strip(); tol = float(p["vfi_tolerance"])14evals = [float(p["eval_y_1"]), float(p["eval_y_2"]), float(p["eval_y_3"])]15data_issues = []16assert spacing == "log"17assert 0 < alpha < 1 and 0 < beta < 11819# ---------------- closed form (guess V = A + B ln y, verified in report) ----------------20ab = alpha * beta21B = 1.0 / (1.0 - ab)22A = (np.log(1.0 - ab) + (ab / (1.0 - ab)) * np.log(ab)) / (1.0 - beta)23c_closed = lambda y: (1.0 - ab) * y24V_closed = lambda y: A + B * np.log(y)25savings_rate = ab # (y - c)/y2627# ---------------- VFI on the prescribed grid (discrete choice of next-period y on the grid) ----28grid = np.exp(np.linspace(np.log(y_min), np.log(y_max), N))29k_needed = grid ** (1.0 / alpha) # investment needed to reach y' = grid[j]30C = grid[:, None] - k_needed[None, :] # consumption for (y_i, y'_j)31feasible = C > 032U = np.full(C.shape, -np.inf)33U[feasible] = np.log(C[feasible])34n_infeasible_pairs = int((~feasible).sum())3536V = np.zeros(N)37n_iter = 038sup_change = np.inf39max_iter = 10000040while sup_change >= tol and n_iter < max_iter:41 V_new = np.max(U + beta * V[None, :], axis=1)42 sup_change = float(np.max(np.abs(V_new - V)))43 V = V_new44 n_iter += 145converged = sup_change < tol46pol_idx = np.argmax(U + beta * V[None, :], axis=1)47c_pol = grid - k_needed[pol_idx]48bound_hits = int(np.sum((pol_idx == 0) | (pol_idx == N - 1)))49if bound_hits:50 data_issues.append({"issue": f"policy hits grid bound at {bound_hits} nodes",51 "handling": "reported; results at such nodes may be distorted by truncation"})5253# interpolation at evaluation points: linear in ln y (grid is uniform in ln y)54lg = np.log(grid)55def interp(vals, y):56 assert y_min <= y <= y_max57 return float(np.interp(np.log(y), lg, vals))58on_node = [bool(np.min(np.abs(lg - np.log(y))) < 1e-12) for y in evals]5960c_num = [interp(c_pol, y) for y in evals]61c_cf = [float(c_closed(y)) for y in evals]62max_err = float(max(abs(a - b) for a, b in zip(c_num, c_cf)))63V_num_1 = interp(V, 1.0)6465# diagnostics (not required keys; printed only)66err_grid_c = float(np.max(np.abs(c_pol - c_closed(grid))))67err_grid_V = float(np.max(np.abs(V - V_closed(grid))))68print(f"grid step in ln y = {lg[1]-lg[0]:.6f}; eval points on nodes: {on_node}")69print(f"VFI iterations={n_iter}, final sup change={sup_change:.3e}, converged={converged}")70print(f"policy bound hits={bound_hits}; infeasible (y,y') pairs={n_infeasible_pairs}")71print(f"max |c_vfi - c_closed| over all grid nodes = {err_grid_c:.3e}; max |V_vfi - V_closed| = {err_grid_V:.3e}")7273if not all(on_node):74 pass # interpolation is allowed by TASK.md; documented in report7576res = {77 "task_id": "F5-conf-04",78 "c_closed_y05": c_cf[0], "c_closed_y1": c_cf[1], "c_closed_y2": c_cf[2],79 "c_num_y05": c_num[0], "c_num_y1": c_num[1], "c_num_y2": c_num[2],80 "max_abs_error_policy": max_err,81 "value_closed_y1": float(V_closed(1.0)),82 "value_num_y1": V_num_1,83 "savings_rate_closed_form": float(savings_rate),84 "contraction_modulus": beta,85 "grid_points_used": int(N),86 "n_iterations": int(n_iter),87 "final_sup_change": float(sup_change),88 "feasibility": "complete" if converged else "partial",89 "feasibility_reason": None if converged else "VFI did not reach tolerance within max_iter",90 "data_issues": data_issues,91 "n_input_rows": {"model.csv": int(len(df))},92 "n_analysed": int(N),93 "random_seed": None,94 "diagnostics": {95 "grid_step_log_y": float(lg[1] - lg[0]),96 "eval_points_on_grid_nodes": on_node,97 "policy_grid_bound_hits": bound_hits,98 "max_abs_policy_error_all_nodes": err_grid_c,99 "max_abs_value_error_all_nodes": err_grid_V,100 "initial_guess": "V0 = 0",101 "choice_variable": "next-period output y' restricted to grid",102 "interpolation": "linear in ln y",103 },104}105os.makedirs(os.path.join(ROOT, "out"), exist_ok=True)106with open(os.path.join(ROOT, "out", "results.json"), "w") as f:107 json.dump(res, f, indent=2)108print(json.dumps({k: v for k, v in res.items() if k != "diagnostics"}, indent=2))109