برابری ریسک سلسله‌مراتبی: خوشه‌بندی و تخصیص بازگشتی سبد

# Hierarchical Risk Parity (HRP)


# Hierarchical Risk Parity (HRP)

**Docker image**: `ml4t`

This notebook demonstrates Hierarchical Risk Parity, a modern portfolio construction
method developed by Marcos López de Prado that reduces sensitivity to noisy covariance
estimates through clustering, quasi-diagonalization, and recursive bisection.

**Learning Objectives**:
- Understand why classical MVO can be fragile in practice (the "Markowitz Curse")
- Implement the three steps of HRP: clustering, quasi-diagonalization, recursive bisection
- Visualize the asset hierarchy with dendrograms
- Run walk-forward backtests comparing HRP to shrinkage MVO and to heuristic allocators, and
  read the ranking against the assets-to-observations ratio that decides when HRP helps

**Book Reference**: Chapter 17, Section 17.6 (Optimizing for stability with Hierarchical
Risk Parity)

**Prerequisites**: `02_mean_variance_optimization`, ETF price data

## 1. Setup and Imports

```python
"""Hierarchical Risk Parity: cluster assets and allocate using inverse variance."""

import hashlib
import warnings

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import plotly.express as px
import plotly.graph_objects as go
import polars as pl
from IPython.display import HTML, display
from ml4t.backtest import (
    BacktestConfig,
    CommissionType,
    DataFeed,
    Engine,
    ExecutionMode,
    Strategy,
)
from ml4t.backtest.config import SlippageType
from ml4t.backtest.execution.rebalancer import RebalanceConfig, TargetWeightExecutor
from ml4t.diagnostic.evaluation import PortfolioAnalysis
from ml4t.diagnostic.visualization import create_portfolio_dashboard
from plotly.subplots import make_subplots
from scipy.cluster.hierarchy import dendrogram, leaves_list, linkage
from scipy.spatial.distance import squareform
from sklearn.covariance import LedoitWolf

from case_studies.utils.backtest_loaders import compute_allocator_metrics
from case_studies.utils.registry.queries import load_prediction_index
from data import load_etfs
from utils.paths import get_case_study_dir, get_output_dir
from utils.reproducibility import set_global_seeds
from utils.style import COLORS, add_message_title, show_plotly_with_alt, show_with_alt
```

```python
# Production defaults; Papermill overrides these values for CI testing
START_DATE = "2010-01-01"
SEED = 42
```

```python
set_global_seeds(SEED)

# Track allocator fallbacks across the walk-forward backtest.
fallback_count: dict[str, int] = {}
```

## 2. Where Mean-Variance Optimization Is Fragile

Mean-variance optimization inverts an estimated covariance matrix. Inversion amplifies the
smallest eigenvalues, which are the ones the sample estimates worst, so small errors in the
inputs can move the weights a long way. The fragility grows with the ratio of assets to
observations: with $N$ assets and $T$ periods the sample covariance is singular once $N > T$
and poorly conditioned well before that.

HRP is one response - never invert. Shrinkage is another - keep inverting, but pull the
estimate toward a well-conditioned target first. This notebook runs both, and section 12 reads
the comparison against the ratio $N/T$. That ratio is one contributor to how hard the estimate
is, not the thing that decides which response to use: conditioning also depends on how the assets
co-move, and a low $N/T$ over highly correlated assets can be worse conditioned than a higher one
over independent ones. A single window over one universe cannot establish when either allocator
helps in general; what section 12 shows is which did better here. See Section 17.6.

## 3. Data Acquisition

The 15 named ETFs form a fixed teaching universe. It is not a point-in-time
membership reconstruction, so the static examples and walk-forward comparison are
conditional on this ex-post universe and must not be read as an unbiased universe-selection test.

```python
# Diversified ETF universe
ETF_UNIVERSE = {
    # US Equity
    "SPY": "S&P 500",
    "QQQ": "NASDAQ 100",
    "IWM": "Russell 2000",
    # International Equity
    "EFA": "EAFE Developed",
    "EEM": "Emerging Markets",
    # Fixed Income
    "AGG": "US Aggregate Bond",
    "TLT": "Long Treasury",
    "HYG": "High Yield",
    # Alternatives
    "GLD": "Gold",
    "VNQ": "Real Estate",
    "DBC": "Commodities",
    # Sectors
    "XLF": "Financials",
    "XLE": "Energy",
    "XLK": "Technology",
    "XLV": "Healthcare",
}

SYMBOLS = list(ETF_UNIVERSE.keys())
END_DATE = "2024-01-01"
```

```python
# Load data from canonical ETFs
print("Loading ETF data...")
etf_data = load_etfs()
etf_filtered = etf_data.filter(
    (pl.col("symbol").is_in(SYMBOLS))
    & (pl.col("timestamp") >= pl.lit(START_DATE).str.to_datetime())
    & (pl.col("timestamp") <= pl.lit(END_DATE).str.to_datetime())
)

close_prices = (
    etf_filtered.select(["timestamp", "symbol", "close"])
    .pivot(on="symbol", index="timestamp", values="close")
    .sort("timestamp")
    .to_pandas()
    .set_index("timestamp")
    .ffill()
    .dropna()
)
returns = close_prices.pct_change().dropna()
print(f"Loaded {len(returns):,} days for {close_prices.shape[1]} ETFs")
```

## 4. The HRP Algorithm

HRP works in three steps:

### Step 1: Tree Clustering
Build a hierarchy of assets based on their correlation structure using
agglomerative clustering.

$$d_{ij} = \sqrt{\frac{1-\rho_{ij}}{2}}$$

### Step 2: Quasi-Diagonalization
Reorder the covariance matrix according to the clustering hierarchy
to group similar assets together.

### Step 3: Recursive Bisection
Allocate risk by recursively splitting the portfolio in half,
with weights inversely proportional to the cluster variance.

```python
def correlation_distance(corr_matrix: np.ndarray) -> np.ndarray:
    """Convert a correlation matrix to a bounded distance matrix."""
    return np.sqrt(np.clip(0.5 * (1 - corr_matrix), 0.0, 1.0))
```

#### Step 1: Hierarchical Clustering

