1"""F5-conf-02: Monte Carlo bias of OLS rho_hat in AR(1) with intercept.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py"""3import json, csv, os4import numpy as np5import pandas as pd67ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))8d = pd.read_csv(os.path.join(ROOT, "data", "design.csv"), dtype=str)9p = dict(zip(d["parameter"].str.strip(), d["value"].str.strip()))10n = int(p["n"]); rho = float(p["rho"]); sigma = float(p["sigma"])11R = int(p["n_replications"]); seed = int(p["seed"])12issues = []13assert p["intercept"].lower() == "yes" and p["initial_condition"].lower() == "stationary"14assert abs(rho) < 11516rng = np.random.default_rng(seed)17# Draw order: y_1 (stationary) for all reps, then innovations e_2..e_n (R x (n-1)).18y = np.empty((R, n))19y[:, 0] = rng.normal(0.0, sigma / np.sqrt(1 - rho**2), size=R)20eps = rng.normal(0.0, sigma, size=(R, n - 1))21for t in range(1, n):22 y[:, t] = rho * y[:, t - 1] + eps[:, t - 1] # c = 02324x = y[:, :-1]; z = y[:, 1:] # regress y_t on (1, y_{t-1}), t=2..n25xd = x - x.mean(axis=1, keepdims=True)26zd = z - z.mean(axis=1, keepdims=True)27rho_hat = (xd * zd).sum(axis=1) / (xd ** 2).sum(axis=1)2829# Cross-check a few replications against lstsq with explicit intercept30for i in range(5):31 X = np.column_stack([np.ones(n - 1), x[i]])32 b = np.linalg.lstsq(X, z[i], rcond=None)[0]33 assert abs(b[1] - rho_hat[i]) < 1e-103435finite = np.isfinite(rho_hat)36rh = rho_hat[finite]37Rk = int(rh.size)38with open(os.path.join(ROOT, "out", "exclusions.csv"), "w", newline="") as f:39 w = csv.writer(f)40 w.writerow(["step", "rule", "n_before", "n_after", "n_removed", "rationale"])41 w.writerow([1, "drop replications with non-finite rho_hat", R, Rk, R - Rk,42 "OLS undefined if regressor has zero variance (probability zero)"])4344mean = float(rh.mean()); sd = float(rh.std(ddof=1))45bias_mc = mean - rho; se = sd / np.sqrt(Rk)46bias_an = -(1 + 3 * rho) / n47res = {48 "task_id": "F5-conf-02",49 "bias_analytical": bias_an,50 "bias_mc": bias_mc,51 "bias_mc_se": float(se),52 "rho_hat_mean": mean,53 "rho_hat_sd_mc": sd,54 "discrepancy_se_units": float((bias_mc - bias_an) / se),55 "se_asymptotic": float(np.sqrt((1 - rho**2) / n)),56 "n_replications": R,57 "feasibility": "complete",58 "feasibility_reason": None,59 "data_issues": issues,60 "n_input_rows": {"design.csv": int(len(d))},61 "n_analysed": Rk,62 "random_seed": seed,63 # supplementary (not required) quantities64 "bias_analytical_neff": -(1 + 3 * rho) / (n - 1),65 "discrepancy_se_units_neff": float((bias_mc + (1 + 3 * rho) / (n - 1)) / se),66 "bias_mc_ci95_low": float(bias_mc - 1.96 * se),67 "bias_mc_ci95_high": float(bias_mc + 1.96 * se),68 "ratio_bias_mc_to_analytical": float(bias_mc / bias_an),69 "rho_hat_median": float(np.median(rh)),70 "share_rho_hat_below_rho": float(np.mean(rh < rho)),71 "se_asymptotic_over_sd_mc": float(np.sqrt((1 - rho**2) / n) / sd),72}73with open(os.path.join(ROOT, "out", "results.json"), "w") as f:74 json.dump(res, f, indent=2)75print(json.dumps(res, indent=2))76