पैनल फ़ीचर: कॉइंटीग्रेशन, गतिशील हेज और क्रॉस-सेक्शनल रैंक

# Panel Features


# Panel Features

**Chapter 9 | Section 9.6**

**Docker image**: `ml4t`

Every feature so far has been computed from one series. This notebook computes features that
need two or more, and they come in two shapes that have almost nothing to do with each other
beyond both requiring a panel.

A **pairwise** feature reads two price series and asks whether a combination of them is
stationary when neither is. That is the cointegration question, and its output is a spread, a
hedge ratio, and a speed of reversion.

A **cross-sectional** feature reads one quantity across many assets on the same date and
replaces its level with its position among them. An annualized volatility of a quarter means
one thing for a Treasury fund and another for an energy fund, and a rank says which of the
two a reader is looking at without needing to know.

**Learning objectives**

- Test a pair for cointegration two ways and read what it means when the two disagree.
- Estimate a hedge ratio that changes over time with a Kalman filter, and see where it
  differs from the single number a full-sample regression gives.
- Estimate a mean-reversion half-life on a first block and use it to size the window a
  trading signal is computed over, so the window length is not chosen with future data.
- Convert a temporal feature into a cross-sectional rank and a benchmark-relative value, and
  say which of the two a downstream model wants.
- Aggregate a per-asset regime measure across a universe without ranking each asset against
  its own future.

**Book reference**

Chapter 9, Section 9.6 (Cross-sectional and panel features).

**Prerequisites**

`04_kalman_filter` for the recursion the hedge ratio uses. `11_hmm_regimes` for the regime
measures the last section aggregates.

## Setup

```python
"""Panel features - cointegrated pairs, cross-sectional ranks and universe aggregates."""

import warnings

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import polars as pl
from IPython.display import display
from ml4t.engineer.features.cross_asset import (
    beta_to_market,
    co_integration_score,
    correlation_regime_indicator,
    rolling_correlation,
)
from sklearn.linear_model import LinearRegression
from statsmodels.tsa.stattools import adfuller, coint
from statsmodels.tsa.vector_ar.vecm import coint_johansen

from data import load_etfs
from utils.style import COLORS, FIGSIZE, show_with_alt

with warnings.catch_warnings():
    # filterpy carries an invalid escape sequence in a docstring, which Python reports the first
    # time it byte-compiles the module and never again. Nothing about the filter is affected.
    warnings.simplefilter("ignore", SyntaxWarning)
    from filterpy.kalman import KalmanFilter
```

```python
START_DATE = "2015-01-01"
END_DATE = "2024-12-31"
MAX_SYMBOLS = 0  # zero loads every symbol in the universe
SIGNIFICANCE = 0.05
TRAIN_FRACTION = 0.5
ENTRY_THRESHOLD = 2.0
MINIMUM_LOOKBACK = 20
MINIMUM_SESSIONS = 252
CORRELATION_WINDOW = 60
COINTEGRATION_WINDOW = 120
FEATURE_WINDOW = 60
REGIME_WINDOW = 252
SESSIONS_PER_YEAR = 252
```

## The universe

The ETF universe, because a pair needs two series whose relationship has a reason behind it
and exchange-traded funds come with the reason written on them. The pair used through the
first half is an energy sector fund against a crude oil fund, and the fallback pair for a
reduced fixture is two broad equity funds.

```python
universe = load_etfs().filter(
    pl.col("timestamp").is_between(pl.lit(START_DATE).str.to_date(), pl.lit(END_DATE).str.to_date())
)
if MAX_SYMBOLS > 0:
    kept = universe["symbol"].unique().sort().to_list()[:MAX_SYMBOLS]
    universe = universe.filter(pl.col("symbol").is_in(kept))

available = set(universe["symbol"].unique().to_list())
DEPENDENT, INDEPENDENT = ("XLE", "USO") if {"XLE", "USO"} <= available else ("SPY", "QQQ")


def pair_frame(dependent: str, independent: str) -> pd.DataFrame:
    """Closing prices of two symbols on the sessions both traded, indexed by date."""
    columns = []
    for symbol, name in ((dependent, "dependent"), (independent, "independent")):
        columns.append(
            universe.filter(pl.col("symbol") == symbol)
            .select(["timestamp", "close"])
            .rename({"close": name})
            .sort("timestamp")
        )
    joined = columns[0].join(columns[1], on="timestamp", how="inner").drop_nulls()
    frame = joined.to_pandas().set_index("timestamp")
    frame.index = pd.DatetimeIndex(frame.index)
    return frame


pair = pair_frame(DEPENDENT, INDEPENDENT)
train_end = int(len(pair) * TRAIN_FRACTION)

print(f"Universe: {len(available)} symbols, {universe.height:,} rows")
print(f"Pair: {DEPENDENT} against {INDEPENDENT}, {len(pair):,} shared sessions")
print(f"Sessions used to estimate before any signal is taken: {train_end:,}")
display(pair.tail(3))
```