```python
def cluster_assets(returns: pd.DataFrame, method: str = "ward") -> np.ndarray:
    """
    Perform hierarchical clustering on assets based on correlation distance.

    Args:
        returns: DataFrame of asset returns
        method: Linkage method ('ward', 'single', 'complete', 'average')

    Returns:
        linkage matrix
    """
    corr = returns.corr().values
    dist = correlation_distance(corr)
    # Convert to condensed distance matrix (upper triangle)
    dist_condensed = squareform(dist, checks=False)
    return linkage(dist_condensed, method=method)
```

#### Step 2: Leaf Ordering

```python
def get_quasi_diagonal_order(link: np.ndarray) -> list[int]:
    """
    Get the quasi-diagonal ordering from linkage matrix.

    This reorders assets so that similar ones are adjacent.
    """
    return list(leaves_list(link))
```

The variance of a sub-portfolio held at inverse-variance weights, which is the number each
split compares its two halves on.

```python
def cluster_variance(cov: np.ndarray, indices: list[int]) -> float:
    """Compute cluster variance using inverse-variance weights within a subset."""
    if len(indices) == 1:
        return float(cov[indices[0], indices[0]])
    c = cov[np.ix_(indices, indices)]
    ivp = 1 / np.diag(c)
    ivp /= ivp.sum()
    return float(np.dot(ivp, np.dot(c, ivp)))
```

#### Step 3: Recursive Bisection Allocation

```python
def recursive_bisection(
    cov: np.ndarray,
    sorted_idx: list[int],
) -> np.ndarray:
    """Allocate weights by recursively splitting ordered clusters."""
    n = len(sorted_idx)
    weights = np.ones(n)
    sorted_position = {original: position for position, original in enumerate(sorted_idx)}
    clusters = [sorted_idx]

    while clusters:
        new_clusters = []
        for cluster in clusters:
            if len(cluster) <= 1:
                continue

            mid = len(cluster) // 2
            left = cluster[:mid]
            right = cluster[mid:]

            left_var = cluster_variance(cov, left)
            right_var = cluster_variance(cov, right)

            alpha = 1 - left_var / (left_var + right_var)
            weights[[sorted_position[i] for i in left]] *= alpha
            weights[[sorted_position[i] for i in right]] *= 1 - alpha

            if len(left) > 1:
                new_clusters.append(left)
            if len(right) > 1:
                new_clusters.append(right)

        clusters = new_clusters

    final_weights = np.zeros(n)
    final_weights[np.asarray(sorted_idx)] = weights

    return final_weights / final_weights.sum()
```

The three steps in order, from a covariance matrix to a weight vector.

```python
def hrp_portfolio(returns: pd.DataFrame) -> np.ndarray:
    """
    Compute HRP portfolio weights.

    Args:
        returns: DataFrame of asset returns

    Returns:
        Array of portfolio weights
    """
    # Step 1: Cluster
    link = cluster_assets(returns)

    # Step 2: Quasi-diagonalize
    sorted_idx = get_quasi_diagonal_order(link)

    # Step 3: Recursive bisection
    cov = returns.cov().values
    weights = recursive_bisection(cov, sorted_idx)

    return weights
```

## 5. Visualizing the Asset Hierarchy

```python
# Compute linkage for visualization
link = cluster_assets(returns)
branch_colors = [COLORS["blue"], COLORS["copper"], COLORS["positive"], COLORS["neutral"]]

# Create dendrogram
fig, ax = plt.subplots(figsize=(14, 8))
dendrogram(
    link,
    labels=[ETF_UNIVERSE.get(s, s) for s in returns.columns],
    leaf_rotation=45,
    leaf_font_size=10,
    link_color_func=lambda node_id: branch_colors[(node_id - len(returns.columns)) % 4],
    ax=ax,
)
ax.set_ylabel("Distance (based on correlation)")
ax.set_xlabel("ETF")
add_message_title(
    ax,
    "Hierarchical clustering of fifteen ETFs on correlation distance",
    subtitle="Ward linkage, fixed teaching universe",
)
show_with_alt(
    fig,
    "Dendrogram of fifteen ETFs on correlation distance, with the bond and gold funds joining low on one branch and the equity and sector funds on another.",
)
```

Read the tree from the leaves upward. Two ETFs joined low down moved together over this
history. The height of a merge is SciPy's Ward distance, which is a monotone transformation of
the increase in within-cluster sum of squares the merge costs, so it rises as the two groups
being joined become less alike. It is not the correlation distance between them, which is what
the leaves were measured on.

Only the leaf order reaches the allocation. Step 2 takes the left-to-right sequence this
tree implies and reorders the covariance matrix by it; step 3 then halves that sequence by
count at every level, without consulting the heights again. The halves it forms therefore need
not be the clusters the tree draws, which is what section 7 traces through the first two
splits.

## 6. Quasi-Diagonal Covariance Matrix

After clustering, we reorder the covariance matrix to reveal the block structure.

```python
# Get quasi-diagonal ordering
sorted_idx = get_quasi_diagonal_order(link)
sorted_symbols = [returns.columns[i] for i in sorted_idx]
sorted_labels = [ETF_UNIVERSE.get(s, s) for s in sorted_symbols]

# Reorder covariance matrix
cov_original = returns.cov()
cov_reordered = cov_original.iloc[sorted_idx, sorted_idx]
cov_colorscale = [
    [0.0, COLORS["negative"]],
    [0.5, COLORS["silver"]],
    [1.0, COLORS["blue"]],
]
```

```python
fig = make_subplots(
    rows=1, cols=2, subplot_titles=["Original Covariance Matrix", "Quasi-Diagonal (Reordered)"]
)
fig.add_trace(
    go.Heatmap(
        z=cov_original.values,
        x=cov_original.columns,
        y=cov_original.index,
        colorscale=cov_colorscale,
        zmid=0,
        showscale=False,
    ),
    row=1,
    col=1,
)
fig.add_trace(
    go.Heatmap(
        z=cov_reordered.values,
        x=sorted_symbols,
        y=sorted_symbols,
        colorscale=cov_colorscale,
        zmid=0,
        showscale=True,
    ),
    row=1,
    col=2,
)
fig.update_layout(
    title="Covariance matrix in the original and the quasi-diagonal ordering",
    height=560,
    width=1100,
    margin=dict(l=90, r=90, b=110, t=100),
)
fig.update_xaxes(tickangle=45, tickfont_size=10, automargin=True)
fig.update_yaxes(tickfont_size=10, automargin=True)
fig.update_xaxes(title_text="ETF ticker", row=1, col=1)
fig.update_xaxes(title_text="ETF ticker", row=1, col=2)
fig.update_yaxes(title_text="ETF ticker", row=1, col=1)
show_plotly_with_alt(
    fig,
    "Two covariance heatmaps side by side, the original ordering on the left and the quasi-diagonal reordering on the right, in which the large values group into blocks along the diagonal.",
)
```

