Árboles de expansión mínima de correlación para diversificar carteras

# Financial Networks for Portfolio Construction


# Financial Networks for Portfolio Construction

**Chapter 23: Knowledge Graphs for Financial AI**

**Docker image**: `ml4t`

This notebook demonstrates using **network theory** for portfolio construction
and risk diagnostics. It forms the investable universe and estimates the
correlation network before a separate evaluation window.

**Learning Objectives**:
- Build correlation networks and extract minimum spanning trees from real price data
- Compute centrality metrics that describe MST structure
- Construct network-diversified portfolios that underweight central nodes
- Run a stylized, order-independent shock-propagation sensitivity

**Book Reference**: Chapter 23, Section 23.5 (Correlations and portfolios in financial networks)

**Prerequisites**: Familiarity with correlation matrices and portfolio construction.
Requires the US equities dataset (3,199 current and delisted securities).

Based on: Konstantinov et al. (2023) - Financial Networks

```python
"""Financial networks for portfolio construction and shock-sensitivity analysis."""

from __future__ import annotations

from logging import getLogger

import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
import polars as pl
from IPython.display import Markdown, display
from matplotlib.colors import LinearSegmentedColormap
from scipy.sparse.csgraph import minimum_spanning_tree

from data import load_us_equities
from utils.reproducibility import set_global_seeds
from utils.style import (
    COLORS,
    FIGSIZE,
    add_message_title,
    format_pct_axis,
    ml4t_palette,
    show_with_alt,
)
```

```python
# Production defaults - Papermill overrides for testing
N_ASSETS = 100
ESTIMATION_DAYS = 504
EVALUATION_DAYS = 252
# Added to degree centrality before inverting, so a zero-degree node gets a finite
# score. Section 6 measures what it does to the weights.
CENTRALITY_FLOOR = 0.01
MAX_WEIGHT = 0.10
# A date is a market date once this many symbols report on it. Section 2 uses it to
# drop the stray provider rows dated to holidays; a subsampled panel needs a lower one.
MIN_SYMBOLS_PER_DATE = 1_000
SEED = 42
```

```python
set_global_seeds(SEED)
getLogger("matplotlib.font_manager").setLevel("ERROR")
```

## 1. Configuration

```python
# Configuration
MIN_DATE = "2015-01-01"

# A cap below 1/N is infeasible: no weight vector on N assets can sum to one while
# every position stays under it, and equal weight is the first thing to violate it.
if MAX_WEIGHT * N_ASSETS < 1:
    raise ValueError(f"MAX_WEIGHT={MAX_WEIGHT} cannot hold {N_ASSETS} assets summing to one")
print(f"Assets: {N_ASSETS}")
print(f"Estimation days: {ESTIMATION_DAYS}")
print(f"Evaluation days: {EVALUATION_DAYS}")
print(f"Min date: {MIN_DATE}")
```

## 2. Load Real Market Data

Load US equities, define a broad-market calendar, and form the universe from
the estimation window only.

```python
# The canonical loader returns `symbol` and `timestamp`.
df = load_us_equities(start_date=MIN_DATE)
df = df.with_columns((pl.col("adj_close") * pl.col("adj_volume")).alias("dollar_volume"))

print(f"Loaded {len(df):,} rows")
print(f"Date range: {df['timestamp'].min()} to {df['timestamp'].max()}")
print(f"Unique symbols: {df['symbol'].n_unique()}")
```

```python
# A handful of provider rows fall on non-market dates. Requiring broad coverage
# removes those rows without consulting future returns for universe membership.
market_dates = (
    df.group_by("timestamp")
    .agg(pl.col("symbol").n_unique().alias("symbols_observed"))
    .filter(pl.col("symbols_observed") >= MIN_SYMBOLS_PER_DATE)
    .sort("timestamp")["timestamp"]
)
required_dates = ESTIMATION_DAYS + EVALUATION_DAYS
if len(market_dates) < required_dates + 1:
    raise ValueError(
        f"Need {required_dates + 1} dates carrying {MIN_SYMBOLS_PER_DATE} symbols or more, "
        f"found {len(market_dates)}"
    )

return_dates = market_dates.tail(required_dates + 1)
analysis_dates = return_dates.tail(required_dates)
estimation_dates = analysis_dates.head(ESTIMATION_DAYS)
evaluation_dates = analysis_dates.tail(EVALUATION_DAYS)
estimation_start, estimation_end = estimation_dates.min(), estimation_dates.max()
evaluation_start, evaluation_end = evaluation_dates.min(), evaluation_dates.max()

print(f"Estimation: {estimation_start} to {estimation_end}")
print(f"Evaluation: {evaluation_start} to {evaluation_end}")
```

