Coverage for src/basanos/math/optimizer.py: 100%
49 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-11-03 18:03 +0000
1"""Correlation-aware risk position optimizer (Basanos).
3This module provides utilities to compute correlation-adjusted risk positions
4from price data and expected-return signals. It relies on volatility-adjusted
5returns to estimate a dynamic correlation matrix (via EWM), applies shrinkage
6towards identity, and solves a normalized linear system per timestamp to
7obtain stable positions.
9Performance characteristics
10---------------------------
11Let *N* be the number of assets and *T* the number of timestamps.
13**Computational complexity**
15+----------------------------------+------------------+--------------------------------------+
16| Operation | Complexity | Bottleneck |
17+==================================+==================+======================================+
18| EWM volatility (``ret_adj``, | O(T·N) | Linear in both T and N; negligible |
19| ``vola``) | | |
20+----------------------------------+------------------+--------------------------------------+
21| EWM correlation (``cor``) | O(T·N²) | ``ewm_covariance`` from |
22| | | ``cvx.linalg`` over all N² pairs |
23+----------------------------------+------------------+--------------------------------------+
24| Linear solve per timestamp | O(N³) | Cholesky / LU per row in |
25| (``cash_position``) | * T solves | ``cash_position`` |
26+----------------------------------+------------------+--------------------------------------+
28**Memory usage** (peak, approximate)
30``ewm_covariance`` from ``cvx.linalg`` processes the input Polars DataFrame
31and returns a dict of covariance matrices. Peak RAM ≈ **O(T · N²)** bytes.
32Typical working sizes on a 16 GB machine:
34+--------+--------------------------+------------------------------------+
35| N | T (daily rows) | Peak memory (approx.) |
36+========+==========================+====================================+
37| 50 | 252 (~1 yr) | ~70 MB |
38+--------+--------------------------+------------------------------------+
39| 100 | 252 (~1 yr) | ~280 MB |
40+--------+--------------------------+------------------------------------+
41| 100 | 2 520 (~10 yr) | ~2.8 GB |
42+--------+--------------------------+------------------------------------+
43| 200 | 2 520 (~10 yr) | ~11 GB |
44+--------+--------------------------+------------------------------------+
45| 500 | 2 520 (~10 yr) | ~70 GB ⚠ exceeds typical RAM |
46+--------+--------------------------+------------------------------------+
48**Practical limits (daily data)**
50* **≤ 150 assets, ≤ 5 years** — well within reach on an 8 GB laptop.
51* **≤ 250 assets, ≤ 10 years** — requires ~11-12 GB; feasible on a 16 GB
52 workstation.
53* **> 500 assets with multi-year history** — peak memory exceeds 16 GB;
54 reduce the time range or switch to a chunked / streaming approach.
55* **> 1 000 assets** — the O(N³) per-solve cost alone makes real-time
56 optimization impractical even with adequate RAM.
58See ``BENCHMARKS.md`` for measured wall-clock timings across representative
59dataset sizes.
61Internal structure
62------------------
63The implementation is split across focused private modules to keep each file
64readable and independently testable:
66* `_config` — `BasanosConfig` and all
67 covariance-mode configuration classes.
68* `_engine_base` — `_BatchCore`, the validated ``prices`` / ``mu`` /
69 ``cfg`` dataclass composing `_CoreDataMixin`, `_SolveMixin` and
70 `_SignalEvaluatorMixin`.
71 `BasanosEngine` subclasses it; the streaming warmup builds it directly.
72* `_engine_validation` — free functions that validate the
73 ``prices`` / ``mu`` / ``cfg`` inputs.
74* `_engine_core` — the `_CoreDataMixin` providing the
75 core data-access properties (``assets``, ``ret_adj``, ``vola``, ``cor``,
76 ``cor_tensor``).
77* `_engine_solve` — private helpers providing the
78 ``_iter_matrices`` and ``_iter_solve`` generators (per-timestamp solve
79 logic).
80* `_engine_diagnostics` — private helpers providing
81 matrix-quality diagnostics (condition number, effective rank, solver
82 residual, signal utilisation).
83* `_engine_performance` — the `_PerformanceMixin` providing the
84 portfolio-Sharpe sweep helpers (``sharpe_at_shrink``,
85 ``sharpe_at_window_factors``, ``naive_sharpe``).
86* `_engine_ic` — private helpers providing signal
87 evaluation metrics (IC, Rank IC, ICIR, and summary statistics).
88* This module — `BasanosEngine`, a single flat class that wires
89 every method together in clearly delimited sections.
90"""
92import dataclasses
94import numpy as np
95import polars as pl
96from jquantstats import Portfolio
98from ._config import (
99 BasanosConfig,
100 CovarianceConfig,
101 CovarianceMode,
102 EwmaShrinkConfig,
103 SlidingWindowConfig,
104)
105from ._config_report import ConfigReport
106from ._engine_base import _BatchCore
107from ._engine_diagnostics import _DiagnosticsMixin
108from ._engine_performance import _PerformanceMixin
110# ---------------------------------------------------------------------------
111# Re-export config symbols so ``from basanos.math.optimizer import …`` keeps
112# working for existing callers.
113# ---------------------------------------------------------------------------
114__all__ = [
115 "BasanosConfig",
116 "BasanosEngine",
117 "CovarianceConfig",
118 "CovarianceMode",
119 "EwmaShrinkConfig",
120 "SlidingWindowConfig",
121]
124@dataclasses.dataclass(frozen=True)
125class BasanosEngine(_BatchCore, _DiagnosticsMixin, _PerformanceMixin):
126 """Engine to compute correlation matrices and optimize risk positions.
128 Encapsulates price data and configuration to build EWM-based
129 correlations, apply shrinkage, and solve for normalized positions.
131 Public methods are organised into clearly delimited sections (some
132 inherited from the private mixin classes):
134 * **Core data access** — `assets`, `ret_adj`, `vola`, `cor`, `cor_tensor`
135 * **Solve / position logic** — `cash_position`, `position_status`,
136 `risk_position`, `position_leverage`, `warmup_state`
137 * **Portfolio and performance** — `portfolio`, `naive_sharpe`,
138 `sharpe_at_shrink`, `sharpe_at_window_factors`
139 * **Matrix diagnostics** — `condition_number`, `effective_rank`,
140 `solver_residual`, `signal_utilisation`
141 * **Signal evaluation** — `ic(h)`, `rank_ic(h)`, `ic_mean(h)`, `ic_std(h)`,
142 `icir(h)`, `rank_ic_mean(h)`, `rank_ic_std(h)` (``h`` defaults to 1)
143 * **Reporting** — `config_report`
145 Data-flow diagram
146 -----------------
148 .. code-block:: text
150 prices (pl.DataFrame)
151 │
152 ├─ vol_adj ──► ret_adj (volatility-adjusted log returns)
153 │ │
154 │ ├─ ewm_covariance ──► cor / cor_tensor
155 │ │ │
156 │ │ └─ shrink2id / FactorModel
157 │ │ │
158 │ vola covariance matrix
159 │ │ │
160 └── mu ──────────┴── _iter_solve ──────────┘
161 │
162 cash_position
163 │
164 ┌────────┴────────┐
165 portfolio diagnostics
166 (Portfolio) (condition_number,
167 effective_rank,
168 solver_residual,
169 signal_utilisation,
170 ic, rank_ic, …)
172 Attributes:
173 prices: Polars DataFrame of price levels per asset over time. Must
174 contain a ``'date'`` column and at least one numeric asset column
175 with strictly positive values that are not monotonically
176 non-decreasing or non-increasing (i.e. they must vary in sign).
177 mu: Polars DataFrame of expected-return signals aligned with *prices*.
178 Must share the same shape and column names as *prices*.
179 cfg: Immutable `BasanosConfig` controlling EWMA half-lives,
180 clipping, shrinkage intensity, and AUM.
182 Examples:
183 Build an engine with two synthetic assets over 30 days and inspect the
184 optimized positions and diagnostic properties.
186 >>> import numpy as np
187 >>> import polars as pl
188 >>> from basanos.math import BasanosConfig, BasanosEngine
189 >>> dates = list(range(30))
190 >>> rng = np.random.default_rng(42)
191 >>> prices = pl.DataFrame({
192 ... "date": dates,
193 ... "A": np.cumprod(1 + rng.normal(0.001, 0.02, 30)) * 100.0,
194 ... "B": np.cumprod(1 + rng.normal(0.001, 0.02, 30)) * 150.0,
195 ... })
196 >>> mu = pl.DataFrame({
197 ... "date": dates,
198 ... "A": rng.normal(0.0, 0.5, 30),
199 ... "B": rng.normal(0.0, 0.5, 30),
200 ... })
201 >>> cfg = BasanosConfig(vola=5, corr=10, clip=2.0, shrink=0.5, aum=1_000_000)
202 >>> engine = BasanosEngine(prices=prices, mu=mu, cfg=cfg)
203 >>> engine.assets
204 ['A', 'B']
205 >>> engine.cash_position.shape
206 (30, 3)
207 >>> engine.position_leverage.columns
208 ['date', 'leverage']
209 """
211 # ------------------------------------------------------------------
212 # Fields and validation — inherited from _BatchCore
213 # ------------------------------------------------------------------
214 # (prices, mu, cfg, __post_init__ -> _validate_inputs)
215 # Defined in _engine_base.py, which BasanosStream also builds on.
217 # ------------------------------------------------------------------
218 # Core data-access properties — inherited from _CoreDataMixin
219 # ------------------------------------------------------------------
220 # (assets, ret_adj, vola, cor, cor_tensor)
221 # Implementations live in _engine_core.py.
223 # ------------------------------------------------------------------
224 # Internal solve helpers — inherited from _SolveMixin
225 # ------------------------------------------------------------------
226 # (_compute_mask, _check_signal, _scale_to_cash, _row_early_check,
227 # _denom_guard_yield, _compute_position, _replay_positions,
228 # _iter_matrices, _iter_solve, warmup_state)
229 # Implementations live in _engine_solve.py; patch targets remain in that
230 # module's namespace, e.g. ``patch("basanos.math._engine_solve.solve")``.
232 # ------------------------------------------------------------------
233 # Position properties
234 # ------------------------------------------------------------------
236 @property
237 def cash_position(self) -> pl.DataFrame:
238 r"""Optimize correlation-aware risk positions for each timestamp.
240 Supports two covariance modes controlled by ``cfg.covariance_config``:
242 * `EwmaShrinkConfig` (default): Computes EWMA correlations, applies
243 linear shrinkage toward the identity, and solves a normalised linear
244 system $C\,x = \mu$ per timestamp via Cholesky / LU.
246 * `SlidingWindowConfig`: At each timestamp uses the
247 ``cfg.covariance_config.window`` most recent vol-adjusted returns to fit a
248 rank-``cfg.covariance_config.n_factors`` factor model via truncated SVD and
249 solves the system via the Woodbury identity at $O(k^3 + kn)$ rather
250 than $O(n^3)$ per step.
252 Non-finite or ill-posed cases yield zero positions for safety.
254 Returns:
255 pl.DataFrame: DataFrame with columns ['date'] + asset names containing
256 the per-timestamp cash positions (risk divided by EWMA volatility).
258 Performance:
259 For ``ewma_shrink``: dominant cost is ``self.cor`` (O(T·N²) time,
260 O(T·N²) memory). The per-timestamp
261 linear solve adds O(N³) per row.
263 For ``sliding_window``: O(T·W·N·k) for sliding SVDs plus
264 O(T·(k³ + kN)) for Woodbury solves. Memory is O(W·N) per step,
265 independent of T.
266 """
267 assets = self.assets
269 # Compute risk positions row-by-row using _replay_positions.
270 prices_num = self.prices.select(assets).to_numpy()
272 risk_pos_np = np.full_like(prices_num, fill_value=np.nan, dtype=float)
273 cash_pos_np = np.full_like(prices_num, fill_value=np.nan, dtype=float)
274 vola_np = self.vola.select(assets).to_numpy()
276 self._replay_positions(risk_pos_np, cash_pos_np, vola_np)
278 # Build Polars DataFrame for cash positions (numeric columns only)
279 cash_position = self.prices.with_columns(
280 [(pl.lit(cash_pos_np[:, i]).alias(asset)) for i, asset in enumerate(assets)]
281 )
283 return cash_position
285 @property
286 def position_status(self) -> pl.DataFrame:
287 """Per-timestamp reason code explaining each `cash_position` row.
289 Labels every row with exactly one of four `SolveStatus`
290 codes (which compare equal to their string equivalents):
292 * ``'warmup'``: Insufficient history for the sliding-window
293 covariance mode (``i + 1 < cfg.covariance_config.window``).
294 Positions are ``NaN`` for all assets at this timestamp.
295 * ``'zero_signal'``: The expected-return vector ``mu`` was
296 all-zeros (or all-NaN) at this timestamp; the optimizer
297 short-circuited and returned zero positions without solving.
298 * ``'degenerate'``: The normalisation denominator was non-finite
299 or below ``cfg.denom_tol``, the Cholesky / Woodbury solve
300 failed, or no asset had a finite price; positions were zeroed
301 for safety.
302 * ``'valid'``: The linear system was solved successfully and
303 positions are non-trivially non-zero.
305 The codes map one-to-one onto the three NaN / zero cases
306 described in the issue and allow downstream consumers (backtests,
307 risk monitors) to distinguish data gaps from signal silence from
308 numerical ill-conditioning without re-inspecting ``mu`` or the
309 engine configuration.
311 Returns:
312 pl.DataFrame: Two-column DataFrame ``{'date': ..., 'status': ...}``
313 with one row per timestamp. The ``status`` column has
314 ``Polars`` dtype ``String``.
315 """
316 statuses = [status for _i, _t, _mask, _pos, status in self._iter_solve()]
317 return pl.DataFrame({"date": self.prices["date"], "status": pl.Series(statuses, dtype=pl.String)})
319 @property
320 def risk_position(self) -> pl.DataFrame:
321 """Risk positions (before EWMA-volatility scaling) at each timestamp.
323 Derives the un-volatility-scaled position by multiplying the cash
324 position by the per-asset EWMA volatility. Equivalently, this is
325 the quantity solved by the correlation-adjusted linear system before
326 dividing by ``vola``.
328 Relationship to other properties::
330 cash_position = risk_position / vola
331 risk_position = cash_position * vola
333 Returns:
334 pl.DataFrame: DataFrame with columns ``['date'] + assets`` where
335 each value is ``cash_position_i * vola_i`` at the given timestamp.
336 """
337 assets = self.assets
338 cp_np = self.cash_position.select(assets).to_numpy()
339 vola_np = self.vola.select(assets).to_numpy()
340 with np.errstate(invalid="ignore"):
341 risk_pos = cp_np * vola_np
342 return self.prices.with_columns([pl.lit(risk_pos[:, i]).alias(asset) for i, asset in enumerate(assets)])
344 @property
345 def position_leverage(self) -> pl.DataFrame:
346 """L1 norm of cash positions (gross leverage) at each timestamp.
348 Sums the absolute values of all asset cash positions at each row.
349 NaN positions are treated as zero (they contribute nothing to gross
350 leverage).
352 Returns:
353 pl.DataFrame: Two-column DataFrame ``{'date': ..., 'leverage': ...}``
354 where ``leverage`` is the L1 norm of the cash-position vector.
355 """
356 assets = self.assets
357 cp_np = self.cash_position.select(assets).to_numpy()
358 leverage = np.nansum(np.abs(cp_np), axis=1)
359 return pl.DataFrame({"date": self.prices["date"], "leverage": pl.Series(leverage, dtype=pl.Float64)})
361 # ------------------------------------------------------------------
362 # Portfolio and performance
363 # ------------------------------------------------------------------
365 @property
366 def portfolio(self) -> Portfolio:
367 """Construct a Portfolio from the optimized cash positions.
369 Converts the computed cash positions into a Portfolio using the
370 configured AUM. The ``cost_per_unit`` from `cfg` is forwarded
371 so that `net_cost_nav` and
372 `position_delta_costs` work out
373 of the box without any further configuration.
375 Returns:
376 Portfolio: Instance built from cash positions with AUM scaling.
377 """
378 cp = self.cash_position
379 assets = [c for c in cp.columns if c != "date" and cp[c].dtype.is_numeric()]
380 scaled = cp.with_columns(pl.col(a) * self.cfg.position_scale for a in assets)
381 return Portfolio.from_cash_position(self.prices, scaled, aum=self.cfg.aum, cost_per_unit=self.cfg.cost_per_unit)
383 # ------------------------------------------------------------------
384 # Performance sweeps — inherited from _PerformanceMixin
385 # ------------------------------------------------------------------
386 # (sharpe_at_shrink, sharpe_at_window_factors, naive_sharpe)
387 # Implementations live in _engine_performance.py.
389 # ------------------------------------------------------------------
390 # Reporting
391 # ------------------------------------------------------------------
393 @property
394 def config_report(self) -> "ConfigReport":
395 """Return a `ConfigReport` facade for this engine.
397 Returns a `ConfigReport` that
398 includes the full **lambda-sweep chart** — an interactive plot of the
399 annualised Sharpe ratio as `shrink` (λ) is swept
400 across [0, 1] — in addition to the parameter table, shrinkage-guidance
401 table, and theory section available from
402 `report`.
404 Returns:
405 basanos.math._config_report.ConfigReport: Report facade with
406 ``to_html()`` and ``save()`` methods.
408 Examples:
409 >>> import numpy as np
410 >>> import polars as pl
411 >>> from basanos.math.optimizer import BasanosConfig, BasanosEngine
412 >>> dates = pl.Series("date", list(range(200)))
413 >>> rng = np.random.default_rng(0)
414 >>> prices = pl.DataFrame({"date": dates, "A": rng.lognormal(size=200), "B": rng.lognormal(size=200)})
415 >>> mu = pl.DataFrame({"date": dates, "A": rng.normal(size=200), "B": rng.normal(size=200)})
416 >>> cfg = BasanosConfig(vola=10, corr=20, clip=3.0, shrink=0.5, aum=1e6)
417 >>> engine = BasanosEngine(prices=prices, mu=mu, cfg=cfg)
418 >>> report = engine.config_report
419 >>> html = report.to_html()
420 >>> "Lambda" in html
421 True
422 """
423 return ConfigReport(config=self.cfg, engine=self)
425 # ------------------------------------------------------------------
426 # Matrix diagnostics — inherited from _DiagnosticsMixin
427 # ------------------------------------------------------------------
428 # (condition_number, effective_rank, solver_residual, signal_utilisation)
429 # Implementations live in _engine_diagnostics.py; patch targets remain in
430 # that module's namespace, e.g.
431 # ``patch("basanos.math._engine_diagnostics.solve")``.
433 # ------------------------------------------------------------------
434 # Signal evaluation — inherited from _SignalEvaluatorMixin
435 # ------------------------------------------------------------------
436 # (_ic_series, ic, rank_ic, ic_mean, ic_std, icir,
437 # rank_ic_mean, rank_ic_std)
438 # Implementations live in _engine_ic.py.