가우시안 혼합 모형으로 요인 수익률의 시장 레짐 발견

# Factor-based regime detection


# Factor-based regime detection

**Chapter 1 · §1.4 Keeping up with changing market regimes**

**Docker image**: `ml4t`

Markets do not behave the same way in every year. The same value strategy that earns a
steady premium through a long expansion can lose money for a decade, and the same
momentum strategy can give back years of gains in a single quarter. A **regime** is a
stretch of time over which the joint behaviour of returns is stable enough to be treated
as one environment: the averages, the volatilities and the correlations hold roughly
still inside it and shift when it ends.

Nobody publishes a regime label. This notebook estimates one from the data with a
**Gaussian mixture model** (GMM), which treats each month's vector of factor returns as
a draw from one of a small number of multivariate normal distributions and infers both
the distributions and which one each month most likely came from. The inputs are monthly
returns to nine factor strategies from AQR's Century of Factor Premia dataset, which
starts in 1927.

## Learning objectives

By the end of this notebook you will be able to:

- Fit a Gaussian mixture model to a panel of monthly factor returns and read off the
  cluster it assigns to each month.
- Choose how many clusters to keep by comparing three criteria that measure different
  things, and say what each one rewards.
- Read a regime timeline against dated market events and judge whether the clusters
  correspond to periods a reader would recognise.
- Measure what each factor earned inside each cluster, and say what that implies about
  which factors diversify each other when conditions deteriorate.
- Explain why a label fitted on the whole sample cannot be used as a trading signal.

## Book reference

Chapter 1, Section 1.4, "Keeping up with changing market regimes". `macro_regimes` is the
companion notebook and estimates regimes from macroeconomic indicators instead.

## Prerequisites

- Monthly return series and standardization (subtracting a mean, dividing by a standard
  deviation) so that series with different scales contribute comparably to a distance.
- The AQR Century of Factor Premia dataset, which `AQRFactorProvider` reads from
  `data/factors/aqr` under `ML4T_DATA_PATH`. Fetch it with
  `uv run python data/factors/aqr_download.py` if it is missing.

## What this notebook establishes, and what it does not

The model is fitted on the whole history and the labels are assigned back over that same
history. What comes out is a description of how the factor space partitions, and a way
of seeing where dated events fall inside that partition. Every month's label is informed
by every other month, including later ones, so the labels describe the sample rather than
forecast it. Building regime labels a strategy could have acted on requires walk-forward
fitting on information available at each date; Chapter 6 onward introduces that machinery
and the case-study chapters apply it.

## Setup

The AQR provider logs its download and cache decisions through `structlog`, which writes
to the notebook rather than through the standard library's logging module. Raising its
threshold to `WARNING` keeps the data-loading cells readable while leaving anything that
reports a problem visible.

```python
"""Factor-based regime detection - Gaussian mixture regimes over AQR factor returns."""

from __future__ import annotations

import logging
from collections.abc import Iterable
from dataclasses import dataclass

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import polars as pl
import structlog
from IPython.display import Markdown, display
from matplotlib.axes import Axes
from matplotlib.collections import LineCollection
from matplotlib.colors import ListedColormap
from ml4t.data.providers import AQRFactorProvider
from ml4t.diagnostic.metrics import sharpe_ratio
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score
from sklearn.mixture import GaussianMixture
from sklearn.preprocessing import StandardScaler

import utils.style as style
from utils.paths import get_output_dir
from utils.reproducibility import set_global_seeds

structlog.configure(wrapper_class=structlog.make_filtering_bound_logger(logging.WARNING))

COLORS = style.COLORS
```

```python
# Production defaults (Papermill overrides these when the test suite runs the notebook)
SEED = 42
```

### Settings, and what each one decides

`SEED` fixes the initialization of every estimator below. A Gaussian mixture is fitted by
expectation-maximization from a random start, so two runs from different starts can land
on different local optima and swap which cluster is numbered 0; fixing the seed is what
makes the figures in this notebook reproduce.

`N_REGIMES_GRID` is the set of cluster counts to fit and compare. Two is the smallest
number that partitions anything. Six is where the search stops: the point of the sweep is
to show how the criteria move as clusters are added, and by six they have already
separated.

`MONTHS_PER_YEAR` annualizes monthly statistics. Multiply a mean monthly return by it and
a monthly standard deviation by its square root.

`ROLLING_VOL_MONTHS` is the width of the moving window used to draw a volatility series.
Twelve months is short enough to move inside a regime and long enough that a single month
does not dominate the estimate.

```python
N_REGIMES_GRID = [2, 3, 4, 5, 6]
MONTHS_PER_YEAR = 12
ROLLING_VOL_MONTHS = 12

OUTPUT_DIR = get_output_dir(1, "factor_regimes")
OUTPUT_DIR.mkdir(parents=True, exist_ok=True)

set_global_seeds(SEED)
```

## Load the factor panel

The AQR Century of Factor Premia dataset holds monthly returns to long-short factor
strategies built inside several asset classes, plus the returns to holding each asset
class outright. The provider returns one row per month and one column per strategy.