```python
# Rank only observations available by the estimation boundary. Requiring a
# complete estimation history is an ex-ante eligibility rule.
avg_dollar_volume = (
    df.filter(pl.col("timestamp").is_in(estimation_dates.implode()))
    .group_by("symbol")
    .agg(pl.col("dollar_volume").mean().alias("avg_dv"), pl.len().alias("n_obs"))
    .filter(pl.col("n_obs") == ESTIMATION_DAYS)
    .sort(["avg_dv", "symbol"], descending=[True, False])
    .head(N_ASSETS)
)
if avg_dollar_volume.height != N_ASSETS:
    raise ValueError(f"Need {N_ASSETS} eligible symbols, found {avg_dollar_volume.height}")

selected_symbols = avg_dollar_volume["symbol"].to_list()
df_filtered = df.filter(
    pl.col("symbol").is_in(selected_symbols) & pl.col("timestamp").is_in(return_dates.implode())
)

print(f"\nSelected {len(selected_symbols)} symbols using estimation-only dollar volume")
print(f"Sample symbols: {selected_symbols[:10]}")
```

## 3. Compute Returns and Correlation Network

Pairwise return correlations are mapped to a distance metric
$$d_{ij}=\sqrt{2\,(1-\rho_{ij})}$$
which lies in $[0, 2]$ and satisfies the triangle inequality, so the minimum
spanning tree built from these distances is well defined.

```python
# Compute daily returns before splitting the windows.
df_filtered = df_filtered.sort(["symbol", "timestamp"])
df_filtered = df_filtered.with_columns(
    (pl.col("adj_close") / pl.col("adj_close").shift(1).over("symbol") - 1).alias("returns")
)

# Pivot to a fixed, estimation-formed column order.
returns_wide = (
    df_filtered.pivot(
        on="symbol",
        index="timestamp",
        values="returns",
    )
    .select(["timestamp", *selected_symbols])
    .sort("timestamp")
)
```

The fixed universe can contain later symbol exits or isolated missing quotes, and
a zero return holds the marked value constant for that date. The diagnostic below
makes that cash and stale-mark assumption explicit, counted over the analysis
dates alone. A return is a difference between two closes, so the first row of the
whole frame is null for every symbol by construction: counting the frame would
report one null per symbol that no exit or missing quote produced, on rows the
analysis never reads.

```python
analysis_frame = returns_wide.filter(pl.col("timestamp").is_in(analysis_dates.implode()))
missing_by_window = analysis_frame.select(pl.exclude("timestamp").is_null().sum()).sum_horizontal()[
    0
]
evaluation_missing = (
    returns_wide.filter(pl.col("timestamp").is_in(evaluation_dates.implode()))
    .select(pl.exclude("timestamp").is_null().sum())
    .sum_horizontal()[0]
)
returns_wide = returns_wide.fill_null(0.0)

estimation_frame = returns_wide.filter(pl.col("timestamp").is_in(estimation_dates.implode()))
evaluation_frame = returns_wide.filter(pl.col("timestamp").is_in(evaluation_dates.implode()))
estimation_returns = estimation_frame.drop("timestamp").to_numpy()
evaluation_returns = evaluation_frame.drop("timestamp").to_numpy()
evaluation_index = evaluation_frame["timestamp"].to_list()

asset_names = selected_symbols
n_assets = len(asset_names)

print(f"Estimation returns: {estimation_returns.shape}")
print(f"Evaluation returns: {evaluation_returns.shape}")
print(
    f"Missing symbol-date returns set to 0 within the analysis window: {missing_by_window} "
    f"of {analysis_frame.height * N_ASSETS:,}"
)
print(f"Missing returns within evaluation: {evaluation_missing}")
print(f"Sample symbols: {asset_names[:5]}")
```

```python
assert estimation_returns.shape == (ESTIMATION_DAYS, N_ASSETS)
assert evaluation_returns.shape == (EVALUATION_DAYS, N_ASSETS)
assert np.isfinite(estimation_returns).all()
assert np.isfinite(evaluation_returns).all()
assert estimation_end < evaluation_start
```

```python
# Calculate the correlation matrix on estimation returns only.
corr_matrix = np.corrcoef(estimation_returns.T)

# Map correlations to MST distances (formula in the markdown above).
distance_matrix = np.sqrt(2 * (1 - np.clip(corr_matrix, -1, 1)))
np.fill_diagonal(distance_matrix, 0)

print(f"Correlation matrix shape: {corr_matrix.shape}")
print(f"Correlation range: [{corr_matrix.min():.3f}, {corr_matrix.max():.3f}]")
print(f"Distance range: [{distance_matrix.min():.3f}, {distance_matrix.max():.3f}]")

# Analyze correlation structure
corr_values = corr_matrix[np.triu_indices(n_assets, k=1)]
print(f"\nAverage correlation: {corr_values.mean():.3f}")
print(f"Median correlation: {np.median(corr_values):.3f}")
```

