1"""F5-conf-04: VFI vs closed form, optimal growth with log utility.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Reads only data/model.csv; deterministic (no randomness)."""4import json, os, csv5import numpy as np67ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))8DATA = os.path.join(ROOT, "data", "model.csv")9OUT = os.path.join(ROOT, "out")1011with open(DATA, newline="") as f:12 rows = list(csv.DictReader(f))13p = {r["parameter"]: r["value"] for r in rows}14alpha = float(p["alpha"]); beta = float(p["beta"])15y_min = float(p["y_min"]); y_max = float(p["y_max"]); N = int(float(p["grid_points"]))16assert p["grid_spacing"].strip() == "log"17tol = float(p["vfi_tolerance"])18eval_y = [float(p["eval_y_1"]), float(p["eval_y_2"]), float(p["eval_y_3"])]1920# ---------------- closed form (guess-and-verify V = A + B ln y) ----------------21ab = alpha * beta22B = 1.0 / (1.0 - ab)23A = (np.log(1.0 - ab) + ab / (1.0 - ab) * np.log(ab)) / (1.0 - beta)24c_closed = lambda y: (1.0 - ab) * np.asarray(y)25V_closed = lambda y: A + B * np.log(y)26savings_rate = ab2728# ---------------- VFI on prescribed grid (discrete choice of y' on grid) --------29grid = np.exp(np.linspace(np.log(y_min), np.log(y_max), N))30k_needed = grid ** (1.0 / alpha) # capital needed to reach y' = grid[j]31C = grid[:, None] - k_needed[None, :] # c(i, j) = y_i - y'_j^(1/alpha)32feasible = C > 033U = np.full_like(C, -np.inf)34U[feasible] = np.log(C[feasible])35assert feasible.any(axis=1).all(), "some grid point has no feasible choice"3637V = np.zeros(N)38n_iter = 039sup_change = np.inf40max_iter = 10000041while sup_change >= tol and n_iter < max_iter:42 W = U + beta * V[None, :]43 V_new = W.max(axis=1)44 sup_change = float(np.max(np.abs(V_new - V)))45 V = V_new46 n_iter += 147converged = sup_change < tol48pol_j = (U + beta * V[None, :]).argmax(axis=1)49c_pol = grid - k_needed[pol_j]5051# boundary diagnostics: does the optimal choice ever hit grid edges?52hits_lower = int(np.sum(pol_j == 0)); hits_upper = int(np.sum(pol_j == N - 1))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_pol, y) for y in eval_y]58c_cl = [float(c_closed(y)) for y in eval_y]59on_grid = [bool(np.any(np.isclose(grid, y, rtol=1e-12, atol=0))) for y in eval_y]60max_err = float(max(abs(a - b) for a, b in zip(c_num, c_cl)))6162# whole-grid diagnostics (for report only)63grid_policy_err = float(np.max(np.abs(c_pol - c_closed(grid))))64grid_value_err = float(np.max(np.abs(V - V_closed(grid))))6566data_issues = []67if not all(on_grid):68 data_issues.append({"issue": "Evaluation points %s are not grid nodes of the log-spaced grid" %69 [y for y, o in zip(eval_y, on_grid) if not o],70 "handling": "VFI policy and value interpolated linearly in ln(y) between adjacent grid nodes"})71data_issues.append({"issue": "Design choice: the choice variable is next-period output y' restricted to the same grid (discrete-choice VFI); consumption c = y - y'^(1/alpha)",72 "handling": "No deviation from prescribed grid/parameters; discretisation error of order of the grid spacing is reported, not corrected"})7374feas = "complete" if converged else "partial"75res = {76 "task_id": "F5-conf-04",77 "c_closed_y05": c_cl[0], "c_closed_y1": c_cl[1], "c_closed_y2": c_cl[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": float(V_closed(1.0)),81 "value_num_y1": interp(V, 1.0),82 "savings_rate_closed_form": float(savings_rate),83 "contraction_modulus": float(beta),84 "grid_points_used": int(N),85 "n_iterations": int(n_iter),86 "final_sup_change": float(sup_change),87 "feasibility": feas,88 "feasibility_reason": None if converged else "VFI did not reach tolerance within max_iter",89 "data_issues": data_issues,90 "n_input_rows": {"model.csv": len(rows)},91 "n_analysed": int(N),92 "random_seed": None,93 # supplementary diagnostics94 "diag_grid_log_step": float(lg[1] - lg[0]),95 "diag_max_abs_policy_error_on_grid": grid_policy_err,96 "diag_max_abs_value_error_on_grid": grid_value_err,97 "diag_policy_hits_lower_bound": hits_lower,98 "diag_policy_hits_upper_bound": hits_upper,99 "diag_value_error_y1": interp(V, 1.0) - float(V_closed(1.0)),100 "diag_closed_form_A": float(A), "diag_closed_form_B": float(B),101 "diag_vfi_error_bound_beta_over_1mbeta_times_final_change": float(beta / (1 - beta) * sup_change),102}103os.makedirs(OUT, exist_ok=True)104with open(os.path.join(OUT, "results.json"), "w") as f:105 json.dump(res, f, indent=2)106with open(os.path.join(OUT, "exclusions.csv"), "w") as f:107 f.write("step,rule,n_before,n_after,n_removed,rationale\n")108 f.write("1,none (all parameter rows used),%d,%d,0,model.csv contains only parameters; nothing excluded\n" % (len(rows), len(rows)))109print(json.dumps({k: v for k, v in res.items() if k != "data_issues"}, indent=1))110