1"""F5-dev-01: Cobb-Douglas Marshallian demands, closed form + independent numerical check.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3"""4import json, os5import numpy as np6import pandas as pd7from scipy.optimize import minimize_scalar, minimize89df = pd.read_csv("data/parameters.csv")10n_rows = int(len(df))11issues = []12P = {}13for _, r in df.iterrows():14 name = str(r["parameter"]).strip()15 if name in P:16 issues.append({"issue": f"duplicate parameter {name}", "handling": "first value kept"})17 continue18 P[name] = float(r["value"])19for k in ["alpha", "p1", "p2", "m"]:20 if k not in P or not np.isfinite(P[k]):21 raise SystemExit(f"missing/invalid parameter {k}")22a, p1, p2, m = P["alpha"], P["p1"], P["p2"], P["m"]23# Assumptions for interior unique optimum24assert 0 < a < 1, "alpha must be in (0,1)"25assert p1 > 0 and p2 > 0 and m > 0, "prices and income must be positive"2627def u(x1, x2):28 return a * np.log(x1) + (1 - a) * np.log(x2)2930# ---- closed form (Lagrangian FOCs: a/x1 = lam p1, (1-a)/x2 = lam p2, budget binds) ----31x1 = a * m / p132x2 = (1 - a) * m / p233v = float(u(x1, x2))34resid = p1 * x1 + p2 * x2 - m35mrs = (a / x1) / ((1 - a) / x2)36price_ratio = p1 / p237share1 = p1 * x1 / m3839# ---- numerical validation (does NOT use closed form) ----40# Monotone utility => budget binds; parametrise by x1 on (0, m/p1), x2 = (m - p1 x1)/p2.41def neg_u_line(t):42 return -u(t, (m - p1 * t) / p2)43eps = 1e-12 * m / p144res1 = minimize_scalar(neg_u_line, bounds=(eps, m / p1 - eps), method="bounded",45 options={"xatol": 1e-12})46x1_num, v_num_line = float(res1.x), float(-res1.fun)47x2_num = (m - p1 * x1_num) / p24849# Fine grid search on budget line as secondary check50grid = np.linspace(m / p1 * 1e-6, m / p1 * (1 - 1e-6), 2_000_001)51gv = u(grid, (m - p1 * grid) / p2)52i = int(np.argmax(gv))53x1_grid, v_grid = float(grid[i]), float(gv[i])5455# Inequality-constrained 2-D optimiser (budget as <=, not imposed as equality), starting off-optimum56cons = [{"type": "ineq", "fun": lambda z: m - p1 * z[0] - p2 * z[1]}]57res2 = minimize(lambda z: -u(z[0], z[1]), x0=[1.0, 1.0], method="SLSQP",58 bounds=[(1e-9, None), (1e-9, None)], constraints=cons,59 options={"ftol": 1e-14, "maxiter": 500})60x_slsqp, v_slsqp = res2.x.tolist(), float(-res2.fun)6162v_num = max(v_num_line, v_grid) # baseline numerical maximum on budget line6364diag = {65 "closed_form": {"x1": x1, "x2": x2, "v": v},66 "bounded_scalar": {"x1": x1_num, "x2": x2_num, "v": v_num_line, "success": bool(res1.success)},67 "grid": {"x1": x1_grid, "v": v_grid, "n_points": int(grid.size),68 "spacing": float(grid[1] - grid[0])},69 "slsqp_inequality": {"x": x_slsqp, "v": v_slsqp, "success": bool(res2.success),70 "budget_slack": float(m - p1 * x_slsqp[0] - p2 * x_slsqp[1])},71 "abs_diff_v_closed_minus_numerical": float(v - v_num),72 "abs_diff_x1_closed_minus_scalar": float(x1 - x1_num),73 "abs_diff_x1_closed_minus_grid": float(x1 - x1_grid),74}7576results = {77 "task_id": "F5-dev-01",78 "x1_demand": x1,79 "x2_demand": x2,80 "indirect_utility": v,81 "budget_identity_residual": float(resid),82 "mrs_at_optimum": float(mrs),83 "price_ratio": float(price_ratio),84 "expenditure_share_good1": float(share1),85 "indirect_utility_numerical": v_num,86 "feasibility": "complete",87 "feasibility_reason": None,88 "data_issues": issues,89 "n_input_rows": {"parameters.csv": n_rows},90 "n_analysed": n_rows,91 "random_seed": None,92}93os.makedirs("out", exist_ok=True)94with open("out/results.json", "w") as f:95 json.dump(results, f, indent=2)96with open("out/validation_details.json", "w") as f:97 json.dump(diag, f, indent=2)98print(json.dumps(results, indent=2))99print(json.dumps(diag, indent=2))100