Dimensionamento de posições por Kelly: otimalidade do crescimento, alavancagem e risco de estimação

# 04_kelly_criterion.py


```py
# ---
# jupyter:
#   jupytext:
#     cell_metadata_filter: tags,-all
#     text_representation:
#       extension: .py
#       format_name: percent
#       format_version: '1.3'
#       jupytext_version: 1.19.3
#   kernelspec:
#     display_name: Python 3 (ipykernel)
#     language: python
#     name: python3
# ---

# %% [markdown]
# # How much to bet, and why the answer is unusable as given
#
# **Docker image**: `ml4t`
#
# ## Purpose
# The previous notebooks ask which assets to hold. This one asks how much of the account to put at
# risk, which is a separate question with a clean answer: the fraction that maximizes the expected
# logarithm of terminal wealth. Kelly derived it for a coin toss with known odds, and it extends to
# continuous returns and to portfolios.
#
# The derivation is exact and the answer, applied to estimated market inputs, is unusable as it
# comes out. Following it through from the coin toss to a multi-asset portfolio shows why, and the
# reason is not that the formula is wrong: it optimizes growth and has no term for ruin, no notion
# of a borrowing limit, and it scales directly with an expected-return estimate that nobody can
# make accurately.
#
# ## Learning objectives
#
# - Derive the optimal bet fraction for a binary wager, and show by simulation what betting more or
#   less than it does to the distribution of outcomes.
# - Extend the result to continuous returns, and say which approximation is being made and when it
#   holds.
# - Compute the multi-asset solution, read the leverage it implies, and work out what adverse move
#   would end the account at that leverage.
# - Recompute the same fraction on rolling windows and judge whether an estimate that unstable can
#   be acted on.
#
# ## Book reference
# Chapter 17, Section 17.4 (Defining baseline allocators).
#
# ## Prerequisites
#
# - `02_mean_variance_optimization`, for covariance estimation and its instability.
# - Basic probability: expectation, variance, and the log of a product.

# %% [markdown]
# ## Imports & Settings

# %%
"""Derive and apply Kelly position sizing under uncertainty."""

from collections.abc import Callable

import numpy as np
import plotly.graph_objects as go
import polars as pl
import sympy
from IPython.display import Markdown, display

# Portfolio analysis
from ml4t.diagnostic.evaluation import PortfolioAnalysis
from plotly.subplots import make_subplots
from scipy.optimize import minimize_scalar
from scipy.stats import binom
from sklearn.covariance import LedoitWolf
from sympy import diff, log, pprint, series, solve, symbols

from data import load_etfs
from utils.reproducibility import set_global_seeds
from utils.style import COLORS, ml4t_palette, show_plotly_with_alt

# %% tags=["parameters"]
# Production defaults; Papermill overrides for CI testing
SEED = 42
RISK_FREE_RATE = 0.02
TRAIN_END = "2019-12-31"
ROLLING_WINDOW_YEARS = 5

# %%
set_global_seeds(SEED)

# %% [markdown]
# ## Part 1: Kelly Criterion for Binary Outcomes
#
# Kelly began by analyzing games with binary outcomes (like a coin toss) with constant
# win-loss probability. A win returns the invested capital plus a payoff defined by the odds.
#
# **Key variables:**
# - $W_n$: Wealth after n bets
# - $b$: Odds (amount won per unit staked)
# - $p$: Probability of winning
# - $f$: Fraction of current wealth to risk

# %% [markdown]
# ### Derivation of the Kelly Formula
#
# After $n$ trials with $W$ wins and $L = n - W$ losses:
#
# $$W_n = W_0 (1 + bf)^W (1 - f)^L$$
#
# The Kelly criterion maximizes the expected growth rate:
#
# $$g(f) = p \log(1 + bf) + (1-p) \log(1 - f)$$
#
# Taking the derivative and setting to zero yields the optimal fraction:
#
# $$f^* = \frac{bp - (1-p)}{b} = \frac{\text{edge}}{\text{odds}}$$

# %%
# Symbolic derivation using SymPy
f, b, p = symbols("f b p", positive=True, real=True)

# Expected growth rate
growth_rate = p * log(1 + b * f) + (1 - p) * log(1 - f)

print("Growth rate function:")
pprint(growth_rate)

# %%
# First derivative (first-order condition)
first_deriv = diff(growth_rate, f)
print("\nFirst derivative:")
pprint(first_deriv)

# %%
# Solve for optimal f
f_star = solve(first_deriv, f)
print("\nOptimal Kelly fraction:")
pprint(f_star[0])

# Simplified form: f* = (bp + p - 1) / b = (bp - q) / b where q = 1 - p
print("\nSimplified: f* = (edge) / (odds) = (bp - q) / b")

# %%
# Verify second-order condition (must be negative for maximum)
second_deriv = diff(first_deriv, f)
print("Second derivative:")
pprint(second_deriv)

# Evaluate at a specific point to verify it's negative
soc_value = second_deriv.subs([(b, 1), (p, 0.6), (f, 0.2)])
print(f"\nSOC at b=1, p=0.6, f=0.2: {float(soc_value):.4f} (negative = maximum)")

# %% [markdown]
# ### Kelly Fraction Examples


# %%
def compute_kelly_fraction(win_prob: float, odds: float = 1.0) -> float:
    """Compute the Kelly fraction for given win probability and odds."""
    edge = odds * win_prob - (1 - win_prob)
    return edge / odds if edge > 0 else 0.0


# %% [markdown]
# The expected growth rate at one bet size: what the wager compounds at per bet in the long run.


# %%
def compute_growth_rate(win_prob: float, fraction: float, odds: float = 1.0) -> float:
    """Compute expected growth rate for given parameters."""
    if fraction <= 0 or fraction >= 1:
        return -np.inf
    return win_prob * np.log(1 + odds * fraction) + (1 - win_prob) * np.log(1 - fraction)


# Example calculations
examples = [
    (0.55, 1.0),  # Slight edge, even odds
    (0.60, 1.0),  # Good edge, even odds
    (0.55, 2.0),  # Slight edge, 2:1 odds
    (0.70, 0.5),  # Strong edge, unfavorable odds
]

pl.DataFrame(
    [
        {
            "win probability": prob,
            "odds won per unit staked": odds,
            "Kelly fraction": compute_kelly_fraction(prob, odds),
            "growth rate per bet": compute_growth_rate(
                prob, compute_kelly_fraction(prob, odds), odds
            ),
        }
        for prob, odds in examples
    ]
)

# %% [markdown]
# ### Growth Rate as Function of Bet Size

# %%
# Compute growth rates for different probabilities and bet fractions
fractions = np.linspace(0.01, 0.99, 100)
probabilities = [0.55, 0.60, 0.65, 0.70, 0.75]

fig = go.Figure()

colors = ml4t_palette(len(probabilities), categorical=True)

for i, prob in enumerate(probabilities):
    kelly = compute_kelly_fraction(prob)
    growth_rates = [compute_growth_rate(prob, f) for f in fractions]

    fig.add_scatter(
        x=fractions,
        y=growth_rates,
        mode="lines",
        name=f"P={prob:.0%}",
        line=dict(color=colors[i % len(colors)], width=2),
    )

    # Mark optimal Kelly fraction
    optimal_growth = compute_growth_rate(prob, kelly)
    fig.add_scatter(
        x=[kelly],
        y=[optimal_growth],
        mode="markers",
        marker=dict(size=10, color=colors[i % len(colors)], symbol="diamond"),
        showlegend=False,
        hovertemplate=f"Kelly={kelly:.1%}<br>Growth={optimal_growth:.4f}",
    )

# %%
fig.update_layout(
    title="Expected growth rate against bet fraction, by win probability",
    xaxis_title="Bet Fraction (f)",
    yaxis_title="Expected Growth Rate",
    xaxis_tickformat=".0%",
    height=500,
    legend=dict(yanchor="top", y=0.99, xanchor="right", x=0.99),
)
show_plotly_with_alt(
    fig,
    "Five curves of expected growth rate against bet fraction, one per win probability, each rising to an interior peak marked with a diamond and then falling away steeply towards a bet fraction of one.",
)

# %% [markdown]
# ### Simulating Wealth Paths


# %%
def simulate_wealth_paths(
    outcomes: np.ndarray,
    fraction: float,
    odds: float = 1.0,
    initial_wealth: float = 100.0,
) -> np.ndarray:
    """Simulate wealth paths for Kelly betting."""
    # Compute growth factors
    growth_factors = np.where(
        outcomes,
        1 + odds * fraction,  # Win
        1 - fraction,  # Loss
    )

    # Cumulative wealth
    wealth = initial_wealth * np.cumprod(growth_factors, axis=0)
    return wealth


# %% [markdown]
# Every panel below is driven by one shared sequence of coin-toss outcomes, so the panels
# differ only in the stake and the comparison between them is paired rather than a comparison
# of two draws. What is plotted is the spread across the simulated paths: the median path, the
# band containing the middle half of them and the band containing the middle nine tenths.


# %%
def add_percentile_bands(
    fig: go.Figure, trials: np.ndarray, wealth_pcts: np.ndarray, col_idx: int, showlegend: bool
) -> None:
    """Add 5-95 and 25-75 percentile bands with median path."""
    bands = [
        (4, None),
        (0, "rgba(10, 22, 40, 0.2)"),
        (3, None),
        (1, "rgba(10, 22, 40, 0.4)"),
    ]
    for idx, fillcolor in bands:
        fig.add_scatter(
            x=trials,
            y=wealth_pcts[idx],
            mode="lines",
            line=dict(width=0),
            fill="tonexty" if fillcolor else None,
            fillcolor=fillcolor,
            showlegend=False,
            row=1,
            col=col_idx,
        )
    fig.add_scatter(
        x=trials,
        y=wealth_pcts[2],
        mode="lines",
        line=dict(color=COLORS["blue"], width=2),
        name="Median",
        showlegend=showlegend,
        row=1,
        col=col_idx,
    )


# %%
# Simulate paths for different Kelly multiples
n_trials = 1000  # bets in one path
n_simulations = 500  # independent paths, which is what the percentile bands are taken across
win_prob = 0.55
odds = 1.0
INITIAL_WEALTH = 100.0
kelly = compute_kelly_fraction(win_prob, odds)
kelly_multiples = [0.25, 0.5, 1.0, 1.5, 2.0]
common_outcomes = np.random.default_rng(SEED).random((n_trials, n_simulations)) < win_prob

print(f"{n_simulations:,} independent paths of {n_trials:,} bets each, one shared outcome draw.")

fig = make_subplots(
    rows=1,
    cols=len(kelly_multiples),
    subplot_titles=[f"{mult:.0%} Kelly" for mult in kelly_multiples],
    shared_yaxes=True,
)

for i, mult in enumerate(kelly_multiples):
    fraction = kelly * mult
    wealth = simulate_wealth_paths(common_outcomes, fraction, odds, INITIAL_WEALTH)
    wealth_pcts = np.percentile(wealth, [5, 25, 50, 75, 95], axis=1)
    add_percentile_bands(fig, np.arange(n_trials), wealth_pcts, col_idx=i + 1, showlegend=(i == 0))

# %%
fig.update_xaxes(title_text="Trial")
fig.update_yaxes(type="log", title_text="Wealth")
fig.update_layout(
    title="Simulated wealth at five multiples of the Kelly fraction",
    height=400,
    showlegend=True,
)
show_plotly_with_alt(
    fig,
    "Five panels of simulated wealth on a logarithmic axis against trial number, at a quarter, a half, one, one and a half, and twice the Kelly fraction, the median path rising and the shaded percentile bands widening as the multiple grows.",
)

print(f"Kelly fraction at a win probability of {win_prob:.0%} and even odds: {kelly:.1%}")
print(f"Terminal wealth after {n_trials:,} bets, from {INITIAL_WEALTH:,.0f} of capital:")

terminal_rows = []
for mult in kelly_multiples:
    terminal = simulate_wealth_paths(common_outcomes, kelly * mult, odds, INITIAL_WEALTH)[-1]
    p05, p50, p95 = np.percentile(terminal, [5, 50, 95])
    terminal_rows.append(
        {
            "multiple of Kelly": f"{mult:.0%}",
            "5th percentile": p05,
            "median": p50,
            "95th percentile": p95,
        }
    )
pl.DataFrame(terminal_rows)

# %% [markdown]
# ### Distribution of Terminal Wealth


# %%
def terminal_wealth_distribution(
    n_trials: int,
    win_prob: float,
    fraction: float,
    odds: float = 1.0,
    initial_wealth: float = 100.0,
) -> pl.DataFrame:
    """Compute theoretical distribution of terminal wealth."""
    rv = binom(n=n_trials, p=win_prob)

    results = []
    for n_wins in range(n_trials + 1):
        n_losses = n_trials - n_wins
        log_wealth = (
            np.log(initial_wealth)
            + n_wins * np.log(1 + odds * fraction)
            + n_losses * np.log(1 - fraction)
        )
        prob = rv.pmf(n_wins)
        results.append(
            {
                "n_wins": n_wins,
                "log_wealth": log_wealth,
                "wealth": np.exp(log_wealth),
                "probability": prob,
            }
        )

    return pl.DataFrame(results)


# %%
# Compare terminal wealth distributions
n_trials = 100
win_prob = 0.55
kelly = compute_kelly_fraction(win_prob)

fig = go.Figure()

for mult in [0.5, 1.0, 1.5, 2.0]:
    dist = terminal_wealth_distribution(n_trials, win_prob, kelly * mult)

    # Compute expected log wealth
    expected_log_wealth = (dist["log_wealth"] * dist["probability"]).sum()

    fig.add_scatter(
        x=dist["log_wealth"].to_list(),
        y=dist["probability"].to_list(),
        mode="lines",
        name=f"{mult:.0%} Kelly (E[log W]={expected_log_wealth:.1f})",
        line=dict(width=2),
    )

fig.update_layout(
    title="Terminal log-wealth distribution, by multiple of the Kelly fraction",
    xaxis_title="Log(Wealth)",
    yaxis_title="Probability",
    height=450,
)
show_plotly_with_alt(
    fig,
    "Four terminal log-wealth distributions after a hundred bets, at a half, one, one and a half and twice Kelly, the full-Kelly curve centred highest and the larger multiples spread wider and further left.",
)

# %% [markdown]
# ## Part 2: Kelly for Continuous Returns (Single Asset)
#
# For continuous return distributions, we use a Taylor expansion of $\log(1+x)$:
#
# $$\log(1+x) \approx x - \frac{x^2}{2} + \frac{x^3}{3} - ...$$
#
# For small returns, the optimal Kelly fraction becomes:
#
# $$f^* \approx \frac{\mu - r_f}{\sigma^2}$$
#
# where $\mu$ is expected return, $r_f$ is risk-free rate, and $\sigma^2$ is variance.

# %% [markdown]
# ### Taylor Expansion Visualization


# %%
def taylor_polynomial(n_terms: int) -> Callable:
    """Return Taylor polynomial approximation of log(1+x)."""
    x = symbols("x")
    expansion = series(log(1 + x), x, x0=0, n=n_terms).removeO()
    return sympy.lambdify(x, expansion, "numpy")


# %%
# Compare Taylor approximations
x_vals = np.linspace(-0.5, 1.0, 200)
true_log = np.log(1 + x_vals)

fig = go.Figure()

fig.add_scatter(
    x=x_vals,
    y=true_log,
    mode="lines",
    name="log(1+x)",
    line=dict(width=3, color=COLORS["blue"]),
)

colors = ml4t_palette(4, categorical=True)
for i, n in enumerate([2, 3, 4, 5]):
    taylor = taylor_polynomial(n)
    approx = taylor(x_vals)
    fig.add_scatter(
        x=x_vals,
        y=approx,
        mode="lines",
        name=f"Taylor (n={n})",
        line=dict(width=2, dash="dash", color=colors[i]),
    )

fig.update_layout(
    title="log(1+x) and its Taylor polynomials of order two to five",
    xaxis_title="x",
    yaxis_title="y",
    yaxis_range=[-2, 1.5],
    height=450,
)
fig.add_vline(x=0, line_dash="dot", line_color=COLORS["neutral"])
fig.add_hline(y=0, line_dash="dot", line_color=COLORS["neutral"])
show_plotly_with_alt(
    fig,
    "The logarithm of one plus x against x, with four Taylor polynomials of increasing order overlaid, each tracking the curve near zero and departing from it further out.",
)

# %% [markdown]
# ### Kelly Fraction for Market Returns


# %%
def kelly_fraction_continuous(
    mean_return: float, std_return: float, risk_free: float = 0.0
) -> float:
    """Kelly fraction from the second-order expansion of log wealth."""
    excess_return = mean_return - risk_free
    return excess_return / (std_return**2)


# %% [markdown]
# The growth rate to maximize, evaluated on the realized return series rather than on an
# assumed distribution - so no normality approximation enters this version.


# %%
def empirical_growth_rate(returns: np.ndarray, fraction: float, risk_free: float = 0.0) -> float:
    """Annualized mean log growth over observed daily returns."""
    daily_risk_free = (1 + risk_free) ** (1 / 252) - 1
    wealth_multiples = 1 + daily_risk_free + fraction * (returns - daily_risk_free)
    if np.any(wealth_multiples <= 0):
        return -np.inf
    return float(np.log(wealth_multiples).mean() * 252)


# %% [markdown]
# A fraction large enough to make wealth negative on some historical return has no logarithm,
# so the search is bounded by the worst return in the sample.


# %%
def optimal_kelly_empirical(returns: np.ndarray, risk_free: float = 0.0) -> tuple[float, float]:
    """Optimize growth while keeping wealth positive for every observed return."""
    daily_risk_free = (1 + risk_free) ** (1 / 252) - 1
    excess_returns = returns - daily_risk_free
    worst_excess_return = float(excess_returns.min())
    if worst_excess_return >= 0:
        raise ValueError("Observed returns do not establish a finite leverage boundary")
    max_observed_safe = (1 + daily_risk_free) / -worst_excess_return
    result = minimize_scalar(
        lambda fraction: -empirical_growth_rate(returns, fraction, risk_free),
        bounds=(0.0, max_observed_safe * (1 - 1e-9)),
        method="bounded",
    )
    return float(result.x), max_observed_safe


# %%
# Load SPY as proxy for S&P 500 from canonical data
etf_data = load_etfs()
spy_data = etf_data.filter(pl.col("symbol") == "SPY").sort("timestamp")
sp500_returns = spy_data.select(
    "timestamp",
    pl.col("close").pct_change().alias("return"),
).drop_nulls()

# Compute annual statistics
annual_ret = float(sp500_returns["return"].mean() * 252)
annual_vol = float(sp500_returns["return"].std() * np.sqrt(252))

print(
    f"SPY statistics ({sp500_returns['timestamp'].min().year}-"
    f"{sp500_returns['timestamp'].max().year}):"
)
print(f"  Annual Return: {annual_ret:.2%}")
print(f"  Annual Volatility: {annual_vol:.2%}")
print(f"  Sharpe Ratio: {(annual_ret - RISK_FREE_RATE) / annual_vol:.2f}")

# Both estimators take the same annual risk-free rate
# so the analytical Taylor approximation and the numerical optimizer are
# compared like-for-like.
kelly_approx = kelly_fraction_continuous(annual_ret, annual_vol, RISK_FREE_RATE)
kelly_empirical, max_observed_safe = optimal_kelly_empirical(
    sp500_returns["return"].to_numpy(), RISK_FREE_RATE
)

print("\nKelly Fractions:")
print(f"  Analytical (Taylor approx): {kelly_approx:.1%}")
print(f"  Empirical log-growth optimum: {kelly_empirical:.1%}")
print(f"  Observed-return leverage boundary: {max_observed_safe:.1%}")

# %% [markdown]
# The empirical optimizer keeps wealth positive for every observed SPY return. That is a
# transparent historical-support constraint, not a guarantee against a worse future return.

# %% [markdown]
# ### Rolling Kelly Fraction
#
# The fraction above is one number computed from the whole history. Recomputing it inside a
# moving window asks the question a position sizer actually faces: what would this formula have
# told someone acting on it at each point in time.
#
# The window length is the trade-off. Kelly divides an expected return by a variance, and the
# expected return is the badly estimated half - its standard error falls with the square root of
# the window, so a short window produces a fraction that swings on noise. Five years is long
# enough that the estimate is not dominated by a single year's returns and short enough that the
# series still moves, which is what makes the instability visible rather than smoothed away.
# Shortening it widens the swings below; lengthening it flattens them without making the
# estimate any more actionable.

# %%
rolling_window = int(ROLLING_WINDOW_YEARS * 252)
# A window of W over N returns yields N - W + 1 estimates. The chart below is about how much
# the fraction moves, so it needs enough of them to move over: a year's worth is the floor.
rolling_estimates = sp500_returns.height - rolling_window + 1
if rolling_estimates < 252:
    raise ValueError(
        f"A {ROLLING_WINDOW_YEARS}-year window over {sp500_returns.height:,} returns leaves "
        f"{rolling_estimates} rolling estimates, fewer than the 252 this chart needs to show "
        "how the fraction moves. Shorten the window or extend the history."
    )

# Rolling mean and std
rolling_stats = sp500_returns.with_columns(
    [
        pl.col("return").rolling_mean(window_size=rolling_window).alias("rolling_mean"),
        pl.col("return").rolling_std(window_size=rolling_window).alias("rolling_std"),
    ]
).drop_nulls()

# Compute Kelly fraction
rolling_stats = rolling_stats.with_columns(
    [
        (
            (pl.col("rolling_mean") * 252 - RISK_FREE_RATE) / (pl.col("rolling_std") ** 2 * 252)
        ).alias("kelly"),
    ]
)

# %%
fig = make_subplots(
    rows=2,
    cols=1,
    shared_xaxes=True,
    vertical_spacing=0.1,
    subplot_titles=[
        f"Annualized return and volatility over a {ROLLING_WINDOW_YEARS}-year window",
        "Kelly fraction those two moments imply",
    ],
)

# Rolling return and volatility
fig.add_scatter(
    x=rolling_stats["timestamp"].to_list(),
    y=(rolling_stats["rolling_mean"] * 252).to_list(),
    name="Annual Return",
    line=dict(color=COLORS["blue"]),
    row=1,
    col=1,
)
fig.add_scatter(
    x=rolling_stats["timestamp"].to_list(),
    y=(rolling_stats["rolling_std"] * np.sqrt(252)).to_list(),
    name="Annual Volatility",
    line=dict(color=COLORS["amber"]),
    row=1,
    col=1,
)

# Kelly fraction
kelly_values = rolling_stats["kelly"].to_list()
fig.add_scatter(
    x=rolling_stats["timestamp"].to_list(),
    y=kelly_values,
    name="Kelly Fraction",
    line=dict(color=COLORS["copper"]),
    row=2,
    col=1,
)
_ = fig.add_hline(y=1.0, line_dash="dash", line_color=COLORS["neutral"], row=2, col=1)

# %%
fig.update_xaxes(title_text="Date", row=2, col=1)
fig.update_yaxes(title_text="Annualized value", tickformat=".0%", row=1, col=1)
fig.update_yaxes(title_text="Kelly fraction", tickformat=".0%", row=2, col=1)
fig.update_layout(
    height=600,
    title="Rolling annualized return and volatility, and the Kelly fraction they imply",
)
show_plotly_with_alt(
    fig,
    "Two stacked panels against date: rolling annualized return and volatility on top, and the Kelly fraction they imply below, which swings across a very wide range with a reference line at one.",
)

print("\nKelly Fraction Statistics:")
print(f"  Mean: {np.mean(kelly_values):.1%}")
print(f"  Std:  {np.std(kelly_values):.1%}")
print(f"  Min:  {np.min(kelly_values):.1%}")
print(f"  Max:  {np.max(kelly_values):.1%}")

# %% [markdown]
# ## Part 3: Kelly for Multiple Assets
#
# For a portfolio of $n$ assets, the optimal Kelly allocation is:
#
# $$\mathbf{f}^* = \Sigma^{-1} \boldsymbol{\mu}$$
#
# where $\Sigma^{-1}$ is the **precision matrix** (inverse covariance) and $\boldsymbol{\mu}$
# is the vector of expected excess returns.
#
# The direction of this vector is the same as the maximum-Sharpe portfolio's under aligned
# assumptions, and the difference is what happens to its length. A mean-variance solution rescales
# its weights to sum to one, which throws the length away and leaves only the mix. Kelly keeps it:
# the magnitude of the vector *is* the leverage, and nothing in the formula bounds it.

# %% [markdown]
# ### Load Multi-Asset Data

# %%
# Load data for a set of ETFs from canonical data
SYMBOLS = ["SPY", "QQQ", "IWM", "EFA", "EEM", "TLT", "GLD", "VNQ"]
START_DATE = "2010-01-01"
END_DATE = "2024-12-01"

# Load from canonical ETF data
multi_etf = 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())
)
prices = (
    multi_etf.select(["timestamp", "symbol", "close"])
    .pivot(on="symbol", index="timestamp", values="close")
    .sort("timestamp")
    .fill_null(strategy="forward")
    .drop_nulls()
)
assets = [column for column in prices.columns if column != "timestamp"]
returns = prices.select(
    "timestamp",
    *[pl.col(symbol).pct_change().alias(symbol) for symbol in assets],
).drop_nulls()
train_returns = returns.filter(pl.col("timestamp") <= pl.lit(TRAIN_END).str.to_datetime())
test_returns = returns.filter(pl.col("timestamp") > pl.lit(TRAIN_END).str.to_datetime())

# %%
if train_returns.is_empty() or test_returns.is_empty():
    raise ValueError("Both train and test windows must contain returns")
if train_returns["timestamp"].max() >= test_returns["timestamp"].min():
    raise ValueError("Train and test windows overlap")

print(
    f"Training: {train_returns['timestamp'].min()} to "
    f"{train_returns['timestamp'].max()} ({train_returns.height:,} returns)"
)
print(
    f"Test: {test_returns['timestamp'].min()} to "
    f"{test_returns['timestamp'].max()} ({test_returns.height:,} returns)"
)

# %%
print(f"Assets with complete history: {len(assets)}")

# %% [markdown]
# This fixed eight-ETF set uses current-vintage histories. It demonstrates allocation mechanics;
# it is not a point-in-time, survivorship-free universe study.

# %%
# Compute statistics
train_matrix = train_returns.select(assets).to_numpy()
annual_total_returns = train_matrix.mean(axis=0) * 252
annual_excess_returns = annual_total_returns - RISK_FREE_RATE
annual_cov = np.cov(train_matrix, rowvar=False, ddof=1) * 252

training_moments = pl.DataFrame(
    {
        "symbol": assets,
        "annual_return": annual_total_returns,
        "annual_excess_return": annual_excess_returns,
    }
)
print(f"Training covariance condition number: {np.linalg.cond(annual_cov):.1f}")
training_moments

# %% [markdown]
# ### Kelly Portfolio Allocation

# %%
# Kelly allocation
kelly_allocation = np.linalg.solve(annual_cov, annual_excess_returns)

print(f"Raw Kelly gross leverage: {np.abs(kelly_allocation).sum():.1%}")
print(f"Raw Kelly net exposure: {kelly_allocation.sum():.1%}")

# %%
fig = go.Figure()

fig.add_bar(
    x=assets,
    y=kelly_allocation,
    name="Kelly (Raw)",
    marker_color=COLORS["blue"],
)

# Add reference line for equal weight
fig.add_hline(
    y=1 / len(assets),
    line_dash="dash",
    line_color=COLORS["neutral"],
    annotation_text=f"Equal Weight: {1 / len(assets):.1%}",
    annotation_position="top left",
)

fig.update_layout(
    title="Raw Kelly allocation per ETF against the equal-weight reference",
    xaxis_title="Asset",
    yaxis_title="Allocation (Can Exceed 100%)",
    yaxis_tickformat=".0%",
    height=400,
)
show_plotly_with_alt(
    fig,
    "Bars of the raw Kelly allocation per ETF, several of them far above the equal-weight reference line and some below zero.",
)

# %% [markdown]
# ### Fractional Kelly for Risk Management
#
# Full Kelly is often too aggressive in practice due to:
# - Estimation error in $\mu$ and $\Sigma$
# - Non-normal return distributions (fat tails)
# - Borrowing constraints
#
# **Half-Kelly** ($f^*/2$) is common in practice, sacrificing some growth for stability.

# %% [markdown]
# The fractional weights below are the raw Kelly solution scaled, with no normalization applied
# afterwards, so each multiple carries genuinely different leverage rather than the same book
# rescaled. Equal weight is included at a gross exposure of exactly one, which is the reference
# every number in the table should be read against.

# %%
kelly_multiples = [0.25, 0.5, 1.0]
test_matrix = test_returns.select(assets).to_numpy()

portfolio_returns = {}
portfolio_gross = {}
for mult in kelly_multiples:
    weights = kelly_allocation * mult
    gross = float(np.abs(weights).sum())

    pf_returns = test_matrix @ weights
    portfolio_returns[f"{mult:.0%} Kelly"] = pf_returns
    portfolio_gross[f"{mult:.0%} Kelly"] = gross

# Add equal weight for comparison (gross = 1.0 by construction)
equal_weights = np.full(len(assets), 1 / len(assets))
portfolio_returns["Equal Weight"] = test_matrix @ equal_weights
portfolio_gross["Equal Weight"] = float(np.abs(equal_weights).sum())

# %% [markdown]
# **Gross exposure** is the sum of the absolute position sizes: a book three times the account
# long and twice it short has five times the account at risk and is one account net long - the
# same gross as a book two and a half times long against the same short, which is the one that
# is flat on net. Its reciprocal is the move that would end
# the account, and that is a worst case rather than a market move - it needs every long to fall
# and every short to rise by that percentage on the same day. Net exposure is far smaller than
# gross, so a uniform market decline of that size does not do it. What the number establishes is
# that at high leverage the worst case sits inside a single ordinary session.

# %%
pl.DataFrame(
    [
        {
            "allocation": name,
            "gross exposure (x account)": gross,
            "adverse move against every leg that empties the account": 1.0 / gross,
        }
        for name, gross in portfolio_gross.items()
    ]
)

# %% [markdown]
# ### What each Kelly multiple did on the test window

# %%
dates = test_returns["timestamp"]
comparison_metrics = []

for name, pf_ret in portfolio_returns.items():
    # The lowest the cumulative wealth path reaches, not the worst single period: the
    # section reads the growth chart off this number, and a chart is a path.
    minimum_wealth_multiple = float(np.min(np.cumprod(1 + pf_ret)))
    if minimum_wealth_multiple <= 0:
        raise ValueError(f"{name} reaches nonpositive wealth in the test window")
    pa = PortfolioAnalysis(
        returns=pl.Series("returns", pf_ret),
        dates=dates,
        risk_free=RISK_FREE_RATE,
        periods_per_year=252,
    )

    metrics = pa.compute_summary_stats()
    comparison_metrics.append(
        {
            "strategy": name,
            "annual_return": metrics.annual_return,
            "annual_vol": metrics.annual_volatility,
            "sharpe": metrics.sharpe_ratio,
            "max_dd": metrics.max_drawdown,
            "gross_exposure": portfolio_gross[name],
            "min_wealth_multiple": minimum_wealth_multiple,
        }
    )

metrics_df = pl.DataFrame(comparison_metrics).sort("gross_exposure", descending=True)
metrics_df

# %% [markdown]
# The rows are ordered by gross exposure rather than by Sharpe ratio, because sorting this table
# by Sharpe puts the allocation that lost the account at the top of it. That is not a quirk of the
# sort: the Sharpe ratio divides a mean daily return by a daily standard deviation, and neither
# term knows that the wealth path they were computed from went to nearly zero and could not have
# been held. Read the `min_wealth_multiple` column, which is the lowest point the compounded path
# reached, against the ratio beside it.
#
# Every path stays strictly above zero wealth over these observations, which is what makes the
# logarithms and the drawdowns well defined. It is a statement about arithmetic, not about
# survival, and it says nothing about a return worse than any in this window.

# %%
# Cumulative returns comparison
fig = go.Figure()

strategy_colors = dict(
    zip(
        portfolio_returns,
        ml4t_palette(len(portfolio_returns), categorical=True),
        strict=True,
    )
)

for name, pf_ret in portfolio_returns.items():
    cum_ret = (1 + pf_ret).cumprod()
    fig.add_scatter(
        x=dates,
        y=cum_ret,
        mode="lines",
        name=name,
        line=dict(color=strategy_colors[name], width=2),
    )

fig.update_layout(
    title="Growth of one dollar on the test window, by allocation",
    xaxis_title="Date",
    yaxis_title="Growth of $1",
    height=500,
    legend=dict(yanchor="top", y=0.99, xanchor="left", x=0.01),
)
show_plotly_with_alt(
    fig,
    "Four growth-of-one-dollar paths on the test window at a quarter, a half and one times Kelly, and at equal weight, the leveraged paths swinging across orders of magnitude while equal weight stays near one.",
)

# %%
fig = go.Figure()
for row in metrics_df.iter_rows(named=True):
    fig.add_scatter(
        x=[row["annual_vol"]],
        y=[row["annual_return"]],
        mode="markers+text",
        name=row["strategy"],
        text=[row["strategy"]],
        textposition="middle right" if row["strategy"] == "Equal Weight" else "top center",
        marker=dict(
            size=max(10, abs(row["sharpe"]) * 15),
            color=strategy_colors[row["strategy"]],
        ),
    )
fig.update_layout(
    title="Annualized return against annualized volatility, by allocation",
    xaxis_title="Annualized volatility",
    yaxis_title="Annualized return",
    xaxis_tickformat=".0%",
    yaxis_tickformat=".0%",
    height=450,
)
show_plotly_with_alt(
    fig,
    "Scatter of the four allocations, annualized volatility against annualized return, marker size scaled by absolute Sharpe ratio, the full-Kelly point far out on the volatility axis and below zero on the return axis.",
)

# %% [markdown]
# ## Part 4: What the derivation assumed
#
# Two different kinds of result sit in this notebook, and only the first is exact. The binary bet's
# $f^* = (bp - q)/b$ is an exact maximizer of expected log wealth for that gamble. The
# mean-over-variance form used for assets is not: it comes from expanding $\log(1 + x)$ and keeping
# terms to the square, so it is an approximation whose accuracy depends on how large the leveraged
# returns entering that logarithm are - not on whether the returns are close to normal. Each result
# also assumed something the market does not supply. Naming those in order separates the parts that
# survive contact with data from the parts that do not.
#
# **The odds are known.** In the binary case $p$ and $b$ are given. In the market case both are
# estimated, the solution divides by the estimated variance and scales with the estimated mean,
# and the rolling section above shows how far that estimate moves on one instrument.
#
# **The bet repeats many times, and the stake is a fraction of current wealth.** Kelly is a
# long-run growth argument: it wins in the limit and says nothing about the horizon a particular
# investor has. Betting a fraction of *current* wealth is what makes ruin impossible in the coin
# toss, because a fraction of a positive number stays positive. That property is lost the moment
# leverage exceeds one, since a position larger than the account can lose more than the account.
#
# **The leveraged returns are small enough for the second-order expansion to hold.** The
# mean-over-variance form drops every term past the square, and what those dropped terms are worth
# grows with the size of the returns fed into the logarithm - which leverage is precisely what
# magnifies. Near-normality does not rescue it: fat tails make the discarded terms matter more, but
# even a well-behaved distribution breaks the approximation once positions are large enough.
#
# **Borrowing is free and unlimited.** There is no funding rate in the objective, no cap on gross
# exposure and no liquidation rule. Part 3 computes what that assumption is worth in this sample.

# %%
# Example: Shrinkage estimator for more stable Kelly allocation
lw = LedoitWolf()
lw.fit(train_matrix)
shrunk_cov = lw.covariance_ * 252

kelly_shrunk = np.linalg.solve(shrunk_cov, annual_excess_returns)

print(f"Total leverage (raw):    {np.abs(kelly_allocation).sum():.1%}")
print(f"Total leverage (shrunk): {np.abs(kelly_shrunk).sum():.1%}")

comparison = pl.DataFrame(
    {
        "symbol": assets,
        "kelly_raw": kelly_allocation,
        "kelly_shrunk": kelly_shrunk,
        "difference": kelly_allocation - kelly_shrunk,
    }
).with_columns(pl.col(pl.Float64).round(3))
comparison

# %% [markdown]
# ## What the test window actually says

# %% [markdown]
# Every number in the paragraph below is read from the table above it rather than typed, so it
# moves with the run instead of describing an earlier one.

# %%
rows = {row["strategy"]: row for row in metrics_df.iter_rows(named=True)}
full, half, quarter = rows["100% Kelly"], rows["50% Kelly"], rows["25% Kelly"]
ruin_move_full = 1.0 / full["gross_exposure"]

display(
    Markdown(
        "The lowest point of the full-Kelly wealth path is what decides how to read the growth "
        f"chart above. At its worst the path is down to {full['min_wealth_multiple']:.1e} times "
        f"the capital it started with, while carrying {full['gross_exposure']:.0f} times that "
        "capital in gross positions - which is the same thing the `max_dd` column says when it "
        f"reads {full['max_dd']:.3f}. The chart's vertical axis is logarithmic, so a fall of that "
        "size still looks like a line on the page; the number is what says the account is gone. "
        "No broker holds a position through it, and the arithmetic above says why: at that gross "
        f"exposure it takes {ruin_move_full:.2%} against every leg at once - every long down that "
        "much and every short up that much - to empty the account. That is a worst case rather "
        "than a market move, since net exposure is a fraction of gross and a uniform decline does "
        f"not do it; what matters is that {full['gross_exposure']:.0f} times leverage puts the "
        "worst case inside a single ordinary session.\n\n"
        f"Halving the fraction does not rescue it: half Kelly bottoms at "
        f"{half['min_wealth_multiple']:.2f} times its starting capital, with a drawdown of "
        f"{abs(half['max_dd']):.1%}, and only at a quarter does the worst point stay above "
        f"{quarter['min_wealth_multiple']:.2f} times. What separates them is not the sign of the "
        "estimates but how much leverage is taken on the strength of them.\n\n"
        "So the curves are not a track record. They are what the formula's answer would have "
        "produced given borrowing that is unlimited, free of interest and never called, and the "
        "value of computing them is precisely that the assumption is visible in the leverage "
        "number rather than hidden."
    )
)

# %% [markdown]
# ## Key takeaways
#
# 1. **Kelly answers a question about growth, not about survival.** Maximizing the expected log of
#    terminal wealth is a defensible objective and it says nothing about the path, so the solution
#    is happy to accept a route through near-total loss on the way to a higher expectation.
# 2. **The formula has no notion of a budget.** Dividing expected returns by a covariance matrix
#    produces a number, and nothing in it is bounded by the capital available. Any use of it in
#    practice is the formula plus a leverage constraint, and the constraint is doing at least as
#    much work as the formula.
# 3. **Fractional Kelly is a cap, not a fix.** Scaling the solution scales the leverage with it,
#    so a quarter of an unimplementable book is a quarter as unimplementable and no more. The
#    reason to use a fraction is that the inputs are estimated and the solution scales with the
#    error in them; the reason it is not sufficient is that a fraction of an unbounded number is
#    still unbounded.
# 4. **The single-asset rolling estimate is the honest picture of the input problem.** The same
#    formula on the same instrument over rolling five-year windows swings across a range no
#    position sizer could act on. That instability is the estimate, not the market.
# 5. **Check the wealth path, not only the return series.** A leveraged return series can produce a
#    respectable annualized figure while passing through a wealth multiple that would have ended
#    the account. The minimum wealth multiple is one line of code and it is what makes the
#    difference visible.
#
# ### Known limitations
#
# - Borrowing is free, unlimited and never margin-called throughout. Introducing any of a funding
#   rate, a leverage cap, or a liquidation rule changes every result in Part 3.
# - Inputs are estimated once on the training window and frozen. A rolling re-estimate would change
#   the position sizes continuously, and the rolling-Kelly section shows by how much.
# - The mean-over-variance form keeps the expansion of the logarithm only to the squared term. It
#   is accurate while the leveraged returns entering that logarithm are small, and leverage is what
#   makes them large, so the approximation is weakest exactly where it is being used hardest.
# - The universe is a small set of funds selected because they exist today.
#
# **Next:** [`05_factor_allocation_evidence`](05_factor_allocation_evidence.ipynb) asks whether
# asset characteristics explain later returns at all. Section 17.4 covers Kelly sizing among the
# baseline allocators.

```