## 4. Minimum Spanning Tree

Extract hierarchical structure from the correlation network.

```python
def compute_mst(distance_matrix: np.ndarray) -> tuple[np.ndarray, list]:
    """Compute minimum spanning tree from distance matrix."""
    # scipy.sparse.csgraph.minimum_spanning_tree
    mst = minimum_spanning_tree(distance_matrix).toarray()

    # Make symmetric
    mst = mst + mst.T

    # Extract edges
    edges = []
    for i in range(len(mst)):
        for j in range(i + 1, len(mst)):
            if mst[i, j] > 0:
                edges.append((i, j, mst[i, j]))

    return mst, edges
```

```python
mst_matrix, mst_edges = compute_mst(distance_matrix)

print(f"MST edges: {len(mst_edges)}")
print(f"Expected edges (N-1): {n_assets - 1}")

# Verify MST properties
total_distance = sum(e[2] for e in mst_edges)
print(f"Total MST distance: {total_distance:.3f}")

# Show some MST connections (most correlated pairs)
sorted_edges = sorted(mst_edges, key=lambda e: e[2])
print("\nMost connected pairs (lowest distance = highest correlation):")
for i, j, d in sorted_edges[:5]:
    rho = 1 - d**2 / 2  # Convert back to correlation
    print(f"  {asset_names[i]} - {asset_names[j]}: ρ={rho:.3f}")
```

## 5. Network Centrality Metrics

Calculate centrality measures that describe estimated MST structure.

### Degree Centrality

Counts the fraction of nodes each asset is connected to. High-degree nodes
are structural hubs in this estimated correlation tree.

```python
def compute_degree_centrality(adjacency: np.ndarray) -> np.ndarray:
    """Compute degree centrality (normalized)."""
    binary_adj = (adjacency > 0).astype(float)
    degrees = binary_adj.sum(axis=1)
    return degrees / (len(degrees) - 1)
```

### Betweenness Centrality

The fraction of MST shortest paths that pass through each node. Bridge nodes
connecting otherwise separate clusters score highest. Computed on the MST with
edge distances via NetworkX.

```python
def build_mst_graph(mst_matrix: np.ndarray, n: int) -> nx.Graph:
    """Build a weighted NetworkX graph from the symmetric MST distance matrix."""
    graph = nx.Graph()
    graph.add_nodes_from(range(n))
    for i in range(len(mst_matrix)):
        for j in range(i + 1, len(mst_matrix)):
            if mst_matrix[i, j] > 0:
                graph.add_edge(i, j, weight=mst_matrix[i, j])
    return graph
```

### Closeness Centrality

Inverse of the average MST shortest-path distance to all other nodes. Assets
near the tree center have short correlation-distance paths to other assets.

```python
# Degree centrality on the MST adjacency; betweenness and closeness via NetworkX
# shortest paths using the MST edge distances.
mst_graph = build_mst_graph(mst_matrix, n_assets)
degree_centrality = compute_degree_centrality(mst_matrix)
betweenness_dict = nx.betweenness_centrality(mst_graph, weight="weight")
closeness_dict = nx.closeness_centrality(mst_graph, distance="weight")
betweenness_centrality = np.array([betweenness_dict[i] for i in range(n_assets)])
closeness_centrality = np.array([closeness_dict[i] for i in range(n_assets)])

# Create centrality DataFrame
centrality_df = pl.DataFrame(
    {
        "symbol": asset_names,
        "degree": degree_centrality,
        "betweenness": betweenness_centrality,
        "closeness": closeness_centrality,
    }
)

print("TOP 10 ASSETS BY CENTRALITY (MST)")

# Sort by composite score
centrality_df = centrality_df.with_columns(
    ((pl.col("degree") + pl.col("betweenness") + pl.col("closeness")) / 3).alias("composite")
)

top_central = centrality_df.sort("composite", descending=True).head(10)
print(top_central)
```

**Interpretation**: High-centrality assets are hubs in the estimated MST.
Centrality describes correlation-network structure; it does not establish
causal influence or systemic importance.

## 6. Network-Based Portfolio Construction

Use centrality metrics to construct diversified portfolios.

### Network-Diversified Weights

Inverse centrality weighting: underweight estimated MST hubs.

