"""Single causal pass over origins producing everything the non-neural baselines and the scorer need: realised targets per horizon, trailing sigma, 19 tree features, HAR regressors, EWMA and seasonal variance forecasts, GARCH(1,1)-t variance forecasts (given per-instrument parameters fitted on training data). Parallel over chunks.""" from __future__ import annotations from concurrent.futures import ProcessPoolExecutor import numpy as np import pandas as pd from ..data.windows import OriginIndex, SymbolGrid, N_DAY, trailing_returns, make_example TREE_FEATURES = ["log_rv30", "log_rv390", "log_rv5d", "ret30", "ret390", "ret_day", "overnight", "relvol30", "range30", "minute_frac", "time_to_close", "weekday", "early_close", "log_sigma", "skew390", "kurt390", "sign_last", "ret5", "rv_ratio_30_390"] HAR_FEATURES = ["log_rv30", "log_rv390", "log_rv5d"] EWMA_LAMBDA = 0.97 def _five_minute_hist(g: SymbolGrid, s: int, m: int, n_days: int = 3) -> np.ndarray: jB = (m + 1) // 5 - 1 out = [] for d in range(max(0, s - n_days), s + 1): last = jB if d == s else 77 if last < 1: continue c = np.asarray(g.flat_b5[d * 78: d * 78 + last + 1, 3]); mm = np.asarray(g.flat_m5[d * 78: d * 78 + last + 1]) > 0 ix = np.flatnonzero(mm) if len(ix) >= 2: out.append(np.diff(c[ix])) return np.concatenate(out) if out else np.zeros(0) def garch_forecast(params: dict, r5: np.ndarray, H: int) -> float: om, al, be = params["omega"], params["alpha"], params["beta"] ab = al + be uncond = om / max(1e-9, 1 - ab) if ab < 0.999 else float(np.var(r5)) if len(r5) else om h = uncond for r in r5[-400:]: h = om + al * r * r + be * h steps = max(1, int(np.ceil(H / 5))) k = np.arange(steps) var = float((uncond + ab ** k * (h - uncond)).sum()) return var * (H / (5 * steps)) def origin_row(g: SymbolGrid, s: int, m: int, horizons: list[int], garch: dict | None) -> dict: ex = make_example(g, s, m, type("W", (), {"L_A": 390, "L_B": 390, "L_C": 260, "H": 60})()) sigma = float(ex["sigma"]); y, ym = ex["y"], ex["y_mask"] > 0 row = {"sym": None, "symbol": g.symbol, "s": s, "m": m, "sigma": sigma, "session_date": g.meta["session_dates"][s], "asset_class": g.meta["asset_class"]} for H in horizons: row[f"cum_{H}"] = float((y[:H] * ym[:H]).sum()); row[f"rv_{H}"] = float(((y[:H] * ym[:H]) ** 2).sum()) row[f"present_{H}"] = float(ym[:H].mean()) r = trailing_returns(g, s, m, 6 * N_DAY) # up to 6 sessions of same-session returns r30, r390, r5 = r[-30:], r[-390:], r[-5:] rv30 = float(np.mean(r30 ** 2)) if len(r30) else sigma ** 2 rv390 = float(np.mean(r390 ** 2)) if len(r390) else sigma ** 2 rv5d = float(np.mean(r ** 2)) if len(r) else sigma ** 2 lc_today = np.asarray(g.flat_lc[s * N_DAY: s * N_DAY + m + 1], dtype=np.float64); ok = ~np.isnan(lc_today) ret_day = float(lc_today[ok][-1] - lc_today[ok][0]) if ok.sum() >= 2 else 0.0 gap = float(g.gap[s]); gap = 0.0 if np.isnan(gap) else gap sl = slice(s * N_DAY + max(0, m - 29), s * N_DAY + m + 1) base = np.asarray(g.flat_base[sl]); mk = np.asarray(g.flat_mask[sl]) > 0; seas = np.asarray(g.flat_seas[sl]) rel = base[mk, 4] - seas[mk]; rel = rel[~np.isnan(rel)] sd390 = float(r390.std()) + 1e-12 if len(r390) > 1 else 1e-12 z390 = (r390 - r390.mean()) / sd390 if len(r390) > 1 else np.zeros(1) row.update({"log_rv30": np.log(rv30 + 1e-12), "log_rv390": np.log(rv390 + 1e-12), "log_rv5d": np.log(rv5d + 1e-12), "ret30": r30.sum() / sigma, "ret390": r390.sum() / sigma, "ret_day": ret_day / sigma, "overnight": gap / sigma, "relvol30": float(rel.mean()) if len(rel) else 0.0, "range30": float((base[mk, 1] - base[mk, 2]).mean()) / sigma if mk.any() else 0.0, "minute_frac": m / 389.0, "time_to_close": (g.n_slots[s] - 1 - m) / 389.0, "weekday": g.weekday[s] / 4.0, "early_close": float(g.early[s]), "log_sigma": np.log(sigma), "skew390": float(np.mean(z390 ** 3)), "kurt390": float(np.mean(z390 ** 4)), "sign_last": float(np.sign(r[-1])) if len(r) else 0.0, "ret5": r5.sum() / sigma, "rv_ratio_30_390": np.log((rv30 + 1e-12) / (rv390 + 1e-12))}) # EWMA (minute scale) variance per minute rr = r[-3 * N_DAY:] if len(rr) >= 30: w = EWMA_LAMBDA ** np.arange(len(rr))[::-1]; v1 = float((w * rr * rr).sum() / w.sum()) else: v1 = sigma ** 2 row["ewma_var1"] = v1 # seasonal profile (past-only squared-return profile for the future slots), level-adjusted by trailing 30-min RV prof = np.asarray(g.flat_seas_r2[s * N_DAY: (s + 1) * N_DAY], dtype=np.float64) prof = np.where(np.isnan(prof), np.nan, prof) ref = np.nanmean(prof[max(0, m - 29): m + 1]) if np.isfinite(prof[max(0, m - 29): m + 1]).any() else np.nan level = float(np.clip(rv30 / ref, 0.25, 4.0)) if (ref is not None and np.isfinite(ref) and ref > 0) else 1.0 for H in horizons: fut = prof[m + 1: m + 1 + H] if np.isfinite(fut).sum() >= H // 2: row[f"seas_var_{H}"] = float(np.nansum(np.where(np.isfinite(fut), fut, np.nanmean(fut))) * level) else: row[f"seas_var_{H}"] = sigma ** 2 * H if garch is not None: row[f"garch_var_{H}"] = garch_forecast(garch, _five_minute_hist(g, s, m), H) if garch is not None: row["garch_nu"] = garch["nu"] return row def _chunk(args) -> list[dict]: index, ks, horizons, garch_params = args rows = [] for k in ks: i, s, m = index.rows[k] r = origin_row(index.grid(int(i)), int(s), int(m), horizons, garch_params.get(int(i)) if garch_params else None) r["sym"] = int(i); r["k"] = int(k) rows.append(r) return rows def origin_frame(index: OriginIndex, horizons: list[int], garch_params: dict | None = None, workers: int = 8, chunk: int = 2000) -> pd.DataFrame: ks = np.arange(len(index)) chunks = [(index, ks[i:i + chunk], horizons, garch_params or {}) for i in range(0, len(ks), chunk)] rows = [] if workers <= 1 or len(chunks) == 1: for c in chunks: rows.extend(_chunk(c)) else: with ProcessPoolExecutor(max_workers=workers) as ex: for part in ex.map(_chunk, chunks): rows.extend(part) df = pd.DataFrame(rows).sort_values("k").reset_index(drop=True) return df