1"""F5-conf-02: Monte Carlo bias of OLS AR(1) coefficient vs Kendall/Marriott-Pope approximation.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3"""4import json, os5import numpy as np6import pandas as pd78d = pd.read_csv("data/design.csv")9p = dict(zip(d["parameter"], d["value"]))10n = int(p["n"]); rho = float(p["rho"]); sigma = float(p["sigma"])11R = int(p["n_replications"]); seed = int(p["seed"])12issues = []13assert str(p["intercept"]).strip().lower() == "yes"14assert str(p["initial_condition"]).strip().lower() == "stationary"15if not (abs(rho) < 1 and n > 2 and R > 1 and sigma > 0):16 issues.append({"issue": "invalid design value", "handling": "none"})1718rng = np.random.default_rng(seed)19# Draw order: y_1 for all replications, then innovations eps_2..eps_n (R x (n-1)).20y = np.empty((R, n))21y[:, 0] = rng.normal(0.0, sigma / np.sqrt(1 - rho**2), size=R)22eps = rng.normal(0.0, sigma, size=(R, n - 1))23for t in range(1, n):24 y[:, t] = rho * y[:, t - 1] + eps[:, t - 1] # c = 02526x = y[:, :-1]; z = y[:, 1:] # regressor y_{t-1}, response y_t, t=2..n27xd = x - x.mean(axis=1, keepdims=True)28zd = z - z.mean(axis=1, keepdims=True)29rho_hat = (xd * zd).sum(1) / (xd**2).sum(1) # OLS slope with intercept3031# sanity check vs. explicit lstsq on first replication32X0 = np.column_stack([np.ones(n - 1), x[0]])33assert np.isclose(np.linalg.lstsq(X0, z[0], rcond=None)[0][1], rho_hat[0])3435ok = np.isfinite(rho_hat)36rh = rho_hat[ok]37mean = float(rh.mean()); sd = float(rh.std(ddof=1)); se = sd / np.sqrt(rh.size)38bias_an = -(1 + 3 * rho) / n39bias_mc = mean - rho40# For reference: same formula using the effective regression sample size n-141bias_an_nm1 = -(1 + 3 * rho) / (n - 1)4243res = {44 "task_id": "F5-conf-02",45 "bias_analytical": bias_an,46 "bias_mc": bias_mc,47 "bias_mc_se": se,48 "rho_hat_mean": mean,49 "rho_hat_sd_mc": sd,50 "discrepancy_se_units": (bias_mc - bias_an) / se,51 "se_asymptotic": float(np.sqrt((1 - rho**2) / n)),52 "n_replications": R,53 "feasibility": "complete",54 "feasibility_reason": None,55 "data_issues": issues,56 "n_input_rows": {"design.csv": int(len(d))},57 "n_analysed": int(rh.size),58 "random_seed": seed,59 # supplementary (not required)60 "bias_analytical_n_minus_1": bias_an_nm1,61 "discrepancy_se_units_n_minus_1": (bias_mc - bias_an_nm1) / se,62 "rho_hat_median": float(np.median(rh)),63 "share_rho_hat_below_rho": float((rh < rho).mean()),64 "relative_discrepancy": (bias_mc - bias_an) / abs(bias_an),65}66os.makedirs("out", exist_ok=True)67with open("out/results.json", "w") as f:68 json.dump(res, f, indent=2)69print(json.dumps(res, indent=2))70