1"""F5-dev-01: Cobb-Douglas consumer. Closed-form Marshallian demands + independent numerical validation.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Reads only data/parameters.csv; writes out/results.json. Fully deterministic (no randomness; random_seed = null).4"""5import json, math, os6import numpy as np7import pandas as pd8from scipy.optimize import minimize_scalar, minimize910DATA = os.path.join("data", "parameters.csv")11OUT = os.path.join("out", "results.json")1213raw = pd.read_csv(DATA)14n_raw = len(raw)15issues = []16if raw["parameter"].duplicated().any():17 issues.append({"issue": "duplicate parameter names", "handling": "first occurrence used"})18par = raw.drop_duplicates("parameter").set_index("parameter")["value"].astype(float)19needed = ["alpha", "p1", "p2", "m"]20missing = [k for k in needed if k not in par.index or not np.isfinite(par[k])]21if missing:22 raise SystemExit(f"missing/non-finite parameters: {missing}")23alpha, p1, p2, m = (float(par[k]) for k in needed)24if not (0 < alpha < 1 and p1 > 0 and p2 > 0 and m > 0):25 raise SystemExit("parameters outside the range where the interior closed form is the unique optimum")26n_analysed = int(sum(raw["parameter"].isin(needed)))272829def u(x1, x2):30 return alpha * np.log(x1) + (1 - alpha) * np.log(x2)313233# ---------- Closed form (Lagrangian FOCs: alpha/x1 = lam p1, (1-alpha)/x2 = lam p2, budget binds => lam = 1/m) ----------34x1 = alpha * m / p135x2 = (1 - alpha) * m / p236v = float(u(x1, x2))37resid = p1 * x1 + p2 * x2 - m38mrs = (alpha / x1) / ((1 - alpha) / x2)39price_ratio = p1 / p240share1 = p1 * x1 / m4142# ---------- Independent numerical validation (does NOT use the closed form) ----------43# Budget line (monotone utility => budget binds): x1 in (0, m/p1), x2 = (m - p1 x1)/p2.44def neg_u_line(x1v):45 return -u(x1v, (m - p1 * x1v) / p2)4647ub = m / p148# (a) dense grid search on the budget line49grid = np.linspace(ub * 1e-6, ub * (1 - 1e-6), 2_000_001)50vals = -neg_u_line(grid)51i = int(np.argmax(vals))52x1_grid, v_grid = float(grid[i]), float(vals[i])53# (b) bounded scalar optimiser (Brent) on the budget line54res_b = minimize_scalar(neg_u_line, bounds=(ub * 1e-9, ub * (1 - 1e-9)), method="bounded",55 options={"xatol": 1e-12})56x1_brent, v_brent = float(res_b.x), float(-res_b.fun)57# (c) 2-D SLSQP with the inequality budget constraint (checks that the budget binds without imposing it)58res_s = minimize(lambda z: -u(z[0], z[1]), x0=np.array([1.0, 1.0]), method="SLSQP",59 bounds=[(1e-9, None), (1e-9, None)],60 constraints=[{"type": "ineq", "fun": lambda z: m - p1 * z[0] - p2 * z[1]}],61 options={"ftol": 1e-14, "maxiter": 1000})62v_slsqp = float(-res_s.fun)6364v_num = max(v_grid, v_brent, v_slsqp) # best utility found numerically65x1_num = x1_brent66x2_num = (m - p1 * x1_num) / p26768results = {69 "task_id": "F5-dev-01",70 "x1_demand": x1,71 "x2_demand": x2,72 "indirect_utility": v,73 "budget_identity_residual": resid,74 "mrs_at_optimum": mrs,75 "price_ratio": price_ratio,76 "expenditure_share_good1": share1,77 "indirect_utility_numerical": v_num,78 "feasibility": "complete",79 "feasibility_reason": None,80 "data_issues": issues,81 "n_input_rows": {"parameters.csv": n_raw},82 "n_analysed": n_analysed,83 "random_seed": None,84 # diagnostics of the numerical validation85 "numerical_x1_brent": x1_num,86 "numerical_x2_brent": x2_num,87 "numerical_x1_grid": x1_grid,88 "numerical_utility_grid": v_grid,89 "numerical_utility_brent": v_brent,90 "numerical_utility_slsqp": v_slsqp,91 "numerical_x_slsqp": [float(res_s.x[0]), float(res_s.x[1])],92 "numerical_slsqp_budget_slack": float(m - p1 * res_s.x[0] - p2 * res_s.x[1]),93 "utility_gap_numerical_minus_closed_form": v_num - v,94 "max_abs_demand_gap_x1": max(abs(x1_num - x1), abs(x1_grid - x1), abs(float(res_s.x[0]) - x1)),95 "max_abs_demand_gap_x2": max(abs(x2_num - x2), abs(float(res_s.x[1]) - x2)),96}97with open(OUT, "w") as f:98 json.dump(results, f, indent=2)99print(f"x1*={x1:.6f} x2*={x2:.6f} v={v:.10f} v_num={v_num:.10f} residual={resid:.3e} MRS={mrs:.6f} p1/p2={price_ratio:.6f}")100