Portfolio_Risk_Management/docs/METHODOLOGY.md

# Methodology

Every formula the services implement, with the conventions and default parameters. Section
numbers follow the project brief ("Portfolio Risk Measurement and Optimization").

## 1. Investment universe (`quantcore.universe`)

* **Candidates.** US: the S&P 500 constituents (Wikipedia). Europe (EU, CH, NO, UK): the STOXX
  Europe 600 constituents, read from the iShares STOXX Europe 600 UCITS ETF holdings file, which
  gives issuer ticker, exchange, listing currency, country and GICS sector. Domiciles outside
  EU/CH/NO/UK (e.g. Bermuda, Israel) are excluded. Both lists are snapshotted into
  `quantcore/data/static/` for offline use.
* **Yahoo symbols.** `<code><suffix>` per exchange (`.SW .L .DE .PA .AS .MI .MC .ST .CO .OL .HE
  .BR .VI .WA .IR .LS .AT`), share classes with a dash (`NOVO-B.CO`, `BRK-B`), London codes lose
  the trailing dot (`RR.` → `RR.L`). All 1 117 symbols were verified against Yahoo, and the
  European prices were cross-checked against the ETF file (after FX) — zero mismatches.
* **Metadata** from `Ticker.info`: long name, domicile, listing currency (GBp → GBP with a price
  scale of 100), exchange, Yahoo sector/industry, market cap, shares outstanding, first trade date.
* **Selection.** One listing per company (share classes collapsed on a normalised name, largest cap
  kept) → rank by market cap in CHF within each region → keep the largest 500 US and 500 European
  → drop stocks whose first price is after `required_history_from` (2005-01-31) or whose coverage
  since then is below 90 % of union-calendar days. The result has fewer than 1 000 names.
* **Market capitalisation** (for index weights) at any date: shares outstanding (today's, held
  constant) × adjusted close × FX to CHF.
* **Price cleaning** (`quantcore.data.cleaning`, applied by the setup and download services and
  re-runnable with `quantcore clean`): non-positive adjusted closes (broken dividend/rights
  adjustments, e.g. ABI.BR before 2008) → NaN; *spike blocks* — at most three days whose entry
  move `|r| > 0.4` is reversed within the block (stale pre-split ticks on exchange holidays, e.g.
  2× ticks in NOVO-B.CO) → NaN. Genuine jumps are not reversed and therefore untouched. Stocks
  whose cleaned series has a gap longer than `max_missing_run` (20 business days) after
  2005-01-31 are excluded from the universe. On the 2026-09-21 data: 1 347 non-positive values in
  2 symbols, 28 spike ticks in 8 symbols; final universe 745 stocks (382 US, 363 European).
* **FX quoting-unit repair** (`clean_fx`, same module and the same two entry points): Yahoo returns
  a cross rate in the wrong unit for whole stretches of days. NOK/CHF came back **100× too small on
  96 observations** — 2016-08-15 to 2017-01-09 in three runs, plus 2022-11-04 — which put the
  measured volatility of the *currency* at 337 % a year instead of 12 %, and with it every Norwegian
  stock (a CHF return is the local return plus the FX return, so EQNR.OL measured 338 % against its
  true 31 %, and its Black–Litterman implied return came out above 100 % a year). A day-on-day
  log-change within 0.15 of `k·ln(10)` is treated as a unit change, accumulated, and the series is
  put back on the scale that most of its observations share (ties go to the scale it ends on).
  The tolerance is tight on purpose: a unit change is *exactly* a power of ten apart from the day's
  own move, while a loose tolerance swallows real events — at 0.5 it flagged ABVX.PA's +510 %
  trial-result day (ratio 6.1) as a decimal shift. A scan of all 1 117 price series with the tight
  rule finds one further case, ROSE.L (11 observations quoted in pence rather than GBP), which is
  not in the investable universe.
* **FX spike ticks** (same function, thresholds an order of magnitude tighter than for prices:
  entry move `|r| > 0.05` reversed within 0.02 over at most three days). Five ticks in twenty-one
  years: EUR, GBP and USD on 2009-02-06 (one bad CHF quote hits every cross that day), SEK on
  2010-09-13, NOK on 2022-10-10. Removing them takes the lag-1 autocorrelation of daily FX
  log-returns to about zero — EUR −0.053 → −0.005, GBP −0.040 → −0.001, USD −0.034 → −0.002 — which
  is what a clean daily close series looks like. The one violent day in the sample, 2015-01-16
  (the SNB abandons the euro floor, every CHF cross moves ~18 %), does not reverse and is kept.

