Portfolio_Risk_Management/docs/optimization_model.tex

% !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}