1"""F5-conf-04: VFI vs closed form for the Brock-Mirman optimal growth model (log utility).2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3"""4import json, os5import numpy as np6import pandas as pd78HERE = os.path.dirname(os.path.abspath(__file__))9ROOT = os.path.dirname(HERE)10df = pd.read_csv(os.path.join(ROOT, "data", "model.csv"))11p = dict(zip(df["parameter"], df["value"]))12alpha = float(p["alpha"]); beta = float(p["beta"])13y_min = float(p["y_min"]); y_max = float(p["y_max"]); N = int(float(p["grid_points"]))14spacing = str(p["grid_spacing"]).strip(); tol = float(p["vfi_tolerance"])15evals = [float(p["eval_y_1"]), float(p["eval_y_2"]), float(p["eval_y_3"])]16data_issues = [{"issue": "Evaluation points y=0.5, 1.0, 2.0 (and y=1 for V) are not nodes of the prescribed 400-point log grid.",17 "handling": "VFI policy and value are linearly interpolated in ln y between adjacent grid nodes (as allowed by TASK.md); grid itself unchanged."},18 {"issue": "VFI implemented as discrete choice of next-period output y' restricted to the same grid (consumption c = y - y'^(1/alpha)); infeasible choices (c<=0) excluded.",19 "handling": "Design choice, not a deviation from prescribed parameters; it introduces a discretisation error in the policy reported as max_abs_error_policy."}]20assert spacing == "log"2122# ---- closed form (guess V = a + b ln y, verified in report) ----23ab = alpha * beta24b = 1.0 / (1.0 - ab)25a = (np.log(1 - ab) + (ab / (1 - ab)) * np.log(ab)) / (1 - beta)26c_closed = lambda y: (1 - ab) * np.asarray(y)27V_closed = lambda y: a + b * np.log(np.asarray(y))2829# ---- VFI: discrete choice of next-period output y' on the same grid ----30grid = np.exp(np.linspace(np.log(y_min), np.log(y_max), N))31k_needed = grid ** (1.0 / alpha) # savings needed to reach y'32C = grid[:, None] - k_needed[None, :] # consumption for (y, y')33feasible = C > 034U = np.where(feasible, np.log(np.where(feasible, C, 1.0)), -np.inf)3536V = np.zeros(N)37n_iter = 038while True:39 RHS = U + beta * V[None, :]40 V_new = RHS.max(axis=1)41 n_iter += 142 diff = np.max(np.abs(V_new - V))43 V = V_new44 if diff < tol:45 break46 if n_iter > 100000:47 raise RuntimeError("VFI did not converge")48pol_idx = (U + beta * V[None, :]).argmax(axis=1)49c_grid = grid - k_needed[pol_idx]5051# check that the optimal choice is not at a grid boundary (truncation)52n_at_lower = int(np.sum(pol_idx == 0)); n_at_upper = int(np.sum(pol_idx == N - 1))5354lg = np.log(grid)55interp = lambda f, y: float(np.interp(np.log(y), lg, f)) # linear in ln y56on_node = [bool(np.any(np.isclose(grid, y, rtol=1e-12, atol=0))) for y in evals]5758c_num = [interp(c_grid, y) for y in evals]59c_cf = [float(c_closed(y)) for y in evals]60max_err = float(max(abs(x - z) for x, z in zip(c_num, c_cf)))61v_num_1 = interp(V, 1.0)62v_cf_1 = float(V_closed(1.0))6364# diagnostics (not required keys)65diag = {66 "eval_points_are_grid_nodes": on_node,67 "policy_at_lower_bound_count": n_at_lower,68 "policy_at_upper_bound_count": n_at_upper,69 "max_abs_policy_error_over_grid": float(np.max(np.abs(c_grid - c_closed(grid)))),70 "max_abs_value_error_over_grid": float(np.max(np.abs(V - V_closed(grid)))),71 "grid_log_step": float(lg[1] - lg[0]),72 "closed_form_a": float(a), "closed_form_b": float(b),73}7475res = {76 "task_id": "F5-conf-04",77 "c_closed_y05": c_cf[0], "c_closed_y1": c_cf[1], "c_closed_y2": c_cf[2],78 "c_num_y05": c_num[0], "c_num_y1": c_num[1], "c_num_y2": c_num[2],79 "max_abs_error_policy": max_err,80 "value_closed_y1": v_cf_1, "value_num_y1": v_num_1,81 "savings_rate_closed_form": float(ab),82 "contraction_modulus": float(beta),83 "grid_points_used": int(N),84 "n_iterations": int(n_iter),85 "final_sup_change": float(diff),86 "feasibility": "complete", "feasibility_reason": None,87 "data_issues": data_issues,88 "n_input_rows": {"model.csv": int(len(df))},89 "n_analysed": int(N),90 "random_seed": None,91 "diagnostics": diag,92}93with open(os.path.join(HERE, "results.json"), "w") as f:94 json.dump(res, f, indent=2)95print(json.dumps(res, indent=2))96