**Why the rates are taken directly against CHF and not crossed through USD.** USD pairs are the most
liquid instruments, so triangulating `CUR/CHF = (USD/CHF) / (USD/CUR)` looks attractive; measured on
this feed it is worse, and the package therefore uses the direct pair wherever it has history (the
USD cross remains the fallback, and is what PLN uses, since `PLNCHF=X` only starts in 2026):

| | direct `CURCHF=X` | crossed through USD |
|---|---|---|
| EUR volatility p.a. / lag-1 AC | 7.49 % / −0.005 | 9.99 % / **−0.219** |
| NOK volatility p.a. / lag-1 AC | 11.45 % / −0.047 | 17.14 % / **−0.330** |
| GBP volatility p.a. | 10.10 % | 10.35 % |
| reversing bad ticks, 2005–2026 | 5 over all currencies | 8, five of them in `EUR=X` alone |

The cross combines two series snapped at different times, and the non-synchronicity shows up as a
strongly negative lag-1 autocorrelation and inflated variance — a portfolio risk figure built on it
would be overstated. It also has its own bad ticks (`NOK=X` falls 33 % on 2020-03-20 and rises 42 %
back on 2020-03-23; `EUR=X` moves −16 % then +14 % on 2008-12-08/09). Triangulation buys only the
removal of stale quotes (SEK 3.9 % of days unchanged, DKK 1.6 %), which costs far less than it adds.

## 2. Prices, FX and log-returns (`quantcore.data`, `quantcore.measurement.returns`)

* Adjusted closes `P_i(t)` in the listing currency on the union calendar of all exchanges, holidays
  forward-filled (≤ 5 days; leading NaNs untouched).
* FX frame `X_c(t)` = CHF per unit of currency `c`, from Yahoo `cCHF=X`, crossed through USD when
  the direct pair lacks history (`c/CHF = (USD/CHF)/(USD/c)`, e.g. PLN).
* Log-returns `r_i(t) = ln P_i(t) − ln P_i(t−1)`; CHF returns add the FX log-return:
  `r_i^CHF(t) = r_i(t) + ln X_{c(i)}(t) − ln X_{c(i)}(t−1)`.
* Weekly returns: Friday-to-Friday closes; monthly: month-end to month-end.
* Horizon `h` in periods: day = 1, week = 5, month = 21 (daily data); at weekly frequency week = 1,
  month = 4. Returns are assumed i.i.d., so `h`-period returns are sums.
* Portfolio log-return with weights `w` fixed at the start of the period (exact aggregation):
  `r_p = ln(1 + Σ_i w_i (e^{r_i} − 1))`.

## 3. Estimators (`quantcore.measurement.estimators`)

Two families, both centred (the corresponding mean is subtracted):

* **Moving window** of `w` observations (default 200 days): sample mean, standard deviation
  (ddof 1), covariance `np.cov`, correlation.
* **Exponentially weighted** with decay `λ` (defaults 0.94 for measurement, 0.99 for the slow
  optimisation covariance): weight of the observation `k` periods back ∝ `λ^k`, normalised to one
  (pandas `ewm(alpha=1−λ, adjust=True)`); covariance
  `Σ = Σ_t ω_t (r_t − m)(r_t − m)'`, `m = Σ_t ω_t r_t`. Effective window ≈ `1/(1−λ)`.
* Streamed sequences `Σ_{t−1}` for every `t` (needed by the innovation filter): the rolling
  estimator recomputes the window, the EWM estimator uses the exact centred recursion
  `m_t = m_{t−1} + α d`, `Σ_t = λ Σ_{t−1} + α d (r_t − m_t)'`, `d = r_t − m_{t−1}`, `α = 1 − λ`.
* Linear shrinkage towards the diagonal `(1−s)Σ + s·diag(Σ)` with `s` either given or `"auto"`:
  the Schäfer–Strimmer / Ledoit–Wolf intensity
  `s* = Σ_{i≠j} Var̂(s_ij) / Σ_{i≠j} s_ij²` (clipped to [0, 1]) computed on the effective estimation
  window (`w` observations, or `3/(1−λ)` for EWM). It is ≈ 0 for a long, strongly correlated
  sample and → 1 when `N` approaches the number of observations — exactly when the sample
  covariance becomes ill-conditioned. `"auto"` is the default of the optimisation covariance
  (≈ 0.05 for 100 stocks / λ = 0.99, ≈ 0.11 for the full universe).
* Annualisation: means × periods per year (252 / 52 / 12), volatilities × √(periods per year).

## 4. Risk measures (`quantcore.measurement.risk`)

Losses are positive. For a log-return `r` over the horizon the loss is `−r` in return space and
`V₀(1 − e^{r})` in CHF (exact, monotone mapping — the CHF VaR/CVaR are obtained by applying it to
the return-space figures, or directly from money P&L for scenario sets).