```python
aqr_raw = AQRFactorProvider().fetch("century_premia")

date_range = aqr_raw.select(
    pl.col("timestamp").min().alias("first"), pl.col("timestamp").max().alias("last")
)
print(f"AQR Century of Factor Premia: {aqr_raw.height} months, {aqr_raw.width - 1} series")
print(f"Covering {date_range['first'][0]:%Y-%m} to {date_range['last'][0]:%Y-%m}")
```

### The nine series used here

Four of them are cross-asset long-short strategies, applied across equity indices, bonds,
currencies and commodities at once:

- **Value** buys what is cheap relative to a fundamental anchor and sells what is
  expensive.
- **Momentum** buys what has risen over the past year and sells what has fallen.
- **Carry** buys what pays a high yield and sells what pays a low one.
- **Defensive** buys low-risk assets and sells high-risk ones.

Two apply the value and momentum definitions inside the US stock market alone. The last
three are the returns to simply holding each asset class - equity indices, fixed income,
commodities - rather than a long-short strategy, and they are what places a month in a
risk environment rather than a style environment.

```python
FACTOR_COLUMNS = [
    "All asset classes Value",
    "All asset classes Momentum",
    "All asset classes Carry",
    "All asset classes Defensive",
    "US Stock Selection Value",
    "US Stock Selection Momentum",
    "Equity indices Market",
    "Fixed income Market",
    "Commodities Market",
]
EQUITY_COLUMN = "Equity indices Market"

missing = [c for c in FACTOR_COLUMNS if c not in aqr_raw.columns]
if missing:
    raise ValueError(
        f"The AQR file is missing {missing}. Re-fetch it with "
        "`uv run python data/factors/aqr_download.py`."
    )

# Short names for the figures and tables below; the full AQR column names do not fit on an
# axis.
SHORT_NAMES = {
    "All asset classes Value": "Value",
    "All asset classes Momentum": "Momentum",
    "All asset classes Carry": "Carry",
    "All asset classes Defensive": "Defensive",
    "US Stock Selection Value": "US Value",
    "US Stock Selection Momentum": "US Momentum",
    "Equity indices Market": "Equity",
    "Fixed income Market": "Bonds",
    "Commodities Market": "Commodities",
}
```

### What is actually in the panel

Clustering treats every month as a point in nine dimensions, so a month is usable only
when all nine series are quoted. The table below reports, for each series, the month it
starts, how many observations it has, and its mean and standard deviation over the raw
file. Read it before anything computed from the panel: it says which series shortens the
usable history and how far apart the series are in scale, which is why they are
standardized before the distance is taken.

```python
coverage = pd.DataFrame(
    [
        {
            "Series": SHORT_NAMES[column],
            "First month": aqr_raw.filter(pl.col(column).is_not_null())["timestamp"].min(),
            "Months quoted": aqr_raw.select(pl.col(column).is_not_null().sum()).item(),
            "Mean (% / month)": aqr_raw.select(pl.col(column).mean()).item() * 100,
            "Std (% / month)": aqr_raw.select(pl.col(column).std()).item() * 100,
        }
        for column in FACTOR_COLUMNS
    ]
).set_index("Series")
coverage["First month"] = pd.to_datetime(coverage["First month"]).dt.strftime("%Y-%m")
coverage.style.format({"Mean (% / month)": "{:.2f}", "Std (% / month)": "{:.2f}"})
```

## Prepare the panel for clustering

Months in which any of the nine series is unquoted are dropped rather than filled. A
forward fill would repeat the previous month's return, and repeating a return is not a
missing observation recovered - it is an invented one, and a mixture model would treat it
as evidence for a cluster.

The nine series are then standardized to zero mean and unit variance. The coverage table
above shows why: the widest series in the panel has several times the monthly standard
deviation of the narrowest, and a Euclidean distance is dominated by whichever coordinate
moves furthest. Standardizing puts every series on the same footing, so the partition
reflects the direction the panel moved in rather than which column happened to be measured
on the largest scale.

```python
factors_pl = aqr_raw.select(["timestamp", *FACTOR_COLUMNS]).sort("timestamp").drop_nulls()

factors_df = factors_pl.select(FACTOR_COLUMNS).to_pandas()
factors_df.index = pd.DatetimeIndex(factors_pl["timestamp"].to_pandas())

factors_scaled = StandardScaler().fit_transform(factors_df)

start_year = factors_df.index.min().year
end_year = factors_df.index.max().year
print(f"Complete months: {len(factors_df)} ({start_year} to {end_year})")
```

## Fit a grid of mixture models

There is no test that returns the true number of regimes, so the count is a modelling
choice made by comparing candidates. The helper below fits one mixture per candidate
count and records three quantities for each.

**BIC**, the Bayesian information criterion, is the fitted log-likelihood penalised by the
number of free parameters times the log of the sample size. **AIC**, the Akaike
information criterion, applies the same idea with a constant penalty of two per parameter.
Both are better when lower, and BIC's penalty grows with the sample, so on a long panel it
is far more reluctant than AIC to buy an extra cluster. A full-covariance mixture spends
parameters quickly - each additional cluster costs a mean vector and a covariance matrix
over nine dimensions - which is what the two criteria are pricing differently.