## Two tests for one question

Two price series are **cointegrated** when some fixed combination of them is stationary even
though neither of them is. That is a stronger statement than correlation and a different one:
two series can move together every day and drift apart without limit over years, and it is
the drift that decides whether a spread between them comes back.

The two standard tests reach the question differently. **Engle-Granger** regresses one series
on the other and tests the residual for a unit root, so it takes one of the two as the
dependent variable and gives a different answer if the roles are swapped. **Johansen** treats
the pair as a system and tests how many stationary combinations exist, which is symmetric in
the two series and extends past two of them.

Both are run on the whole sample here, and that is the correct thing for what they are being
used for. Whether a pair is worth trading is a question asked before any trading, out of the
history available at that point; the answer is not a feature and does not enter a design
matrix. What must not use the whole sample is anything the signal below reads.

```python
statistic, p_value, _ = coint(pair["dependent"], pair["independent"], trend="c")[:3]

johansen = coint_johansen(pair[["dependent", "independent"]].to_numpy(), det_order=0, k_ar_diff=1)
trace_statistics = johansen.lr1
critical_values = johansen.cvt[:, 1]
johansen_rejects = bool(trace_statistics[0] > critical_values[0])
johansen_vector = johansen.evec[:, 0]

display(
    pd.DataFrame(
        [
            {
                "test": "Engle-Granger",
                "statistic": statistic,
                "threshold": np.nan,
                "p value": p_value,
                "rejects no cointegration": bool(p_value < SIGNIFICANCE),
            },
            {
                "test": "Johansen, no stationary combination",
                "statistic": trace_statistics[0],
                "threshold": critical_values[0],
                "p value": np.nan,
                "rejects no cointegration": johansen_rejects,
            },
            {
                "test": "Johansen, at most one",
                "statistic": trace_statistics[1],
                "threshold": critical_values[1],
                "p value": np.nan,
                "rejects no cointegration": bool(trace_statistics[1] > critical_values[1]),
            },
        ]
    ).set_index("test")
)

print(
    f"Hedge ratio implied by the leading Johansen vector: {-johansen_vector[1] / johansen_vector[0]:.4f}"
)
```

This pair is the lesson rather than the demonstration. An energy sector fund and a crude oil
fund are related by something a reader can state in a sentence, and neither test calls them
cointegrated over this sample. The reason is visible in the funds: one holds equity in companies
whose earnings depend on the oil price, the other holds futures and pays to roll them, and a
roll cost accumulates as a drift that no fixed combination removes.

The implied hedge ratio printed above is a good place to see what a rejected test means. It is
far from the regression's, and on a pair where no stationary combination is found the leading
eigenvector is not estimating one: it is the least badly behaved direction in a system that has
no well behaved one, and its ratio is not a number to hold a position on.

The mechanics below are worth working through anyway, because every step still computes and
every step's output says something about a pair that fails.

## Two hedge ratios

The spread needs a ratio: how many units of one series to hold against one unit of the other.
A regression over the whole sample gives one number and a filter gives a series, and the
difference between them is not a matter of preference.

The regression is fitted on everything, which makes its ratio unavailable at every date
except the last. It is reported because it is the number the cointegration test used, and its
spread is what the stationarity test below is run on.

The **Kalman filter** treats the intercept and the ratio as a two-dimensional state following
a random walk and updates it one session at a time. Its estimate at session $t$ has seen
sessions up to $t$ and no more, which is what a hedge ratio has to be if a position is taken
on it. The two parameters that decide its behavior are the measurement noise and the process
noise, and their ratio is the whole tuning: a larger process noise lets the ratio move faster
and tracks a genuine structural change sooner, at the cost of chasing noise.

```python
MEASUREMENT_NOISE = 1e-3
PROCESS_NOISE = 1e-5

regression = LinearRegression().fit(pair[["independent"]], pair["dependent"])
static_ratio = float(regression.coef_[0])
pair["spread_static"] = pair["dependent"] - static_ratio * pair["independent"]

adf_statistic, adf_p_value = adfuller(pair["spread_static"], autolag="AIC")[:2]


def kalman_hedge_ratio(
    dependent: np.ndarray, independent: np.ndarray
) -> tuple[np.ndarray, np.ndarray]:
    """Filter an intercept and a hedge ratio forward, one session at a time."""
    filter_ = KalmanFilter(dim_x=2, dim_z=1)
    filter_.F = np.eye(2)
    filter_.x = np.array([[0.0], [1.0]])
    filter_.P = np.eye(2)
    filter_.R = np.array([[MEASUREMENT_NOISE]])
    filter_.Q = np.eye(2) * PROCESS_NOISE

    intercepts = np.empty(len(dependent))
    ratios = np.empty(len(dependent))
    for step in range(len(dependent)):
        filter_.H = np.array([[1.0, independent[step]]])
        filter_.predict()
        filter_.update(np.array([[dependent[step]]]))
        intercepts[step] = filter_.x[0, 0]
        ratios[step] = filter_.x[1, 0]
    return ratios, intercepts


pair["hedge_ratio"], _ = kalman_hedge_ratio(
    pair["dependent"].to_numpy(), pair["independent"].to_numpy()
)
pair["spread"] = pair["dependent"] - pair["hedge_ratio"] * pair["independent"]

print(f"Hedge ratio from a regression on the whole sample: {static_ratio:.4f}")
print(
    f"Its spread, augmented Dickey-Fuller: statistic {adf_statistic:.4f}, p value {adf_p_value:.4f}"
)
print(
    f"Filtered hedge ratio: from {pair['hedge_ratio'].min():.4f} to "
    f"{pair['hedge_ratio'].max():.4f}, ending at {pair['hedge_ratio'].iloc[-1]:.4f}"
)
```

