1#!/usr/bin/env python32"""Generates the synthetic 1-minute OHLCV fixtures used by CI (no vendor data) and the reference values3(numpy/statsmodels/exchange_calendars/arch) consumed by tests/reference. Deterministic (seed 20261102).45Fixture: 2 symbols × 5 sessions around Thanksgiving 2024 (2024-11-25 … 2024-11-29, incl. the early close on 11-296and the holiday 11-28), extended hours 04:00–19:59 ET with regular bars 09:30–16:00 inclusive; ~3 % of regular7minutes randomly absent (zero-volume minutes are absent at the source); the 16:00 bar is present on every session.8"""9import json, math, os, sys10import numpy as np11import pandas as pd1213root = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))14out_dir = os.path.join(root, "datasets", "fixtures"); os.makedirs(out_dir, exist_ok=True)15rng = np.random.default_rng(20261102)16sessions = ["2024-11-25", "2024-11-26", "2024-11-27", "2024-11-29", "2024-12-02"]17early = {"2024-11-29"}18rows = []19truth = {}20for sym, p0, vol in [("SYN1", 100.0, 0.0008), ("SYN2", 25.0, 0.0015)]:21 price = p022 truth[sym] = {}23 for s in sessions:24 close_min = 13 * 60 if s in early else 16 * 6025 ext_end = 17 * 60 if s in early else 20 * 6026 minutes = list(range(4 * 60, ext_end))27 day = pd.Timestamp(s, tz="America/New_York")28 sess_rv = 0.0; n_reg = 0; closes = []29 for m in minutes:30 regular = 9 * 60 + 30 <= m <= close_min31 # drop ~3 % of regular minutes and ~40 % of extended minutes, but never the open/close prints32 if m not in (9 * 60 + 30, close_min):33 if regular and rng.random() < 0.03: continue34 if not regular and rng.random() < 0.40: continue35 r = rng.normal(0, vol * (1.0 if regular else 0.5))36 o = price; c = price * math.exp(r); h = max(o, c) * (1 + abs(rng.normal(0, 1e-4))); l = min(o, c) * (1 - abs(rng.normal(0, 1e-4)))37 v = float(rng.integers(100, 5000)) if regular else float(rng.integers(0, 300) + 1)38 ts = (day + pd.Timedelta(minutes=m)).tz_convert("UTC")39 rows.append((sym, ts.strftime("%Y-%m-%dT%H:%M:%SZ"), round(o, 4), round(h, 4), round(l, 4), round(c, 4), v))40 price = c41 if regular:42 closes.append(round(c, 4)); n_reg += 1 # truth from the rounded prices written to the CSV43 lr = np.diff(np.log(np.array(closes)))44 truth[sym][s] = {"n_regular_bars": n_reg, "rv": float(np.sum(lr ** 2)), "bpv": float((math.pi / 2) * (len(lr) / (len(lr) - 1)) * np.sum(np.abs(lr[1:]) * np.abs(lr[:-1]))),45 "rsv_neg": float(np.sum(lr[lr < 0] ** 2)), "rsv_pos": float(np.sum(lr[lr > 0] ** 2)), "rq": float(len(lr) / 3 * np.sum(lr ** 4)),46 "expected_bars": 211 if s in early else 391}47df = pd.DataFrame(rows, columns=["symbol", "timestamp_utc", "open", "high", "low", "close", "volume"])48df.to_csv(os.path.join(out_dir, "synthetic_bars_1m.csv"), index=False)49json.dump(truth, open(os.path.join(out_dir, "synthetic_bars_1m.expected.json"), "w"), indent=1)50print("wrote", len(df), "bars for", df.symbol.nunique(), "symbols;", df.groupby("symbol").size().to_dict())5152# --- reference values from established software -------------------------------------------------------------53ref = {}54try:55 import exchange_calendars as xc56 cal = xc.get_calendar("XNYS", start="2000-01-03", end="2026-12-31")57 sess = cal.sessions_in_range("2000-01-03", "2026-12-31")58 ref["xnys"] = {"n_sessions_2000_2026": int(len(sess)), "sessions_2024": [str(d.date()) for d in cal.sessions_in_range("2024-01-01", "2024-12-31")],59 "early_closes_2000_2026": [str(d.date()) for d in cal.early_closes if "2000-01-03" <= str(d.date()) <= "2026-12-31"],60 "holidays_2000_2026": [str(d.date()) for d in pd.bdate_range("2000-01-03", "2026-12-31") if d not in sess],61 "exchange_calendars_version": xc.__version__}62except Exception as e:63 ref["xnys"] = {"error": str(e)}64import statsmodels.api as sm65from scipy import stats as sps66X = np.column_stack([np.ones(500), rng.normal(size=500), rng.normal(size=500) * 2 + 1]); y = 0.5 + 1.5 * X[:, 1] - 0.7 * X[:, 2] + rng.normal(size=500) * (1 + 0.5 * np.abs(X[:, 1]))67ols = sm.OLS(y, X).fit(); hc0 = sm.OLS(y, X).fit(cov_type="HC0"); hc3 = sm.OLS(y, X).fit(cov_type="HC3"); hac = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": 5, "use_correction": False})68groups = np.arange(500) // 25; cl = sm.OLS(y, X).fit(cov_type="cluster", cov_kwds={"groups": groups})69np.savetxt(os.path.join(out_dir, "ols_design.csv"), np.column_stack([y, X[:, 1], X[:, 2]]), delimiter=",", header="y,x1,x2", comments="", fmt="%.17g")70ref["ols"] = {"params": ols.params.tolist(), "bse": ols.bse.tolist(), "pvalues": ols.pvalues.tolist(), "rsquared": float(ols.rsquared), "ssr": float(ols.ssr),71 "hc0_bse": hc0.bse.tolist(), "hc3_bse": hc3.bse.tolist(), "hac5_bse": hac.bse.tolist(), "cluster_bse": cl.bse.tolist(), "statsmodels_version": sm.__version__,72 "t_pvalues": {"1.96_inf": float(2 * sps.norm.sf(1.96)), "2.0_10": float(2 * sps.t.sf(2.0, 10)), "-3.0_30": float(2 * sps.t.sf(3.0, 30))}}73try:74 from arch import arch_model75 import arch76 s2 = 0.04 / (1 - 0.05 - 0.9); e = 0.0; r = np.empty(3000); g = np.random.default_rng(7)77 for t in range(3000):78 s2 = 0.04 + 0.05 * e * e + 0.9 * s2; e = math.sqrt(s2) * g.standard_normal(); r[t] = 0.01 + e79 np.savetxt(os.path.join(out_dir, "garch_returns.csv"), r, header="r", comments="", fmt="%.17g")80 out = {"arch_version": arch.__version__}81 for name, kw in [("garch", dict(vol="GARCH", p=1, q=1)), ("gjr", dict(vol="GARCH", p=1, o=1, q=1)), ("egarch", dict(vol="EGARCH", p=1, o=1, q=1))]:82 for dist in ["normal", "t"]:83 am = arch_model(r, mean="Constant", dist=dist, rescale=False, **kw); res = am.fit(disp="off")84 out[f"{name}.{dist}"] = {"params": {k: float(v) for k, v in res.params.items()}, "std_err": {k: float(v) for k, v in res.std_err.items()}, "loglik": float(res.loglikelihood),85 "converged": bool(res.convergence_flag == 0), "forecast_h5": res.forecast(horizon=5, reindex=False).variance.values[-1].tolist() if name != "egarch" else res.forecast(horizon=5, method="simulation", simulations=20000, reindex=False).variance.values[-1].tolist()}86 ref["garch"] = out87except Exception as e:88 ref["garch"] = {"error": repr(e)}89json.dump(ref, open(os.path.join(out_dir, "reference_values.json"), "w"), indent=1)90print("reference values written:", list(ref.keys()))9192# --- synthetic Arrow IPC page + Parquet in the FinData bars schema v1 (one regular session, 391 bars) --------------93import pyarrow as pa, pyarrow.ipc as ipc, pyarrow.parquet as pq94day = pd.Timestamp("2024-11-27", tz="America/New_York"); n = 39195ts = [(day + pd.Timedelta(minutes=9 * 60 + 30 + i)).tz_convert("UTC").to_pydatetime() for i in range(n)]96g = np.random.default_rng(3); c = 100 * np.exp(np.cumsum(g.normal(0, 5e-4, n))); o = np.r_[100.0, c[:-1]]97schema = pa.schema([("instrument_id", pa.string()), ("symbol_as_reported", pa.string()), ("venue", pa.string()), ("asset_class", pa.string()), ("interval", pa.string()),98 ("timestamp_utc", pa.timestamp("us", tz="UTC")), ("session_date", pa.date32()), ("open", pa.float64()), ("high", pa.float64()), ("low", pa.float64()), ("close", pa.float64()),99 ("volume", pa.float64()), ("currency", pa.string()), ("trade_count", pa.int32()), ("vwap", pa.float64()), ("adjustment_type", pa.string()), ("source_id", pa.string()),100 ("ingested_at", pa.timestamp("us", tz="UTC")), ("quality_flags", pa.list_(pa.string())), ("revision_id", pa.int32())],101 metadata={"findata.schema": "bars", "findata.schema_version": "1", "findata.timestamp_convention": "bar_open"})102tbl = pa.table({"instrument_id": ["synthetic:stock:SYN1"] * n, "symbol_as_reported": ["SYN1"] * n, "venue": [None] * n, "asset_class": ["stock"] * n, "interval": ["1m"] * n,103 "timestamp_utc": pa.array(ts, pa.timestamp("us", tz="UTC")), "session_date": [day.date()] * n, "open": o, "high": np.maximum(o, c) * 1.0001, "low": np.minimum(o, c) * 0.9999, "close": c,104 "volume": g.integers(100, 1000, n).astype(float), "currency": ["USD"] * n, "trade_count": pa.array([None] * n, pa.int32()), "vwap": pa.array([None] * n, pa.float64()),105 "adjustment_type": ["raw"] * n, "source_id": ["synthetic"] * n, "ingested_at": pa.array([ts[0]] * n, pa.timestamp("us", tz="UTC")), "quality_flags": [[]] * n, "revision_id": pa.array([1] * n, pa.int32())}, schema=schema)106with open(os.path.join(out_dir, "synthetic_bars_page.arrow"), "wb") as fh:107 with ipc.new_stream(fh, schema) as w:108 w.write_table(tbl)109pq.write_table(tbl, os.path.join(out_dir, "synthetic_bars_page.parquet"))110print("wrote synthetic_bars_page.arrow/.parquet", tbl.num_rows, "rows")111112# --- HAR reference: simulated daily RV (log-AR(1)), HAR-RV h=1 and h=5 (average target) with HAC, log-HAR, HARQ -----113g = np.random.default_rng(40); T = 1500; lv = np.empty(T); lv[0] = np.log(1e-4)114for t in range(1, T): lv[t] = 0.3 * np.log(1e-4) + 0.7 * lv[t - 1] + 0.4 * g.standard_normal()115rv = np.exp(lv); rq = rv ** 2 * (1 + 0.5 * g.random(T))116np.savetxt(os.path.join(out_dir, "har_rv.csv"), np.column_stack([rv, rq]), delimiter=",", header="rv,rq", comments="", fmt="%.17g")117def har_design(x, lags=(1, 5, 22), h=1, log=False, harq=None, target="average"):118 z = np.log(x) if log else x; L = max(lags); rows = []; ys = []119 sqrt_rq_mean = np.sqrt(np.mean(np.maximum(harq, 0.0))) if harq is not None else 0.0120 for t in range(L - 1, len(x) - h):121 row = [1.0] + [z[t - l + 1:t + 1].mean() for l in lags]122 if harq is not None: row.append((np.sqrt(max(harq[t], 0.0)) - sqrt_rq_mean) * x[t])123 if target == "average": yv = x[t + 1:t + h + 1].mean(); yv = np.log(yv) if log else yv124 else: yv = z[t + h]125 rows.append(row); ys.append(yv)126 return np.array(rows), np.array(ys)127har = {}128for name, kw in [("h1", dict(h=1)), ("h5", dict(h=5)), ("h22", dict(h=22)), ("log_h1", dict(h=1, log=True)), ("harq_h1", dict(h=1, harq=rq))]:129 X, y = har_design(rv, **kw); n = len(y); h = kw["h"]; L = max(h - 1, int(np.floor(4 * (n / 100) ** (2 / 9))))130 f = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": L, "use_correction": False})131 har[name] = {"params": f.params.tolist(), "hac_bse": f.bse.tolist(), "hac_lags": L, "n": n, "ssr": float(f.ssr), "rsquared": float(f.rsquared)}132ref["har"] = har; json.dump(ref, open(os.path.join(out_dir, "reference_values.json"), "w"), indent=1)133print("HAR references:", {k: v["n"] for k, v in har.items()})134135# --- VAR(2) reference: statsmodels VAR, orthogonalised FEVD, generalized FEVD (Pesaran–Shin) by hand, DY totals, LW ----136from statsmodels.tsa.api import VAR137g = np.random.default_rng(77); N = 4; T = 1200138A1 = np.array([[0.5, 0.1, 0.0, 0.05], [0.0, 0.4, 0.15, 0.0], [0.1, 0.0, 0.3, 0.1], [0.05, 0.05, 0.0, 0.45]]); A2 = 0.15 * np.eye(N)139L = np.linalg.cholesky(np.array([[1.0, 0.3, 0.1, 0.0], [0.3, 1.0, 0.2, 0.1], [0.1, 0.2, 1.0, 0.3], [0.0, 0.1, 0.3, 1.0]]))140Y = np.zeros((T, N))141for t in range(2, T): Y[t] = 0.1 + A1 @ Y[t - 1] + A2 @ Y[t - 2] + L @ g.standard_normal(N)142np.savetxt(os.path.join(out_dir, "var_panel.csv"), Y, delimiter=",", header=",".join(f"y{i}" for i in range(N)), comments="", fmt="%.17g")143res = VAR(Y).fit(2, trend="c"); H = 10144Phi = [np.eye(N)]145coefs = [res.coefs[i] for i in range(2)]146for h in range(1, H): Phi.append(sum(coefs[j - 1] @ Phi[h - j] for j in range(1, min(h, 2) + 1)))147S = res.sigma_u148num = np.zeros((N, N)); den = np.zeros(N)149for h in range(H):150 PS = Phi[h] @ S151 for i in range(N):152 den[i] += PS[i] @ Phi[h][i]153 for j in range(N): num[i, j] += PS[i, j] ** 2 / S[j, j]154th = num / den[:, None]; th = th / th.sum(axis=1, keepdims=True)155dy_total = 100 * (th.sum() - np.trace(th)) / N; dy_from = 100 * (th.sum(axis=1) - np.diag(th)); dy_to = 100 * (th.sum(axis=0) - np.diag(th))156sel = VAR(Y).select_order(4)157try:158 from sklearn.covariance import ledoit_wolf as sk_lw159 lw_cov, lw_s = sk_lw(Y)160except Exception:161 Xc = Y - Y.mean(axis=0); n = Xc.shape[0]; Sc = Xc.T @ Xc / n; mu = np.trace(Sc) / N162 d2 = np.sum((Sc - mu * np.eye(N)) ** 2) / N; b = sum(np.sum((np.outer(x, x) - Sc) ** 2) for x in Xc) / (n * n * N); lw_s = min(1.0, b / d2); lw_cov = (1 - lw_s) * Sc + lw_s * mu * np.eye(N)163ref = json.load(open(os.path.join(out_dir, "reference_values.json")))164ref["var"] = {"nobs": int(res.nobs), "params": res.params.tolist(), "sigma_u": res.sigma_u.tolist(), "llf": float(res.llf), "fevd10": res.fevd(H).decomp[:, -1, :].tolist(),165 "gfevd10": th.tolist(), "dy_total10": float(dy_total), "dy_from10": dy_from.tolist(), "dy_to10": dy_to.tolist(), "bic_order": int(sel.bic), "aic_order": int(sel.aic),166 "lw_shrinkage": float(lw_s), "lw_cov": np.asarray(lw_cov).tolist(), "statsmodels_version": sm.__version__}167json.dump(ref, open(os.path.join(out_dir, "reference_values.json"), "w"), indent=1)168print("VAR reference: nobs", res.nobs, "DY total", round(dy_total, 3), "bic order", int(sel.bic), "LW shrinkage", round(float(lw_s), 4))169