## 7. Compute HRP Weights

```python
# Compute HRP weights
hrp_weights = hrp_portfolio(returns)

# Create weights DataFrame
weights_df = pd.DataFrame(
    {
        "Symbol": returns.columns,
        "Name": [ETF_UNIVERSE.get(s, s) for s in returns.columns],
        "HRP Weight": hrp_weights,
    }
)
weights_df = weights_df.sort_values("HRP Weight", ascending=False)
```

```python
# Visualization
fig = px.bar(
    weights_df,
    x="Name",
    y="HRP Weight",
    title="HRP weight per ETF, largest first",
    color="HRP Weight",
    color_continuous_scale=[COLORS["silver_muted"], COLORS["blue"]],
)
fig.update_layout(
    height=400,
    xaxis_tickangle=45,
    xaxis_title="ETF",
    yaxis_title="Portfolio weight",
    yaxis_tickformat=".0%",
)
show_plotly_with_alt(
    fig,
    "Bars of HRP weight per ETF, sorted from largest to smallest, with one bond fund holding the majority of the portfolio and the rest falling away sharply.",
)
```

### Where That Concentration Comes From

HRP is introduced as the answer to the concentrated portfolios that mean-variance optimization
produces, so a single dominant weight deserves an explanation rather than a shrug. It follows
from the algorithm as specified.

Recursive bisection splits the *ordered list* in half by count at every level. It does not cut
the tree at the linkage's own cluster boundaries, so the halves it forms need not be the
clusters the dendrogram shows. Each half then receives capital inversely to its cluster
variance. Two splits are enough to see the effect.

```python
def trace_bisection(cov: np.ndarray, sorted_idx: list[int], columns, levels: int = 2) -> None:
    """Print the halves, their annualized cluster variances, and the resulting split share."""
    clusters, depth = [sorted_idx], 0
    while clusters and depth < levels:
        next_clusters = []
        for cluster in clusters:
            if len(cluster) <= 1:
                continue
            mid = len(cluster) // 2
            left, right = cluster[:mid], cluster[mid:]
            left_var, right_var = cluster_variance(cov, left), cluster_variance(cov, right)
            alpha = 1 - left_var / (left_var + right_var)
            print(f"  split at depth {depth}:")
            print(
                f"    {[columns[i] for i in left]}\n"
                f"      annualized cluster variance {left_var * 252:.4f} -> {alpha:.1%} of the branch"
            )
            print(
                f"    {[columns[i] for i in right]}\n"
                f"      annualized cluster variance {right_var * 252:.4f} -> {1 - alpha:.1%} of the branch"
            )
            if len(left) > 1:
                next_clusters.append(left)
            if len(right) > 1:
                next_clusters.append(right)
        clusters = next_clusters
        depth += 1


print("First two levels of the recursive bisection:")
trace_bisection(returns.cov().values, get_quasi_diagonal_order(link), list(returns.columns))
print()
print(f"Weight in the single largest holding: {weights_df['HRP Weight'].max():.1%}")
print(f"Weight in the three largest: {weights_df['HRP Weight'].nlargest(3).sum():.1%}")
```

The low-variance half takes most of the capital at both splits, and the shares multiply. Bonds and gold
end up holding most of the portfolio, and within them the lowest-volatility holding takes most
of what is left.

This is inverse-variance allocation doing exactly what it is defined to do, not a bug. But it
means HRP is not concentration-free: it concentrates on the *low-variance* assets rather than on
whichever assets the covariance inverse happens to favor. Whether that is an improvement depends
on whether low realized variance in the estimation window is a better guide to the future than
the inverse covariance is - a question the walk-forward comparison below can address and this
static example cannot.

## 8. Comparison: HRP vs Other Methods

```python
def inverse_volatility_weights(returns: pd.DataFrame) -> np.ndarray:
    """Inverse volatility (risk parity) weights."""
    vols = returns.std().values
    inv_vols = 1 / vols
    return inv_vols / inv_vols.sum()
```

#### Equal-Weight Baseline

```python
def equal_weights(n: int) -> np.ndarray:
    """Equal weights."""
    return np.ones(n) / n
```

#### Minimum-Variance Baseline

The long-only projection starts from
$w = \Sigma^{-1}\mathbf{1}/(\mathbf{1}^{\top}\Sigma^{-1}\mathbf{1})$,
clips negative weights, and renormalizes the remaining capital.

```python
def minimum_variance_weights(returns: pd.DataFrame) -> np.ndarray:
    """Minimum variance portfolio (no expected returns input)."""
    cov = returns.cov().values
    n = len(returns.columns)

    try:
        cov_inv = np.linalg.inv(cov)
        ones = np.ones(n)
        weights = cov_inv @ ones
        weights /= weights.sum()
        # Long-only projection: clip negatives, then renormalize so weights sum to 1
        weights = np.clip(weights, 0, None)
        total = weights.sum()
        return weights / total if total > 0 else equal_weights(n)
    except np.linalg.LinAlgError:
        return equal_weights(n)
```

#### Ledoit-Wolf Shrinkage Baseline

```python
def mvo_shrinkage_weights(returns: pd.DataFrame) -> np.ndarray:
    """MVO with Ledoit-Wolf shrinkage."""
    # Use Ledoit-Wolf for robust covariance
    lw = LedoitWolf().fit(returns)
    cov_shrunk = lw.covariance_

    # Minimum variance with shrunk covariance
    n = len(returns.columns)
    try:
        cov_inv = np.linalg.inv(cov_shrunk)
        ones = np.ones(n)
        weights = cov_inv @ ones
        weights /= weights.sum()
        # Long-only projection: clip negatives, then renormalize so weights sum to 1
        weights = np.clip(weights, 0, None)
        total = weights.sum()
        return weights / total if total > 0 else equal_weights(n)
    except np.linalg.LinAlgError:
        return equal_weights(n)
```

```python
# Compute all weights
n_assets = len(returns.columns)

all_weights = pd.DataFrame(
    {
        "Symbol": returns.columns,
        "Name": [ETF_UNIVERSE.get(s, s) for s in returns.columns],
        "Equal": equal_weights(n_assets),
        "Inv Vol": inverse_volatility_weights(returns),
        "Min Var (LW)": mvo_shrinkage_weights(returns),
        "HRP": hrp_weights,
    }
)
```