## How fast the spread comes back

Fit an autoregression of order one to the spread and the coefficient says what fraction of a
deviation is still there one session later. Written as a change,

$$\Delta s_t = \alpha + \phi\, s_{t-1} + \varepsilon_t$$

a negative $\phi$ means the spread is pulled back toward its level, and the **half-life** is
the number of sessions over which half of a deviation is expected to have gone:
$-\ln 2 / \ln(1 + \phi)$. A non-negative $\phi$ says the spread is not being pulled anywhere
and there is no half-life to report.

The half-life is used below to choose the window the trading signal's mean and standard
deviation are computed over, which is why it is estimated on the first `TRAIN_FRACTION` of the
sample only. A window length chosen from the whole sample is a small leak and an easy one to
miss, because a length does not look like a parameter.

```python
def estimate_half_life(spread: pd.Series) -> float:
    """Sessions over which half a deviation decays, from an autoregression of order one."""
    level = spread.dropna()
    change = level.diff().dropna()
    lagged = level.shift(1).dropna().loc[change.index]

    coefficient = float(
        LinearRegression().fit(lagged.to_numpy().reshape(-1, 1), change.to_numpy()).coef_[0]
    )
    if coefficient >= 0.0:
        return float("inf")
    return float(-np.log(2.0) / np.log1p(coefficient))


half_life_train = estimate_half_life(pair["spread"].iloc[:train_end])
half_life_full = estimate_half_life(pair["spread"])
LOOKBACK = max(int(2 * half_life_train), MINIMUM_LOOKBACK)

print(f"Half-life on the first block: {half_life_train:.1f} sessions")
print(f"Half-life on the whole sample, for comparison only: {half_life_full:.1f} sessions")
print(f"Signal window, twice the first block's half-life: {LOOKBACK} sessions")
```

The two half-lives differ, and the gap between them is how far the window length would have
moved. Neither is more correct as an estimate. The first is the only one available at the
moment the window has to be chosen.

## From spread to position

The signal is the spread measured against its own recent history: subtract a rolling mean,
divide by a rolling standard deviation, and the result is in units of its own variability. A
spread `ENTRY_THRESHOLD` standard deviations below its recent mean is an entry on the long
side, meaning long the dependent series and short `hedge_ratio` units of the other, and the
position is held until the spread crosses back through the mean.

Both the mean and the standard deviation are rolling and both use the window chosen above, so
every value is available on its own session. The position is stepped forward one session at a
time rather than computed with a vector operation, because it depends on its own previous
value: an entry that has not yet reverted stays open.

```python
rolling_mean = pair["spread"].rolling(LOOKBACK).mean()
rolling_deviation = pair["spread"].rolling(LOOKBACK).std()
pair["spread_mean"] = rolling_mean
pair["upper_band"] = rolling_mean + ENTRY_THRESHOLD * rolling_deviation
pair["lower_band"] = rolling_mean - ENTRY_THRESHOLD * rolling_deviation
pair["z_score"] = (pair["spread"] - rolling_mean) / rolling_deviation

scores = pair["z_score"].to_numpy()
position = np.zeros(len(pair), dtype=np.int64)
held = 0
for step in range(1, len(pair)):
    if step < train_end:
        continue
    if scores[step] < -ENTRY_THRESHOLD:
        held = 1
    elif scores[step] > ENTRY_THRESHOLD:
        held = -1
    elif held * scores[step] > 0:
        held = 0
    position[step] = held

pair["position"] = position

held_positions = pair["position"].to_numpy()
previous_positions = np.concatenate([[0], held_positions[:-1]])
entries = int(((held_positions != previous_positions) & (held_positions != 0)).sum())

print(f"Sessions with a position open: {int((held_positions != 0).sum()):,}")
print(f"Entries, counting a reversal as a new one: {entries:,}")
```

