Formulating Sharpe Ratio Objective Function in Gurobi
Ongoing
Dear Gurobi Community,
I am currently working on a portfolio optimization problem using Gurobi and encountered a challenge in formulating the objective function to optimize the Sharpe Ratio more precisely this equation:
$$
\frac{\mathbb{E}_t\left[R_{i, t+1}\right]-R_{f, t}}{\sqrt{\operatorname{Var}_t\left[R_{i, t+1}\right]}}
$$
Here is a snippet of the code I have implemented so far for the minimum variance portfolio:
m = gp.Model('portfolio')
x = m.addMVar(len(stocks))
portfolio_risk = x @ sigma @ x
m.setObjective(portfolio_risk, GRB.MINIMIZE)
z = m.addVars(len(stocks), vtype=gp.GRB.BINARY, name="z")
m.addConstr(sum(z[i] for i in range(len(stocks))) == 15, "cardinality_constraint")
m.write('portfolio_selection_optimization.lp')
m.addConstr(sum(x[i] for i in range(len(stocks))) == 1, "budget_constraint") m.write('portfolio_selection_optimization.lp')
for i in range(len(stocks)):
m.addConstr(x[i] <= 0.13 * z[i], f"linking_constraint_{i}")
for i in range(len(stocks)):
m.addConstr(0.03 * z[i] <= x[i], f"lower_bound_constraint_{i}")
r = np.mean(delta)
m.addConstr(sum(delta[i] * x[i] for i in range(len(stocks))) >= r, "target_return_constraint")
m.write('portfolio_selection_optimization.lp')
m.optimize()
if m.status == GRB.OPTIMAL:
for i in range(len(stocks)):
print(f"x[{i}] = {x[i].x}")
I would greatly appreciate any guidance on how to modify this code to set the Sharpe Ratio as the objective function in Gurobi.
Thank you in advance for your assistance.
Best regards,
Marcell Frisch
-
Hi Marcell,
We have an optimod on this:
Cheers,
David0 -
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 -
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 -
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,
David0 -
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: floatAnd the efficient frontier can be built by varying the target return right?
Thank you in advance for your answer!
Cheers,
Marcell
0 -
Yes, that should work. Varying the target return should produce the efficient frontier.
0 -
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 -
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 -
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 memoryHow 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 -
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 -
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,
Raik0 -
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.
Comments
12 comments