% !TEX encoding = UTF-8 Unicode % --------------------------------------------------------------------------------------------- % quantcore -- documentation of the portfolio optimisation model % From the Rockafellar-Uryasev mean-CVaR LP in portfolio weights to the mixed-integer % programme in whole shares with rebalancing and transaction costs, as implemented in % quantcore/optimization/problem.py. % Build: latexmk -pdf docs/optimization_model.tex % --------------------------------------------------------------------------------------------- \documentclass[11pt,a4paper]{article} \usepackage[utf8]{inputenc} \usepackage[T1]{fontenc} \usepackage[british]{babel} \usepackage{lmodern} \usepackage{microtype} \usepackage[margin=2.4cm]{geometry} \usepackage{amsmath,amssymb,amsthm,mathtools} \usepackage{booktabs,array,longtable} \usepackage{enumitem} \usepackage{xcolor} \usepackage[colorlinks=true,linkcolor=blue!50!black,citecolor=blue!50!black,urlcolor=blue!50!black]{hyperref} \newtheorem{theorem}{Theorem}[section] \newtheorem{proposition}[theorem]{Proposition} \newtheorem{lemma}[theorem]{Lemma} \theoremstyle{definition} \newtheorem{definition}[theorem]{Definition} \newtheorem{example}[theorem]{Example} \theoremstyle{remark} \newtheorem{remark}[theorem]{Remark} \DeclareMathOperator{\CVaR}{CVaR} \DeclareMathOperator{\VaR}{VaR} \DeclareMathOperator{\E}{\mathbb{E}} \DeclareMathOperator*{\argmin}{arg\,min} \newcommand{\R}{\mathbb{R}} \newcommand{\Z}{\mathbb{Z}} \newcommand{\pos}[1]{\left[#1\right]^{+}} \newcommand{\code}[1]{\texttt{#1}} \newcommand{\TC}{\mathit{TC}} \setlength{\parskip}{4pt} \setlength{\parindent}{0pt} \title{Mean--CVaR Portfolio Optimisation in Whole Shares\\[1ex] \large From the Rockafellar--Uryasev linear programme in portfolio weights to the mixed-integer programme with rebalancing and transaction costs implemented in \code{quantcore}} \author{AISE502 -- AI-Augmented Portfolio Intelligence Platform\\ \small Deterministic quant core, module \code{quantcore.optimization}} \date{September 2026} \begin{document} \maketitle \begin{abstract} \noindent The portfolio optimisation service of the platform selects a portfolio of stocks such that the expected profit over a horizon is maximal while the Conditional Value-at-Risk (expected shortfall) of the loss over that horizon stays below an absolute risk budget in CHF. Positions are whole numbers of shares, sector and single-stock limits apply, short positions are not allowed, and an existing portfolio can be rebalanced at proportional and fixed transaction costs. This note states the classical scenario-based mean--CVaR linear programme of Rockafellar and Uryasev in portfolio weights, and derives from it, step by step, the mixed-integer linear programme (MILP) that \code{quantcore/optimization/problem.py} builds and that is solved on NVIDIA cuOpt (with HiGHS as fallback). Every constraint of the implementation is traced back to its origin, a worked example shows why rounding a continuous solution is not an option, a section on scenario sets quantifies the sampling error of the approximation and derives the variance-reduction, importance-sampling and row-generation techniques in the code together with the exact (mean--variance-equivalent) solution of the Gaussian case, and the last section lists the complete model with the variable and row names used in the code. \end{abstract} \tableofcontents % ============================================================================================ \section{Setting and notation} % ============================================================================================ There are $N$ candidate stocks $i=1,\dots,N$ with current prices $P_i>0$ in the base currency (CHF). The investor has cash $c^0\ge 0$ and, in the rebalancing case, existing holdings $n^0_i\in\Z_{\ge 0}$ shares; the total wealth is \begin{equation} W \;=\; c^0 + \sum_{i=1}^N P_i\, n^0_i . \label{eq:wealth} \end{equation} The uncertainty is described by $S$ equiprobable scenarios of simple returns over the horizon, \begin{equation} R_{s,i} \;=\; e^{\,r_{s,i}} - 1 , \qquad s = 1,\dots,S, \label{eq:scenarios} \end{equation} where the log-returns $r_{s,i}$ (in CHF) are simulated by one of the stock-price models of the platform (geometric Brownian motion or empirical innovations) with the expected returns of the Black--Litterman model as drift; see \code{docs/METHODOLOGY.md}, Sections~6--7. Each scenario has probability $1/S$. The confidence level of the risk measure is $\beta\in(0,1)$ (0.95 or 0.99). Throughout, \emph{losses are positive numbers}: a loss of $L$ means the portfolio is worth $L$ less than at the start of the horizon. Vectors are columns, $\pos{z}=\max(z,0)$, and $\mathbf{1}$ is the vector of ones. % ============================================================================================ \section{Conditional Value-at-Risk and its linear representation} % ============================================================================================ Let $L(x,Y)$ be the loss of a decision $x\in X\subseteq\R^N$ under the random outcome $Y$, with distribution function $\Psi(x,\zeta)=\mathbb{P}\{L(x,Y)\le\zeta\}$. \begin{definition}[VaR and CVaR] The \emph{Value-at-Risk} at level $\beta$ is the $\beta$-quantile of the loss, $\VaR_\beta(x)=\min\{\zeta\in\R:\Psi(x,\zeta)\ge\beta\}$. The \emph{Conditional Value-at-Risk} (expected shortfall) is the mean of the $\beta$-tail of the loss distribution, \begin{equation} \CVaR_\beta(x) \;=\; \frac{1}{1-\beta}\int_{\beta}^{1} \VaR_{\gamma}(x)\,d\gamma , \label{eq:cvar-def} \end{equation} which for a continuous loss distribution equals $\E\big[L(x,Y)\,\big|\,L(x,Y)\ge\VaR_\beta(x)\big]$. For distributions with atoms (in particular for a finite set of scenarios) the definition \eqref{eq:cvar-def} splits the atom at the VaR appropriately \cite{RockafellarUryasev2002, AcerbiTasche2002}. \end{definition} CVaR is a coherent risk measure (monotone, sub-additive, positively homogeneous, translation invariant), which VaR is not. Its value for optimisation purposes rests on the following representation. \begin{theorem}[Rockafellar--Uryasev \cite{RockafellarUryasev2000,RockafellarUryasev2002}] \label{thm:ru} Define the auxiliary function \begin{equation} F_\beta(x,\zeta) \;=\; \zeta + \frac{1}{1-\beta}\,\E\Big[\pos{L(x,Y)-\zeta}\Big]. \label{eq:F} \end{equation} Then $F_\beta(x,\cdot)$ is finite and convex, and \begin{equation} \CVaR_\beta(x) \;=\; \min_{\zeta\in\R} F_\beta(x,\zeta), \qquad \VaR_\beta(x) \;=\; \text{left endpoint of } \argmin_{\zeta} F_\beta(x,\zeta). \label{eq:ru} \end{equation} If $L(x,Y)$ is convex in $x$, then $F_\beta$ is jointly convex in $(x,\zeta)$ and \begin{equation} \min_{x\in X}\CVaR_\beta(x) \;=\; \min_{(x,\zeta)\in X\times\R} F_\beta(x,\zeta), \label{eq:ru-joint} \end{equation} with $(x^*,\zeta^*)$ minimising the right-hand side if and only if $x^*$ minimises the left-hand side and $\zeta^*\in\argmin_\zeta F_\beta(x^*,\zeta)$. \end{theorem} The theorem removes the need to know the VaR before computing the CVaR: both come out of one minimisation over the scalar $\zeta$, and $\zeta$ can simply be added to the decision variables. A second consequence is used for risk \emph{constraints} \cite{KrokhmalPalmquistUryasev2002}: for any budget $\omega$, \begin{equation} \CVaR_\beta(x)\le\omega \quad\Longleftrightarrow\quad \exists\,\zeta:\; F_\beta(x,\zeta)\le\omega , \label{eq:cvar-constraint} \end{equation} because the minimum over $\zeta$ is at most $\omega$ exactly when some $\zeta$ attains at most $\omega$. \paragraph{Scenario approximation and linearisation.} With $S$ equiprobable scenarios the expectation in \eqref{eq:F} becomes an average, \begin{equation} \tilde F_\beta(x,\zeta) \;=\; \zeta + \frac{1}{(1-\beta)S}\sum_{s=1}^{S} \pos{L_s(x)-\zeta}, \qquad L_s(x)=L(x,Y_s). \label{eq:F-scen} \end{equation} The positive parts are linearised with auxiliary variables $u_s\ge 0$: since $u_s\ge L_s(x)-\zeta$ and $u_s\ge 0$ imply $u_s\ge\pos{L_s(x)-\zeta}$, and the objective/constraint pushes $u_s$ down, \begin{equation} \min_\zeta \tilde F_\beta(x,\zeta) \;=\; \min_{\zeta,\,u}\Big\{\zeta + \frac{1}{(1-\beta)S}\sum_{s} u_s \;:\; u_s \ge L_s(x)-\zeta,\; u_s\ge 0,\; s=1,\dots,S\Big\}. \label{eq:lin} \end{equation} Whenever $L_s(x)$ is linear in $x$, \eqref{eq:lin} is a linear programme in $(x,\zeta,u)$. \begin{remark}[The discrete tail and the meaning of $\zeta$] \label{rem:discrete} For $S$ scenarios and $k=(1-\beta)S$, the minimum of \eqref{eq:F-scen} equals the average of the $k$ largest losses, where a fractional $k$ gives the $(\lfloor k\rfloor+1)$-th largest loss the fractional weight $k-\lfloor k\rfloor$. This is exactly what \code{quantcore.measurement.risk.scenario\_cvar} computes, and it is the value the LP returns. When $k$ is an integer, $\tilde F_\beta(x,\cdot)$ is flat between the $(k{+}1)$-th and the $k$-th largest loss, so the minimiser $\zeta$ is \emph{not unique}: any point of that interval is optimal, and the VaR is its left endpoint. Consequently the $\zeta$ (called $\alpha$ in the code) reported by a solver need not equal the VaR, in particular when the CVaR constraint is not tight; the implementation therefore recomputes VaR and CVaR of the decoded portfolio from the scenarios rather than reading them off the LP variables. \end{remark} % ============================================================================================ \section{The classical model in portfolio weights} % ============================================================================================ The textbook formulation \cite{RockafellarUryasev2000,KrokhmalPalmquistUryasev2002} describes a portfolio by weights $w\in\R^N$, the fractions of wealth invested in each stock. Its loss over the horizon, relative to the initial wealth, is minus the portfolio return, \begin{equation} L_s(w) \;=\; -\sum_{i=1}^N R_{s,i}\,w_i \;=\; -R_s^{\top} w , \label{eq:loss-w} \end{equation} which is linear in $w$, and the expected return is $\bar R^{\top}w$ with $\bar R_i=\frac1S\sum_s R_{s,i}$. With a risk budget $\omega$ expressed as a fraction of wealth (e.g.\ ``the 95\,\% CVaR over one week may not exceed 4\,\% of the portfolio value''), the \emph{mean--CVaR linear programme} reads \begin{subequations} \label{eq:lp-w} \begin{alignat}{3} \max_{w,\zeta,u}\quad & \bar R^{\top} w \label{eq:lp-w-obj}\\ \text{s.t.}\quad & u_s \;\ge\; -R_s^{\top} w - \zeta, \quad u_s\ge 0, & \quad & s=1,\dots,S, \label{eq:lp-w-u}\\ & \zeta + \frac{1}{(1-\beta)S}\sum_{s=1}^S u_s \;\le\; \omega, \label{eq:lp-w-cvar}\\ & \mathbf{1}^{\top} w \;\le\; 1, \qquad w \;\ge\; 0, \label{eq:lp-w-budget}\\ & w_i \;\le\; \bar w, & & i=1,\dots,N, \label{eq:lp-w-single}\\ & \sum_{i\in G_k} w_i \;\le\; \bar w_k, & & k=1,\dots,K. \label{eq:lp-w-sector} \end{alignat} \end{subequations} Here \eqref{eq:lp-w-u}--\eqref{eq:lp-w-cvar} are \eqref{eq:lin} and \eqref{eq:cvar-constraint} combined, \eqref{eq:lp-w-budget} is the budget with a non-negative cash residual $1-\mathbf{1}^\top w$ (earning nothing over the horizon) and the no-short-sale condition, \eqref{eq:lp-w-single} caps each stock at the fraction $\bar w$ (10\,\% in the platform) and \eqref{eq:lp-w-sector} caps every GICS sector $G_k$ at $\bar w_k$. \paragraph{Three equivalent objectives.} Krokhmal, Palmquist and Uryasev \cite{KrokhmalPalmquistUryasev2002} show that the three problems \begin{equation} \text{(a)}\;\max\,\bar R^{\top}w \text{ s.t. } \CVaR_\beta\le\omega, \quad \text{(b)}\;\min\,\CVaR_\beta \text{ s.t. } \bar R^{\top}w\ge\rho, \quad \text{(c)}\;\max\,\bar R^{\top}w-\lambda\CVaR_\beta \label{eq:three} \end{equation} trace out the same efficient frontier as $\omega$, $\rho$ and $\lambda\ge 0$ vary. All three are linear after the substitution \eqref{eq:lin}; the implementation offers them as the objectives \code{max\_return}, \code{min\_cvar} and \code{utility}. The platform's primary mode is (a), because the risk-management specification is an \emph{absolute} loss level for a period. \paragraph{Properties.} Problem \eqref{eq:lp-w} has $N+1+S$ variables and $S+1+N+K+1$ constraints (the weight caps are variable bounds). Its constraint matrix is dense in the $S$ scenario rows: $S(N+2)$ non-zeros. It is a linear programme, so every local optimum is global, the solution is a vertex, and it can be solved for thousands of scenarios and hundreds of assets in seconds. % ============================================================================================ \section{From weights to numbers of shares} \label{sec:shares} % ============================================================================================ Weights are a convenient abstraction, but a real portfolio consists of whole shares: one cannot buy 2.66 shares of a stock, and for a stock priced at CHF~2\,000 the difference between 2 and 3 shares is 10\,\% of a CHF~10\,000 portfolio. The platform therefore optimises directly over the integer numbers of shares. \subsection{Money instead of fractions} Let $n_i\in\Z_{\ge0}$ be the number of shares of stock $i$ held over the horizon. The value of the position is $P_i n_i$, and the weight of the classical model is recovered as \begin{equation} w_i \;=\; \frac{P_i\,n_i}{W}. \label{eq:w-n} \end{equation} Substituting \eqref{eq:w-n} into the relative loss \eqref{eq:loss-w} and multiplying by $W$ turns it into a loss \emph{in money}: \begin{equation} W L_s(w) \;=\; -\sum_i R_{s,i}\,P_i\,n_i \;=\; -\sum_i g_{s,i}\,n_i, \qquad g_{s,i} := P_i R_{s,i}, \label{eq:loss-n} \end{equation} where $g_{s,i}$ is the CHF gain of one share of stock $i$ in scenario $s$. Because CVaR is positively homogeneous, $\CVaR_\beta(W L)=W\CVaR_\beta(L)$: the CVaR of the money loss is the money amount the risk manager specifies (``the 95\,\% CVaR over a week may not exceed CHF~40\,000''), and the relative budget $\omega$ of \eqref{eq:lp-w-cvar} becomes the absolute budget $\Lambda=\omega W$. The expected return $\bar R^\top w$ becomes the expected profit \begin{equation} \E[\text{P\&L}] \;=\; \sum_i \bar g_i\,n_i, \qquad \bar g_i=\frac1S\sum_s g_{s,i}=P_i\bar R_i . \label{eq:epnl} \end{equation} Both the objective and the risk constraint use the \emph{same} scenario set, so the trade-off they describe is internally consistent (the drift of the scenarios is the expected-return model). \subsection{Budget, cash and the constraints} Multiplying each constraint of \eqref{eq:lp-w} by $W$ and substituting \eqref{eq:w-n} gives: \begin{itemize}[nosep] \item budget: $\sum_i P_i n_i \le W$ --- the purchases may not exceed the available cash, what is not invested stays in cash; \item single stock: $P_i n_i\le \bar w W$; sector: $\sum_{i\in G_k}P_i n_i\le \bar w_k W$; \item no short sales: $n_i\ge 0$. \end{itemize} The single-stock cap deserves a closer look. As a linear constraint it is a row of the model, but for an integer variable it is equivalent to the \emph{variable bound} \begin{equation} n_i \;\le\; \bar n_i := \Big\lfloor \frac{\bar w\,W}{P_i} \Big\rfloor , \label{eq:nbar} \end{equation} which is strictly tighter than the linear row whenever $\bar w W/P_i$ is not an integer: the rounding down is free information that a MILP solver exploits (tighter LP relaxation, smaller branch-and-bound tree). The implementation uses \eqref{eq:nbar} with $\bar w$ replaced by $\min(\bar w, f_{\max})$ if a maximum invested fraction $f_{\max}<1$ is given, and never leaves the cap to a row. \subsection{The mean--CVaR MILP in shares} Putting the pieces together, the new-portfolio problem (all wealth in cash, $n^0=0$, $W=c^0$) reads \begin{subequations} \label{eq:milp-n} \begin{alignat}{3} \max_{n,\zeta,u}\quad & \sum_{i} \bar g_i\, n_i \label{eq:milp-n-obj}\\ \text{s.t.}\quad & \sum_i g_{s,i}\, n_i + \zeta + u_s \;\ge\; 0, \quad u_s\ge 0, & \quad & s=1,\dots,S, \label{eq:milp-n-u}\\ & \zeta + \frac{1}{(1-\beta)S}\sum_{s} u_s \;\le\; \Lambda, \label{eq:milp-n-cvar}\\ & \sum_i P_i\, n_i \;\le\; c^0, \label{eq:milp-n-budget}\\ & \sum_{i\in G_k} P_i\, n_i \;\le\; \bar w_k\, W, & & k=1,\dots,K, \label{eq:milp-n-sector}\\ & n_i \in \{0,1,\dots,\bar n_i\}, & & i=1,\dots,N, \label{eq:milp-n-int}\\ & \zeta\in\R . \nonumber \end{alignat} \end{subequations} Row \eqref{eq:milp-n-u} is \eqref{eq:lp-w-u} written as ``$\ge 0$'' with the money loss \eqref{eq:loss-n} inserted ($u_s\ge L_s-\zeta \Leftrightarrow \sum_i g_{s,i}n_i+\zeta+u_s\ge0$). Apart from the integrality \eqref{eq:milp-n-int} the problem is the LP \eqref{eq:lp-w} scaled by $W$; with the integrality it is a mixed-integer linear programme, and the LP obtained by dropping \eqref{eq:milp-n-int} is its \emph{LP relaxation}. \subsection{Why the integrality cannot be handled by rounding} The optimum of the LP relaxation is a vertex, and its coordinates $n_i$ are generically fractional. Rounding them is not a substitute for solving the MILP: \begin{itemize}[nosep] \item rounding \emph{up} (or to the nearest integer) can violate the budget, a cap, or the CVaR constraint --- the risk budget is exceeded; \item rounding \emph{down} is always feasible (all constraints are ``$\le$'' in $n$, and the scenario rows are satisfied by re-optimising $\zeta,u$) but can be far from optimal, because the ``lost'' fractions can often be re-invested in other, cheaper shares; \item the true integer optimum can have a different \emph{structure} than the relaxed one, not just rounded coordinates. \end{itemize} The gap is largest when $W/P_i$ is small, i.e.\ for expensive stocks or small portfolios; for a CHF~1\,000\,000 portfolio and stocks below CHF~1\,000 it is negligible, but the 10\,\% cap keeps each position to at most $\bar w W/P_i$ shares, and a few stocks in the universe trade above CHF~10\,000 (Lindt \& Spr\"ungli partizipation certificates above CHF~100\,000). \begin{example}[Two stocks, eight scenarios] \label{ex:rounding} Prices $P=(572,\ 2208)$, cash $W=c^0=10\,000$, $\beta=0.75$ (so $k=(1-\beta)S=2$: the CVaR is the mean of the two worst scenarios), risk budget $\Lambda=250$, caps $\bar w=0.7$, no costs. The simple returns of the eight scenarios are \[ R=\begin{pmatrix} 0.005&-0.082\\ 0.014&-0.052\\ -0.069&-0.020\\ -0.002&0.055\\ -0.009&0.083\\ 0.048&0.064\\ 0.044&0.028\\ -0.011&0.042 \end{pmatrix}, \qquad \bar R=(0.0025,\ 0.01475). \] Stock 2 has six times the expected return of stock 1 but a heavier left tail. \begin{center} \small \begin{tabular}{@{}lcccccl@{}} \toprule portfolio & $n_1$ & $n_2$ & invested & $\E[\text{P\&L}]$ & $\CVaR_{0.75}$ & feasible \\ \midrule LP relaxation ($n\in\R^2_{\ge0}$) & 2.6605 & 1.7876 & 5\,468.90 & 62.02 & 250.00 & yes, CVaR binding \\ nearest rounding & 3 & 2 & 6\,132 & 69.43 & 280.13 & \textbf{no}: $\CVaR>\Lambda$ \\ rounding down & 2 & 1 & 3\,352 & 35.43 & 149.22 & yes, $-43\,\%$ vs.\ LP \\ MILP optimum & 7 & 1 & 6\,212 & 42.58 & 240.74 & yes \\ \bottomrule \end{tabular} \end{center} Rounding to the nearest integers breaks the risk budget (the CVaR limit, not the cash: 6\,132 of 10\,000 are invested). The relaxed solution buys 1.79 shares of the attractive stock~2; the second whole share of it is too large a block of risk for the budget, so the integer optimum keeps one share of stock~2 and fills the remaining risk budget with \emph{seven} shares of stock~1 --- a portfolio that no rounding of $(2.66, 1.79)$ produces. The enumeration of all feasible integer portfolios confirms $(7,1)$ as the optimum; the integrality gap $62.02-42.58$ is 31\,\% of the relaxed value. (The numbers were computed with \code{MeanCVaRProblem} and HiGHS and are pinned by \code{tests/test\_problem.py::test\_documentation\_example\_rounding\_fails}.) \end{example} % ============================================================================================ \section{Rebalancing an existing portfolio with transaction costs} % ============================================================================================ When holdings $n^0\ne 0$ exist, the decision is no longer ``which shares to buy'' but ``which shares to buy and which to sell'', and every trade costs money. Weights cannot express this naturally (a weight change of $+2\,\%$ says nothing about the number of shares traded, and the cost depends on the traded value), whereas in share numbers it is immediate. \subsection{Trades, linking and the cash balance} Introduce the integer trade variables $b_i\ge0$ (shares bought) and $s_i\ge0$ (shares sold), linked to the final position by \begin{equation} n_i \;=\; n^0_i + b_i - s_i, \qquad s_i\le n^0_i \;(\text{implied by } n_i\ge0), \label{eq:link} \end{equation} and the transaction cost with proportional rates $c^b, c^s$ (10 basis points by default) and an optional fixed cost $F$ per traded stock, \begin{equation} \TC(b,s,t) \;=\; \sum_i \big(c^b P_i\, b_i + c^s P_i\, s_i\big) + F\sum_i t_i, \qquad b_i + s_i \le M_i\, t_i,\; t_i\in\{0,1\}, \label{eq:tc} \end{equation} where $t_i$ indicates whether stock $i$ is traded at all and $M_i=\max(\bar n_i, n^0_i)+1$ is a valid big-M (no trade can exceed it). The fixed cost makes small trades unattractive, which is what it should do; it is the only place where binary variables enter the cost model. The budget \eqref{eq:milp-n-budget} becomes a \emph{cash balance}: sales bring cash in, purchases and costs take it out, and the cash at the end of the rebalancing may not be negative, \begin{equation} c^{\text{end}} \;=\; c^0 + \sum_i P_i s_i - \sum_i P_i b_i - \TC \;\ge\; 0 \quad\Longleftrightarrow\quad \sum_i P_i b_i - \sum_i P_i s_i + \TC \;\le\; c^0 . \label{eq:cash} \end{equation} \begin{remark}[Buying and selling the same stock] With positive proportional costs an optimal solution never has $b_i>0$ and $s_i>0$ simultaneously: reducing both by one share leaves $n_i$ unchanged and saves $(c^b+c^s)P_i$ in cost and cash. With zero costs the decomposition of a net trade into $b_i$ and $s_i$ is not unique but harmless. The decoded solution reports the net trade $b_i-s_i$. \end{remark} \subsection{Transaction costs enter the loss} The costs are paid at the start of the horizon and are therefore a deterministic part of the loss relative to the initial wealth $W$ of \eqref{eq:wealth}: \begin{equation} L_s(n,b,s,t) \;=\; \TC(b,s,t) - \sum_i g_{s,i}\, n_i . \label{eq:loss-reb} \end{equation} By translation invariance, $\CVaR_\beta(\TC-G)=\TC+\CVaR_\beta(-G)$: the costs shift the whole loss distribution and consume part of the risk budget, exactly as they reduce the expected profit, \begin{equation} \E[\text{P\&L}] \;=\; \sum_i \bar g_i\, n_i - \TC(b,s,t). \label{eq:epnl-reb} \end{equation} The scenario rows \eqref{eq:milp-n-u} therefore become \begin{equation} \sum_i g_{s,i}\, n_i - \TC(b,s,t) + \zeta + u_s \;\ge\; 0, \qquad s=1,\dots,S . \label{eq:u-reb} \end{equation} \subsection{Limits on the rebalancing itself} Two constraints of the specification concern the trading rather than the resulting portfolio: a cap on the total transaction cost and a cap on the turnover (the traded value as a fraction of wealth), \begin{equation} \TC(b,s,t)\;\le\;\TC_{\max}, \qquad \sum_i P_i\,(b_i+s_i)\;\le\;\tau\,W . \label{eq:turnover} \end{equation} Both are optional. The single-stock and sector caps apply to the \emph{final} positions $n=n^0+b-s$: a legacy position above the cap has to be sold down to it, and if a turnover cap prevents that, the problem is (correctly) infeasible. \subsection{The new portfolio as a special case} With $n^0=0$ nothing can be sold, so $s=0$ and $b=n$; the linking rows \eqref{eq:link} become trivial and are dropped, the costs become $\TC=c^b\sum_i P_i n_i+F\sum_i t_i$, and the cash balance \eqref{eq:cash} reduces to $\sum_i P_i(1+c^b)\,n_i+F\sum_i t_i\le c^0$. This is \eqref{eq:milp-n} with proportional costs, i.e.\ mode \code{new} of the implementation. \subsection{Optional cardinality and minimum position size} A portfolio of 30 stocks is easier to monitor than one of 120 with many tiny positions; the integer model can express both a maximum number of positions and a minimum value per open position with one binary $y_i$ per stock: \begin{equation} n_i \;\le\; \max(\bar n_i,1)\, y_i, \qquad P_i n_i \;\ge\; v_{\min}\, y_i, \qquad \sum_i y_i\le K_{\max}, \qquad y_i\in\{0,1\}. \label{eq:card} \end{equation} These constraints are optional (\code{max\_positions}, \code{min\_position\_value}); they make the problem combinatorially harder and are usually not needed --- the CVaR constraint and the integrality already produce concentrated portfolios. % ============================================================================================ \section{Scenario sets: sampling error, variance reduction and large sets} \label{sec:sampling} % ============================================================================================ The MILP sees the return distribution only through the scenario matrix $g_{s,i}=P_iR_{s,i}$ and the probabilities of its rows. It is a \emph{sample-average approximation} (SAA) of the problem one would like to solve --- the one with the true distribution of the horizon returns \cite{ShapiroDentchevaRuszczynski2009} --- and its solution inherits the sampling error of the scenario set. This section quantifies that error on the real universe, derives the three remedies implemented in \code{quantcore.models.sampling} and \code{quantcore.optimization}, and closes with the Gaussian special case, in which the exact solution is available and the mean--CVaR problem coincides with the mean--variance problem. \subsection{How noisy is a plain Monte-Carlo set?} \label{sec:mc-noise} Let the scenarios be $S$ independent draws from the model with mean $\mu$ (per asset, over the horizon) and covariance $\Sigma$. The sample mean of asset $i$ has standard error $\sigma_i/\sqrt S$. With $S=2\,000$ and weekly volatilities of 3--11\,\% this is 0.07--0.25\,\% per week, i.e.\ 3--13\,\% per year --- the same order as the expected returns themselves (5--12\,\% p.a.\ from the Black--Litterman model). The optimiser cannot tell signal from noise: it overweights the assets whose sample mean happens to be high, the classical ``error maximisation'' of \cite{Michaud1989}. The tail is only the second problem: the CVaR at 95\,\% is estimated from the $(1-\beta)S=100$ worst scenarios, with a relative standard error of about $0.37/\sqrt{(1-\beta)S}\approx 2\,\%$ for a Gaussian tail. Table~\ref{tab:noise} measures the consequences on the 30 largest stocks of the universe (CHF 1 million, weekly horizon, $\CVaR_{0.95}\le 5\,\%$ of wealth, eight seeds per method). The reference is the exact optimum of Section~\ref{sec:gaussian}, so ``true'' figures are exact. The comparison is made in portfolio weights (the continuous relaxation) so that integrality neither hides nor causes noise; the integer results are the same to within a percent. \begin{table}[ht] \centering \footnotesize \setlength{\tabcolsep}{4pt} \begin{tabular}{@{}lrrrrr@{}} \toprule $S=2\,048$, $\beta=0.95$ & plain MC & antithetic & MC + mm & Sobol' & Sobol' + mm \\ \midrule in-sample $\E[\text{P\&L}]$, mean $\pm$ sd & $2\,791\pm552$ & $2\,228\pm19$ & $2\,221\pm10$ & $2\,221\pm8$ & $2\,220\pm5$ \\ true $\E[\text{P\&L}]$ of the chosen portfolio & $1\,991\pm119$ & $2\,228\pm19$ & $2\,221\pm10$ & $2\,220\pm5$ & $2\,220\pm5$ \\ loss versus the optimum ($2\,219$) & $-228$ & $-9$ & $-2$ & $-0.6$ & $-0.4$ \\ true CVaR over the limit, mean $\pm$ sd & $-0.5\pm5.6\,\%$ & $+1.8\pm2.8\,\%$ & $+0.6\pm1.7\,\%$ & $+0.2\pm0.7\,\%$ & $+0.15\pm0.7\,\%$ \\ turnover between two seeds & $101\,\%$ & $12.6\,\%$ & $11.7\,\%$ & $7.2\,\%$ & $5.8\,\%$ \\ turnover to the optimum & $75\,\%$ & $8.9\,\%$ & $8.3\,\%$ & $5.0\,\%$ & $4.1\,\%$ \\ \bottomrule \end{tabular} \caption{Sampling noise of the mean--CVaR optimisation (30 stocks, eight seeds, CHF per week; ``turnover'' is $\sum_i|w_i-w_i'|$). ``mm'' = moment matching.} \label{tab:noise} \end{table} With plain Monte Carlo two seeds give essentially disjoint portfolios (turnover 101\,\%), the in-sample expected profit is 26\,\% too high and the true expected profit 10\,\% below the optimum. Fixing the first two moments of the set removes the illusion entirely; Sobol' points halve the remaining tail noise. The residual 4--6\,\% turnover is the flatness of the optimum (many portfolios within 0.02\,\% of the optimal objective), not sampling. Increasing $S$ does not rescue plain Monte Carlo: at $S=8\,192$ its turnover between seeds is still 76\,\% (the mean error decays like $1/\sqrt S$), while Sobol' + mm is at 4.0\,\%. \subsection{Quasi-Monte Carlo and moment matching} \label{sec:qmc} \paragraph{Sobol' points.} Instead of pseudo-random $z_s\sim N(0,I)$ the implementation uses scrambled Sobol' points \cite{Owen1998}: $S$ (rounded up to a power of two) points $u_s\in(0,1)^N$ that fill the unit cube far more evenly than random points, mapped through the inverse normal distribution function, $z_s=\Phi^{-1}(u_s)$. The scrambling is seeded, so the set is reproducible. The first coordinates of a Sobol' sequence are the best distributed, so the draws are coloured with the eigen-root $A=V\Lambda^{1/2}$ of $\Sigma$ with eigenvalues in descending order: the first coordinate drives the market factor, which dominates every long-only portfolio's loss \cite[ch.~5]{Glasserman2004}. \paragraph{Moment matching.} Whatever the base sample, its first two moments can be made \emph{exact} \cite{HoylandWallace2001}. With probabilities $p_s$ (equal unless importance sampling is used, Section~\ref{sec:is}), let $m=\sum_sp_sz_s$ and $C=\sum_sp_s(z_s-m)(z_s-m)^\top$; then \begin{equation} \tilde z_s \;=\; C^{-1/2}(z_s-m), \qquad \sum_sp_s\tilde z_s=0,\quad \sum_sp_s\tilde z_s\tilde z_s^\top=I, \label{eq:mm} \end{equation} and the scenarios $r_s=h m_0+\sqrt h\,A\tilde z_s$ have mean $hm_0$ and covariance $h\Sigma$ exactly, for any $S>N$. Only the \emph{shape} of the tail remains sampled. For the empirical innovation model the same transform is applied to the resampled sums of innovations (target covariance $hI$), so that both models reproduce the latest covariance exactly. \paragraph{The objective uses the model mean.} Once the mean is exact there is no reason to estimate it from the scenarios at all. The implementation puts the model's expected simple return $\mu_i$ (the Black--Litterman posterior, converted to the horizon) directly into the objective, $\bar g_i=P_i\mu_i$, and uses the scenarios only in the CVaR rows. This also keeps the objective exact under importance sampling, where a tail-focused sample estimates the mean \emph{worse} than a plain one. The convention behind $\mu_i$ matters: Black--Litterman returns are arithmetic expectations, so the log-drift of the scenarios is $\ln(1+h\mu_i)/h-\sigma_i^2/2$ (It\^o correction) and $\E[e^{r_i}-1]=h\mu_i$; using $\mu_i$ as the log-drift would add $\sigma_i^2/2$ to the expected simple return --- 5\,\% p.a.\ for a typical stock and 34\,\% for the most volatile one in the universe --- and tilt the optimiser towards volatile names. \subsection{Importance sampling for the tail} \label{sec:is} Only the scenarios in the tail of the loss matter for the CVaR: at $\beta=0.99$, 1\,\% of them. Importance sampling draws more scenarios where they matter and corrects with likelihood ratios \cite{GlassermanHeidelbergerShahabuddin2000}. Three things are specific to the optimisation setting. \emph{The tail depends on the decision.} The loss of portfolio $v=P\circ n$ is $L=-v^\top r=-v^\top hm_0-(\sqrt hA^\top v)^\top z$, so its tail is the half-space $\{z:\,d^\top z\ge c\}$ with the unit direction $d=-A^\top v/\|A^\top v\|$. Because all long-only, capped, sector-limited candidates load mainly on the market factor, one direction serves all of them. The default takes $v$ from the \emph{cap-weighted portfolio of the candidate stocks} (\code{importance\_sampling="market"}), which depends on the universe and the covariance but on none of the parameters a user varies; the alternative \code{"adaptive"} takes it from a first-pass solution of the problem itself, which is marginally sharper at $\beta=0.99$ but makes the scenario set depend on the risk limit and the constraints (Section~\ref{sec:onesetup}). \emph{Bounded likelihood ratios.} A share $\varepsilon$ (default one half) of the Sobol' points is shifted, $z_s\mapsto z_s+\theta$ with $\theta=z_\beta d$, which centres the loss of $v^{(1)}$ at its VaR; the rest stay nominal. The proposal is the defensive mixture $g=(1-\varepsilon)\varphi+\varepsilon\varphi(\cdot-\theta)$ \cite{Hesterberg1995}, whose likelihood ratios \begin{equation} \ell_s=\frac{\varphi(z_s)}{g(z_s)}=\frac{1}{(1-\varepsilon)+\varepsilon\,e^{\theta^\top z_s-\|\theta\|^2/2}} \;\le\;\frac{1}{1-\varepsilon} \label{eq:lr} \end{equation} are bounded, so a candidate whose tail points elsewhere is never wildly mis-estimated (the optimiser would exploit any such region). The scenario probabilities are the self-normalised ratios $p_s=\ell_s/\sum_{s'}\ell_{s'}$, and the moment matching \eqref{eq:mm} is applied with these weights. \emph{Weighted scenarios in the MILP.} Theorem~\ref{thm:ru} holds for any discrete distribution: with probabilities $p_s$ the CVaR row \eqref{eq:milp-n-cvar} becomes \begin{equation} \zeta+\frac{1}{1-\beta}\sum_sp_su_s\;\le\;\Lambda, \label{eq:cvar-weighted} \end{equation} and nothing else changes; the decoded VaR and CVaR are computed from the weighted sample in the same way (the tail is the set of largest losses whose probability adds up to $1-\beta$, the last one taken fractionally). For the empirical innovation model the tilt acts on the days: the resampling probabilities are $q_t\propto\exp(\vartheta\,d^\top\varepsilon_t)$ with $\vartheta$ chosen so that the tilted mean of $d^\top\varepsilon$ equals $z_\beta/\sqrt h$ per day (a one-dimensional root), the mixture draws all $h$ days of a scenario either from $q$ or uniformly, and the likelihood ratio of a scenario is $1/\big((1-\varepsilon)+\varepsilon\prod_kTq_{t_k}\big)$. \begin{table}[ht] \centering \footnotesize \setlength{\tabcolsep}{3.2pt} \begin{tabular}{@{}lrrrrrrrr@{}} \toprule & \multicolumn{2}{c}{$0.95$, $2\,048$} & \multicolumn{2}{c}{$0.95$, $8\,192$} & \multicolumn{2}{c}{$0.99$, $2\,048$} & \multicolumn{2}{c}{$0.99$, $8\,192$} \\ \cmidrule(lr){2-3}\cmidrule(lr){4-5}\cmidrule(lr){6-7}\cmidrule(lr){8-9} $\beta$, $S$ & S+mm & + IS & S+mm & + IS & S+mm & + IS & S+mm & + IS \\ \midrule scenarios beyond the VaR & 101 & 560 & 410 & 2\,249 & 19 & 512 & 80 & 2\,075 \\ true CVaR over the limit & 0.15\,\% & 0.21\,\% & 0.08\,\% & 0.01\,\% & 2.45\,\% & 0.50\,\% & 0.75\,\% & 0.03\,\% \\ sd of the true CVaR (\% of limit) & 0.74 & 0.30 & 0.36 & 0.15 & 1.67 & 0.40 & 0.79 & 0.20 \\ loss versus the optimum (CHF) & 0.4 & 1.1 & 0.3 & 0.0 & 10.7 & 2.7 & 3.6 & 0.1 \\ turnover between two seeds & 5.9\,\% & 3.6\,\% & 4.0\,\% & 1.6\,\% & 29\,\% & 11.5\,\% & 15.4\,\% & 5.7\,\% \\ \bottomrule \end{tabular} \caption{Importance sampling (``+ IS'': mixture $\varepsilon=0.5$, two-stage tilt) versus moment-matched Sobol' points (``S+mm'') at equal $S$; 30 stocks, eight seeds, exact Gaussian reference.} \label{tab:is} \end{table} At equal $S$ the tail noise falls by a factor 2--2.5 at 95\,\% and 3--4 at 99\,\%; equivalently, 2\,048 tilted scenarios are as good as 8\,192 plain ones. The gain is largest exactly where it is needed, at $\beta=0.99$, where a plain set of 2\,048 scenarios has 19 tail points. The tilted LP is slower per scenario (a quarter of its rows are active instead of $1-\beta$), but the market tilt is on by default (\code{importance\_sampling="auto"} $=$ \code{"market"}) because it costs nothing in accuracy at 95\,\% and halves the noise at 99\,\%. Table~\ref{tab:is-market} compares the two directions. Under the Gaussian model importance sampling is unnecessary --- the exact solution of Section~\ref{sec:gaussian} needs no scenarios at all. \begin{table}[ht] \centering \footnotesize \setlength{\tabcolsep}{4pt} \begin{tabular}{@{}lrrrrrr@{}} \toprule & \multicolumn{3}{c}{$\beta=0.95$} & \multicolumn{3}{c}{$\beta=0.99$} \\ \cmidrule(lr){2-4}\cmidrule(lr){5-7} $S=2\,048$ & no tilt & market & adaptive & no tilt & market & adaptive \\ \midrule true CVaR over the limit & 0.14\,\% & 0.37\,\% & 0.20\,\% & 2.43\,\% & 1.37\,\% & 0.50\,\% \\ sd of the true CVaR (\% of limit) & 0.74 & 0.72 & 0.30 & 1.69 & 0.61 & 0.42 \\ turnover between two seeds & 5.7\,\% & 3.2\,\% & 3.6\,\% & 28.5\,\% & 17.4\,\% & 11.6\,\% \\ turnover to the optimum & 4.1\,\% & 2.6\,\% & 2.4\,\% & 22.4\,\% & 13.0\,\% & 8.7\,\% \\ depends on the risk parameters & no & no & \textbf{yes} & no & no & \textbf{yes} \\ \bottomrule \end{tabular} \caption{Tilt direction: the cap-weighted candidate portfolio (``market'', two blocks at $z_{0.95}$ and $z_{0.99}$) versus a first-pass solution (``adaptive'', one block at $z_\beta$); 30 stocks, eight seeds.} \label{tab:is-market} \end{table} \emph{How far does one direction carry?} The tilt enriches the tail of the \emph{reference} portfolio; another portfolio benefits in proportion to how much its loss is correlated with the reference's. Table~\ref{tab:alignment} measures this for portfolios from market-like to deliberately orthogonal. Long-only portfolios under single-stock and sector caps are market-like in practice --- the optimiser's solutions come out at an alignment of about $0.90$ --- and the reference is the cap-weighted portfolio \emph{of the candidate stocks}, so a defensive universe gets a defensive reference. Where the alignment is nevertheless low the tilt does not help and costs a little precision (the nominal half of the mixture bounds the damage); the optimisation result reports the realised alignment (\code{tilt\_alignment}) and warns below $0.3$, where \code{importance\_sampling="adaptive"} or \code{"off"} is the better choice. \begin{table}[ht] \centering \footnotesize \setlength{\tabcolsep}{4pt} \begin{tabular}{@{}lrrrrrr@{}} \toprule & & \multicolumn{2}{c}{scenarios beyond the VaR} & \multicolumn{3}{c}{noise of the CVaR estimate} \\ \cmidrule(lr){3-4}\cmidrule(lr){5-7} portfolio & alignment & plain & tilted & plain & tilted & reduction \\ \midrule cap-weighted (the reference) & 1.00 & 41 & 789 & 0.99\,\% & 0.19\,\% & 5.3$\times$ \\ equal weight & 0.89 & 41 & 644 & 0.79\,\% & 0.30\,\% & 2.6$\times$ \\ mean--CVaR optimum & 0.90 & 41 & 630 & 1.39\,\% & 0.61\,\% & 2.3$\times$ \\ five highest-beta stocks & 0.68 & 41 & 368 & 0.22\,\% & 0.34\,\% & 0.6$\times$ \\ minimum variance & 0.69 & 41 & 350 & 1.34\,\% & 1.25\,\% & 1.1$\times$ \\ five lowest-beta stocks & $-0.09$ & 41 & 34 & 1.66\,\% & 1.96\,\% & 0.8$\times$ \\ three least market-like stocks & $-0.15$ & 41 & 31 & 1.93\,\% & 1.95\,\% & 1.0$\times$ \\ \bottomrule \end{tabular} \caption{Reach of one tilt direction: ``alignment'' is the correlation between the portfolio's loss and the reference's; $S=4\,096$, $\beta=0.99$, 12 seeds, 30 stocks. Noise is the standard deviation of the CVaR estimate across seeds, relative to its mean.} \label{tab:alignment} \end{table} \subsection{One scenario set for a family of runs} \label{sec:onesetup} A portfolio proposal is rarely produced once: a risk limit is varied, a constraint tightened, a view added, two confidence levels compared. If each of those runs drew its own scenarios, the differences between the results would mix the effect of the parameter with the sampling noise of a new set --- precisely the noise Sections~\ref{sec:qmc} and \ref{sec:is} remove. The implementation therefore makes the scenario set a function of the \emph{model} only. Its \emph{scenario key} (reported with every result) contains \begin{center} \footnotesize \begin{tabular}{@{}p{7.2cm}p{7.2cm}@{}} \toprule in the key (changing it draws new scenarios) & not in the key (the set is reused) \\ \midrule candidate stocks, data as-of date, frequency and horizon, covariance model, scenario model, number of scenarios, sampling scheme, moment matching, seed, tilt mode and mixture, innovation history and re-standardisation & CVaR limit (absolute or as a fraction), confidence level $\beta$, objective, risk aversion, target return, all constraints, transaction costs, cash, existing holdings, expected-return model and views, solver and time limit \\ \bottomrule \end{tabular} \end{center} Two properties make this possible. First, the tilt of Section~\ref{sec:is} uses the cap-weighted candidate portfolio and blocks at both usual confidence levels, so it does not look at the optimisation parameters. Second, the drift enters \emph{after} the random part: scenarios are generated with zero expected return and the expected-return model is added as a per-asset constant, $r_s\mapsto r_s+hm$, by \code{ScenarioSet.with\_drift}. Changing a view therefore shifts every scenario by a constant and leaves the shocks --- and hence the whole risk picture --- untouched; the result reports a fingerprint of the shocks (\code{shock\_hash}) next to the one of the scenarios (\code{hash}). Sets are cached per key inside the service, so a family of runs also pays the expensive innovation filter only once (measured: 14\,s for the first run, 0\,s afterwards). Only \code{importance\_sampling="adaptive"} breaks the property, since its tilt comes from a first-pass solve; the key then lists the parameters of that solve, so runs that differ in them are not silently compared on different scenarios. \subsection{How many scenarios, and how many assets?} \label{sec:howmany} \paragraph{The count has to grow with the candidate set.} The approximation searches an $N$-dimensional space, so a fixed $S$ buys less precision as $N$ grows. Measured against the exact optimum of Section~\ref{sec:gaussian} (GBM model, four seeds): \begin{center} \footnotesize \begin{tabular}{@{}rrrrrr@{}} \toprule $N$ & $S$ & $S/N$ & turnover between two seeds & turnover to the optimum & sd of $\E[\text{P\&L}]$ \\ \midrule 30 & 512 & 17 & 9.6\,\% & 7.2\,\% & 17.6 CHF \\ 30 & 2\,048 & 68 & 3.5\,\% & 4.2\,\% & 2.8 \\ 30 & 8\,192 & 273 & 2.6\,\% & 3.8\,\% & 1.1 \\ 100 & 512 & 5 & 26.5\,\% & 21.5\,\% & 15.5 \\ 100 & 2\,048 & 20 & 12.5\,\% & 10.0\,\% & 2.8 \\ 100 & 8\,192 & 82 & 6.2\,\% & 6.2\,\% & 2.0 \\ 300 & 512 & 1.7 & 22.7\,\% & 23.4\,\% & 18.5 \\ 300 & 2\,048 & 6.8 & 9.0\,\% & 6.8\,\% & 4.2 \\ \bottomrule \end{tabular} \end{center} Two readings of the same table matter. The \emph{reported figures} are stable everywhere: from 2\,000 scenarios on, the expected profit of the proposal varies by less than $0.1\,\%$ between seeds. The \emph{composition} is what moves, and at fixed $S$ it moves roughly in proportion to $N$ (3.5\,\% at $N=30$ against 12.5\,\% at $N=100$, both at $S=2\,048$) --- with hundreds of candidates and a 10\,\% cap there are many portfolios within a fraction of a percent of the optimum, and the sample decides between them. About 64 scenarios per candidate keeps the disagreement near 5\,\%; as a default that is too expensive ($N=300$ would mean 32\,768 scenarios and minutes per solve), so the implementation draws 16 per candidate within $[4\,096,\,16\,384]$ and says in the result when the ratio is low. Where the composition itself has to be reproducible --- a backtest, a teaching example --- either raise $S$ or use \code{scenario\_model="gaussian\_exact"}, which has no sampling error at all. \paragraph{Cost at the top end, and which solver.} A model with $N=300$ and $S=32\,768$ has 9.9\,M non-zeros: building it takes 3.3\,s, encoding it as msgpack 0.6\,s, and the 108\,MB body transfers to the cuOpt server at about 10\,MB/s. Transport is therefore not the constraint; proving optimality is. Table~\ref{tab:solvers} compares the two back ends with both running on the Spark, end to end through the service (scenario generation, build, solve, decode). \begin{table}[ht] \centering \footnotesize \begin{tabular}{@{}lrr@{}} \toprule case & cuOpt (GPU) & HiGHS (CPU, 4 threads) \\ \midrule $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) & 350.7\,s (optimal) \\ $N=300$, $S=8\,192$, row generation & 139\,s & 64\,s \\ \bottomrule \end{tabular} \caption{The same models on the same machine, weekly horizon, $\beta=0.95$, 600\,s solver time limit.} \label{tab:solvers} \end{table} HiGHS is at least as fast in every case and twice as fast on the two largest. This is not a statement about GPU optimisation in general. cuOpt's published speedups concern \emph{linear} programmes solved with its PDLP first-order method; for mixed-integer problems its documentation states that the GPU runs primal heuristics while \emph{branch and bound runs on the CPU}, that the MIP solver is in beta, and that proving feasible solutions optimal ``remains under active development''. The models here have an easy LP relaxation and a small integrality gap, so nearly all the work is the optimality proof --- the part that is not accelerated --- whereas HiGHS warm-starts its dual simplex across branch-and-bound nodes and, with row generation, solves a sequence of small problems in-process while cuOpt pays a fixed cost per request. HiGHS is a demanding reference: in the Mittelmann benchmark of 7 July 2026 it solves 179 of 240 instances at a scaled shifted geometric mean of 5.44, against 1.72 for the commercial COPT and 9.93 for SCIP. Three caveats: cuOpt's time limit covers only its internal solve, it ignores the MIP-gap settings sent over REST and may therefore prove optimality to a tighter tolerance, and this is one problem family on one machine rather than a benchmark. \subsection{Tail-row generation for large scenario sets} \label{sec:rowgen} The MILP has one row and one variable $u_s$ per scenario; with $N+2$ non-zeros per row the matrix grows linearly in $S$, but the simplex work does not: HiGHS' dual simplex on the LP relaxation scales roughly like $S^{1.9}$ and the branch-and-bound like $S^{1.65}$, cuOpt like $S^{1.2}$--$S^{1.6}$ (Table~\ref{tab:scaling}). Going from 2\,000 to 100\,000 scenarios costs about 45$\times$ on cuOpt and more than 200$\times$ on HiGHS. \begin{table}[ht] \centering \small \begin{tabular}{@{}rrrrrrr@{}} \toprule $S$ & non-zeros & HiGHS LP & HiGHS MILP & cuOpt MILP & payload & row generation (HiGHS) \\ \midrule \multicolumn{7}{@{}l}{$N=30$} \\ 2\,048 & 68\,k & 0.03\,s & 1.5\,s & 5.1\,s & 1.7\,MB & -- \\ 8\,192 & 270\,k & 0.4\,s & 14.5\,s & 10.8\,s & 6.6\,MB & 5.4\,s (2 iterations) \\ 32\,768 & 1.1\,M & 6.0\,s & 166\,s & 66.6\,s & 27\,MB & 25\,s (3 iterations) \\ 100\,000 & 3.3\,M & 56.5\,s & $>300$\,s & 220\,s & 82\,MB & 158\,s (2 iterations) \\ \multicolumn{7}{@{}l}{$N=100$} \\ 2\,048 & 211\,k & 0.1\,s & 9.5\,s & 7.0\,s & 5\,MB & \\ 8\,192 & 844\,k & 1.2\,s & 50\,s & 11.5\,s & 21\,MB & \\ 32\,768 & 3.4\,M & 16.5\,s & 232\,s & 101\,s & 82\,MB & \\ \bottomrule \end{tabular} \caption{Solve times versus the number of scenarios (MIP gap $10^{-4}$, $\beta=0.95$). The HiGHS column was measured on a development laptop whose SciPy bundles HiGHS 1.2.0 and the cuOpt column on the DGX Spark, so the two columns show the growth with $S$ for each solver but are not a like-for-like comparison; Table~\ref{tab:solvers} compares them on the same machine.} \label{tab:scaling} \end{table} At the optimum only the scenarios in the tail have $u_s>0$; every other scenario row is slack. \emph{Row generation} \cite{KunziBayMayer2006} exploits this: \begin{enumerate}[nosep] \item Start with a working set $\mathcal W$: the $2(1-\beta)S$ worst scenarios of a reference portfolio (the current holdings, or an equal-weight allocation). \item Solve the MILP with the scenario rows of $\mathcal W$ only, keeping the full-set probabilities $p_s$ in \eqref{eq:cvar-weighted} (omitted scenarios contribute $u_s=0$) and the full-set objective. \item Compute the losses $L_s$ of all $S$ scenarios for the solution $(n^\ast,\zeta^\ast)$; add every $s\notin\mathcal W$ with $L_s>\zeta^\ast$ to $\mathcal W$ and repeat; stop when there is none. \end{enumerate} The restricted problem drops constraints, so it is a relaxation and its objective bounds the full one from above; when no omitted scenario is violated the restricted solution is feasible for the full problem with $u_s=0$ on the omitted rows, hence optimal (within the MIP gap). The working set stays at about 11\,\% of $S$ at $\beta=0.95$ (2--3\,\% at $\beta=0.99$), two or three iterations suffice, and the speed-up over the full MILP is 6--7$\times$ from 32\,768 scenarios on; only the working-set rows are sent to the solver. The default (\code{row\_generation="auto"}) uses it from 8\,192 scenarios on. \subsection{The Gaussian special case: mean--variance equivalence and an exact solution} \label{sec:gaussian} \begin{proposition} \label{prop:elliptic} Let the horizon returns be elliptically distributed with mean $m$ and dispersion $C$, e.g.\ $R\sim N(m,C)$. Then for every portfolio $\CVaR_\beta(L(w))=-m^\top w+k_\beta\,\sigma(w)$ with $\sigma(w)=\sqrt{w^\top Cw}$ and a constant $k_\beta$ ($k_\beta=\varphi(z_\beta)/(1-\beta)$ in the normal case, $2.063$ at 95\,\% and $2.665$ at 99\,\%). Consequently the solution of ``maximise $m^\top w$ subject to $\CVaR_\beta\le\Lambda$'' lies on the mean--variance efficient frontier, and it is also the solution of ``maximise $m^\top w$ subject to $\sigma(w)\le\bar\sigma$'' with the implicit bound $\bar\sigma=(\Lambda+m^\top w^\ast)/k_\beta$. \end{proposition} \begin{proof} For elliptical returns the loss of every portfolio is a location--scale transform of one and the same symmetric distribution, so every positively homogeneous, translation-invariant risk measure is an affine function of $\sigma(w)$ \cite{EmbrechtsMcNeilStraumann2002}. Let $w^\ast$ be mean--CVaR optimal with mean $m^\ast$. Any feasible $w$ with $m^\top w\ge m^\ast$ and $\sigma(w)\le\bar\sigma$ has $\CVaR_\beta\le-m^\ast+k_\beta\bar\sigma=\Lambda$, hence is mean--CVaR feasible and cannot beat $m^\ast$: the two problems have the same optimal value, and the same solutions up to ties. The same argument applies over the integer feasible set. \end{proof} The variance bound is only known \emph{ex post} (it depends on $m^\ast$), so the mean--CVaR problem with a Gaussian model is solved directly. In whole shares it is \begin{equation} \max_n\; m^\top v-\TC \quad\text{s.t.}\quad g(n):=\TC-m^\top v+k_\beta\sqrt{v^\top Cv}\;\le\;\Lambda,\qquad v=P\circ n, \label{eq:misocp} \end{equation} plus the linear constraints of Section~\ref{sec:shares}: a convex mixed-integer second-order-cone programme (a mixed-integer QP in the utility form), which neither cuOpt nor HiGHS solves directly. The function $g$ is convex, so it can be outer-approximated by tangent cuts on the same MILP \cite{Kelley1960,DuranGrossmann1986}. At an exposure $v_k$ with $\sigma_k=\sqrt{v_k^\top Cv_k}>0$ the gradient of $\sigma$ is $Cv_k/\sigma_k$, and because $(Cv_k)^\top v_k/\sigma_k=\sigma_k$ the tangent simplifies to \begin{equation} g(n)\;\ge\;\TC+\Big(k_\beta\frac{Cv_k}{\sigma_k}-m\Big)^{\!\top}v . \label{eq:cut} \end{equation} With an epigraph variable $t$ (\code{cvar} in the code) the MILP carries the rows $\TC+(k_\beta Cv_k/\sigma_k-m)^\top(P\circ n)-t\le0$ for the cuts collected so far and $t\le\Lambda$ (or $t$ in the objective for \code{min\_cvar} and \code{utility}). The loop starts with a cut at a reference portfolio, solves the MILP, compares the true $g(n^\ast)$ with $t^\ast$, and adds a cut at $v^\ast$ until they agree (relative tolerance $10^{-6}$). Each cut removes the current infeasible point and the integer set is finite, so the loop terminates at the exact optimum within the MIP gap. On the universe it takes 20--70 cuts (each a MILP with $N$ integers and a few dozen rows, milliseconds on HiGHS) and reproduces Gurobi's direct MISOCP solution to the cent. This is \code{scenario\_model="gaussian\_exact"}: no scenarios, no sampling error, fully deterministic. \begin{remark} The GBM scenario model is lognormal in simple returns ($R=e^r-1$ with $r$ normal), the exact mode is normal in simple returns. At weekly volatilities of 4--5\,\% the lognormal left tail is 2--3\,\% thinner, so the scenario optimum takes a few percent more exposure; the sharp test of the sampling machinery therefore uses Gaussian scenarios linearly, where moment-matched Sobol' points reproduce the exact optimum within 1\,\% and plain Monte Carlo does not (\code{tests/test\_rowgen\_gaussian.py}). For the empirical innovation model no closed form exists --- that is the point of using it --- and the equivalence serves only as a calibration of the sampling mechanism. \end{remark} \subsection{Rank of the covariance estimate} \label{sec:rank} An estimator with $T_{\mathrm{eff}}$ effective observations cannot support more than $T_{\mathrm{eff}}$ assets: a 200-day window with $N=300$ gives a sample covariance of rank 200, and a portfolio in its null space has \emph{zero} estimated variance --- an optimiser maximising return under a risk cap will find it. The shrinkage of the measurement note repairs exactly this, and its data-driven intensity rises with $N/T$ on its own: \begin{center} \footnotesize \begin{tabular}{@{}lrrrrr@{}} \toprule estimator, $N$ & $T_{\mathrm{eff}}$ & intensity & $\lambda_{\min}/\overline{\sigma^2}$ & condition & zero eigenvalues \\ \midrule rolling 200, $N=300$ & 200 & 0 & $-3.4\cdot10^{-15}$ & $4\cdot10^{298}$ & 101 \\ rolling 200, $N=300$ & 200 & 0.080 & 0.023 & $3.0\cdot10^{3}$ & 0 \\ rolling 200, $N=745$ & 200 & 0 & $-2.0\cdot10^{-14}$ & $6\cdot10^{298}$ & 546 \\ rolling 200, $N=745$ & 200 & 0.128 & 0.025 & $4.3\cdot10^{3}$ & 0 \\ EWM $\lambda=0.99$, $N=745$ & 300 & 0 & $2.2\cdot10^{-5}$ & $5.5\cdot10^{6}$ & 44 (numerically) \\ EWM $\lambda=0.99$, $N=745$ & 300 & 0.114 & 0.023 & $4.8\cdot10^{3}$ & 0 \\ \bottomrule \end{tabular} \end{center} With \code{shrinkage="auto"} --- the default everywhere --- the smallest eigenvalue never falls below 2\,\% of the average variance, whatever the number of candidates. Without it the degeneracy is real and silent: a null-space portfolio comes out at variance $-1.6\cdot10^{-20}$ and the Sobol' draws span only 222 of 300 dimensions. \code{CovarianceModel.conditioning} reports the figures and the optimisation warns when the estimate is singular or ill-conditioned for the size of the candidate set. Two things soften the problem in this particular model, and neither should be relied on. An EWM estimator over the full history is formally full rank (every observation keeps a positive weight), so its degeneracy is numerical rather than exact; and no-short-selling with a 10\,\% cap makes null directions hard to reach, since exploiting one generally requires a short leg. Measured across $N=30\ldots745$, the optimised portfolio's out-of-sample-to-in-sample volatility ratio exceeds the cap-weighted portfolio's by a factor of only 0.89--1.18, with or without shrinkage --- the optimiser is not systematically mining the estimate. The exactly singular case does bite elsewhere, though: with a 200-day window and $N\ge300$ the outer approximation of Section~\ref{sec:gaussian} stops converging and returns a time-limited solution. \subsection{Reproducibility} \label{sec:reproducibility} Two runs of the same request use the same scenario set, and so do two runs that differ only in the parameters listed on the right of Section~\ref{sec:onesetup}: the base sample depends only on the seed (scrambled Sobol' points and the PCG64 generator are deterministic across platforms), the transformations are deterministic, and the result reports a fingerprint of the final matrix and its probabilities (\code{scenario\_hash}). The fingerprint is exact on one platform; between BLAS implementations (laptop versus the DGX Spark) the GBM sets agree to $10^{-11}$ and keep the same hash, while the empirical model's innovation filter --- thousands of eigendecompositions --- can differ at the $10^{-10}$ level, which changes the hash but none of the reported figures. Whether the \emph{solution} is the same depends on the solver. HiGHS' branch-and-bound is deterministic: identical inputs give identical shares (\code{reproducibility="exact"}) as long as the HiGHS build is the same --- across SciPy versions the search may end on a different optimal vertex among ties, at the same objective. cuOpt combines GPU heuristics with a time-limited search and, in version 26.08, ignores the MIP-gap settings sent over its REST interface; it proves optimality within its default relative gap of $10^{-4}$, so the objective, expected profit and CVaR of two runs agree to that gap (about 0.2 CHF per week on a CHF 1 million portfolio) while near-tied share counts can differ by a few shares (\code{reproducibility="within\_gap"}). A result that hit the time limit (\code{FeasibleFound}) is not reproducible at all. The Gaussian exact mode runs on HiGHS by default and is exact in this sense. % ============================================================================================ \section{The complete model as implemented} % ============================================================================================ \code{MeanCVaRProblem.build()} assembles the following MILP from a \code{ProblemSpec}. The variable vector is laid out as $x=(n,\;[b,\;s],\;\alpha,\;u,\;[y],\;[t])$; blocks in brackets exist only in rebalancing mode ($b,s$), with cardinality/minimum-position rules ($y$) or with a fixed cost per trade ($t$). The CVaR level $\zeta$ is called $\alpha$ in the code. \begin{center} \small \begin{tabular}{@{}llll@{}} \toprule block & name in the code & type, bounds & meaning \\ \midrule $n_i$ & \code{n[symbol]} & integer, $[0,\bar n_i]$ & final number of shares \\ $b_i$ & \code{buy[symbol]} & integer, $[0,\bar n_i]$ & shares bought (rebalancing) \\ $s_i$ & \code{sell[symbol]} & integer, $[0,n^0_i]$ & shares sold (rebalancing) \\ $\alpha$ & \code{alpha} & continuous, free & the $\zeta$ of Theorem~\ref{thm:ru} \\ $u_s$ & \code{u[s]} & continuous, $[0,\infty)$ & excess loss of scenario $s$ over $\alpha$ \\ $t_{\mathrm{cvar}}$ & \code{cvar} & continuous, free & Gaussian CVaR epigraph (Section~\ref{sec:gaussian}) \\ $y_i$ & \code{y[symbol]} & binary & position open (cardinality / minimum position) \\ $t_i$ & \code{t[symbol]} & binary & stock traded (fixed cost) \\ \bottomrule \end{tabular} \end{center} With $\bar n_i=\lfloor \min(\bar w,f_{\max})\,W/P_i\rfloor$, $g_{s,i}=P_iR_{s,i}$, $\bar g_i=P_i\mu_i$ (the model's expected simple return over the horizon; the probability-weighted scenario mean if no model mean is supplied), $\kappa_s=\frac{p_s}{1-\beta}$ with $p_s=1/S$ unless the scenarios carry importance-sampling probabilities, and $\TC=\sum_i(c^bP_ib_i+c^sP_is_i)+F\sum_it_i$ (in mode \code{new}: $c^b\sum_iP_in_i+F\sum_it_i$), the rows are: \begin{center} \footnotesize \begin{tabular}{@{}lp{6.1cm}p{2.9cm}l@{}} \toprule row name & constraint & count & origin \\ \midrule \code{cvar\_scen\_s} & $\sum_i g_{s,i}n_i - \TC + \alpha + u_s \ge 0$ & $S$ & \eqref{eq:lin}, \eqref{eq:u-reb} \\ \code{cvar\_limit} & $\alpha + \sum_s \kappa_s u_s \le \Lambda$ \newline (Gaussian: $t_{\mathrm{cvar}}\le\Lambda$) & 1, if a limit is given & \eqref{eq:cvar-constraint}, \eqref{eq:cvar-weighted} \\ \code{cut\_j} & $\TC+(k_\beta Cv_j/\sigma_j-m)^{\top}(P\circ n)$ \newline $\quad-\,t_{\mathrm{cvar}}\le0$ & one per cut, \newline Gaussian mode & \eqref{eq:cut} \\ \code{budget} & $\sum_iP_ib_i-\sum_iP_is_i+\TC\le c^0$ \newline (new: $\sum_iP_in_i+\TC\le c^0$) & 1 & \eqref{eq:cash} \\ \code{link\_i} & $n_i-b_i+s_i=n^0_i$ & $N$, rebalancing & \eqref{eq:link} \\ \code{invested} & $f_{\min}W\le\sum_iP_in_i\le f_{\max}W$ & $\le1$, if bounds given & invested fraction \\ \code{sector\_k} & $\sum_{i\in G_k}P_in_i\le\bar w_kW$ & $K$, binding sectors & \eqref{eq:milp-n-sector} \\ \code{max\_cost} & $\TC\le\TC_{\max}$ & $\le1$ & \eqref{eq:turnover} \\ \code{turnover} & $\sum_iP_i(b_i+s_i)\le\tau W$ & $\le1$, rebalancing & \eqref{eq:turnover} \\ \code{open\_i} & $n_i-\max(\bar n_i,1)\,y_i\le0$ & $N$, optional & \eqref{eq:card} \\ \code{minpos\_i} & $P_in_i-v_{\min}y_i\ge0$ & $N$, optional & \eqref{eq:card} \\ \code{cardinality} & $\sum_iy_i\le K_{\max}$ & $\le1$, optional & \eqref{eq:card} \\ \code{trade\_i} & $b_i+s_i-M_it_i\le0$ \newline (new: $n_i-M_it_i\le0$) & $N$, fixed cost only & \eqref{eq:tc} \\ \code{target\_pnl} & $\sum_i\bar g_in_i-\TC\ge\rho$ & 1, \code{min\_cvar} only & \eqref{eq:three}(b) \\ \bottomrule \end{tabular} \end{center} The objective is one of \begin{align} \code{max\_return}:&\quad \max\; \sum_i\bar g_i n_i - \TC \quad\text{subject to \code{cvar\_limit}},\\ \code{min\_cvar}:&\quad \min\; \alpha+\sum_s\kappa_s u_s \quad\text{subject to \code{target\_pnl}},\\ \code{utility}:&\quad \max\; \sum_i\bar g_i n_i - \TC - \lambda\Big(\alpha+\sum_s\kappa_s u_s\Big), \end{align} the three forms of \eqref{eq:three} (in Gaussian mode $\alpha+\sum_s\kappa_su_s$ is replaced by $t_{\mathrm{cvar}}$). All coefficients are in CHF; the scenario rows carry $N(+2N)+2$ non-zeros each, so the matrix has about $S(N+2)$ non-zeros (420\,k for $N=100$, $S=4\,096$; 3\,M for the full universe of 745 stocks with $S=4\,096$). Duplicate index pairs that arise when a variable appears in several terms of a row (e.g.\ $n_i$ in the gain and in the cost of mode \code{new}) are summed when the sparse matrix is assembled. \paragraph{The two solvers.} The same \code{MILPModel} is accepted by both back ends. \emph{cuOpt} is NVIDIA's GPU optimisation library, reached over HTTP as a self-hosted server (version 26.08 here) running as a container on the DGX Spark; it is the only component that uses the GPU. \emph{HiGHS} is the open-source LP/MILP solver of ERGO-Code at the University of Edinburgh (MIT licence, \url{https://highs.dev}) \cite{HuangfuHall2018}; it is not installed or started separately but ships compiled inside SciPy, so \code{scipy.optimize.milp} runs it \emph{in the calling process on the CPU} --- in production inside the same container as the service (SciPy 1.18.1, HiGHS 1.12.0, four threads). It serves as the fallback when cuOpt is unreachable, as the solver of the small cut MILPs of the Gaussian mode, and as the oracle of the test suite. \paragraph{Solving and decoding.} The model is sent to the cuOpt server as an LP/MILP request (CSR matrix, bounds, variable types, a mandatory time limit; msgpack-encoded) and polled until a solution arrives; if the server is unreachable the same \code{MILPModel} is handed to HiGHS. From 8\,192 scenarios on the service solves by row generation (Section~\ref{sec:rowgen}), i.e.\ a sequence of such MILPs on the working set; in Gaussian mode by the outer-approximation loop (Section~\ref{sec:gaussian}). A status \code{Optimal} means proven optimality within the solver's relative gap; \code{FeasibleFound} means the time limit stopped the search at an incumbent whose \code{mip\_gap} is reported. The result states the scenario fingerprint and the level of reproducibility (Section~\ref{sec:reproducibility}). From the solution vector, \code{MeanCVaRProblem.decode()} rounds the integer blocks (they are integral up to solver tolerance), rebuilds trades, costs, cash, position values and weights, and \emph{recomputes} the expected profit (from the model means), VaR and CVaR of the decoded portfolio from the weighted scenario set (Remark~\ref{rem:discrete}) or, in Gaussian mode, from the closed form; these recomputed figures are what the service reports, together with the provenance of every model input (covariance model, scenario model and seed, expected returns and views, constraints, costs, solver statistics). \paragraph{What the model does not do.} It assumes the scenario set represents the horizon distribution and that all trades execute at today's prices; it holds the portfolio unchanged over the horizon (no intra-horizon rebalancing); cash earns nothing over the horizon (the CHF risk-free rate is zero); it is single-period. These are deliberate simplifications of a teaching platform, not limitations of the formulation --- multi-period extensions replace the scenario rows by a scenario tree. % ============================================================================================ \begin{thebibliography}{99} \bibitem{RockafellarUryasev2000} R.\,T.~Rockafellar and S.~Uryasev, \emph{Optimization of conditional value-at-risk}, Journal of Risk 2(3), 21--41, 2000. \bibitem{RockafellarUryasev2002} R.\,T.~Rockafellar and S.~Uryasev, \emph{Conditional value-at-risk for general loss distributions}, Journal of Banking \& Finance 26(7), 1443--1471, 2002. \bibitem{KrokhmalPalmquistUryasev2002} P.~Krokhmal, J.~Palmquist and S.~Uryasev, \emph{Portfolio optimization with conditional value-at-risk objective and constraints}, Journal of Risk 4(2), 43--68, 2002. \bibitem{AcerbiTasche2002} C.~Acerbi and D.~Tasche, \emph{On the coherence of expected shortfall}, Journal of Banking \& Finance 26(7), 1487--1503, 2002. \bibitem{LoboFazelBoyd2007} M.\,S.~Lobo, M.~Fazel and S.~Boyd, \emph{Portfolio optimization with linear and fixed transaction costs}, Annals of Operations Research 152, 341--365, 2007. \bibitem{KellererMansiniSperanza2000} H.~Kellerer, R.~Mansini and M.\,G.~Speranza, \emph{Selecting portfolios with fixed costs and minimum transaction lots}, Annals of Operations Research 99, 287--304, 2000. \bibitem{HuangfuHall2018} Q.~Huangfu and J.\,A.\,J.~Hall, \emph{Parallelizing the dual revised simplex method}, Mathematical Programming Computation 10, 119--142, 2018. \bibitem{ShapiroDentchevaRuszczynski2009} A.~Shapiro, D.~Dentcheva and A.~Ruszczy\'nski, \emph{Lectures on Stochastic Programming: Modeling and Theory}, SIAM, Philadelphia, 2009. \bibitem{Michaud1989} R.\,O.~Michaud, \emph{The Markowitz optimization enigma: is `optimized' optimal?}, Financial Analysts Journal 45(1), 31--42, 1989. \bibitem{Owen1998} A.\,B.~Owen, \emph{Scrambling Sobol' and Niederreiter--Xing points}, Journal of Complexity 14(4), 466--489, 1998. \bibitem{Glasserman2004} P.~Glasserman, \emph{Monte Carlo Methods in Financial Engineering}, Springer, New York, 2004. \bibitem{HoylandWallace2001} K.~H\o yland and S.\,W.~Wallace, \emph{Generating scenario trees for multistage decision problems}, Management Science 47(2), 295--307, 2001. \bibitem{GlassermanHeidelbergerShahabuddin2000} P.~Glasserman, P.~Heidelberger and P.~Shahabuddin, \emph{Variance reduction techniques for estimating value-at-risk}, Management Science 46(10), 1349--1364, 2000. \bibitem{Hesterberg1995} T.~Hesterberg, \emph{Weighted average importance sampling and defensive mixture distributions}, Technometrics 37(2), 185--194, 1995. \bibitem{KunziBayMayer2006} A.~K\"unzi-Bay and J.~Mayer, \emph{Computational aspects of minimizing conditional value-at-risk}, Computational Management Science 3(1), 3--27, 2006. \bibitem{EmbrechtsMcNeilStraumann2002} P.~Embrechts, A.\,J.~McNeil and D.~Straumann, \emph{Correlation and dependence in risk management: properties and pitfalls}, in: M.~Dempster (ed.), Risk Management: Value at Risk and Beyond, Cambridge University Press, 176--223, 2002. \bibitem{Kelley1960} J.\,E.~Kelley, \emph{The cutting-plane method for solving convex programs}, Journal of the SIAM 8(4), 703--712, 1960. \bibitem{DuranGrossmann1986} M.\,A.~Duran and I.\,E.~Grossmann, \emph{An outer-approximation algorithm for a class of mixed-integer nonlinear programs}, Mathematical Programming 36, 307--339, 1986. \bibitem{MansiniOgryczakSperanza2014} R.~Mansini, W.~Ogryczak and M.\,G.~Speranza, \emph{Twenty years of linear programming based portfolio optimization}, European Journal of Operational Research 234(2), 518--535, 2014. \end{thebibliography} \end{document}