```python
fig, axes = plt.subplots(4, 1, figsize=FIGSIZE["grid_3x2"], sharex=True)

ax = axes[0]
ax.plot(pair.index, pair["dependent"], linewidth=0.8, color=COLORS["blue"], label=DEPENDENT)
ax.plot(pair.index, pair["independent"], linewidth=0.8, color=COLORS["copper"], label=INDEPENDENT)
ax.set_ylabel("US dollars")
ax.set_title("Two funds with a reason to move together", fontsize=9)
ax.legend(fontsize=7)

ax = axes[1]
ax.plot(pair.index, pair["hedge_ratio"], linewidth=0.9, color=COLORS["blue"], label="filtered")
ax.axhline(static_ratio, color=COLORS["amber"], linestyle="--", linewidth=0.9, label="whole sample")
ax.set_ylabel("Units held short")
ax.set_title("The filtered hedge ratio moves; the regression cannot", fontsize=9)
ax.legend(fontsize=7)

ax = axes[2]
ax.plot(pair.index, pair["spread"], linewidth=0.8, color=COLORS["blue"])
ax.plot(pair.index, pair["spread_mean"], linewidth=0.8, color=COLORS["amber"], linestyle="--")
ax.fill_between(
    pair.index, pair["lower_band"], pair["upper_band"], alpha=0.2, color=COLORS["silver_muted"]
)
ax.set_ylabel("US dollars")
ax.set_title("The spread against a band of its own recent variability", fontsize=9)

ax = axes[3]
ax.plot(pair.index, pair["z_score"], linewidth=0.8, color=COLORS["blue"])
for level in (-ENTRY_THRESHOLD, ENTRY_THRESHOLD):
    ax.axhline(level, color=COLORS["negative"], linestyle="--", linewidth=0.7)
ax.axhline(0, color=COLORS["recede"], linewidth=0.7)
ax.axvline(pair.index[train_end], color=COLORS["neutral"], linestyle=":", linewidth=0.9)
ax.set_ylabel("Standard deviations")
ax.set_xlabel("Session")
ax.set_title("Entries at the dashed lines, exits at the crossing", fontsize=9)

fig.suptitle(f"{DEPENDENT} against {INDEPENDENT}: spread, band and signal")
show_with_alt(
    fig,
    "Four stacked panels for one pair of funds. The top plots both price series. The second "
    "plots the filtered hedge ratio as a moving line against a flat dashed line for the "
    "whole-sample regression, and the two are far apart for long stretches. The third plots the "
    "spread with a shaded band around a dashed rolling mean. The bottom plots the spread's "
    "standardized deviation with dashed entry lines above and below zero and a dotted vertical "
    "line where the estimation block ends.",
)
```

## What the position would have returned

What was held over session $t$ is one share of the dependent fund against $\beta_{t-1}$ shares
of the other, because $\beta_{t-1}$ is the ratio the filter had produced by the previous close.
Its profit in dollars is the two price changes weighted by those holdings:

$$\Delta P_t - \beta_{t-1}\, \Delta Q_t$$

and to express it as a return it is divided by what holding both legs committed, which is the
value of one share of the first plus $|\beta_{t-1}|$ shares of the second.

That is deliberately not the first difference of the `spread` column. The spread at $t$ is
defined with $\beta_t$, so differencing it mixes the session's price moves with the change in
the ratio itself, and a holder cannot rebalance at a ratio they only learn at that session's
close. The two agree only where the filtered ratio is flat.

It is also not the difference of the two funds' returns. That quantity is the profit of equal
dollar amounts in the two legs, which is a hedge in money rather than in shares, and it uses no
estimated ratio at all. Which of the two a strategy wants is a choice; reporting one while
having estimated the other is not.

Three limits, all structural rather than fixable here. No transaction costs and no slippage are
charged, and a spread strategy trades often enough to be sensitive to both. The short leg is
assumed available at no borrowing cost. And the pair was chosen by looking at the same history,
which the screening section takes up.

```python
evaluated = pair.index >= pair.index[train_end]

profit = pair["position"].shift(1) * (
    pair["dependent"].diff() - pair["hedge_ratio"].shift(1) * pair["independent"].diff()
)
committed = pair["dependent"].shift(1) + pair["hedge_ratio"].shift(1).abs() * pair[
    "independent"
].shift(1)
pair["strategy_return"] = (profit / committed).fillna(0.0)

realized = pair.loc[evaluated, "strategy_return"]
curve = (1.0 + realized).cumprod()

display(
    pd.DataFrame(
        [
            {
                "sessions evaluated": int(len(realized)),
                "sessions with a position": int((pair.loc[evaluated, "position"] != 0).sum()),
                "total return": float(curve.iloc[-1] - 1.0),
                "annualized mean return": float(realized.mean() * SESSIONS_PER_YEAR),
                "annualized volatility": float(realized.std() * np.sqrt(SESSIONS_PER_YEAR)),
                "deepest drawdown": float((curve / curve.cummax() - 1.0).min()),
            }
        ]
    ).T.rename(columns={0: DEPENDENT + " against " + INDEPENDENT})
)
```