The **silhouette score** measures something else entirely. For each month it compares the
average distance to the other months in its own cluster against the average distance to
the months of the nearest other cluster, and reports the normalized difference; the score
is the average over months. It runs from -1 to 1, it is better when higher, and a negative
value means the typical month sits closer to a different cluster than to its own. It is a
geometric statement about separation and knows nothing about likelihood.

K-means is fitted alongside for the same counts. It partitions on distance to a centroid
with no covariance and no probabilities, so putting its silhouette beside the mixture's
says how much of the separation is coming from the shape the mixture is allowed to fit.

The last two columns count things none of the three criteria look at: how often consecutive
months carry different labels, and how many months the emptiest cluster holds. A criterion
can only say how well a partition fits. Whether the partition holds still long enough to be
called a regime, and whether every cluster in it has enough months to say anything about,
are separate questions, and these two counts are the cheapest way to ask them.

```python
@dataclass(frozen=True)
class GmmFitResult:
    """One fitted mixture and the three quantities used to compare it against the others."""

    model: GaussianMixture
    labels: np.ndarray
    probabilities: np.ndarray
    bic: float
    aic: float
    silhouette: float
```

```python
def fit_gmm_grid(
    x: np.ndarray,
    n_components_list: Iterable[int],
    random_state: int = SEED,
) -> dict[int, GmmFitResult]:
    """Fit one full-covariance mixture per requested cluster count."""
    results: dict[int, GmmFitResult] = {}
    for n in n_components_list:
        model = GaussianMixture(
            n_components=n,
            covariance_type="full",
            random_state=random_state,
            n_init=10,
            reg_covar=1e-6,
        )
        model.fit(x)
        labels = model.predict(x)
        results[n] = GmmFitResult(
            model=model,
            labels=labels,
            probabilities=model.predict_proba(x),
            bic=float(model.bic(x)),
            aic=float(model.aic(x)),
            silhouette=float(silhouette_score(x, labels)),
        )
    return results
```

```python
gmm_grid = fit_gmm_grid(factors_scaled, N_REGIMES_GRID)

kmeans_labels = {
    n: KMeans(n_clusters=n, random_state=SEED, n_init=10).fit_predict(factors_scaled)
    for n in N_REGIMES_GRID
}


def switch_count(labels: np.ndarray) -> int:
    """How many times consecutive months carry different labels."""
    return int(np.sum(np.diff(labels) != 0))


selection = pd.DataFrame(
    {
        "Clusters": N_REGIMES_GRID,
        "BIC": [gmm_grid[n].bic for n in N_REGIMES_GRID],
        "AIC": [gmm_grid[n].aic for n in N_REGIMES_GRID],
        "Silhouette (mixture)": [gmm_grid[n].silhouette for n in N_REGIMES_GRID],
        "Silhouette (k-means)": [
            silhouette_score(factors_scaled, kmeans_labels[n]) for n in N_REGIMES_GRID
        ],
        "Switches": [switch_count(gmm_grid[n].labels) for n in N_REGIMES_GRID],
        "Smallest cluster (months)": [
            int(np.bincount(gmm_grid[n].labels).min()) for n in N_REGIMES_GRID
        ],
    }
).set_index("Clusters")
selection.style.format(
    {
        "BIC": "{:,.0f}",
        "AIC": "{:,.0f}",
        "Silhouette (mixture)": "{:.3f}",
        "Silhouette (k-means)": "{:.3f}",
    }
)
```

### Reading the three criteria against each other

The bars below draw three of the table's columns. Both information criteria are drawn on a
zoomed vertical axis, because the differences between candidate models are small relative
to their level and would be invisible on a scale that started at zero; the silhouette
panel is drawn on its natural scale, which includes zero.

```python
def plot_selection_criterion(ax: Axes, values: list[float], label: str, best: str) -> None:
    """Draw one criterion as bars, highlighting the count it favours."""
    bars = ax.bar(N_REGIMES_GRID, values, color=COLORS["recede"])
    chosen = int(np.argmin(values)) if best == "min" else int(np.argmax(values))
    bars[chosen].set_color(COLORS["amber"])
    ax.set_xlabel("Clusters")
    ax.set_ylabel(label)
    ax.set_xticks(N_REGIMES_GRID)
```

```python
fig, axes = plt.subplots(1, 3, figsize=style.FIGSIZE["triple_h_tall"])

plot_selection_criterion(axes[0], selection["BIC"].tolist(), "BIC (lower is better)", "min")
plot_selection_criterion(axes[1], selection["AIC"].tolist(), "AIC (lower is better)", "min")
plot_selection_criterion(
    axes[2], selection["Silhouette (mixture)"].tolist(), "Silhouette (higher is better)", "max"
)
for ax, values in ((axes[0], selection["BIC"]), (axes[1], selection["AIC"])):
    span = values.max() - values.min()
    ax.set_ylim(values.min() - 0.15 * span, values.max() + 0.15 * span)

style.add_message_title(
    axes[0],
    "BIC, AIC and silhouette do not pick the same number of regimes",
    subtitle="Highlighted bar is the count each criterion favours; note the zoomed axes on BIC and AIC",
)
style.show_with_alt(
    fig,
    "Three bar charts over cluster counts two to six. BIC rises with the count, AIC falls "
    "with it, and the silhouette score falls sharply after the smallest count and hovers "
    "near zero at the larger ones.",
)
```

