"""Base Convex Optimization estimator."""
# Copyright (c) 2023-2026
# Author: Hugo Delatte <hugo.delatte@skfoliolabs.com>
# SPDX-License-Identifier: BSD-3-Clause
# The optimization features are derived
# from Riskfolio-Lib, Copyright (c) 2020-2023, Dany Cajas, Licensed under BSD 3 clause.
from __future__ import annotations
import functools
import warnings
from abc import ABC, abstractmethod
from collections.abc import Callable
from enum import auto
from typing import Any, TypeVar
import cvxpy as cp
import cvxpy.constraints.constraint as cpc
import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as sla
import sklearn.utils.metadata_routing as skm
from cvxpy.reductions.solvers.defines import MI_SOLVERS
import skfolio.typing as skt
from skfolio._constants import (
_MANAGEMENT_FEES,
_TRANSACTION_COSTS,
)
from skfolio.exceptions import ConvexOptimizationError, OptimizationError
from skfolio.measures import RiskMeasure, owa_gmd_weights
from skfolio.optimization._base import BaseOptimization
from skfolio.prior import BasePrior, ReturnDistribution
from skfolio.typing import AnyArray, ArrayLike, BoolArray, FloatArray, StrArray
from skfolio.uncertainty_set import (
BaseCovarianceUncertaintySet,
BaseMuUncertaintySet,
CompactCovarianceUncertaintySet,
UncertaintySet,
)
from skfolio.utils.equations import equations_to_matrix, group_cardinalities_to_matrix
from skfolio.utils.tools import (
AutoEnum,
_get_liquidation_turnover_and_cost,
_validate_positive_real,
cache_method,
input_to_array,
)
INSTALLED_SOLVERS = cp.installed_solvers()
# Conic residuals can represent an equality-constrained zero as a tiny positive
# value. Dividing homogeneous weights by such a value is unsafe.
_MIN_HOMOGENIZATION_FACTOR = 1e-12
_ResultT = TypeVar("_ResultT")
[docs]
class ObjectiveFunction(AutoEnum):
r"""Enumeration of objective functions.
Attributes
----------
MINIMIZE_RISK : str
Minimize the risk measure.
MAXIMIZE_RETURN : str
Maximize the expected return.
MAXIMIZE_UTILITY : str
Maximize the utility :math:`w^T\mu - \lambda \times risk(w)`.
MAXIMIZE_RATIO : str
Maximize the ratio :math:`\frac{w^T\mu - R_{f}}{risk(w)}`.
"""
MINIMIZE_RISK = auto()
MAXIMIZE_RETURN = auto()
MAXIMIZE_UTILITY = auto()
MAXIMIZE_RATIO = auto()
[docs]
class ConvexOptimization(BaseOptimization, ABC):
r"""Base class for all convex optimization estimators in skfolio.
All risk measures that have a convex formulation are defined in class methods with
naming convention: `_{risk_measure}_risk`. That naming convention is used for
dynamic lookup.
CVX expressions that are shared among multiple risk measures are cached in a
dictionary named `_cvx_cache`.
This is to avoid cvx expression duplication and improve performance and convergence.
Parameters
----------
risk_measure : RiskMeasure, default=RiskMeasure.VARIANCE
:class:`~skfolio.meta.RiskMeasure` of the optimization.
Can be any of:
* VARIANCE
* SEMI_VARIANCE
* STANDARD_DEVIATION
* SEMI_DEVIATION
* MEAN_ABSOLUTE_DEVIATION
* FIRST_LOWER_PARTIAL_MOMENT
* CVAR
* EVAR
* WORST_REALIZATION
* CDAR
* MAX_DRAWDOWN
* AVERAGE_DRAWDOWN
* EDAR
* ULCER_INDEX
* GINI_MEAN_DIFFERENCE_RATIO
The default is `RiskMeasure.VARIANCE`.
prior_estimator : BasePrior, optional
:ref:`Prior estimator <prior>`.
The prior estimator is used to estimate the :class:`~skfolio.prior.ReturnDistribution`
containing estimates of expected asset returns, covariance matrix,
returns and Cholesky decomposition of the covariance.
The default (`None`) is to use :class:`~skfolio.prior.EmpiricalPrior`.
min_weights : float | dict[str, float] | array-like of shape (n_assets, ) | None, default=0.0
Minimum assets weights (weights lower bounds).
If a float is provided, it is applied to each asset.
`None` is equivalent to `-np.Inf` (no lower bound).
If a dictionary is provided, its (key/value) pair must be the
(asset name/asset minimum weight) and the input `X` of the `fit` method must
be a DataFrame with the assets names in columns.
When using a dictionary, assets values that are not provided are assigned
a minimum weight of `0.0`.
The default value is `0.0` (no short selling).
Example:
* `min_weights = 0` --> long only portfolio (no short selling).
* `min_weights = None` --> no lower bound (same as `-np.Inf`).
* `min_weights = -2` --> each weight must be above -200%.
* `min_weights = {"SX5E": 0, "SPX": -2}`
* `min_weights = [0, -2]`
max_weights : float | dict[str, float] | array-like of shape (n_assets, ) | None, default=1.0
Maximum assets weights (weights upper bounds).
If a float is provided, it is applied to each asset.
`None` is equivalent to `+np.Inf` (no upper bound).
If a dictionary is provided, its (key/value) pair must be the
(asset name/asset maximum weight) and the input `X` of the `fit` method must
be a DataFrame with the assets names in columns.
When using a dictionary, assets values that are not provided are assigned
a minimum weight of `1.0`.
The default value is `1.0` (each asset is below 100%).
Example:
* `max_weights = 0` --> no long position (short only portfolio).
* `max_weights = None` --> no upper bound.
* `max_weights = 2` --> each weight must be below 200%.
* `max_weights = {"SX5E": 1, "SPX": 2}`
* `max_weights = [1, 2]`
budget : float | None, default=1.0
Investment budget. It is the sum of long positions and short positions (sum of
all weights). `None` means no budget constraints.
The default value is `1.0` (fully invested portfolio).
For example:
* `budget = 1` --> fully invested portfolio.
* `budget = 0` --> market neutral portfolio.
* `budget = None` --> no constraints on the sum of weights.
min_budget : float, optional
Minimum budget. It is the lower bound of the sum of long and short positions
(sum of all weights). If provided, you must set `budget=None`.
The default (`None`) means no minimum budget constraint.
max_budget : float, optional
Maximum budget. It is the upper bound of the sum of long and short positions
(sum of all weights). If provided, you must set `budget=None`.
The default (`None`) means no maximum budget constraint.
max_short : float, optional
Maximum short position. The short position is defined as the sum of negative
weights (in absolute term).
The default (`None`) means no maximum short position.
max_long : float, optional
Maximum long position. The long position is defined as the sum of positive
weights.
The default (`None`) means no maximum long position.
cardinality : int, optional
Specifies the cardinality constraint to limit the number of invested assets
(non-zero weights). This feature requires a mixed-integer solver. For an
open-source option, we recommend using SCIP by setting `solver="SCIP"`.
To install it, use: `pip install cvxpy[SCIP]`. For commercial solvers,
supported options include MOSEK, GUROBI, or CPLEX.
group_cardinalities : dict[str, int], optional
A dictionary specifying cardinality constraints for specific groups of assets.
The keys represent group names (strings), and the values specify the maximum
number of assets allowed in each group. You must provide the groups using the
`groups` parameter. This requires a mixed-integer solver (see `cardinality`
for more details).
threshold_long : float | dict[str, float] | array-like of shape (n_assets, ), optional
Specifies the minimum weight threshold for assets in the portfolio to be
considered as a long position. Assets with weights below this threshold
will not be included as part of the portfolio's long positions. This
constraint can help eliminate insignificant allocations.
This requires a mixed-integer solver (see `cardinality` for more details).
It follows the same format as `min_weights` and `max_weights`.
threshold_short : float | dict[str, float] | array-like of shape (n_assets, ), optional
Specifies the maximum weight threshold for assets in the portfolio to be
considered as a short position. Assets with weights above this threshold
will not be included as part of the portfolio's short positions. This
constraint can help control the magnitude of short positions.
This requires a mixed-integer solver (see `cardinality` for more details).
It follows the same format as `min_weights` and `max_weights`.
transaction_costs : float | dict[str, float] | array-like of shape (n_assets, ), default=0.0
Transaction costs of the assets. It is used to add linear transaction costs to
the optimization problem:
.. math:: total\_cost = \sum_{i=1}^{N} c_{i} \times |w_{i} - w\_prev_{i}|
with :math:`c_{i}` the transaction cost of asset i, :math:`w_{i}` its weight
and :math:`w\_prev_{i}` its previous weight (defined in `previous_weights`).
The float :math:`total\_cost` is impacting the portfolio expected return in the optimization:
.. math:: expected\_return = \mu^{T} \cdot w - total\_cost
with :math:`\mu` the vector of assets' expected returns and :math:`w` the
vector of assets weights.
For positions in `previous_weights` whose assets are no longer in the
investment universe, transaction costs are calculated assuming full
liquidation. These costs are included in both the optimization and
`Portfolio.total_cost`. For assets absent from `X`, `transaction_costs`
must be a single rate applied to all assets or a dictionary keyed by asset name.
If a float is provided, it is applied to each asset.
If a dictionary is provided, its (key/value) pair must be the
(asset name/asset cost) and the input `X` of the `fit` method must be a
DataFrame with the assets names in columns.
The default value is `0.0`.
.. warning::
Based on the above formula, the periodicity of the transaction costs
must match the periodicity of :math:`\mu`. For example, if the input
`X` is composed of **daily** returns, the `transaction_costs` need to be
expressed as **daily** costs. A transaction cost is paid once per
rebalancing while a position earns its expected return on every period it
is held, so the one-off cost is converted by dividing it by the expected
investment duration (e.g. `0.001 / 21` for a 10 bps cost with daily
returns and a one-month expected holding period).
(See :ref:`Periodicity Convention <periodicity_convention>`)
management_fees : float | dict[str, float] | array-like of shape (n_assets, ), default=0.0
Management fees of the assets. It is used to add linear management fees to the
optimization problem:
.. math:: total\_fee = \sum_{i=1}^{N} f_{i} \times w_{i}
with :math:`f_{i}` the management fee of asset i and :math:`w_{i}` its weight.
The float :math:`total\_fee` is impacting the portfolio expected return in the optimization:
.. math:: expected\_return = \mu^{T} \cdot w - total\_fee
with :math:`\mu` the vector of assets' expected returns and :math:`w` the vector
of assets weights.
If a float is provided, it is applied to each asset.
If a dictionary is provided, its (key/value) pair must be the
(asset name/asset fee) and the input `X` of the `fit` method must be a
DataFrame with the assets names in columns.
The default value is `0.0`.
.. warning::
Based on the above formula, the periodicity of the management fees
must match the periodicity of :math:`\mu`. For example, if the input
`X` is composed of **daily** returns, the `management_fees` need to be
expressed in **daily** fees. Unlike transaction costs, management fees
accrue with holding time, so a stated annual fee converts directly to the
return periodicity (e.g. `0.02 / 252` for a 2% annual fee on daily
returns).
.. note::
Another approach is to directly impact the management fees to the input `X`
in order to express the returns net of fees. However, when estimating the
:math:`\mu` parameter using for example Shrinkage estimators, this approach
would mix a deterministic value with an uncertain one leading to unwanted
bias in the management fees.
previous_weights : float | dict[str, float] | array-like of shape (n_assets, ), optional
Previous weights of the assets. Previous weights are used to compute the
portfolio cost and the portfolio turnover.
For named positions in assets absent from `X`, these calculations assume
full liquidation.
If a float is provided, it is applied to each asset.
If a dictionary is provided, its (key/value) pair must be the
(asset name/asset previous weight) and the input `X` of the `fit` method must
be a DataFrame with the assets names in columns.
The default (`None`) means no previous weights.
Additionally, when `fallback="previous_weights"`, failures will fall back to
these weights if provided.
l1_coef : float, default=0.0
L1 regularization coefficient.
It is used to penalize the objective function by the L1 norm:
.. math:: l1\_coef \times \Vert w \Vert_{1} = l1\_coef \times \sum_{i=1}^{N} |w_{i}|
Increasing this coefficient will reduce the number of non-zero weights
(cardinality). It tends to increase robustness (out-of-sample stability) but
reduces diversification.
The default value is `0.0`.
l2_coef : float, default=0.0
L2 regularization coefficient.
It is used to penalize the objective function by the L2 norm:
.. math:: l2\_coef \times \Vert w \Vert_{2}^{2} = l2\_coef \times \sum_{i=1}^{N} w_{i}^2
It tends to increase robustness (out-of-sample stability).
The default value is `0.0`.
mu_uncertainty_set_estimator : BaseMuUncertaintySet, optional
:ref:`Mu Uncertainty set estimator <uncertainty_set_estimator>`.
If provided, the expected asset returns are modelled with a norm-ball
uncertainty set. It is called worst-case optimization and is a class of robust
optimization. It reduces the instability that arises from the estimation errors
of the expected returns.
The worst-case portfolio expected return is:
.. math:: w^T\hat{\mu} - \kappa_{\mu}\lVert L_{\mu}^Tw\rVert_{q}
with :math:`\kappa` the radius of the uncertainty set (confidence region),
:math:`L` its linear geometry map and :math:`q` the dual norm. For an
ellipsoidal set with shape matrix :math:`S`, :math:`L` is a square-root factor
satisfying :math:`S = L L^T` and :math:`q` is :math:`2`.
The default (`None`) means that no uncertainty set is used.
covariance_uncertainty_set_estimator : BaseCovarianceUncertaintySet, optional
:ref:`Covariance Uncertainty set estimator <uncertainty_set_estimator>`.
If provided, covariance estimation uncertainty is included in the optimized
variance. This approach is known as worst-case optimization, a form of robust
optimization. It reduces sensitivity to covariance estimation errors.
Covariance uncertainty is applied when `risk_measure=RiskMeasure.VARIANCE` or
when `max_variance` is set.
The default (`None`) means that no uncertainty set is used.
linear_constraints : array-like of shape (n_constraints,), optional
Linear constraints on portfolio weights or factor exposures.
Constraint names can reference:
* Asset names: individual asset weights (e.g. `"SPX"`, `"AAPL"`)
* Group names: sums of weights in groups defined by `groups`
* Factor names: portfolio factor exposure (requires factor model prior)
* Factor families: sum of portfolio exposures to all factors in one family
Supported equation patterns include:
* `"name <= value"` or `"name >= value"`
* `"name == value"`
* `"a * name1 + b * name2 <= c * name3 + d"`
For example:
* `"SPX >= 0.10"` --> SPX weight >= 10%
* `"SX5E + SPX >= 0.2"` --> sum of SX5E and SPX weights >= 20%
* `"US == 0.7"` --> sum of weights in US group == 70%
* `"Equity == 3 * Bond"` --> sum of weights in Equity group == 3x sum of weights in Bond group
* `"Momentum <= 0.30"` --> portfolio Momentum exposure <= 30%
* `"style <= 0.50"` --> sum of all style factor exposures (Momentum, Value, Size, etc.) <= 50%
Factor constraints require a prior estimator (e.g.
:class:`~skfolio.prior.TimeSeriesFactorModel`,
:class:`~skfolio.prior.CharacteristicsFactorModel`)
that provides `loading_matrix`, `factor_names` and optionally `factor_families`
in its :class:`~skfolio.prior.FactorModel`.
Asset, group, factor, and factor family names must be unique.
groups : dict[str, list[str]] or array-like of shape (n_groups, n_assets), optional
The assets groups referenced in `linear_constraints`.
If a dictionary is provided, its (key/value) pair must be the
(asset name/asset groups) and the input `X` of the `fit` method must be a
DataFrame with the assets names in columns.
For example:
* `groups = {"SX5E": ["Equity", "Europe"], "SPX": ["Equity", "US"], "TLT": ["Bond", "US"]}`
* `groups = [["Equity", "Equity", "Bond"], ["Europe", "US", "US"]]`
left_inequality : array-like of shape (n_constraints, n_assets), optional
Left inequality matrix :math:`A` of the linear
constraint :math:`A \cdot w \leq b`.
right_inequality : array-like of shape (n_constraints, ), optional
Right inequality vector :math:`b` of the linear
constraint :math:`A \cdot w \leq b`.
risk_free_rate : float, default=0.0
Risk-free rate, expressed in the same frequency as the returns `X` (for
example, :math:`0.04 / 252` for a 4% annual rate with daily returns).
The default value is `0.0`.
min_acceptable_return : float, optional
The minimum acceptable return used to distinguish "downside" and "upside"
returns for the computation of lower partial moments:
* First Lower Partial Moment
* Semi-Variance
* Semi-Deviation
The default (`None`) is to use the mean.
cvar_beta : float, default=0.95
CVaR (Conditional Value at Risk) confidence level in [0, 1].
At 1, CVaR equals the largest loss with positive sample weight.
evar_beta : float, default=0.95
EVaR (Entropic Value at Risk) confidence level in [0, 1].
At 1, EVaR equals the largest loss with positive sample weight.
cdar_beta : float, default=0.95
CDaR (Conditional Drawdown at Risk) confidence level in [0, 1].
At 1, CDaR equals the largest drawdown magnitude with positive sample weight.
edar_beta : float, default=0.95
EDaR (Entropic Drawdown at Risk) confidence level in [0, 1].
At 1, EDaR equals the largest drawdown magnitude with positive sample weight.
add_objective : Callable[[cp.Variable], cp.Expression], optional
Add a custom objective to the existing objective expression.
It is a function that must take as argument the weights `w` and returns a
CVXPY expression.
add_constraints : Callable[[cp.Variable], cp.Expression | list[cp.Expression]], optional
Add a custom constraint or a list of constraints to the existing constraints.
It must be a function taking the CVXPY weight variable `w` as its first
positional argument and, optionally, the estimator instance as its second.
It must return a CVXPY expression or a list of CVXPY expressions, evaluated
when `fit` is called.
For example, to require an effective number of assets of at least 20:
>>> import cvxpy as cp
>>> from skfolio.optimization import MeanRisk
>>> model = MeanRisk(add_constraints=lambda w: cp.sum_squares(w) <= 1 / 20)
The optional second argument gives access to the estimator's attributes,
including quantities estimated during `fit`. For example, to cap each
position size in risk units at 20 bps, using the volatilities estimated
by the prior:
>>> import numpy as np
>>> def position_risk_cap(w, model):
... covariance = model.prior_estimator_.return_distribution_.covariance
... vols = np.sqrt(np.diag(covariance))
... return cp.multiply(vols, w) <= 0.002
>>> model = MeanRisk(add_constraints=position_risk_cap)
overwrite_expected_return : Callable[[cp.Variable], cp.Expression], optional
Overwrite the expected return :math:`\mu \cdot w` with a custom CVXPY
expression. It must be a function taking the CVXPY weight variable `w` as
its first positional argument and, optionally, the estimator instance as
its second. It must return a concave CVXPY expression, evaluated when
`fit` is called. The custom expression replaces the expected return in the
objective function and in the constraints where the expected return is
used.
For example, to adjust the expected return for volatility drag,
approximating the portfolio geometric mean return:
>>> import cvxpy as cp
>>> from skfolio.optimization import MeanRisk
>>> def geometric_expected_return(w, model):
... dist = model.prior_estimator_.return_distribution_
... return dist.mu @ w - 0.5 * cp.quad_form(w, dist.covariance)
>>> model = MeanRisk(overwrite_expected_return=geometric_expected_return)
solver : str, default="CLARABEL"
The solver to use. The default is "CLARABEL" which is written in Rust and has
better numerical stability and performance than ECOS and SCS. Cvxpy will replace
its default solver "ECOS" by "CLARABEL" in future releases.
For more details about available solvers, check the CVXPY documentation:
https://www.cvxpy.org/tutorial/advanced/index.html#choosing-a-solver
solver_params : dict, optional
Solver parameters. For example, `solver_params=dict(verbose=True)`.
The default (`None`) is to use `{"tol_gap_abs": 1e-9, "tol_gap_rel": 1e-9}`
for the solver "CLARABEL" and the CVXPY default otherwise.
For more details about solver arguments, check the CVXPY documentation:
https://www.cvxpy.org/tutorial/solvers
scale_objective : float, optional
Scale each objective element by this value.
It can be used to increase the optimization accuracies in specific cases.
The default (`None`) is set depending on the problem.
scale_constraints : float, optional
Scale each constraint element by this value.
It can be used to increase the optimization accuracies in specific cases.
The default (`None`) is set depending on the problem.
save_problem : bool, default=False
If this is set to True, the CVXPY Problem is saved in `problem_`.
The default is `False`.
portfolio_params : dict, optional
Portfolio parameters forwarded to the resulting `Portfolio` in `predict`.
If not provided and if available on the estimator, the following
attributes are propagated to the portfolio by default: `name`,
`transaction_costs`, `management_fees`, `previous_weights` and `risk_free_rate`.
fallback : BaseOptimization | "previous_weights" | list[BaseOptimization | "previous_weights"], optional
Fallback estimator or a list of estimators to try, in order, when the primary
optimization raises during `fit`. Alternatively, use `"previous_weights"` (alone
or in a list) to fall back to the estimator's `previous_weights`. When a
fallback succeeds, its fitted `weights_` are copied back to the primary
estimator so that `fit` still returns the original instance. For traceability,
`fallback_` stores the successful estimator (or the string `"previous_weights"`)
and `fallback_chain_` stores each attempt with the associated outcome. See
:ref:`optimization_fallbacks`.
raise_on_failure : bool, default=True
Controls error handling when fitting fails and no fallback succeeds. If True,
the estimator raises the final error. If False, the estimator emits a warning
and sets `weights_` to None, so subsequent calls to `predict` return a
:class:`~skfolio.portfolio.FailedPortfolio`. During `fit`, `raise_on_failure`
applies to any fitting error, including errors raised by the prior estimator.
See :ref:`optimization_failure_handling`. When computing multiple portfolios,
setting `raise_on_failure=False` preserves successful allocations and records
each failure separately. See :ref:`optimization_multiple_results`.
Attributes
----------
weights_ : ndarray of shape (n_assets,) or (n_optimizations, n_assets)
Weights of the assets.
problem_values_ : dict[str, float] | list[dict[str, float]] of size n_optimizations
Expression values retrieved from the CVXPY problem.
prior_estimator_ : BasePrior
Fitted `prior_estimator`.
mu_uncertainty_set_estimator_ : BaseMuUncertaintySet
Fitted `mu_uncertainty_set_estimator` if provided.
covariance_uncertainty_set_estimator_ : BaseCovarianceUncertaintySet
Fitted `covariance_uncertainty_set_estimator` if provided.
problem_: cvxpy.Problem
CVXPY problem used for the optimization. Only when `save_problem` is set to
`True`.
fallback_ : BaseOptimization | "previous_weights" | None
The fallback estimator instance, or the string `"previous_weights"`, that
produced the final result. `None` if no fallback was used.
fallback_chain_ : list[tuple[str, str]] | None
Sequence describing the optimization fallback attempts. Each element is a
pair `(estimator_repr, outcome)` where `estimator_repr` is the string
representation of the primary estimator or a fallback (e.g. `"EqualWeighted()"`,
`"previous_weights"`), and `outcome` is `"success"` if that step produced
a valid solution, otherwise the stringified error message. For successful
fits without any fallback, this is `None`.
error_ : str | list[str | None] | None
For a single portfolio, this is the recorded error message, or None after a
successful allocation or fallback. For multiple portfolios, it is a list with
one entry per row of `weights_`, containing an error message for each failed
portfolio and None for each successful portfolio.
Notes
-----
All estimators should specify all parameters as explicit keyword arguments in
`__init__` (no `*args` or `**kwargs`), following scikit-learn conventions.
"""
_solver_params: dict
_scale_objective: cp.Constant
_scale_constraints: cp.Constant
_cvx_cache: dict
problem_: cp.Problem
problem_values_: dict[str, float] | list[dict[str, float] | None] | None
prior_estimator_: BasePrior
mu_uncertainty_set_estimator_: BaseMuUncertaintySet
covariance_uncertainty_set_estimator_: BaseCovarianceUncertaintySet
@abstractmethod
def __init__(
self,
risk_measure: RiskMeasure = RiskMeasure.VARIANCE,
prior_estimator: BasePrior | None = None,
min_weights: skt.MultiInput | None = 0.0,
max_weights: skt.MultiInput | None = 1.0,
budget: float | None = 1.0,
min_budget: float | None = None,
max_budget: float | None = None,
max_short: float | None = None,
max_long: float | None = None,
cardinality: int | None = None,
group_cardinalities: dict[str, int] | None = None,
threshold_long: skt.MultiInput | None = None,
threshold_short: skt.MultiInput | None = None,
transaction_costs: skt.MultiInput = 0.0,
management_fees: skt.MultiInput = 0.0,
previous_weights: skt.MultiInput | None = None,
target_weights: skt.MultiInput | None = None,
groups: skt.Groups | None = None,
linear_constraints: skt.LinearConstraints | None = None,
left_inequality: skt.Inequality | None = None,
right_inequality: skt.Inequality | None = None,
l1_coef: float = 0.0,
l2_coef: float = 0.0,
mu_uncertainty_set_estimator: BaseMuUncertaintySet | None = None,
covariance_uncertainty_set_estimator: (
BaseCovarianceUncertaintySet | None
) = None,
risk_free_rate: float = 0.0,
min_acceptable_return: skt.Target | None = None,
cvar_beta: float = 0.95,
evar_beta: float = 0.95,
cdar_beta: float = 0.95,
edar_beta: float = 0.95,
solver: str = "CLARABEL",
solver_params: dict | None = None,
scale_objective: float | None = None,
scale_constraints: float | None = None,
save_problem: bool = False,
add_objective: skt.ExpressionFunction | None = None,
add_constraints: skt.ConstraintFunction | None = None,
overwrite_expected_return: skt.ExpressionFunction | None = None,
portfolio_params: dict | None = None,
fallback: skt.Fallback = None,
raise_on_failure: bool = True,
) -> None:
super().__init__(
previous_weights=previous_weights,
portfolio_params=portfolio_params,
fallback=fallback,
raise_on_failure=raise_on_failure,
)
if risk_measure.is_annualized:
warnings.warn(
f"The annualized risk measure {risk_measure} will be converted"
f"to its non-annualized version {risk_measure.non_annualized_measure}",
stacklevel=2,
)
risk_measure = risk_measure.non_annualized_measure
self.risk_measure = risk_measure
self.prior_estimator = prior_estimator
self.mu_uncertainty_set_estimator = mu_uncertainty_set_estimator
self.covariance_uncertainty_set_estimator = covariance_uncertainty_set_estimator
self.min_weights = min_weights
self.max_weights = max_weights
self.budget = budget
self.min_budget = min_budget
self.max_budget = max_budget
self.max_short = max_short
self.max_long = max_long
self.cardinality = cardinality
self.group_cardinalities = group_cardinalities
self.threshold_long = threshold_long
self.threshold_short = threshold_short
self.min_acceptable_return = min_acceptable_return
self.transaction_costs = transaction_costs
self.management_fees = management_fees
self.target_weights = target_weights
self.groups = groups
self.linear_constraints = linear_constraints
self.left_inequality = left_inequality
self.right_inequality = right_inequality
self.l1_coef = l1_coef
self.l2_coef = l2_coef
self.risk_free_rate = risk_free_rate
self.add_objective = add_objective
self.add_constraints = add_constraints
self.overwrite_expected_return = overwrite_expected_return
self.solver = solver
self.solver_params = solver_params
self.save_problem = save_problem
self.scale_objective = scale_objective
self.scale_constraints = scale_constraints
self.cvar_beta = cvar_beta
self.evar_beta = evar_beta
self.cdar_beta = cdar_beta
self.edar_beta = edar_beta
self._clear_models_cache()
def _call_custom_func(
self,
func: Callable[..., _ResultT],
w: cp.Variable,
name: str = "custom_func",
) -> _ResultT:
"""Call a user specific function, infer arguments and perform validation.
Parameters
----------
func : Callable[[cvxpy Variable, any], cvxpy Expression]
The custom function. Must have one or two positional arguments.
The first argument is the CVXPY weight variable `w` and the second is
the reference to the class itself.
w : cvxpy Variable
The CVXPY Variable representing assets weights.
Returns
-------
result : cvxpy Expression | list[cvxpy Expression]
Result of calling the custom function.
"""
try:
func_code = func.__code__ # ty: ignore[unresolved-attribute]
except AttributeError as err:
raise ValueError("Custom functions is invalid") from err
if func_code.co_argcount == 1:
args = (w,)
elif func_code.co_argcount == 2:
args = (w, self)
else:
raise ValueError(
"Custom functions must have 1 or 2 positional arguments, got"
f" {func_code.co_argcount}"
)
try:
return func(*args)
except Exception as err:
raise TypeError(
f"Error while calling {name}. "
f"{name} must be a function taking as argument "
"the weight variable OR the weight variable and the estimator object."
) from err
def _clear_models_cache(self) -> None:
"""Clear the cache of CVX models."""
self._cvx_cache = {}
def _get_weight_constraints(
self,
n_assets: int,
w: cp.Variable,
factor: skt.Factor,
allow_negative_weights: bool = True,
return_distribution: ReturnDistribution | None = None,
) -> list[cpc.Constraint]:
"""Compute weight constraints from input parameters.
Parameters
----------
n_assets : int
Number of investable assets.
w : cvxpy Variable
The CVXPY Variable representing assets weights.
factor : cvxpy Variable | cvxpy Constant
Cvxpy variable or constant.
allow_negative_weights : bool, default=True
Whether to allow negative weights.
return_distribution : ReturnDistribution, optional
The return distribution containing the factor model data (loading_matrix,
factor_names, factor_families) for factor constraints.
Returns
-------
constraints : list[cvxpy Constraint]
The list of weight constraints.
"""
constraints = []
# Clean and convert to array
min_weights = self.min_weights
max_weights = self.max_weights
threshold_long = self.threshold_long
threshold_short = self.threshold_short
groups = self.groups
assets_names = getattr(self, "feature_names_in_", None)
investable_mask = getattr(self, "investable_mask_", None)
if min_weights is not None:
min_weights = self._clean_input(
min_weights,
n_assets=n_assets,
fill_value=0,
name="min_weights",
)
if max_weights is not None:
max_weights = self._clean_input(
max_weights,
n_assets=n_assets,
fill_value=1,
name="max_weights",
)
if threshold_long is not None:
threshold_long = self._clean_input(
threshold_long,
n_assets=n_assets,
fill_value=0,
name="threshold_long",
)
if np.all(threshold_long == 0):
threshold_long = None
if threshold_short is not None:
threshold_short = self._clean_input(
threshold_short,
n_assets=n_assets,
fill_value=0,
name="threshold_short",
)
if np.all(threshold_short == 0):
threshold_short = None
if groups is not None:
groups = input_to_array(
items=groups,
n_assets=n_assets,
fill_value="",
dim=2,
assets_names=assets_names,
investable_mask=investable_mask,
name="groups",
)
is_mip = (
(self.cardinality is not None and self.cardinality < n_assets)
or (self.group_cardinalities is not None)
or threshold_long is not None
or threshold_short is not None
)
if is_mip and self.solver not in MI_SOLVERS:
raise ValueError(
"You are using constraints that require a mixed-integer solver and "
f"{self.solver} doesn't support MIP problems. For an open-source "
"option, we recommend using SCIP by setting `solver='SCIP'`. "
"To install it, use: `pip install cvxpy[SCIP]`. For commercial "
"solvers, supported options include MOSEK, GUROBI, or CPLEX."
)
constraints += self._weight_and_exposure_constraints(
w=w,
factor=factor,
min_weights=min_weights,
max_weights=max_weights,
allow_negative_weights=allow_negative_weights,
)
constraints += self._budget_constraints(w=w, factor=factor)
if is_mip:
constraints += self._mixed_integer_constraints(
n_assets=n_assets,
w=w,
factor=factor,
min_weights=min_weights,
max_weights=max_weights,
threshold_long=threshold_long,
threshold_short=threshold_short,
groups=groups,
)
constraints += self._linear_constraints(
w=w,
factor=factor,
groups=groups,
assets_names=assets_names,
investable_mask=investable_mask,
return_distribution=return_distribution,
)
constraints += self._matrix_inequality_constraints(
n_assets=n_assets, w=w, factor=factor, investable_mask=investable_mask
)
return constraints
def _weight_and_exposure_constraints(
self,
w: cp.Variable,
factor: skt.Factor,
min_weights: float | FloatArray | None,
max_weights: float | FloatArray | None,
allow_negative_weights: bool,
) -> list[cpc.Constraint]:
"""Constrain individual asset weights and total long and short exposure.
Weight bounds must already be aligned with the investable assets.
"""
constraints = []
if min_weights is not None:
if not allow_negative_weights and np.any(min_weights < 0):
raise ValueError(
f"{self.__class__.__name__} must have non negative `min_weights` "
f"constraint otherwise the problem becomes non-convex."
)
constraints.append(
w * self._scale_constraints
>= min_weights * factor * self._scale_constraints
)
if max_weights is not None:
constraints.append(
w * self._scale_constraints
<= max_weights * factor * self._scale_constraints
)
if self.max_long is not None:
max_long = float(self.max_long)
if max_long <= 0:
raise ValueError("`max_long` must be strictly positive")
constraints.append(
cp.sum(cp.pos(w)) * self._scale_constraints
<= max_long * factor * self._scale_constraints
)
if self.max_short is not None:
max_short = float(self.max_short)
if max_short <= 0:
raise ValueError("`max_short` must be strictly positive")
constraints.append(
cp.sum(cp.neg(w)) * self._scale_constraints
<= max_short * factor * self._scale_constraints
)
return constraints
def _budget_constraints(
self, w: cp.Variable, factor: skt.Factor
) -> list[cpc.Constraint]:
"""Constrain the sum of weights using the configured budget or limits."""
constraints = []
if self.min_budget is not None:
constraints.append(
cp.sum(w) * self._scale_constraints
>= float(self.min_budget) * factor * self._scale_constraints
)
if self.max_budget is not None:
constraints.append(
cp.sum(w) * self._scale_constraints
<= float(self.max_budget) * factor * self._scale_constraints
)
if self.budget is not None:
if self.max_budget is not None:
raise ValueError(
"`max_budget`and `budget` cannot be provided at the same time"
)
if self.min_budget is not None:
raise ValueError(
"`min_budget`and `budget` cannot be provided at the same time"
)
constraints.append(
cp.sum(w) * self._scale_constraints
== float(self.budget) * factor * self._scale_constraints
)
return constraints
def _mixed_integer_constraints(
self,
n_assets: int,
w: cp.Variable,
factor: skt.Factor,
min_weights: float | FloatArray | None,
max_weights: float | FloatArray | None,
threshold_long: float | FloatArray | None,
threshold_short: float | FloatArray | None,
groups: AnyArray | None,
) -> list[cpc.Constraint]:
"""Build cardinality and position-threshold constraints.
Bounds, thresholds, and groups must already be aligned with the investable
assets. All-zero thresholds must be converted to `None`. The caller checks
that the solver supports mixed-integer problems.
"""
if max_weights is None or min_weights is None:
raise ValueError(
"'max_weights' and 'min_weights' must be provided with cardinality "
"constraint"
)
is_short = np.any(min_weights < 0)
if np.all(min_weights > 0):
raise ValueError(
"Cardinality and Threshold constraint can only be applied "
"if 'min_weights' are not all strictly positive (you allow some "
"weights to be 0)"
)
if self.group_cardinalities is not None and groups is None:
raise ValueError(
"When 'group_cardinalities' is provided, you must also "
"also provide 'groups'"
)
if threshold_long is not None and threshold_short is None and is_short:
raise ValueError(
"When 'threshold_long' is provided and 'min_weights' can be negative "
"(short positions are allowed), then 'threshold_short' must also be "
"provided"
)
if threshold_short is not None and threshold_long is None:
raise ValueError(
"When 'threshold_short' is provided, 'threshold_long' must also be "
"provided"
)
if threshold_short is not None and is_short:
return _mip_weight_constraints_threshold_short(
n_assets=n_assets,
w=w,
factor=factor,
scale_constraints=self._scale_constraints,
cardinality=self.cardinality,
group_cardinalities=self.group_cardinalities,
max_weights=max_weights,
groups=groups,
min_weights=min_weights,
threshold_long=threshold_long, # ty: ignore[invalid-argument-type]
threshold_short=threshold_short,
)
return _mip_weight_constraints_no_short_threshold(
n_assets=n_assets,
w=w,
factor=factor,
scale_constraints=self._scale_constraints,
cardinality=self.cardinality,
group_cardinalities=self.group_cardinalities,
max_weights=max_weights,
groups=groups,
min_weights=min_weights,
threshold_long=threshold_long,
)
def _linear_constraints(
self,
w: cp.Variable,
factor: skt.Factor,
groups: AnyArray | None,
assets_names: StrArray | None,
investable_mask: BoolArray | None,
return_distribution: ReturnDistribution | None,
) -> list[cpc.Constraint]:
"""Build equalities and inequalities from `linear_constraints`.
Supplied groups must already be aligned with the investable assets.
Factor constraints use the factor model from `return_distribution`.
"""
constraints = []
if self.linear_constraints is not None:
if groups is None:
if assets_names is None:
raise ValueError(
"If `linear_constraints` is provided you must provide either"
" `groups` or `X` as a DataFrame with asset names in columns"
)
if investable_mask is None:
groups = np.asarray([assets_names])
else:
groups = np.asarray([assets_names[investable_mask]])
# Extract factor info from return_distribution for factor constraints
loading_matrix = None
factor_groups = None
if return_distribution is not None:
factor_model = return_distribution.factor_model
if factor_model is not None:
loading_matrix = factor_model.loading_matrix
if factor_model.factor_families is not None:
factor_groups = np.array(
[
factor_model.factor_names,
factor_model.factor_families,
]
)
else:
factor_groups = np.array([factor_model.factor_names])
a_eq, b_eq, a_ineq, b_ineq = equations_to_matrix(
groups=groups,
equations=self.linear_constraints,
raise_if_group_missing=False,
loading_matrix=loading_matrix,
factor_groups=factor_groups,
)
if len(a_eq) != 0:
constraints.append(
a_eq @ w * self._scale_constraints
- b_eq * factor * self._scale_constraints
== 0
)
if len(a_ineq) != 0:
constraints.append(
a_ineq @ w * self._scale_constraints
- b_ineq * factor * self._scale_constraints
<= 0
)
return constraints
def _matrix_inequality_constraints(
self,
n_assets: int,
w: cp.Variable,
factor: skt.Factor,
investable_mask: BoolArray | None,
) -> list[cpc.Constraint]:
"""Build `left_inequality @ w <= right_inequality * factor` constraints.
When `investable_mask` is provided, the matrix must include all original
assets, including non-investable ones.
"""
constraints = []
if self.left_inequality is not None and self.right_inequality is not None:
left_inequality = np.asarray(self.left_inequality)
right_inequality = np.asarray(self.right_inequality)
if left_inequality.ndim != 2:
raise ValueError(
f"`left_inequality` must be a 2D array, got {left_inequality.ndim}D"
" array"
)
if right_inequality.ndim != 1:
raise ValueError(
"`right_inequality` must be a 1D array, got"
f" {right_inequality.ndim}D array"
)
if investable_mask is not None:
n_total_assets = len(investable_mask)
if left_inequality.shape[1] != n_total_assets:
raise ValueError(
"`left_inequality` must be of shape (n_inequalities, n_total_assets) "
f"with n_total_assets={n_total_assets}, got {left_inequality.shape[1]}"
)
left_inequality = left_inequality[:, investable_mask]
if left_inequality.shape[1] != n_assets:
raise ValueError(
"`left_inequality` must be of shape (n_inequalities, n_assets) "
f"with n_assets={n_assets}, got {left_inequality.shape[1]}"
)
if left_inequality.shape[0] != right_inequality.shape[0]:
raise ValueError(
"`left_inequality` and `right_inequality` must have same number of"
f" rows (i.e. n_inequalities) , got {left_inequality.shape[0]} and"
f" {right_inequality.shape[0]}"
)
constraints.append(
left_inequality @ w * self._scale_constraints
- right_inequality * factor * self._scale_constraints
<= 0
)
return constraints
def _set_solver_params(self, default: dict | None) -> None:
"""Set the solver params by saving its value in `_solver_params`.
It uses `solver` if provided otherwise it uses the `default` solver.
Parameters
----------
default : str
The default solver params to use when `solver_params` is `None`.
"""
if self.solver_params is None:
self._solver_params = default if default is not None else {}
else:
self._solver_params = self.solver_params
def _set_scale_objective(self, default: float) -> None:
"""Set the objective scale by saving its value in `_scale_objective`.
It uses `scale_objective` if provided otherwise it uses the `default` scale.
Parameters
----------
default : float
The default objective scale to use when `scale_objective` is `None`.
"""
if self.scale_objective is None:
self._scale_objective = cp.Constant(default)
else:
self._scale_objective = cp.Constant(self.scale_objective)
def _set_scale_constraints(self, default: float) -> None:
"""Set the constraints scale by saving its value in `_scale_constraints`.
It uses `scale_constraints` if provided otherwise it uses the `default` scale.
Parameters
----------
default : float
The default constraints scale to use when `scale_constraints` is `None`.
"""
if self.scale_constraints is None:
self._scale_constraints = cp.Constant(default)
else:
self._scale_constraints = cp.Constant(self.scale_constraints)
def _get_custom_objective(self, w: cp.Variable) -> cp.Expression:
"""Return the CVXPY expression evaluated by calling the `add_objective`
function if provided, otherwise returns the CVXPY constant `0`.
Parameters
----------
w : cvxpy Variable
The CVXPY Variable representing assets weights.
Returns
-------
expression : cvxpy Expression
The CVXPY expression evaluated by calling the `add_objective`
function if provided, otherwise returns the CVXPY constant `0`.
"""
if self.add_objective is None:
return cp.Constant(0)
return self._call_custom_func(
func=self.add_objective, w=w, name="add_objective"
)
def _get_custom_constraints(self, w: cp.Variable) -> list[cpc.Constraint]:
"""Return the list of CVXPY expressions evaluated by calling the
`add_constraint`s function if provided, otherwise returns an empty list.
Parameters
----------
w : cvxpy Variable
The CVXPY Variable representing assets weights.
Returns
-------
expressions : list of cvxpy Expression
The list of CVXPY expressions evaluated by calling the
`add_constraints` function if provided, otherwise returns an empty list.
"""
if self.add_constraints is None:
return []
constraints = self._call_custom_func(
func=self.add_constraints, w=w, name="add_constraint"
)
if isinstance(constraints, list):
return constraints
return [constraints]
@cache_method("_cvx_cache")
def _cvx_expected_return(
self, return_distribution: ReturnDistribution, w: cp.Variable
) -> cp.Expression:
"""Expected Return expression."""
if self.overwrite_expected_return is None:
expected_return = return_distribution.mu @ w
else:
expected_return = self._call_custom_func(
func=self.overwrite_expected_return,
w=w,
name="overwrite_expected_return",
)
return expected_return
# Model reused among multiple risk measure
def _solve_problem(
self,
problem: cp.Problem,
w: cp.Variable,
factor: skt.Factor,
parameters_values: skt.ParametersValues | None = None,
expressions: dict[str, cp.Expression] | None = None,
) -> None:
"""Solve the CVXPY problem and store `weights_` and `problem_values_`.
Reuse the problem across parameter values when solving a sweep. Save it in
`problem_` when `save_problem` is enabled.
Parameters
----------
problem : cvxpy Problem
Optimization problem containing the objective and constraints.
w : cvxpy Variable of shape (n_assets,)
Weights of the investable assets in optimization units.
factor : cvxpy Variable | cvxpy Constant
Scalar homogenization factor. Final portfolio weights are `w / factor`.
parameters_values : list[tuple[cvxpy Parameter, float | ndarray]], optional
Parameters and the values to assign before solving. Arrays are paired by
index and must be nonempty with the same length. Scalar values are reused
for every solve. None gives a single solve with the current parameters.
expressions : dict[str, cvxpy Expression] | None, optional
CVXPY expressions in reporting units, evaluated after each solve and saved
in `problem_values_`. The `"objective"` entry contains the solver objective
divided by `_scale_objective`.
Raises
------
ValueError
If the solver is unavailable or parameter arrays are empty or have
different lengths.
ConvexOptimizationError
If a single solve fails, all targets in a sweep fail, or a target fails
while `raise_on_failure` is True.
Warns
-----
UserWarning
If a solution may be inaccurate or a sweep continues after a failed target.
Notes
-----
With `raise_on_failure=False`, a failed target in a sweep has NaN weights,
None for its reported values, and an error message in `error_`. The saved
problem retains the state left by the last solve attempt.
"""
if self.solver not in INSTALLED_SOLVERS:
raise ValueError(f"The solver {self.solver} is not installed.")
if parameters_values is None:
parameters_values = []
expressions = {
**(expressions or {}),
"objective": problem.objective.expr / self._scale_objective,
}
# If the parameter value is a list, each element is the parameter value of
# a distinct optimization. Therefore, each list must have same length.
values_by_parameter = [
(parameter, np.asarray(values, dtype=float))
for parameter, values in parameters_values
]
sizes = [len(values) for _, values in values_by_parameter if values.ndim != 0]
n_optimizations = sizes[0] if sizes else 1
if n_optimizations == 0 or any(size != n_optimizations for size in sizes):
raise ValueError(
"Array values in `parameters_values` must be nonempty and have the same length."
)
# Scalar parameter values will be used in each optimization, therefore we
# broadcast them.
values_by_parameter = [
(
parameter,
np.full(n_optimizations, values) if values.ndim == 0 else values,
)
for parameter, values in values_by_parameter
]
self.problem_values_ = None
try:
if n_optimizations == 1:
for parameter, values in values_by_parameter:
parameter.value = values[0]
weights, self.problem_values_ = _solve(
w=w,
factor=factor,
expressions=expressions,
problem=problem,
solver=self.solver,
solver_params=self._solver_params,
)
self.weights_ = self._expand_weights_to_full_universe(weights=weights)
else:
all_weights = []
all_problem_values = []
all_errors = []
with warnings.catch_warnings():
warnings.simplefilter("once", UserWarning)
for i in range(n_optimizations):
for parameter, values in values_by_parameter:
parameter.value = values[i]
try:
weights, problem_values = _solve(
w=w,
factor=factor,
expressions=expressions,
problem=problem,
solver=self.solver,
solver_params=self._solver_params,
)
error = None
except OptimizationError as optimization_error:
if self.raise_on_failure:
raise
error = str(optimization_error)
warnings.warn(error, stacklevel=2)
problem_values = None
weights = np.full(w.shape, np.nan, dtype=float)
all_problem_values.append(problem_values)
all_weights.append(weights)
all_errors.append(error)
all_weights = np.array(all_weights, dtype=float)
failed = np.array([error is not None for error in all_errors])
if failed.all():
raise ConvexOptimizationError(
f"All {n_optimizations} optimizations failed, with last optimization error {all_errors[-1]}"
)
self.weights_ = self._expand_weights_to_full_universe(
weights=all_weights
)
# Failed portfolios remain all-NaN, including unavailable assets.
self.weights_[failed] = np.nan
self.problem_values_ = all_problem_values
self.error_ = all_errors
finally:
if self.save_problem:
self.problem_ = problem
self._clear_models_cache()
@cache_method("_cvx_cache")
def _cvx_mu_uncertainty_set(
self, mu_uncertainty_set: UncertaintySet, w: cp.Variable
) -> cp.Expression:
"""Uncertainty Set expression of expected returns.
Parameters
----------
mu_uncertainty_set : UncertaintySet
The uncertainty set model of expected returns.
w : cvxpy Variable
The CVXPY Variable representing assets weights.
Returns
-------
expression : cvxpy Expression
The CVXPY Expression of the uncertainty set of expected returns.
"""
# cvxpy annotates `p` as int or str but accepts any float, including inf.
return mu_uncertainty_set.radius * cp.pnorm(
mu_uncertainty_set.geometry.T @ w,
mu_uncertainty_set.dual_norm, # ty: ignore[invalid-argument-type]
)
@cache_method("_cvx_cache")
def _cvx_regularization(self, w: cp.Variable) -> cp.Expression:
"""L1 and L2 regularization expression.
Parameters
----------
w : cvxpy Variable
The CVXPY Variable representing assets weights.
Returns
-------
expression : cvxpy Expression
The CVXPY Expression of L1 and L2 regularization.
"""
# Norm L1
if self.l1_coef is None or self.l1_coef == 0:
l1_reg = cp.Constant(0)
else:
l1_reg = cp.Constant(self.l1_coef) * cp.norm(w, 1)
# Norm L2
if self.l2_coef is None or self.l2_coef == 0:
l2_reg = cp.Constant(0)
else:
l2_reg = self.l2_coef * cp.sum_squares(w)
regularization = l1_reg + l2_reg
return regularization
@cache_method("_cvx_cache")
def _cvx_transaction_cost(
self,
return_distribution: ReturnDistribution,
w: cp.Variable,
factor: skt.Factor,
) -> cp.Expression:
"""Transaction cost expression.
Parameters
----------
return_distribution : ReturnDistribution
asset returns distribution DataModel.
w : cvxpy Variable
The CVXPY Variable representing assets weights.
factor : cvxpy Variable | cvxpy Constant
Additional variable used for the optimization of some objective function
like the ratio maximization.
Returns
-------
expression : cvxpy Expression
The CVXPY Expression of transaction cost.
"""
n_assets = return_distribution.returns.shape[1]
transaction_costs = self._clean_input(
self.transaction_costs,
n_assets=n_assets,
fill_value=0,
name=_TRANSACTION_COSTS,
)
if np.all(transaction_costs == 0):
cost = cp.Constant(0)
else:
previous_weights = self._clean_previous_weights(n_assets=n_assets)
if np.isscalar(transaction_costs):
cost = transaction_costs * cp.norm(previous_weights * factor - w, 1)
else:
cost = cp.norm(
cp.multiply(transaction_costs, (previous_weights * factor - w)), 1
)
_, liquidation_cost = _get_liquidation_turnover_and_cost(
previous_weights=self.previous_weights,
transaction_costs=self.transaction_costs,
assets_names=getattr(self, "feature_names_in_", None),
investable_mask=getattr(self, "investable_mask_", None),
)
if liquidation_cost:
cost += liquidation_cost * factor
return cost
@cache_method("_cvx_cache")
def _cvx_management_fee(
self, return_distribution: ReturnDistribution, w: cp.Variable
) -> cp.Expression:
"""Management fee expression.
Parameters
----------
return_distribution : ReturnDistribution
asset returns distribution DataModel.
w : cvxpy Variable
The CVXPY Variable representing assets weights.
Returns
-------
expression : cvxpy Expression
The CVXPY Expression of management fee .
"""
n_assets = return_distribution.returns.shape[1]
management_fees = self._clean_input(
self.management_fees,
n_assets=n_assets,
fill_value=0,
name=_MANAGEMENT_FEES,
)
if np.all(management_fees == 0):
return cp.Constant(0)
if np.isscalar(management_fees):
management_fees *= np.ones(n_assets)
return management_fees @ w
@cache_method("_cvx_cache")
def _cvx_returns(
self, return_distribution: ReturnDistribution, w: cp.Variable
) -> cp.Expression:
"""Expression of the portfolio returns series.
Parameters
----------
return_distribution : ReturnDistribution
asset returns distribution DataModel.
w : cvxpy Variable
The CVXPY Variable representing assets weights.
Returns
-------
expression : cvxpy Expression
The CVXPY Expression the portfolio returns series.
"""
returns = return_distribution.returns @ w
return returns
@cache_method("_cvx_cache")
def _turnover(
self, n_assets: int, w: cp.Variable, factor: skt.Factor
) -> cp.Expression:
"""Per-asset turnover in the investable optimization universe.
Parameters
----------
n_assets : int
The number of investable assets.
w : cvxpy Variable
The CVXPY Variable representing assets weights.
factor : cvxpy Variable | cvxpy Constant
Additional variable used for the optimization of some objective function
like the ratio maximization.
Returns
-------
expression : cvxpy Expression
The CVXPY Expression the portfolio turnover.
"""
previous_weights = self._clean_previous_weights(n_assets=n_assets)
turnover = cp.abs(w - previous_weights * factor)
return turnover
@cache_method("_cvx_cache")
def _cvx_min_acceptable_return(
self,
return_distribution: ReturnDistribution,
w: cp.Variable,
min_acceptable_return: skt.Target | None = None,
) -> cp.Expression:
"""Expression of the portfolio Minimum Acceptable Returns.
Parameters
----------
return_distribution : ReturnDistribution
asset returns distribution DataModel..
w : cvxpy Variable
The CVXPY Variable representing assets weights.
min_acceptable_return : float | ndarray of shape (n_assets,)
The minimum acceptable return used to distinguish "downside" and "upside"
returns for the computation of lower partial moments.
Returns
-------
expression : cvxpy Expression
The CVXPY Expression the portfolio Minimum Acceptable Returns.
"""
if min_acceptable_return is None:
min_acceptable_return = return_distribution.mu
if isinstance(
min_acceptable_return, np.ndarray
) and min_acceptable_return.shape != (len(min_acceptable_return), 1):
min_acceptable_return = min_acceptable_return[np.newaxis, :]
mar = (return_distribution.returns - min_acceptable_return) @ w
return mar
@cache_method("_cvx_cache")
def __cvx_drawdown(
self,
return_distribution: ReturnDistribution,
w: cp.Variable,
factor: skt.Factor,
) -> tuple[cp.Variable, list[cpc.Constraint]]:
"""Expression of the portfolio drawdown.
Parameters
----------
return_distribution : ReturnDistribution
asset returns distribution DataModel.
w : cvxpy Variable
The CVXPY Variable representing assets weights.
factor : cvxpy Variable | cvxpy Constant
Additional variable used for the optimization of some objective function
like the ratio maximization.
Returns
-------
expression : cvxpy Expression
The CVXPY Expression the portfolio drawdown.
"""
n_observations = return_distribution.returns.shape[0]
ptf_returns = self._cvx_returns(return_distribution=return_distribution, w=w)
ptf_transaction_cost = self._cvx_transaction_cost(
return_distribution=return_distribution, w=w, factor=factor
)
ptf_management_fee = self._cvx_management_fee(
return_distribution=return_distribution, w=w
)
v = cp.Variable(n_observations + 1)
constraints = [
v[1:] * self._scale_constraints
>= v[:-1] * self._scale_constraints
- ptf_returns * self._scale_constraints
+ ptf_transaction_cost * self._scale_constraints
+ ptf_management_fee * self._scale_constraints,
v[1:] * self._scale_constraints >= 0,
v[0] * self._scale_constraints == 0,
]
return v, constraints
def _cvx_drawdown(
self,
return_distribution: ReturnDistribution,
w: cp.Variable,
factor: skt.Factor,
) -> tuple[cp.Variable, list[cpc.Constraint]]:
"""Expression of the portfolio drawdown.
Wrapper around __cvx_drawdown to avoid re-adding the constraints when they
have already been included in the problem.
Parameters
----------
return_distribution : ReturnDistribution
asset returns distribution DataModel.
w : cvxpy Variable
The CVXPY Variable representing assets weights.
factor : cvxpy Variable | cvxpy Constant
Additional variable used for the optimization of some objective function
like the ratio maximization.
Returns
-------
expression : cvxpy Expression
The CVXPY Expression the portfolio drawdown.
"""
if "__cvx_drawdown" in self._cvx_cache:
v, _ = self.__cvx_drawdown(
return_distribution=return_distribution, w=w, factor=factor
)
return v, []
return self.__cvx_drawdown(
return_distribution=return_distribution, w=w, factor=factor
)
def _tracking_error(
self,
return_distribution: ReturnDistribution,
w: cp.Variable,
y: FloatArray,
factor: skt.Factor,
) -> cp.Expression:
"""Expression of the portfolio tracking error.
Parameters
----------
return_distribution : ReturnDistribution
asset returns distribution DataModel.
w : cvxpy Variable
The CVXPY Variable representing assets weights.
y : ndarray of shape (n_observations,)
Benchmark for the tracking error computation.
factor : cvxpy Variable | cvxpy Constant
Additional variable used for the optimization of some objective function
like the ratio maximization.
Returns
-------
expression : cvxpy Expression
The CVXPY Expression the portfolio tracking error.
"""
n_observations = return_distribution.returns.shape[0]
ptf_returns = self._cvx_returns(return_distribution=return_distribution, w=w)
tracking_error = cp.norm(ptf_returns - y * factor, "fro") / cp.sqrt(
n_observations - 1
)
return tracking_error
# Risk Measures risk models
# They need to be named f'_{risk_measure}_risk' as they are loaded dynamically in
# mean_risk_optimization()
def _mean_absolute_deviation_risk(
self,
return_distribution: ReturnDistribution,
w: cp.Variable,
min_acceptable_return: skt.Target,
) -> skt.RiskResult:
"""Expression and Constraints of the Mean Absolute Deviation risk measure.
Parameters
----------
return_distribution : ReturnDistribution
asset returns distribution DataModel.
w : cvxpy Variable
The CVXPY Variable representing assets weights.
min_acceptable_return : float | ndarray of shape (n_assets,)
The minimum acceptable return used to distinguish "downside" and "upside"
returns for the computation of lower partial moments.
Returns
-------
expression : tuple[cvxpy Expression , list[cvxpy Expression]]
CVXPY expression and constraints of the mean absolute deviation risk
measure.
"""
n_observations = return_distribution.returns.shape[0]
ptf_min_acceptable_return = self._cvx_min_acceptable_return(
return_distribution=return_distribution,
w=w,
min_acceptable_return=min_acceptable_return,
)
v = cp.Variable(n_observations, nonneg=True)
if return_distribution.sample_weight is None:
risk = 2 * cp.sum(v) / n_observations
else:
risk = 2 * cp.sum(cp.multiply(return_distribution.sample_weight, v))
constraints = [
ptf_min_acceptable_return * self._scale_constraints
>= -v * self._scale_constraints
]
return risk, constraints
def _first_lower_partial_moment_risk(
self,
return_distribution: ReturnDistribution,
w: cp.Variable,
min_acceptable_return: skt.Target,
factor: skt.Factor,
) -> skt.RiskResult:
"""Expression and Constraints of the First Lower Partial Moment risk measure.
Parameters
----------
return_distribution : ReturnDistribution
asset returns distribution DataModel.
w : cvxpy Variable
The CVXPY Variable representing assets weights.
min_acceptable_return : float | ndarray of shape (n_assets,)
The minimum acceptable return used to distinguish "downside" and "upside"
returns for the computation of lower partial moments.
factor : cvxpy Variable | cvxpy Constant
Additional variable used for the optimization of some objective function
like the ratio maximization.
Returns
-------
expression : tuple[cvxpy Expression , list[cvxpy Expression]]
CVXPY expression and constraints of the first lower partial moment risk
measure.
"""
n_observations = return_distribution.returns.shape[0]
ptf_min_acceptable_return = self._cvx_min_acceptable_return(
return_distribution=return_distribution,
w=w,
min_acceptable_return=min_acceptable_return,
)
v = cp.Variable(n_observations, nonneg=True)
if return_distribution.sample_weight is None:
risk = cp.sum(v) / n_observations
else:
risk = cp.sum(cp.multiply(return_distribution.sample_weight, v))
constraints = [
self.risk_free_rate * factor * self._scale_constraints
- ptf_min_acceptable_return * self._scale_constraints
<= v * self._scale_constraints
]
return risk, constraints
def _standard_deviation_risk(
self, return_distribution: ReturnDistribution, w: cp.Variable
) -> skt.RiskResult:
"""Expression and constraints of the standard deviation risk measure.
Parameters
----------
return_distribution : ReturnDistribution
Asset return distribution.
w : cvxpy Variable
CVXPY variable representing the asset weights.
Returns
-------
expression : tuple[cvxpy Expression, list[cvxpy Expression]]
CVXPY expression and constraints of the standard deviation risk measure.
"""
# nonneg=True instead of a separate constraint improves DCP analysis.
risk = cp.Variable(nonneg=True)
scale = self._scale_constraints
covariance_sqrt = return_distribution.covariance_sqrt
terms = [component.T @ w for component in covariance_sqrt.components]
if covariance_sqrt.diagonal is not None:
terms.append(cp.multiply(covariance_sqrt.diagonal, w))
constraints = [cp.SOC(risk * scale, cp.hstack(terms) * scale)]
return risk, constraints
def _variance_risk(
self, return_distribution: ReturnDistribution, w: cp.Variable
) -> skt.RiskResult:
"""Expression and constraints of the variance risk measure.
Parameters
----------
return_distribution : ReturnDistribution
Asset return distribution.
w : cvxpy Variable
CVXPY variable representing the asset weights.
Returns
-------
expression : tuple[cvxpy Expression, list[cvxpy Expression]]
CVXPY expression and constraints of the variance risk measure.
"""
risk, constraints = self._standard_deviation_risk(
return_distribution=return_distribution, w=w
)
risk = cp.square(risk)
return risk, constraints
def _worst_case_variance_risk(
self,
return_distribution: ReturnDistribution,
covariance_uncertainty_set: UncertaintySet | CompactCovarianceUncertaintySet,
w: cp.Variable,
factor: skt.Factor,
) -> skt.RiskResult:
r"""Expression and constraints of the worst-case variance.
A :class:`~skfolio.uncertainty_set.CompactCovarianceUncertaintySet`
adds the reduced quadratic penalty
.. math::
\kappa \min_z \lVert Cw - Qz \rVert_2^2
to the nominal variance. A generic
:class:`~skfolio.uncertainty_set.UncertaintySet` uses a lifted semidefinite
formulation.
Parameters
----------
return_distribution : ReturnDistribution
Asset return distribution.
covariance_uncertainty_set : UncertaintySet | CompactCovarianceUncertaintySet
Fitted :ref:`covariance uncertainty set <uncertainty_set_estimator>`.
w : cvxpy Variable
CVXPY variable representing the asset weights.
factor : cvxpy Variable | cvxpy Constant
Homogenization factor used by the generic lifted formulation for ratio
optimization.
Returns
-------
expression : tuple[cvxpy Expression, list[cvxpy Expression]]
CVXPY expression and constraints of the worst-case variance.
"""
if isinstance(covariance_uncertainty_set, CompactCovarianceUncertaintySet):
# The compact representation avoids lifted matrix variables.
risk, constraints = self._variance_risk(return_distribution, w)
residual = cp.multiply(covariance_uncertainty_set.metric_sqrt, w)
rank = covariance_uncertainty_set.basis.shape[1]
if rank > 0:
projection_coefficients = cp.Variable(rank)
residual -= covariance_uncertainty_set.basis @ projection_coefficients
risk += covariance_uncertainty_set.radius * cp.sum_squares(residual)
return risk, constraints
# Generic covariance uncertainty uses a lifted semidefinite formulation.
n_assets = return_distribution.returns.shape[1]
x = cp.Variable((n_assets, n_assets), symmetric=True)
y = cp.Variable((n_assets, n_assets), symmetric=True)
w_reshaped = cp.reshape(w, (n_assets, 1), order="F")
factor_reshaped = cp.reshape(factor, (1, 1), order="F")
z1 = cp.vstack([x, w_reshaped.T])
z2 = cp.vstack([w_reshaped, factor_reshaped])
# cvxpy annotates `p` as int or str but accepts any float, including inf.
risk = covariance_uncertainty_set.radius * cp.pnorm(
covariance_uncertainty_set.geometry.T
@ (cp.vec(x, order="F") + cp.vec(y, order="F")),
covariance_uncertainty_set.dual_norm, # ty: ignore[invalid-argument-type]
) + cp.trace(return_distribution.covariance @ (x + y))
constraints = [
cp.hstack([z1, z2]) * self._scale_constraints >> 0,
y * self._scale_constraints >> 0,
]
return risk, constraints
def _semi_variance_risk(
self,
return_distribution: ReturnDistribution,
w: cp.Variable,
min_acceptable_return: skt.Target | None = None,
) -> skt.RiskResult:
"""Expression and Constraints of the Semi Variance risk measure.
Parameters
----------
return_distribution : ReturnDistribution
asset returns distribution DataModel.
w : cvxpy Variable
The CVXPY Variable representing assets weights.
min_acceptable_return : float | ndarray of shape (n_assets,)
The minimum acceptable return used to distinguish "downside" and "upside"
returns for the computation of lower partial moments.
Returns
-------
expression : tuple[cvxpy Expression , list[cvxpy Expression]]
CVXPY expression and constraints of the semi variance risk measure.
"""
n_observations = return_distribution.returns.shape[0]
ptf_min_acceptable_return = self._cvx_min_acceptable_return(
return_distribution=return_distribution,
w=w,
min_acceptable_return=min_acceptable_return,
)
v = cp.Variable(n_observations, nonneg=True)
if return_distribution.sample_weight is None:
risk = cp.sum_squares(v) / (n_observations - 1)
else:
sample_weight = return_distribution.sample_weight
# Match the sample correction used by the numerical semi-variance.
correction = 1 - sample_weight @ sample_weight
_validate_positive_real(
float(correction), "semi-variance correction 1 - sum(sample_weight**2)"
)
risk = cp.sum_squares(cp.multiply(np.sqrt(sample_weight), v)) / correction
constraints = [
ptf_min_acceptable_return * self._scale_constraints
>= -v * self._scale_constraints
]
return risk, constraints
def _semi_deviation_risk(
self,
return_distribution: ReturnDistribution,
w: cp.Variable,
min_acceptable_return: skt.Target | None = None,
) -> skt.RiskResult:
"""Expression and Constraints of the Semi Standard Deviation risk measure.
Parameters
----------
return_distribution : ReturnDistribution
asset returns distribution DataModel.
w : cvxpy Variable
The CVXPY Variable representing assets weights.
min_acceptable_return : float | ndarray of shape (n_assets,)
The minimum acceptable return used to distinguish "downside" and "upside"
returns for the computation of lower partial moments.
Returns
-------
expression : tuple[cvxpy Expression , list[cvxpy Expression]]
CVXPY expression and constraints of the semi standard deviation risk measure.
"""
n_observations = return_distribution.returns.shape[0]
ptf_min_acceptable_return = self._cvx_min_acceptable_return(
return_distribution=return_distribution,
w=w,
min_acceptable_return=min_acceptable_return,
)
v = cp.Variable(n_observations, nonneg=True)
if return_distribution.sample_weight is None:
risk = cp.norm(v, 2) / np.sqrt(n_observations - 1)
else:
sample_weight = return_distribution.sample_weight
# Match the sample correction used by the numerical semi-deviation.
correction = 1 - sample_weight @ sample_weight
_validate_positive_real(
float(correction), "semi-deviation correction 1 - sum(sample_weight**2)"
)
risk = cp.norm(cp.multiply(np.sqrt(sample_weight), v), 2) / np.sqrt(
correction
)
constraints = [
ptf_min_acceptable_return * self._scale_constraints
>= -v * self._scale_constraints
]
return risk, constraints
def _fourth_central_moment_risk(
self, w: cp.Variable, factor: skt.Factor
) -> skt.RiskResult:
"""Fourth central moment risk, not supported in convex optimization."""
raise NotImplementedError
def _fourth_lower_partial_moment_risk(
self, w: cp.Variable, factor: skt.Factor
) -> skt.RiskResult:
"""Fourth lower partial moment risk, not supported in convex optimization."""
raise NotImplementedError
def _worst_realization_risk(
self,
return_distribution: ReturnDistribution,
w: cp.Variable,
factor: skt.Factor,
) -> skt.RiskResult:
"""Expression and Constraints of the Worst Realization risk measure.
Parameters
----------
return_distribution : ReturnDistribution
asset returns distribution DataModel.
w : cvxpy Variable
The CVXPY Variable representing assets weights.
factor : cvxpy Variable | cvxpy Constant
Additional variable used for the optimization of some objective function
like the ratio maximization.
Returns
-------
expression : tuple[cvxpy Expression , list[cvxpy Expression]]
CVXPY expression and constraints of the worst realization risk measure.
"""
ptf_returns = self._cvx_returns(return_distribution=return_distribution, w=w)
ptf_transaction_cost = self._cvx_transaction_cost(
return_distribution=return_distribution, w=w, factor=factor
)
ptf_management_fee = self._cvx_management_fee(
return_distribution=return_distribution, w=w
)
v = cp.Variable()
risk = v
constraints = [
-ptf_returns * self._scale_constraints
+ ptf_transaction_cost * self._scale_constraints
+ ptf_management_fee * self._scale_constraints
<= v * self._scale_constraints
]
return risk, constraints
def _cvar_risk(
self,
return_distribution: ReturnDistribution,
w: cp.Variable,
factor: skt.Factor,
) -> skt.RiskResult:
"""Expression and Constraints of the CVaR risk measure.
Parameters
----------
return_distribution : ReturnDistribution
asset returns distribution DataModel.
w : cvxpy Variable
The CVXPY Variable representing assets weights.
factor : cvxpy Variable | cvxpy Constant
Additional variable used for the optimization of some objective function
like the ratio maximization.
Returns
-------
expression : tuple[cvxpy Expression , list[cvxpy Expression]]
CVXPY expression and constraints of the CVaR risk measure.
"""
n_observations = return_distribution.returns.shape[0]
ptf_returns = self._cvx_returns(return_distribution=return_distribution, w=w)
ptf_transaction_cost = self._cvx_transaction_cost(
return_distribution=return_distribution, w=w, factor=factor
)
ptf_management_fee = self._cvx_management_fee(
return_distribution=return_distribution, w=w
)
if self.cvar_beta == 1:
return _cvx_worst_loss(
-ptf_returns + ptf_transaction_cost + ptf_management_fee,
sample_weight=return_distribution.sample_weight,
), []
alpha = cp.Variable()
v = cp.Variable(n_observations, nonneg=True)
if return_distribution.sample_weight is None:
risk = alpha + cp.sum(v) / (n_observations * (1 - self.cvar_beta))
else:
risk = alpha + cp.sum(cp.multiply(return_distribution.sample_weight, v)) / (
1 - self.cvar_beta
)
constraints = [
ptf_returns * self._scale_constraints
- ptf_transaction_cost * self._scale_constraints
- ptf_management_fee * self._scale_constraints
+ alpha * self._scale_constraints
+ v * self._scale_constraints
>= 0
]
return risk, constraints
def _evar_risk(
self,
return_distribution: ReturnDistribution,
w: cp.Variable,
factor: skt.Factor,
) -> skt.RiskResult:
"""Expression and Constraints of the EVaR risk measure.
Parameters
----------
return_distribution : ReturnDistribution
asset returns distribution DataModel.
w : cvxpy Variable
The CVXPY Variable representing assets weights.
factor : cvxpy Variable | cvxpy Constant
Additional variable used for the optimization of some objective function
like the ratio maximization.
Returns
-------
expression : tuple[cvxpy Expression , list[cvxpy Expression]]
CVXPY expression and constraints of the EVaR risk measure.
"""
ptf_returns = self._cvx_returns(return_distribution=return_distribution, w=w)
ptf_transaction_cost = self._cvx_transaction_cost(
return_distribution=return_distribution, w=w, factor=factor
)
ptf_management_fee = self._cvx_management_fee(
return_distribution=return_distribution, w=w
)
# We don't include the transaction_cost in the constraint otherwise the problem
# is not DCP
if not isinstance(ptf_transaction_cost, cp.Constant):
warnings.warn(
"The EVaR problem will be relaxed by removing the transaction costs"
" from the Cone constraint to keep the problem DCP. The solution may"
" not be accurate.",
stacklevel=2,
)
return _cvx_entropic_risk(
losses=-ptf_returns + ptf_management_fee,
beta=self.evar_beta,
sample_weight=return_distribution.sample_weight,
scale_constraints=self._scale_constraints,
)
def _max_drawdown_risk(
self,
return_distribution: ReturnDistribution,
w: cp.Variable,
factor: skt.Factor,
) -> skt.RiskResult:
"""Expression and Constraints of the Maximum Drawdown risk measure.
Parameters
----------
return_distribution : ReturnDistribution
asset returns distribution DataModel.
w : cvxpy Variable
The CVXPY Variable representing assets weights.
factor : cvxpy Variable | cvxpy Constant
Additional variable used for the optimization of some objective function
like the ratio maximization.
Returns
-------
expression : tuple[cvxpy Expression , list[cvxpy Expression]]
CVXPY expression and constraints of the maximum drawdown risk measure.
"""
v, constraints = self._cvx_drawdown(
return_distribution=return_distribution, w=w, factor=factor
)
u = cp.Variable()
risk = u
constraints += [u * self._scale_constraints >= v[1:] * self._scale_constraints]
return risk, constraints
def _average_drawdown_risk(
self,
return_distribution: ReturnDistribution,
w: cp.Variable,
factor: skt.Factor,
) -> skt.RiskResult:
"""Expression and Constraints of the Average Drawdown risk measure.
Parameters
----------
return_distribution : ReturnDistribution
asset returns distribution DataModel.
w : cvxpy Variable
The CVXPY Variable representing assets weights.
factor : cvxpy Variable | cvxpy Constant
Additional variable used for the optimization of some objective function
like the ratio maximization.
Returns
-------
expression : tuple[cvxpy Expression , list[cvxpy Expression]]
CVXPY expression and constraints of the average drawdown risk measure.
"""
n_observations = return_distribution.returns.shape[0]
v, constraints = self._cvx_drawdown(
return_distribution=return_distribution, w=w, factor=factor
)
if return_distribution.sample_weight is None:
risk = cp.sum(v[1:]) / n_observations
else:
risk = return_distribution.sample_weight @ v[1:]
return risk, constraints
def _cdar_risk(
self,
return_distribution: ReturnDistribution,
w: cp.Variable,
factor: skt.Factor,
) -> skt.RiskResult:
"""Expression and Constraints of the CDaR risk measure.
Parameters
----------
return_distribution : ReturnDistribution
asset returns distribution DataModel.
w : cvxpy Variable
The CVXPY Variable representing assets weights.
factor : cvxpy Variable | cvxpy Constant
Additional variable used for the optimization of some objective function
like the ratio maximization.
Returns
-------
expression : tuple[cvxpy Expression , list[cvxpy Expression]]
CVXPY expression and constraints of the CDaR risk measure.
"""
n_observations = return_distribution.returns.shape[0]
v, constraints = self._cvx_drawdown(
return_distribution=return_distribution, w=w, factor=factor
)
if self.cdar_beta == 1:
return _cvx_worst_loss(
v[1:], sample_weight=return_distribution.sample_weight
), constraints
alpha = cp.Variable()
z = cp.Variable(n_observations, nonneg=True)
if return_distribution.sample_weight is None:
tail_average = cp.sum(z) / n_observations
else:
tail_average = return_distribution.sample_weight @ z
risk = alpha + tail_average / (1 - self.cdar_beta)
constraints += [
z * self._scale_constraints
>= v[1:] * self._scale_constraints - alpha * self._scale_constraints
]
return risk, constraints
def _edar_risk(
self,
return_distribution: ReturnDistribution,
w: cp.Variable,
factor: skt.Factor,
) -> skt.RiskResult:
"""Expression and Constraints of the EDaR risk measure.
Parameters
----------
return_distribution : ReturnDistribution
asset returns distribution DataModel.
w : cvxpy Variable
The CVXPY Variable representing assets weights.
factor : cvxpy Variable | cvxpy Constant
Additional variable used for the optimization of some objective function
like the ratio maximization.
Returns
-------
expression : tuple[cvxpy Expression , list[cvxpy Expression]]
CVXPY expression and constraints of the EDaR risk measure.
"""
v, constraints = self._cvx_drawdown(
return_distribution=return_distribution, w=w, factor=factor
)
risk, risk_constraints = _cvx_entropic_risk(
losses=v[1:],
beta=self.edar_beta,
sample_weight=return_distribution.sample_weight,
scale_constraints=self._scale_constraints,
)
return risk, constraints + risk_constraints
def _ulcer_index_risk(
self,
return_distribution: ReturnDistribution,
w: cp.Variable,
factor: skt.Factor,
) -> skt.RiskResult:
"""Expression and Constraints of the Ulcer Index risk measure.
Parameters
----------
return_distribution : ReturnDistribution
asset returns distribution DataModel.
w : cvxpy Variable
The CVXPY Variable representing assets weights.
factor : cvxpy Variable | cvxpy Constant
Additional variable used for the optimization of some objective function
like the ratio maximization.
Returns
-------
expression : tuple[cvxpy Expression , list[cvxpy Expression]]
CVXPY expression and constraints of the Ulcer Index risk measure.
"""
v, constraints = self._cvx_drawdown(
return_distribution=return_distribution, w=w, factor=factor
)
n_observations = return_distribution.returns.shape[0]
if return_distribution.sample_weight is None:
risk = cp.norm(v[1:], 2) / np.sqrt(n_observations)
else:
risk = cp.norm(
cp.multiply(np.sqrt(return_distribution.sample_weight), v[1:]), 2
)
return risk, constraints
def _gini_mean_difference_risk(
self,
return_distribution: ReturnDistribution,
w: cp.Expression,
) -> skt.RiskResult:
r"""Build the Gini mean difference risk expression and its constraints.
A sparse sorting network represents an upper bound on portfolio GMD with
:math:`O(T \log^2 T)` variables and constraints for :math:`T` observations.
Minimizing the risk expression makes this bound exact.
Parameters
----------
return_distribution : ReturnDistribution
Fitted distribution containing asset returns with at least two
observations.
w : cvxpy Expression of shape (n_assets,)
Portfolio weights in optimization units. For ratio objectives, these
are the weights before division by the homogenization factor. To model
risk relative to a target, subtract target weights in those same units.
Returns
-------
risk : cvxpy Expression
Scalar upper bound on the GMD of `return_distribution.returns @ w`.
constraints : list[cvxpy Constraint]
Linear constraints linking portfolio returns to the sorting network.
Raises
------
ValueError
If the return distribution contains fewer than two observations.
Notes
-----
Asset returns are centered and scaled to improve numerical conditioning.
The risk expression restores the original return units. When GMD is only
upper bounded, the expression can exceed the empirical GMD of the final
weights. Use `_gmd_reported_risk` to report that empirical value.
"""
returns = return_distribution.returns
n_observations = returns.shape[0]
equality, inequality, risk_coefficients = _gmd_sorting_network(n_observations)
z = cp.Variable(equality.shape[1], name="gmd_network")
# GMD is translation invariant. Centering returns and omitting costs and
# fees remove observation-independent shifts without changing the risk.
centered_returns = returns - returns.mean(axis=0)
# Keep network coefficients near one, including when ratio optimization
# produces small homogeneous weights. Restore return units in the risk.
return_scale = np.max(np.abs(centered_returns))
if return_scale == 0:
return_scale = 1.0
constraints = [
z[:n_observations] * self._scale_constraints
== (centered_returns / return_scale) @ w * self._scale_constraints,
equality @ z * self._scale_constraints == 0,
inequality @ z * self._scale_constraints <= 0,
]
return return_scale * (risk_coefficients @ z), constraints
[docs]
@abstractmethod
def fit(
self, X: ArrayLike, y: ArrayLike | None = None, **fit_params: Any
) -> ConvexOptimization:
"""Fit the Convex Optimization estimator.
Parameters
----------
X : array-like of shape (n_observations, n_assets)
Price returns of the assets.
y : array-like of shape (n_observations, n_targets), optional
Price returns of factors or a target benchmark.
The default is `None`.
**fit_params : dict
Parameters to pass to the underlying estimators.
Only available if `enable_metadata_routing=True`, which can be
set by using `sklearn.set_config(enable_metadata_routing=True)`.
See :ref:`Metadata Routing User Guide <metadata_routing>` for
more details.
Returns
-------
self : ConvexOptimization
Fitted estimator.
"""
...
def _mip_weight_constraints_no_short_threshold(
n_assets: int,
w: cp.Variable,
factor: skt.Factor,
scale_constraints: cp.Constant,
cardinality: int | None,
group_cardinalities: dict[str, int] | None,
max_weights: float | FloatArray,
groups: AnyArray | None,
min_weights: float | FloatArray,
threshold_long: float | FloatArray | None,
) -> list[cpc.Constraint]:
"""
Create a list of MIP constraints for cardinality and threshold conditions
when no short threshold is present. This only requires the creation of a single
boolean variable array.
"""
constraints = []
is_short = np.any(min_weights < 0)
is_invested_bool = cp.Variable(n_assets, boolean=True)
if cardinality is not None and cardinality < n_assets:
constraints.append(cp.sum(is_invested_bool) <= cardinality)
if group_cardinalities is not None:
a_card, b_card = group_cardinalities_to_matrix(
groups=groups,
group_cardinalities=group_cardinalities,
raise_if_group_missing=False,
)
constraints.append(a_card @ is_invested_bool - b_card <= 0)
if isinstance(factor, cp.Variable):
is_invested_factor = cp.Variable(n_assets, nonneg=True)
# We want (w <= cp.multiply(is_invested_short_bool, max_weights) * factor
# but this is not DCP. So we introduce another variable and set
# constraint to ensure its value is equal to is_invested_short_bool * factor
M = 1e3
# Big M method to activate or deactivate constraints
# In the ratio homogenization procedure, the factor has been calibrated
# to be around 0.1-10. By using M=1e3, we ensure that M is large enough while
# not too large for improved MIP convergence.
constraints += [
is_invested_factor <= factor,
is_invested_factor <= M * is_invested_bool,
is_invested_factor >= factor - M * (1 - is_invested_bool),
]
is_invested = is_invested_factor
else:
is_invested = is_invested_bool
if threshold_long is not None:
constraints.append(
w * scale_constraints
>= cp.multiply(is_invested, threshold_long) * scale_constraints
)
constraints.append(
w * scale_constraints
<= cp.multiply(is_invested, max_weights) * scale_constraints
)
if is_short:
constraints.append(
w * scale_constraints
>= cp.multiply(is_invested, min_weights) * scale_constraints
)
return constraints
def _mip_weight_constraints_threshold_short(
n_assets: int,
w: cp.Variable,
factor: skt.Factor,
scale_constraints: cp.Constant,
max_weights: float | FloatArray,
min_weights: float | FloatArray,
threshold_long: float | FloatArray,
threshold_short: float | FloatArray,
cardinality: int | None,
group_cardinalities: dict[str, int] | None,
groups: AnyArray | None,
) -> list[cpc.Constraint]:
"""
Create a list of MIP constraints for cardinality and threshold constraints
when a short threshold is allowed. This requires the creation of two boolean
variable arrays, one for long positions and one for short positions.
"""
constraints = []
is_invested_short_bool = cp.Variable(n_assets, boolean=True)
is_invested_long_bool = cp.Variable(n_assets, boolean=True)
is_invested_bool = is_invested_short_bool + is_invested_long_bool
if cardinality is not None and cardinality < n_assets:
constraints.append(cp.sum(is_invested_bool) <= cardinality)
if group_cardinalities is not None:
a_card, b_card = group_cardinalities_to_matrix(
groups=groups,
group_cardinalities=group_cardinalities,
raise_if_group_missing=False,
)
constraints.append(a_card @ is_invested_bool - b_card <= 0)
M = 1e3
# Big M method to activate or deactivate constraints
# In the ratio homogenization procedure, the factor has been calibrated
# to be around 0.1-10. By using M=1e3, we ensure that M is large enough while
# not too large for improved MIP convergence.
if isinstance(factor, cp.Variable):
is_invested_short_factor = cp.Variable(n_assets, nonneg=True)
is_invested_long_factor = cp.Variable(n_assets, nonneg=True)
# We want (w <= cp.multiply(is_invested_short_bool, max_weights) * factor
# but this is not DCP. So we introduce another variable and set
# constraint to ensure its value is equal to is_invested_short_bool * factor
constraints += [
is_invested_short_factor <= factor,
is_invested_long_factor <= factor,
is_invested_short_factor <= M * is_invested_short_bool,
is_invested_long_factor <= M * is_invested_long_bool,
is_invested_short_factor >= factor - M * (1 - is_invested_short_bool),
is_invested_long_factor >= factor - M * (1 - is_invested_long_bool),
]
is_invested_short = is_invested_short_factor
is_invested_long = is_invested_long_factor
else:
is_invested_short = is_invested_short_bool
is_invested_long = is_invested_long_bool
constraints += [
is_invested_bool <= 1.0,
w * scale_constraints
<= cp.multiply(is_invested_long, max_weights) * scale_constraints,
w * scale_constraints
>= cp.multiply(is_invested_short, min_weights) * scale_constraints,
# Apply threshold_long if is_invested_long == 1,
# unrestricted if is_invested_long == 0
w * scale_constraints
>= cp.multiply(is_invested_long, threshold_long) * scale_constraints
- M * (1 - is_invested_long_bool) * scale_constraints,
# # Apply threshold_short if is_invested_short == 1,
# # unrestricted if is_invested_short == 0
w * scale_constraints
<= cp.multiply(is_invested_short, threshold_short) * scale_constraints
+ M * (1 - is_invested_short_bool) * scale_constraints,
]
return constraints
def _solve(
w: cp.Variable,
factor: skt.Factor,
expressions: dict[str, cp.Expression],
problem: cp.Problem,
solver: str,
solver_params: dict,
) -> tuple[FloatArray, dict[str, float]]:
"""Solve `problem` and return normalized weights and reported values.
Parameters
----------
w : cvxpy Variable of shape (n_assets,)
Portfolio weights in optimization units.
factor : cvxpy Variable | cvxpy Constant
Scalar normalization factor. Its solved value must be finite and greater
than 1e-12.
expressions : dict[str, cvxpy Expression]
Named scalar expressions in reporting units, evaluated after solving.
problem : cvxpy Problem
Optimization problem containing `w` and its supporting constraints.
solver : str
Name of the CVXPY solver to use.
solver_params : dict
Additional keyword arguments passed to `problem.solve`.
Returns
-------
weights : ndarray of shape (n_assets,)
Final portfolio weights, computed as `w.value / factor.value`.
problem_values : dict[str, float]
Values of the supplied reporting expressions.
Raises
------
ConvexOptimizationError
If solving fails, weights are unavailable, the factor is invalid, or the
normalized weights or reported values are non-finite.
Warns
-----
UserWarning
If the solver returns a usable solution with a status other than `optimal`.
"""
try:
# We suppress cvxpy warning as it is redundant with our warning
with warnings.catch_warnings():
warnings.simplefilter("ignore")
problem.solve(solver=solver, **solver_params)
if w.value is None:
raise cp.SolverError("No solution found")
except (cp.SolverError, sla.ArpackNoConvergence) as solver_error:
params_string = " ".join([f"{p.value:0g}" for p in problem.parameters()])
if len(params_string) != 0:
params_string = f" with parameters {params_string}"
error = (
f"Solver '{solver}' failed{params_string}. Try another"
" solver, or solve with solver_params=dict(verbose=True) for more"
" information"
)
raise ConvexOptimizationError(error) from solver_error
factor_value = factor.value
if (
factor_value is None
or not np.isfinite(factor_value)
or factor_value <= _MIN_HOMOGENIZATION_FACTOR
):
raise ConvexOptimizationError(
"The optimization returned an invalid homogeneous normalization factor"
)
factor_value = float(factor_value)
with np.errstate(over="ignore", invalid="ignore"):
weights = w.value / factor_value
if not np.isfinite(weights).all():
raise ConvexOptimizationError("Allocation produced non-finite weights.")
with np.errstate(over="ignore", invalid="ignore"):
problem_values = {
name: expression.value for name, expression in expressions.items()
}
for name, value in problem_values.items():
if not np.isfinite(value).all():
raise ConvexOptimizationError(
f"Optimization produced a non-finite '{name}' value."
)
weights = np.array(weights, dtype=float)
if not problem.status == cp.OPTIMAL:
warnings.warn(
"Solution may be inaccurate. Try changing the solver params or the"
" scale. For more details, set `solver_params=dict(verbose=True)`",
stacklevel=2,
)
return weights, problem_values # ty: ignore[invalid-return-type]
def _cvx_worst_loss(
losses: cp.Expression, *, sample_weight: FloatArray | None
) -> cp.Expression:
"""Return the largest loss with positive probability.
Parameters
----------
losses : cvxpy Expression of shape (n_observations,)
Losses or drawdown magnitudes for all observations.
sample_weight : ndarray of shape (n_observations,) or None
Validated scenario probabilities. None gives every observation equal weight.
Returns
-------
risk : cvxpy Expression
Maximum over observations with positive probability.
"""
if sample_weight is not None:
positive = sample_weight > 0
if not positive.all():
losses = losses[positive]
return cp.max(losses)
def _cvx_entropic_risk(
losses: cp.Expression,
*,
beta: float,
sample_weight: FloatArray | None,
scale_constraints: cp.Constant,
) -> skt.RiskResult:
"""Create the CVXPY expression and constraints for EVaR.
Pass portfolio losses to model EVaR, or nonnegative drawdown magnitudes
to model EDaR.
Parameters
----------
losses : cvxpy Expression of shape (n_observations,)
Loss at each observation. For EDaR, pass nonnegative drawdown
magnitudes computed from the full return path, including dates
whose sample weight is zero.
beta : float
Confidence level between 0 and 1, inclusive. At zero, the risk equals
the mean loss. At one, it equals the largest loss with positive weight.
sample_weight : ndarray of shape (n_observations,) or None
Finite, nonnegative probabilities from the return distribution,
summing to one. If None, use equal probabilities.
scale_constraints : cvxpy Constant
Positive factor used to scale the constraints for the solver.
Returns
-------
expression : tuple[cvxpy Expression, list[cvxpy Expression]]
Risk expression and its constraints.
Notes
-----
Nonuniform probabilities enter the exponential cone as logarithms.
This avoids using very small probabilities as constraint coefficients,
which solvers can treat as zero. Only exactly zero probabilities are
excluded from the losses.
"""
if beta == 1:
return _cvx_worst_loss(losses, sample_weight=sample_weight), []
if sample_weight is not None:
positive = sample_weight > 0
if not positive.all():
losses = losses[positive]
sample_weight = sample_weight[positive]
if np.all(sample_weight == sample_weight[0]):
sample_weight = None
n_observations = losses.shape[0]
if beta == 0:
if sample_weight is None:
risk = cp.sum(losses) / n_observations
else:
risk = sample_weight @ losses
return risk, []
location = cp.Variable()
temperature = cp.Variable(nonneg=True)
exponential = cp.Variable(n_observations)
cone_losses = losses - location
if sample_weight is None:
# Keep the existing equal-weight formulation and its solver scaling.
log_tail = np.log(n_observations * (1 - beta))
else:
cone_losses += temperature * np.log(sample_weight)
log_tail = np.log1p(-beta)
risk = location - temperature * log_tail
constraints = [
cp.sum(exponential) * scale_constraints <= temperature * scale_constraints,
cp.constraints.ExpCone(
cone_losses * scale_constraints,
np.ones(n_observations) * temperature * scale_constraints,
exponential * scale_constraints,
),
]
return risk, constraints
def _gmd_reported_risk(returns: FloatArray, w: cp.Expression) -> cp.Expression:
"""Build an empirical GMD expression for evaluation after solving.
The expression computes GMD directly from portfolio returns when its value
is read after solving.
Parameters
----------
returns : ndarray of shape (n_observations, n_assets)
Asset returns with at least two observations.
w : cvxpy Expression of shape (n_assets,)
Portfolio weights after division by the homogenization factor. Subtract
target weights to report the GMD relative to a target portfolio.
Returns
-------
risk : cvxpy Expression
Scalar empirical GMD expression for reporting after solving.
"""
centered_returns = returns - returns.mean(axis=0)
return cp.dotsort(centered_returns @ w, owa_gmd_weights(len(returns)))
def _odd_even_merge_sort_comparators(n: int) -> list[tuple[int, int]]:
"""Generate the compare-and-swap operations for a Batcher sorting network.
Parameters
----------
n : int
Number of input values to sort. Counts that are not powers of two are
supported by omitting comparisons involving positions beyond the input.
Returns
-------
comparators : list[tuple[int, int]]
Pairs of zero-based positions in execution order. For each pair `(i, j)`,
place the smaller value at position `i` and the larger value at position
`j`. Applying all pairs sorts the input in ascending order.
"""
comparators = []
p = 1
while p < n:
k = p
while k >= 1:
for j in range(k % p, n - k, 2 * k):
for i in range(min(k, n - j - k)):
if (i + j) // (2 * p) == (i + j + k) // (2 * p):
comparators.append((i + j, i + j + k))
k //= 2
p *= 2
return comparators
@functools.lru_cache(maxsize=8)
def _gmd_sorting_network(
n_observations: int,
) -> tuple[sp.csr_matrix, sp.csr_matrix, FloatArray]:
r"""Build cached sparse matrices for the GMD sorting network.
For fixed input observations, minimizing the network's linear risk expression
under these constraints gives their empirical GMD.
Parameters
----------
n_observations : int
Number of input observations. Must be at least two.
Returns
-------
equality : scipy.sparse.csr_matrix of shape (n_comparators, n_variables)
Coefficients of the comparator equalities `equality @ z == 0`.
inequality : scipy.sparse.csr_matrix of shape (2 * n_comparators, n_variables)
Coefficients of the comparator inequalities `inequality @ z <= 0`.
risk_coefficients : ndarray of shape (n_variables,)
Ordered weighted average coefficients on the final network outputs. The
risk expression is `risk_coefficients @ z`.
Raises
------
ValueError
If `n_observations` is less than two.
Notes
-----
The variable vector `z` starts with `n_observations` inputs, followed by two
output variables for each comparator. Thus, `n_variables` equals
`n_observations + 2 * n_comparators`. The caller binds the inputs to portfolio
returns. The returned matrices and coefficients are cached and read-only.
Each comparator with inputs `u` and `v` has outputs `lo` and `hi` constrained
by `lo + hi = u + v`, `hi >= u` and `hi >= v`. These constraints relax sorting,
so feasible auxiliary values may be unsorted. Minimizing the weighted final
outputs yields the exact GMD.
The formulation uses :math:`O(T \log^2 T)` variables and constraints for
:math:`T` observations. It is the linear programming dual of Goemans' (2015)
permutahedron formulation:
https://doi.org/10.1007/s10107-014-0757-1
"""
if n_observations < 2:
raise ValueError("GMD requires at least two observations")
comparators = _odd_even_merge_sort_comparators(n_observations)
n_comparators = len(comparators)
n_variables = n_observations + 2 * n_comparators
source = np.arange(n_observations)
lo = n_observations + 2 * np.arange(n_comparators)
hi = lo + 1
u = np.empty(n_comparators, dtype=int)
v = np.empty(n_comparators, dtype=int)
for k, (i, j) in enumerate(comparators):
u[k], v[k] = source[i], source[j]
source[i], source[j] = lo[k], hi[k]
rows = np.repeat(np.arange(n_comparators), 4)
cols = np.column_stack([lo, hi, u, v]).ravel()
vals = np.tile([1.0, 1.0, -1.0, -1.0], n_comparators)
equality = sp.csr_matrix((vals, (rows, cols)), shape=(n_comparators, n_variables))
rows = np.repeat(np.arange(2 * n_comparators), 2)
cols = np.column_stack([u, hi, v, hi]).ravel()
vals = np.tile([1.0, -1.0], 2 * n_comparators)
inequality = sp.csr_matrix(
(vals, (rows, cols)), shape=(2 * n_comparators, n_variables)
)
risk_coefficients = np.zeros(n_variables)
risk_coefficients[source] = owa_gmd_weights(n_observations)
# Cached arrays are shared by independent fits. Keep them immutable.
for matrix in (equality, inequality):
for values in (matrix.data, matrix.indices, matrix.indptr):
values.flags.writeable = False
risk_coefficients.flags.writeable = False
return equality, inequality, risk_coefficients