The entry count is the number to read, and it is the section's conclusion. A half-life of
months sets a window of about a year, a band that wide is crossed rarely, and the result is a
position that opens a couple of times across the evaluated years and is then held for more than
half of their sessions. That is not
a mean-reversion strategy; it is a slow directional bet on a spread, which is what the
cointegration tests said would happen when they declined to find a stationary combination.

The return is therefore not evidence about the method. It is one draw from as many trades as the
entry count reports. The
drawdown is the one figure here worth carrying: it is what a position of that duration exposed
a holder to while waiting for a reversion the spread had no mechanism to deliver.

## Screening, and what a screen costs

Six candidate pairs, each with a stated reason to be related, run through the same two tests.
The table is the point of the section: a reason is not evidence, and the two tests do not
always agree.

The screen also carries a cost that the table cannot show. Testing six pairs at the same
level and reporting whichever passes makes the reported level wrong, because the chance that
at least one of six independent tests rejects a true null is far above the level of any one of
them. What follows from that is not a correction to apply here but a rule for reading: a pair
that passes a screen has cleared a lower bar than a pair tested on its own, and its spread
needs to hold up on data the screen did not see.

```python
CANDIDATE_PAIRS = [
    ("GLD", "SLV", "gold against silver"),
    ("XLE", "USO", "energy equity against crude oil"),
    ("QQQ", "SMH", "the Nasdaq 100 against semiconductors"),
    ("SPY", "VTI", "the S&P 500 against the total market"),
    ("TLT", "IEF", "long against intermediate Treasuries"),
    ("EEM", "VWO", "two emerging market funds"),
]

screen_rows = []
for dependent, independent, description in CANDIDATE_PAIRS:
    if not {dependent, independent} <= available:
        continue
    candidate = pair_frame(dependent, independent)
    if len(candidate) < MINIMUM_SESSIONS:
        continue

    candidate_p_value = coint(candidate["dependent"], candidate["independent"], trend="c")[1]
    candidate_johansen = coint_johansen(
        candidate[["dependent", "independent"]].to_numpy(), det_order=0, k_ar_diff=1
    )
    ratio = float(
        LinearRegression().fit(candidate[["independent"]], candidate["dependent"]).coef_[0]
    )
    spread = candidate["dependent"] - ratio * candidate["independent"]

    screen_rows.append(
        {
            "pair": f"{dependent} / {independent}",
            "why it might hold": description,
            "sessions": len(candidate),
            "Engle-Granger p value": candidate_p_value,
            "Johansen rejects": bool(candidate_johansen.lr1[0] > candidate_johansen.cvt[0, 1]),
            "hedge ratio": ratio,
            "half-life in sessions": estimate_half_life(spread),
        }
    )

screen = pd.DataFrame(screen_rows).set_index("pair")
display(screen)

print(
    f"Pairs both tests reject no cointegration for: {int(((screen['Engle-Granger p value'] < SIGNIFICANCE) & screen['Johansen rejects']).sum())}"
)
print(
    f"Pairs exactly one test rejects for: {int(((screen['Engle-Granger p value'] < SIGNIFICANCE) ^ screen['Johansen rejects']).sum())}"
)
```

Read the two counts printed under the table first. Over this sample no pair clears the
Engle-Granger test at the stated level, and one pair is called cointegrated by Johansen alone.
Six pairs with a stated reason, one weak signal: that is the result, and it is more useful than
a table of successes would have been, because the reasons were not bad ones.

The half-life column decides something different from whether the relationship exists. A
half-life of a few sessions leaves room to enter and exit inside the reversion. Every half-life
here is measured in months, and the shortest of them still holds a position through a quarter in
which the relationship may have changed for reasons no test on past data can see. A pair needs
both: a spread that reverts, and reversion fast enough to trade.

## The same features as library expressions

The pairwise work above is written out with statsmodels and numpy, one pair at a time.
`ml4t-engineer` supplies the rolling versions as Polars expressions, which is the form a
pipeline over many pairs wants: `rolling_correlation` and `beta_to_market` read return
columns, `co_integration_score` reads price columns, and `correlation_regime_indicator` reads
a correlation column and emits flags for its level.

What these are not is a replacement for the tests above. A rolling score computed over a
window is a description of that window, and cointegration is a statement about a relationship
holding over a long sample. The two answer different questions.