```python
def project_capped_weights(raw_weights: np.ndarray, max_weight: float) -> np.ndarray:
    """Project nonnegative scores to sum to one without violating a weight cap."""
    if max_weight <= 0 or max_weight * len(raw_weights) < 1:
        raise ValueError("max_weight is infeasible for this portfolio size")
    scores = np.asarray(raw_weights, dtype=float)
    if (scores < 0).any() or not np.isfinite(scores).all() or scores.sum() <= 0:
        raise ValueError("raw_weights must be finite, nonnegative, and nonzero")

    weights = np.zeros_like(scores)
    free = np.ones(len(scores), dtype=bool)
    remaining = 1.0
    while free.any():
        free_scores = scores[free]
        if free_scores.sum() > 0:
            proposal = free_scores / free_scores.sum() * remaining
        else:
            proposal = np.full(free.sum(), remaining / free.sum())
        capped = proposal > max_weight
        if not capped.any():
            weights[free] = proposal
            break
        free_idx = np.flatnonzero(free)
        weights[free_idx[capped]] = max_weight
        free[free_idx[capped]] = False
        remaining = 1.0 - weights.sum()
    return weights
```

The portfolio rule converts inverse centrality scores into feasible weights
with the capped-simplex projection above.

```python
def network_diversified_weights(
    centrality: np.ndarray,
    max_weight: float = MAX_WEIGHT,
    floor: float = CENTRALITY_FLOOR,
) -> np.ndarray:
    """
    Construct portfolio weights using inverse centrality.

    High-centrality assets are estimated MST hubs, so the rule
    underweights them by construction.
    """
    inverse_cent = 1 / (centrality + floor)
    return project_capped_weights(inverse_cent, max_weight)
```

### Benchmark Portfolios

Equal weight and inverse-volatility as baselines for comparison.

```python
def equal_weight_portfolio(n: int) -> np.ndarray:
    """Simple equal weight portfolio."""
    return np.ones(n) / n
```

### Inverse-Volatility Weights

Weight each asset by the inverse of its volatility - a simple risk-based
benchmark. (A constrained minimum-variance optimizer would also use the full
covariance structure; this notebook uses the lighter inverse-volatility rule.)

```python
def inverse_volatility_weights(returns: np.ndarray, max_weight: float = MAX_WEIGHT) -> np.ndarray:
    """Weight assets by inverse estimation-window volatility."""
    cov = np.cov(returns.T)
    vols = np.sqrt(np.diag(cov))
    inv_vol = 1 / (vols + 1e-8)
    return project_capped_weights(inv_vol, max_weight)
```

```python
# Construct different portfolios
w_equal = equal_weight_portfolio(n_assets)
w_inv_vol = inverse_volatility_weights(estimation_returns)
w_network_div = network_diversified_weights(degree_centrality)

# Portfolio statistics
portfolios = {
    "Equal Weight": w_equal,
    "Inverse Volatility": w_inv_vol,
    "Network Diversified": w_network_div,
}

print("PORTFOLIO CONCENTRATION ANALYSIS")
```

```python
concentration_rows = []
for name, weights in portfolios.items():
    hhi = np.sum(weights**2)
    effective_n = 1 / hhi
    max_weight = weights.max()
    concentration_rows.append(
        {
            "portfolio": name,
            "hhi": hhi,
            "effective_names": effective_n,
            "max_weight": max_weight,
        }
    )
    assert np.isclose(weights.sum(), 1.0)
    assert weights.min() >= 0
    assert weights.max() <= MAX_WEIGHT + 1e-12

concentration = pl.DataFrame(concentration_rows)
print(concentration)
```

### Concentration Finding

Equal weight is not a comparison here, it is the ceiling. The effective-name
count is the reciprocal of the Herfindahl index, which any weight vector on N
assets maximises by putting 1/N in each, so equal weight scores exactly N and no
rule can beat it. Its largest position is 1/N, which for this universe never
approaches the cap either. Reporting a network portfolio "versus equal weight"
on this measure can only report how far short it falls.

The number that carries information is the shortfall, and the rule to compare
against is one that also trades off something: inverse volatility, which is
already computed above.

```python
equal_concentration = concentration.filter(pl.col("portfolio") == "Equal Weight").row(0, named=True)
concentration = concentration.with_columns(
    (pl.col("effective_names") / float(n_assets)).alias("share_of_ceiling")
)
if not np.isclose(equal_concentration["effective_names"], n_assets):
    raise RuntimeError(
        f"equal weight scored {equal_concentration['effective_names']:.4f} effective names "
        f"rather than the {n_assets} its construction guarantees"
    )
print(concentration)

network_concentration = concentration.filter(pl.col("portfolio") == "Network Diversified").row(
    0, named=True
)
inv_vol_concentration = concentration.filter(pl.col("portfolio") == "Inverse Volatility").row(
    0, named=True
)
display(
    Markdown(
        f"**Finding:** equal weight reaches the ceiling of {n_assets} effective names "
        f"by construction. Network diversification holds "
        f"{network_concentration['effective_names']:.1f}, "
        f"{network_concentration['share_of_ceiling']:.0%} of it, against "
        f"{inv_vol_concentration['effective_names']:.1f} "
        f"({inv_vol_concentration['share_of_ceiling']:.0%}) for inverse volatility. "
        f"The largest network position is {network_concentration['max_weight']:.1%} "
        f"against a {MAX_WEIGHT:.0%} cap."
    )
)
```