* **Gaussian (GBM)**: with per-period mean `μ` and volatility `σ` from either estimator,
  `μ_h = hμ`, `σ_h = √h σ`,
  `VaR_β = −μ_h + z_β σ_h`, `CVaR_β = −μ_h + σ_h φ(z_β)/(1−β)`.
  Reference: `σ = 1 %`, `μ = 0`: VaR₉₅ = 1.645 %, CVaR₉₅ = 2.063 %, VaR₉₉ = 2.326 %, CVaR₉₉ = 2.665 %.
* **Historical**: `VaR_β = −q_{1−β}(r_h)`, `CVaR_β = −E[r_h | r_h ≤ q_{1−β}]` on the last
  `max(5w, 750)` observations aggregated to overlapping `h`-period sums.
* **Time series** (`risk_series`): rolling / EWM volatilities and the Gaussian VaR/CVaR built from
  them, plus rolling historical quantiles and tail means.
* **Portfolio** (`portfolio_risk`): holdings in shares (+ cash) → value weights `w` at today's CHF
  prices; Gaussian with `σ_p² = w'Σw`, historical on the realised portfolio return, or on a scenario
  set (`method="scenarios"`).
* **Scenario CVaR** (`scenario_cvar`): with `S` equiprobable losses sorted descending and
  `k = (1−β)S`, `CVaR_β = (Σ_{s≤⌊k⌋} L_(s) + (k−⌊k⌋) L_(⌊k⌋+1)) / k` — exactly what the
  Rockafellar–Uryasev LP reproduces; `VaR_β = L_(⌈k⌉)`.

## 5. Performance (`quantcore.measurement.performance`)