```python
LIBRARY_PAIR = ("GLD", "SLV") if {"GLD", "SLV"} <= available else (DEPENDENT, INDEPENDENT)

library = (
    pair_frame(*LIBRARY_PAIR)
    .pipe(pl.from_pandas, include_index=True)
    .rename({"dependent": "first_close", "independent": "second_close"})
    .with_columns(
        first_return=pl.col("first_close").pct_change(),
        second_return=pl.col("second_close").pct_change(),
    )
    .drop_nulls()
    .with_columns(
        correlation=rolling_correlation("first_return", "second_return", window=CORRELATION_WINDOW),
        beta=beta_to_market("second_return", "first_return", window=CORRELATION_WINDOW),
        cointegration_score=co_integration_score(
            "first_close", "second_close", window=COINTEGRATION_WINDOW
        ),
    )
)
library = library.with_columns(**correlation_regime_indicator("correlation"))

print(f"Library pair: {LIBRARY_PAIR[0]} against {LIBRARY_PAIR[1]}, {library.height:,} sessions")
display(
    library.select(["correlation", "beta", "cointegration_score"])
    .describe()
    .filter(pl.col("statistic").is_in(["mean", "std", "min", "max"]))
)
display(library.select(pl.col("^corr_regime.*$")).mean())
```

## A level replaced by a position

The second half of the notebook changes what a feature is computed over. Nothing here is a
rolling window: the operation reads one date across many assets and returns each asset's
standing among them on that date.

Three versions of the same idea, and they are not interchangeable. A **rank** is an integer
position and is bounded by the number of assets, so it says nothing about how far apart the
assets are. A **percentile** divides the rank by the count, which makes a universe of ten
comparable with one of five hundred. A **z-score** subtracts the cross-sectional mean and
divides by the cross-sectional standard deviation, which keeps the distances and therefore
keeps the outliers, and on a date when one asset moved five times as far as the rest it is the
only one of the three that says so.

None of the three uses a future date. That is worth stating because the causality question
looks the same as everywhere else in the chapter and has a different answer: a cross-sectional
operation reads across assets at one date, so it cannot look forward in time whatever else it
does.

```python
RANK_SYMBOLS = ["SPY", "QQQ", "IWM", "EFA", "EEM", "TLT", "GLD", "XLE", "XLF", "XLV"]
selected = [symbol for symbol in RANK_SYMBOLS if symbol in available]

panel = (
    universe.filter(pl.col("symbol").is_in(selected))
    .select(["timestamp", "symbol", "close"])
    .sort(["symbol", "timestamp"])
    .with_columns(returns=pl.col("close").pct_change().over("symbol"))
    .with_columns(
        momentum=pl.col("returns").rolling_mean(FEATURE_WINDOW).over("symbol"),
        volatility=pl.col("returns").rolling_std(FEATURE_WINDOW).over("symbol")
        * np.sqrt(SESSIONS_PER_YEAR),
    )
    .drop_nulls()
    .with_columns(
        momentum_rank=pl.col("momentum").rank().over("timestamp"),
        volatility_rank=pl.col("volatility").rank().over("timestamp"),
        momentum_percentile=(pl.col("momentum").rank().over("timestamp") - 1)
        / (pl.col("momentum").count().over("timestamp") - 1),
        momentum_z_score=(pl.col("momentum") - pl.col("momentum").mean().over("timestamp"))
        / pl.col("momentum").std().over("timestamp"),
    )
)

print(f"Panel: {panel['symbol'].n_unique()} symbols, {panel.height:,} rows")
snapshot_date = panel["timestamp"].unique().sort().to_list()[-REGIME_WINDOW]
display(
    panel.filter(pl.col("timestamp") == snapshot_date)
    .select(["symbol", "momentum", "momentum_rank", "momentum_percentile", "momentum_z_score"])
    .sort("momentum_rank")
)
```

Read the snapshot's last two columns against each other. The percentiles are evenly spaced by
construction, one step per asset, whatever the momentum values were. The z-scores are not, and
where they bunch together the percentile has manufactured a distinction between assets that
were nearly identical. That is the trade: a rank is robust because it discards the distances,
and it is misleading for the same reason.

```python
fig, axes = plt.subplots(2, 1, figsize=FIGSIZE["dual_v"], sharex=True)
HIGHLIGHT = [symbol for symbol in ("SPY", "GLD", "XLE") if symbol in selected]
palette = [COLORS["blue"], COLORS["amber"], COLORS["copper"]]

for symbol, color in zip(HIGHLIGHT, palette, strict=False):
    series = panel.filter(pl.col("symbol") == symbol).to_pandas()
    axes[0].plot(
        series["timestamp"],
        series["momentum_percentile"],
        linewidth=0.8,
        color=color,
        label=symbol,
    )
    axes[1].plot(
        series["timestamp"], series["volatility_rank"], linewidth=0.8, color=color, label=symbol
    )

axes[0].axhline(0.5, color=COLORS["recede"], linestyle="--", linewidth=0.7)
axes[0].set_ylabel("Percentile")
axes[0].set_title("Momentum measured against the rest of the universe", fontsize=9)
axes[0].legend(fontsize=7)

axes[1].set_ylabel("Rank")
axes[1].set_xlabel("Session")
axes[1].set_title("Volatility rank, where one is the calmest of the universe", fontsize=9)

fig.suptitle("A cross-sectional feature is bounded whatever the market does")
show_with_alt(
    fig,
    "Two stacked panels tracking three funds. The top plots each fund's momentum percentile "
    "against a dashed line at one half; the lines travel the full range and cross often. The "
    "bottom plots each fund's volatility rank, which steps between integer levels and holds "
    "for months at a time.",
)
```