The concentration table is determined entirely by weights fixed at the
estimation boundary. What it shows is how much dispersion each rule gives up in
exchange for the thing it targets, not whether it disperses capital better than
a rule that does nothing else.

### What the Centrality Floor Does

`CENTRALITY_FLOOR` is added to every centrality before inverting, so it decides
how much the rule can separate the least connected assets from each other. A leaf
in this tree has centrality 1/(N-1), which is the same order as the floor, so the
floor roughly halves a leaf's raw score before the projection. Removing it
entirely would divide by zero for any isolated node; shrinking it sharpens the
rule until the cap absorbs the difference.

```python
floor_rows = []
for candidate_floor in (CENTRALITY_FLOOR / 10, CENTRALITY_FLOOR, CENTRALITY_FLOOR * 10):
    candidate_weights = network_diversified_weights(degree_centrality, floor=candidate_floor)
    floor_rows.append(
        {
            "floor": candidate_floor,
            "effective_names": 1 / np.sum(candidate_weights**2),
            "max_weight": candidate_weights.max(),
            "weight_ratio_leaf_to_hub": float(
                candidate_weights.max() / candidate_weights[int(np.argmax(degree_centrality))]
            ),
        }
    )
floor_sensitivity = pl.DataFrame(floor_rows)
print("Sensitivity of the weights to CENTRALITY_FLOOR:")
print(floor_sensitivity)
```

## 7. Stylized Shock-Propagation Sensitivity

This deterministic exercise diffuses a shock across correlations above a
threshold. It is not a causal contagion model or an estimate of portfolio loss.

```python
def simulate_shock_diffusion(
    corr_matrix: np.ndarray,
    shock_asset: int,
    shock_magnitude: float = -0.10,
    contagion_threshold: float = 0.5,
    attenuation: float = 0.5,
    max_rounds: int = 5,
) -> dict:
    """Diffuse a shock synchronously through a thresholded correlation graph."""
    n = len(corr_matrix)
    impacted = np.zeros(n, dtype=bool)
    impacts = np.zeros(n)
    impacted[shock_asset] = True
    impacts[shock_asset] = shock_magnitude
    frontier = np.array([shock_asset])
    rounds = 0

    while len(frontier) and rounds < max_rounds:
        proposals = np.zeros(n)
        for source in frontier:
            eligible = (~impacted) & (corr_matrix[source] > contagion_threshold)
            propagated = impacts[source] * corr_matrix[source] * attenuation
            proposals[eligible] = np.minimum(proposals[eligible], propagated[eligible])
        frontier = np.flatnonzero(proposals < 0)
        impacts[frontier] = proposals[frontier]
        impacted[frontier] = True
        rounds += 1

    return {
        "shocked_asset": shock_asset,
        "n_impacted": int(impacted.sum()),
        "impacts": impacts,
        "propagation_rounds": rounds,
    }
```

The scenarios rank assets by MST degree and the diffusion runs over correlations
above a threshold, and the MST is built from those same correlations. So "the
high-degree node reaches more assets" restates how the tree was constructed
unless the threshold binds somewhere the tree does not. The count below says how
many pairs clear it at all.

```python
CONTAGION_THRESHOLD = 0.5
pairs_above_threshold = int((corr_values > CONTAGION_THRESHOLD).sum())
print(
    f"Pairs with correlation above {CONTAGION_THRESHOLD}: {pairs_above_threshold:,} of "
    f"{len(corr_values):,} ({pairs_above_threshold / len(corr_values):.1%})"
)
```

```python
# Compare shocks from high- and low-degree MST nodes.
print("STYLIZED SHOCK-DIFFUSION SENSITIVITY")

# Find highest and lowest centrality assets. argmin runs over positive-degree
# nodes only, then maps back to the original asset index.
high_central_idx = int(np.argmax(degree_centrality))
positive_degree_idx = np.where(degree_centrality > 0)[0]
low_central_idx = int(positive_degree_idx[np.argmin(degree_centrality[positive_degree_idx])])

# Every leaf of a spanning tree shares the minimum degree, so the low-centrality
# pick is one of many tied assets rather than a distinguished one. Say how many.
min_centrality = degree_centrality[low_central_idx]
tied_at_minimum = int(np.isclose(degree_centrality, min_centrality).sum())
print(
    f"Least connected asset: {asset_names[low_central_idx]}, one of {tied_at_minimum} "
    f"tied at centrality {min_centrality:.4f}"
)

scenarios = [
    ("High centrality shock", high_central_idx),
    ("Low centrality shock", low_central_idx),
]

shock_rows = []
for name, idx in scenarios:
    result = simulate_shock_diffusion(corr_matrix, idx, contagion_threshold=CONTAGION_THRESHOLD)
    shock_rows.append(
        {
            "scenario": name,
            "symbol": asset_names[idx],
            "assets_reached": result["n_impacted"],
            "rounds": result["propagation_rounds"],
            "equal_weight_sensitivity": float(result["impacts"] @ w_equal),
            "network_weight_sensitivity": float(result["impacts"] @ w_network_div),
        }
    )

shock_summary = pl.DataFrame(shock_rows)
print(shock_summary)
```

