1"""F5-conf-04: VFI vs closed form, optimal growth with log utility.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Deterministic; no random numbers used (random_seed = null)."""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"]); spacing = str(p["grid_spacing"]).strip()14evals = [float(p["eval_y_1"]), float(p["eval_y_2"]), float(p["eval_y_3"])]15data_issues = []16assert spacing == "log"1718# ---------- closed form (guess V = A + B ln y, verified in report) ----------19ab = alpha * beta20B = 1.0 / (1.0 - ab)21A = (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: A + B * np.log(y)24savings_rate = ab # (y - c)/y2526# ---------- VFI on prescribed grid ----------27grid = np.exp(np.linspace(np.log(y_min), np.log(y_max), N))28# choice = next-period output y' on grid node j; requires capital k = y'^(1/alpha) < y29k_needed = grid ** (1.0 / alpha) # (N,)30C = grid[:, None] - k_needed[None, :] # c for (current i, next j)31feasible = C > 032U = np.full_like(C, -np.inf)33U[feasible] = np.log(C[feasible])34if not feasible.any(axis=1).all():35 data_issues.append({"issue": "some grid states have no feasible choice", "handling": "none"})36V = 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 += 145pol_j = np.argmax(U + beta * V[None, :], axis=1)46c_grid = grid - k_needed[pol_j]47converged = sup_change < tol48# corner checks: does optimal y' hit grid bounds?49hits_lower = bool(np.any(pol_j == 0)); hits_upper = bool(np.any(pol_j == N - 1))50if hits_lower or hits_upper:51 data_issues.append({"issue": f"optimal next-period output hits grid bound (lower={hits_lower}, upper={hits_upper}) for some states",52 "handling": "reported; policy at evaluation points checked to be interior"})5354# interpolation at evaluation points: linear in ln y (grid is uniform in ln y)55lg = np.log(grid)56interp = lambda f, y: float(np.interp(np.log(y), lg, f))57c_num = [interp(c_grid, y) for y in evals]58c_cl = [float(c_closed(y)) for y in evals]59max_err = float(max(abs(a - b) for a, b in zip(c_num, c_cl)))60V_num1 = interp(V, 1.0)6162# diagnostics (not required keys)63grid_err_c = float(np.max(np.abs(c_grid - c_closed(grid))))64grid_err_V = float(np.max(np.abs(V - V_closed(grid))))65on_node = [bool(np.min(np.abs(lg - np.log(y))) < 1e-12) for y in evals]6667res = {68 "task_id": "F5-conf-04",69 "c_closed_y05": c_cl[0], "c_closed_y1": c_cl[1], "c_closed_y2": c_cl[2],70 "c_num_y05": c_num[0], "c_num_y1": c_num[1], "c_num_y2": c_num[2],71 "max_abs_error_policy": max_err,72 "value_closed_y1": float(V_closed(1.0)),73 "value_num_y1": V_num1,74 "savings_rate_closed_form": float(savings_rate),75 "contraction_modulus": beta,76 "grid_points_used": int(N),77 "n_iterations": int(n_iter),78 "final_sup_change": sup_change,79 "feasibility": "complete" if converged else "partial",80 "feasibility_reason": None if converged else "VFI did not reach tolerance within max_iter",81 "data_issues": data_issues,82 "n_input_rows": {"model.csv": int(len(df))},83 "n_analysed": int(N),84 "random_seed": None,85 "diagnostics": {86 "eval_points_are_grid_nodes": on_node,87 "max_abs_policy_error_on_grid": grid_err_c,88 "max_abs_value_error_on_grid": grid_err_V,89 "value_error_bound_from_tolerance": beta * tol / (1 - beta),90 "log_grid_step": float(lg[1] - lg[0]),91 "closed_form_A": float(A), "closed_form_B": float(B),92 "policy_hits_lower_bound": hits_lower, "policy_hits_upper_bound": hits_upper,93 "interpolation": "linear in ln(y) on grid nodes",94 },95}96os.makedirs(os.path.join(ROOT, "out"), exist_ok=True)97with open(os.path.join(ROOT, "out", "results.json"), "w") as f:98 json.dump(res, f, indent=2)99print(json.dumps(res, indent=2))100