1"""F5-conf-02: Monte Carlo bias of OLS AR(1) coefficient vs Marriott-Pope/Kendall approximation.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py"""3import json, os4import numpy as np5import pandas as pd67d = pd.read_csv("data/design.csv")8p = dict(zip(d["parameter"], d["value"]))9n = int(p["n"]); rho = float(p["rho"]); sigma = float(p["sigma"])10R = int(p["n_replications"]); seed = int(p["seed"])11assert str(p["intercept"]).strip().lower() == "yes"12assert str(p["initial_condition"]).strip().lower() == "stationary"13c = 0.01415rng = np.random.default_rng(seed)16# One (R, n) standard-normal draw: column 0 -> stationary y_1, columns 1..n-1 -> innovations.17z = rng.standard_normal((R, n))18y = np.empty((R, n))19y[:, 0] = z[:, 0] * sigma / np.sqrt(1 - rho**2)20eps = sigma * z[:, 1:]21for t in range(1, n):22 y[:, t] = c + rho * y[:, t - 1] + eps[:, t - 1]2324# OLS of y_t on (1, y_{t-1}), t = 2..n -> slope = cov/var (demeaned within each sample)25x = y[:, :-1]; yy = y[:, 1:]26xd = x - x.mean(axis=1, keepdims=True); yd = yy - yy.mean(axis=1, keepdims=True)27rho_hat = (xd * yd).sum(axis=1) / (xd**2).sum(axis=1)2829ok = np.isfinite(rho_hat)30rho_hat = rho_hat[ok]31m = float(rho_hat.mean()); sd = float(rho_hat.std(ddof=1))32bias_mc = m - rho; se = sd / np.sqrt(len(rho_hat))33bias_an = -(1 + 3 * rho) / n34res = {35 "task_id": "F5-conf-02",36 "bias_analytical": bias_an,37 "bias_mc": bias_mc,38 "bias_mc_se": se,39 "rho_hat_mean": m,40 "rho_hat_sd_mc": sd,41 "discrepancy_se_units": (bias_mc - bias_an) / se,42 "se_asymptotic": float(np.sqrt((1 - rho**2) / n)),43 "n_replications": R,44 "feasibility": "complete",45 "feasibility_reason": None,46 "data_issues": [],47 "n_input_rows": {"design.csv": int(len(d))},48 "n_analysed": int(ok.sum()),49 "random_seed": seed,50 # supplementary51 "bias_analytical_n_minus_1": -(1 + 3 * rho) / (n - 1),52 "bias_mc_ci95": [bias_mc - 1.96 * se, bias_mc + 1.96 * se],53 "rho_hat_median": float(np.median(rho_hat)),54 "share_rho_hat_below_rho": float((rho_hat < rho).mean()),55}56os.makedirs("out", exist_ok=True)57with open("out/results.json", "w") as f:58 json.dump(res, f, indent=2)59print(json.dumps(res, indent=2))60