1"""F5-dev-01: Cobb-Douglas consumer. Closed-form Marshallian demands + independent2numerical validation. Reads only data/parameters.csv; writes out/results.json and3out/exclusions.csv. Deterministic: no randomness is used (random_seed = null)."""4import json, os, hashlib5import numpy as np6import pandas as pd7from scipy.optimize import minimize_scalar, minimize89ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))10DATA = os.path.join(ROOT, "data", "parameters.csv")11OUT = os.path.join(ROOT, "out")1213# ---------------- data audit ----------------14raw = pd.read_csv(DATA)15n_raw = len(raw)16issues = []17required = ["alpha", "p1", "p2", "m"]18dups = raw["parameter"][raw["parameter"].duplicated()].tolist()19if dups:20 issues.append({"issue": f"duplicate parameter rows: {dups}", "handling": "fatal; not resolved"})21missing = [k for k in required if k not in set(raw["parameter"])]22if missing:23 issues.append({"issue": f"missing parameters: {missing}", "handling": "fatal"})24vals = pd.to_numeric(raw["value"], errors="coerce")25if vals.isna().any():26 issues.append({"issue": "non-numeric values", "handling": "fatal"})27extra = [k for k in raw["parameter"] if k not in required]28if extra:29 issues.append({"issue": f"unexpected parameters {extra}", "handling": "ignored"})30used = raw[raw["parameter"].isin(required)]31pd.DataFrame([{"step": 1, "rule": "keep rows whose parameter is in {alpha,p1,p2,m}",32 "n_before": n_raw, "n_after": len(used), "n_removed": n_raw - len(used),33 "rationale": "only these parameters enter the model"}]).to_csv(34 os.path.join(OUT, "exclusions.csv"), index=False)3536P = dict(zip(used["parameter"], pd.to_numeric(used["value"]).astype(float)))37a, p1, p2, m = P["alpha"], P["p1"], P["p2"], P["m"]38valid = (0 < a < 1) and p1 > 0 and p2 > 0 and m > 039if not valid:40 issues.append({"issue": "parameters outside 0<alpha<1, p>0, m>0", "handling": "closed form not reported"})4142def u(x1, x2):43 return a * np.log(x1) + (1 - a) * np.log(x2)4445# ---------------- independent numerical maximisation (no closed form) ----------------46# (i) 1-D: substitute the budget line x2 = (m - p1 x1)/p2, maximise over x1 in (0, m/p1)47# with a bounded Brent search (strict monotonicity => budget binds; checked below in (iii)).48# (ii) dense grid on the budget line (1,000,001 points) as a derivative-free check.49# (iii) 2-D SLSQP with the inequality budget constraint (does not impose binding).50def neg_u_line(x1):51 return -u(x1, (m - p1 * x1) / p2)52eps = 1e-12 * m / p153r1 = minimize_scalar(neg_u_line, bounds=(eps, m / p1 - eps), method="bounded",54 options={"xatol": 1e-12})55x1_num, v_num = float(r1.x), float(-r1.fun)56x2_num = (m - p1 * x1_num) / p257grid = np.linspace(0, m / p1, 1_000_001)[1:-1]58vg = u(grid, (m - p1 * grid) / p2)59ig = int(np.argmax(vg)); x1_grid, v_grid = float(grid[ig]), float(vg[ig])60x0 = np.array([0.5 * m / (2 * p1), 0.5 * m / (2 * p2)]) # arbitrary feasible start61r3 = minimize(lambda x: -u(x[0], x[1]), x0, method="SLSQP",62 bounds=[(1e-9, None), (1e-9, None)],63 constraints=[{"type": "ineq", "fun": lambda x: m - p1 * x[0] - p2 * x[1]}],64 options={"ftol": 1e-14, "maxiter": 1000})65x_slsqp, v_slsqp = r3.x.astype(float), float(-r3.fun)6667# ---------------- closed form (derived from FOCs; see report) ----------------68x1 = a * m / p169x2 = (1 - a) * m / p270v = float(u(x1, x2))71resid = p1 * x1 + p2 * x2 - m72mrs = (a / x1) / ((1 - a) / x2)7374results = {75 "task_id": "F5-dev-01",76 "x1_demand": x1 if valid else None,77 "x2_demand": x2 if valid else None,78 "indirect_utility": v if valid else None,79 "budget_identity_residual": resid if valid else None,80 "mrs_at_optimum": mrs if valid else None,81 "price_ratio": p1 / p2,82 "expenditure_share_good1": p1 * x1 / m if valid else None,83 "indirect_utility_numerical": v_num,84 "feasibility": "complete" if valid and not missing and not dups else "infeasible",85 "feasibility_reason": None if valid else "parameters violate assumptions",86 "data_issues": issues,87 "n_input_rows": {"parameters.csv": int(n_raw)},88 "n_analysed": int(len(used)),89 "random_seed": None,90 "validation_details": {91 "input_sha256": hashlib.sha256(open(DATA, "rb").read()).hexdigest(),92 "brent_x1": x1_num, "brent_x2": x2_num, "brent_converged": bool(r1.success),93 "grid_points": int(grid.size), "grid_x1": x1_grid, "grid_utility": v_grid,94 "grid_step_x1": float(grid[1] - grid[0]),95 "slsqp_x1": float(x_slsqp[0]), "slsqp_x2": float(x_slsqp[1]),96 "slsqp_utility": v_slsqp, "slsqp_success": bool(r3.success),97 "slsqp_budget_slack": float(m - p1 * x_slsqp[0] - p2 * x_slsqp[1]),98 "abs_diff_utility_closed_vs_brent": abs(v - v_num),99 "abs_diff_x1_closed_vs_brent": abs(x1 - x1_num),100 "abs_diff_utility_closed_vs_grid": abs(v - v_grid),101 "abs_diff_utility_closed_vs_slsqp": abs(v - v_slsqp),102 "mrs_minus_price_ratio": mrs - p1 / p2,103 },104}105with open(os.path.join(OUT, "results.json"), "w") as f:106 json.dump(results, f, indent=2, sort_keys=False)107print(json.dumps({k: results[k] for k in ["x1_demand", "x2_demand", "indirect_utility",108 "indirect_utility_numerical", "budget_identity_residual", "mrs_at_optimum", "price_ratio"]}, indent=1))109