```python
# Visualization: Weight comparison
fig = go.Figure()

# One colour per allocator, in the order the four are always listed. Equal weight is the
# benchmark, so it takes the neutral grey; HRP is the subject of the notebook and takes the
# primary colour. The same four are reused for the equity curves in section 9.
methods = ["Equal", "Inv Vol", "Min Var (LW)", "HRP"]
ALLOCATOR_COLORS = [COLORS["neutral"], COLORS["amber"], COLORS["copper"], COLORS["blue"]]

for method, color in zip(methods, ALLOCATOR_COLORS, strict=True):
    fig.add_trace(
        go.Bar(
            name=method,
            x=all_weights["Name"],
            y=all_weights[method],
            marker_color=color,
        )
    )

fig.update_layout(
    title="Portfolio weight per ETF under four allocators",
    barmode="group",
    xaxis_tickangle=45,
    xaxis_title="ETF",
    yaxis_title="Portfolio weight",
    yaxis_tickformat=".0%",
    height=500,
    legend=dict(orientation="h", yanchor="bottom", y=1.02),
)
show_plotly_with_alt(
    fig,
    "Grouped bars of portfolio weight per ETF for equal weight, inverse volatility, shrinkage "
    "minimum variance and HRP. The two variance-based allocators put the bulk of the "
    "portfolio in the aggregate bond fund, while the equal-weight bars sit at the same "
    "height across the universe.",
)
```

## 9. Walk-Forward Backtest with ML Predictions

The previous sections demonstrated *how* HRP allocates capital given a
covariance matrix. We now apply it inside a realistic trading strategy:

1. Load the ETF case-study GBM walk-forward predictions (Chapter 12).
2. At each rebalance date, use the latest available prediction to *select* the
   top-$N$ assets.
3. Apply HRP (and the other methods) to the selected subset using only
   historical returns for covariance estimation.

Asset selection is driven by a walk-forward validation signal, while the
allocator only sees information available at the rebalance date. The target
becomes effective on the following bar.

Asset selection uses the ETF case study's best-IC GBM walk-forward validation predictions. HRP
itself needs only a covariance matrix, so the predictions are not part of the allocator; they
decide *which* assets each allocator sizes, and every allocator sees the same selection.

The highest-validation-IC GBM checkpoint is resolved at runtime from
`case_studies/etfs/run_log/registry.db`. No `prediction_hash` is baked in, so re-running the
GBM sweep changes which predictions feed the comparison without an edit here.

The registry ranks and evaluates on validation, so what follows demonstrates allocation
behavior on a sample the selection step has already seen. It is not an out-of-sample estimate.

```python
etf_case_dir = get_case_study_dir("etfs")
etf_registry_path = etf_case_dir / "run_log" / "registry.db"
etf_registry_sha256 = hashlib.sha256(etf_registry_path.read_bytes()).hexdigest()
prediction_index = load_prediction_index(
    "etfs",
    label="fwd_ret_21d",
    split="validation",
    family="gbm",
    case_dir=etf_case_dir,
)
if prediction_index.is_empty():
    raise RuntimeError("No registered ETF GBM validation predictions are available")
best_ic = prediction_index["ic_mean"][0]
if best_ic is None or not np.isfinite(best_ic):
    raise RuntimeError("The leading ETF GBM validation prediction has no finite IC")
leaders = prediction_index.filter((pl.col("ic_mean") - best_ic).abs() <= 1e-12)
if leaders.height != 1:
    hashes = leaders["prediction_hash"].to_list()
    raise RuntimeError(f"Ambiguous best-IC ETF GBM predictions: {hashes}")
best_gbm = leaders.row(0, named=True)
ETF_GBM_PRED_HASH = best_gbm["prediction_hash"]
```

Validate the selected prediction parquet and record immutable input hashes.

```python
PRED_PATH = etf_case_dir / "run_log" / "predictions" / ETF_GBM_PRED_HASH / "predictions.parquet"
if not PRED_PATH.exists():
    raise FileNotFoundError(
        f"Resolved best-IC GBM hash {ETF_GBM_PRED_HASH} (config={best_gbm['config_name']}) "
        f"has no predictions parquet at {PRED_PATH}. Re-run case_studies/etfs/07_gbm.py."
    )
print(
    f"Resolved best-IC GBM: hash={ETF_GBM_PRED_HASH}, config={best_gbm['config_name']}, "
    f"IC={best_gbm['ic_mean']:.4f}"
)
print(f"ETF registry SHA-256: {etf_registry_sha256}")
print(f"Prediction parquet SHA-256: {hashlib.sha256(PRED_PATH.read_bytes()).hexdigest()}")

upstream_preds = pl.read_parquet(PRED_PATH).select("timestamp", "symbol", "prediction")
print(
    f"Loaded ETF GBM predictions: {upstream_preds.height:,} rows, "
    f"{upstream_preds['symbol'].n_unique()} symbols, "
    f"{upstream_preds['timestamp'].min()} to {upstream_preds['timestamp'].max()}"
)
```

### Asset Selection Rule

At each rebalance date we take the most recent prediction available for each
symbol (as-of join), restrict the ranking to the fixed teaching universe, and
select the top $N$. If no prediction is available
yet for the date (e.g., before the first walk-forward fold), the rebalance
is skipped by returning an empty list.

```python
def select_top_assets(
    date: pd.Timestamp,
    full_returns: pd.DataFrame,
    predictions: pl.DataFrame,
    top_n: int,
) -> list[str]:
    """Select top-N assets by latest available GBM prediction at `date`.

    Polars `group_by` does not guarantee row order, and ties at the top-N
    cutoff would otherwise produce non-deterministic selections across runs.
    We secondary-sort by symbol so the selection is reproducible.
    """
    as_of = predictions.filter(pl.col("timestamp") <= date)
    if as_of.is_empty():
        return []
    latest = (
        as_of.sort("timestamp")
        .group_by("symbol")
        .agg(pl.col("prediction").last())
        .filter(pl.col("symbol").is_in(full_returns.columns.to_list()))
        .filter(pl.col("prediction").is_finite())
        .sort(["prediction", "symbol"], descending=[True, False])
    )
    return latest.head(top_n)["symbol"].to_list()
```

