quanttoolbox.portfolio

portfolio.risk_budgeting

Python alternatives

Keep — the single largest and most-tested piece of custom work in the whole port (75+ original MATLAB files consolidated). riskparityportfolio covers basic risk parity but not the box-constrained/general-linear-constrained/VaR-ES/target-matching breadth here. Genuinely the strongest "keep" case in the whole library.

quanttoolbox.portfolio.risk_budgeting

Risk budgeting / risk parity portfolio construction.

Ported from QuantToolBox/rpb/{compute_risk_contribution,compute_rc_sd, compute_rc_vol,compute_rb_sd,compute_rb_sd_ccd,compute_rb_sd_newton, lagrange_rb_sd,compute_erc_portfolio}.m and QuantToolBox/crb/compute_rb_sd_bc_admm1*.m (+ ~40 near-duplicate ADMM variants), and QuantToolBox/mloapa/compute_ERC_{ADMM,CCD}.m.

Consolidation notes -- this is the single largest simplification in the whole port:

The original has ~75 files across rpb/ and crb/ implementing risk budgeting under different combinations of {solver algorithm} x {constraint type}. The solver algorithms (CCD, Newton, ADMM+Newton, ADMM+CCD, ADMM+QP, ADMM+fmincon, bisection-on-lambda) are all just different numerical routes to the same mathematical solution -- they don't change what problem is being solved, only how. This module ports:

• risk_contribution -- the risk-decomposition arithmetic shared by every variant (consolidates compute_risk_contribution/compute_rc_sd/ compute_rc_vol).
• solve_unconstrained -- CCD (default, Roncalli's cyclical coordinate descent -- the standard reference algorithm) or Newton, for the classic budget-constrained-only (sum(x)=1) risk budgeting problem.
• solve_box_constrained -- ADMM (Newton-based x-update, closed-form box-projection z-update, bisection on the budget-constraint Lagrange multiplier lambda) for box-constrained [x_minus, x_plus] risk budgeting. This one function replaces the entire ~40-file family of crb/compute_rb_sd_bc_admm{1,2,3,4}(_lambda).m and crb/compute_rb_sd_bc_ccd(_lambda).m variants, which differ only in which inner solver (Newton/CCD/QP/fmincon) handles the ADMM x-update or whether lambda is searched via bisection vs. passed directly -- all converge to the same box-constrained solution.
• solve_constrained -- the same ADMM structure, generalized to any combination of extra linear equality/inequality constraints and box bounds (the z-update projects onto their intersection via optim.proximal.proximal_linear_constraints's Dykstra algorithm instead of a closed-form box clip). Consolidates the crb/compute_rb_sd_constrained_admm*.m family. solve_box_constrained is now a thin wrapper around this for the box-only case.
• solve_unconstrained_var/solve_unconstrained_es -- Value-at-Risk / Expected-Shortfall risk budgeting under a Gaussian assumption, by mapping the confidence level alpha to the equivalent standard-deviation multiplier c and reusing solve_unconstrained directly (consolidates compute_rb_var.m/compute_rb_es.m/ compute_rc_var.m/compute_rc_es.m -- each just a 5-line wrapper in the original, computing a different scalar multiplier).
• erc_portfolio -- Equal Risk Contribution, the b=1/n special case (consolidates compute_erc_portfolio and mloapa/compute_ERC_{ADMM,CCD}.m, which are just the unconstrained solver called with equal budgets).
• risk_budgeting_target -- bisects the risk-aversion parameter c to hit a target active return or active volatility (relative to an optional benchmark), covering the "mu-problem"/"sigma-problem" modes of compute_risk_parity_portfolio.m/compute_risk_parity_portfolio_bounds.m. Works with any of the three solvers above via a solver= argument ("unconstrained"/"box"/"constrained"), consolidating both original functions (and their four small _return/_volatility objective-function helper files) into one ~80-line implementation built on quanttoolbox.optim.bisection -- versus the originals' ~550 lines combined.

MATLAB's global RB_CCD_*/RB_Newton_*/RB_ADMM_* blocks are replaced by quanttoolbox.config.{CCDConfig,NewtonConfig,ADMMConfig}.

erc_portfolio(cov_matrix, x0=None, method='ccd')

Equal Risk Contribution portfolio: risk budgeting with all budgets equal (b = 1/n), pure volatility risk measure.

Original: rpb/compute_erc_portfolio.m (+ mloapa/compute_ERC_{ADMM,CCD}.m, equivalent alternative solvers for the same problem)

Source code in src/quanttoolbox/portfolio/risk_budgeting.py

def erc_portfolio(
    cov_matrix: np.ndarray, x0: np.ndarray | None = None, method: str = "ccd"
) -> RiskBudgetingResult:
    """Equal Risk Contribution portfolio: risk budgeting with all budgets
    equal (b = 1/n), pure volatility risk measure.

    Original: rpb/compute_erc_portfolio.m (+ mloapa/compute_ERC_{ADMM,CCD}.m,
    equivalent alternative solvers for the same problem)
    """
    n = np.asarray(cov_matrix).shape[0]
    return solve_unconstrained(
        cov_matrix, b=np.full(n, 1.0 / n), mu=0.0, c=1.0, x0=x0, method=method
    )

risk_budgeting_frontier(cov_matrix, b, mu, c_values, x0=None)

Evaluate the risk budgeting solution at each risk-aversion value in c_values (the "gamma-problem" mode of compute_risk_parity_portfolio.m).

See module docstring for target-matching modes not ported here.

Original: rpb/compute_risk_parity_portfolio.m (gamma-problem branch)

Source code in src/quanttoolbox/portfolio/risk_budgeting.py

def risk_budgeting_frontier(
    cov_matrix: np.ndarray,
    b: np.ndarray | None,
    mu: np.ndarray,
    c_values: np.ndarray,
    x0: np.ndarray | None = None,
) -> list[RiskBudgetingResult]:
    """Evaluate the risk budgeting solution at each risk-aversion value in
    c_values (the "gamma-problem" mode of compute_risk_parity_portfolio.m).

    See module docstring for target-matching modes not ported here.

    Original: rpb/compute_risk_parity_portfolio.m (gamma-problem branch)
    """
    cov_matrix = np.asarray(cov_matrix, dtype=float)
    return [
        solve_unconstrained(cov_matrix, b=b, mu=mu, c=float(c), x0=x0, method="ccd")
        for c in np.atleast_1d(c_values)
    ]

risk_budgeting_target(cov_matrix, mu, target, target_type='return', b=None, x_benchmark=None, x0=None, c_min=1.0, c_max=100.0, solver='unconstrained', bisection_config=None, **solver_kwargs)

Find the risk budgeting portfolio (relative to an optional benchmark x_benchmark) whose active return or active volatility matches target, by bisecting the risk-aversion parameter c between c_min and c_max.

target_type="return" (default): bisect to hit a target active return (the "mu-problem" mode). target_type="volatility": bisect to hit a target active volatility (the "sigma-problem" mode).

solver="unconstrained" (default), "box", or "constrained" selects which underlying risk-budgeting solver to use at each trial c; solver_kwargs are forwarded to it (e.g. x_minus/x_plus for "box", a_eq/b_eq/c_ineq/d_ineq/x_minus/x_plus for "constrained").

Returns a result with weights all-NaN and converged=False if target lies outside the achievable range [c_min, c_max] maps to.

Original: rpb/{compute_risk_parity_portfolio, compute_risk_parity_portfolio_bounds}.m ("mu-problem"/"sigma-problem" branches, consolidated with their _return/_volatility objective helper files -- see module docstring)

Source code in src/quanttoolbox/portfolio/risk_budgeting.py

def risk_budgeting_target(
    cov_matrix: np.ndarray,
    mu: np.ndarray,
    target: float,
    target_type: str = "return",
    b: np.ndarray | None = None,
    x_benchmark: np.ndarray | None = None,
    x0: np.ndarray | None = None,
    c_min: float = 1.0,
    c_max: float = 100.0,
    solver: str = "unconstrained",
    bisection_config: BisectionConfig | None = None,
    **solver_kwargs,
) -> RiskBudgetingTargetResult:
    """Find the risk budgeting portfolio (relative to an optional benchmark
    x_benchmark) whose active return or active volatility matches `target`,
    by bisecting the risk-aversion parameter c between c_min and c_max.

    target_type="return" (default): bisect to hit a target active return
    (the "mu-problem" mode). target_type="volatility": bisect to hit a
    target active volatility (the "sigma-problem" mode).

    solver="unconstrained" (default), "box", or "constrained" selects
    which underlying risk-budgeting solver to use at each trial c;
    solver_kwargs are forwarded to it (e.g. x_minus/x_plus for "box",
    a_eq/b_eq/c_ineq/d_ineq/x_minus/x_plus for "constrained").

    Returns a result with weights all-NaN and converged=False if `target`
    lies outside the achievable range [c_min, c_max] maps to.

    Original: rpb/{compute_risk_parity_portfolio,
    compute_risk_parity_portfolio_bounds}.m ("mu-problem"/"sigma-problem"
    branches, consolidated with their _return/_volatility objective
    helper files -- see module docstring)
    """
    cov_matrix = np.asarray(cov_matrix, dtype=float)
    n = cov_matrix.shape[0]
    mu = np.asarray(mu, dtype=float).flatten()
    x_benchmark_arr = (
        np.zeros(n) if x_benchmark is None else np.asarray(x_benchmark, dtype=float).flatten()
    )
    bisect_cfg = bisection_config or BisectionConfig()

    _, mu_min, sigma_min, _ = _solve_rb_at_c(
        cov_matrix, b, mu, c_min, x0, x_benchmark_arr, solver, solver_kwargs
    )
    _, mu_max, sigma_max, _ = _solve_rb_at_c(
        cov_matrix, b, mu, c_max, x0, x_benchmark_arr, solver, solver_kwargs
    )

    if target_type == "volatility":
        lo, hi = min(sigma_min, sigma_max), max(sigma_min, sigma_max)

        def objective(c: np.ndarray) -> np.ndarray:
            _, _, sigma_x, _ = _solve_rb_at_c(
                cov_matrix, b, mu, float(c), x0, x_benchmark_arr, solver, solver_kwargs
            )
            return np.array(sigma_x - target)

    elif target_type == "return":
        lo, hi = min(mu_min, mu_max), max(mu_min, mu_max)

        def objective(c: np.ndarray) -> np.ndarray:
            _, mu_x, _, _ = _solve_rb_at_c(
                cov_matrix, b, mu, float(c), x0, x_benchmark_arr, solver, solver_kwargs
            )
            return np.array(mu_x - target)

    else:
        raise ValueError(f"risk_budgeting_target: unknown target_type '{target_type}'")

    if target < lo or target > hi:
        nan_arr = np.full(n, np.nan)
        return RiskBudgetingTargetResult(
            weights=nan_arr,
            active_return=np.nan,
            active_volatility=np.nan,
            c=np.nan,
            converged=False,
        )

    c_star = bisection(objective, c_min, c_max, bisect_cfg)
    if np.isnan(c_star):
        nan_arr = np.full(n, np.nan)
        return RiskBudgetingTargetResult(
            weights=nan_arr,
            active_return=np.nan,
            active_volatility=np.nan,
            c=np.nan,
            converged=False,
        )

    x, mu_x, sigma_x, converged = _solve_rb_at_c(
        cov_matrix, b, mu, float(c_star), x0, x_benchmark_arr, solver, solver_kwargs
    )
    return RiskBudgetingTargetResult(
        weights=x,
        active_return=mu_x,
        active_volatility=sigma_x,
        c=float(c_star),
        converged=converged,
    )

risk_budgeting_target_frontier(cov_matrix, mu, targets, target_type='return', **kwargs)

Evaluate risk_budgeting_target at each value in targets.

Original: rpb/{compute_risk_parity_portfolio, compute_risk_parity_portfolio_bounds}.m (looping over multiple targets)

Source code in src/quanttoolbox/portfolio/risk_budgeting.py

def risk_budgeting_target_frontier(
    cov_matrix: np.ndarray,
    mu: np.ndarray,
    targets: np.ndarray,
    target_type: str = "return",
    **kwargs,
) -> list[RiskBudgetingTargetResult]:
    """Evaluate ``risk_budgeting_target`` at each value in `targets`.

    Original: rpb/{compute_risk_parity_portfolio,
    compute_risk_parity_portfolio_bounds}.m (looping over multiple targets)
    """
    return [
        risk_budgeting_target(cov_matrix, mu, float(t), target_type=target_type, **kwargs)
        for t in np.atleast_1d(targets)
    ]

risk_contribution(x, cov_matrix, mu=0.0, c=1.0)

Decompose portfolio risk = -x'mu + csqrt(x'Covx) into each asset's marginal and total risk contribution.

With mu=0, c=1 this is the pure volatility risk measure (consolidates compute_risk_contribution.m / compute_rc_vol.m). With mu != 0 it's the standard-deviation-based measure net of expected return (compute_rc_sd.m).

Original: rpb/{compute_risk_contribution,compute_rc_sd,compute_rc_vol}.m

Source code in src/quanttoolbox/portfolio/risk_budgeting.py

def risk_contribution(
    x: np.ndarray, cov_matrix: np.ndarray, mu: np.ndarray | float = 0.0, c: float = 1.0
) -> RiskContribution:
    """Decompose portfolio risk = -x'mu + c*sqrt(x'Cov*x) into each asset's
    marginal and total risk contribution.

    With mu=0, c=1 this is the pure volatility risk measure (consolidates
    compute_risk_contribution.m / compute_rc_vol.m). With mu != 0 it's the
    standard-deviation-based measure net of expected return
    (compute_rc_sd.m).

    Original: rpb/{compute_risk_contribution,compute_rc_sd,compute_rc_vol}.m
    """
    x = np.asarray(x, dtype=float).flatten()
    cov_matrix = np.asarray(cov_matrix, dtype=float)
    n = x.shape[0]
    mu = np.zeros(n) if np.isscalar(mu) and mu == 0 else np.asarray(mu, dtype=float).flatten()

    sigma = np.sqrt(x @ cov_matrix @ x)
    risk = -x @ mu + c * sigma
    mr = -mu + c * (cov_matrix @ x) / sigma
    rc = x * mr
    prc = rc / np.sum(rc)

    return RiskContribution(
        risk=risk, marginal_risk=mr, risk_contribution=rc, pct_risk_contribution=prc
    )

risk_contribution_es(x, cov_matrix, mu, alpha)

Risk contribution decomposition under the (Gaussian) ES risk measure.

Original: rpb/compute_rc_es.m

Source code in src/quanttoolbox/portfolio/risk_budgeting.py

def risk_contribution_es(
    x: np.ndarray, cov_matrix: np.ndarray, mu: np.ndarray | float, alpha: float
) -> RiskContribution:
    """Risk contribution decomposition under the (Gaussian) ES risk measure.

    Original: rpb/compute_rc_es.m
    """
    return risk_contribution(x, cov_matrix, mu, c=_es_multiplier(alpha))

risk_contribution_var(x, cov_matrix, mu, alpha)

Risk contribution decomposition under the (Gaussian) VaR risk measure.

Original: rpb/compute_rc_var.m

Source code in src/quanttoolbox/portfolio/risk_budgeting.py

def risk_contribution_var(
    x: np.ndarray, cov_matrix: np.ndarray, mu: np.ndarray | float, alpha: float
) -> RiskContribution:
    """Risk contribution decomposition under the (Gaussian) VaR risk measure.

    Original: rpb/compute_rc_var.m
    """
    return risk_contribution(x, cov_matrix, mu, c=_var_multiplier(alpha))

solve_box_constrained(cov_matrix, x_minus, x_plus, b=None, mu=0.0, c=1.0, x0=None, lambda_bracket=None, admm_config=None, newton_config=None)

Solve risk budgeting subject to box constraints x_minus <= x <= x_plus (in addition to the sum(x)=1 budget constraint). Thin wrapper around solve_constrained for the box-only case (uses the cheaper closed-form box projection directly, no Dykstra iteration needed).

Original: crb/compute_rb_sd_bc_admm1*.m (+ ~40 equivalent solver variants -- see module docstring)

Source code in src/quanttoolbox/portfolio/risk_budgeting.py

def solve_box_constrained(
    cov_matrix: np.ndarray,
    x_minus: np.ndarray | float,
    x_plus: np.ndarray | float,
    b: np.ndarray | None = None,
    mu: np.ndarray | float = 0.0,
    c: float = 1.0,
    x0: np.ndarray | None = None,
    lambda_bracket: tuple[float, float] | None = None,
    admm_config: ADMMConfig | None = None,
    newton_config: NewtonConfig | None = None,
) -> RiskBudgetingResult:
    """Solve risk budgeting subject to box constraints x_minus <= x <= x_plus
    (in addition to the sum(x)=1 budget constraint). Thin wrapper around
    ``solve_constrained`` for the box-only case (uses the cheaper
    closed-form box projection directly, no Dykstra iteration needed).

    Original: crb/compute_rb_sd_bc_admm1*.m (+ ~40 equivalent solver
    variants -- see module docstring)
    """
    return solve_constrained(
        cov_matrix,
        b=b,
        mu=mu,
        c=c,
        x0=x0,
        x_minus=x_minus,
        x_plus=x_plus,
        lambda_bracket=lambda_bracket,
        admm_config=admm_config,
        newton_config=newton_config,
    )

solve_constrained(cov_matrix, b=None, mu=0.0, c=1.0, x0=None, a_eq=None, b_eq=None, c_ineq=None, d_ineq=None, x_minus=None, x_plus=None, lambda_bracket=None, admm_config=None, newton_config=None, proximal_config=None)

Solve risk budgeting subject to any combination of extra linear equality constraints (a_eq @ x == b_eq), inequality constraints (c_ineq @ x <= d_ineq), and box bounds [x_minus, x_plus] -- in addition to the implicit sum(x)=1 budget constraint (enforced, as in solve_box_constrained, via the outer bisection on lambda, not via a_eq). Pass None for any constraint set you don't need.

Same ADMM structure as solve_box_constrained, generalized so the z-update projects onto the intersection of whichever constraint sets are given (via optim.proximal.proximal_linear_constraints's Dykstra alternating-projection algorithm) instead of just a box. solve_box_constrained is now a thin wrapper around this function for the box-only case (using the cheaper closed-form box projection directly, no Dykstra iteration needed, when no other constraints are given).

Original: crb/compute_rb_sd_constrained_admm1*.m (+ equivalent solver variants, same consolidation as solve_box_constrained -- see module docstring)

Note: Dykstra's alternating projection (used here whenever a_eq or c_ineq is given alongside box bounds) can, for constraint sets that meet at a sharp corner, exit on a temporary numerical plateau before reaching the true intersection point -- see the caveat in optim.proximal.proximal_linear_constraints's docstring. If a result looks suspicious, verify constraint satisfaction directly.

Source code in src/quanttoolbox/portfolio/risk_budgeting.py

def solve_constrained(
    cov_matrix: np.ndarray,
    b: np.ndarray | None = None,
    mu: np.ndarray | float = 0.0,
    c: float = 1.0,
    x0: np.ndarray | None = None,
    a_eq: np.ndarray | None = None,
    b_eq: np.ndarray | None = None,
    c_ineq: np.ndarray | None = None,
    d_ineq: np.ndarray | None = None,
    x_minus: np.ndarray | float | None = None,
    x_plus: np.ndarray | float | None = None,
    lambda_bracket: tuple[float, float] | None = None,
    admm_config: ADMMConfig | None = None,
    newton_config: NewtonConfig | None = None,
    proximal_config: ProximalConfig | None = None,
) -> RiskBudgetingResult:
    """Solve risk budgeting subject to any combination of extra linear
    equality constraints (a_eq @ x == b_eq), inequality constraints
    (c_ineq @ x <= d_ineq), and box bounds [x_minus, x_plus] -- in addition
    to the implicit sum(x)=1 budget constraint (enforced, as in
    ``solve_box_constrained``, via the outer bisection on lambda, not via
    a_eq). Pass None for any constraint set you don't need.

    Same ADMM structure as ``solve_box_constrained``, generalized so the
    z-update projects onto the *intersection* of whichever constraint sets
    are given (via ``optim.proximal.proximal_linear_constraints``'s
    Dykstra alternating-projection algorithm) instead of just a box.
    ``solve_box_constrained`` is now a thin wrapper around this function
    for the box-only case (using the cheaper closed-form box projection
    directly, no Dykstra iteration needed, when no other constraints are
    given).

    Original: crb/compute_rb_sd_constrained_admm1*.m (+ equivalent solver
    variants, same consolidation as ``solve_box_constrained`` -- see
    module docstring)

    Note: Dykstra's alternating projection (used here whenever a_eq or
    c_ineq is given alongside box bounds) can, for constraint sets that
    meet at a sharp corner, exit on a temporary numerical plateau before
    reaching the true intersection point -- see the caveat in
    ``optim.proximal.proximal_linear_constraints``'s docstring. If a
    result looks suspicious, verify constraint satisfaction directly.
    """
    cov_matrix = np.asarray(cov_matrix, dtype=float)
    n = cov_matrix.shape[0]
    mu_arr = np.zeros(n) if np.isscalar(mu) and mu == 0 else np.asarray(mu, dtype=float).flatten()
    b_arr = np.full(n, 1.0 / n) if b is None else np.asarray(b, dtype=float).flatten()
    b_arr = b_arr / np.sum(b_arr)
    x0_arr = np.full(n, 1.0 / n) if x0 is None else np.asarray(x0, dtype=float).flatten()

    admm_cfg = admm_config or ADMMConfig()
    newton_cfg = newton_config or NewtonConfig()
    prox_cfg = proximal_config or ProximalConfig()

    has_extra_constraints = a_eq is not None or c_ineq is not None
    box_only = not has_extra_constraints and (x_minus is not None or x_plus is not None)

    def _project(v: np.ndarray) -> np.ndarray:
        if has_extra_constraints:
            x_out, _ = proximal_linear_constraints(
                v,
                a_eq=a_eq,
                b_eq=b_eq,
                c_ineq=c_ineq,
                d_ineq=d_ineq,
                lb=x_minus,
                ub=x_plus,
                config=prox_cfg,
            )
            return x_out
        if box_only:
            return proximal_bounds(v, x_minus, x_plus)
        return v  # no constraints beyond the implicit budget constraint

    def _admm_at_lambda(lambda_: float) -> tuple[np.ndarray, bool, int]:
        x = x0_arr.copy()
        z = x.copy()
        z0 = z.copy()
        u = np.zeros(n)
        varphi = admm_cfg.varphi
        converged = False

        n_iter = 1
        while n_iter < admm_cfg.max_iters:
            v_x = z - u
            x, _, _ = _solve_newton(
                x, mu_arr, cov_matrix, c, b_arr, newton_cfg, lambda_, varphi, v_x
            )

            v_z = x + u
            z = _project(v_z)

            r = x - z
            s = varphi * (z - z0)
            u = u + r

            cvg = max(np.sum((x - x0_arr) ** 2), np.sum((x - z) ** 2), np.sum((z - z0) ** 2))
            if cvg <= admm_cfg.tol:
                converged = True
                break

            if admm_cfg.varphi_method == 2:
                primal_error = np.sum(r * r)
                dual_error = np.sum(s * s)
                if primal_error > admm_cfg.tau_primal * dual_error:
                    varphi = varphi * admm_cfg.tau_primal
                    u = u / admm_cfg.tau_primal
                elif dual_error > admm_cfg.tau_dual * primal_error:
                    varphi = varphi / admm_cfg.tau_dual
                    u = u * admm_cfg.tau_dual

            x0_arr[:] = x
            z0 = z.copy()
            n_iter += 1

        return x, converged, n_iter

    def _budget_gap(lambda_: np.ndarray) -> np.ndarray:
        lam = float(lambda_)
        x, _, _ = _admm_at_lambda(lam)
        return np.array(np.sum(x) - 1.0)

    if lambda_bracket is None:
        lambda_star = _auto_bracket_and_solve(_budget_gap)
    else:
        lambda_star = bisection(_budget_gap, lambda_bracket[0], lambda_bracket[1])

    if np.isnan(lambda_star):
        n_arr = np.full(n, np.nan)
        return RiskBudgetingResult(
            weights=n_arr,
            risk=np.nan,
            marginal_risk=n_arr,
            risk_contribution=n_arr,
            pct_risk_contribution=n_arr,
            converged=False,
            n_iters=0,
        )

    x, converged, n_iters = _admm_at_lambda(float(lambda_star))
    rc = risk_contribution(x, cov_matrix, mu_arr, c)
    return RiskBudgetingResult(
        weights=x,
        risk=rc.risk,
        marginal_risk=rc.marginal_risk,
        risk_contribution=rc.risk_contribution,
        pct_risk_contribution=rc.pct_risk_contribution,
        converged=converged,
        n_iters=n_iters,
    )

solve_unconstrained(cov_matrix, b=None, mu=0.0, c=1.0, x0=None, method='ccd', config=None)

Solve the classic (budget-constraint-only, sum(x)=1) risk budgeting problem: find weights x such that each asset's risk contribution matches its target budget b.

method="ccd" (default): Roncalli's cyclical coordinate descent. method="newton": Newton's method on the Lagrangian stationary conditions.

Original: rpb/{compute_rb_sd,compute_rb_sd_ccd,compute_rb_sd_newton}.m

Source code in src/quanttoolbox/portfolio/risk_budgeting.py

def solve_unconstrained(
    cov_matrix: np.ndarray,
    b: np.ndarray | None = None,
    mu: np.ndarray | float = 0.0,
    c: float = 1.0,
    x0: np.ndarray | None = None,
    method: str = "ccd",
    config: CCDConfig | NewtonConfig | None = None,
) -> RiskBudgetingResult:
    """Solve the classic (budget-constraint-only, sum(x)=1) risk budgeting
    problem: find weights x such that each asset's risk contribution
    matches its target budget b.

    method="ccd" (default): Roncalli's cyclical coordinate descent.
    method="newton": Newton's method on the Lagrangian stationary conditions.

    Original: rpb/{compute_rb_sd,compute_rb_sd_ccd,compute_rb_sd_newton}.m
    """
    cov_matrix = np.asarray(cov_matrix, dtype=float)
    n = cov_matrix.shape[0]
    mu_arr = np.zeros(n) if np.isscalar(mu) and mu == 0 else np.asarray(mu, dtype=float).flatten()
    b_arr = np.full(n, 1.0 / n) if b is None else np.asarray(b, dtype=float).flatten()
    b_arr = b_arr / np.sum(b_arr)
    x0_arr = np.full(n, 1.0 / n) if x0 is None else np.asarray(x0, dtype=float).flatten()

    if method == "ccd":
        cfg = config if isinstance(config, CCDConfig) else CCDConfig()
        x, converged, n_iters = _solve_ccd(x0_arr, mu_arr, cov_matrix, c, b_arr, cfg)
    elif method == "newton":
        cfg = config if isinstance(config, NewtonConfig) else NewtonConfig()
        x, converged, n_iters = _solve_newton(x0_arr, mu_arr, cov_matrix, c, b_arr, cfg)
        x = x / np.sum(x)
    else:
        raise ValueError(f"solve_unconstrained: unknown method '{method}' (use 'ccd' or 'newton')")

    rc = risk_contribution(x, cov_matrix, mu_arr, c)
    return RiskBudgetingResult(
        weights=x,
        risk=rc.risk,
        marginal_risk=rc.marginal_risk,
        risk_contribution=rc.risk_contribution,
        pct_risk_contribution=rc.pct_risk_contribution,
        converged=converged,
        n_iters=n_iters,
    )

solve_unconstrained_es(cov_matrix, alpha, b=None, mu=0.0, x0=None, method='ccd', config=None)

Risk budgeting under a (Gaussian) Expected Shortfall risk measure at confidence level alpha. See solve_unconstrained_var.

Original: rpb/{compute_rb_es,compute_rc_es}.m

Source code in src/quanttoolbox/portfolio/risk_budgeting.py

def solve_unconstrained_es(
    cov_matrix: np.ndarray,
    alpha: float,
    b: np.ndarray | None = None,
    mu: np.ndarray | float = 0.0,
    x0: np.ndarray | None = None,
    method: str = "ccd",
    config: CCDConfig | NewtonConfig | None = None,
) -> RiskBudgetingResult:
    """Risk budgeting under a (Gaussian) Expected Shortfall risk measure at
    confidence level alpha. See ``solve_unconstrained_var``.

    Original: rpb/{compute_rb_es,compute_rc_es}.m
    """
    return solve_unconstrained(
        cov_matrix, b=b, mu=mu, c=_es_multiplier(alpha), x0=x0, method=method, config=config
    )

solve_unconstrained_var(cov_matrix, alpha, b=None, mu=0.0, x0=None, method='ccd', config=None)

Risk budgeting under a (Gaussian) Value-at-Risk risk measure at confidence level alpha, rather than plain volatility -- solved by mapping alpha to the equivalent standard-deviation multiplier c and reusing solve_unconstrained directly (VaR is just a rescaled standard-deviation measure under normality).

Original: rpb/{compute_rb_var,compute_rc_vol}.m (VaR variant)

Source code in src/quanttoolbox/portfolio/risk_budgeting.py

def solve_unconstrained_var(
    cov_matrix: np.ndarray,
    alpha: float,
    b: np.ndarray | None = None,
    mu: np.ndarray | float = 0.0,
    x0: np.ndarray | None = None,
    method: str = "ccd",
    config: CCDConfig | NewtonConfig | None = None,
) -> RiskBudgetingResult:
    """Risk budgeting under a (Gaussian) Value-at-Risk risk measure at
    confidence level alpha, rather than plain volatility -- solved by
    mapping alpha to the equivalent standard-deviation multiplier c and
    reusing ``solve_unconstrained`` directly (VaR is just a rescaled
    standard-deviation measure under normality).

    Original: rpb/{compute_rb_var,compute_rc_vol}.m (VaR variant)
    """
    return solve_unconstrained(
        cov_matrix, b=b, mu=mu, c=_var_multiplier(alpha), x0=x0, method=method, config=config
    )

Examples

Equal risk contribution and box-constrained risk budgeting — rpb/test_box1.py

"""Translated from Examples/rpb/test_box1.m -- ERC and box-constrained
("C-ERC") risk budgeting portfolios at progressively wider bounds around
a starting position."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.portfolio.risk_budgeting import erc_portfolio, solve_box_constrained
from quanttoolbox.stats.moments import corr_to_cov

x0 = np.array([0.29, 0.25, 0.23, 0.18, 0.05])
sigma = np.array([0.20, 0.20, 0.25, 0.15, 0.25])
rho = xpnd(
    np.array(
        [1.00, 0.40, 1.00, 0.70, 0.75, 1.00, 0.60, 0.55, 0.90, 1.00, 0.70, 0.60, 0.70, 0.65, 1.00]
    ),
    method=1,
)
cov_matrix = corr_to_cov(sigma, rho)

r1 = erc_portfolio(cov_matrix)
print("ERC weights:", np.round(r1.weights, 4))

for delta in [0.02, 0.07, 0.20]:
    x_minus, x_plus = x0 - delta, x0 + delta
    r = solve_box_constrained(cov_matrix, x_minus=x_minus, x_plus=x_plus, x0=x0)
    print(f"box (delta={delta}) weights:", np.round(r.weights, 4), "converged:", r.converged)

Risk budgeting toward unequal target budgets — rpb/test_erc3.py

"""Translated from Examples/rpb/test_erc3.m -- Example 17 (page 123),
Roncalli (2013). RB portfolio with unequal target budgets."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.portfolio.risk_budgeting import risk_contribution, solve_unconstrained
from quanttoolbox.stats.moments import corr_to_cov

sigma = np.array([0.15, 0.20, 0.30, 0.10])
rho = xpnd(np.array([1.00, 0.50, 1.00, 0.00, 0.20, 1.00, -0.10, 0.40, 0.70, 1.00]), method=1)
cov_matrix = corr_to_cov(sigma, rho)
x = np.full(4, 0.25)

rc = risk_contribution(x, cov_matrix)
print("equal-weight risk contribution:", rc.risk, np.round(100 * rc.pct_risk_contribution, 2))

b = np.array([0.20, 0.20, 0.30, 0.30])
r = solve_unconstrained(cov_matrix, b=b, method="ccd")
print("RB weights (target budgets 20/20/30/30):", np.round(r.weights, 4))
print("converged:", r.converged, "n_iters:", r.n_iters)

Risk contribution decomposition and equal-budget risk budgeting — rpb/test_erc2.py

"""Translated from Examples/rpb/test_erc2.m -- Example 7 (page 80),
Roncalli (2013). Risk contribution decomposition + risk budgeting at
various target budgets."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.portfolio.risk_budgeting import risk_contribution, solve_unconstrained
from quanttoolbox.stats.moments import corr_to_cov

sigma = np.array([0.30, 0.20, 0.15])
rho = xpnd(np.array([1.00, 0.80, 1.00, 0.50, 0.30, 1.00]), method=1)
cov_matrix = corr_to_cov(sigma, rho)
x = np.array([0.50, 0.20, 0.30])

rc = risk_contribution(x, cov_matrix)
print("fixed-x risk contribution:", rc.risk, np.round(100 * rc.pct_risk_contribution, 2))

r_equal = solve_unconstrained(cov_matrix, b=np.full(3, 1 / 3), method="ccd")
print("equal-budget RB weights:", np.round(r_equal.weights, 4))

r_custom = solve_unconstrained(cov_matrix, b=x, method="ccd")
print("target-x-as-budget RB weights:", np.round(r_custom.weights, 4))

portfolio.mean_variance, portfolio.black_litterman, portfolio.tracking_error

Python alternatives

Hybrid: PyPortfolioOpt is mature and actively maintained, covering mean-variance, Black-Litterman, CVaR, and discrete allocation — genuinely worth using directly for standard portfolio construction. Keep these modules for tighter integration with the rest of this codebase (shared solve_qp, direct ridge/lasso penalty composability).

quanttoolbox.portfolio.mean_variance

Mean-variance, minimum-variance, and most-diversified portfolio construction.

Ported from QuantToolBox/rpb/{compute_mvo_portfolio,compute_minvar_portfolio, compute_mdp_portfolio,compute_mdp_objective_function}.m and QuantToolBox/mloapa/{compute_MDP_ADMM,compute_MinVar_ADMM1, compute_MinVar_ADMM2}.m.

Consolidation notes:

• mvo_portfolio/minvar_portfolio route through quanttoolbox.optim.quadprog.solve_qp (see that module's docstring for why this eliminates the original's separate quadprog_lasso/ridge/etc. variable-splitting machinery); minvar_portfolio is literally mvo_portfolio with mu=0.
• mvo_frontier evaluates at a list of risk-aversion (gamma) values, covering the original's "gamma-problem" mode (problem=0). mvo_target_portfolio covers the "mu-problem"/"sigma-problem" target-matching modes (problem=1/problem=2): for each target expected-return or target volatility, it bisects on gamma (via quanttoolbox.optim.bisection.bisection) until mvo_portfolio's achieved return/volatility hits that target, exactly mirroring compute_mvo_portfolio.m's three-branch structure (boundary cases at gamma=0/gamma=100, bisection in between over gamma in [0, 10] by default -- both bounds preserved as separate, independently configurable parameters, matching a real quirk in the original: the achievability check against the "infinite risk aversion" case uses gamma=100, but the interior bisection search is only ever bracketed to [0, 10], so a target only reachable at a gamma between 10 and 100 will come back as unreachable (NaN), exactly as the MATLAB source does. See docs/matlab_bugs_found.md for a related, genuine bug this surfaced in one of the original's own example scripts.)
• mdp_portfolio (most diversified portfolio) has a genuinely nonlinear objective (log(portfolio vol) - log(weighted-avg individual vol)), so it uses scipy.optimize.minimize (SLSQP) rather than solve_qp -- matching the original's use of fmincon rather than quadprog. mloapa/compute_MDP_ADMM.m (an alternative ADMM-based MDP solver) is not separately ported since SLSQP solves the same small nonlinear program directly and reliably.
• mloapa/compute_MinVar_ADMM{1,2}.m are alternative ADMM solvers for the same minimum-variance QP that solve_qp already solves directly; not separately ported.

mdp_portfolio(cov_matrix, a_eq=None, b_eq=None, c_ineq=None, d_ineq=None, lb=None, ub=None, x0=None)

Most Diversified Portfolio: maximize the diversification ratio (weighted-average individual volatility / portfolio volatility), i.e. minimize log(portfolio vol) - log(weighted-average individual vol).

Original: rpb/{compute_mdp_portfolio,compute_mdp_objective_function}.m

Source code in src/quanttoolbox/portfolio/mean_variance.py

def mdp_portfolio(
    cov_matrix: np.ndarray,
    a_eq: np.ndarray | None = None,
    b_eq: np.ndarray | None = None,
    c_ineq: np.ndarray | None = None,
    d_ineq: np.ndarray | None = None,
    lb: np.ndarray | float | None = None,
    ub: np.ndarray | float | None = None,
    x0: np.ndarray | None = None,
) -> MDPResult:
    """Most Diversified Portfolio: maximize the diversification ratio
    (weighted-average individual volatility / portfolio volatility), i.e.
    minimize log(portfolio vol) - log(weighted-average individual vol).

    Original: rpb/{compute_mdp_portfolio,compute_mdp_objective_function}.m
    """
    cov_matrix = np.asarray(cov_matrix, dtype=float)
    n = cov_matrix.shape[0]
    sigma = np.sqrt(np.diag(cov_matrix))

    a_eq = np.ones((1, n)) if a_eq is None else np.asarray(a_eq, dtype=float)
    b_eq = np.array([1.0]) if b_eq is None else np.asarray(b_eq, dtype=float)
    lb_arr = (
        np.full(n, -100.0)
        if lb is None
        else np.full(n, lb)
        if np.isscalar(lb)
        else np.asarray(lb, dtype=float)
    )
    ub_arr = (
        np.full(n, 100.0)
        if ub is None
        else np.full(n, ub)
        if np.isscalar(ub)
        else np.asarray(ub, dtype=float)
    )
    x0_arr = np.full(n, 1.0 / n) if x0 is None else np.asarray(x0, dtype=float)

    def objective(x: np.ndarray) -> float:
        port_vol = np.sqrt(x @ cov_matrix @ x)
        weighted_avg_vol = x @ sigma
        return float(np.log(port_vol) - np.log(weighted_avg_vol))

    constraints = [LinearConstraint(a_eq, b_eq, b_eq)]
    if c_ineq is not None:
        constraints.append(
            LinearConstraint(
                np.asarray(c_ineq, dtype=float), -np.inf, np.asarray(d_ineq, dtype=float)
            )
        )

    result = minimize(
        objective,
        x0_arr,
        method="SLSQP",
        bounds=list(zip(lb_arr, ub_arr, strict=True)),
        constraints=constraints,
        options={"maxiter": 1000, "ftol": 1e-12},
    )

    x = result.x
    x = np.where(np.abs(x) < 1e-10, 0.0, x)
    sigma_x = float(np.sqrt(x @ cov_matrix @ x))
    dr_x = float((x @ sigma) / sigma_x)

    return MDPResult(
        weights=x, volatility=sigma_x, diversification_ratio=dr_x, converged=result.success
    )

minvar_portfolio(cov_matrix, a_eq=None, b_eq=None, c_ineq=None, d_ineq=None, lb=None, ub=None)

Minimum-variance portfolio (mean-variance optimal with gamma=0).

Original: rpb/compute_minvar_portfolio.m

Source code in src/quanttoolbox/portfolio/mean_variance.py

def minvar_portfolio(
    cov_matrix: np.ndarray,
    a_eq: np.ndarray | None = None,
    b_eq: np.ndarray | None = None,
    c_ineq: np.ndarray | None = None,
    d_ineq: np.ndarray | None = None,
    lb: np.ndarray | float | None = None,
    ub: np.ndarray | float | None = None,
) -> PortfolioResult:
    """Minimum-variance portfolio (mean-variance optimal with gamma=0).

    Original: rpb/compute_minvar_portfolio.m
    """
    cov_matrix = np.asarray(cov_matrix, dtype=float)
    n = cov_matrix.shape[0]
    return mvo_portfolio(
        np.zeros(n),
        cov_matrix,
        gamma=0.0,
        a_eq=a_eq,
        b_eq=b_eq,
        c_ineq=c_ineq,
        d_ineq=d_ineq,
        lb=lb,
        ub=ub,
    )

mvo_frontier(mu, cov_matrix, gamma_values, **kwargs)

Evaluate the mean-variance frontier at each risk-aversion value in gamma_values (the "gamma-problem" mode of compute_mvo_portfolio.m).

See mvo_target_portfolio for the mu-problem/sigma-problem target-matching modes.

Source code in src/quanttoolbox/portfolio/mean_variance.py

def mvo_frontier(
    mu: np.ndarray,
    cov_matrix: np.ndarray,
    gamma_values: np.ndarray,
    **kwargs,
) -> list[PortfolioResult]:
    """Evaluate the mean-variance frontier at each risk-aversion value in
    gamma_values (the "gamma-problem" mode of compute_mvo_portfolio.m).

    See ``mvo_target_portfolio`` for the mu-problem/sigma-problem
    target-matching modes.
    """
    return [
        mvo_portfolio(mu, cov_matrix, gamma=float(g), **kwargs) for g in np.atleast_1d(gamma_values)
    ]

mvo_portfolio(mu, cov_matrix, gamma=1.0, a_eq=None, b_eq=None, c_ineq=None, d_ineq=None, lb=None, ub=None, ridge_penalty=None, lasso_penalty=None)

Mean-variance optimal portfolio: maximize gammamu'x - 0.5x'Cov*x subject to a budget constraint (sum(x)=1 by default) and any other given constraints.

gamma=0 gives the minimum-variance portfolio; larger gamma weights expected return more heavily relative to risk.

Original: rpb/compute_mvo_portfolio.m (gamma-problem branch)

Source code in src/quanttoolbox/portfolio/mean_variance.py

def mvo_portfolio(
    mu: np.ndarray,
    cov_matrix: np.ndarray,
    gamma: float = 1.0,
    a_eq: np.ndarray | None = None,
    b_eq: np.ndarray | None = None,
    c_ineq: np.ndarray | None = None,
    d_ineq: np.ndarray | None = None,
    lb: np.ndarray | float | None = None,
    ub: np.ndarray | float | None = None,
    ridge_penalty: tuple[np.ndarray | float, np.ndarray] | None = None,
    lasso_penalty: tuple[np.ndarray | float, np.ndarray] | None = None,
) -> PortfolioResult:
    """Mean-variance optimal portfolio: maximize gamma*mu'x - 0.5*x'Cov*x
    subject to a budget constraint (sum(x)=1 by default) and any other
    given constraints.

    gamma=0 gives the minimum-variance portfolio; larger gamma weights
    expected return more heavily relative to risk.

    Original: rpb/compute_mvo_portfolio.m (gamma-problem branch)
    """
    mu = np.asarray(mu, dtype=float).flatten()
    cov_matrix = np.asarray(cov_matrix, dtype=float)

    x = solve_qp(
        cov_matrix,
        gamma * mu,
        a_eq=a_eq,
        b_eq=b_eq,
        c_ineq=c_ineq,
        d_ineq=d_ineq,
        lb=lb,
        ub=ub,
        ridge_penalty=ridge_penalty,
        lasso_penalty=lasso_penalty,
        default_budget_constraint=a_eq is None,
    )
    return PortfolioResult(
        weights=x, expected_return=float(x @ mu), volatility=float(np.sqrt(x @ cov_matrix @ x))
    )

mvo_target_portfolio(mu, cov_matrix, targets, problem='sigma', a_eq=None, b_eq=None, c_ineq=None, d_ineq=None, lb=None, ub=None, gamma_bracket=(0.0, 10.0), gamma_max=100.0)

Target-matching mean-variance portfolios: for each value in targets, find the risk-aversion gamma whose mvo_portfolio solution achieves that target expected return (problem="mu") or that target volatility (problem="sigma"), via bisection on gamma.

Original: rpb/compute_mvo_portfolio.m (mu-problem/problem=1 and sigma-problem/problem=2 branches), via compute_mvo_portfolio_return.m/compute_mvo_portfolio_volatility.m as the bisection objective.

For each target, first the gamma=0 and gamma=gamma_max solutions are used to bracket what's achievable:

• problem="sigma": a target below the gamma=0 volatility is unreachable (NaN weights); a target at or above the gamma=gamma_max volatility returns that portfolio directly; otherwise gamma is bisected within gamma_bracket.
• problem="mu": a target at or below the gamma=0 return returns that portfolio directly; a target above the gamma=gamma_max return is unreachable (NaN weights); a target exactly at the gamma=gamma_max return returns that portfolio directly; otherwise gamma is bisected within gamma_bracket.

Note the original's own quirk, preserved here: the boundary checks use gamma_max (100 by default) but the bisection search itself is only ever bracketed to gamma_bracket ((0, 10) by default) -- so a target only reachable at a gamma between the two is misreported as unreachable. See the module docstring and docs/matlab_bugs_found.md for more.

Source code in src/quanttoolbox/portfolio/mean_variance.py

def mvo_target_portfolio(
    mu: np.ndarray,
    cov_matrix: np.ndarray,
    targets: np.ndarray | float,
    problem: str = "sigma",
    a_eq: np.ndarray | None = None,
    b_eq: np.ndarray | None = None,
    c_ineq: np.ndarray | None = None,
    d_ineq: np.ndarray | None = None,
    lb: np.ndarray | float | None = None,
    ub: np.ndarray | float | None = None,
    gamma_bracket: tuple[float, float] = (0.0, 10.0),
    gamma_max: float = 100.0,
) -> list[MVOTargetResult]:
    """Target-matching mean-variance portfolios: for each value in
    ``targets``, find the risk-aversion gamma whose ``mvo_portfolio``
    solution achieves that target expected return (``problem="mu"``) or
    that target volatility (``problem="sigma"``), via bisection on gamma.

    Original: rpb/compute_mvo_portfolio.m (mu-problem/problem=1 and
    sigma-problem/problem=2 branches), via
    compute_mvo_portfolio_return.m/compute_mvo_portfolio_volatility.m as
    the bisection objective.

    For each target, first the gamma=0 and gamma=``gamma_max`` solutions
    are used to bracket what's achievable:

    - ``problem="sigma"``: a target below the gamma=0 volatility is
      unreachable (NaN weights); a target at or above the gamma=gamma_max
      volatility returns that portfolio directly; otherwise gamma is
      bisected within ``gamma_bracket``.
    - ``problem="mu"``: a target at or below the gamma=0 return returns
      that portfolio directly; a target above the gamma=gamma_max return
      is unreachable (NaN weights); a target exactly at the gamma=gamma_max
      return returns that portfolio directly; otherwise gamma is bisected
      within ``gamma_bracket``.

    Note the original's own quirk, preserved here: the boundary checks use
    ``gamma_max`` (100 by default) but the bisection search itself is only
    ever bracketed to ``gamma_bracket`` ((0, 10) by default) -- so a target
    only reachable at a gamma between the two is misreported as
    unreachable. See the module docstring and
    ``docs/matlab_bugs_found.md`` for more.
    """
    if problem not in ("mu", "sigma"):
        raise ValueError('problem must be "mu" or "sigma"')

    mu = np.asarray(mu, dtype=float).flatten()
    cov_matrix = np.asarray(cov_matrix, dtype=float)
    kwargs = dict(a_eq=a_eq, b_eq=b_eq, c_ineq=c_ineq, d_ineq=d_ineq, lb=lb, ub=ub)

    r_min = mvo_portfolio(mu, cov_matrix, gamma=0.0, **kwargs)
    r_max = mvo_portfolio(mu, cov_matrix, gamma=gamma_max, **kwargs)

    def nan_result() -> MVOTargetResult:
        return MVOTargetResult(
            weights=np.full_like(mu, np.nan),
            expected_return=np.nan,
            volatility=np.nan,
            gamma=np.nan,
        )

    def as_target_result(r: PortfolioResult, gamma: float) -> MVOTargetResult:
        return MVOTargetResult(
            weights=r.weights,
            expected_return=r.expected_return,
            volatility=r.volatility,
            gamma=gamma,
        )

    def achieved(gamma: float) -> float:
        r = mvo_portfolio(mu, cov_matrix, gamma=float(gamma), **kwargs)
        return r.expected_return if problem == "mu" else r.volatility

    results = []
    for target in np.atleast_1d(np.asarray(targets, dtype=float)):
        target = float(target)

        if problem == "sigma":
            if target < r_min.volatility:
                results.append(nan_result())
                continue
            if target == r_min.volatility:
                results.append(as_target_result(r_min, 0.0))
                continue
            if target >= r_max.volatility:
                results.append(as_target_result(r_max, np.inf))
                continue
        else:  # mu-problem
            if target <= r_min.expected_return:
                results.append(as_target_result(r_min, 0.0))
                continue
            if target > r_max.expected_return:
                results.append(nan_result())
                continue
            if target == r_max.expected_return:
                results.append(as_target_result(r_max, np.inf))
                continue

        def objective(gamma: float, target: float = target) -> float:
            return achieved(gamma) - target

        gamma_star = bisection(objective, gamma_bracket[0], gamma_bracket[1])
        if np.isnan(gamma_star):
            results.append(nan_result())
        else:
            results.append(
                as_target_result(
                    mvo_portfolio(mu, cov_matrix, gamma=float(gamma_star), **kwargs), gamma_star
                )
            )

    return results

quanttoolbox.portfolio.black_litterman

Black-Litterman implied returns and posterior (view-updated) moments.

Ported from QuantToolBox/rpb/{implied_risk_premia, compute_Black_Litterman_moments}.m

black_litterman_moments(mu_tilde, gamma_matrix, p, q, omega)

Combine an equilibrium prior (mu_tilde, gamma_matrix) with investor views (P @ mu = Q, view uncertainty Omega) into posterior expected returns and covariance, via the standard Black-Litterman conditional normal update.

Original: rpb/compute_Black_Litterman_moments.m

Parameters:

Name	Type	Description	Default
mu_tilde	(n,) prior (equilibrium) expected returns.		required
gamma_matrix	(n, n) prior covariance of expected returns (often	tau * Sigma, where Sigma is the asset covariance and tau is a small scalar).	required
p	(k, n) view matrix, one row per view.		required
q	(k,) view target values.		required
omega	(k, k) view uncertainty covariance.		required

Source code in src/quanttoolbox/portfolio/black_litterman.py

def black_litterman_moments(
    mu_tilde: np.ndarray,
    gamma_matrix: np.ndarray,
    p: np.ndarray,
    q: np.ndarray,
    omega: np.ndarray,
) -> BlackLittermanMoments:
    """Combine an equilibrium prior (mu_tilde, gamma_matrix) with investor
    views (P @ mu = Q, view uncertainty Omega) into posterior expected
    returns and covariance, via the standard Black-Litterman conditional
    normal update.

    Original: rpb/compute_Black_Litterman_moments.m

    Parameters
    ----------
    mu_tilde : (n,) prior (equilibrium) expected returns.
    gamma_matrix : (n, n) prior covariance of expected returns (often
        tau * Sigma, where Sigma is the asset covariance and tau is a
        small scalar).
    p : (k, n) view matrix, one row per view.
    q : (k,) view target values.
    omega : (k, k) view uncertainty covariance.
    """
    mu_tilde = np.asarray(mu_tilde, dtype=float).flatten()
    gamma_matrix = np.asarray(gamma_matrix, dtype=float)
    p = np.asarray(p, dtype=float)
    q = np.asarray(q, dtype=float).flatten()
    omega = np.asarray(omega, dtype=float)

    pg = p @ gamma_matrix
    inner = np.linalg.inv(omega + pg @ p.T)

    mu_bar = mu_tilde + gamma_matrix @ p.T @ inner @ (q - p @ mu_tilde)
    sigma_bar = gamma_matrix - gamma_matrix @ p.T @ inner @ pg

    return BlackLittermanMoments(mu_bar=mu_bar, sigma_bar=sigma_bar)

implied_risk_premia(x, cov_matrix, sharpe_ratio)

Back out the implied excess-return vector (and risk-aversion parameters) consistent with a given market-cap-weighted portfolio x achieving Sharpe ratio sharpe_ratio -- the standard first step of the Black-Litterman framework's equilibrium prior.

Original: rpb/implied_risk_premia.m

Source code in src/quanttoolbox/portfolio/black_litterman.py

def implied_risk_premia(
    x: np.ndarray, cov_matrix: np.ndarray, sharpe_ratio: float
) -> ImpliedRiskPremia:
    """Back out the implied excess-return vector (and risk-aversion
    parameters) consistent with a given market-cap-weighted portfolio x
    achieving Sharpe ratio sharpe_ratio -- the standard first step of the
    Black-Litterman framework's equilibrium prior.

    Original: rpb/implied_risk_premia.m
    """
    x = np.asarray(x, dtype=float).flatten()
    cov_matrix = np.asarray(cov_matrix, dtype=float)

    sigma_x = float(np.sqrt(x @ cov_matrix @ x))
    phi = sharpe_ratio / sigma_x
    gamma = sigma_x / sharpe_ratio
    pi = sharpe_ratio * (cov_matrix @ x) / sigma_x

    return ImpliedRiskPremia(pi=pi, phi=phi, gamma=gamma)

quanttoolbox.portfolio.tracking_error

Tracking-error-optimized portfolio construction relative to a benchmark.

Ported from QuantToolBox/rpb/{compute_te_portfolio, compute_minimum_te_portfolio,compute_te_portfolio_mixed_norm}.m.

Consolidation notes:

• All three original functions solve the same underlying QP (minimize active variance minus a return-tilt term) via quadprog; here they route through quanttoolbox.optim.quadprog.solve_qp. te_portfolio_mixed_norm folds naturally into te_portfolio as optional ridge_penalty/lasso_penalty arguments (already supported directly by solve_qp) rather than being a separate function.
• minimum_te_portfolio is te_portfolio with gamma=0 (pure tracking-error minimization, no return tilt).
• te_frontier covers the "gamma-problem" mode (evaluate at a list of gamma values, problem=0). te_target_portfolio covers the "mu-problem"/"sigma-problem" target-matching modes (problem=1/ problem=2), mirroring compute_te_portfolio.m exactly the same way mean_variance.mvo_target_portfolio mirrors compute_mvo_portfolio.m -- see that function's docstring for the bisection-bracket/boundary-check quirk both share.

minimum_te_portfolio(x_benchmark, cov_matrix, a_eq=None, b_eq=None, c_ineq=None, d_ineq=None, lb=None, ub=None, on_infeasible='raise')

Minimum-tracking-error portfolio (te_portfolio with gamma=0, no return tilt).

Original: rpb/compute_minimum_te_portfolio.m

on_infeasible : "raise" (default) lets solve_qp's RuntimeError on infeasibility/non-convergence propagate unchanged -- existing behavior. "nearest" catches that RuntimeError and instead solves the same tracking-error QP via scipy.optimize.minimize (SLSQP), starting from the benchmark weights and using the same equality/inequality/box-constraint handling (including pruning rank-deficient equality rows, and defaulting missing lb/ub to [0, 1]) as HSF-Notebooks chapter 11h's local min_te_portfolio, returning the "nearest" near-feasible point SLSQP settles at instead of raising. Not ported from the MATLAB HSF toolbox -- no .m file in hfs-archive implements this fallback. Promoted from HSF-Notebooks chapter 11h.

Source code in src/quanttoolbox/portfolio/tracking_error.py

def minimum_te_portfolio(
    x_benchmark: np.ndarray,
    cov_matrix: np.ndarray,
    a_eq: np.ndarray | None = None,
    b_eq: np.ndarray | None = None,
    c_ineq: np.ndarray | None = None,
    d_ineq: np.ndarray | None = None,
    lb: np.ndarray | float | None = None,
    ub: np.ndarray | float | None = None,
    on_infeasible: str = "raise",
) -> TrackingErrorResult:
    """Minimum-tracking-error portfolio (te_portfolio with gamma=0, no
    return tilt).

    Original: rpb/compute_minimum_te_portfolio.m

    on_infeasible : "raise" (default) lets ``solve_qp``'s ``RuntimeError``
        on infeasibility/non-convergence propagate unchanged -- existing
        behavior. "nearest" catches that ``RuntimeError`` and instead
        solves the same tracking-error QP via ``scipy.optimize.minimize``
        (SLSQP), starting from the benchmark weights and using the same
        equality/inequality/box-constraint handling (including pruning
        rank-deficient equality rows, and defaulting missing lb/ub to
        [0, 1]) as HSF-Notebooks chapter 11h's local ``min_te_portfolio``,
        returning the "nearest" near-feasible point SLSQP settles at
        instead of raising. Not ported from the MATLAB HSF toolbox -- no
        .m file in hfs-archive implements this fallback. Promoted from
        HSF-Notebooks chapter 11h.
    """
    if on_infeasible not in ("raise", "nearest"):
        raise ValueError('on_infeasible must be "raise" or "nearest"')

    x_benchmark = np.asarray(x_benchmark, dtype=float).flatten()
    n = x_benchmark.shape[0]
    cov_matrix = np.asarray(cov_matrix, dtype=float)

    try:
        return te_portfolio(
            x_benchmark,
            np.zeros(n),
            cov_matrix,
            gamma=0.0,
            a_eq=a_eq,
            b_eq=b_eq,
            c_ineq=c_ineq,
            d_ineq=d_ineq,
            lb=lb,
            ub=ub,
        )
    except RuntimeError:
        if on_infeasible != "nearest":
            raise
        return _minimum_te_portfolio_nearest(
            x_benchmark, cov_matrix, a_eq, b_eq, c_ineq, d_ineq, lb, ub
        )

solve_mixed_norm_te_portfolio(x_benchmark, carbon_intensity, md, dts, sector, reduction, varphi_as, varphi_md, varphi_dts, lb=None, ub=None, x0=None)

L1 mixed active-share/MD/DTS bond tracking-error objective under a budget constraint and a carbon-intensity-reduction target, solved two ways: direct nonsmooth minimization (SLSQP) and an epigraph LP reformulation (linprog).

Minimizes R_Mix(w) = varphi_as*R_AS(w) + varphi_md*R_MD(w) + varphi_dts*R_DTS(w) where R_AS(w) = 0.5*||w-x_benchmark||_1, R_MD(w) = ||C_MD@(w-x_benchmark)||_1 and R_DTS(w) = ||C_DTS@(w-x_benchmark)||_1 for the per-sector aggregation matrices C_MD/C_DTS built from sector, subject to sum(w) == 1 and carbon_intensity@w <= (1-reduction) * carbon_intensity@x_benchmark. Because C_MD/C_DTS are not diagonal, these L1 terms are not expressible via solve_qp's elementwise lasso_penalty and need this dedicated two-path implementation instead.

Not ported from the MATLAB HSF toolbox -- no .m file in hfs-archive implements this. Promoted from HSF-Notebooks chapter 11b/16f.

Source code in src/quanttoolbox/portfolio/tracking_error.py

def solve_mixed_norm_te_portfolio(
    x_benchmark: np.ndarray,
    carbon_intensity: np.ndarray,
    md: np.ndarray,
    dts: np.ndarray,
    sector: np.ndarray,
    reduction: float,
    varphi_as: float,
    varphi_md: float,
    varphi_dts: float,
    lb: np.ndarray | float | None = None,
    ub: np.ndarray | float | None = None,
    x0: np.ndarray | None = None,
) -> MixedNormTEResult:
    """L1 mixed active-share/MD/DTS bond tracking-error objective under a
    budget constraint and a carbon-intensity-reduction target, solved two
    ways: direct nonsmooth minimization (SLSQP) and an epigraph LP
    reformulation (``linprog``).

    Minimizes ``R_Mix(w) = varphi_as*R_AS(w) + varphi_md*R_MD(w) +
    varphi_dts*R_DTS(w)`` where ``R_AS(w) = 0.5*||w-x_benchmark||_1``,
    ``R_MD(w) = ||C_MD@(w-x_benchmark)||_1`` and ``R_DTS(w) =
    ||C_DTS@(w-x_benchmark)||_1`` for the per-sector aggregation matrices
    ``C_MD``/``C_DTS`` built from ``sector``, subject to
    ``sum(w) == 1`` and ``carbon_intensity@w <= (1-reduction) *
    carbon_intensity@x_benchmark``. Because ``C_MD``/``C_DTS`` are not
    diagonal, these L1 terms are not expressible via ``solve_qp``'s
    elementwise ``lasso_penalty`` and need this dedicated two-path
    implementation instead.

    Not ported from the MATLAB HSF toolbox -- no .m file in hfs-archive
    implements this. Promoted from HSF-Notebooks chapter 11b/16f.
    """
    x_benchmark = np.asarray(x_benchmark, dtype=float).flatten()
    carbon_intensity = np.asarray(carbon_intensity, dtype=float).flatten()
    md = np.asarray(md, dtype=float).flatten()
    dts = np.asarray(dts, dtype=float).flatten()
    sector = np.asarray(sector)
    n = len(x_benchmark)
    sectors = np.unique(sector)
    n_sectors = len(sectors)
    s_ji = (sector[None, :] == sectors[:, None]).astype(float)  # (n_sectors, n)
    md_star = (s_ji * md) @ x_benchmark  # per-sector MD@benchmark
    dts_star = (s_ji * dts) @ x_benchmark
    ci_b = carbon_intensity @ x_benchmark

    def r_as(w: np.ndarray) -> float:
        return 0.5 * np.abs(w - x_benchmark).sum()

    def r_md(w: np.ndarray) -> float:
        return np.abs((s_ji * md) @ (w - x_benchmark)).sum()

    def r_dts(w: np.ndarray) -> float:
        return np.abs((s_ji * dts) @ (w - x_benchmark)).sum()

    def r_mix(w: np.ndarray) -> float:
        return varphi_as * r_as(w) + varphi_md * r_md(w) + varphi_dts * r_dts(w)

    if lb is None:
        lb_arr = np.zeros(n)
    elif np.isscalar(lb):
        lb_arr = np.full(n, float(lb))
    else:
        lb_arr = np.asarray(lb, dtype=float)
    if ub is None:
        ub_arr = np.ones(n)
    elif np.isscalar(ub):
        ub_arr = np.full(n, float(ub))
    else:
        ub_arr = np.asarray(ub, dtype=float)
    x0_arr = x_benchmark.copy() if x0 is None else np.asarray(x0, dtype=float).flatten()

    # -- direct nonsmooth minimization --
    cons = [
        {"type": "eq", "fun": lambda w: w.sum() - 1.0},
        {"type": "ineq", "fun": lambda w: (1 - reduction) * ci_b - carbon_intensity @ w},
    ]
    res = minimize(
        r_mix,
        x0_arr,
        method="SLSQP",
        bounds=list(zip(lb_arr, ub_arr, strict=True)),
        constraints=cons,
        options={"maxiter": 1000, "ftol": 1e-15},
    )
    w_slsqp = res.x

    # -- epigraph LP reformulation --
    # variables: [w (n), tau_w (n), tau_md (n_sectors), tau_dts (n_sectors)]
    c = np.concatenate(
        [
            np.zeros(n),
            0.5 * varphi_as * np.ones(n),
            varphi_md * np.ones(n_sectors),
            varphi_dts * np.ones(n_sectors),
        ]
    )
    i_n, i_ns = np.eye(n), np.eye(n_sectors)
    z1, z2, z3 = (
        np.zeros((n, n_sectors)),
        np.zeros((n_sectors, n)),
        np.zeros((n_sectors, n_sectors)),
    )
    c_md = s_ji * md
    c_dts = s_ji * dts
    a_ub = np.block(
        [
            [i_n, -i_n, z1, z1],
            [-i_n, -i_n, z1, z1],
            [c_md, z2, -i_ns, z3],
            [-c_md, z2, -i_ns, z3],
            [c_dts, z2, z3, -i_ns],
            [-c_dts, z2, z3, -i_ns],
        ]
    )
    b_ub_lp = np.concatenate([x_benchmark, -x_benchmark, md_star, -md_star, dts_star, -dts_star])
    a_ub = np.vstack([a_ub, np.concatenate([carbon_intensity, np.zeros(n + 2 * n_sectors)])])
    b_ub_lp = np.concatenate([b_ub_lp, [(1 - reduction) * ci_b]])
    a_eq = np.concatenate([np.ones(n), np.zeros(n + 2 * n_sectors)]).reshape(1, -1)
    b_eq = np.array([1.0])
    bounds_lp = list(zip(lb_arr, ub_arr, strict=True)) + [(0, 1e5)] * (n + 2 * n_sectors)

    res_lp = linprog(
        c, A_ub=a_ub, b_ub=b_ub_lp, A_eq=a_eq, b_eq=b_eq, bounds=bounds_lp, method="highs"
    )
    w_lp = res_lp.x[:n]

    def metrics(w: np.ndarray) -> MixedNormTEMetrics:
        ci_w = w @ carbon_intensity
        return MixedNormTEMetrics(
            active_share=r_as(w),
            md=float(w @ md),
            dts=float(w @ dts),
            r_as=r_as(w),
            r_md=r_md(w),
            r_dts=r_dts(w),
            ci=float(ci_w),
            reduction=float(1 - ci_w / ci_b),
            r_mix=r_mix(w),
        )

    return MixedNormTEResult(
        weights_slsqp=w_slsqp,
        weights_lp=w_lp,
        metrics_slsqp=metrics(w_slsqp),
        metrics_lp=metrics(w_lp),
        metrics_benchmark=metrics(x_benchmark),
        max_abs_diff=float(np.max(np.abs(w_slsqp - w_lp))),
    )

te_frontier(x_benchmark, mu, cov_matrix, gamma_values, **kwargs)

Evaluate the tracking-error frontier at each risk-aversion value in gamma_values (the "gamma-problem" mode of compute_te_portfolio.m).

Source code in src/quanttoolbox/portfolio/tracking_error.py

def te_frontier(
    x_benchmark: np.ndarray,
    mu: np.ndarray,
    cov_matrix: np.ndarray,
    gamma_values: np.ndarray,
    **kwargs,
) -> list[TrackingErrorResult]:
    """Evaluate the tracking-error frontier at each risk-aversion value in
    gamma_values (the "gamma-problem" mode of compute_te_portfolio.m).
    """
    return [
        te_portfolio(x_benchmark, mu, cov_matrix, gamma=float(g), **kwargs)
        for g in np.atleast_1d(gamma_values)
    ]

te_portfolio(x_benchmark, mu, cov_matrix, gamma=1.0, a_eq=None, b_eq=None, c_ineq=None, d_ineq=None, lb=None, ub=None, ridge_penalty=None, lasso_penalty=None)

Tracking-error-optimal portfolio: minimize 0.5(x-x_b)'Cov(x-x_b) - gamma*mu'x, i.e. trade off tracking error against active return.

gamma=0 gives pure tracking-error minimization (no return tilt).

Original: rpb/compute_te_portfolio.m (gamma-problem branch), also covers compute_te_portfolio_mixed_norm.m via ridge_penalty/lasso_penalty.

Source code in src/quanttoolbox/portfolio/tracking_error.py

def te_portfolio(
    x_benchmark: np.ndarray,
    mu: np.ndarray,
    cov_matrix: np.ndarray,
    gamma: float = 1.0,
    a_eq: np.ndarray | None = None,
    b_eq: np.ndarray | None = None,
    c_ineq: np.ndarray | None = None,
    d_ineq: np.ndarray | None = None,
    lb: np.ndarray | float | None = None,
    ub: np.ndarray | float | None = None,
    ridge_penalty: tuple[np.ndarray | float, np.ndarray] | None = None,
    lasso_penalty: tuple[np.ndarray | float, np.ndarray] | None = None,
) -> TrackingErrorResult:
    """Tracking-error-optimal portfolio: minimize 0.5*(x-x_b)'Cov*(x-x_b) -
    gamma*mu'x, i.e. trade off tracking error against active return.

    gamma=0 gives pure tracking-error minimization (no return tilt).

    Original: rpb/compute_te_portfolio.m (gamma-problem branch), also
    covers compute_te_portfolio_mixed_norm.m via ridge_penalty/lasso_penalty.
    """
    x_benchmark = np.asarray(x_benchmark, dtype=float).flatten()
    mu = np.asarray(mu, dtype=float).flatten()
    cov_matrix = np.asarray(cov_matrix, dtype=float)

    r_lin = gamma * mu + cov_matrix @ x_benchmark
    x = solve_qp(
        cov_matrix,
        r_lin,
        a_eq=a_eq,
        b_eq=b_eq,
        c_ineq=c_ineq,
        d_ineq=d_ineq,
        lb=lb,
        ub=ub,
        ridge_penalty=ridge_penalty,
        lasso_penalty=lasso_penalty,
        default_budget_constraint=a_eq is None,
    )

    active = x - x_benchmark
    return TrackingErrorResult(
        weights=x,
        active_return=float(active @ mu),
        tracking_error=float(np.sqrt(active @ cov_matrix @ active)),
    )

te_target_portfolio(x_benchmark, mu, cov_matrix, targets, problem='sigma', a_eq=None, b_eq=None, c_ineq=None, d_ineq=None, lb=None, ub=None, gamma_bracket=(0.0, 10.0), gamma_max=100.0)

Target-matching tracking-error portfolios: for each value in targets, find the risk-aversion gamma whose te_portfolio solution achieves that target active return (problem="mu") or that target tracking error (problem="sigma"), via bisection on gamma.

Original: rpb/compute_te_portfolio.m (mu-problem/problem=1 and sigma-problem/problem=2 branches), via compute_te_portfolio_return.m/compute_te_portfolio_volatility.m as the bisection objective.

Same boundary-case and bisection-bracket structure as mean_variance.mvo_target_portfolio -- see that function's docstring for the full description of the achievability checks and the gamma_max/gamma_bracket quirk they share.

Source code in src/quanttoolbox/portfolio/tracking_error.py

def te_target_portfolio(
    x_benchmark: np.ndarray,
    mu: np.ndarray,
    cov_matrix: np.ndarray,
    targets: np.ndarray | float,
    problem: str = "sigma",
    a_eq: np.ndarray | None = None,
    b_eq: np.ndarray | None = None,
    c_ineq: np.ndarray | None = None,
    d_ineq: np.ndarray | None = None,
    lb: np.ndarray | float | None = None,
    ub: np.ndarray | float | None = None,
    gamma_bracket: tuple[float, float] = (0.0, 10.0),
    gamma_max: float = 100.0,
) -> list[TETargetResult]:
    """Target-matching tracking-error portfolios: for each value in
    ``targets``, find the risk-aversion gamma whose ``te_portfolio``
    solution achieves that target active return (``problem="mu"``) or
    that target tracking error (``problem="sigma"``), via bisection on
    gamma.

    Original: rpb/compute_te_portfolio.m (mu-problem/problem=1 and
    sigma-problem/problem=2 branches), via
    compute_te_portfolio_return.m/compute_te_portfolio_volatility.m as the
    bisection objective.

    Same boundary-case and bisection-bracket structure as
    ``mean_variance.mvo_target_portfolio`` -- see that function's
    docstring for the full description of the achievability checks and
    the gamma_max/gamma_bracket quirk they share.
    """
    if problem not in ("mu", "sigma"):
        raise ValueError('problem must be "mu" or "sigma"')

    x_benchmark = np.asarray(x_benchmark, dtype=float).flatten()
    mu = np.asarray(mu, dtype=float).flatten()
    cov_matrix = np.asarray(cov_matrix, dtype=float)
    kwargs = dict(a_eq=a_eq, b_eq=b_eq, c_ineq=c_ineq, d_ineq=d_ineq, lb=lb, ub=ub)

    r_min = te_portfolio(x_benchmark, mu, cov_matrix, gamma=0.0, **kwargs)
    r_max = te_portfolio(x_benchmark, mu, cov_matrix, gamma=gamma_max, **kwargs)

    def nan_result() -> TETargetResult:
        return TETargetResult(
            weights=np.full_like(x_benchmark, np.nan),
            active_return=np.nan,
            tracking_error=np.nan,
            gamma=np.nan,
        )

    def as_target_result(r: TrackingErrorResult, gamma: float) -> TETargetResult:
        return TETargetResult(
            weights=r.weights,
            active_return=r.active_return,
            tracking_error=r.tracking_error,
            gamma=gamma,
        )

    def achieved(gamma: float) -> float:
        r = te_portfolio(x_benchmark, mu, cov_matrix, gamma=float(gamma), **kwargs)
        return r.active_return if problem == "mu" else r.tracking_error

    results = []
    for target in np.atleast_1d(np.asarray(targets, dtype=float)):
        target = float(target)

        if problem == "sigma":
            if target < r_min.tracking_error:
                results.append(nan_result())
                continue
            if target == r_min.tracking_error:
                results.append(as_target_result(r_min, 0.0))
                continue
            if target >= r_max.tracking_error:
                results.append(as_target_result(r_max, np.inf))
                continue
        else:  # mu-problem
            if target <= r_min.active_return:
                results.append(as_target_result(r_min, 0.0))
                continue
            if target > r_max.active_return:
                results.append(nan_result())
                continue
            if target == r_max.active_return:
                results.append(as_target_result(r_max, np.inf))
                continue

        def objective(gamma: float, target: float = target) -> float:
            return achieved(gamma) - target

        gamma_star = bisection(objective, gamma_bracket[0], gamma_bracket[1])
        if np.isnan(gamma_star):
            results.append(nan_result())
        else:
            results.append(
                as_target_result(
                    te_portfolio(x_benchmark, mu, cov_matrix, gamma=float(gamma_star), **kwargs),
                    gamma_star,
                )
            )

    return results

Examples

Black-Litterman view matched to six tracking-error targets — rpb/test_bl4.py

"""Translated from Examples/rpb/test_bl4.m -- Roncalli [2013], "Introduction
to Risk Parity and Budgeting", page 24: a single Black-Litterman view
(same P/Q/Omega/tau as test_bl3.py's scenario 1) turned into a *tracking-error*
frontier against the x0 benchmark -- the sigma-problem branch of
`compute_te_portfolio.m`, now covered by `te_target_portfolio`. Six target
tracking-error levels (effectively 0% through 5%) show how active return and
the information ratio (alpha/te, here genuinely dimensionless -- unlike
test_bl3.py's alpha/sigma_x convention) scale with the amount of tracking
error taken."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.portfolio.black_litterman import black_litterman_moments, implied_risk_premia
from quanttoolbox.portfolio.tracking_error import te_target_portfolio
from quanttoolbox.stats.moments import corr_to_cov

mu = np.array([0.05, 0.06, 0.08, 0.06])
sigma = np.array([0.15, 0.20, 0.25, 0.30])
rho = xpnd(np.array([1.00, 0.10, 1.00, 0.40, 0.70, 1.00, 0.50, 0.40, 0.80, 1.00]), method=1)
cov_matrix = corr_to_cov(sigma, rho)

x0 = np.array([0.40, 0.30, 0.20, 0.10])
r = 0.03
irp = implied_risk_premia(x0, cov_matrix, 0.25)
mu_tilde = r + irp.pi

p_matrix = np.array([[1, 0, 0, 0], [0, 1, -1, 0]], dtype=float)
q = np.array([0.04, -0.01])
omega = np.diag([0.10**2, 0.05**2])
bl = black_litterman_moments(mu_tilde, 1.0 * cov_matrix, p_matrix, q, omega)

te_targets = np.array([1e-5, 0.01, 0.02, 0.03, 0.04, 0.05])
results = te_target_portfolio(
    x0, bl.mu_bar, cov_matrix, te_targets, problem="sigma", lb=0.0, ub=1.0
)
for target, res in zip(te_targets, results, strict=False):
    ir = res.active_return / res.tracking_error if res.tracking_error > 0 else float("nan")
    print(
        f"target_te={100 * target:5.2f}  gamma={res.gamma:6.3f}  w={np.round(100 * res.weights, 2)}  "
        f"alpha={100 * res.active_return:6.3f}  te={100 * res.tracking_error:5.2f}  IR={ir:.3f}"
    )

Black-Litterman view sensitivity across five scenarios — rpb/test_bl2.py

"""Translated from Examples/rpb/test_bl2.m -- Black-Litterman sensitivity
analysis across 6 scenarios (base case + 5 view/uncertainty/tau variants),
each solved as a fixed-risk-aversion MVO portfolio."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.portfolio.black_litterman import black_litterman_moments, implied_risk_premia
from quanttoolbox.portfolio.mean_variance import mvo_portfolio
from quanttoolbox.stats.moments import corr_to_cov

sigma = np.array([0.15, 0.20, 0.25, 0.30])
rho = xpnd(np.array([1.00, 0.10, 1.00, 0.40, 0.70, 1.00, 0.50, 0.40, 0.80, 1.00]), method=1)
cov_matrix = corr_to_cov(sigma, rho)

x0 = np.array([0.40, 0.30, 0.20, 0.10])
r = 0.03
irp = implied_risk_premia(x0, cov_matrix, sharpe_ratio=0.25)
mu_tilde = r + irp.pi
gamma0 = irp.gamma

scenarios = [
    dict(
        P=np.array([[1, 0, 0, 0], [0, 1, -1, 0]], dtype=float),
        Q=np.array([0.04, -0.01]),
        Omega=np.diag([0.10**2, 0.05**2]),
        tau=1,
    ),
    dict(
        P=np.array([[1, 0, 0, 0], [0, 1, -1, 0]], dtype=float),
        Q=np.array([0.07, -0.01]),
        Omega=np.diag([0.10**2, 0.05**2]),
        tau=1,
    ),
    dict(
        P=np.array([[1, 0, 0, 0], [0, 1, -1, 0]], dtype=float),
        Q=np.array([0.04, -0.01]),
        Omega=np.diag([0.20**2, 0.20**2]),
        tau=1,
    ),
    dict(
        P=np.array([[1, 0, 0, 0], [0, 1, -1, 0]], dtype=float),
        Q=np.array([0.04, -0.01]),
        Omega=np.diag([0.10**2, 0.05**2]),
        tau=0.10,
    ),
    dict(
        P=np.array([[1, 0, 0, 0], [0, 1, -1, 0]], dtype=float),
        Q=np.array([0.04, -0.01]),
        Omega=np.diag([0.10**2, 0.05**2]),
        tau=0.01,
    ),
]

results = [
    dict(weights=x0, mu=x0 @ mu_tilde, sigma=np.sqrt(x0 @ cov_matrix @ x0), alpha=0.0, te=0.0)
]
for s in scenarios:
    bl = black_litterman_moments(mu_tilde, s["tau"] * cov_matrix, s["P"], s["Q"], s["Omega"])
    mvo = mvo_portfolio(bl.mu_bar, cov_matrix, gamma=gamma0, lb=0.0, ub=1.0)
    alpha = (mvo.weights - x0) @ bl.mu_bar
    te = np.sqrt((mvo.weights - x0) @ cov_matrix @ (mvo.weights - x0))
    results.append(
        dict(weights=mvo.weights, mu=mvo.expected_return, sigma=mvo.volatility, alpha=alpha, te=te)
    )

for i, r_ in enumerate(results):
    print(
        f"scenario {i}: weights={np.round(r_['weights'],4)} mu={round(r_['mu'],5)} te={round(r_['te'],5)}"
    )

Black-Litterman views matched to a benchmark's volatility — rpb/test_bl3.py

"""Translated from Examples/rpb/test_bl3.m -- Roncalli [2013], "Introduction
to Risk Parity and Budgeting", page 24: five Black-Litterman view scenarios
(varying the view itself, its confidence Omega, and the scaling tau), each
turned into a portfolio matched to the *same* target volatility as the
original x0 benchmark -- the sigma-problem branch of
`compute_mvo_portfolio.m`, now covered by `mvo_target_portfolio`. For each
scenario this reports the resulting weights, expected return, volatility
(pinned to sigma0 by construction), active return (alpha) and information
ratio relative to x0.

Note on the displayed IR: the original computes
`IR_x(i) = alpha_x(i)/sigma_x(i)` (active return over the portfolio's *own*
volatility, not tracking error -- both are dimensionless ratios), and then
the whole results matrix is multiplied by 100 for percent-style display
before formatting. Since IR_x is already a ratio, that display convention
ends up printing 100*IR rather than IR itself; this translation keeps that
same convention (labelled `100*IR`) so the printed numbers match the
original table (Roncalli page 26) directly."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.portfolio.black_litterman import black_litterman_moments, implied_risk_premia
from quanttoolbox.portfolio.mean_variance import mvo_target_portfolio
from quanttoolbox.stats.moments import corr_to_cov

mu = np.array([0.05, 0.06, 0.08, 0.06])
sigma = np.array([0.15, 0.20, 0.25, 0.30])
rho = xpnd(np.array([1.00, 0.10, 1.00, 0.40, 0.70, 1.00, 0.50, 0.40, 0.80, 1.00]), method=1)
cov_matrix = corr_to_cov(sigma, rho)

x0 = np.array([0.40, 0.30, 0.20, 0.10])
sigma0 = np.sqrt(x0 @ cov_matrix @ x0)
r = 0.03
irp = implied_risk_premia(x0, cov_matrix, 0.25)
mu_tilde = r + irp.pi

print(f"x0={np.round(100 * x0, 2)}  mu0={100 * (x0 @ mu_tilde):.2f}  sigma0={100 * sigma0:.2f}")

# (P, Q, Omega, tau) for the 5 view scenarios
p_matrix = np.array([[1, 0, 0, 0], [0, 1, -1, 0]], dtype=float)
scenarios = [
    (np.array([0.04, -0.01]), np.diag([0.10**2, 0.05**2]), 1.0),
    (np.array([0.07, -0.01]), np.diag([0.10**2, 0.05**2]), 1.0),
    (np.array([0.04, -0.01]), np.diag([0.20**2, 0.20**2]), 1.0),
    (np.array([0.04, -0.01]), np.diag([0.10**2, 0.05**2]), 0.10),
    (np.array([0.04, -0.01]), np.diag([0.10**2, 0.05**2]), 0.01),
]

for i, (q, omega, tau) in enumerate(scenarios, start=1):
    bl = black_litterman_moments(mu_tilde, tau * cov_matrix, p_matrix, q, omega)
    res = mvo_target_portfolio(bl.mu_bar, cov_matrix, sigma0, problem="sigma", lb=0.0, ub=1.0)[0]
    alpha = (res.weights - x0) @ bl.mu_bar
    te = np.sqrt((res.weights - x0) @ cov_matrix @ (res.weights - x0))
    ir = alpha / res.volatility
    print(
        f"scenario {i}: gamma={res.gamma:.3f}  w={np.round(100 * res.weights, 2)}  "
        f"mu={100 * res.expected_return:.2f}  sigma={100 * res.volatility:.2f}  "
        f"alpha={100 * alpha:.2f}  te={100 * te:.2f}  100*IR={100 * ir:.2f}"
    )

Full-view Black-Litterman with tracking-error target-matching — rpb/test_bl5.py

"""Translated from Examples/rpb/test_bl5.m -- a second, independent 4-asset
Black-Litterman + tracking-error worked example (no book page cited in the
original). Every asset gets a direct view (`P=I`, `Q=mu`, full confidence
`Omega=Sigma_epsilon=covMatrix`) -- an unusually strong-conviction setup
that pulls the posterior most of the way from the equilibrium prior toward
the raw sample means.

Four blocks, each comparing the resulting active portfolio (x) against the
x0 benchmark:

1. gamma-problem MVO (mu_BL, Sigma_BL) at the implied risk-aversion gamma0.
2. sigma-problem TE-target=1% (mu_BL, Sigma_BL), i.e. `te_target_portfolio`.
3. gamma-problem MVO (mu_BL, covMatrix) -- same posterior mean, but
   optimized against the *original* (pre-view) covariance rather than the
   BL posterior covariance Sigma_BL.
4. sigma-problem TE-target=1% (mu_BL, covMatrix).

The original's final block repeats block 2 verbatim (mu_BL, Sigma_BL again)
-- included here too for a faithful translation, and it does reproduce
identical numbers, as expected."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.portfolio.black_litterman import black_litterman_moments, implied_risk_premia
from quanttoolbox.portfolio.mean_variance import mvo_portfolio
from quanttoolbox.portfolio.tracking_error import te_target_portfolio
from quanttoolbox.stats.moments import corr_to_cov

mu = np.array([0.03, 0.03, 0.08, 0.07])
sigma = np.array([0.06, 0.07, 0.18, 0.17])
rho = xpnd(np.array([1.00, 0.50, 1.00, -0.40, -0.40, 1.00, -0.40, -0.40, 0.80, 1.00]), method=1)
cov_matrix = corr_to_cov(sigma, rho)
print("Sigma_hat =\n", np.round(100 * cov_matrix, 2))

x0 = np.array([0.40, 0.40, 0.10, 0.10])
r = 0.02
irp = implied_risk_premia(x0, cov_matrix, 0.25)
mu_tilde = r + irp.pi
print("tilde(pi) =", np.round(100 * irp.pi, 2))
print("tilde(mu) =", np.round(100 * mu_tilde, 2))
print("gamma0 =", round(irp.gamma, 4))

p_matrix = np.eye(4)
bl = black_litterman_moments(mu_tilde, cov_matrix, p_matrix, mu, cov_matrix)
print("mu(BL) =", np.round(100 * bl.mu_bar, 2))
print("Sigma(BL) =\n", np.round(100 * bl.sigma_bar, 2))


def report(label: str, weights: np.ndarray, cov_for_te: np.ndarray) -> None:
    active = weights - x0
    te = np.sqrt(active @ cov_for_te @ active)
    print(f"{label}: x0={np.round(x0, 4)}  x={np.round(weights, 4)}  te={100 * te:.2f}")


gamma0 = irp.gamma

res = mvo_portfolio(
    bl.mu_bar,
    bl.sigma_bar,
    gamma=gamma0,
    a_eq=np.ones((1, 4)),
    b_eq=np.array([1.0]),
    lb=0.0,
    ub=1.0,
)
report("1. MVO(mu_BL, Sigma_BL) @ gamma0", res.weights, bl.sigma_bar)

te_res = te_target_portfolio(
    x0,
    bl.mu_bar,
    bl.sigma_bar,
    0.01,
    problem="sigma",
    a_eq=np.ones((1, 4)),
    b_eq=np.array([1.0]),
    lb=0.0,
    ub=1.0,
)[0]
report("2. TE-target 1% (mu_BL, Sigma_BL)", te_res.weights, bl.sigma_bar)

res2 = mvo_portfolio(
    bl.mu_bar, cov_matrix, gamma=gamma0, a_eq=np.ones((1, 4)), b_eq=np.array([1.0]), lb=0.0, ub=1.0
)
report("3. MVO(mu_BL, covMatrix) @ gamma0", res2.weights, cov_matrix)

te_res2 = te_target_portfolio(
    x0,
    bl.mu_bar,
    cov_matrix,
    0.01,
    problem="sigma",
    a_eq=np.ones((1, 4)),
    b_eq=np.array([1.0]),
    lb=0.0,
    ub=1.0,
)[0]
report("4. TE-target 1% (mu_BL, covMatrix)", te_res2.weights, cov_matrix)

te_res3 = te_target_portfolio(
    x0,
    bl.mu_bar,
    bl.sigma_bar,
    0.01,
    problem="sigma",
    a_eq=np.ones((1, 4)),
    b_eq=np.array([1.0]),
    lb=0.0,
    ub=1.0,
)[0]
report("5. TE-target 1% (mu_BL, Sigma_BL), repeated", te_res3.weights, bl.sigma_bar)

Mean-variance frontier: risk-aversion, return, and volatility targets — rpb/test_mvo2.py

"""Translated from Examples/rpb/test_mvo2.m -- Roncalli [2013], "Introduction
to Risk Parity and Budgeting", Example 1 (pages 7-8): the same 4-asset
mean-variance problem evaluated three ways -- the gamma-problem (pick a
risk-aversion, solve directly), the mu-problem (pick a target expected
return, bisect on gamma to hit it), and the sigma-problem (pick a target
volatility, bisect on gamma to hit it). All three route through
`compute_mvo_portfolio.m`'s three branches; here that's
`mvo_frontier`/`mvo_target_portfolio`.

The original passes `lb=0, ub=0` (MATLAB's "use the default -100/100 wide
bounds" sentinel); passed through explicitly here as `lb=-100, ub=100`."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.portfolio.mean_variance import mvo_frontier, mvo_target_portfolio
from quanttoolbox.stats.moments import corr_to_cov

mu = np.array([0.05, 0.06, 0.08, 0.06])
sigma = np.array([0.15, 0.20, 0.25, 0.30])
rho = xpnd(np.array([1.00, 0.10, 1.00, 0.40, 0.70, 1.00, 0.50, 0.40, 0.80, 1.00]), method=1)
cov_matrix = corr_to_cov(sigma, rho)

print("1. gamma-problem (page 7)")
gamma_values = np.array([0.00, 0.20, 0.50, 1.00, 2.00, 5.00])
results = mvo_frontier(mu, cov_matrix, gamma_values, lb=-100.0, ub=100.0)
for g, r in zip(gamma_values, results, strict=False):
    print(
        f"  gamma={g:5.2f}  mu={100 * r.expected_return:6.2f}  sigma={100 * r.volatility:6.2f}  "
        f"w={np.round(100 * r.weights, 2)}"
    )

print("\n2. mu-problem (page 8)")
mu_targets = np.array([5.00, 6.00, 7.00, 8.00, 9.00]) / 100
mu_results = mvo_target_portfolio(mu, cov_matrix, mu_targets, problem="mu", lb=-100.0, ub=100.0)
for target, r in zip(mu_targets, mu_results, strict=False):
    print(
        f"  target_mu={100 * target:5.2f}  gamma={r.gamma:6.3f}  mu={100 * r.expected_return:6.2f}  "
        f"sigma={100 * r.volatility:6.2f}  w={np.round(100 * r.weights, 2)}"
    )

print("\n3. sigma-problem (page 8)")
sigma_targets = np.array([15.00, 20.00, 25.00, 30.00, 35.00]) / 100
sigma_results = mvo_target_portfolio(
    mu, cov_matrix, sigma_targets, problem="sigma", lb=-100.0, ub=100.0
)
for target, r in zip(sigma_targets, sigma_results, strict=False):
    print(
        f"  target_sigma={100 * target:5.2f}  gamma={r.gamma:6.3f}  mu={100 * r.expected_return:6.2f}  "
        f"sigma={100 * r.volatility:6.2f}  w={np.round(100 * r.weights, 2)}"
    )

Mean-variance optimization plus ridge/lasso-penalized portfolios — rpb/test_lasso1.py

"""Translated from Examples/rpb/test_lasso1.m -- Roncalli [2013],
"Introduction to Risk Parity and Budgeting", Example 1 (page 53):
compares an unconstrained gamma-problem MVO portfolio against the same
problem with a budget (sum-to-1) constraint, and ridge/lasso-penalized
variants (toward 0 and toward equal-weight) of the unconstrained
problem.

`compute_mvo_portfolio`/`quadprog_ridge`/`quadprog_lasso` map to
`mvo_portfolio`/`solve_qp(..., ridge_penalty=...)`/`solve_qp(...,
lasso_penalty=...)` as established in test_lasso3.py/test_lasso5.py; the
original's `0, 0` sentinel arguments for "no equality constraint" /
"no bounds" map to simply omitting those keyword arguments."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.optim.quadprog import solve_qp
from quanttoolbox.portfolio.mean_variance import mvo_portfolio
from quanttoolbox.stats.moments import corr_to_cov

mu = np.array([0.05, 0.06, 0.08, 0.06])
sigma = np.array([0.15, 0.20, 0.25, 0.30])
rho = xpnd(np.array([1.00, 0.10, 1.00, 0.40, 0.70, 1.00, 0.50, 0.40, 0.80, 1.00]), method=1)
cov_matrix = corr_to_cov(sigma, rho)
n = 4
x0 = np.full(n, 1 / n)

# Case gamma-problem
gamma_x = 0.5

x1 = mvo_portfolio(mu, cov_matrix, gamma=gamma_x).weights
x2 = mvo_portfolio(
    mu, cov_matrix, gamma=gamma_x, a_eq=np.ones((1, n)), b_eq=np.array([1.0])
).weights

lambda_ridge = 0.03
s_ridge = lambda_ridge * np.eye(n)
x3 = solve_qp(cov_matrix, gamma_x * mu, ridge_penalty=(s_ridge, np.zeros(n)))
x4 = solve_qp(cov_matrix, gamma_x * mu, ridge_penalty=(s_ridge, x0))

lambda_lasso = 0.03 / 2
s_lasso = lambda_lasso * np.ones(n)
x5 = solve_qp(cov_matrix, gamma_x * mu, lasso_penalty=(s_lasso, np.zeros(n)))
x6 = solve_qp(cov_matrix, gamma_x * mu, lasso_penalty=(s_lasso, x0))

results = 100 * np.column_stack([x1, x2, x3, x4, x5, x6])
print("            x1      x2      x3      x4      x5      x6")
print(np.round(results, 2))

Minimum-variance portfolio under general linear constraints — rpb/test_minvar2.py

"""Translated from Examples/rpb/test_minvar2.m -- minimum-variance
portfolio with general linear equality/inequality constraints and box
bounds."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.portfolio.mean_variance import minvar_portfolio
from quanttoolbox.stats.moments import corr_to_cov

sigma = np.array([0.15, 0.20, 0.25, 0.30])
rho = xpnd(np.array([1.00, 0.10, 1.00, 0.40, 0.70, 1.00, 0.50, 0.40, 0.80, 1.00]), method=1)
cov_matrix = corr_to_cov(sigma, rho)

a_eq = np.array([[1.0, 1.0, 1.0, 1.0], [1.0, 1.0, 0.0, 0.0]])
b_eq = np.array([1.0, 0.0])
c_ineq = np.array([[0.0, 0.0, 0.0, -1.0]])
d_ineq = np.array([-0.90])

r = minvar_portfolio(
    cov_matrix, a_eq=a_eq, b_eq=b_eq, c_ineq=c_ineq, d_ineq=d_ineq, lb=-1.50, ub=2.00
)
print("weights:", np.round(r.weights, 3))
print("volatility:", round(r.volatility, 5))

Mixed ridge+lasso penalties toward two different target vectors — rpb/test_lasso5.py

"""Translated from Examples/rpb/test_lasso5.m -- ridge/lasso/mixed
portfolios with two different ridge/lasso target vectors (equal-weight
and a custom 20/20/30/30 target)."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.optim.quadprog import solve_qp
from quanttoolbox.portfolio.mean_variance import mvo_portfolio
from quanttoolbox.stats.moments import corr_to_cov

mu = np.array([0.05, 0.06, 0.08, 0.06])
sigma = np.array([0.15, 0.20, 0.25, 0.30])
rho = xpnd(np.array([1.00, 0.10, 1.00, 0.40, 0.70, 1.00, 0.50, 0.40, 0.80, 1.00]), method=1)
cov_matrix = corr_to_cov(sigma, rho)
n = 4
a_eq, b_eq = np.ones((1, n)), np.array([1.0])
gamma_x = 0.5

y1 = np.full(n, 1 / n)
y2 = np.array([0.20, 0.20, 0.30, 0.30])
S_ridge = np.diag(np.diag(cov_matrix))
lambda_lasso = 0.005

x1 = mvo_portfolio(mu, cov_matrix, gamma=gamma_x, a_eq=a_eq, b_eq=b_eq, lb=0.0, ub=1.0).weights
x2 = solve_qp(
    cov_matrix, gamma_x * mu, a_eq=a_eq, b_eq=b_eq, lb=0.0, ub=1.0, ridge_penalty=(S_ridge, y1)
)
x4 = solve_qp(
    cov_matrix, gamma_x * mu, a_eq=a_eq, b_eq=b_eq, lb=0.0, ub=1.0, lasso_penalty=(lambda_lasso, y1)
)
x6 = solve_qp(
    cov_matrix,
    gamma_x * mu,
    a_eq=a_eq,
    b_eq=b_eq,
    lb=0.0,
    ub=1.0,
    ridge_penalty=(S_ridge, y1),
    lasso_penalty=(lambda_lasso, y1),
)
# mixed with DIFFERENT targets for ridge (toward y1) vs lasso (toward y2)
x8 = solve_qp(
    cov_matrix,
    gamma_x * mu,
    a_eq=a_eq,
    b_eq=b_eq,
    lb=0.0,
    ub=1.0,
    ridge_penalty=(S_ridge, y1),
    lasso_penalty=(lambda_lasso, y2),
)

for name, x in [
    ("MVO", x1),
    ("Ridge->y1", x2),
    ("Lasso->y1", x4),
    ("Mixed(ridge->y1,lasso->y1)", x6),
    ("Mixed(ridge->y1,lasso->y2)", x8),
]:
    print(f"{name}: {np.round(x, 4)}")

Raw QP solve vs. te_portfolio, plus a gamma-recovery check — rpb/test_bl6.py

"""Translated from Examples/rpb/test_bl6.m -- same 4-asset setup as
test_bl5.py, demonstrating that a raw QP call and `te_portfolio` at a fixed
gamma agree, then cross-checking that fixed-gamma solution against
`te_target_portfolio`'s sigma-problem mode targeting the tracking error that
fixed-gamma solution actually achieves.

The original's raw `quadprog(H, f, ..., A, B, lb, ub, x0, options)` call
with `H=covMatrix`, `f=-gamma0*mu - covMatrix*x0` is exactly what
`te_portfolio(x0, mu, covMatrix, gamma=gamma0, ...)` computes internally
(see tracking_error.py's `r_lin = gamma*mu + cov_matrix @ x_benchmark`) --
so it's translated directly as a `te_portfolio` call rather than a raw
`solve_qp` call, and the two are verified to agree below."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.optim.quadprog import solve_qp
from quanttoolbox.portfolio.tracking_error import te_portfolio, te_target_portfolio
from quanttoolbox.stats.moments import corr_to_cov

mu = np.array([0.03, 0.03, 0.08, 0.07])
sigma = np.array([0.06, 0.07, 0.18, 0.17])
rho = xpnd(np.array([1.00, 0.50, 1.00, -0.40, -0.40, 1.00, -0.40, -0.40, 0.80, 1.00]), method=1)
cov_matrix = corr_to_cov(sigma, rho)
print("Sigma_hat =\n", np.round(100 * cov_matrix, 2))

x0 = np.array([0.40, 0.40, 0.10, 0.10])
gamma0 = 0.0422

# Raw quadprog(H, f, [], [], A, B, lb, ub, x0, options) call from the
# original, translated literally via solve_qp.
x_raw = solve_qp(
    cov_matrix,
    gamma0 * mu + cov_matrix @ x0,
    a_eq=np.ones((1, 4)),
    b_eq=np.array([1.0]),
    lb=0.0,
    ub=1.0,
)
print("x (raw quadprog) =", np.round(x_raw, 6))

# Same QP via te_portfolio -- should match x_raw exactly.
res = te_portfolio(
    x0, mu, cov_matrix, gamma=gamma0, a_eq=np.ones((1, 4)), b_eq=np.array([1.0]), lb=0.0, ub=1.0
)
te_x = np.sqrt((res.weights - x0) @ cov_matrix @ (res.weights - x0))
ir_x = res.active_return / te_x
print("x0 =", x0, " x (te_portfolio) =", np.round(res.weights, 6))
print(f"te(x) = {100 * te_x:.2f}   IR = {ir_x:.4f}")

# Cross-check: te_target_portfolio's sigma-problem, targeting the tracking
# error the fixed-gamma solution above actually achieves, should recover
# gamma0 (and the same portfolio) via bisection.
target_res = te_target_portfolio(
    x0,
    mu,
    cov_matrix,
    te_x,
    problem="sigma",
    a_eq=np.ones((1, 4)),
    b_eq=np.array([1.0]),
    lb=0.0,
    ub=1.0,
)[0]
print(
    f"cross-check: gamma from te_target_portfolio={target_res.gamma:.4f} "
    f"(vs gamma0={gamma0}), weights match={np.allclose(target_res.weights, res.weights, atol=1e-4)}"
)

Ridge, lasso, and mixed-norm penalized portfolios — rpb/test_lasso3.py

"""Translated from Examples/rpb/test_lasso3.m (byte-identical to
test_lasso4.m) -- MVO/ridge/lasso/mixed-penalty portfolios via the
consolidated solve_qp (replacing the original's separate quadprog_ridge/
quadprog_lasso/quadprog_mixed calls -- see optim/quadprog.py docstring)."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.optim.quadprog import solve_qp
from quanttoolbox.portfolio.mean_variance import mvo_portfolio
from quanttoolbox.stats.moments import corr_to_cov

mu = np.array([0.05, 0.06, 0.08, 0.06])
sigma = np.array([0.15, 0.20, 0.25, 0.30])
rho = xpnd(np.array([1.00, 0.10, 1.00, 0.40, 0.70, 1.00, 0.50, 0.40, 0.80, 1.00]), method=1)
cov_matrix = corr_to_cov(sigma, rho)
n = 4
x0 = np.full(n, 1 / n)
gamma_x = 0.5
a_eq, b_eq = np.ones((1, n)), np.array([1.0])

# gamma-problem (plain MVO)
r1 = mvo_portfolio(mu, cov_matrix, gamma=gamma_x, a_eq=a_eq, b_eq=b_eq, lb=0.0, ub=1.0)
x1 = r1.weights

# ridge-problem: penalize toward zero, scaled by each asset's own variance
S_ridge = np.diag(np.diag(cov_matrix))
x2 = solve_qp(
    cov_matrix, gamma_x * mu, a_eq=a_eq, b_eq=b_eq, lb=0.0, ub=1.0, ridge_penalty=(S_ridge, x0)
)

# lasso-problem: L1 penalty toward equal weight
lambda_lasso = 0.005
x3 = solve_qp(
    cov_matrix, gamma_x * mu, a_eq=a_eq, b_eq=b_eq, lb=0.0, ub=1.0, lasso_penalty=(lambda_lasso, x0)
)

# mixed-problem: both ridge and lasso, toward various targets
x4 = solve_qp(
    cov_matrix,
    gamma_x * mu,
    a_eq=a_eq,
    b_eq=b_eq,
    lb=0.0,
    ub=1.0,
    ridge_penalty=(S_ridge, np.zeros(n)),
    lasso_penalty=(lambda_lasso, np.zeros(n)),
)
x5 = solve_qp(
    cov_matrix,
    gamma_x * mu,
    a_eq=a_eq,
    b_eq=b_eq,
    lb=0.0,
    ub=1.0,
    ridge_penalty=(S_ridge, x0),
    lasso_penalty=(lambda_lasso, np.zeros(n)),
)
x6 = solve_qp(
    cov_matrix,
    gamma_x * mu,
    a_eq=a_eq,
    b_eq=b_eq,
    lb=0.0,
    ub=1.0,
    ridge_penalty=(S_ridge, np.zeros(n)),
    lasso_penalty=(lambda_lasso, x0),
)
x7 = solve_qp(
    cov_matrix,
    gamma_x * mu,
    a_eq=a_eq,
    b_eq=b_eq,
    lb=0.0,
    ub=1.0,
    ridge_penalty=(S_ridge, x0),
    lasso_penalty=(lambda_lasso, x0),
)

for name, x in [
    ("MVO", x1),
    ("Ridge", x2),
    ("Lasso", x3),
    ("Mixed(0,0)", x4),
    ("Mixed(EW,0)", x5),
    ("Mixed(0,EW)", x6),
    ("Mixed(EW,EW)", x7),
]:
    print(f"{name}: {np.round(x, 4)}")

Ridge/lasso penalty sweep, points from a 250-point scan — rpb/test_lasso2.py

"""Translated from Examples/rpb/test_lasso2.m -- Roncalli [2013],
"Introduction to Risk Parity and Budgeting", Example 1 (page 53): the
original sweeps 250 ridge/lasso penalty strengths and plots the resulting
weight paths (4 subplots: ridge toward zero with a budget constraint,
ridge toward equal-weight with no budget constraint, and the same two
cases for lasso). This translates the *numeric core* at a handful of
representative penalty strengths spanning each sweep's range, rather than
the full 250-point plot -- same "numeric core, plot dropped" convention
used throughout this port, and the same
`solve_qp(..., ridge_penalty=...)`/`solve_qp(..., lasso_penalty=...)`
machinery test_lasso1.py/test_lasso3.py already establish for
quadprog_ridge/quadprog_lasso (see optim/quadprog.py's docstring for why
the variable-splitting original collapses into these two keyword
arguments).

Not cross-verified against Octave: the original's raw
`quadprog_ridge`/`quadprog_lasso` variable-splitting formulation turned
out to be numerically fragile for the *lasso* branches specifically under
Octave's free `quadprog` (returns an infeasible/degenerate all-zero
result, retcode=-3, even for the small, well-conditioned inputs used
here), which made an apples-to-apples comparison unreliable in this
environment; `solve_qp`'s native (variable-splitting-free) L1 handling via
cvxpy does not have that problem and returns sensible, well-conditioned
weights throughout, consistent with test_lasso1.py's already-established
output."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.optim.quadprog import solve_qp
from quanttoolbox.stats.moments import corr_to_cov

mu = np.array([0.05, 0.06, 0.08, 0.06])
sigma = np.array([0.15, 0.20, 0.25, 0.30])
rho = xpnd(np.array([1.00, 0.10, 1.00, 0.40, 0.70, 1.00, 0.50, 0.40, 0.80, 1.00]), method=1)
cov_matrix = corr_to_cov(sigma, rho)
n = 4
x0 = np.full(n, 1 / n)
gamma_x = 0.5
a_eq, b_eq = np.ones((1, n)), np.array([1.0])

lambda_ridge_grid = np.array([0.0, 0.1, 0.2, 0.3, 0.4, 0.5])
lambda_lasso_grid = np.array([0.0, 0.5, 1.0, 1.5, 2.0, 2.5]) / 100

print("Ridge (Static, toward zero, budget constrained):")
for lam in lambda_ridge_grid:
    x = solve_qp(
        cov_matrix,
        gamma_x * mu,
        a_eq=a_eq,
        b_eq=b_eq,
        lb=-100.0,
        ub=100.0,
        ridge_penalty=(lam * np.eye(n), np.zeros(n)),
    )
    print(f"  lambda={lam:.2f}  w={np.round(100 * x, 2)}")

print("Ridge (Dynamic, toward equal-weight, no budget constraint):")
for lam in lambda_ridge_grid:
    x = solve_qp(cov_matrix, gamma_x * mu, lb=-100.0, ub=100.0, ridge_penalty=(lam * np.eye(n), x0))
    print(f"  lambda={lam:.2f}  w={np.round(100 * x, 2)}")

print("Lasso (Static, toward zero, budget constrained):")
for lam in lambda_lasso_grid:
    x = solve_qp(
        cov_matrix,
        gamma_x * mu,
        a_eq=a_eq,
        b_eq=b_eq,
        lb=-100.0,
        ub=100.0,
        lasso_penalty=(lam * np.ones(n), np.zeros(n)),
    )
    print(f"  lambda={100 * lam:.2f}  w={np.round(100 * x, 2)}")

print("Lasso (Dynamic, toward equal-weight, no budget constraint):")
for lam in lambda_lasso_grid:
    x = solve_qp(
        cov_matrix, gamma_x * mu, lb=-100.0, ub=100.0, lasso_penalty=(lam * np.ones(n), x0)
    )
    print(f"  lambda={100 * lam:.2f}  w={np.round(100 * x, 2)}")

Volatility-target mean-variance under three weight-bound configurations — rpb/test_mvo3.py

"""Translated from Examples/rpb/test_mvo3.m -- Roncalli [2013], "Introduction
to Risk Parity and Budgeting", Example 1 (page 10): the same sigma-problem
target-matching as test_mvo2.py's section 3, now under three different
weight-bound configurations (unconstrained-ish, long-only, long-only capped
at 40% per asset), showing how the achievable target range narrows as
constraints tighten.

**A genuine bug in the original**, found while cross-checking this example
against Octave: `test_mvo3.m` (unlike `test_mvo2.m` and every other example
in this cluster) never calls `init_global`, which is what sets the global
`BISECTION_Tol` that `bisection.m`'s convergence loop depends on. With
`BISECTION_Tol` undefined, `while max(abs(a-b)) > BISECTION_Tol` compares
against an empty value, which both MATLAB and Octave treat as false --
so the loop runs zero iterations and `bisection` silently returns the raw
bracket midpoint `(0+10)/2 = 5.0` for every target, regardless of what the
target actually is. `compute_mvo_portfolio` then reports this as a
successful solve (retcode=1) even though the resulting portfolio doesn't
come close to the requested volatility target. Running `test_mvo3.m`
standalone in a fresh MATLAB/Octave session reproduces this; running it
right after `test_mvo2.m` in the same session "accidentally" works, because
`clear` (unlike `clear all`) doesn't clear globals, so `test_mvo2.m`'s
earlier `init_global` call leaves `BISECTION_Tol` set. See
`docs/matlab_bugs_found.md` for the full writeup. This translation produces
the *correct* target-matching numbers (as `test_mvo3.m` would if it called
`init_global` like its neighbors do), since `mvo_target_portfolio`'s Python
`bisection` always has a well-defined tolerance."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.portfolio.mean_variance import mvo_target_portfolio
from quanttoolbox.stats.moments import corr_to_cov

mu = np.array([0.05, 0.06, 0.08, 0.06])
sigma = np.array([0.15, 0.20, 0.25, 0.30])
rho = xpnd(np.array([1.00, 0.10, 1.00, 0.40, 0.70, 1.00, 0.50, 0.40, 0.80, 1.00]), method=1)
cov_matrix = corr_to_cov(sigma, rho)

sigma_targets = np.array([15.00, 20.00]) / 100

bound_configs = {
    "x1 (lb=-100, ub=100)": (-100.0, 100.0),
    "x2 (lb=0, ub=100, long-only)": (0.0, 100.0),
    "x3 (lb=0, ub=40, long-only, capped)": (0.0, 0.40),
}

for label, (lb, ub) in bound_configs.items():
    print(f"{label}:")
    results = mvo_target_portfolio(mu, cov_matrix, sigma_targets, problem="sigma", lb=lb, ub=ub)
    for target, r in zip(sigma_targets, results, strict=False):
        print(
            f"  target_sigma={100 * target:5.2f}  gamma={r.gamma:6.3f}  mu={100 * r.expected_return:6.2f}  "
            f"sigma={100 * r.volatility:6.2f}  w={np.round(100 * r.weights, 2)}"
        )

portfolio.erc_mdp

Thin re-export convenience module — see risk_budgeting (ERC) and mean_variance (MDP) for the actual implementations and alternatives discussion.

quanttoolbox.portfolio.erc_mdp

Convenience re-exports for Equal Risk Contribution and Most Diversified Portfolio construction.

Ported from QuantToolBox/rpb/compute_erc_portfolio.m, QuantToolBox/mloapa/compute_{ERC_ADMM,ERC_CCD,MDP_ADMM}.m.

These are both already implemented in full elsewhere in this package -- ERC in portfolio.risk_budgeting (it's the b=1/n special case of risk budgeting) and MDP in portfolio.mean_variance (it needs its own nonlinear solve, unrelated to risk budgeting's machinery). This module just re-exports both under the names matching the original MATLAB toolbox's dedicated rpb/mloapa entry points, so callers used to those names don't need to know the underlying module split.

erc_portfolio(cov_matrix, x0=None, method='ccd')

Equal Risk Contribution portfolio: risk budgeting with all budgets equal (b = 1/n), pure volatility risk measure.

Original: rpb/compute_erc_portfolio.m (+ mloapa/compute_ERC_{ADMM,CCD}.m, equivalent alternative solvers for the same problem)

Source code in src/quanttoolbox/portfolio/risk_budgeting.py

def erc_portfolio(
    cov_matrix: np.ndarray, x0: np.ndarray | None = None, method: str = "ccd"
) -> RiskBudgetingResult:
    """Equal Risk Contribution portfolio: risk budgeting with all budgets
    equal (b = 1/n), pure volatility risk measure.

    Original: rpb/compute_erc_portfolio.m (+ mloapa/compute_ERC_{ADMM,CCD}.m,
    equivalent alternative solvers for the same problem)
    """
    n = np.asarray(cov_matrix).shape[0]
    return solve_unconstrained(
        cov_matrix, b=np.full(n, 1.0 / n), mu=0.0, c=1.0, x0=x0, method=method
    )

mdp_portfolio(cov_matrix, a_eq=None, b_eq=None, c_ineq=None, d_ineq=None, lb=None, ub=None, x0=None)

Most Diversified Portfolio: maximize the diversification ratio (weighted-average individual volatility / portfolio volatility), i.e. minimize log(portfolio vol) - log(weighted-average individual vol).

Original: rpb/{compute_mdp_portfolio,compute_mdp_objective_function}.m

Source code in src/quanttoolbox/portfolio/mean_variance.py

def mdp_portfolio(
    cov_matrix: np.ndarray,
    a_eq: np.ndarray | None = None,
    b_eq: np.ndarray | None = None,
    c_ineq: np.ndarray | None = None,
    d_ineq: np.ndarray | None = None,
    lb: np.ndarray | float | None = None,
    ub: np.ndarray | float | None = None,
    x0: np.ndarray | None = None,
) -> MDPResult:
    """Most Diversified Portfolio: maximize the diversification ratio
    (weighted-average individual volatility / portfolio volatility), i.e.
    minimize log(portfolio vol) - log(weighted-average individual vol).

    Original: rpb/{compute_mdp_portfolio,compute_mdp_objective_function}.m
    """
    cov_matrix = np.asarray(cov_matrix, dtype=float)
    n = cov_matrix.shape[0]
    sigma = np.sqrt(np.diag(cov_matrix))

    a_eq = np.ones((1, n)) if a_eq is None else np.asarray(a_eq, dtype=float)
    b_eq = np.array([1.0]) if b_eq is None else np.asarray(b_eq, dtype=float)
    lb_arr = (
        np.full(n, -100.0)
        if lb is None
        else np.full(n, lb)
        if np.isscalar(lb)
        else np.asarray(lb, dtype=float)
    )
    ub_arr = (
        np.full(n, 100.0)
        if ub is None
        else np.full(n, ub)
        if np.isscalar(ub)
        else np.asarray(ub, dtype=float)
    )
    x0_arr = np.full(n, 1.0 / n) if x0 is None else np.asarray(x0, dtype=float)

    def objective(x: np.ndarray) -> float:
        port_vol = np.sqrt(x @ cov_matrix @ x)
        weighted_avg_vol = x @ sigma
        return float(np.log(port_vol) - np.log(weighted_avg_vol))

    constraints = [LinearConstraint(a_eq, b_eq, b_eq)]
    if c_ineq is not None:
        constraints.append(
            LinearConstraint(
                np.asarray(c_ineq, dtype=float), -np.inf, np.asarray(d_ineq, dtype=float)
            )
        )

    result = minimize(
        objective,
        x0_arr,
        method="SLSQP",
        bounds=list(zip(lb_arr, ub_arr, strict=True)),
        constraints=constraints,
        options={"maxiter": 1000, "ftol": 1e-12},
    )

    x = result.x
    x = np.where(np.abs(x) < 1e-10, 0.0, x)
    sigma_x = float(np.sqrt(x @ cov_matrix @ x))
    dr_x = float((x @ sigma) / sigma_x)

    return MDPResult(
        weights=x, volatility=sigma_x, diversification_ratio=dr_x, converged=result.success
    )