```python
bic_choice = int(selection["BIC"].idxmin())
aic_choice = int(selection["AIC"].idxmin())
silhouette_choice = int(selection["Silhouette (mixture)"].idxmax())
display(
    Markdown(
        f"BIC is lowest at **{bic_choice} clusters** and the silhouette score is highest at "
        f"**{silhouette_choice}**, while AIC keeps falling out to **{aic_choice}**, the "
        "largest count fitted."
    )
)
```

When the criteria disagree, the choice is made on what the model is for. AIC targets
predictive likelihood and keeps buying clusters as long as they raise it, and the last two
columns say what it is buying at the larger counts: a partition that changes label more
often, containing at least one cluster with too few months to describe. A sliver of the
sample with its own covariance improves the fit and is not a regime anyone can name.

BIC's sample-scaled penalty and the silhouette's separation test both ask for clusters
distinct enough to survive being described, which is what a regime has to be if it is going
to appear in a risk report. The rest of this notebook works with the two-cluster model on
that basis, and the timelines in the next section show what the larger counts do with the
clusters they add.

```python
N_REGIMES_SELECTED = 2
# fit_gmm_grid fitted exactly the counts in N_REGIMES_GRID, and the swimlane and label
# cells below read gmm_grid[N_REGIMES_SELECTED]. Narrowing the grid without moving this
# constant reaches them as a KeyError.
assert N_REGIMES_SELECTED in N_REGIMES_GRID, (
    f"N_REGIMES_SELECTED={N_REGIMES_SELECTED} was never fitted; N_REGIMES_GRID={N_REGIMES_GRID}"
)
```

## What each cluster count partitions

Each band below is one fitted model. A month is shaded in the row of the cluster it was
assigned to, so a solid horizontal stretch is a period the model held in one regime and a
speckled stretch is one it switched through repeatedly. Comparing the bands shows what the
extra clusters buy: whether an added cluster carves out a recognisable era, or whether it
scatters across the whole century.

```python
def decade_ticks(dates: pd.DatetimeIndex) -> tuple[list[int], list[str]]:
    """Positions and labels for the first month of each decade in a monthly index."""
    years = dates.year.tolist()
    ticks = [
        j for j, year in enumerate(years) if year % 10 == 0 and (j == 0 or years[j - 1] != year)
    ]
    return ticks, [str(years[j]) for j in ticks]
```

```python
def draw_swimlanes(ax: Axes, rows: np.ndarray, n_rows: int, cmap: ListedColormap) -> None:
    """Shade each month in the row given by *rows*, one row per regime."""
    matrix = np.zeros((n_rows, len(rows)))
    matrix[rows, np.arange(len(rows))] = 1
    ax.imshow(
        matrix,
        aspect="auto",
        cmap=cmap,
        vmin=0,
        vmax=1,
        extent=(0, len(rows), -0.5, n_rows - 0.5),
        interpolation="nearest",
        origin="lower",
    )
    for boundary in range(1, n_rows):
        ax.axhline(y=boundary - 0.5, color=COLORS["bg_light"], linewidth=1.5)
```

```python
lane_cmap = ListedColormap([COLORS["bg_light"], COLORS["blue"]])
ticks, tick_labels = decade_ticks(factors_df.index)

fig, axes = plt.subplots(
    len(N_REGIMES_GRID),
    1,
    figsize=(style.PAGE_WIDTH, 1.0 * len(N_REGIMES_GRID) + 1.0),
    sharex=True,
)
for ax, n in zip(axes, N_REGIMES_GRID):
    draw_swimlanes(ax, gmm_grid[n].labels, n, lane_cmap)
    ax.set_yticks([])
    ax.set_ylabel(f"K = {n}", rotation=0, ha="right", va="center")

axes[-1].set_xticks(ticks)
axes[-1].set_xticklabels(tick_labels)
axes[-1].set_xlabel("Year")
style.add_message_title(
    axes[0],
    "Adding clusters makes the timeline switch more often, not cleaner",
    subtitle="One panel per fitted model; a month is shaded in the row of the cluster it was assigned",
)
style.show_with_alt(
    fig,
    "Five stacked panels, one per cluster count. The two-cluster panel shows one row "
    "occupied most of the time and the other taking short scattered stretches. Each panel "
    "below it is more finely speckled than the one above, with no row holding a long "
    "uninterrupted stretch by the largest count.",
)
```

No band settles down as clusters are added. An added cluster does not carve out an era the
earlier models had merged; it takes months from across the whole century, so each panel is
more finely speckled than the one above it. A mixture has no notion of time and nothing in
it prefers a month to keep its neighbour's label, so there is no mechanism by which more
clusters could produce longer stretches. That absence is the reason the last section of
the section on how long the regimes last points at a hidden Markov model, which supplies
exactly the mechanism this one lacks.