The rank in the lower panel holds for long stretches and then steps, which is a different
statistical object from the series it was computed from: a rolling standard deviation moves
every session, and its rank moves only when the ordering changes. A downstream model reading
the rank sees a step function, and one reading the raw volatility sees a continuous series.
Which is wanted depends on whether the question is how volatile an asset is or which assets
are the volatile ones.

## A level replaced by a difference from a benchmark

The other way to place a feature in context is to subtract a benchmark's version of it, which
keeps the units. Momentum minus the market's momentum is a momentum, in the same units, with
the part every asset shared removed. Volatility divided by the market's volatility is a
ratio, and above one says the asset moved more than the market did.

How the benchmark enters the feature decides which of the two operations to use. A momentum is
a sum of returns and the market component is additive in it, so subtraction removes it. A
volatility is a scale and the market's stress multiplies it, so a ratio removes it and a
difference does not.

```python
BENCHMARK = "SPY" if "SPY" in selected else selected[0]
COMPARED = "XLE" if "XLE" in selected else selected[-1]

benchmark = (
    panel.filter(pl.col("symbol") == BENCHMARK)
    .select(["timestamp", "momentum", "volatility"])
    .rename({"momentum": "market_momentum", "volatility": "market_volatility"})
)

relative = panel.join(benchmark, on="timestamp", how="inner").with_columns(
    momentum_less_market=pl.col("momentum") - pl.col("market_momentum"),
    volatility_over_market=pl.col("volatility") / pl.col("market_volatility"),
)

compared = relative.filter(pl.col("symbol") == COMPARED).to_pandas()

fig, axes = plt.subplots(2, 1, figsize=FIGSIZE["dual_v"], sharex=True)

axes[0].plot(
    compared["timestamp"],
    compared["momentum"],
    linewidth=0.8,
    color=COLORS["blue"],
    label="its own",
)
axes[0].plot(
    compared["timestamp"],
    compared["momentum_less_market"],
    linewidth=0.8,
    color=COLORS["copper"],
    label="less the market's",
)
axes[0].axhline(0, color=COLORS["recede"], linestyle="--", linewidth=0.7)
axes[0].set_ylabel("Mean daily return")
axes[0].set_title("Subtracting the market leaves what was specific to the fund", fontsize=9)
axes[0].legend(fontsize=7)

axes[1].plot(
    compared["timestamp"],
    compared["volatility_over_market"],
    linewidth=0.8,
    color=COLORS["blue"],
)
axes[1].axhline(1.0, color=COLORS["recede"], linestyle="--", linewidth=0.7)
axes[1].set_ylabel("Ratio")
axes[1].set_xlabel("Session")
axes[1].set_title("Volatility as a multiple of the market's, not a difference", fontsize=9)

fig.suptitle(f"{COMPARED} against {BENCHMARK}, in units that survive the comparison")
show_with_alt(
    fig,
    f"Two stacked panels for {COMPARED}. The top plots its own momentum and its momentum less "
    f"{BENCHMARK}'s, both against a dashed zero line; the two diverge most where the market "
    "itself moved. The bottom plots the ratio of its volatility to the market's against a "
    "dashed line at one, and the ratio stays above one for most of the sample.",
)
```

The lower panel is the one that changes a reading. A ratio near one during a market-wide
selloff says this fund was no more volatile than everything else, which the absolute
volatility could not have said: that number was high because every number was high.

## One measure aggregated across the universe

The last construction turns a per-asset regime measure into universe-level features. Three
aggregates, each answering a different question about the same date:

- The **mean** stress across assets, which is the level.
- The **breadth**, the fraction of assets above a threshold, which distinguishes a few assets
  in trouble from all of them.
- The **dispersion**, the standard deviation across assets, which is low when everything moves
  together and high when the universe has split.

The per-asset measure here is each asset's volatility ranked within a trailing window of
`REGIME_WINDOW` sessions, a stand-in for the filtered regime probability that
`11_hmm_regimes` produces. `rolling_rank` is what makes it a stand-in rather than a leak:
ranking each asset's volatility over its whole history would place today's value against
values from years that had not happened, and the resulting series would look like a regime
indicator while being unavailable on every date but the last.