Cumulative return `e^{Σr} − 1`; annualised log-return `ppy · mean(r)` and its geometric version
`e^{ppy·mean} − 1`; annualised volatility; Sharpe `(ann. return − r_f)/ann. vol` with the CHF
risk-free rate (`QUANTCORE_RISK_FREE_RATE`, default 0); Sortino (downside deviation); drawdown of
the cumulative value path, max drawdown, Calmar; skewness, excess kurtosis, hit ratio, best/worst
period. Rolling (window) and EWM (λ) versions as time series. Portfolios in shares: constant-mix
(today's value weights) or buy-and-hold.

## 6. Stock-price models (`quantcore.models.simulation`)

Both assume i.i.d. one-period log-returns; `Σ_T` is the latest matrix of the chosen
`CovarianceModel` (e.g. 200-day equally weighted, or EWM λ = 0.99).

* **Geometric Brownian motion**: `r_h ~ N(hm, hΣ_T)`, sampled with `A z`, `A A' = Σ_T`
  (eigen-root with eigenvalues descending by default, so the first draw coordinate drives the market
  factor; symmetric and Cholesky roots optional).
* **Empirical innovations** (filtered historical simulation):
  `ε_t = Σ_{t−1}^{−1/2} r_t` with `Σ_{t−1}` estimated from data up to `t−1` (streamed sequence,
  symmetric root `V Λ^{−1/2} V'` — Cholesky optional for large `N`); forward simulation
  `r(j | T) = hm + Σ_T^{1/2} Σ_{k=1..h} ε_{j_k}` with the dates `j_k` drawn with replacement from
  the history (scrambled Sobol' points on the market-factor ranks by default, `"mc"` uniform or
  `"balanced"` optional). Innovations are centred so the scenario mean equals the supplied drift.
  **Re-standardisation**: a lagged volatility filter leaves the innovations with variance above
  one (the forecast is too low right after every shock — 1.1 for 10 stocks, 1.5 for 100 with
  λ = 0.99), so the centred sample is rescaled before resampling: `standardize="diagonal"`
  (default) divides each asset's innovations by their standard deviation, `"full"` whitens with
  the innovations' sample covariance (needs `T ≥ 3N`). The scenario covariance then reproduces
  `hΣ_T` (per-asset volatility within ±2 %, exactly with full whitening) while skewness, fat
  tails and non-linear co-movements of the data are preserved. Innovations beyond ±8 standard
  deviations are winsorised (a numerical guard affecting ≈ 0.2 % of the values); ill-conditioned
  `Σ_{t−1}` for large `N` is handled by the automatic shrinkage of Section 3.
* **Sampling and moment matching** (`quantcore.models.sampling`): the optimiser is a sample-average
  approximation and the sample mean of a plain Monte-Carlo set is off by `σ/√S` per asset — 3–13 %
  p.a. for `S = 2000` — which the optimiser maximises (measured: two seeds give portfolios differing
  by 100 % of wealth). Therefore: scrambled Sobol' points (`S` rounded up to a power of two),
  **moment matching** `z̃ = C^{-1/2}(z − z̄)` so that the (probability-weighted) scenario mean and
  covariance are *exactly* `hm` and `hΣ_T`, and the model expected returns in the optimiser's
  objective. Only the shape of the tail remains sampled.
* **Drift convention**: the expected-return model delivers *arithmetic* expected returns `μ`
  (per period, annual / 252). The log-drift is the Itô-corrected `m_i = ln(1 + hμ_i)/h − σ_i²/2`,
  so that `E[e^{r_h} − 1] = hμ` exactly (`drift_convention="arithmetic"`); `"log"` uses `μ` as the
  log-drift and adds `σ²/2` to the expected simple return (5 % p.a. for a typical stock, 34 % for
  the most volatile one — an unintended tilt to volatile stocks).
* **Importance sampling** (`importance_sampling="auto"` = `"market"`): the loss direction `d` of
  the *cap-weighted candidate portfolio*; two blocks of 25 % of the scenarios are shifted along it
  by `z_0.95` and `z_0.99` (GBM: `z ↦ z + c_k d`; empirical: dates drawn with exponentially tilted
  probabilities `q_t ∝ exp(ϑ_k d'ε_t)`), so one set serves both confidence levels. Every scenario
  carries the self-normalised likelihood ratio of the defensive mixture (bounded by
  `1/(1 − Σ mixture)`) as its probability `p_s`. Tail noise falls 2–4×; 2 048 tilted scenarios ≈
  8 192 plain ones at β = 0.99. The shift is driven by the *covariance*: in return space it is
  `−c_k Σv/√(v'Σv)`, i.e. each asset moves by its covariance with the reference portfolio (beta),
  not by its own volatility — on the top 30 the correlation of the shift with beta is −1.00 versus
  −0.93 with own volatility, and a negatively correlated stock is shifted *up*. The weighted
  tilted set still reproduces `hΣ` and all correlations exactly (the unweighted one does not). The
  gain a given portfolio gets scales with the correlation between its loss and the reference's
  (`tilt_alignment` in the result: 5.3× at 1.00, 2.3× at 0.90 — the usual value for the optimiser's
  solutions — and nothing below ~0.3, where a warning is raised). `"adaptive"` tilts towards a first-pass solution (sharper at 99 %,
  but the set then depends on the risk parameters); `"off"` disables the tilt.
* **One scenario set per family of runs**: the scenario *key* (`OptimizationService.scenario_key`,
  reported with every result) contains only model inputs — candidates, as-of date, frequency and
  horizon, covariance model, scenario model, `n_scenarios`, sampling, moment matching, seed, tilt
  mode, innovation settings. The CVaR limit, confidence, objective, constraints, costs, cash,
  holdings and the expected-return model are *not* in it: varying them reuses the same scenarios
  (cached in the service), so differences between runs are the effect of the parameter and not
  sampling noise. The drift is applied after the random part (`ScenarioSet.with_drift`), so a new
  view shifts the scenarios without touching the shocks (`shock_hash` is unchanged).
* `ScenarioSet` reports `expected_simple` (the expectation the optimiser uses), `probabilities`,
  the effective sample size and a SHA-256 fingerprint `hash` of the matrix.

## 7. Expected returns — Black–Litterman with Idzorek confidence (`quantcore.models.black_litterman`)

Weekly CHF log-returns of the candidate stocks over the last 10 years (`bl_years`):

* `Σ` = 10-year sample covariance of weekly returns; `w` = current market-cap weights of the
  candidates; `μ_M` = average weekly log-return of the **cap-weighted broad market index** of the
  whole universe (weights drift with prices between rebalancings, today's shares held constant);
  `r_f` = CHF risk-free rate per week.
* Risk aversion `δ = (μ_M − r_f)/(w'Σw)`; equilibrium excess returns `π = δ Σ w` (so the market
  portfolio earns exactly `μ_M − r_f`).
* Views `P μ = Q + ε`, `ε ~ N(0, Ω)`, Ω diagonal. A view is `{assets: {ticker: weight}, q, confidence}`;
  `q` is the annual total return of the view portfolio (converted to weekly excess internally;
  the risk-free part cancels for relative views).
* **Idzorek confidence** `C_k ∈ (0, 1]`: the view's uncertainty is the `ω_k` for which the
  unconstrained posterior weights move the fraction `C_k` of the way from `w` to the weights implied
  by the view at 100 % confidence. For a diagonal Ω that calibration has the closed form
  `ω_k = (1 − C_k)/C_k · p_k (τΣ) p_k'`; the tilt-matching minimisation of the paper is implemented
  as `idzorek_omega_numeric` and the tests show both agree to 1e-4 (Idzorek 2005; the closed form is
  due to Walters 2014).
* Posterior `μ_BL = π + τΣP'(PτΣP' + Ω)^{-1}(Q − Pπ)`,
  `Σ_BL = Σ + τΣ − τΣP'(PτΣP' + Ω)^{-1}PτΣ`; `τ` default 0.05.
* Sentiment → view (`view_from_sentiment`): an assessment `s ∈ [−1, 1]` with degree of belief `C`
  on stock `i` becomes the absolute view `q = π_i^{ann} + r_f + s · strength · σ_i^{ann}`,
  `strength` = 0.25 by default (a fully positive view = +¼ annual volatility above equilibrium).
* Output: annual and per-period posterior total returns (used as scenario drift), Ω, δ, `μ_M`.
  Alternatives: `expected_returns="historical"` (10-year means) or `"zero"`.

## 8. Mean–CVaR optimisation as a MILP (`quantcore.optimization.problem`)

Scenarios `R_{s,i}` = simple horizon returns `e^{r_{s,i}} − 1` in CHF (Section 6, with the drift of
Section 7), prices `P_i`, wealth `W = cash₀ + Σ_i n⁰_i P_i`.

Variables: shares `n_i ∈ ℤ₊` (no shorts); rebalancing adds buys `b_i`, sells `s_i ∈ ℤ₊` with
`n_i = n⁰_i + b_i − s_i`; `α`, `u_s ≥ 0` bounded by the worst gain/loss attainable within the share
caps (free variables make the CVaR objective numerically hard — cuOpt's concurrent method aborts on
such models; the client then retries once with its dual simplex); binaries `y_i` (cardinality / minimum position) and
`t_i` (fixed cost per trade) only when those options are used.

Transaction costs `TC = Σ_i (c_b P_i b_i + c_s P_i s_i) + F Σ_i t_i` (defaults 10 bp, no fixed cost).

Loss in scenario `s`: `L_s = TC − Σ_i P_i R_{s,i} n_i`. CVaR (Rockafellar–Uryasev):

```
u_s ≥ L_s − α                                (S rows)
CVaR_β = α + 1/(1−β) Σ_s p_s u_s             (p_s = 1/S unless importance sampling is used)
```

Objectives: `max_return` — maximise `E[P&L] = Σ_i P_i μ_i n_i − TC` (model expected simple returns
`μ_i` over the horizon; the probability-weighted scenario mean if none are given) s.t. `CVaR_β ≤ L`
(the absolute CHF risk budget for the period, e.g. "CVaR₉₅ over a week ≤ 40 000"); `min_cvar`
s.t. `E[P&L] ≥ target`; `utility` — maximise `E[P&L] − λ·CVaR_β`.

Constraints: budget `cash₀ + Σ P_i s_i − Σ P_i b_i − TC ≥ 0` (cash may stay uninvested); single-stock
`P_i n_i ≤ w_max W` (10 %, implemented as tight variable bounds `n_i ≤ ⌊w_max W / P_i⌋`); sector
`Σ_{i∈k} P_i n_i ≤ w_k W` (default 25 %, per-sector dict possible); `TC ≤ TC_max`; turnover
`Σ P_i (b_i + s_i) ≤ τ W`; min/max invested fraction; cardinality `Σ y_i ≤ K`, `n_i ≤ M y_i`;
minimum position `P_i n_i ≥ v_min y_i`; `b_i + s_i ≤ M t_i`.

The reported expected P&L (from the model means), VaR and CVaR of the decoded integer portfolio are
recomputed from the weighted scenario set and agree with the solver's values (`tests/test_problem.py`,
including a brute-force optimality check). Note that `α` is only *a* minimiser of the RU function —
it equals the VaR only when the CVaR constraint is tight.

Sizes: `S + N(+2N) + 1` variables, `S + N + K + O(1)` rows, `≈ S·(N+2)` non-zeros; with the default
`N = 100`, `S = 4096` about 420 k non-zeros, solved by cuOpt in ~10 s. Solve times grow faster than
linearly in `S` (HiGHS LP ≈ `S^1.9`, cuOpt ≈ `S^1.2–1.6`): 2 000 → 100 000 scenarios costs ≈ 45× on
cuOpt.

**Tail-row generation** (`quantcore.optimization.rowgen`, `row_generation="auto"` from 8 192
scenarios on): only tail scenarios have `u_s > 0` at the optimum, so the MILP is solved on a working
set (the `2(1−β)S` worst scenarios of a reference portfolio), the losses of all scenarios are
checked, violated ones (`L_s > α*`) are added, and the loop stops when none is violated — the
restricted problem is a relaxation, so the final solution is optimal for the full set (within the
MIP gap). Measured: working sets ≈ 11 % of `S`, 2–3 iterations, 6–7× faster at 32 k–100 k scenarios.

**Exact Gaussian mode** (`scenario_model="gaussian_exact"`, `quantcore.optimization.gaussian_exact`):
with normal horizon returns `R ~ N(m, C)` the CVaR is `TC − m'v + k_β √(v'Cv)` (`v = P∘n`,
`k_β = φ(z_β)/(1−β)`); the mean–CVaR problem then coincides with the mean–variance problem (the
CVaR bound corresponds to an implicit variance bound `σ̄ = (Λ + m'v*)/k_β`). In integer shares this is
a convex MISOCP, solved by outer approximation: tangent cuts `t ≥ TC + (k_β C v_k/σ_k − m)'v` on the
same MILP with an epigraph variable `t` (`cvar`), a new cut at each solution until the true CVaR equals
`t` (20–70 cuts; identical to Gurobi's direct MISOCP solution to the cent). No scenarios, no sampling
error; the reference for the sampling tests.

## 9. Solvers (`quantcore.optimization.solvers`)

Both solvers receive the *same* `MILPModel` (objective, CSR matrix, row and variable bounds,
variable types); `solve_with_fallback` picks between them and the result always names the one that
produced it.

**cuOpt** — NVIDIA's GPU optimisation library, used through its self-hosted REST server
(`nvidia/cuopt`, version 26.08 here) running as the `cuopt` container on the DGX Spark:
`POST /cuopt/request` with a msgpack body, then poll `GET /cuopt/solution/{id}`; `"inf"/"ninf"` for
infinite bounds, `variable_types`, a mandatory `time_limit`. Statuses `Optimal` / `FeasibleFound`
(time limit — check `mip_gap`) / `Infeasible`. Docs: <https://docs.nvidia.com/cuopt/>. It is the
only part of the pipeline that uses the GPU.

**HiGHS** — the open-source LP/MILP/QP solver of ERGO-Code, University of Edinburgh (MIT licence,
<https://highs.dev>, source <https://github.com/ERGO-Code/HiGHS>, reference: Huangfu & Hall,
*Parallelizing the dual revised simplex method*, Math. Prog. Comp. 10, 119–142, 2018). quantcore
neither installs nor starts it separately: **HiGHS ships compiled inside SciPy**, and
`scipy.optimize.milp()` is a thin wrapper around it. It therefore runs **in-process, on the CPU, in
whichever process calls quantcore** — no network, no daemon, no separate container — and its version
follows the installed SciPy:

| where | Python | SciPy | HiGHS | notes |
|---|---|---|---|---|
| `quantcore` container on the Spark (production) | 3.12.14 | 1.18.1 | 1.12.0 (`4f96ee8`) | aarch64 Linux, 4 threads (`OMP/OPENBLAS/MKL_NUM_THREADS=4`); numpy 2.5.3, pandas 3.0.6 |
| developer machine used for this work | 3.11.11 | 1.13.1 | 1.2.0 | arm64 macOS; numpy 1.26.4, pandas 2.1.4 |

The container tracks the newest releases (nothing is pinned beyond the lower bounds in
`pyproject.toml`), so a rebuild can move it a long way from the developer machine — the build of
2026-09-23 jumped it to pandas 3, where the month-end resampling alias is `'ME'` and the old `'M'`
is gone (`quantcore.measurement.returns.MONTH_END` resolves this at import time, and is the reason
monthly frequency works on both).

Verify it in any environment with
`python -c "import scipy.optimize as so, numpy as np, scipy.sparse as sp; print(so.milp(np.array([-3.,-2.]), constraints=so.LinearConstraint(sp.csr_matrix([[2.,3.],[4.,1.]]), [-np.inf]*2, [12.,10.]), integrality=np.ones(2), bounds=so.Bounds([0,0],[np.inf]*2), options={'disp': True}))"`,
which prints the HiGHS banner with its version.

HiGHS is used in four places: `solver="highs"`; as the fallback when cuOpt is unreachable or errors
(`solver="auto"`); for the cut MILPs of the exact Gaussian mode, which are small and numerous; and as
the oracle in the tests (brute-force optimality checks).

**Scenario count and candidate count**: `n_scenarios` defaults to 16 per candidate within
[4 096, 16 384]. The reported figures (expected P&L, CVaR) are stable at any size from ~2 000 on
(sd < 0.1 %), but the *composition* moves between near-equivalent portfolios in proportion to N/S —
3.5 % of wealth between seeds at N = 30, S = 2 048 against 12.5 % at N = 100 — and ~64 per candidate
is needed to hold it near 5 %. The result warns below 8 per candidate.

**Covariance rank**: an estimator with `T_eff` effective observations (rolling: the window; EWM:
`3/(1−λ)`) cannot support more than `T_eff` assets. A 200-day window with N = 300 gives 101 zero
eigenvalues and a portfolio with *negative* estimated variance; `shrinkage="auto"` (the default)
raises the intensity with N/T (0.033 at N = 30 → 0.128 at N = 745) and keeps the smallest eigenvalue
above 2 % of the average variance in every case. `CovarianceModel.conditioning()` reports the
figures and the optimisation warns when the estimate is singular or ill-conditioned for the size of
the candidate set. Long-only with a 10 % cap limits the damage in any case: the optimised portfolio's
out-of-sample/in-sample volatility ratio exceeds the cap-weighted portfolio's by only 0.89–1.18×
across N = 30…745, with or without shrinkage.

**cuOpt transport**: the model travels as one msgpack body (9.9 M non-zeros = 108 MB at N = 300,
S = 32 768), accepted at ~10 MB/s, so transport is not the constraint. The client uses a payload-scaled
timeout and a *fresh connection per attempt* — retrying on a pooled connection appends a second body and
the server rejects it with "unpack received extra data" (422) — and `solver="auto"` falls back to HiGHS
on failure. cuOpt's concurrent method also aborts on some CVaR-objective models; the client retries once
with its dual simplex (`method=2`). With row generation only the working set is sent.

**Which models go to cuOpt.** cuOpt 26.08 aborts on one specific combination: the CVaR in the
*objective* (`objective="min_cvar"`) **and** probability-weighted scenarios. The cause is the scaling
of the objective vector. In the min-CVaR formulation the coefficient on each excess `u_s` is
`p_s/(1−β)`, so the scenario probabilities *are* the objective coefficients; importance sampling
spreads them over about six orders of magnitude. Measured on a 4 137-variable / 4 106-row model
(40 candidates, 4 096 scenarios, β = 99 %):

| scenarios | objective `|c|` range | ratio | cuOpt |
|---|---|---|---|
| importance-sampled | 1.1e-06 … 1.0 | 8.8e5 | aborts in ~5 s |
| uniform (`importance_sampling="off"`) | 2.4e-02 … 1.0 | 4.1e1 | **Optimal in 3.8 s**, 5 runs out of 5, matching HiGHS to 7e-5 |

Relaxing the integrality does not help — the pure LP aborts too, so it is the LP and not the branch
and bound — and multiplying the objective by a constant does not either, since that leaves the ratio
unchanged. In the `max_return` formulation the same probabilities sit in a *constraint row* instead
(the objective is then just 40 expected returns, range 1e3), where cuOpt's row scaling copes, and it
solves at either weighting. Its methods are not interchangeable either: `method=1` (PDLP) returned the
weighted optimum once in 4.9 s and then aborted three times on the identical model, while `method=2`
and `method=3` always abort.

So `solver="auto"` diverts exactly that combination to HiGHS (1.7 s for that model) and the result
says why — otherwise a green cuOpt health pill next to `solver: "highs"` looks like the GPU being
ignored. Everything else keeps the GPU path. An explicit `solver="cuopt"` is always honoured and
returns cuOpt's own answer, abort included; an abort under `"auto"` is retried down
`CUOPT_RETRY_METHODS = (1, 2)` before the fallback.

**cuOpt vs HiGHS, both on the Spark** (end-to-end through the API: scenario generation, model build,
solve, decode; GBM model, weekly, β = 0.95, no tilt, 600 s solver time limit):

| case | cuOpt (GPU) | HiGHS (CPU, 4 threads) |
|---|---|---|
| N = 30, S = 8 192 | 15.6 s | 14.5 s |
| N = 30, S = 32 768 | 171.4 s | 171.0 s |
| N = 100, S = 4 096 | 10.0 s | 8.9 s |
| N = 100, S = 32 768 | 606.9 s (time limit, `FeasibleFound`) | 350.7 s (`Optimal`) |
| N = 300, S = 8 192 (row generation) | 139 s | 64 s |

HiGHS is at least as fast in every case measured and 2× faster on the two largest. That is not a
statement about GPUs in general: cuOpt's published speedups are for *LP* with its PDLP first-order
method, while [its own documentation](https://docs.nvidia.com/cuopt/user-guide/latest/milp-features.html)
states that for MILP the GPU runs primal heuristics and **branch and bound runs on the CPU**, that the
MIP solver is beta, and that "proving feasible solutions optimal remains under active development".
Our models have an easy LP relaxation and a small integrality gap, so the work sits almost entirely in
proving optimality — the part cuOpt does not accelerate — while HiGHS' dual simplex warm-starts across
branch-and-bound nodes. HiGHS is also a strong reference point: in the [Mittelmann MILP benchmark](https://plato.asu.edu/ftp/milp.html)
(7 July 2026) HiGHS 1.15.1 solves 179 of 240 instances at a scaled shifted geometric mean of 5.44,
against 1.72 for COPT 8.0.3 and 9.93 for SCIP. Caveats: cuOpt's `time_limit` covers only its internal
solve, it ignores the MIP-gap settings sent over REST (so it may prove optimality to a tighter
tolerance), and this is one problem family on one machine, not a benchmark.

**Reproducibility**: the scenario set depends only on the seed (`scenario_hash` in every result), and
the hash is identical for two runs **in the same environment**. Across environments it is not, and the
reason is the fingerprint rather than the mathematics. Measured between the developer Mac (numpy 1.26.4,
pandas 2.1.4, scipy 1.13.1) and the Spark container (numpy 2.5.3, pandas 3.0.6, scipy 1.18.1) on the
same 100 stocks and the same stored prices: scrambled Sobol' points and the numpy RNG stream are
bit-identical, the estimated covariance agrees to `max |dC| = 1.7e-18` (entries ~5e-4), its
eigenvalues to 1.4e-17 and the columns of `pca_root` to 4e-13 relative — no sign flips, no rotation.
The scenario matrices therefore agree to 13–14 significant digits. But `scenario_hash` fingerprints
`np.round(scenarios, 12)`, and on a 4096 × 100 matrix of weekly log-returns a *relative* perturbation
of 1e-15 already pushes ~3 entries of 409 600 across a rounding boundary (1e-14 pushes ~47), which
changes the digest. The figures are unaffected: the same request optimised on both machines gave
E[P&L] 1 888.65 against 1 888.76 and CVaR 39 999.71 against 39 999.82 on a CHF 1 M portfolio.

The portable identity is the **scenario key**, not the hash: it is built from the request and the
model description, so two environments that agree on the inputs produce the same `key` (verified:
`7162939ddf75ff74` on both machines for the run above) even where the digest of the realised matrix
differs. Read `scenario_hash` as "these two results used bit-identical scenarios" — a statement about
one machine — and `scenario_key` as "these two results were asked for the same scenarios".
HiGHS' branch-and-bound is deterministic → identical shares for the same inputs *and the same HiGHS
build* (`reproducibility="exact"`); across SciPy versions (1.2.0 on the developer machine against
1.12.0 in the container) the search can end on a different optimal vertex among ties, with the same
objective. cuOpt 26.08
ignores the MIP-gap settings sent over REST and stops within its default relative gap of 1e-4: the
objective, expected P&L and CVaR of two runs agree to that gap (≈ 0.2 CHF/week on CHF 1 M) while
near-tied share counts can differ by a few shares (`reproducibility="within_gap"`). The Gaussian exact
mode runs on HiGHS by default.

## 10. Default parameters

| Parameter | Default | Where |
|---|---|---|
| history start / required first price | 2004-12-01 / ≤ 2005-01-31, coverage ≥ 90 % | `Settings` |
| base currency, risk-free | CHF, 0 % | `Settings` |
| measurement window / λ | 200 periods / 0.94 | `RiskService` |
| optimisation covariance | EWM λ = 0.99 (≈ 100-day effective window) | `OptimizationRequest.covariance` |
| scenarios | empirical innovations, 4 096 scrambled Sobol' draws on the market-factor ranks, moment matched, arithmetic drift, last 2 500 days, diagonal re-standardisation, seed 42 | `OptimizationRequest` |
| importance sampling / row generation | `auto` = market tilt (blocks at z₀.₉₅, z₀.₉₉, 25 % each); row generation from 8 192 scenarios | `OptimizationRequest` |
| MIP gap | 1e-4 relative (HiGHS; cuOpt uses its own default) | `Settings.mip_rel_gap` |
| Black–Litterman | 10 years weekly, τ = 0.05, sentiment strength 0.25 | `OptimizationRequest`, `view_from_sentiment` |
| constraints | 10 % per stock, 25 % per sector, no shorts, 10 bp costs | `OptimizationConstraints`, `TransactionCosts` |
| CVaR limit | explicit CHF amount or fraction (5 % of wealth if nothing is given) | `OptimizationRequest` |
| cuOpt time limit | 120 s (`QUANTCORE_CUOPT_TIME_LIMIT`) | `Settings` |