## The two-cluster model: risk-on and risk-off

The mixture returns cluster numbers, not names, and the numbering carries no meaning - a
different seed can swap it. To describe the two clusters they have to be identified by
something outside the model, and the natural choice is the series that says what the
market as a whole was doing: the average return to holding equity indices inside each
cluster. **Risk-on** and **risk-off** are the market's ordinary terms for the environments
that result, meaning periods when investors were adding exposure to risky assets and
periods when they were shedding it.

Naming the clusters this way uses the whole sample, exactly as fitting them did. It is a
description of what the partition turned out to be, not a rule that could have been
applied in 1929.

```python
labels_2 = gmm_grid[N_REGIMES_SELECTED].labels
equity_returns = factors_df[EQUITY_COLUMN]

mean_equity_by_cluster = equity_returns.groupby(labels_2).mean()
risk_on_cluster = int(mean_equity_by_cluster.idxmax())
risk_off_cluster = int(mean_equity_by_cluster.idxmin())

regime_name = pd.Series(
    np.where(labels_2 == risk_on_cluster, "Risk-on", "Risk-off"),
    index=factors_df.index,
    name="regime",
)
```

### The timeline against dated events

The upper panel is the same swim-lane band as above, restricted to the two-cluster model,
with seven dated episodes marked. The lower panel is the cumulative product of the equity
index returns on a logarithmic scale, drawn in the colour of the regime each month was
assigned to. A logarithmic scale is what makes the 1930s comparable to the 2010s: equal
vertical distances are equal percentage moves, so a fall of a given percentage takes up the
same space on the axis wherever in the century it happens.

```python
HISTORICAL_EVENTS = {
    1929: ("Great Crash", "top"),
    1937: ("Recession", "bottom"),
    1973: ("Oil crisis", "top"),
    1987: ("Black Monday", "bottom"),
    2000: ("Dot-com peak", "top"),
    2008: ("Global financial crisis", "bottom"),
    2020: ("COVID-19", "top"),
}
```

```python
def mark_events(ax: Axes, years: list[int]) -> None:
    """Draw a rule and a label at the first month of each dated episode."""
    for event_year, (event_label, side) in HISTORICAL_EVENTS.items():
        position = next((j for j, year in enumerate(years) if year == event_year), None)
        if position is None:
            continue
        ax.axvline(x=position, color=COLORS["neutral"], alpha=0.6, linewidth=1.0)
        ax.annotate(
            event_label,
            xy=(position, 1.7 if side == "top" else -0.7),
            ha="center",
            va="bottom" if side == "top" else "top",
            fontsize=7,
            color=COLORS["neutral"],
        )
```

```python
def plot_by_regime(ax: Axes, values: np.ndarray, in_risk_on: np.ndarray) -> None:
    """Draw one line whose every segment takes the colour of the month it ends in.

    The segment from month t-1 to month t is the move month t produced, so month t's regime
    is the one that colours it. Masking each regime into its own line instead would drop
    every segment that spans a transition, and would draw nothing at all for a one-month
    episode, which is most of them here.
    """
    y = np.asarray(values, dtype=float)
    points = np.column_stack([np.arange(len(y)), y]).reshape(-1, 1, 2)
    segments = np.concatenate([points[:-1], points[1:]], axis=1)
    drawable = np.isfinite(segments[:, :, 1]).all(axis=1)
    colours = np.where(in_risk_on[1:], COLORS["recede"], COLORS["blue"])
    ax.add_collection(LineCollection(segments[drawable], colors=colours[drawable], linewidths=1.2))
    finite = y[np.isfinite(y)]
    ax.set_xlim(0, len(y) - 1)
    ax.set_ylim(finite.min() * 0.9, finite.max() * 1.1)
```

```python
fig, (ax_band, ax_equity) = plt.subplots(
    2, 1, figsize=style.FIGSIZE["dual_v"], height_ratios=[1, 2], sharex=True
)

rows = np.where(labels_2 == risk_on_cluster, 1, 0)
draw_swimlanes(ax_band, rows, 2, lane_cmap)
ax_band.set_yticks([0, 1])
ax_band.set_yticklabels(["Risk-off", "Risk-on"])
ax_band.set_ylim(-0.5, 1.5)
mark_events(ax_band, factors_df.index.year.tolist())

cumulative = (1 + equity_returns).cumprod()
in_risk_on = labels_2 == risk_on_cluster
plot_by_regime(ax_equity, cumulative.to_numpy(), in_risk_on)
ax_equity.set_yscale("log")
ax_equity.set_ylabel("Growth of 1 unit invested (log scale)")
ax_equity.set_xlabel("Year")
ax_equity.set_xticks(ticks)
ax_equity.set_xticklabels(tick_labels)
ax_equity.grid(True, alpha=0.3)

style.add_message_title(
    ax_band,
    "Risk-off months come in short bursts, not in long bear markets",
    subtitle="Upper: assigned regime. Lower: cumulative equity index return, coloured by regime",
)
style.show_with_alt(
    fig,
    "A two-row regime band above a rising cumulative return curve on a log scale. The "
    "risk-off row is occupied in many short bursts spread across the whole century rather "
    "than in a few long blocks, and the curve is drawn dark through those same stretches.",
)
```