```python
STRESS_THRESHOLD = 0.7

stress = (
    panel.select(["timestamp", "symbol", "volatility"])
    .sort(["symbol", "timestamp"])
    .with_columns(
        stress=pl.col("volatility").rolling_rank(window_size=REGIME_WINDOW).over("symbol")
        / REGIME_WINDOW
    )
    .drop_nulls()
)

aggregates = (
    stress.group_by("timestamp")
    .agg(
        mean_stress=pl.col("stress").mean(),
        breadth=(pl.col("stress") > STRESS_THRESHOLD).mean(),
        dispersion=pl.col("stress").std(),
        assets=pl.len(),
    )
    .sort("timestamp")
)

benchmark_level = (
    universe.filter(pl.col("symbol") == BENCHMARK)
    .select(["timestamp", "close"])
    .sort("timestamp")
    .with_columns(cumulative=(pl.col("close") / pl.col("close").first() - 1.0) * 100)
)

aggregated = aggregates.join(
    benchmark_level.select(["timestamp", "cumulative"]), on="timestamp", how="inner"
).to_pandas()

print(
    f"Aggregate rows: {len(aggregated):,}, assets per date: {aggregates['assets'].min()} to {aggregates['assets'].max()}"
)
display(aggregates.select(["mean_stress", "breadth", "dispersion"]).describe())
```

```python
fig, axes = plt.subplots(3, 1, figsize=FIGSIZE["grid_3x2"], sharex=True)

axes[0].plot(
    aggregated["timestamp"], aggregated["cumulative"], linewidth=0.8, color=COLORS["neutral"]
)
axes[0].set_ylabel("Percent")
axes[0].set_title(f"{BENCHMARK} cumulative return, for context", fontsize=9)

axes[1].fill_between(
    aggregated["timestamp"], 0, aggregated["breadth"], alpha=0.4, color=COLORS["copper"]
)
axes[1].set_ylabel("Fraction")
axes[1].set_title("Breadth: how much of the universe is in its own high-stress range", fontsize=9)

axes[2].plot(aggregated["timestamp"], aggregated["dispersion"], linewidth=0.8, color=COLORS["blue"])
axes[2].set_ylabel("Standard deviation")
axes[2].set_xlabel("Session")
axes[2].set_title("Dispersion: whether the universe agrees about the stress", fontsize=9)

fig.suptitle("Three aggregates of one per-asset measure, each saying something different")
show_with_alt(
    fig,
    f"Three stacked panels. The top plots {BENCHMARK}'s cumulative return. The middle fills the "
    "fraction of the universe in its own high-stress range, which rises to near one during the "
    "sharp declines above it. The bottom plots the cross-sectional dispersion of the stress "
    "measure, which falls when the breadth is at its highest.",
)
```

Breadth and dispersion move against each other at the extremes, and that is the reason to
carry both. When everything is stressed the breadth is near one and there is nothing left for
the assets to disagree about, so the dispersion falls to near zero. A high dispersion is a
universe where some assets are under stress and others are not, which is a different
environment and not a milder version of the same one.

Both are computed from a trailing rank, so both are available on the session they describe.
What neither is, is a regime probability: a rank says where today sits among the last
`REGIME_WINDOW` sessions of the same asset, and a fitted model says how likely a state is
given everything observed. The aggregation mechanics are the same either way, which is what
this section is for.

## Takeaways

1. **Cointegration is a claim about drift, not about co-movement.** Two funds can track each
   other daily and still fail both tests, and one of the two funds carrying a roll cost is
   enough to do it. A stated economic reason is where a screen starts and not evidence: six
   pairs with good reasons produced no Engle-Granger rejection at the level used here.
2. **Engle-Granger and Johansen can disagree, and the disagreement is information.** One takes
   a dependent variable and one treats the pair as a system. A pair only one of them rejects
   for is a pair whose relationship depends on how the question was asked.
3. **A hedge ratio used to hold a position has to be filtered.** A regression over the whole
   sample gives one number that was unavailable on every date but the last. The Kalman
   recursion gives a series, and its process noise is the choice between tracking a real
   change and chasing noise.
4. **A window length is a parameter and leaks like one.** The half-life that sizes the signal
   window is estimated on the first block only, because a length chosen from the whole sample
   is a leak that does not look like one.
5. **Price the position that was held, not the column that was plotted.** The profit is the two
   price changes weighted by the previous session's holdings. Differencing a spread whose ratio
   moves adds the ratio's own change to it, and differencing the two funds' returns prices equal
   dollar amounts in each leg, which is a hedge in money and not the share ratio estimated.
6. **Rank, percentile and z-score discard different things.** Ranks and percentiles throw away
   the distances between assets, which is what makes them robust and what makes them invent
   distinctions between assets that were nearly identical.
7. **A cross-sectional operation cannot look forward, and a rank over an asset's own history
   can.** The first reads one date across assets. The second is the same word applied along
   time, where it needs a trailing window like every other temporal feature.

**Previous**: `13_regime_as_feature` puts a regime column into a downstream model.
**Next**: `case_study_temporal_summary` collects what this chapter's notebooks produced.
![notebook output](figures/p1_1.png)
![notebook output](figures/p1_2.png)
![notebook output](figures/p1_3.png)
![notebook output](figures/p1_4.png)