### Guarding against a covariance the allocator cannot use

Validate each allocator's long-only weights before mapping them to the full universe. A singular
covariance or invalid weight vector triggers a visible equal-weight fallback and increments the
method's fallback counter.

```python
def allocate_selected_assets(
    full_columns: pd.Index,
    top_assets: list[str],
    selected_returns: pd.DataFrame,
    allocation_fn,
    method: str = "allocator",
    date: pd.Timestamp | None = None,
) -> pd.Series:
    """Map valid selected weights to the full universe, with a reported equal-weight fallback."""
    weights = pd.Series(0.0, index=full_columns)
    try:
        alloc_weights = np.asarray(allocation_fn(selected_returns), dtype=float)
        if alloc_weights.shape != (len(top_assets),):
            raise ValueError(
                f"allocator returned shape {alloc_weights.shape}, expected {(len(top_assets),)}"
            )
        if not np.isfinite(alloc_weights).all() or (alloc_weights < 0).any():
            raise ValueError("allocator returned non-finite or negative weights")
        total = alloc_weights.sum()
        if total <= 0:
            raise ValueError("allocator returned non-positive total weight")
        alloc_weights = alloc_weights / total
        for asset, weight in zip(top_assets, alloc_weights, strict=False):
            weights[asset] = weight
    except (np.linalg.LinAlgError, ValueError) as exc:
        warnings.warn(
            f"{method} singular/invalid for {len(top_assets)} assets at {date}: {exc}",
            RuntimeWarning,
            stacklevel=2,
        )
        fallback_count[method] = fallback_count.get(method, 0) + 1
        for asset in top_assets:
            weights[asset] = 1.0 / len(top_assets)
    return weights
```

### Rebalance Decision

Build a target from trailing returns and the latest registered prediction. Returning
`None` leaves the prior target unchanged and prevents a pre-signal baseline from entering results.

```python
def rebalance_target(
    returns: pd.DataFrame,
    predictions: pl.DataFrame,
    date: pd.Timestamp,
    allocation_fn,
    top_n: int,
    lookback: int,
    method: str,
    min_history: int,
) -> pd.Series | None:
    """Return the target decided at `date`, or None when history/signal is unavailable."""
    loc = returns.index.get_loc(date)
    hist_returns = returns.iloc[max(0, loc - lookback) : loc]
    if len(hist_returns) < min_history:
        return None
    top_assets = select_top_assets(date, returns, predictions, top_n)
    if len(top_assets) < 2:
        return None
    return allocate_selected_assets(
        returns.columns,
        top_assets,
        hist_returns[top_assets],
        allocation_fn,
        method=method,
        date=date,
    )
```

### Walk-Forward Backtest Engine

```python
def walk_forward_backtest(
    returns: pd.DataFrame,
    allocation_fn,
    predictions: pl.DataFrame,
    lookback: int,
    rebalance_freq: str,
    top_n: int,
    min_history: int,
    method: str = "allocator",
) -> tuple[pd.DataFrame, pd.DataFrame]:
    """Hold positions between month-end decisions and let weights drift with returns."""
    rebalance_dates = set(returns.groupby(returns.index.to_period(rebalance_freq)).tail(1).index)
    portfolio_returns, target_history = [], []
    current_weights = pending_target = None
    for date in returns.index:
        # A target decided at the prior close becomes the next session's opening allocation.
        if pending_target is not None:
            current_weights = pending_target
            pending_target = None
        if current_weights is not None:
            asset_returns = returns.loc[date]
            daily_ret = float((asset_returns * current_weights).sum())
            portfolio_returns.append({"date": date, "return": daily_ret})
            end_values = current_weights * (1.0 + asset_returns)
            current_weights = end_values / end_values.sum()
        if date in rebalance_dates:
            target = rebalance_target(
                returns, predictions, date, allocation_fn, top_n, lookback, method, min_history
            )
            if target is not None:
                pending_target = target
                target_history.append(
                    {"date": date, **{f"w_{s}": target[s] for s in returns.columns}}
                )

    return (
        pd.DataFrame(portfolio_returns).set_index("date"),
        pd.DataFrame(target_history).set_index("date"),
    )
```

Four settings decide what the backtest does, and every allocator gets the same four.
`LOOKBACK` is how much trailing history each covariance is estimated on - one year of
sessions, which with five selected names leaves the sample covariance comfortably
over-determined. `TOP_N` is how many of the fifteen the prediction picks each month, and it is
the number that decides whether the ratio of assets to observations is anywhere near the
regime HRP is designed for. `REBALANCE_FREQ` is month-end, so the targets are decided twelve
times a year and held in between. `MIN_HISTORY` is the shortest trailing window an allocator
is allowed to size from at all; below it the rebalance is skipped rather than estimated on too
few rows.

```python
LOOKBACK = 252
TOP_N = 5
REBALANCE_FREQ = "M"
MIN_HISTORY = 60

allocation_methods = {
    "Equal Weight": lambda r: equal_weights(len(r.columns)),
    "Inverse Volatility": inverse_volatility_weights,
    "Min Variance (LW)": mvo_shrinkage_weights,
    "HRP": hrp_portfolio,
}

results = {}
weights_all = {}

print("Running walk-forward backtests...")
for name, alloc_fn in allocation_methods.items():
    print(f"  {name}...")
    ret_df, w_df = walk_forward_backtest(
        returns=returns,
        allocation_fn=alloc_fn,
        predictions=upstream_preds,
        lookback=LOOKBACK,
        rebalance_freq=REBALANCE_FREQ,
        top_n=TOP_N,
        method=name,
        min_history=MIN_HISTORY,
    )
    results[name] = ret_df["return"]
    weights_all[name] = w_df

print("All backtests complete")

# Combine into DataFrame
portfolio_returns = pd.DataFrame(results)
```

### Gross Backtest Metrics

The vectorized comparison isolates allocation effects and is gross of implementation costs.
Turnover is reported beside performance. The execution-aware bridge below applies its declared
commission and slippage assumptions and uses next-bar fills.

```python
def _weights_wide_to_long(weights_df: pd.DataFrame) -> pl.DataFrame:
    """Convert dense target weights to canonical timestamp/symbol/weight format."""
    wide = weights_df.reset_index()
    ts_col = wide.columns[0]
    wide = wide.rename(columns={ts_col: "timestamp"})
    long = wide.melt(id_vars=["timestamp"], var_name="symbol", value_name="weight")
    long["symbol"] = long["symbol"].str.replace("w_", "", regex=False)
    return pl.from_pandas(long[["timestamp", "symbol", "weight"]]).sort(["timestamp", "symbol"])
```

