ساخت ویژگی‌های علّی نوسان، طیفی و امضای مسیر

# NASDAQ-100 microstructure: features a model has to be fitted to produce


# 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).

```python
"""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")
```

```python
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 = ""
```

## 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.

```python
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)
```

```python
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.")
```

## 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.

```python
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
```

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.

```python
_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")
```

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.

```python
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}",
            ],
        }
    )
)
```

### 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.

```python
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}")
```

## 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.

## 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.

```python
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"])
)
```

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.

```python
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.")
```

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.

```python
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."
)
```

```python
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
        ]
    )
```

**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.

```python
# 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.",
)
```

## 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.

```python
_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
```

### 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$.

```python
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}
```

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.

```python
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))])
```

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.

```python
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
```

```python
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."
)
```

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.

```python
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"),
    )
)
```

**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.

```python
_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
```

**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.

```python
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):,}")
```

```python
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.",
)
```

### 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.

```python
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,
    }
```

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.

```python
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"],
        }
    )
```

```python
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.")
```

### 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.

```python
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()])
```

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.

```python
_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.")
```

```python
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
```

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.

```python
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)
```

```python
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.")
```

### 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.

```python
_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 truncated at the holdout boundary: something reads ahead"
)
print(
    f"{_full.height:,} rows on {len(_check_symbols)} symbols recomputed from a panel ending "
    f"{HOLDOUT_START.date()}; every value identical."
)
# The minute panel's own decision times, kept because the coverage check in Section E is a
# comparison between what the panel offered and what the artifact carries. Taken from `df`
# rather than from the feature frame, which is the point: the two differ by exactly the
# warm-up rows, and a check that read the feature frame on both sides would be comparing a
# frame against itself.
PANEL_TS = set(df["timestamp"].unique().to_list())
del _full, _truncated, _truncated_parts, df
```

## D. Fit stability across folds

Only one of the three procedures has parameters, and it is refitted every bar, so the question
"do the parameters move as the window rolls" is a question about a distribution rather than
about three numbers per fold. **Figure F3** shows that distribution, one box per component pe