### Shock-Sensitivity Finding

The comparison below is tied to the executed scenarios and remains explicitly
bounded to this stylized diffusion rule.

```python
high_shock, low_shock = shock_rows
same_reach = high_shock["assets_reached"] == low_shock["assets_reached"]
display(
    Markdown(
        f"**Finding:** the high-degree {high_shock['symbol']} shock and the "
        f"low-degree {low_shock['symbol']} shock reach "
        + (
            f"the same {high_shock['assets_reached']} assets"
            if same_reach
            else f"{high_shock['assets_reached']} and {low_shock['assets_reached']} assets"
        )
        + f". Their sensitivities still differ, "
        f"{high_shock['equal_weight_sensitivity']:.1%} against "
        f"{low_shock['equal_weight_sensitivity']:.1%} under equal weights, so the "
        f"difference is in how large an impact each asset receives rather than in "
        f"how many receive one. Holding the scenario fixed and changing the weight "
        f"vector is the separate comparison: the {high_shock['symbol']} shock moves "
        f"from {high_shock['equal_weight_sensitivity']:.1%} to "
        f"{high_shock['network_weight_sensitivity']:.1%} on the network weights."
    )
)
```

Reach does not separate the two scenarios. The diffusion runs five rounds over
the thresholded graph counted above, which is dense enough that a shock starting
anywhere inside its connected part covers the same assets. What the starting node
changes is the path, and the path is what sets the size of each impact: an
impact is the initial shock multiplied by the correlations along the route that
reached the asset and by the attenuation at every hop, so two shocks that touch
the same assets deliver very different amounts to them.

That makes the scenario difference a statement about impact magnitudes under one
fixed weight vector, not about the weights. The weights are the other comparison,
and it runs within a scenario: the same shock evaluated on equal weights and on
network weights. Reading a reach difference as evidence that MST centrality
identifies systemic assets would in any case be reading the tree's construction
back out of the correlation matrix it was built from.

The table reports portfolio-weighted sensitivities, not a sum of hypothetical
asset losses. Synchronous updates make the result independent of loop order.

## 8. Performance Backtest

Apply the weights fixed at the estimation boundary to the later evaluation
window. `returns @ weights` holds the target weights every day, which is daily
rebalancing rather than buy and hold: the drift a held portfolio accumulates is
traded away each session. Returns are gross of trading costs and financing, and
daily rebalancing to a fixed target is exactly where those costs would arise.

```python
def backtest_portfolio(returns: np.ndarray, weights: np.ndarray) -> dict:
    """Gross evaluation metrics for a portfolio rebalanced daily to fixed weights."""
    portfolio_returns = returns @ weights

    # Annualized metrics
    ann_return = portfolio_returns.mean() * 252
    ann_vol = portfolio_returns.std() * np.sqrt(252)
    sharpe = ann_return / ann_vol if ann_vol > 0 else 0

    # Drawdown
    cumulative = (1 + portfolio_returns).cumprod()
    running_max = np.maximum.accumulate(np.concatenate(([1.0], cumulative)))[1:]
    drawdown = (cumulative - running_max) / running_max
    max_drawdown = drawdown.min()

    return {
        "ann_return": ann_return,
        "ann_volatility": ann_vol,
        "sharpe_ratio": sharpe,
        "max_drawdown": max_drawdown,
        "portfolio_returns": portfolio_returns,
        "growth": cumulative,
        "drawdown": drawdown,
    }
```

```python
print("BACKTEST RESULTS")

results = {}
result_rows = []
for name, weights in portfolios.items():
    metrics = backtest_portfolio(evaluation_returns, weights)
    results[name] = metrics
    result_rows.append(
        {
            "portfolio": name,
            "annualized_return": metrics["ann_return"],
            "annualized_volatility": metrics["ann_volatility"],
            "sharpe_ratio": metrics["sharpe_ratio"],
            "max_drawdown": metrics["max_drawdown"],
        }
    )
    assert all(np.isfinite(value) for key, value in result_rows[-1].items() if key != "portfolio")

performance = pl.DataFrame(result_rows)
print(performance)
```

### Evaluation Finding

This output-derived comparison reports the later window without turning one
historical split into a general performance claim.