Compute allocator metrics for each method. The dense month-end target state preserves both new
positions and explicit zero-weight liquidations before the per-symbol turnover difference.

```python
allocator_metrics: dict[str, dict] = {}
metrics_list = []
for name in portfolio_returns.columns:
    returns_arr = portfolio_returns[name].to_numpy()
    monthly_targets = weights_all[name]
    weights_long = _weights_wide_to_long(monthly_targets)
    m = compute_allocator_metrics(
        pl.Series("returns", returns_arr),
        weights_df=weights_long,
        ann_factor=np.sqrt(252),
    )
    allocator_metrics[name] = m
    finite = returns_arr[np.isfinite(returns_arr)]
    annual_vol = float(np.nanstd(finite, ddof=1) * np.sqrt(252)) if finite.size > 1 else 0.0
    metrics_list.append(
        {
            "Method": name,
            "Annual Return": m["annual_return"],
            "Annual Vol": annual_vol,
            "Sharpe Ratio": m["sharpe"],
            "Max Drawdown": m["max_drawdown"],
            "Calmar Ratio": m["calmar"],
            "Avg Turnover": m["avg_turnover"],
        }
    )
```

```python
metrics_df = pd.DataFrame(metrics_list)
metrics_df = metrics_df[
    [
        "Method",
        "Annual Return",
        "Annual Vol",
        "Sharpe Ratio",
        "Max Drawdown",
        "Calmar Ratio",
        "Avg Turnover",
    ]
]
metrics_df.round(4)
```

### Scheduled Weight Strategy

This adapter converts the notebook's target-weight schedule into the
event-driven strategy interface expected by `ml4t-backtest`.

```python
class ScheduledWeightStrategy(Strategy):
    def __init__(self, weights_long: pl.DataFrame, allow_short: bool):
        self.executor = TargetWeightExecutor(
            config=RebalanceConfig(
                min_trade_value=100.0,
                min_weight_change=0.001,
                allow_fractional=True,
                allow_short=allow_short,
            )
        )
        self._targets_by_ts: dict[pd.Timestamp, dict[str, float]] = {}
        for row in weights_long.iter_rows(named=True):
            ts = pd.Timestamp(row["timestamp"]).tz_localize(None)
            self._targets_by_ts.setdefault(ts, {})
            self._targets_by_ts[ts][str(row["symbol"])] = float(row["weight"])

    def on_data(self, timestamp, data, context, broker):
        ts = pd.Timestamp(timestamp).tz_localize(None)
        targets = self._targets_by_ts.get(ts)
        if not targets:
            return
        targets = {asset: weight for asset, weight in targets.items() if asset in data}
        if targets:
            self.executor.execute(targets, data, broker)
```

### Restating the schedule in the form the engine reads

Convert the wide monthly target schedule into the long-form price and target tables required by
the execution engine. Zero targets remain present so liquidations are explicit.

```python
bridge_method = "HRP" if "HRP" in weights_all else next(iter(weights_all))
weights_long = _weights_wide_to_long(weights_all[bridge_method]).with_columns(
    pl.col("timestamp").cast(pl.Datetime("us"))
)
allow_short_engine = (
    weights_long.filter(pl.col("weight") < 0).height > 0 if not weights_long.is_empty() else False
)

prices_panel = pl.from_pandas(close_prices.reset_index())
ts_col = prices_panel.columns[0]
if ts_col != "timestamp":
    prices_panel = prices_panel.rename({ts_col: "timestamp"})
prices_long = (
    prices_panel.unpivot(index="timestamp", variable_name="symbol", value_name="close")
    .with_columns(
        [
            pl.col("timestamp").cast(pl.Datetime("us")),
            pl.col("close").alias("open"),
            pl.col("close").alias("high"),
            pl.col("close").alias("low"),
            pl.lit(1_000_000).alias("volume"),
        ]
    )
    .sort(["timestamp", "symbol"])
)
```

### Run the Execution-Aware Backtest

Replay the same target weights with commissions, slippage, and next-bar fills. Synthetic
OHLC fields equal the daily close, so this is an accounting and timing bridge rather than an
intraday fill-quality model.

```python
engine = Engine(
    feed=DataFeed(prices_df=prices_long),
    strategy=ScheduledWeightStrategy(weights_long, allow_short=allow_short_engine),
    config=BacktestConfig(
        initial_cash=100_000.0,
        execution_mode=ExecutionMode.NEXT_BAR,
        commission_type=CommissionType.PERCENTAGE,
        commission_rate=0.0005,
        slippage_type=SlippageType.PERCENTAGE,
        slippage_rate=0.0005,
        allow_short_selling=allow_short_engine,
    ),
)
```

### Compare Engine and Vectorized Returns

The bridge checks whether the event-driven backtest preserves the ranking and
risk shape implied by the vectorized walk-forward simulation.

```python
engine_daily = (
    engine.run()
    .to_daily_pnl()
    .select(
        pl.col("date").cast(pl.Datetime("us")).alias("date"),
        pl.col("return_pct").alias("engine_return"),
    )
)
vectorized_daily = pl.DataFrame(
    {
        "date": pl.Series(portfolio_returns.index.to_list()).cast(pl.Datetime("us")),
        "vectorized_return": portfolio_returns[bridge_method].to_numpy(),
    }
)

bridge = (
    vectorized_daily.join(engine_daily, on="date", how="inner")
    .drop_nulls(["vectorized_return", "engine_return"])
    .sort("date")
)

vec_summary = PortfolioAnalysis(returns=bridge["vectorized_return"], periods_per_year=252)
eng_summary = PortfolioAnalysis(returns=bridge["engine_return"], periods_per_year=252)
vec_stats = vec_summary.compute_summary_stats()
eng_stats = eng_summary.compute_summary_stats()

print(f"\nExecution bridge ({bridge_method} walk-forward):")
print(
    f"  Vectorized Sharpe={vec_stats.sharpe_ratio:.3f}, Engine Sharpe={eng_stats.sharpe_ratio:.3f}"
)
print(f"  Vectorized MaxDD={vec_stats.max_drawdown:.2%}, Engine MaxDD={eng_stats.max_drawdown:.2%}")
```