## Volatility inside each regime

A regime label is only worth carrying if it corresponds to something a risk manager would
act on, and the first thing to check is whether the two environments differ in how much
the market moves. The series below is the standard deviation of the equity index return
over a trailing twelve-month window, annualized by multiplying by the square root of
twelve, drawn in the colour of the regime assigned to each month. The dashed rules are the
average of that rolling series within each regime.

```python
rolling_vol = equity_returns.rolling(ROLLING_VOL_MONTHS).std() * np.sqrt(MONTHS_PER_YEAR)

rolling_vol_by_regime = (
    rolling_vol.groupby(regime_name)
    .agg(["mean", "median", "max"])
    .rename(columns={"mean": "Mean (%)", "median": "Median (%)", "max": "Max (%)"})
    .mul(100)
    .loc[["Risk-on", "Risk-off"]]
)
rolling_vol_by_regime.index.name = "Regime"
rolling_vol_by_regime.style.format("{:.1f}")
```

```python
fig, (ax_band, ax_vol) = plt.subplots(
    2, 1, figsize=style.FIGSIZE["dual_v"], height_ratios=[1, 3], sharex=True
)

draw_swimlanes(ax_band, rows, 2, lane_cmap)
ax_band.set_yticks([0, 1])
ax_band.set_yticklabels(["Risk-off", "Risk-on"])
ax_band.set_ylim(-0.5, 1.5)

plot_by_regime(ax_vol, rolling_vol.to_numpy() * 100, in_risk_on)
for regime, colour in (("Risk-on", COLORS["copper"]), ("Risk-off", COLORS["amber"])):
    ax_vol.axhline(
        rolling_vol_by_regime.loc[regime, "Mean (%)"],
        color=colour,
        linestyle="--",
        linewidth=1.2,
        label=f"{regime} average",
    )
ax_vol.set_ylabel(f"Trailing {ROLLING_VOL_MONTHS}-month volatility, annualized (%)")
ax_vol.set_xlabel("Year")
ax_vol.set_xticks(ticks)
ax_vol.set_xticklabels(tick_labels)
ax_vol.legend(loc="upper right")
ax_vol.grid(True, alpha=0.3)

style.add_message_title(
    ax_band,
    "The risk-off cluster sits on the higher-volatility stretches of the century",
    subtitle=f"Trailing {ROLLING_VOL_MONTHS}-month standard deviation of the equity index, "
    f"annualized by the square root of {MONTHS_PER_YEAR}",
)
style.show_with_alt(
    fig,
    "A two-row regime band above a rolling volatility series in percent. The dark risk-off "
    "segments concentrate on the peaks of the volatility series, and the risk-off average "
    "rule sits above the risk-on one.",
)
```

### Figure 1.5 inputs

The print version of Figure 1.5 is drawn by
`book/01_process_is_edge/figures/scripts/generate_figure_1_5_factor_regimes_volatility.py`
in the book repository. It reads the arrays written below, so the book build renders the
figure from this notebook's fit rather than re-estimating the mixture.

```python
ARTIFACT_DIR = OUTPUT_DIR / "figure_1_5"
ARTIFACT_DIR.mkdir(parents=True, exist_ok=True)
np.savez(
    ARTIFACT_DIR / "inputs.npz",
    dates=factors_df.index.astype("datetime64[ns]").astype("int64"),
    labels_2=np.asarray(labels_2, dtype=np.int64),
    good_regime=np.int64(risk_on_cluster),
    rolling_vol=rolling_vol.to_numpy(dtype=float),
    risk_on_mean_vol=float(rolling_vol_by_regime.loc["Risk-on", "Mean (%)"]) / 100,
    risk_off_mean_vol=float(rolling_vol_by_regime.loc["Risk-off", "Mean (%)"]) / 100,
    start_year=np.int64(start_year),
    end_year=np.int64(end_year),
)
```

## What each factor earned in each regime

The volatility check says the two environments differ in how much the market moved. The
more useful question for portfolio construction is whether the strategies behave
differently in them, because a strategy that earns its premium in both is a diversifier
and one that only earns it in the calm environment is leverage on the market in disguise.

Each bar below is the mean monthly return of one series over the months assigned to one
regime, multiplied by twelve. This is an average within a set of months, not the return of
a strategy that switched between them - no rule available in real time would have produced
this partition.

```python
annualized_by_regime = (
    factors_df.groupby(regime_name).mean().mul(MONTHS_PER_YEAR * 100).rename(columns=SHORT_NAMES)
).loc[["Risk-on", "Risk-off"]]
```

