docs: score frontier portfolios with jquantstats

book/marimo/notebooks/cla.py

@@ -4,6 +4,8 @@
 #     "marimo==0.14.13",
 #     "numpy==2.3.0",
 #     "plotly==6.7.0",
+#     "polars==1.44.1",
+#     "jquantstats==0.11.0",
 #     "cvx-linalg>=0.9.3",
 #     "cvxcla"
 # ]
@@ -22,9 +24,15 @@
 with app.setup:
     import marimo as mo
     import numpy as np
+    import polars as pl
+    from jquantstats import Data
 
     from cvxcla import CLA
 
+    # Trading days in a year, and the length of the simulated history.
+    PERIODS = 260
+    HISTORY = 4 * PERIODS
+
 
 @app.cell
 def _():
@@ -33,6 +41,13 @@ def _():
     # The Critical Line Algorithm
     We compute an efficient frontier using the critical line algorithm (cla).
     The method was introduced by Harry M Markowitz in 1956.
+
+    Rather than invent a mean vector and a covariance matrix out of thin air, we
+    simulate a return history, estimate both from it, and trace the frontier of
+    the estimated problem. That way every portfolio on the frontier has a
+    *realised* return series too, which we hand to
+    [jQuantStats](https://github.com/jebel-quant/jquantstats) at the bottom of
+    this notebook.
     """
     )
     return
@@ -47,40 +62,149 @@ def _():
 
 
 @app.function(hide_code=True)
-def cla(n):
+def business_dates(periods):
+    """Build a Monday-to-Friday date index of the given length."""
+    calendar = pl.date_range(pl.date(2015, 1, 1), pl.date(2035, 1, 1), interval="1d", eager=True)
+    return calendar.filter(calendar.dt.weekday() <= 5).head(periods)
+
+
+@app.function(hide_code=True)
+def simulate(n, seed=42):
+    """Simulate a daily return history for n assets.
+
+    The assets differ in their true drift and load on a handful of common
+    factors, so the estimated problem below has a genuinely tilted frontier
+    instead of one driven purely by estimation noise.
+
+    Args:
+        n (int): Number of assets.
+        seed (int): Seed for the random generator, so the notebook is reproducible.
+
+    Returns:
+        polars.DataFrame: A frame with a ``date`` column and one return column
+            per asset, with ``HISTORY`` rows.
+    """
+    rng = np.random.default_rng(seed)
+
+    # True annual drifts spread across the assets, expressed per day.
+    drift = np.linspace(0.02, 0.20, n) / PERIODS
+    # A low-rank common factor structure plus idiosyncratic noise.
+    k = max(2, n // 10)
+    exposures = rng.standard_normal((n, k)) * 0.4
+    factors = rng.standard_normal((HISTORY, k)) * 0.01
+    idiosyncratic = rng.standard_normal((HISTORY, n)) * 0.01
+
+    returns = drift + factors @ exposures.T + idiosyncratic
+    columns = [f"asset_{i:03d}" for i in range(n)]
+    return pl.DataFrame({"date": business_dates(HISTORY), **dict(zip(columns, returns.T, strict=True))})
+
+
+@app.function(hide_code=True)
+def cla(returns):
     """Compute using the Critical Line Algorithm (CLA) an efficient frontier.
 
+    The mean vector and the covariance matrix are the sample estimates taken
+    from the simulated return history. The portfolios are long-only, capped at
+    100% per name, and fully invested.
+
     Args:
-        n (int): The dimension size of the mean vector, lower and upper bounds
-            arrays, and covariance matrix used in the computation.
+        returns (polars.DataFrame): Return history with a leading ``date`` column.
 
     Returns:
-        numpy.ndarray: The efficient frontier generated by the CLA based on the
-            provided parameters.
+        cvxcla.types.Frontier: The efficient frontier of the estimated problem.
     """
-    mean = np.random.randn(n)
-    lower_bounds = np.zeros(n)
-    upper_bounds = np.ones(n)
-
-    factor = np.random.randn(n, n)
-    covariance = factor @ factor.T
-
-    f1 = CLA(
-        mean=mean,
-        covariance=covariance,
-        lower_bounds=lower_bounds,
-        upper_bounds=upper_bounds,
-        a=np.ones((1, len(mean))),
+    matrix = returns.drop("date").to_numpy()
+    n = matrix.shape[1]
+
+    return CLA(
+        mean=matrix.mean(axis=0),
+        covariance=np.cov(matrix, rowvar=False),
+        lower_bounds=np.zeros(n),
+        upper_bounds=np.ones(n),
+        a=np.ones((1, n)),
         b=np.ones(1),
     ).frontier
-    return f1
 
 
 @app.cell
 def _(slider):
-    frontier = cla(slider.value)
-    frontier.interpolate(2).plot(volatility=True, markers=True)
-    frontier.plot()
+    returns = simulate(slider.value)
+    frontier = cla(returns)
+    mo.md(f"The frontier of the estimated problem has **{len(frontier)}** turning points.")
+    return frontier, returns
+
+
+@app.cell
+def _(frontier):
+    frontier.plot(volatility=True, markers=True)
+    return
+
+
+@app.cell
+def _():
+    mo.md(
+        r"""
+    ## From weights to a track record
+
+    The frontier is a set of weight vectors. Applied to the return history they
+    generated, each one becomes a return series, and a return series is what
+    jQuantStats analyses. We look at three portfolios: the maximum-Sharpe point
+    on the frontier, the minimum-variance point, and equal weight as a
+    reference.
+
+    These are *in-sample* numbers -- the same history produced the estimates the
+    optimiser used -- so read the table as a description of the frontier, not as
+    a backtest. `experiments/frontier_stats.py` runs the out-of-sample version
+    on real S&P 500 data.
+    """
+    )
+    return
+
+
+@app.function(hide_code=True)
+def track_records(frontier, returns):
+    """Turn frontier portfolios into a jQuantStats `Data` object.
+
+    Args:
+        frontier (cvxcla.types.Frontier): The traced efficient frontier.
+        returns (polars.DataFrame): The return history the frontier was estimated on.
+
+    Returns:
+        jquantstats.Data: The realised return series of the maximum-Sharpe, the
+            minimum-variance and the equal-weight portfolio.
+    """
+    matrix = returns.drop("date").to_numpy()
+    n = matrix.shape[1]
+
+    _, max_sharpe = frontier.max_sharpe
+    # Frontier order runs from maximum return towards minimum variance, but read
+    # the minimiser off the variance vector rather than relying on that order.
+    min_variance = frontier.weights[int(np.argmin(frontier.variance))]
+
+    portfolios = {
+        "max_sharpe": max_sharpe,
+        "min_variance": min_variance,
+        "equal_weight": np.full(n, 1.0 / n),
+    }
+    series = {name: matrix @ weights for name, weights in portfolios.items()}
+    return Data.from_returns(pl.DataFrame({"date": returns["date"], **series}))
+
+
+@app.cell
+def _(frontier, returns):
+    data = track_records(frontier, returns)
+    return (data,)
+
+
+@app.cell
+def _(data):
+    mo.ui.table(data.stats.summary(), selection=None)
+    return
+
+
+@app.cell
+def _(data):
+    data.plots.returns(title="Cumulative return of three frontier portfolios")
     return
 
 

book/marimo/notebooks/factor.py

@@ -4,6 +4,8 @@
 #     "marimo==0.14.13",
 #     "numpy==2.3.0",
 #     "plotly==6.7.0",
+#     "polars==1.44.1",
+#     "jquantstats==0.11.0",
 #     "cvx-linalg>=0.9.3",
 #     "cvxcla"
 # ]
@@ -22,6 +24,8 @@
 with app.setup:
     import marimo as mo
     import numpy as np
+    import polars as pl
+    from jquantstats import Data
 
     from cvxcla import CLA, FactorCovariance
 
@@ -42,6 +46,9 @@ def _():
        which is exactly what `FactorCovariance` solves via the Woodbury
        identity in $O(nk)$ memory.
     4. Hand it to `CLA` and plot the frontier.
+    5. Check the operator against the data: the volatility the Woodbury
+       quadratic form predicts for a frontier portfolio has to equal the
+       volatility its realised return series actually shows.
     """
     )
     return