```python
# Equity curves
cumulative = (1 + portfolio_returns).cumprod()

fig = go.Figure()

for col, color in zip(cumulative.columns, ALLOCATOR_COLORS, strict=True):
    fig.add_trace(
        go.Scatter(
            x=cumulative.index,
            y=cumulative[col],
            mode="lines",
            name=col,
            line=dict(color=color, width=2 if col == "HRP" else 1.5),
        )
    )

fig.add_hline(y=1.0, line_dash="dot", line_color=COLORS["neutral"])

fig.update_layout(
    title="Four allocators over the same monthly selection, gross of costs",
    xaxis_title="Date",
    yaxis_title="Growth of $1",
    height=500,
)
show_plotly_with_alt(
    fig,
    "Four growth-of-one-dollar paths from the walk-forward backtest, one per allocation "
    "method. They track each other until about 2018 and then separate into two pairs, equal "
    "weight and inverse volatility above, shrinkage minimum variance and HRP below, each "
    "falling sharply in early 2020.",
)
```

## 10. The HRP series on its own terms

The table above compares allocators to each other. The dashboard below looks at one of them -
HRP - the way an investor would: cumulative growth against SPY, the underwater curve, rolling
Sharpe, and the monthly return grid. Its metrics come from the same return series as the table,
so the two agree to rounding; what it adds is the path behind the summary statistics.

```python
hrp_returns_series = portfolio_returns["HRP"].dropna()
spy_returns_series = returns["SPY"].reindex(hrp_returns_series.index).dropna()
common_idx = hrp_returns_series.index.intersection(spy_returns_series.index)

hrp_analysis = PortfolioAnalysis(
    returns=hrp_returns_series.loc[common_idx].to_numpy(),
    benchmark=spy_returns_series.loc[common_idx].to_numpy(),
    dates=common_idx,
    risk_free=0.0,
    periods_per_year=252,
)

hrp_tear_sheet = create_portfolio_dashboard(hrp_analysis)

for dashboard_figure in hrp_tear_sheet.figures.values():
    dashboard_figure.update_layout(
        paper_bgcolor=COLORS["bg_light"],
        plot_bgcolor=COLORS["bg_light"],
        font_color=COLORS["neutral"],
    )

rolling_beta_figure = hrp_tear_sheet.figures["Rolling Beta"]
rolling_beta_figure.update_layout(margin=dict(l=60, r=90, t=40, b=40))
_ = rolling_beta_figure.update_annotations(x=0.995, xanchor="right")
```

`hrp_tear_sheet.show()` would render the metrics block and then loop `fig.show()` over the
nine figures, which publishes each PNG with no alt text and leaves a screen reader with
nothing. Displaying the same content a figure at a time is what lets each one carry a
description of what it plots.

```python
DASHBOARD_ALT = {
    "Cumulative Returns": (
        "Cumulative return of the HRP portfolio and of the SPY benchmark against date, the "
        "benchmark drawn dashed and finishing well above the portfolio."
    ),
    "Drawdown": (
        "The HRP portfolio's underwater curve against date, filled to zero, showing the "
        "percentage below its own running peak, with the deepest point marked in March 2020."
    ),
    "Rolling Sharpe Ratio": (
        "Two lines of rolling Sharpe ratio against date, over sixty-three and two hundred and "
        "fifty-two sessions, the shorter window swinging more widely than the longer one."
    ),
    "Rolling Volatility": (
        "Three lines of annualized rolling volatility against date, over twenty-one, "
        "sixty-three and two hundred and fifty-two sessions."
    ),
    "Rolling Beta": (
        "Rolling beta of the HRP portfolio against SPY, plotted against date and shaded down "
        "to zero, with a dashed reference line at the market's own beta."
    ),
    "Annual Returns": (
        "Bars of the HRP portfolio's annual return by calendar year, with the benchmark's "
        "annual return marked as points and a line at zero."
    ),
    "Monthly Returns Heatmap": (
        "Heatmap of monthly return, years down the vertical axis and calendar months across "
        "with a compounded annual column at the right, each cell labelled with its return and "
        "coloured from red for losses to green for gains."
    ),
    "Returns Distribution": (
        "Histogram of the HRP portfolio's daily returns with a fitted normal density drawn "
        "over it and vertical reference lines in the left tail."
    ),
    "Top Drawdowns": (
        "Horizontal bars of drawdown depth for the five deepest episodes, one bar per episode, "
        "ordered deepest at the top and annotated with the depth reached."
    ),
}

display(HTML(f"<pre>{hrp_tear_sheet.metrics.summary()}</pre>"))
for figure_name, dashboard_figure in hrp_tear_sheet.figures.items():
    show_plotly_with_alt(dashboard_figure, DASHBOARD_ALT[figure_name])
```

```python
# HTML delivery: the same content as a single self-contained file.
output_dir = get_output_dir(17, "hrp")
output_dir.mkdir(parents=True, exist_ok=True)
hrp_tear_sheet_path = output_dir / "hrp_tear_sheet.html"
hrp_tear_sheet.save_html(hrp_tear_sheet_path, include_plotlyjs="cdn")
print(f"HRP tear sheet saved: {output_dir.name}/{hrp_tear_sheet_path.name}")
print(f"  Figures embedded: {list(hrp_tear_sheet.figures.keys())}")
```

## 11. Weight Evolution in Strategy Context

Analyze how the submitted HRP targets change over the walk-forward backtest. The portfolio holds
each allocation until the next month-end decision, so realized weights drift between targets.

```python
# Plot HRP weight evolution from actual backtest
hrp_weights_hist = weights_all.get("HRP", pd.DataFrame())

if not hrp_weights_hist.empty:
    # Get weight columns
    weight_cols = [c for c in hrp_weights_hist.columns if c.startswith("w_")]

    fig = go.Figure()

    # Show top assets by average weight
    avg_weights = hrp_weights_hist[weight_cols].mean().sort_values(ascending=False)
    top_cols = avg_weights.head(6).index.tolist()

    for col in top_cols:
        symbol = col.replace("w_", "")
        fig.add_trace(
            go.Scatter(
                x=hrp_weights_hist.index,
                y=hrp_weights_hist[col],
                mode="lines",
                name=ETF_UNIVERSE.get(symbol, symbol),
                line_shape="hv",
            )
        )

    fig.update_layout(
        title="Monthly HRP target weight for the six largest average holdings",
        xaxis_title="Date",
        yaxis_title="Weight",
        height=450,
    )
    show_plotly_with_alt(
        fig,
        "Step lines of the six largest average HRP target weights against date, each holding flat between month-end decisions and jumping when the selected set changes.",
    )
else:
    print("No HRP weights available")
```

