1"""F5-dev-01: Cobb-Douglas consumer - closed form + independent numerical validation.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Deterministic; no randomness used (random_seed = null)."""4import json, math, os5import numpy as np6import pandas as pd7from scipy.optimize import minimize_scalar, minimize89DATA = "data/parameters.csv"10OUT = "out"1112raw = pd.read_csv(DATA)13n_raw = len(raw)14issues = []15dups = raw["parameter"][raw["parameter"].duplicated()].tolist()16if dups:17 issues.append({"issue": f"duplicate parameters {dups}", "handling": "first occurrence used"})18p = raw.drop_duplicates("parameter").set_index("parameter")["value"].astype(float).to_dict()19required = ["alpha", "p1", "p2", "m"]20missing = [k for k in required if k not in p or not math.isfinite(p[k])]21if missing:22 raise SystemExit(f"missing/non-finite parameters: {missing}")23a, p1, p2, m = (p[k] for k in required)24assert 0 < a < 1 and p1 > 0 and p2 > 0 and m > 0, "parameters outside closed-form validity range"2526def u(x1, x2):27 return a * np.log(x1) + (1 - a) * np.log(x2)2829# ---- closed form (Lagrangian FOCs: a/x1 = lam p1, (1-a)/x2 = lam p2, budget binds -> lam = 1/m)30x1 = a * m / p131x2 = (1 - a) * m / p232v = float(u(x1, x2))33resid = p1 * x1 + p2 * x2 - m34mrs = (a / x1) / ((1 - a) / x2)35price_ratio = p1 / p236share1 = p1 * x1 / m3738# ---- numerical validation (does NOT use the closed form)39# (a) bounded 1-D maximisation along the budget line x2 = (m - p1 x1)/p240f = lambda t: -u(t, (m - p1 * t) / p2)41eps = 1e-12 * m / p142r1 = minimize_scalar(f, bounds=(eps, m / p1 - eps), method="bounded",43 options={"xatol": 1e-12, "maxiter": 1000})44# (b) dense grid on the budget line45grid = np.linspace(m / p1 * 1e-6, m / p1 * (1 - 1e-6), 2_000_001)46ug = u(grid, (m - p1 * grid) / p2)47ig = int(np.argmax(ug))48# (c) 2-D SLSQP with inequality budget constraint (starts at an arbitrary feasible point)49r3 = minimize(lambda z: -u(z[0], z[1]), x0=[0.1 * m / p1, 0.1 * m / p2], method="SLSQP",50 bounds=[(1e-9, None), (1e-9, None)],51 constraints=[{"type": "ineq", "fun": lambda z: m - p1 * z[0] - p2 * z[1]}],52 options={"ftol": 1e-12, "maxiter": 1000})53v_num = float(max(-r1.fun, ug[ig], -r3.fun))54validation = {55 "bounded_scalar": {"x1": float(r1.x), "x2": float((m - p1 * r1.x) / p2), "u": float(-r1.fun), "success": bool(r1.success)},56 "grid": {"n_points": int(grid.size), "x1": float(grid[ig]), "x2": float((m - p1 * grid[ig]) / p2), "u": float(ug[ig])},57 "slsqp_2d": {"x1": float(r3.x[0]), "x2": float(r3.x[1]), "u": float(-r3.fun), "success": bool(r3.success),58 "budget_slack": float(m - p1 * r3.x[0] - p2 * r3.x[1])},59 "abs_diff_v_closed_vs_numerical": abs(v - v_num),60}6162results = {63 "task_id": "F5-dev-01",64 "x1_demand": x1, "x2_demand": x2, "indirect_utility": v,65 "budget_identity_residual": resid, "mrs_at_optimum": mrs, "price_ratio": price_ratio,66 "expenditure_share_good1": share1, "indirect_utility_numerical": v_num,67 "feasibility": "complete", "feasibility_reason": None, "data_issues": issues,68 "n_input_rows": {"parameters.csv": n_raw}, "n_analysed": len(p), "random_seed": None,69 "numerical_validation": validation,70}71os.makedirs(OUT, exist_ok=True)72with open(os.path.join(OUT, "results.json"), "w") as fh:73 json.dump(results, fh, indent=2, sort_keys=False)74pd.DataFrame([{"step": "load", "rule": "read data/parameters.csv; drop duplicate parameter names; require finite alpha,p1,p2,m",75 "n_before": n_raw, "n_after": len(p), "n_removed": n_raw - len(p),76 "rationale": "no rows removed unless duplicated"}]).to_csv(os.path.join(OUT, "exclusions.csv"), index=False)77print(json.dumps(results, indent=2))78