Skip to main content

Formulating Sharpe Ratio Objective Function in Gurobi

Ongoing

Comments

12 comments

  • Eli Towle
    Gurobi Staff Gurobi Staff

    Constraints in your original model (expressed in terms of variables \( x \) representing portfolio weight) can be translated to the OptiMod model's variable space (reformulated to use variables \( y \)) using the relationship \( x_i = y_i / \sum_{j=1}^n y_j \) for all assets \( i \). For example, building off of the Sharpe ratio OptiMod model code:

    z = model.addMVar(mu.size, vtype=GRB.BINARY, name="z")

    # Variable ysum equal to sum of y variables
    ysum = model.addMVar(shape=(1,), name="ysum")
    model.addConstr(ysum == y.sum(), name="set_ysum")

    # Minimum/maximum holdings
    # lb * z <= x <= ub * z
    # lb * z <= y / ysum <= ub * z
    # lb * z * ysum <= y <= ub * z * ysum
    model.addConstr(0.03 * z * ysum <= y, name="minimum_holding")
    model.addConstr(y <= 0.43 * z * ysum, name="maximum_holding")

    # Cardinality constraint
    model.addConstr(z.sum() == 3, name="cardinality_constraint")

    # Target return
    # mu @ x >= r
    # mu @ (y / ysum) >= r
    # mu @ y >= r * ysum
    r = np.mean(mu)
    model.addConstr(mu @ y >= r * ysum, name="target_return_rate")
    1
  • David Torres Sanchez
    Gurobi Staff Gurobi Staff

    Hi Marcell,

    We have an optimod on this:

    Cheers, 
    David

    0
  • Marcell Frisch
    First Comment
    First Question

    Hi David, 

    The problem is I have 482 assets and I have these constraints :
    how many assets should be in the portfolio,
    how much should be the minimum and the maximum weight can be assigned to one asset,
    weights should add up to one, 
    the returns * weights should be higher or equal to the target return.

    the problem that in the optimod is just doing this on the whole dataset and there is no cardinality constraints but in my case, I need to do an optimization where I have to choose those assets that maximise the shape ratio while not breaking these constraints.

    Cheers, 

    Marcell

    0
  • David Torres Sanchez
    Gurobi Staff Gurobi Staff

    Thanks, I see the difference now.
    Can you not use the same reformulation trick (as in the drop down "Mathematical Model")?
    And re-write the constraints to use \(y\) variables instead?

    Cheers, 
    David

    0
  • Marcell Frisch
    First Comment
    First Question

    Dear David,

    So if I understand correctly I should just put this part:

    z = model.addMVar(mu.size, vtype=GRB.BINARY, name="z")
    
    # Variable ysum equal to sum of y variables
    ysum = model.addMVar(shape=(1,), name="ysum")
    model.addConstr(ysum == y.sum(), name="set_ysum")
    
    # Minimum/maximum holdings
    # lb * z <= x <= ub * z
    # lb * z <= y / ysum <= ub * z
    # lb * z * ysum <= y <= ub * z * ysum
    model.addConstr(0.03 * z * ysum <= y, name="minimum_holding")
    model.addConstr(y <= 0.43 * z * ysum, name="maximum_holding")
    
    # Cardinality constraint
    model.addConstr(z.sum() == 3, name="cardinality_constraint")
    
    # Target return
    # mu @ x >= r
    # mu @ (y / ysum) >= r
    # mu @ y >= r * ysum
    r = np.mean(mu)
    model.addConstr(mu @ y >= r * ysum, name="target_return_rate")

    into the function like this and it should work right?:


    import math
    from dataclasses import dataclass
    from typing import Union

    import gurobipy as gp
    import numpy as np
    import pandas as pd
    from gurobipy import GRB

    from gurobi_optimods.utils import optimod


    @optimod()
    def max_sharpe_ratio(cov_matrix, mu, rf_rate=0, *, create_env):
        """
        Solve the problem of finding a portfolio that maximizes the
        Sharpe ratio.

        Parameters
        ----------
        cov_matrix : ndarray or DataFrame
            2D positive-semidefinite variance-covariance matrix :math:`\Sigma`
        mu : ndarray or Series
            Expected return rates :math:`\mu`
        rf_rate : float >= 0, optional
            Non-negative risk-free rate of return (defaults to ``0``)

        Returns
        -------
        result : SharpeRatioResult
            A data class representing the portfolio that maximizes the Sharpe
            ratio:

            * ``result.x``: The relative portfolio allocations :math:`x`
            * ``result.sharpe_ratio``: The Sharpe ratio of the portfolio
            * ``result.ret``: The estimated return :math:`\mu^T x` of the portfolio
            * ``result.risk``: The estimated risk :math:`x^T \Sigma x` of the
              portfolio
        """
        indices = None

        if isinstance(cov_matrix, pd.DataFrame):
            indices = cov_matrix.index
            cov_matrix = cov_matrix.to_numpy()
        elif not isinstance(cov_matrix, np.ndarray):
            raise TypeError(f"Unknown covariance matrix type: {type(cov_matrix)}")

        if cov_matrix.ndim != 2:
            raise ValueError(
                f"Covariance matrix should be in 2 dimensions, not {cov_matrix.ndim}"
            )

        if isinstance(mu, pd.Series):
            if indices is None:
                indices = mu.index
            elif not mu.index.equals(indices):
                raise ValueError("Indices of cov_matrix and mu are misaligned")
            mu = mu.to_numpy()
        elif not isinstance(mu, np.ndarray):
            raise TypeError(f"Unknown return rates type: {type(mu)}")

        if mu.ndim != 1:
            raise ValueError(f"Return rates should be in 1 dimension, not {mu.ndim}")

        if not isinstance(rf_rate, float) and not isinstance(rf_rate, int):
            raise TypeError(f"Unknown risk-free return rate type: {type(rf_rate)}")
        elif rf_rate < 0:
            raise ValueError("Risk-free return rate should be non-negative")
        elif (mu < rf_rate).all():
            raise ValueError(
                f"No expected returns are greater than risk-free return rate of {rf_rate}"
            )

        result = _max_sharpe_ratio_numpy(cov_matrix, mu, rf_rate, create_env)

        if indices is not None:
            result.x = pd.Series(data=result.x, index=indices)

        return result


    def _max_sharpe_ratio_numpy(cov_matrix, mu, rf_rate, create_env):
        with create_env() as env, gp.Model("sharpe_ratio", env=env) as model:
            y = model.addMVar(mu.size, name="y")
            model.addConstr((mu - rf_rate) @ y == 1)
            model.setObjective(y @ cov_matrix @ y, sense=GRB.MINIMIZE)
            z = model.addMVar(mu.size, vtype=GRB.BINARY, name="z")

    # Variable ysum equal to sum of y variables
            ysum = model.addMVar(shape=(1,), name="ysum")
            model.addConstr(ysum == y.sum(), name="set_ysum")

    # Minimum/maximum holdings
    # lb * z <= x <= ub * z
    # lb * z <= y / ysum <= ub * z
    # lb * z * ysum <= y <= ub * z * ysum
            model.addConstr(0.03 * z * ysum <= y, name="minimum_holding")
            model.addConstr(y <= 0.13 * z * ysum, name="maximum_holding")

    # Cardinality constraint
            model.addConstr(z.sum() == 15, name="cardinality_constraint")

    # Target return
    # mu @ x >= r
    # mu @ (y / ysum) >= r
    # mu @ y >= r * ysum
            r = np.mean(mu)
            model.addConstr(mu @ y >= r * ysum, name="target_return_rate")
            
            model.optimize()

            # Translate solution to original variable space
            x = y.X / y.X.sum()
            ret = mu @ x
            risk = x @ cov_matrix @ x
            sharpe_ratio = (ret - rf_rate) / math.sqrt(risk)
            return SharpeRatioResult(x, sharpe_ratio, ret, risk)


    @dataclass
    class SharpeRatioResult:
        """
        Data class representing the portfolio that maximizes the Sharpe ratio.


        Attributes
        ----------
        x : ndarray or Series
            The relative portfolio allocations :math:`x`
        sharpe_ratio : float
            The Sharpe ratio of the portfolio
        ret : float
            The estimated return :math:`\mu^T x` of the portfolio
        risk : float
            The estimated risk :math:`x^T \Sigma x` of the portfolio
        """

        x: Union[np.ndarray, pd.Series]
        sharpe_ratio: float
        ret: float
      risk: float

    And the efficient frontier can be built by varying the target return right?

    Thank you in advance for your answer! 

    Cheers,

    Marcell

    0
  • Eli Towle
    Gurobi Staff Gurobi Staff

    Yes, that should work. Varying the target return should produce the efficient frontier.

    0
  • Raik Becker
    First Comment

    Hi!

    I have a similar problem and I am making use of the transformed representation. All assets have a minimum and maximum holding. In total I have three assets, say A1, A2 and A3. Now this is different in my case:

    * The asset A1 is special and can only be taken into the portfolio when A2 has passed a certain threshold. That is easy to implement with the integers z and ysum (as defined above).

    * A1 should not count into the weights (x in the original problem in OptiMods) constraint, while its mean and covariance should be considered as usual. So the sum of the weights of A2 and A3 needs to be one, while the sum of A1, A2 and A3 can be even larger (1.15 due to the maximum holding of A1 is 0.15). I struggle bringing this constraint into the transformed formulation. Could you help?

    0
  • Eli Towle
    Gurobi Staff Gurobi Staff

    My initial thought is that you don't need a special constraint. In the original problem, \( x \) represents what proportion of your budget is allocated to each asset. By definition, these numbers should sum to 1. Maybe another way of looking at it: if you invested $0.15 in A1, $0.40 in A2, and $0.60 in A3, this would be represented in the original model as \(x_1=0.13043\), \(x_2=0.34783\), and \(x_3=0.52174\).

    Do the results not make sense without a special constraint? It would be helpful to see a small example showing that.

    0
  • Marcell Frisch
    First Comment
    First Question

    Dear Eli Towle,

    The code works perfectly for my basic cases but when I try to run it for 843 assets it runs for a very long time and I get this error message:


    Explored 5308982 nodes (64204304 simplex iterations) in 36602.53 seconds (7673.94 work units)
    Thread count was 8 (of 8 available processors)

    Solution count 1: 19.9949 

    Solve interrupted (error code 10001)
    Best objective 1.999494801962e+01, best bound 1.935112130712e+01, gap 3.2199%
    Traceback (most recent call last):

      Cell In[8], line 1
        result_US_india_full=max_sharpe_ratio(sigma, mu)

      File ~\anaconda3\Lib\site-packages\gurobi_optimods\utils.py:137 in optimod_decorated
        return func(*args, create_env=create_env, **kwargs)

      Cell In[1], line 74 in max_sharpe_ratio
        result = _max_sharpe_ratio_numpy(cov_matrix, mu, rf_rate, create_env)

      Cell In[1], line 110 in _max_sharpe_ratio_numpy
      model.optimize()

    File src\\gurobipy\\model.pxi:893 in gurobipy.Model.optimize GurobiError: Out of memory

    How can it be solved for this problem? I saw parallelization could solve the problem but the problem is for me it is taking even more time in that case, or is it possible to just set the gap higher or how can I solve this problem?

    Thank you in advance for your answer and help!

    Best regards,

    Marcell Frisch

    0
  • Eli Towle
    Gurobi Staff Gurobi Staff

    Your machine is running out of memory. The article How do I avoid an out-of-memory condition? contains ideas for reducing memory usage.

    By default, Gurobi will not stop solving this model until it reaches a 0.01% optimality gap. Your log output shows that the optimality gap is 3.2199% after the machine runs out of memory 36602 seconds into the solve. You could set the MIPGap parameter to direct Gurobi to terminate at a larger relative optimality gap.

    0
  • Raik Becker
    First Comment

    Hi Eli Towle,
    I haven't had the time to work on this project but now I would like to get back to your answer. There is one case where my special asset A1 is causing issues.

    Here is the code I am running.

    import math
    from gurobipy import Model, GRB
    import pandas as pd
    # Inputs
    mu_dict = {
        "A1": 1000,
        "A2": 50,
        "A3": 1000,
        }

    mu = pd.Series(mu_dict).to_numpy()

    cov_matrix = pd.DataFrame(
        data={
            "A1": [1, 0, 0],
            "A2": [0, 1, 0],
            "A3": [0, 0, 1],
            },
    ).to_numpy()
    # Model
    model = Model("sharpe_ratio")

    y = model.addMVar(mu.size, name="y")
    z = model.addMVar(mu.size, vtype=GRB.BINARY, name="z")
    model.addConstr(mu @ y == 1)

    # Minimum/maximum holding
    ysum = model.addMVar(shape=(1,), name="ysum")
    model.addConstr(ysum == y.sum(), name="set_ysum")

    # A1
    model.addConstr(0.1 * z[0] * ysum <= y[0])
    model.addConstr(y[0] <= 0.3 * z[0] * ysum)
    # A2
    model.addConstr(0.4 * z[1] * ysum <= y[1])
    model.addConstr(y[1] <= 1 * z[1] * ysum)
    # A2
    model.addConstr(0.1 * z[2] * ysum <= y[2])
    model.addConstr(y[2] <= 0.3 * z[2] * ysum)

    # In order to invest into A1, A2 needs to have passed
    # a certain threshold (0.4) and the amount the A2 is
    # above that threshold can be put into A1.
    model.addConstr(0.4 * z[1] * ysum + y[0] <= y[1])

    # Objective function
    model.setObjective(
        y @ cov_matrix @ y,
        sense=GRB.MINIMIZE,
    )

    model.optimize()

    x = y.X / y.X.sum()
    ret = mu @ x
    risk = x @ cov_matrix @ x
    ratio = ret / math.sqrt(risk)
    print(x)

    This prints out [0.15 0.55 0.3]. So A1 gets only 15% because of the threshold constraint (see code) with A2. However, A2 should become 0.7, A3 0.3 (maximum holding) and A1 should get 0.3 as well. The total sum is then 1 = A2 + A3 - A1.

    Formulating this in the original formulation is easy but I I haven't find a way in the transformed version. Could you please have a look?


    Cheers,
    Raik

    0
  • Eli Towle
    Gurobi Staff Gurobi Staff

    Sorry for the delayed reply. I was out of office the last few weeks.

    The Background: Mathematical Model section of the Maximum Sharpe Ratio OptiMod documentation describes the original Sharpe ratio model:

    $$\begin{alignat*}{2}\max_x\ && \frac{\mu^\top x}{\sqrt{x^\top \Sigma x}} & \\ \textrm{s.t.}\ && \sum_{i=1}^n x_i ={} &1 \qquad \textbf{(1)} \\ && x \geq {} &0\end{alignat*}$$

    and its reformulation:

    $$\begin{alignat*}{2} \max_y\ && \frac{1}{\sqrt{y^\top \Sigma y}} & \\ \textrm{s.t.}\enspace && \mu^\top y = {} &1 \qquad \textbf{(2)} \\ && y \geq {} &0. \end{alignat*}$$

    The section describes how models \(\textbf{(1)}\) and \(\textbf{(2)}\) are equivalent. In particular, a solution \(\bar{x}\) of \(\textbf{(1)}\) can be mapped to a solution \(\bar{y}\) of \(\textbf{(2)}\) of equal objective value with the transformation \(\bar{y}_i := \bar{x}_i / \mu^\top \bar{x}\). Similarly, a solution \(\bar{y}\) of \(\textbf{(2)}\) can be mapped to a solution \(\bar{x}\) of \(\textbf{(1)}\) of equal objective value with the transformation \(\bar{x}_i := \bar{y}_i / \sum_{j=1}^n \bar{y}_j\).

    Your problem differs slightly from \(\textbf{(1)}\) because the first asset should not be included in the weight summation. That is, the weight summation should start at \(2\) instead of \(1\):

    $$\begin{alignat*}{2}\max_x\ && \frac{\mu^\top x}{\sqrt{x^\top \Sigma x}} & \\ \textrm{s.t.}\ && \sum_{i=2}^n x_i ={} &1 \qquad \textbf{(3)} \\ && x \geq {} &0.\end{alignat*}$$

    I believe we can show that your model \(\textbf{(3)}\) is also equivalent to the reformulated model \(\textbf{(2)}\). The transformation from \(\textbf{(3)}\) to \(\textbf{(2)}\) is the same as before:

    $$\begin{align*}\bar{y}_i := \frac{\bar{x}_i}{\mu^\top \bar{x}}.\end{align*}$$

    But the transformation from \(\textbf{(2)}\) to \(\textbf{(3)}\) is different, with the summation now ignoring \(\bar{y}_1\):

    $$\begin{align*}\bar{x}_i := \frac{\bar{y}_i}{\sum_{j=2}^n \bar{y}_j}.\end{align*}$$

    If this is correct, your code needs two small updates to reflect this new transformation. First, update the constraint defining the \(\texttt{ysum}\) variable to not include the first asset:

    model.addConstr(ysum == y[1:].sum(), name="set_ysum")

    Then, update the transformation back to the original variable space:

    x = y.X / ysum.X

    After these changes, the code prints the solution [0.3, 0.7, 0.3].

    It's probably worth mentioning that the Sharpe ratio is invariant to the amount of money you invest in a portfolio with a fixed proportion of assets. That is, a portfolio with

    $$\begin{align*}\textrm{A1} &= 0.3 \\ \textrm{A2} &= 0.7 \\ \textrm{A3} &= 0.3\end{align*}$$

    yields the same Sharpe ratio as a portfolio with these values normalized so they sum to 1:

    $$\begin{align*}\textrm{A1} &\approx 0.230769 \\ \textrm{A2} &\approx 0.538462 \\ \textrm{A3} &\approx 0.230769.\end{align*}$$

    0

Please sign in to leave a comment.