1"""F5-conf-02: finite-sample bias of OLS AR(1) coefficient (with intercept).2Reads data/design.csv only; writes out/results.json deterministically.3Run from workspace root: ~/wp29/venv/bin/python out/analysis.py4"""5import json, os6import numpy as np7import pandas as pd89ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))10design = pd.read_csv(os.path.join(ROOT, "data", "design.csv"))11p = dict(zip(design["parameter"], design["value"]))1213n = int(p["n"]); rho = float(p["rho"]); sigma = float(p["sigma"])14intercept = str(p["intercept"]).strip().lower() == "yes"15init = str(p["initial_condition"]).strip().lower()16R = int(p["n_replications"]); seed = int(p["seed"])17data_issues = []18assert intercept and init == "stationary" and abs(rho) < 1 and n >= 3 and R > 11920# --- Baseline (closed-form) quantities ---21bias_analytical = -(1.0 + 3.0 * rho) / n22se_asymptotic = float(np.sqrt((1.0 - rho ** 2) / n))23# sensitivity: same formula using effective regression sample size n-124bias_analytical_neff = -(1.0 + 3.0 * rho) / (n - 1)2526# --- Monte Carlo (vectorised) ---27# Draw order: first R stationary initial values, then an R x (n-1) innovation matrix.28rng = np.random.default_rng(seed)29y = np.empty((R, n))30y[:, 0] = rng.standard_normal(R) * sigma / np.sqrt(1.0 - rho ** 2)31eps = rng.standard_normal((R, n - 1)) * sigma32c = 0.033for t in range(1, n):34 y[:, t] = c + rho * y[:, t - 1] + eps[:, t - 1]35x = y[:, :-1]; z = y[:, 1:] # regress y_t (t=2..n) on (1, y_{t-1})36xd = x - x.mean(axis=1, keepdims=True)37zd = z - z.mean(axis=1, keepdims=True)38rho_hat = (xd * zd).sum(axis=1) / (xd ** 2).sum(axis=1)39finite = np.isfinite(rho_hat)40n_nonfinite = int((~finite).sum())41if n_nonfinite:42 data_issues.append({"issue": f"{n_nonfinite} non-finite rho_hat", "handling": "excluded"})43rh = rho_hat[finite]44R_used = int(rh.size)4546# cross-check one replication against numpy lstsq47X0 = np.column_stack([np.ones(n - 1), x[0]])48b0 = np.linalg.lstsq(X0, z[0], rcond=None)[0][1]49assert abs(b0 - rho_hat[0]) < 1e-105051rho_hat_mean = float(rh.mean())52rho_hat_sd = float(rh.std(ddof=1))53bias_mc = rho_hat_mean - rho54bias_mc_se = rho_hat_sd / np.sqrt(R_used)55disc = (bias_mc - bias_analytical) / bias_mc_se5657res = {58 "task_id": "F5-conf-02",59 "bias_analytical": bias_analytical,60 "bias_mc": float(bias_mc),61 "bias_mc_se": float(bias_mc_se),62 "rho_hat_mean": rho_hat_mean,63 "rho_hat_sd_mc": rho_hat_sd,64 "discrepancy_se_units": float(disc),65 "se_asymptotic": se_asymptotic,66 "n_replications": R,67 "feasibility": "complete",68 "feasibility_reason": None,69 "data_issues": data_issues,70 "n_input_rows": {"design.csv": int(len(design))},71 "n_analysed": R_used,72 "random_seed": seed,73 # supplementary keys74 "baseline_bias_analytical": bias_analytical,75 "baseline_se_asymptotic": se_asymptotic,76 "bias_analytical_neff": bias_analytical_neff,77 "discrepancy_se_units_neff": float((bias_mc - bias_analytical_neff) / bias_mc_se),78 "bias_mc_minus_analytical": float(bias_mc - bias_analytical),79 "ratio_bias_mc_to_analytical": float(bias_mc / bias_analytical),80 "rho_hat_median": float(np.median(rh)),81 "share_rho_hat_below_rho": float((rh < rho).mean()),82 "ratio_sd_mc_to_se_asymptotic": float(rho_hat_sd / se_asymptotic),83 "n_regression_obs": n - 1,84 "n_nonfinite_rho_hat": n_nonfinite,85}86os.makedirs(os.path.join(ROOT, "out"), exist_ok=True)87with open(os.path.join(ROOT, "out", "results.json"), "w") as f:88 json.dump(res, f, indent=2)89print(json.dumps(res, indent=2))90