```python
fig, ax = plt.subplots(figsize=style.FIGSIZE["single_tall"])

x = np.arange(annualized_by_regime.shape[1])
width = 0.38
ax.bar(
    x - width / 2,
    annualized_by_regime.loc["Risk-on"],
    width,
    label="Risk-on",
    color=COLORS["recede"],
)
ax.bar(
    x + width / 2,
    annualized_by_regime.loc["Risk-off"],
    width,
    label="Risk-off",
    color=COLORS["blue"],
)
ax.set_xticks(x)
ax.set_xticklabels(annualized_by_regime.columns, rotation=30, ha="right")
ax.set_ylabel("Mean return within the regime, annualized (%)")
style.zero_line(ax)
ax.legend(loc="upper right")
ax.grid(axis="y", alpha=0.3)

style.add_message_title(
    ax,
    "Value and bonds hold up in the risk-off months; carry and defensive do not",
    subtitle="Mean monthly return within each regime times twelve, by series",
)
style.show_with_alt(
    fig,
    "Paired bars per factor. Value and bonds have taller risk-off bars than risk-on bars; "
    "momentum, equity, carry and defensive have taller risk-on bars, with carry and "
    "defensive falling below zero in the risk-off months.",
)
```

```python
annualized_by_regime.T.style.format("{:+.1f}")
```

The two strategies a calm market reads as safe income - carry, which collects a yield
spread, and defensive, which is long low-risk assets - are the ones that stop paying when
conditions deteriorate. Value and bonds go the other way and earn more.

A strategy's average return over a full sample therefore says nothing about whether it
will be there in the environment a portfolio needs it in. Finding out means conditioning
on the environment and looking, which is what makes a regime label worth estimating even
when nothing can act on it.

## Regime statistics

The table below characterises each regime through the equity index alone. Every statistic
is computed on one stream: the monthly returns of the months the model assigned to that
regime, concatenated in date order with the months in the other regime removed.

For the mean, the standard deviation and the Sharpe ratio that is an ordinary conditional
statistic - the average of a set of months does not depend on their order. **Maximum
drawdown** is different, because it is a path statistic: it compounds the stream, tracks
its running high-water mark, and reports the largest percentage fall from that mark to a
later point. Compounding a stream whose gaps have been closed up is not the same as
compounding the index, so the number below describes an investor who held equities during
one regime and held nothing during the other. The section after the table shows how far
apart the two readings are.

```python
def drawdown_path(returns: pd.Series) -> pd.Series:
    """Fractional decline from the running high-water mark of the compounded stream."""
    wealth = (1 + returns).cumprod()
    return wealth / wealth.cummax() - 1
```

```python
regime_statistics = pd.DataFrame(
    [
        {
            "Regime": regime,
            "Months": int(len(returns)),
            "Share of sample": len(returns) / len(equity_returns),
            "Return, annualized": returns.mean() * MONTHS_PER_YEAR,
            "Volatility, annualized": returns.std() * np.sqrt(MONTHS_PER_YEAR),
            "Sharpe ratio": sharpe_ratio(returns, periods_per_year=MONTHS_PER_YEAR),
            "Max drawdown": float(drawdown_path(returns).min()),
        }
        for regime, returns in (
            (name, equity_returns[regime_name == name]) for name in ("Risk-on", "Risk-off")
        )
    ]
).set_index("Regime")
regime_statistics.style.format(
    {
        "Share of sample": "{:.1%}",
        "Return, annualized": "{:+.1%}",
        "Volatility, annualized": "{:.1%}",
        "Sharpe ratio": "{:.2f}",
        "Max drawdown": "{:.1%}",
    }
)
```

```python
pooled_vol_ratio = (
    regime_statistics.loc["Risk-off", "Volatility, annualized"]
    / regime_statistics.loc["Risk-on", "Volatility, annualized"]
)
rolling_vol_ratio = (
    rolling_vol_by_regime.loc["Risk-off", "Mean (%)"]
    / rolling_vol_by_regime.loc["Risk-on", "Mean (%)"]
)
display(
    Markdown(
        f"Risk-off volatility is **{pooled_vol_ratio:.2f}x** risk-on volatility measured on "
        f"the pooled monthly returns, against **{rolling_vol_ratio:.2f}x** measured as the "
        "average of the trailing twelve-month series."
    )
)
```

The two ratios answer different questions and a reader should expect them to differ. The
rolling series smooths each month against the eleven around it, so a single extreme month
raises twelve windows a little; the pooled standard deviation lets it enter once, at full
size. The pooled figure is the one that describes what a portfolio holding through the
regime experiences, and the rolling figure is the one a monitoring dashboard would show.
Reporting either alone would be defensible; reporting one and calling it "the" volatility
ratio would not.

### What the restricted drawdown actually measures

The drawdown reported for the risk-off regime is large, and it is worth reading where it
came from before quoting it. The cell below finds the high-water mark of the restricted
stream and the month it fell furthest below it, and puts the drawdown of the index itself
over the same span beside it.

