1"""F5-conf-04: VFI vs closed-form policy, optimal growth with log utility.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Deterministic; no randomness 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"]))13spacing = str(p["grid_spacing"]).strip(); tol = float(p["vfi_tolerance"])14eval_y = [float(p["eval_y_1"]), float(p["eval_y_2"]), float(p["eval_y_3"])]15assert spacing == "log" and 0 < alpha < 1 and 0 < beta < 1 and 0 < y_min < y_max and N > 116data_issues = []1718# ---- Closed form (guess-and-verify V = A + B ln y): B = 1/(1-ab), c = (1-ab) y19ab = 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) * y23V_closed = lambda y: A + B * np.log(y)24savings_rate = ab # baseline quantity2526# ---- VFI on prescribed grid: state y_i, choice y'_j on the same grid,27# consumption c_ij = y_i - (y'_j)^(1/alpha) (capital needed to produce y'_j); infeasible if c<=0.28grid = np.exp(np.linspace(np.log(y_min), np.log(y_max), N))29k_needed = grid ** (1.0 / alpha)30C = grid[:, None] - k_needed[None, :]31U = np.where(C > 0, np.log(np.where(C > 0, C, 1.0)), -np.inf)32V = np.zeros(N)33n_it, sup = 0, np.inf34while sup >= tol:35 Vn = np.max(U + beta * V[None, :], axis=1)36 sup = float(np.max(np.abs(Vn - V)))37 V = Vn; n_it += 138 if n_it > 100000:39 raise RuntimeError("VFI did not converge")40jstar = np.argmax(U + beta * V[None, :], axis=1)41c_pol = grid - k_needed[jstar]4243# check policy is interior w.r.t. grid bounds (no corner at the top/bottom of the choice grid)44n_low = int(np.sum(jstar == 0)); n_high = int(np.sum(jstar == N - 1))45if n_low or n_high:46 data_issues.append({"issue": f"VFI choice hits grid bound at {n_low} (low) / {n_high} (high) states",47 "handling": "reported; policy near those states is truncated by the grid"})4849lg = np.log(grid)50interp = lambda f, y: float(np.interp(np.log(y), lg, f)) # linear in ln y between nodes51for y in eval_y:52 if not np.any(np.isclose(grid, y, rtol=0, atol=1e-12)):53 pass # off-node: linear interpolation in ln y (allowed by the task)54c_num = [interp(c_pol, y) for y in eval_y]55c_cf = [float(c_closed(y)) for y in eval_y]56max_err = float(max(abs(a - b) for a, b in zip(c_num, c_cf)))5758res = {59 "task_id": "F5-conf-04",60 "c_closed_y05": c_cf[0], "c_closed_y1": c_cf[1], "c_closed_y2": c_cf[2],61 "c_num_y05": c_num[0], "c_num_y1": c_num[1], "c_num_y2": c_num[2],62 "max_abs_error_policy": max_err,63 "value_closed_y1": float(V_closed(1.0)),64 "value_num_y1": interp(V, 1.0),65 "savings_rate_closed_form": float(savings_rate),66 "contraction_modulus": float(beta),67 "grid_points_used": int(N),68 "n_iterations": int(n_it),69 "final_sup_change": float(sup),70 "feasibility": "complete",71 "feasibility_reason": None,72 "data_issues": data_issues,73 "n_input_rows": {"model.csv": int(len(df))},74 "n_analysed": int(N),75 "random_seed": None,76 # supplementary diagnostics (not required keys)77 "diag_eval_points_on_grid": [bool(np.any(np.isclose(grid, y, rtol=0, atol=1e-12))) for y in eval_y],78 "diag_grid_log_step": float(lg[1] - lg[0]),79 "diag_policy_max_abs_error_on_grid": float(np.max(np.abs(c_pol - c_closed(grid)))),80 "diag_value_max_abs_error_on_grid": float(np.max(np.abs(V - V_closed(grid)))),81 "diag_value_tolerance_bound": float(beta * tol / (1 - beta)),82 "diag_value_abs_error_y1": abs(interp(V, 1.0) - float(V_closed(1.0))),83 "diag_closed_form_A": float(A), "diag_closed_form_B": float(B),84 "diag_choice_at_lower_bound": n_low, "diag_choice_at_upper_bound": n_high,85}86os.makedirs(os.path.join(ROOT, "out"), exist_ok=True)87with open(os.path.join(ROOT, "out", "results.json"), "w") as f:88 json.dump(res, f, indent=2)89print(json.dumps(res, indent=2))90