@@ -94,7 +101,7 @@ def _(n_slider):
     returns = simulate_returns(rng, t=2 * n, n=n)
     covariance = clip_covariance(returns)
     mo.md(f"Kept **{covariance.k}** factors out of {n} sample eigenvalues.")
-    return covariance, n, rng
+    return covariance, n, returns, rng
 
 
 @app.cell
@@ -117,5 +124,88 @@ def _(frontier):
     return
 
 
+@app.cell
+def _():
+    mo.md(
+        r"""
+    ## Does the operator agree with the data?
+
+    `FactorCovariance` never forms the $n \times n$ matrix, so the volatilities
+    plotted above come out of the Woodbury identity rather than a dense
+    quadratic form. That is worth checking against the returns themselves.
+
+    Below, each frontier portfolio is applied to the simulated history to give a
+    realised return series, which
+    [jQuantStats](https://github.com/jebel-quant/jquantstats) measures. The
+    predicted column is $\sqrt{w^\top \Sigma w}$ from the operator; the
+    realised column is the sample standard deviation of the series, and their
+    ratio sits within a few percent of 1 across the whole slider range. Not
+    exactly 1: clipping deliberately discards the part of the sample spectrum
+    it calls noise, so the cleaned $\Sigma$ is *not* the sample covariance of
+    this history. A few percent is that discarded noise. An order-of-magnitude
+    gap, or one that widened with $n$, would instead point at the low-rank
+    solve.
+
+    Risk is all we ask of this table. The means fed to the CLA are a synthetic
+    forecast unrelated to the simulation, and `simulate_returns` draws from a
+    standard normal rather than at a realistic return magnitude, so the
+    return-, Sharpe- and drawdown-based metrics jQuantStats also offers would
+    be measuring the simulation's conventions rather than the frontier.
+    `experiments/frontier_stats.py` reads those metrics off real S&P 500 data
+    instead.
+    """
+    )
+    return
+
+
+@app.function(hide_code=True)
+def score(frontier, returns):
+    """Compare each portfolio's predicted volatility with its realised one.
+
+    Args:
+        frontier (cvxcla.types.Frontier): The traced efficient frontier.
+        returns (numpy.ndarray): The t x n simulated return history.
+
+    Returns:
+        polars.DataFrame: One row per portfolio with the volatility the
+            covariance operator predicts, the volatility jQuantStats measures on
+            the realised series, and their ratio.
+    """
+    t, n = returns.shape
+    # Frontier order runs from maximum return towards minimum variance, but read
+    # the minimiser off the variance vector rather than relying on that order.
+    min_variance = frontier.weights[int(np.argmin(frontier.variance))]
+    portfolios = {
+        "max_sharpe": frontier.max_sharpe[1],
+        "min_variance": min_variance,
+        "equal_weight": np.full(n, 1.0 / n),
+    }
+
+    calendar = pl.date_range(pl.date(2015, 1, 1), pl.date(2045, 1, 1), interval="1d", eager=True)
+    dates = calendar.filter(calendar.dt.weekday() <= 5).head(t)
+    series = {name: returns @ weights for name, weights in portfolios.items()}
+    data = Data.from_returns(pl.DataFrame({"date": dates, **series}))
+
+    # Per-observation volatility, so it is comparable with the frontier's own:
+    # the simulation carries no annualisation convention.
+    realised = data.stats.volatility(annualize=False)
+    # The Woodbury operator's own quadratic form -- no dense matrix is formed.
+    predicted = {name: float(np.sqrt(w @ frontier.covariance.matvec(w))) for name, w in portfolios.items()}
+    return pl.DataFrame(
+        {
+            "portfolio": list(portfolios),
+            "predicted volatility": [predicted[name] for name in portfolios],
+            "realised volatility": [realised[name] for name in portfolios],
+            "ratio": [realised[name] / predicted[name] for name in portfolios],
+        }
+    )
+
+
+@app.cell
+def _(frontier, returns):
+    mo.ui.table(score(frontier, returns), selection=None)
+    return
+
+
 if __name__ == "__main__":
     app.run()