```python
best_result = performance.sort("sharpe_ratio", descending=True).row(0, named=True)
network_result = performance.filter(pl.col("portfolio") == "Network Diversified").row(0, named=True)
equal_result = performance.filter(pl.col("portfolio") == "Equal Weight").row(0, named=True)
display(
    Markdown(
        f"**Finding:** {best_result['portfolio']} has the highest evaluation Sharpe "
        f"({best_result['sharpe_ratio']:.2f}). Network diversified records "
        f"{network_result['sharpe_ratio']:.2f} with a {network_result['max_drawdown']:.1%} "
        f"maximum drawdown, versus {equal_result['sharpe_ratio']:.2f} and "
        f"{equal_result['max_drawdown']:.1%} for equal weight."
    )
)
```

**Interpretation**: This is a single, later 252-session evaluation, not a model
selection exercise. The comparison is descriptive and does not establish that
one weighting rule will outperform in another period.

### Growth and Drawdown Comparison

Compare gross growth of one dollar and the associated drawdown paths.

```python
leader = performance.sort("sharpe_ratio", descending=True)["portfolio"][0]
strategy_colors = dict(zip(portfolios, ml4t_palette(3, categorical=True), strict=True))
fig, (ax_growth, ax_drawdown) = plt.subplots(2, 1, figsize=FIGSIZE["dual_v"], sharex=True)
for name, metrics in results.items():
    ax_growth.plot(evaluation_index, metrics["growth"], label=name, color=strategy_colors[name])
    ax_drawdown.plot(evaluation_index, metrics["drawdown"], color=strategy_colors[name])
    ax_drawdown.fill_between(
        evaluation_index, metrics["drawdown"], 0, color=strategy_colors[name], alpha=0.08
    )
ax_growth.set_ylabel("Growth of $1")
ax_drawdown.set_ylabel("Drawdown (%)")
ax_drawdown.set_xlabel("Evaluation date")
ax_drawdown.set_ylim(top=0)
format_pct_axis(ax_drawdown)
ax_growth.legend(loc="upper left")
add_message_title(
    ax_growth,
    "Growth of a dollar and drawdown, three weighting rules",
    subtitle=(
        f"{evaluation_start} to {evaluation_end}; gross returns, rebalanced daily to "
        "weights fixed at the estimation boundary"
    ),
)
show_with_alt(
    fig,
    f"Two stacked panels sharing a date axis over the evaluation window. The upper "
    f"panel plots growth of one dollar for {', '.join(portfolios)}, ending at "
    + ", ".join(f"{results[name]['growth'][-1]:.2f} for {name.lower()}" for name in portfolios)
    + f". The lower panel plots their drawdown paths below zero, with a maximum "
    f"drawdown across the three of "
    f"{abs(min(results[name]['max_drawdown'] for name in portfolios)) * 100:.1f} "
    "percent, in the final weeks of the window.",
)
```

### Centrality vs Portfolio Weight

The network-diversified rule mechanically underweights high-degree MST nodes.

```python
fig, ax = plt.subplots(figsize=FIGSIZE["single"])
ax.scatter(
    degree_centrality,
    w_network_div * 100,
    alpha=0.6,
    color=COLORS["blue"],
    s=30,
    label="Network Diversified",
)
ax.scatter(
    degree_centrality,
    np.full(n_assets, 100 / n_assets),
    alpha=0.4,
    color=COLORS["copper"],
    s=20,
    marker="x",
    label="Equal Weight",
)
# Label only the highest-degree hub to preserve the pattern without collisions.
hub_idx = int(np.argmax(degree_centrality))
ax.annotate(
    asset_names[hub_idx],
    (degree_centrality[hub_idx], w_network_div[hub_idx] * 100),
    xytext=(-4, 8),
    textcoords="offset points",
    fontsize=8,
    ha="right",
)
ax.set_xlabel("Degree Centrality")
ax.set_ylabel("Portfolio Weight (%)")
add_message_title(
    ax,
    "Portfolio weight against MST degree centrality",
    subtitle=f"weights estimated through {estimation_end}; cap verified at {MAX_WEIGHT:.0%}",
)
ax.legend()
show_with_alt(
    fig,
    f"A scatter of portfolio weight in percent against MST degree centrality for "
    f"{n_assets} assets. The network-diversified points trace a hyperbola: the "
    f"least connected assets sit at a weight of "
    f"{w_network_div.max() * 100:.1f} percent and the labelled hub "
    f"{asset_names[hub_idx]} at {w_network_div[hub_idx] * 100:.1f} percent. The "
    f"equal-weight crosses sit on a flat line at {100 / n_assets:.1f} percent.",
)
```

The inverse relationship is a direct consequence of
$w_i \propto 1/(c_i + 0.01)$. It visualizes the rule; it is not independent
evidence that centrality predicts future returns or losses.

## 9. Summary Statistics