```python
risk_off_returns = equity_returns[regime_name == "Risk-off"]
risk_off_drawdown = drawdown_path(risk_off_returns)
trough_month = risk_off_drawdown.idxmin()
peak_month = (1 + risk_off_returns).cumprod().loc[:trough_month].idxmax()

restricted_drawdown = float(risk_off_drawdown.min())
index_over_span = equity_returns.loc[peak_month:trough_month]
index_drawdown = float(drawdown_path(index_over_span).min())
display(
    Markdown(
        f"The restricted stream peaks in **{peak_month:%B %Y}** and reaches its low in "
        f"**{trough_month:%B %Y}**, a fall of "
        f"**{abs(restricted_drawdown):.1%}** spread over "
        f"**{trough_month.year - peak_month.year} years**. Over the same span the equity "
        f"index itself fell at most **{abs(index_drawdown):.1%}** from its own high-water "
        "mark."
    )
)
```

Nobody lived through the first of those two numbers. The risk-off months are selected by
the mixture's assignment, not by the sign of the equity return, so the stream holds gains
as well as losses - its average monthly return is positive. What it does not hold is the
months in between, and closing those gaps puts a loss from one decade immediately beside a
loss from the next, with none of the recovery that separated them. Compounded, decades of
separate episodes read as one uninterrupted decline.

The statistic is not wrong; it is the correct drawdown of the object it was computed on.
But the object is a construction of this notebook, and a reader who takes it for a crash
has been misled by a table that never said which object it meant.

That is the general form: a conditional statistic inherits the assumptions of the
conditioning. Averages survive it, because a mean does not care in what order its terms
arrive. Anything that reads a path - a drawdown, a run length, a time to recovery, a
stop-loss - is a different statistic once the path has been cut up, and reporting it
without saying so is how a number that means one thing gets read as another.

## How long the regimes last

A regime label is useful to a strategy only if it persists long enough to act on. An
**episode** below is a maximal run of consecutive months carrying the same label, and the
table reports how many episodes each regime had and how long they ran.

```python
episode_id = (labels_2 != np.roll(labels_2, 1)).cumsum()
episode_lengths = (
    pd.DataFrame({"regime": regime_name.to_numpy(), "episode": episode_id})
    .groupby(["regime", "episode"])
    .size()
    .rename("months")
    .reset_index()
)

duration = (
    episode_lengths.groupby("regime")["months"]
    .agg(["mean", "max", "count"])
    .rename(
        columns={"mean": "Mean length (months)", "max": "Longest (months)", "count": "Episodes"}
    )
    .loc[["Risk-on", "Risk-off"]]
)
duration.index.name = "Regime"
duration.style.format({"Mean length (months)": "{:.1f}"})
```

```python
switches = int(np.sum(np.diff(labels_2) != 0))
months_between_switches = len(labels_2) / switches
display(
    Markdown(
        f"The two-cluster model changes regime **{switches}** times over "
        f"{end_year - start_year} years, about once every "
        f"**{months_between_switches:.1f} months**."
    )
)
```

That is the finding that decides what this model is for. A label that changes several
times a year cannot drive a reallocation: the turnover would be paid on every switch,
including the ones the model reverses a month later, and the transaction costs would
arrive whether or not the regime call was right. The same label is useful as a lens on
history and as one input to a risk report, which is the use Section 1.4 argues for.

A model built to be acted on would be built differently. It would carry a transition
structure that makes staying in a regime more likely than leaving it - a hidden Markov
model is the standard choice - and it would be fitted forward in time so that each
month's label used only the months before it. `09_model_based_features/11_hmm_regimes`
builds the first, and `09_model_based_features/13_regime_as_feature` turns the result into
a feature a model downstream can read.

## Key takeaways

- **The number of clusters is a decision, not an estimate.** BIC, AIC and the silhouette
  score reward different things, so they will disagree; pick the one whose objective
  matches what the labels are for, and say which you picked and why.
- **Cluster numbers carry no meaning until something outside the model names them.** Here
  the average equity return within each cluster does the naming, and a different seed can
  renumber the clusters without changing the partition.
- **Condition on the environment before trusting a factor's average.** A strategy's
  full-sample premium can be entirely earned in one environment, which is invisible in the
  pooled number and decisive for whether it diversifies anything.
- **State how a conditional statistic was built.** Restricting a return stream to the
  months in a regime and concatenating them produces a drawdown no investor experienced
  unless they held only during that regime; the same number is either informative or
  misleading depending on whether that construction is stated.
- **Fitting on the whole sample buys description, not prediction.** Every label here is
  informed by the months after it, so the labels cannot be used as features. Walk-forward
  fitting is what makes them usable, and Chapter 6 onward introduces it.

**Known limitations.** The mixture has no notion of time: it can assign January 1932 and
January 2019 to the same cluster and has no mechanism preferring a month to keep the
previous month's label, which is why the episodes are short. The AQR series are monthly, so
nothing here resolves anything faster than a month. And the panel is a century long, over
which the composition of the asset classes and the microstructure of the markets changed
considerably; the model treats 1927 and 2024 as draws from the same set of distributions.

**Next**: `macro_regimes` asks the same question of macroeconomic indicators rather than
factor returns, and checks the resulting clusters against realized equity volatility.
![notebook output](figures/p1_1.png)
![notebook output](figures/p1_2.png)
![notebook output](figures/p1_3.png)
![notebook output](figures/p1_4.png)
![notebook output](figures/p1_5.png)