Building Causal HAR, Spectral, and Path Signature Features

# 04_model_based_features.py


```py
# ---
# jupyter:
#   jupytext:
#     cell_metadata_filter: tags,-all
#     text_representation:
#       extension: .py
#       format_name: percent
#       format_version: '1.3'
#       jupytext_version: 1.19.3
#   kernelspec:
#     display_name: Python 3 (ipykernel)
#     language: python
#     name: python3
# ---

# %% [markdown]
# # NASDAQ-100 microstructure: features a model has to be fitted to produce
#
# **Chapter 9: Model-Based Feature Extraction**
#
# The previous stage built features that are arithmetic on past bars - a spread, a realized
# variance, a share of volume. This one builds features that only exist once a model has been
# estimated: the number a reader gets depends on parameters fitted from data, so the window
# those parameters came from is part of what the feature knows. That is the whole subject of
# this notebook, and Section A is where it is argued.
#
# Three procedures run on the minute panel, in increasing distance from ordinary arithmetic:
#
# | Procedure | What it produces | What is estimated |
# |---|---|---|
# | HAR regression | a forecast of the next few minutes' variance, and the error in the last one | three regression coefficients, refitted every bar |
# | Rolling Fourier transform | how much of the recent activity sits at which frequency | nothing; a fixed transform of the window |
# | Depth-2 path signature | the order in which price, order flow and trade intensity moved | nothing; a fixed transform of the window |
#
# **What you will be able to do after working through it**
#
# - Split a volatility forecast into components measured over different lengths of history, and
#   refit it as time passes so that no coefficient is ever estimated from a bar the forecast is
#   supposed to precede.
# - Turn a rolling window of volume into a small set of numbers describing how repetitive the
#   recent activity has been, and read those numbers back.
# - Summarise a short stretch of price and order flow by which one moved first, in a form a
#   model can use as a column.
# - Write out a feature table keyed by bar rather than by fold, and check that it covers every
#   window every configured prediction target will ask for.
# - Measure whether any of it ranks the cross-section, with a standard error that accounts for
#   the dependence between neighbouring observations of a time series.
#
# **What it reads and what it writes**
#
# - Reads the AlgoSeek minute archive through `load_nasdaq100_bars`, and the label files under
#   `labels/` for their timestamps alone - those decide the walk-forward windows.
# - Writes `features/model_based.parquet`, one row per (`timestamp`, `symbol`) and no `fold`
#   column - Section E says why - with a `.digest.json` sidecar beside it recording what was
#   written and what it was built from.
#
# **Prerequisites**: [`02_labels`](02_labels.ipynb) must have run, because its output defines the
# windows here. [`03_financial_features`](03_financial_features.ipynb) need not have run; its
# output is joined against for a coverage count and is never read for a value.
#
# **The same methods, taught one at a time**:
# [`09_har_rough_volatility`](../../09_model_based_features/09_har_rough_volatility.ipynb),
# [`05_spectral_features`](../../09_model_based_features/05_spectral_features.ipynb),
# [`06_path_signatures`](../../09_model_based_features/06_path_signatures.ipynb).

# %%
"""NASDAQ-100 Microstructure: Model-Based Features (Ch9)."""

import warnings
from pathlib import Path

import numpy as np
import pandas as pd
import plotly.graph_objects as go
import polars as pl
import yaml
from IPython.display import display
from ml4t.diagnostic.evaluation.stats import benjamini_hochberg_fdr
from ml4t.diagnostic.metrics import compute_ic_hac_stats, cross_sectional_ic_series
from ml4t.diagnostic.splitters.calendar import TradingCalendar

from case_studies.utils.artifact_digest import value_digest
from case_studies.utils.temporal import walk_forward_feature, write_model_based
from data import load_nasdaq100_bars
from utils.artifact_specs import resolve_label_buffer, resolve_label_horizon
from utils.cv_splits import generate_cv_splits, load_evaluation_config
from utils.data_quality import top_entities
from utils.paths import get_case_study_dir
from utils.style import COLORS, show_plotly_with_alt

warnings.filterwarnings("ignore")

# %% tags=["parameters"]
CASE_STUDY_ID = "nasdaq100_microstructure"
START_DATE = "2020-01-01"
END_DATE = "2021-12-31"
MAX_SYMBOLS = 0
# Where the artifact is written. Empty is the production path: beside the case study's other
# features. A reduced run sets it to a throwaway directory so that a five-symbol preview cannot
# replace the 115-symbol artifact every registered training run is pinned to by file hash.
# Reads are never redirected - a preview resolves its labels and its panel from the committed
# artifacts, so it is a run of this notebook rather than of a smaller case study.
WORKSPACE = ""

# %% [markdown]
# ## Configuration
#
# Everything the run depends on is bound from `config/setup.yaml` rather than typed here, so a
# change to the case study's declared window or label set reaches this notebook without an edit.
# The estimation windows below are this notebook's own model specification and come from the same
# file, under `model_based`: a window typed at a call site is a number a reader has to find in the
# code to know what was fitted.

# %%
CASE_DIR = get_case_study_dir(CASE_STUDY_ID)
FEATURES_DIR = CASE_DIR / "features"
LABELS_DIR = CASE_DIR / "labels"
OUTPUT_DIR = Path(WORKSPACE) / "features" if WORKSPACE else FEATURES_DIR
OUTPUT_DIR.mkdir(parents=True, exist_ok=True)

SETUP = yaml.safe_load((CASE_DIR / "config" / "setup.yaml").read_text())
EVAL_CFG = load_evaluation_config(CASE_STUDY_ID)
PRIMARY_LABEL = SETUP["labels"]["primary"]
CONFIGURED_LABELS = [PRIMARY_LABEL, *SETUP["labels"].get("variants", [])]
UNIVERSE = sorted(SETUP["universe"]["symbols"])
CALENDAR = EVAL_CFG["calendar"]
HOLDOUT_START = pd.Timestamp(EVAL_CFG["holdout_start"])
HOLDOUT_END = pd.Timestamp(EVAL_CFG["holdout_end"])
# The config states the holdout's last date; parsed as a timestamp it is that date at midnight,
# which is before every intraday bar of the session it names. Anything comparing a bar against the
# end of the holdout uses the exclusive bound instead.
HOLDOUT_END_EXCLUSIVE = HOLDOUT_END + pd.Timedelta(days=1)

# The bar is one minute, so a horizon written as a duration converts to a bar count exactly.
BAR = pd.Timedelta(minutes=1)
# The horizon the label measures, not the buffer that purges it. The buffer is one bar
# wider, because the label reads a quote one bar past its own horizon, and sampling IC
# on the buffer would space observations by an interval no label spans.
LABEL_HORIZON = SETUP["labels"]["horizons"][PRIMARY_LABEL]
LABEL_HORIZON_BARS = int(pd.Timedelta(LABEL_HORIZON) // BAR)
IC_SAMPLE_STEP = LABEL_HORIZON_BARS

# The model specification, read from `setup.yaml::model_based`. The three lengths of history the
# HAR regression averages squared returns over, the history each of its refits reads, how often it
# refits, the floor on a refit's usable observations, and the window each spectrum and each
# signature path spans.
MODEL_BASED = SETUP["model_based"]
HAR_COMPONENTS = tuple(MODEL_BASED["har"]["components"])
HAR_FIT_WINDOW = MODEL_BASED["har"]["fit_window"]
HAR_REFIT_EVERY = MODEL_BASED["har"]["refit_every"]
HAR_MIN_TRAIN_OBS = MODEL_BASED["har"]["min_train_obs"]
# The refit window is the burn-in: the walk's first estimate reads `HAR_FIT_WINDOW` bars and the
# bar after them is the first one any parameters speak for.
HAR_BURNIN = HAR_FIT_WINDOW + 1
FFT_WINDOW = MODEL_BASED["spectrum"]["window"]
FFT_LOW_FREQ_PERIOD = MODEL_BASED["spectrum"]["low_frequency_period"]
SIG_WINDOW = MODEL_BASED["signature"]["window"]

# The exchange's regular session. Bars outside it are not part of any decision.
OPEN_HOUR, OPEN_MINUTE, CLOSE_HOUR = 9, 30, 16

# Every table below is short enough to read whole, and a frame shown as ten rows and an
# ellipsis is a table the reader has to take on trust.
pl.Config.set_tbl_rows(40)

# %%
print(f"Sample: {START_DATE} to {END_DATE}, minute bars on the {CALENDAR} calendar.")
print(
    f"The holdout runs {HOLDOUT_START.date()} to {HOLDOUT_END.date()}. Nothing is fitted on it "
    "and no number printed below is measured over it; features are still written for it, from "
    "windows that end before it, so that a model scored there has inputs."
)
print(
    f"Predictions are made for {PRIMARY_LABEL}, the return over the next {LABEL_HORIZON} "
    f"({LABEL_HORIZON_BARS} bars). The case study also configures "
    f"{', '.join(CONFIGURED_LABELS[1:])}, whose horizons differ, and each of them asks for a "
    "different walk-forward split - which is why Section B resolves one per label."
)
print(
    f"The HAR regression reads {HAR_COMPONENTS[0]}, {HAR_COMPONENTS[1]} and "
    f"{HAR_COMPONENTS[2]} minutes of squared returns and is refitted on the trailing "
    f"{HAR_FIT_WINDOW} bars: long enough for the four coefficients to be identified, short "
    "enough that the fit follows the day rather than the quarter."
)
print(
    f"Each spectrum spans {FFT_WINDOW} bars, which resolves periods up to an hour, and each "
    f"signature path spans {SIG_WINDOW} bars, which is the horizon over which order flow and "
    "price are expected to lead one another."
)
if MAX_SYMBOLS:
    print(f"Universe limited to the {MAX_SYMBOLS} symbols with the most bars.")

# %% [markdown]
# ## The panel this notebook reads
#
# The archive is one row per symbol and minute, carrying the closing quote on each side, the
# traded volume split by where each trade printed against the prevailing quote, and the volume
# reported away from the exchanges. Only the columns the three procedures consume are kept; the
# rest of the archive's sixty-odd columns would triple the memory this notebook holds and
# nothing here reads them.

# %%
READ = [
    "timestamp",
    "symbol",
    "close_bid_price",
    "close_ask_price",
    "volume",
    "finra_volume",
    "total_trades",
    "trade_at_bid",
    "trade_at_bid_mid",
    "trade_at_mid_ask",
    "trade_at_ask",
]

# `lazy=True` is what makes the projection above worth writing. `include_microstructure`
# returns the raw archive schema with no projection of its own, so collecting first and
# selecting afterwards reads all sixty columns into memory to keep eleven. Deferring the
# collect pushes the projection, the universe filter and the session-hours filter into the
# parquet scan, so only the columns and rows this notebook reads are ever materialized.
_lf = load_nasdaq100_bars(
    start_date=START_DATE,
    end_date=str(END_DATE),
    include_microstructure=True,
    symbols=UNIVERSE,
    lazy=True,
).select(READ)

if MAX_SYMBOLS:
    # `top_entities` and not a local sort. Every symbol on this panel quotes on the same
    # padded minute grid, so the row counts a reduction sorts on are equal for every name
    # that traded the whole window - and a descending sort over a group of equal counts
    # returns whatever order the group-by produced, which polars does not fix and which is
    # not stable between two runs of this notebook. Two runs of the same reduced
    # configuration measured here on 2026-09-05 chose {AAPL, MSFT, TSLA} and {AMZN, FB,
    # GOOG}. The shared rule breaks the tie on the symbol name, which is what makes a
    # reduced 04 read the same universe that the reduced 02, 03 and 05 wrote.
    top_syms = top_entities(_lf, MAX_SYMBOLS)
    _lf = _lf.filter(pl.col("symbol").is_in(pl.Series("symbol", top_syms).implode()))
    print(f"Restricted to {MAX_SYMBOLS} symbols: {top_syms}")

_hour, _minute = pl.col("timestamp").dt.hour(), pl.col("timestamp").dt.minute()
df = (
    _lf.filter(
        ((_hour > OPEN_HOUR) | ((_hour == OPEN_HOUR) & (_minute >= OPEN_MINUTE)))
        & (_hour < CLOSE_HOUR)
    )
    .with_columns(pl.col("timestamp").dt.date().alias("session_date"))
    .collect()
)
del _lf

# %% [markdown]
# The vendor emits a padded 390-bar grid on every date, including the afternoons the exchange
# closes early. A bar after the close carries the last quote forward and no
# position could have been taken on it, so a feature computed for that minute is a feature for
# a time at which no decision existed. The session's real length comes from the exchange
# calendar, and the padding is dropped before anything is built, so nothing below describes a
# minute at which the exchange was not open.

# %%
_schedule = TradingCalendar(CALENDAR).calendar.schedule(start_date=START_DATE, end_date=END_DATE)
sessions = pl.DataFrame(
    {
        "session_date": [d.date() for d in _schedule.index],
        "session_bars": (
            (_schedule["market_close"] - _schedule["market_open"]).dt.total_seconds() // 60
        ).astype("int32"),
    }
)
SHORT = sessions.filter(pl.col("session_bars") < sessions["session_bars"].max())
_unscheduled = set(df["session_date"].unique()) - set(sessions["session_date"])
assert not _unscheduled, f"{len(_unscheduled)} session dates are not on the {CALENDAR} calendar"

_minute_of_day = (
    pl.col("timestamp").dt.hour().cast(pl.Int32) * 60
    + pl.col("timestamp").dt.minute().cast(pl.Int32)
    - (OPEN_HOUR * 60 + OPEN_MINUTE)
)
_padded = df.height
df = (
    df.join(sessions, on="session_date", how="inner")
    .filter(_minute_of_day < pl.col("session_bars"))
    .drop("session_bars")
    .sort(["symbol", "timestamp"])
)
print(f"{sessions.height} scheduled sessions, {SHORT.height} of them early closes")
print(f"{_padded - df.height:,} padded bars dropped past the scheduled close")

# %% [markdown]
# What is in the panel, before anything is computed from it: how many names, how much history,
# and how much of a session an average name actually trades through.

# %%
sessions_seen = df.group_by("session_date").agg(pl.col("symbol").n_unique().alias("symbols"))
symbol_sessions = df.group_by("session_date", "symbol").agg(pl.len().alias("bars"))
display(
    pl.DataFrame(
        {
            "quantity": [
                "symbols",
                "sessions",
                "minute bars",
                "first bar",
                "last bar",
                "median bars a name trades in a session",
                "median names quoting per session",
            ],
            "value": [
                f"{df['symbol'].n_unique():,}",
                f"{df['session_date'].n_unique():,}",
                f"{df.height:,}",
                str(df["timestamp"].min()),
                str(df["timestamp"].max()),
                f"{symbol_sessions['bars'].median():,.0f}",
                f"{sessions_seen['symbols'].median():,.0f}",
            ],
        }
    )
)

# %% [markdown]
# ### The series the three procedures read
#
# The **mid price** is the average of the two sides of the quote, which is the price neither
# side of the market has paid a spread to reach. Its one-minute log change is the return series
# every volatility quantity below is built from, and it restarts at each session open so that
# the overnight move never enters a one-minute return.
#
# **Signed volume** is the volume that printed on the ask side minus the volume that printed on
# the bid side: positive when buyers were the ones crossing the spread. Divided by the volume it
# is counted over it becomes a share between -1 and 1, comparable across a heavily traded name
# and a thin one. That denominator has to be the traded volume on **both** venues, because the
# trade-location buckets the numerator comes from count every trade in the bar, including the
# ones reported to the FINRA trade reporting facility rather than to an exchange; `volume` alone
# counts the exchange prints and would make the share exceed one whenever off-exchange activity
# was heavy. The assertion below is what keeps that true rather than assumed.
#
# The **trade count** is the third dimension of the signature path in Section C.3, taken from
# the archive unchanged: how many separate trades made up the bar's volume, which distinguishes
# one large print from a hundred small ones.

# %%
TRADED_VOLUME = (pl.col("volume") + pl.col("finra_volume")).clip(lower_bound=1)
signed_vol = (pl.col("trade_at_ask") + pl.col("trade_at_mid_ask")) - (
    pl.col("trade_at_bid") + pl.col("trade_at_bid_mid")
)
mid = (pl.col("close_bid_price") + pl.col("close_ask_price")) / 2

df = df.with_columns(mid_close=mid)
df = df.filter(pl.col("mid_close").is_not_null() & (pl.col("mid_close") > 0))

group_cols = ["symbol", "session_date"]
df = df.with_columns(
    r1m=(pl.col("mid_close").log() - pl.col("mid_close").log().shift(1).over(group_cols)),
    signed_vol=signed_vol,
    signed_vol_share=(signed_vol / TRADED_VOLUME),
    bar_of_day=pl.col("timestamp").rank("ordinal").over(group_cols).cast(pl.Int32) - 1,
)

_worst = df["signed_vol_share"].abs().max()
assert _worst <= 1.0 + 1e-9, (
    f"signed volume share reaches {_worst:.3f}: the denominator does not cover the numerator"
)
print(f"Signed volume share stays inside [-1, 1]; the largest magnitude is {_worst:.3f}.")

# The digest of the panel as consumed, taken here so that it describes exactly the rows the
# three procedures read. It goes into the artifact's sidecar at the end of Section E.
RAW_DIGEST = value_digest(df.select(READ))
print(f"Minute panel digest: {RAW_DIGEST}")

# %% [markdown]
# ## A. Why a fitted feature is different
#
# A financial feature is a function of past bars. Written down, it is arithmetic: take the last
# thirty midpoints, take their standard deviation, that is the number. Two people with the same
# thirty bars get the same answer, and the answer for 10:45 does not change when 10:46 arrives.
#
# A model-based feature is a function of *parameters that were estimated from* bars. The HAR
# forecast for 10:45 is a weighted sum of three realized variances, and the weights came from a
# regression on some stretch of history. Change that stretch and the forecast changes, even
# though the three variances did not. So the feature's information set is not just the window it
# reads - it is that window **plus** every bar the parameters were estimated from.
#
# This is where look-ahead gets into a feature block without anyone writing anything obviously
# wrong. Fit the regression once on the whole sample and the weights carry information from the
# end of the sample into a forecast made at the beginning; every forecast is then partly a
# summary of what happened afterwards. The forecast will look good, and the backtest built on it
# will look better, and neither result was available to anyone at the time.
#
# The discipline that removes it is to make the estimation window part of the feature's
# definition and then keep that window behind the bar being described. Here it is kept behind by
# construction: the HAR is refitted at every bar on the immediately preceding stretch, so its
# weights at 10:45 were estimated from bars ending at 10:44 and the question of which fold the
# fit belonged to does not arise. The Fourier transform and the path signature estimate nothing
# at all - they are fixed transforms of a trailing window, and they are here because a reader
# who has met the hazard on the HAR should see what the same discipline costs when there is no
# parameter to place.
#
# Two consequences run through the rest of the notebook. First, the feature values do not depend
# on the walk-forward split at all, which is why the artifact written in Section E carries no
# fold column: there is nothing for a fold tag to distinguish, and the label being fitted selects
# its own rows by its own boundaries. A case study whose model is fitted once per fold cannot do
# that, and Section E says what changes there. Second, nothing protects a *diagnostic* the same
# way:
# a number printed about the features is as capable of reading the holdout as a fitted parameter
# is. That is why the folds are resolved next, before anything is computed or printed.

# %% [markdown]
# ## B. The fold contract
#
# A walk-forward split cuts the history into a training window and the validation window that
# follows it, with a gap between them wide enough that the outcome of the last training decision
# has already been realized before the validation window opens. The width of that gap is the
# horizon of the thing being predicted, so **each prediction target gets its own split**: a
# five-minute return seals five minutes and a sixty-minute return seals an hour, and the two
# disagree about where the training window ends and where validation runs to.
#
# This case study configures more than one target and their horizons differ, so this notebook
# resolves one split per target rather than one for the notebook. Each is derived from that
# target's own label file, which is the same
# frame `load_modeling_dataset` uses downstream: fold boundaries are positions in a timestamp
# index, so deriving them from a different frame - the price panel, say, or a feature frame with
# its warm-up rows removed - moves every boundary by however many timestamps the two indexes
# differ by, and the artifact then answers a question about folds nobody downstream is asking.

# %%
label_splits: dict[str, list[dict]] = {}
label_timeline_digest: dict[str, str] = {}
for label in CONFIGURED_LABELS:
    label_path = LABELS_DIR / f"{label}.parquet"
    if not label_path.exists():
        raise FileNotFoundError(f"{label} is configured but not built - run 02_labels.py first.")
    label_ts = pl.scan_parquet(label_path).select("timestamp").unique().collect()
    label_timeline_digest[label] = value_digest(label_ts)
    label_splits[label] = generate_cv_splits(
        label_ts,
        case_study_id=CASE_STUDY_ID,
        label_buffer=resolve_label_buffer(CASE_STUDY_ID, label, SETUP),
        outcome_horizon=resolve_label_horizon(CASE_STUDY_ID, label, SETUP),
        date_col="timestamp",
    )

splits = label_splits[PRIMARY_LABEL]
N_FOLDS = len(splits)
assert all(len(s) == N_FOLDS for s in label_splits.values()), (
    "the configured labels do not agree on how many folds there are"
)

display(
    pl.DataFrame(
        [
            {
                "label": label,
                "seals": resolve_label_buffer(CASE_STUDY_ID, label, SETUP),
                "fold": s["fold"],
                "train_start": s["train_start"],
                "train_end": s["train_end"],
                "val_start": s["val_start"],
                "val_end": s["val_end"],
            }
            for label, label_split in label_splits.items()
            for s in label_split
        ]
    ).sort(["fold", "label"])
)

# %% [markdown]
# The next cell executes the contract rather than describing it. The first check is that no
# training window runs into the validation window it is scored against. The second is the one
# that binds a supervised quantity: a validation bar at $t$ carries an outcome that resolves at
# $t + h$, so the last validation bar a target may use is $h$ before the holdout opens, not the
# bar before it. Both are checked for every configured target, because a split that holds for
# the fifteen-minute return can fail for the sixty-minute one.

# %%
for label, label_split in label_splits.items():
    seal = pd.Timedelta(resolve_label_buffer(CASE_STUDY_ID, label, SETUP))
    for s in label_split:
        assert pd.Timestamp(s["train_end"]) < pd.Timestamp(s["val_start"]), (
            f"{label} fold {s['fold']}: training window runs into its own validation window"
        )
        assert pd.Timestamp(s["val_end"]) + seal <= HOLDOUT_START, (
            f"{label} fold {s['fold']}: a validation outcome resolves inside the holdout"
        )
print(f"The contract holds for {N_FOLDS} folds on each of {len(label_splits)} targets.")

# %% [markdown]
# The four targets disagree about where each fold starts and ends, by minutes, because their
# horizons differ by minutes. The span a fold needs covered therefore runs from the earliest
# training start to the latest validation end across the targets, and minutes are exactly what
# a coverage check downstream is counting.
#
# These spans are what the artifact has to reach, not a tag it carries. Section E writes one row
# per bar with no fold column, so nothing here decides which rows a model reads - the boundaries
# of the label being fitted do. What the spans are used for is the coverage check in Section E
# and the figure below.

# %%
fold_window = {
    s["fold"]: (
        min(
            pd.Timestamp(x["train_start"])
            for sp in label_splits.values()
            for x in sp
            if x["fold"] == s["fold"]
        ),
        max(
            pd.Timestamp(x["val_end"])
            for sp in label_splits.values()
            for x in sp
            if x["fold"] == s["fold"]
        ),
    )
    for s in splits
}
for fold, (start, end) in sorted(fold_window.items()):
    print(f"  Fold {fold} needs {start} .. {end}")
print(
    f"  Fold {N_FOLDS} needs every bar from {min(w[0] for w in fold_window.values())} to "
    f"{HOLDOUT_END.date()}, and is trained on everything before the holdout opens."
)


# %%
def validation_rows(frame: pl.DataFrame) -> pl.DataFrame:
    """Restrict a frame to the validation windows of the primary target's folds.

    Every diagnostic in this notebook goes through this function. The feature frame carries no
    fold column of its own, so a readout built straight from it spans whatever the frame spans,
    holdout included. The primary target's windows are the right ones here because that is the
    target the readouts in Sections C, D and F are about.
    """
    return pl.concat(
        [
            frame.filter(
                (pl.col("timestamp") >= pd.Timestamp(s["val_start"]))
                & (pl.col("timestamp") <= pd.Timestamp(s["val_end"]))
            ).with_columns(pl.lit(s["fold"], dtype=pl.Int32).alias("fold"))
            for s in splits
        ]
    )


# %% [markdown]
# **Figure F1** draws the geometry the artifact has to cover. Each fold is a training span and
# the validation span that follows it; the top row is the holdout, whose training bars all lie
# before the holdout opens so that a model scored on the holdout has features for it without any
# of them having been built from it. The point to read off the figure is that no bar of any
# training span lies to the right of the rule.
#
# The artifact spans the union of everything drawn here, in one row per bar. The figure is a
# picture of what will be asked of it, not of how it is laid out.

# %%
# Sorted by fold id, so the row order is a property of this cell rather than of the order
# `generate_cv_splits` happens to return. Plotly lays a categorical axis out in order of first
# appearance, bottom upwards, so without the sort the bottom row is whichever fold the splits
# list puts first - and that is exactly what the fold renumbering changes. The rendered figure
# is unchanged today; what changes is that the description below stays true when the numbering
# moves.
spans = [
    (f"Fold {s['fold']}", kind, pd.Timestamp(s[f"{key}_start"]), pd.Timestamp(s[f"{key}_end"]))
    for s in sorted(splits, key=lambda s: s["fold"])
    for kind, key in (("Training bars", "train"), ("Validation bars", "val"))
]
spans += [
    (
        f"Fold {N_FOLDS}",
        "Training bars",
        min(w[0] for w in fold_window.values()),
        HOLDOUT_START,
    ),
    (f"Fold {N_FOLDS}", "Holdout bars", HOLDOUT_START, HOLDOUT_END_EXCLUSIVE),
]
span_colors = {
    "Training bars": COLORS["blue"],
    "Validation bars": COLORS["amber"],
    "Holdout bars": COLORS["recede"],
}

fig = go.Figure()
seen = set()
for row, kind, start, end in spans:
    fig.add_trace(
        go.Scatter(
            x=[start.isoformat(), end.isoformat()],
            y=[row, row],
            mode="lines",
            line={"width": 16, "color": span_colors[kind]},
            name=kind,
            legendgroup=kind,
            showlegend=kind not in seen,
        )
    )
    seen.add(kind)

fig.add_vrect(
    x0=HOLDOUT_START.isoformat(),
    x1=HOLDOUT_END_EXCLUSIVE.isoformat(),
    fillcolor=COLORS["recede"],
    opacity=0.12,
    line_width=0,
    layer="below",
)
fig.add_vline(x=HOLDOUT_START.isoformat(), line_dash="dash", line_color=COLORS["negative"])
fig.update_layout(
    title=(
        "Every fold trains left of the validation span it is scored on"
        "<br><sup>The dashed rule is where the holdout opens and the shaded region is the"
        "<br>holdout itself. The last fold is the one written so that a model scored on the"
        "<br>holdout has features there; its training bars all predate the rule. Spans overlap"
        "<br>between folds because the fold tag selects rows rather than changing values.</sup>"
    ),
    xaxis_title="Session",
    yaxis_title="",
    height=460,
    margin={"l": 90, "t": 140},
)
show_plotly_with_alt(
    fig,
    "Horizontal timeline with one row per fold on a session axis, lowest-numbered fold at the "
    "bottom. Every row below the top is a validation fold: a long dark navy training bar "
    "followed by the shorter amber validation bar it is scored on. The validation windows sit "
    "at different points along the axis and the training bars overlap between rows, because a "
    "fold tag selects rows rather than changing values. A dashed red rule marks where the "
    "holdout opens and a shaded band to its right is the holdout itself. The top row is the "
    "extra fold written for the holdout and has no validation bar: its training bar runs from "
    "the left edge up to the rule and its light grey holdout bar sits inside the band, later "
    "than every validation window. No training bar of any row crosses the rule.",
)

# %% [markdown]
# ## C. One section per model: what it infers, and why it cannot see ahead
#
# Each of the three procedures reads a trailing window and writes a value for the bar at the end
# of it. Those windows are counted in bars within a symbol, and a symbol's bars are the
# sessions laid end to end, so a window that is longer than the distance back to the session
# open reaches across the overnight gap. The return series itself does not - `r1m` is null at
# each open and enters the windows below as a zero - but the aggregation is not restarted, so a
# bar early in the session is described partly by yesterday afternoon.
#
# The share of rows this affects is the window length over the session length, which is worth
# measuring rather than asserting: `bar_of_day` is a row's position in its own session, so a row
# with `bar_of_day` below the window is one whose window crosses the gap. A production system
# would bound each window by the session. The approximation is kept here because it makes the cost of the shortcut visible and because it is the cost that
# grows with the window - which is the reason the fit window is the shortest one that identifies
# the regression rather than the longest one available.

# %%
_bod = validation_rows(df.select("timestamp", "bar_of_day"))
display(
    pl.DataFrame(
        {
            "window": ["signature path", "spectrum", "HAR components", "HAR fit window"],
            "bars": [SIG_WINDOW, FFT_WINDOW, HAR_COMPONENTS[-1], HAR_FIT_WINDOW],
            "share of rows reaching across a session gap": [
                f"{_bod.select((pl.col('bar_of_day') < w).mean()).item():.1%}"
                for w in (SIG_WINDOW, FFT_WINDOW, HAR_COMPONENTS[-1], HAR_FIT_WINDOW)
            ],
        }
    )
)
del _bod

# %% [markdown]
# ### C.1 HAR: a variance forecast built from three lengths of history
#
# Realized volatility is persistent, and it is persistent at more than one time scale at once:
# what happened in the last five minutes, the last quarter of an hour and the last hour all say
# something, and they do not say the same thing. The heterogeneous autoregressive model (Corsi,
# 2009) is the simplest way to use all three - a linear regression of the next period's realized
# variance on the realized variance measured over each of those three lengths:
#
# $$RV_{t+1}^{(5)} = c + \beta_5 \, RV_t^{(5)} + \beta_{15} \, RV_t^{(15)} + \beta_{60} \, RV_t^{(60)} + \varepsilon_{t+1}$$
#
# Corsi's original components are a day, a week and a month, on the argument that different
# participants look at different lengths of history. On a minute grid the same argument gives
# minutes, quarter hours and hours, which is what the three components here are.
#
# Two features come out. The **forecast** is the fitted right-hand side, the model's statement
# about the variance of the next few minutes. The **residual** is what the last such statement
# got wrong - realized variance minus the forecast made for it - and it is the more interesting
# of the two, because a large positive residual is variance that arrived without the recent past
# implying it: a news arrival, a liquidity event, something the persistence did not contain.
#
# Realized variance at horizon $w$ is the average squared one-minute return over the $w$ bars
# **before** $t$, so the value at $t$ never includes the bar at $t$.


# %%
def build_har_features_intraday(
    r1m: np.ndarray, components: tuple[int, int, int] = HAR_COMPONENTS
) -> dict[str, np.ndarray]:
    """Realized variance at each of the three HAR horizons.

    Every window ends at ``t`` exclusive, so the value at ``t`` is a function of bars strictly
    before ``t``.
    """
    n = len(r1m)
    r2 = r1m**2
    window_5, window_15, window_60 = components

    rv_5 = np.full(n, np.nan)
    rv_15 = np.full(n, np.nan)
    rv_60 = np.full(n, np.nan)

    for t in range(window_60, n):
        rv_5[t] = np.mean(r2[t - window_5 : t])
        rv_15[t] = np.mean(r2[t - window_15 : t])
        rv_60[t] = np.mean(r2[t - window_60 : t])

    return {"rv_5m": rv_5, "rv_15m": rv_15, "rv_60m": rv_60}


# %% [markdown]
# The regression is refitted at every bar on the immediately preceding stretch of history. That
# is the discipline Section A described, applied at the finest cadence available: the
# coefficients used to describe bar $t$ come from a regression whose last observation is bar
# $t-1$, so there is no window in which a parameter and the bar it describes share information.
# A refit that reads fewer than `HAR_MIN_TRAIN_OBS` usable observations is skipped rather than
# fitted on whatever survived, and the bar keeps a null.
#
# The schedule that walks the fit forward is not written here. `walk_forward_feature` in
# `case_studies/utils/temporal.py` owns it for every case study that fits anything, and it is
# what makes the two channels a fitted feature carries - which bars a value is computed from,
# and which bars its parameters came from - end at or before the bar being described. What this
# notebook supplies is the two halves that are specific to a HAR: how one refit window is
# turned into coefficients, and how those coefficients produce a bar's forecast.
#
# `apply_scope="block"` is the second half of that. A GARCH or a Kalman filter has to be run
# from the start of the series to reach the bar it is describing, so the harness hands it every
# row up to the block. A HAR forecast is a dot product of the bar's own three regressors with
# the coefficients, and nothing earlier enters it - so the harness hands it the block alone.
# At one refit per bar over 174,000 bars the difference is not stylistic: the prefix form would
# ask for 1.5e10 rows per symbol to keep 174,000 of them.


# %%
def fit_har_window(block: np.ndarray) -> np.ndarray:
    """Least squares of the next bar's short-horizon variance on the three components.

    ``block`` is one refit window, ``(fit_window, 4)``: the three regressors at each bar,
    followed by the value being regressed on them, which is the short-horizon variance of the
    bar after. Carrying the target as a column of the same array is what keeps the fit inside
    the window the schedule handed over - a target read from outside it would be an
    observation the schedule never sanctioned.
    """
    regressors = np.column_stack([np.ones(len(block)), block[:, :3]])
    target = block[:, 3]
    usable = np.isfinite(target) & np.all(np.isfinite(regressors), axis=1)
    if usable.sum() < HAR_MIN_TRAIN_OBS:
        raise ValueError(
            f"{usable.sum()} usable observations in the refit window, "
            f"below the declared {HAR_MIN_TRAIN_OBS}"
        )
    return np.linalg.lstsq(regressors[usable], target[usable], rcond=None)[0]


def apply_har(beta: np.ndarray, block: np.ndarray) -> np.ndarray:
    """The forecast for each bar of the block, and the coefficients that produced it.

    The coefficients travel with the forecast because Section D reads them: they are what says
    which of the three horizons the fit is leaning on, and they exist per bar for the same
    reason the forecast does.
    """
    regressors = np.column_stack([np.ones(len(block)), block[:, :3]])
    forecast = np.where(np.all(np.isfinite(regressors), axis=1), regressors @ beta, np.nan)
    return np.column_stack([forecast, np.tile(beta[1:], (len(block), 1))])


# %% [markdown]
# One symbol at a time: build the three regressors, roll the fit across them, and keep the
# coefficients as well as the features. The coefficients are thinned to one per symbol-session
# before they leave the function, because Section D asks how the fit moves over months and
# twenty million rows of it would answer that no better than fifty thousand.


# %%
def compute_har_per_symbol(
    symbol_df: pl.DataFrame,
) -> tuple[pl.DataFrame, pl.DataFrame]:
    """HAR features and rolling coefficients for one symbol's sessions."""
    r1m = symbol_df["r1m"].to_numpy().copy()
    r1m = np.nan_to_num(r1m, nan=0.0)

    har_regs = build_har_features_intraday(r1m)
    rv_5 = har_regs["rv_5m"]
    # The fourth column is the target: the short-horizon variance of the bar after this one.
    # It is a column of the same array rather than a second argument because the walk hands
    # `fit_har_window` a slice and nothing else, which is what bounds the fit.
    series = np.column_stack(
        [rv_5, har_regs["rv_15m"], har_regs["rv_60m"], np.append(rv_5[1:], np.nan)]
    )

    emitted = walk_forward_feature(
        series,
        # One symbol's bars, already sorted: the caller below partitions by symbol and sorts by
        # timestamp, and this is what makes the harness refuse the call rather than trust it.
        timestamps=symbol_df["timestamp"],
        burnin=HAR_BURNIN,
        refit_every=HAR_REFIT_EVERY,
        window=HAR_FIT_WINDOW,
        fit=fit_har_window,
        apply=apply_har,
        apply_scope="block",
        n_features=4,
        # A window too thin to identify the regression is a statement about that stretch of the
        # symbol's history, not a failure of the notebook. The bar keeps a null and the walk
        # carries on to the next refit.
        on_fit_error="skip",
    )
    har_forecast, har_betas = emitted[:, 0], emitted[:, 1:]

    # The residual is what the PREVIOUS bar's forecast got wrong, so it is read off the emitted
    # column rather than computed inside the walk: `har_residual[t] = rv_5[t] - forecast[t-1]`.
    # It exists only where both forecasts do - a bar whose own refit was skipped has no
    # coefficients to be judged against, which is the rule the loop this replaces applied.
    previous_forecast = np.concatenate([[np.nan], har_forecast[:-1]])
    har_residual = np.where(
        np.isfinite(har_forecast) & np.isfinite(previous_forecast),
        rv_5 - previous_forecast,
        np.nan,
    )

    features = pl.DataFrame(
        {
            "timestamp": symbol_df["timestamp"],
            "symbol": symbol_df["symbol"],
            "har_rv5_pred": har_forecast,
            "har_residual": har_residual,
        }
    )

    betas = (
        pl.DataFrame(
            {
                "timestamp": symbol_df["timestamp"],
                "symbol": symbol_df["symbol"],
                "session_date": symbol_df["session_date"],
                "bar_of_day": symbol_df["bar_of_day"],
                "beta_5": har_betas[:, 0],
                "beta_15": har_betas[:, 1],
                "beta_60": har_betas[:, 2],
            }
        )
        .filter(pl.col("bar_of_day") == pl.col("bar_of_day").max().over("session_date"))
        .drop("bar_of_day")
        .with_columns(pl.col("^beta_.*$").fill_nan(None))
        .drop_nulls(["beta_5", "beta_15", "beta_60"])
    )

    return features, betas


# %%
symbols = df["symbol"].unique().sort().to_list()
har_results = []
beta_results = []

for i, sym in enumerate(symbols):
    sym_df = df.filter(pl.col("symbol") == sym).sort("timestamp")
    result, betas = compute_har_per_symbol(sym_df)
    har_results.append(result)
    beta_results.append(betas)
    if (i + 1) % 20 == 0 or (i + 1) == len(symbols):
        print(f"  HAR: {i + 1}/{len(symbols)} symbols processed")

har_df = pl.concat(har_results)
har_beta_df = pl.concat(beta_results)
del har_results, beta_results

for c in ["har_rv5_pred", "har_residual"]:
    har_df = har_df.with_columns(pl.col(c).fill_nan(None))

print(
    f"HAR: {har_df['har_rv5_pred'].drop_nulls().len():,} forecasts on {har_df.height:,} bars, "
    f"from {har_beta_df.height:,} retained fits."
)

# %% [markdown]
# There is no single representative set of coefficients to quote: the model is refitted every
# bar, so the object is the distribution of those fits. The medians below say which of the three
# horizons the fit leans on, and their sum says how much of a variance shock the model expects
# to still be there next period. They are taken over validation rows, like every other readout
# here.

# %%
display(
    validation_rows(har_beta_df).select(
        pl.col("beta_5", "beta_15", "beta_60").median().round(4).name.suffix("_median"),
        (pl.col("beta_5") + pl.col("beta_15") + pl.col("beta_60"))
        .median()
        .round(4)
        .alias("persistence_median"),
        pl.len().alias("fits"),
    )
)

# %% [markdown]
# **The forecast is an unconstrained linear extrapolation, and it shows.** The HAR is a
# regression on a variance with nothing in it that keeps a prediction non-negative. When a
# symbol's realized variance jumps well outside the range the trailing window was fitted on -
# a single-name event, an earnings gap, a halt - the fit extrapolates and the forecast can land
# below zero. The distribution below is reported as quantiles rather than as a mean and a
# standard deviation, because those two are set by a handful of such rows and describe nothing
# a reader can use. Modelling the logarithm of the variance is the standard remedy; leaving the
# forecast unconstrained here is what makes the failure mode visible instead of quietly clipped.

# %%
_fc = validation_rows(har_df.select("timestamp", "har_rv5_pred"))["har_rv5_pred"].drop_nulls()
display(
    pl.DataFrame(
        {
            "statistic": [
                "rows",
                "1st percentile",
                "median",
                "99th percentile",
                "minimum",
                "share below zero",
            ],
            "value": [
                f"{len(_fc):,}",
                f"{_fc.quantile(0.01):.3e}",
                f"{_fc.median():.3e}",
                f"{_fc.quantile(0.99):.3e}",
                f"{_fc.min():.3e}",
                f"{(_fc < 0).mean():.2%}",
            ],
        }
    )
)
del _fc

# %% [markdown]
# **Figure F2** shows what the fitted model actually inferred. The forecast is drawn against the
# realized variance it was forecasting, on validation rows only, with the boundaries between
# consecutive validation windows marked. Both series are the cross-sectional median across symbols
# within a session, because one symbol's minute-level realized variance is far too noisy to read
# over a year.
#
# The realized series is not fetched from anywhere: it is reconstructed from what this notebook
# emits, since `har_residual[t] = rv_5[t] - har_forecast[t-1]` rearranges to
# `rv_5[t] = har_forecast[t-1] + har_residual[t]`. Reading the two columns back that way is also
# a check that they mean what the docstring says they mean.

# %%
har_view = (
    validation_rows(har_df.select("timestamp", "symbol", "har_rv5_pred", "har_residual"))
    .sort(["symbol", "timestamp"])
    .with_columns(realized=pl.col("har_rv5_pred").shift(1).over("symbol") + pl.col("har_residual"))
    .with_columns(session=pl.col("timestamp").dt.date())
    .group_by("session")
    .agg(
        pl.col("har_rv5_pred").median().alias("forecast"),
        pl.col("realized").median().alias("realized"),
    )
    .sort("session")
)
print(f"Validation sessions drawn: {len(har_view):,}")

# %%
fig = go.Figure()
fig.add_trace(
    go.Scatter(
        x=har_view["session"],
        y=har_view["realized"],
        mode="lines",
        name="Realized 5-bar variance",
        line={"color": COLORS["blue"], "width": 1.5},
    )
)
fig.add_trace(
    go.Scatter(
        x=har_view["session"],
        y=har_view["forecast"],
        mode="lines",
        name="HAR forecast",
        line={"color": COLORS["amber"], "width": 2},
    )
)
for s in sorted(splits, key=lambda s: pd.Timestamp(s["val_start"]))[1:]:
    fig.add_vline(
        x=pd.Timestamp(s["val_start"]).isoformat(),
        line_dash="dot",
        line_color=COLORS["neutral"],
    )
fig.update_layout(
    title=(
        "The HAR forecast tracks realized variance closely and overshoots its peaks"
        "<br><sup>Cross-sectional median across symbols per session, validation rows only."
        "<br>Each dotted rule is a boundary between consecutive validation windows. Both series"
        "<br>are means of squared one-minute log returns.</sup>"
    ),
    xaxis_title="Session",
    yaxis_title="Mean squared 1-minute log return",
    height=440,
    margin={"t": 120},
)
show_plotly_with_alt(
    fig,
    "Two lines over validation sessions on a shared axis of mean squared one-minute log return. "
    "A thin dark navy line is realized 5-bar variance, spiky, with occasional tall isolated "
    "peaks. A thicker amber line is the HAR forecast, tracking the same path closely and rising "
    "above the navy line at every peak. A single dotted vertical rule near the middle is the "
    "boundary between the two consecutive validation windows.",
)

# %% [markdown]
# ### C.2 A rolling spectrum of volume and of variance
#
# Intraday activity is repetitive. Volume is heavy at the open, thins through the middle of the
# day and picks up into the close, and that shape repeats every session; volatility clusters at
# its own frequencies. A Fourier transform of a trailing window is a way of asking how much of
# the recent activity sits at which repetition rate, and it produces conditioning features - not
# a prediction of direction, but a description of what kind of hour this is.
#
# Four numbers come out of each window. **Spectral energy** is how much variation there is in
# total once the level is removed. The **dominant period** is the repetition length carrying the
# most of it, in bars. **Spectral entropy** is how evenly the variation is spread across
# frequencies: low when one rhythm dominates, high when the window is closer to noise. The
# **low-frequency ratio** is the share sitting at periods longer than `FFT_LOW_FREQ_PERIOD`
# bars, which separates a slow drift in activity from minute-to-minute churn.
#
# Each window ends at $t$ exclusive, so nothing at or after the bar being described enters its
# own spectrum.


# %%
def rolling_fft_features(
    signal: np.ndarray,
    window: int = FFT_WINDOW,
    low_frequency_period: int = FFT_LOW_FREQ_PERIOD,
) -> dict[str, np.ndarray]:
    """Four descriptions of the power spectrum of each trailing window of *signal*."""
    n = len(signal)
    spectral_energy = np.full(n, np.nan)
    dominant_period = np.full(n, np.nan)
    spectral_entropy = np.full(n, np.nan)
    low_freq_ratio = np.full(n, np.nan)

    for t in range(window, n):
        segment = signal[t - window : t]

        if np.all(np.isnan(segment)) or np.nanstd(segment) < 1e-12:
            continue

        seg_clean = np.nan_to_num(segment, nan=0.0)
        seg_clean = seg_clean - seg_clean.mean()

        fft_vals = np.fft.rfft(seg_clean)
        power = np.abs(fft_vals) ** 2
        freqs = np.fft.rfftfreq(window)

        total_power = np.sum(power[1:])
        if total_power <= 0:
            continue

        spectral_energy[t] = total_power

        dom_idx = np.argmax(power[1:]) + 1
        if freqs[dom_idx] > 0:
            dominant_period[t] = 1.0 / freqs[dom_idx]

        p_norm = power[1:] / total_power
        p_norm = p_norm[p_norm > 0]
        spectral_entropy[t] = -np.sum(p_norm * np.log(p_norm))

        low_mask = freqs[1:] < (1.0 / low_frequency_period)
        if low_mask.any():
            low_freq_ratio[t] = np.sum(power[1:][low_mask]) / total_power

    return {
        "spectral_energy": spectral_energy,
        "dominant_period": dominant_period,
        "spectral_entropy": spectral_entropy,
        "low_freq_ratio": low_freq_ratio,
    }


# %% [markdown]
# The transform is run twice per symbol, on two different signals. Volume answers how structured
# the recent activity pattern was; squared returns answer the same question about volatility.
# Volume is passed through a logarithm first, because raw share counts span several orders of
# magnitude within a session and a spectrum of them is dominated by the largest few bars.


# %%
def compute_fft_per_symbol(
    symbol_df: pl.DataFrame,
    window: int = FFT_WINDOW,
) -> pl.DataFrame:
    """Spectral descriptions of trailing volume and squared-return windows for one symbol."""
    vol_raw = symbol_df["volume"].to_numpy().astype(float)
    vol_signal = np.log1p(np.clip(vol_raw, 0, None))

    r1m = symbol_df["r1m"].to_numpy().copy()
    r2_signal = np.nan_to_num(r1m, nan=0.0) ** 2

    vol_fft = rolling_fft_features(vol_signal, window=window)
    r2_fft = rolling_fft_features(r2_signal, window=window)

    return pl.DataFrame(
        {
            "timestamp": symbol_df["timestamp"],
            "symbol": symbol_df["symbol"],
            "vol_spectral_energy": vol_fft["spectral_energy"],
            "vol_dominant_period": vol_fft["dominant_period"],
            "vol_spectral_entropy": vol_fft["spectral_entropy"],
            "vol_low_freq_ratio": vol_fft["low_freq_ratio"],
            "rv_spectral_energy": r2_fft["spectral_energy"],
            "rv_dominant_period": r2_fft["dominant_period"],
            "rv_spectral_entropy": r2_fft["spectral_entropy"],
            "rv_low_freq_ratio": r2_fft["low_freq_ratio"],
        }
    )


# %%
fft_results = []

for i, sym in enumerate(symbols):
    sym_df = df.filter(pl.col("symbol") == sym).sort("timestamp")
    fft_results.append(compute_fft_per_symbol(sym_df, window=FFT_WINDOW))
    if (i + 1) % 20 == 0 or (i + 1) == len(symbols):
        print(f"  FFT: {i + 1}/{len(symbols)} symbols processed")

fft_df = pl.concat(fft_results)
del fft_results

fft_feature_cols = [c for c in fft_df.columns if c not in ["timestamp", "symbol"]]
for c in fft_feature_cols:
    fft_df = fft_df.with_columns(pl.col(c).fill_nan(None))

print(f"FFT: {fft_df['vol_spectral_energy'].drop_nulls().len():,} of {fft_df.height:,} bars.")

# %% [markdown]
# ### C.3 Path signatures: which moved first, price or flow
#
# The microstructure question this case study is built around is whether order flow leads price
# or price leads order flow. The distinction matters: flow arriving before a move is what an
# informed trade looks like, and a move arriving before flow is what liquidity chasing a price
# looks like, and they imply opposite things about whether the next minute continues.
#
# A **path signature** is a way of summarising a multi-dimensional path so that the order in
# which its dimensions moved is still readable off the summary. Take the three series - price,
# signed volume share, trade count - over a trailing window and treat them as one path through
# three dimensions. The signature is a sequence of iterated integrals of that path. Truncated at
# depth two it is $d + d^2$ numbers for a $d$-dimensional path: $d$ net displacements, one per
# dimension, and $d^2$ cross terms.
#
# The cross terms are the point. The term $S^{i,j}$ accumulates movement in dimension $j$
# weighted by how far dimension $i$ has already travelled, so it is large when $i$ moved first.
# `sig2_svs_ret` large means flow moved before price; `sig2_ret_svs` large means price moved
# first. Their difference is the asymmetry the question is about, and neither a correlation nor
# a lagged regression puts it in one number this way.
#
# Depth two has a closed form, so no library is needed.


# %%
def compute_depth2_signature(path: np.ndarray) -> np.ndarray:
    """The depth-2 truncated signature of a path of shape ``(T, d)``.

    Returns ``d`` net displacements followed by the ``d * d`` iterated integrals
    ``S^{i,j} = int int_{s<t} dX^i_s dX^j_t``, flattened row-major.

    The path is piecewise linear between samples, so a pair of increments contributes to
    ``S^{i,j}`` in two ways: whole earlier segments, ``dX^i_s dX^j_t`` for ``s < t``, and the
    half of each segment that lies below its own diagonal, ``0.5 dX^i_t dX^j_t``. Dropping the
    second is the difference between the signature and a strictly-lagged double sum, and it is
    visible in the diagonal: the identity below fails without it, and ``S^{i,i}`` goes negative
    whenever the increments partly cancel.
    """
    T, d = path.shape
    increments = np.diff(path, axis=0)

    sig1 = path[-1] - path[0]

    sig2 = np.zeros((d, d))
    cumsum = np.zeros(d)
    for t in range(len(increments)):
        sig2 += np.outer(cumsum, increments[t]) + 0.5 * np.outer(increments[t], increments[t])
        cumsum += increments[t]

    return np.concatenate([sig1, sig2.ravel()])


# %% [markdown]
# Two identities hold for the depth-2 signature of any path, whatever the path is, so they are
# what says the implementation computes a signature rather than something that resembles one.
# The diagonal is fixed by the net displacement alone, $S^{i,i} = \frac{1}{2}(\Delta X^i)^2$,
# which also makes it non-negative; and the shuffle relation
# $S^{i,j} + S^{j,i} = \Delta X^i \Delta X^j$ says the symmetric part carries no information
# beyond depth one, which is why the *antisymmetric* part is the feature worth reading.

# %%
_rng = np.random.default_rng(0)
for _trial in range(20):
    _p = np.cumsum(_rng.standard_normal((30, 3)), axis=0)
    _sig = compute_depth2_signature(_p)
    _s1, _s2 = _sig[:3], _sig[3:].reshape(3, 3)
    assert np.allclose(np.diag(_s2), 0.5 * _s1**2), (
        "depth-2 diagonal is not half the squared net move"
    )
    assert np.allclose(_s2 + _s2.T, np.outer(_s1, _s1)), "depth-2 shuffle identity fails"
print("Depth-2 signature identities hold on 20 random 3-dimensional paths.")


# %%
def _window_normalize(x: np.ndarray) -> np.ndarray:
    """Centre and scale a window by its own mean and standard deviation."""
    s = np.std(x)
    if s < 1e-12:
        return x - np.mean(x)
    return (x - np.mean(x)) / s


# %% [markdown]
# The three dimensions arrive on wildly different scales - a log return near $10^{-4}$, a share
# between minus one and one, a trade count in the hundreds - and a signature of the raw path
# would be a description of those scales rather than of the path's shape. Each window is
# therefore standardised by its **own** mean and standard deviation before the path is built.
# That is what keeps the signature free of any quantity computed over the whole sample: the
# scale each window is put on comes from the window, never from a constant estimated across the
# symbol's history, which would be exactly the leak Section A is about.


# %%
def compute_signatures_per_symbol(
    symbol_df: pl.DataFrame,
    window: int = SIG_WINDOW,
) -> pl.DataFrame:
    """Rolling depth-2 signatures of the (return, signed volume share, trades) path."""
    r1m = symbol_df["r1m"].to_numpy().copy()
    svs = symbol_df["signed_vol_share"].to_numpy().copy()
    trades = symbol_df["total_trades"].to_numpy().astype(float).copy()

    r1m = np.nan_to_num(r1m, nan=0.0)
    svs = np.nan_to_num(svs, nan=0.0)
    trades = np.nan_to_num(trades, nan=0.0)

    n = len(r1m)
    d = 3
    n_features = d + d * d

    sig_features = np.full((n, n_features), np.nan)

    for t in range(window, n):
        seg_r = np.cumsum(_window_normalize(r1m[t - window : t]))
        seg_svs = np.cumsum(_window_normalize(svs[t - window : t]))
        seg_trades = np.cumsum(_window_normalize(trades[t - window : t]))

        path = np.column_stack([seg_r, seg_svs, seg_trades])

        if np.all(np.abs(np.diff(path, axis=0)) < 1e-12):
            continue

        sig_features[t] = compute_depth2_signature(path)

    dims = ["ret", "svs", "trd"]
    col_names = [f"sig1_{name}" for name in dims]
    col_names += [f"sig2_{name_i}_{name_j}" for name_i in dims for name_j in dims]

    result = {
        "timestamp": symbol_df["timestamp"],
        "symbol": symbol_df["symbol"],
    }
    for k, col_name in enumerate(col_names):
        result[col_name] = sig_features[:, k]

    return pl.DataFrame(result)


# %%
sig_results = []

for i, sym in enumerate(symbols):
    sym_df = df.filter(pl.col("symbol") == sym).sort("timestamp")
    sig_results.append(compute_signatures_per_symbol(sym_df, window=SIG_WINDOW))
    if (i + 1) % 20 == 0 or (i + 1) == len(symbols):
        print(f"  Signatures: {i + 1}/{len(symbols)} symbols processed")

sig_df = pl.concat(sig_results)
del sig_results

sig_feature_cols = [c for c in sig_df.columns if c not in ["timestamp", "symbol"]]
for c in sig_feature_cols:
    sig_df = sig_df.with_columns(pl.col(c).fill_nan(None))

print(f"Signatures: {sig_df['sig1_ret'].drop_nulls().len():,} of {sig_df.height:,} bars.")

# %% [markdown]
# ### C.4 Withholding the future changes nothing
#
# Every claim made so far about these three procedures is a claim that a value at $t$ is a
# function of bars before $t$. A notebook cannot establish that by agreeing with itself, so the
# check below computes the features a second time from a different input: the same code on a
# panel that stops at the holdout boundary. If any window, any fit or any normalisation reached
# forward, truncating the panel would move the values on the rows the two runs share.
#
# It runs on three symbols rather than the whole universe because the property is a property of
# the code, not of the sample, and three symbols is enough for a difference to appear. Exact
# equality is the bar - not a tolerance - because these are the same arithmetic on the same
# bars.

# %%
_check_symbols = symbols[:3]
_full = (
    har_df.join(fft_df, on=["timestamp", "symbol"], how="inner")
    .join(sig_df, on=["timestamp", "symbol"], how="inner")
    .filter(pl.col("symbol").is_in(_check_symbols) & (pl.col("timestamp") < HOLDOUT_START))
    .sort(["symbol", "timestamp"])
)
_truncated_parts = []
for sym in _check_symbols:
    _sym = df.filter((pl.col("symbol") == sym) & (pl.col("timestamp") < HOLDOUT_START)).sort(
        "timestamp"
    )
    _h, _ = compute_har_per_symbol(_sym)
    _truncated_parts.append(
        _h.join(compute_fft_per_symbol(_sym, window=FFT_WINDOW), on=["timestamp", "symbol"]).join(
            compute_signatures_per_symbol(_sym, window=SIG_WINDOW), on=["timestamp", "symbol"]
        )
    )
_truncated = (
    pl.concat(_truncated_parts)
    .with_columns(
        [pl.col(c).fill_nan(None) for c in _full.columns if c not in ("timestamp", "symbol")]
    )
    .sort(["symbol", "timestamp"])
    .select(_full.columns)
)
assert _truncated.equals(_full), (
    "a feature moved when the panel was trunc