1"""F5-conf-02: Monte Carlo finite-sample bias of OLS AR(1) estimator vs Marriott-Pope/Kendall approximation.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py"""3import json, os4import numpy as np, pandas as pd56d = pd.read_csv("data/design.csv")7n_rows = len(d)8p = dict(zip(d["parameter"].astype(str).str.strip(), d["value"].astype(str).str.strip()))9issues = []10n = int(p["n"]); rho = float(p["rho"]); sigma = float(p["sigma"])11R = int(p["n_replications"]); seed = int(p["seed"])12intercept = p["intercept"].lower() in ("yes", "true", "1")13init = p["initial_condition"].lower()14assert intercept and init == "stationary" and abs(rho) < 1 and sigma > 0 and n > 21516rng = np.random.default_rng(seed)17# Draw all shocks as one (R, n) standard-normal array; column 0 -> y_1 from the stationary law.18z = rng.standard_normal((R, n))19y = np.empty((R, n))20y[:, 0] = z[:, 0] * sigma / np.sqrt(1 - rho**2)21for t in range(1, n):22 y[:, t] = rho * y[:, t - 1] + sigma * z[:, t] # c = 023# OLS of y_t on (1, y_{t-1}), t = 2..n (n-1 usable observations)24x = y[:, :-1]; yy = y[:, 1:]25xc = x - x.mean(axis=1, keepdims=True); yc = yy - yy.mean(axis=1, keepdims=True)26rho_hat = (xc * yc).sum(axis=1) / (xc**2).sum(axis=1)27assert np.all(np.isfinite(rho_hat))2829mean = float(rho_hat.mean()); sd = float(rho_hat.std(ddof=1))30bias_mc = mean - rho; se = sd / np.sqrt(R)31bias_an = -(1 + 3 * rho) / n32res = {33 "task_id": "F5-conf-02",34 "bias_analytical": bias_an,35 "bias_mc": bias_mc,36 "bias_mc_se": se,37 "rho_hat_mean": mean,38 "rho_hat_sd_mc": sd,39 "discrepancy_se_units": (bias_mc - bias_an) / se,40 "se_asymptotic": float(np.sqrt((1 - rho**2) / n)),41 "n_replications": R,42 "feasibility": "complete",43 "feasibility_reason": None,44 "data_issues": issues,45 "n_input_rows": {"design.csv": n_rows},46 "n_analysed": int(rho_hat.size),47 "random_seed": seed,48 # supplementary (not required) diagnostics49 "supplementary": {50 "bias_analytical_using_n_minus_1": -(1 + 3 * rho) / (n - 1),51 "median_rho_hat": float(np.median(rho_hat)),52 "share_rho_hat_below_rho": float((rho_hat < rho).mean()),53 "bias_mc_ci95": [bias_mc - 1.96 * se, bias_mc + 1.96 * se],54 "ratio_bias_mc_to_analytical": bias_mc / bias_an,55 "numpy_version": np.__version__,56 },57}58os.makedirs("out", exist_ok=True)59with open("out/results.json", "w") as f:60 json.dump(res, f, indent=2)61print(json.dumps(res, indent=2))62