```python
summary_stats = {
    "Data source": "Wiki Prices (real data)",
    "Assets analyzed": n_assets,
    "Estimation range": f"{estimation_start} to {estimation_end}",
    "Evaluation range": f"{evaluation_start} to {evaluation_end}",
    "Estimation days": ESTIMATION_DAYS,
    "Evaluation days": EVALUATION_DAYS,
    "MST edges": len(mst_edges),
    "Avg correlation": f"{corr_values.mean():.3f}",
    "Median correlation": f"{np.median(corr_values):.3f}",
    "Network portfolio HHI": f"{np.sum(w_network_div**2):.4f}",
    "Network portfolio Sharpe": f"{results['Network Diversified']['sharpe_ratio']:.2f}",
    "Equal weight Sharpe": f"{results['Equal Weight']['sharpe_ratio']:.2f}",
}

print("NOTEBOOK RESULTS FOR CHAPTER 23.5")

for key, value in summary_stats.items():
    print(f"  {key:35s}: {value}")
```

## 10. Visualization: Estimated MST

Visualize the estimation-window tree while labeling only its most central
nodes so the topology remains legible.

```python
# A seeded spring layout and one hub label expose structure without a label cloud.
positions = nx.spring_layout(mst_graph, seed=SEED, weight="weight")

centrality_cmap = LinearSegmentedColormap.from_list(
    "ml4t_centrality", [COLORS["silver_muted"], COLORS["blue"]]
)
fig, ax = plt.subplots(figsize=FIGSIZE["single_tall"])
nx.draw_networkx_edges(mst_graph, positions, ax=ax, edge_color=COLORS["silver_muted"], width=0.8)
nodes = nx.draw_networkx_nodes(
    mst_graph,
    positions,
    ax=ax,
    node_size=70 + 1_500 * degree_centrality,
    node_color=degree_centrality,
    cmap=centrality_cmap,
    edgecolors=COLORS["blue"],
    linewidths=0.4,
)
hub_x, hub_y = positions[hub_idx]
ax.annotate(
    asset_names[hub_idx],
    (hub_x, hub_y),
    xytext=(6, 6),
    textcoords="offset points",
    fontsize=8,
)
fig.colorbar(nodes, ax=ax, shrink=0.75, label="MST degree centrality")
add_message_title(
    ax,
    "The estimation-window minimum spanning tree",
    subtitle="node size and colour encode degree centrality; only the top hub is labelled",
)
ax.set_axis_off()
show_with_alt(
    fig,
    f"A spring-layout network diagram of the {n_assets}-asset minimum spanning tree "
    f"with {len(mst_edges)} edges. Node size and colour encode degree centrality, "
    f"which runs from {degree_centrality.min():.3f} to {degree_centrality.max():.3f}. "
    f"The largest and darkest node, labelled {asset_names[hub_idx]}, sits among "
    f"chains of small pale leaf nodes.",
)
```

Node size and color encode degree centrality in the estimation-window MST.
Larger nodes are the hubs the inverse-centrality rule mechanically underweights.

## 11. Verification

```python
print("NOTEBOOK EXECUTION COMPLETE")
print("Data: Wiki Prices (real market data)")
print(f"Assets: {n_assets}")
print(f"MST edges: {len(mst_edges)}")
print(f"Centrality range: [{degree_centrality.min():.3f}, {degree_centrality.max():.3f}]")
print(f"Estimation end: {estimation_end}")
print(f"Evaluation start: {evaluation_start}")
print(f"Network portfolio Sharpe: {results['Network Diversified']['sharpe_ratio']:.2f}")
print("Network-based portfolio evaluation complete.")
```

## Key Takeaways

1. **Correlation distance**: Converting correlations to distances via
   $d = \sqrt{2(1 - \rho)}$ satisfies the triangle inequality and enables
   MST construction that reveals hierarchical market structure.
2. **Centrality as a graph property**: Degree, betweenness, and closeness
   centrality describe how connected each asset is in the estimated MST. They
   do not establish causal influence or systemic importance.
3. **Separated evaluation**: Universe membership, the network, volatility,
   and all three weight vectors are fixed before the later 252-session gross
   evaluation. The output tables report the resulting concentration and risk.
4. **Shock sensitivity**: Synchronous diffusion makes the stylized scenario
   order-independent, and portfolio-weighted sensitivity avoids presenting a
   sum of hypothetical asset impacts as a portfolio loss.
5. **MST visualization**: The MST projects the dense correlation matrix into
   a tree that makes cluster structure and bridge nodes visually inspectable.
   This single-window comparison is descriptive, gross of costs, and not a
   claim of future outperformance.

**Next**: See `06_gnn_feature_engineering.py` for using GNN embeddings
derived from these networks as features for ML models.
**Book**: Chapter 23.5 discusses financial network theory and portfolio
applications in detail.
![notebook output](figures/p1_1.png)
![notebook output](figures/p1_2.png)
![notebook output](figures/p1_3.png)