The step lines show submitted targets only. Actual portfolio weights drift with asset returns
between month-end decisions; neither the vectorized path nor the execution engine resets them
during the month.

## 12. Reading the Result

The comparison above is the notebook's evidence. What follows reads it rather than restating
what HRP is supposed to achieve, and it reads two columns rather than one: the Sharpe ratio the
allocators are ranked on, and the turnover each of them needed to get there.

```python
ranked = metrics_df.sort_values("Sharpe Ratio", ascending=False).reset_index(drop=True)
ranked.index = pd.RangeIndex(1, len(ranked) + 1, name="Rank by Sharpe")

print(f"Assets in the teaching universe: {n_assets}")
print(f"Estimation window each allocator sizes from: {LOOKBACK} sessions")
print(f"Assets selected each month: {TOP_N}")
ranked[["Method", "Sharpe Ratio", "Annual Vol", "Avg Turnover"]].round(4)
```

Read the ranking against what each allocator had to estimate, because that is the axis the four
differ on. Equal weight estimates nothing. Inverse volatility estimates one variance per asset.
HRP reads the correlations to build the tree and order the assets, and each split then compares
its two halves on their cluster variances, which `cluster_variance` computes as
$w^\top \Sigma w$ over the block - so the off-diagonal entries reach the weights. Minimum
variance with Ledoit-Wolf shrinkage reads and inverts the whole matrix.

Whether estimating more pays depends on how good the estimate is, and the ratio that decides
that is not favourable to the elaborate methods here. The argument for HRP is that inverting a
noisy covariance matrix amplifies estimation error, and it has force when the estimate is badly
under-determined - when the number of assets approaches or exceeds the number of observations.
This comparison is the opposite case: five selected assets estimated over 252 daily
observations, where the sample covariance is well conditioned and its inverse is not dominated
by noise. The regime HRP was designed for is not the regime tested here, so wherever it lands
in this table, the table is not evidence about that regime.

What HRP never does is invert the matrix, and that is the property it is chosen for: it returns
weights whatever the conditioning, where inversion has no unique solution at all once the
assets outnumber the observations.

The turnover column carries a second reading, and it is the one that survives a different
sample. Clustering is often described as producing more stable allocations, and the column is
where that claim would be tested - but it measures two things at once. Every row pays for the
monthly re-selection of five names out of fifteen, and every row then pays for whatever its
own sizing rule does with the names it keeps. The equal-weight row shows what the schedule
costs an allocator that makes no sizing decision at all, which is a reference point and not a
quantity the other three sit above: what a replacement costs depends on the weights being
replaced, so a concentrated allocator's re-selection cost is a different number rather than
the equal-weight one plus a margin.

### What is true regardless of the ranking

These are properties of the algorithm, and they hold wherever it lands in a given sample:

- It never inverts a covariance matrix, so it returns weights for any $N$ and $T$, including
  $N > T$ where minimum-variance optimization has no unique solution at all.
- It needs no expected-return forecast. Every input is second-moment.
- Its weights are determined by the correlation hierarchy and the variances, so they can be
  traced back to a specific split, as section 7 traced them.

### What this notebook cannot tell you

One universe, one sample, one estimation window, and a selection step evaluated on validation
data. A ranking under those conditions is a fact about this run. The regime where HRP was
designed to help - many assets relative to observations - is not the regime tested here, so this
result is not evidence against HRP in that regime either.

### How often an allocator had to fall back

When the allocation function raises (e.g., singular covariance), the
walk-forward backtest falls back to equal weights and increments a counter.

```python
print(f"Assets in universe: {n_assets}")
print(
    f"Invested backtest period: {portfolio_returns.index[0].date()} "
    f"to {portfolio_returns.index[-1].date()}"
)
print(f"Asset selection: ETF GBM walk-forward predictions (hash {ETF_GBM_PRED_HASH})")
if fallback_count:
    print("Allocator fallbacks (singular/invalid covariance):")
    for method_name, count in sorted(fallback_count.items()):
        print(f"  {method_name}: {count}")
else:
    print("Allocator fallbacks: none triggered")
```

## Key Takeaways

1. **HRP replaces inversion with ordering and recursive splits.** It still depends on estimated
   correlations and variances, so it reduces the exposure to estimation error rather than
   removing it.
2. **Avoiding inversion pays off when the estimate is under-determined.** With five assets and
   252 observations it is not, so nothing in this comparison exercises the property HRP is
   chosen for. The case for it is strongest when the number of assets approaches or exceeds the
   sample length, and that regime is not tested here.
3. **HRP concentrates too, on the low-variance assets.** Section 7 traces a majority weight in a
   single bond fund to two successive inverse-variance splits. Risk-based is not the same thing
   as diversified.
4. **Turnover measures re-selection and sizing together, and no column separates them.** Every
   allocator pays for replacing names the monthly signal drops, and then pays again for
   whatever its own sizing rule does with the names it keeps. Equal weight shows what the
   schedule costs an allocator that makes no sizing decision, which is a reference point and
   not a quantity the others sit above: what a replacement costs depends on the weights being
   replaced. So HRP's turnover against equal weight's is suggestive about stability and is not
   a measurement of it.
5. **The comparison is conditional in three ways.** A fixed 15-ETF ex-post universe, gross of
   costs in the vectorized path, and a selection signal ranked on validation data. It shows
   allocation behavior, not an out-of-sample estimate.
6. **Timing and costs need separate evidence.** Targets use returns through each decision date
   and become effective on the next bar; the execution bridge shows what commissions and
   slippage do to the vectorized result.

**Next**: [`07_conformal_position_sizing`](07_conformal_position_sizing.ipynb) sizes positions
from the width of a prediction interval rather than from a covariance matrix.
![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)
![notebook output](figures/p1_6.png)
![notebook output](figures/p1_7.png)
![notebook output](figures/p1_8.png)
![notebook output](figures/p1_9.png)
![notebook output](figures/p1_10.png)
![notebook output](figures/p1_11.png)
![notebook output](figures/p1_12.png)
![notebook output](figures/p1_13.png)
![notebook output](figures/p1_14.png)
![notebook output](figures/p1_15.png)