The source files for all examples can be found in /examples.

Portfolio Optimisation

We consider a single-period Markowitz portfolio optimisation example. Assume that we have a portfolio with $n$ assets at the beginning of time period $t$. Given some forecasts on risks and expected returns we try to find the optimal trade vector that rebalances the portfolio to achieve a good balance between expected risk (variance) $x^\top \Sigma x$ and returns $\mu^\top x$. In it's most simple form we want to solve:

\[\begin{array}{ll} \text{maximize} & \mu^\top x - \gamma (x^\top \Sigma x)\\ \text{subject to} & 1^\top x = d + 1^\top x^0 \\ & x \geq 0, \end{array}\]

with variable $x \in \mathbf{R}^n$, $\mu$ forecasted (expected) returns, $\gamma > 0 $ risk aversion parameter. $x^0_i$ represents the initial investment in asset $i$ and $d$ represents the cash reserve. Consequently, the equality constraint tells us that the sum of the new allocation vector $x$ has to equal the initial allocation plus the cash reserve. Furthermore, the covariance matrix of our risk model is given by $\Sigma \in \mathbf{S}_+^n$. Here we assume a factor risk model, i.e. we can write $\Sigma$ as $\Sigma = D + F F^\top$ where $D$ is diagonal and the factor matrix $F$ has a lower rank $k < n$. This approach allows us to reduce the number of nonzeros in the problem. Furthermore, note that we don't consider shortselling in this example. Let's generate some problem data:

using LinearAlgebra, SparseArrays, Random, COSMO, JuMP, Test

# generate the data
rng = Random.MersenneTwister(1)
k = 5; # number of factors
n = k * 10; # number of assets
D = spdiagm(0 => rand(rng, n) .* sqrt(k))
F = sprandn(rng, n, k, 0.5); # factor loading matrix
μ = (3 .+ 9. * rand(rng, n)) / 100. # expected returns between 3% - 12%
γ = 1.0; # risk aversion parameter
d = 1 # we are starting from all cash
x0 = zeros(n);

We can now write the problem as a QP:

\[\begin{array}{ll} \text{minimize} & x^\top D x + y^\top y - \gamma^{-1} \mu^\top x \\ \text{subject to} & y = F^\top x \\ & 1^\top x = d + 1^\top x^0 \\ & x \geq 0. \end{array}\]

Before considering other effects, let's create the model in JuMP and solve it using COSMO:

model = JuMP.Model(COSMO.Optimizer);
@variable(model, x[1:n]);
@variable(model, y[1:k]);
@objective(model, Min,  x' * D * x + y' * y - 1/γ * μ' * x);
@constraint(model, y .== F' * x);
@constraint(model, sum(x) == d + sum(x0));
@constraint(model, x .>= 0);
JuMP.optimize!(model)

------------------------------------------------------------------
          COSMO v0.8.9 - A Quadratic Objective Conic Solver
                         Michael Garstka
                University of Oxford, 2017 - 2022
------------------------------------------------------------------

Problem:  x ∈ R^{55},
          constraints: A ∈ R^{56x55} (228 nnz),
          matrix size to factor: 111x111,
          Floating-point precision: Float64
Sets:     Nonnegatives of dim: 50
          ZeroSet of dim: 6
Settings: ϵ_abs = 1.0e-05, ϵ_rel = 1.0e-05,
          ϵ_prim_inf = 1.0e-04, ϵ_dual_inf = 1.0e-04,
          ρ = 0.1, σ = 1e-06, α = 1.6,
          max_iter = 5000,
          scaling iter = 10 (on),
          check termination every 25 iter,
          check infeasibility every 40 iter,
          KKT system solver: QDLDL
Acc:      Anderson Type2{QRDecomp},
          Memory size = 15, RestartedMemory,
          Safeguarded: true, tol: 2.0
Setup Time: 0.13ms

Iter:	Objective:	Primal Res:	Dual Res:	Rho:
1	-1.0930e-01	6.0222e-01	6.5188e-02	1.0000e-01
25	-7.9519e-02	3.5102e-09	1.9787e-09	1.0000e-01

------------------------------------------------------------------
>>> Results
Status: Solved
Iterations: 25
Optimal objective: -0.07952
Runtime: 0.001s (0.73ms)

After solving the problem, we can calculate the expected return and risk $\sigma= \sqrt{x^{* \top} \Sigma x^*}$:

x_opt = JuMP.value.(x);
y_opt = JuMP.value.(y);
expected_return_basic = dot(μ, x_opt)

0.09960959937666913

expected_risk_basic = sqrt(dot(y_opt, y_opt))

0.013048755030479095

Using standard deviation in the model

It is pointed out in [1] that above problem formulation can lead to numerical problems, e.g. if $\Sigma$ is not strictly positive semidefinite. Another option is to formulate the risk constraint in terms of the standard deviation $\|M^\top x \|$ where $M M^\top = D + F F^\top$ and bound it using a second-order cone constraint:

\[\begin{array}{ll} \text{minimize} & - \mu^\top x \\ \text{subject to} & \|M^\top x\| \leq \gamma \\ & 1^\top x = d + 1^\top x^0 \\ & x \geq 0. \end{array}\]

Mt = [D.^0.5; F']
model = JuMP.Model(COSMO.Optimizer);
@variable(model, x[1:n]);
@objective(model, Min, - μ' * x);
@constraint(model,  [γ; Mt * x] in SecondOrderCone()); # ||M'x|| <= γ
@constraint(model, sum(x) == d + sum(x0));
@constraint(model, x .>= 0);
JuMP.optimize!(model)

------------------------------------------------------------------
          COSMO v0.8.9 - A Quadratic Objective Conic Solver
                         Michael Garstka
                University of Oxford, 2017 - 2022
------------------------------------------------------------------

Problem:  x ∈ R^{50},
          constraints: A ∈ R^{107x50} (273 nnz),
          matrix size to factor: 157x157,
          Floating-point precision: Float64
Sets:     SecondOrderCone of dim: 56
          Nonnegatives of dim: 50
          ZeroSet of dim: 1
Settings: ϵ_abs = 1.0e-05, ϵ_rel = 1.0e-05,
          ϵ_prim_inf = 1.0e-04, ϵ_dual_inf = 1.0e-04,
          ρ = 0.1, σ = 1e-06, α = 1.6,
          max_iter = 5000,
          scaling iter = 10 (on),
          check termination every 25 iter,
          check infeasibility every 40 iter,
          KKT system solver: QDLDL
Acc:      Anderson Type2{QRDecomp},
          Memory size = 15, RestartedMemory,
          Safeguarded: true, tol: 2.0
Setup Time: 0.14ms

Iter:	Objective:	Primal Res:	Dual Res:	Rho:
1	-2.9234e-01	6.0111e-01	6.3843e-02	1.0000e-01
25	-1.1531e-01	6.3586e-03	2.8492e-03	1.0000e-01
50	-1.2033e-01	1.9181e-03	1.2572e-03	1.0000e-01
75	-1.1842e-01	2.7093e-04	1.1518e-03	1.0000e-01
100	-1.1871e-01	4.3383e-04	1.2079e-03	1.1116e-02
125	-1.1965e-01	3.5900e-02	1.3409e-03	1.1116e-02
150	-1.1953e-01	4.1597e-03	1.2176e-03	1.1116e-02
175	-1.1930e-01	6.3593e-02	1.5866e-03	1.1116e-02
200	-1.1842e-01	2.4249e-02	9.1097e-04	1.1116e-02
225	-1.1438e-01	3.1034e-02	1.0748e-03	1.1116e-02
250	-1.1795e-01	1.2000e-03	1.1944e-04	1.1116e-02
275	-1.1913e-01	2.6930e-04	1.8974e-05	1.1116e-02
300	-1.1882e-01	6.7647e-05	1.3779e-06	1.1116e-02
325	-1.1887e-01	7.0063e-06	4.1292e-07	1.1116e-02

------------------------------------------------------------------
>>> Results
Status: Solved
Iterations: 326 (incl. 1 safeguarding iter)
Optimal objective: -0.1189
Runtime: 0.005s (5.47ms)

Note that the result is different from the example above because $\gamma$ scales the problem in a different way. Here it can be seen as an upper bound on the standard deviation of the portfolio.

x_opt = JuMP.value.(x);
expected_return = dot(μ, x_opt)

0.11886638687579054

Let us verify that the bound holds:

@test norm(Mt * x_opt) <= γ

Test Passed

Pareto-optimal front

The above portfolio optimisation approach yields the optimal expected return for a given level of risk. The result is obviously impacted by the risk aversion $\gamma$ parameter. To visualise the trade-off and present the investor with an efficient Pareto optimal portfolio for their risk appetite we can compute the optimal portfolio for many choices of $\gamma$ and plot the corresponding risk-return trade-off curve.

gammas = [ 0.001, 0.01, 0.1,  0.5,  1., 3., 10, 100, 1000]
risks = zeros(length(gammas))
returns = zeros(length(gammas))
model = JuMP.Model(optimizer_with_attributes(COSMO.Optimizer, "verbose" => false));
@variable(model, x[1:n]);
@variable(model, y[1:k]);
@objective(model, Min,  x' * D * x + y' * y - 1/γ * μ' * x);
@constraint(model, y .== F' * x);
@constraint(model, sum(x) == d + sum(x0));
@constraint(model, x .>= 0);

# solve the same problem for different values of γ
for (k, gamma) in enumerate(gammas)
    coeff = - 1/gamma * μ
    JuMP.set_objective_coefficient.(model, x, coeff)
    JuMP.optimize!(model)
    local x_opt = JuMP.value.(x);
    local y_opt = JuMP.value.(y);
    returns[k] = dot(μ, x_opt)
    risks[k] = sqrt(dot(y_opt, y_opt))
end

We can now plot the risk-return trade-off curve:

using Plots
Plots.plot(risks, returns, xlabel = "Standard deviation (risk)", ylabel = "Expected return", title = "Risk-return trade-off for efficient portfolios", legend = false)

Transaction costs

In the model above we assume that trading the assets is free and does not impact the market. However, this is clearly not the case in reality. To make the example more realistic consider the following cost $c_j$ associated with the trade $δ_j = x_j - x_j^0$:

\[c_j(\delta_j) = a_j |\delta_j| + b_j |\delta_j|^{3/2},\]

where the first term models the bid-ask spread and broker fees for asset $j$. The second term models the impact on the market that our trade has. This is obviously only a factor if the volume of our trade is significant. The constant $b_j$ is a function of the total volume traded in the considered time periode and the price volatility of the asset and has to be estimated by the trader. To make this example simple we consider the same coefficients $a$ and $b$ for every asset. The $|\delta_j|^{3/2}$ term can be easily modeled using a power cone constraint $\mathcal{K}_{pow} = \{(x, y, z) \mid x^\alpha y^{(1-\alpha)} \geq |z|, x \geq 0, y \geq 0, 0 \leq \alpha \leq 1 \}$. In fact this can be used to model any market impact function with exponent greater than 1. We can write the total transaction cost $a^\top s + b^\top t$ where $s_j$ bounds the absolute value of $\delta_j$ and $t_{j}$ is used to bound the term $|x_j - x_j^0|^{3/2} \leq t_{j}$ using a power cone formulation: $(t_{j}, 1, x_j - x_j^0) \in \mathcal{K}_{pow}(2/3)$.

a = 1e-3
b = 1e-1
γ = 1.0;
model = JuMP.Model(optimizer_with_attributes(COSMO.Optimizer, "eps_abs" => 1e-5, "eps_rel" => 1e-5));
@variable(model, x[1:n]);
@variable(model, y[1:k]);
@variable(model, t[1:n]);
@variable(model, s[1:n]);
@objective(model, Min, x' * D * x + y' * y - 1/γ * μ' * x);
@constraint(model, y .== F' * x);
@constraint(model, x .>= 0);

# transaction costs
@constraint(model, sum(x) + a * sum(s) + b * sum(t) == d + sum(x0) );
@constraint(model, [i = 1:n], x[i] - x0[i] <= s[i]); # model the absolute value with slack variable s
@constraint(model, [i = 1:n], x0[i] - x[i] <= s[i]);
@constraint(model, [i = 1:n], [t[i], 1, x[i] - x0[i]] in MOI.PowerCone(2/3));
JuMP.optimize!(model)

------------------------------------------------------------------
          COSMO v0.8.9 - A Quadratic Objective Conic Solver
                         Michael Garstka
                University of Oxford, 2017 - 2022
------------------------------------------------------------------

Problem:  x ∈ R^{155},
          constraints: A ∈ R^{306x155} (628 nnz),
          matrix size to factor: 461x461,
          Floating-point precision: Float64
Sets:     Box of dim: 100
          Nonnegatives of dim: 50
          ZeroSet of dim: 6
          PowerCone of dim: 3
          PowerCone of dim: 3
          ... and 48 more
Settings: ϵ_abs = 1.0e-05, ϵ_rel = 1.0e-05,
          ϵ_prim_inf = 1.0e-04, ϵ_dual_inf = 1.0e-04,
          ρ = 0.1, σ = 1e-06, α = 1.6,
          max_iter = 5000,
          scaling iter = 10 (on),
          check termination every 25 iter,
          check infeasibility every 40 iter,
          KKT system solver: QDLDL
Acc:      Anderson Type2{QRDecomp},
          Memory size = 15, RestartedMemory,
          Safeguarded: true, tol: 2.0
Setup Time: 3.32ms

Iter:	Objective:	Primal Res:	Dual Res:	Rho:
1	-1.5637e-01	1.1398e+00	6.3250e-02	1.0000e-01
25	-7.7854e-02	2.0863e-03	5.8259e-05	1.0000e-01
50	-7.7862e-02	4.1404e-05	3.8814e-07	1.0000e-01
75	-7.7862e-02	1.7501e-06	3.3927e-08	1.0000e-01

------------------------------------------------------------------
>>> Results
Status: Solved
Iterations: 80 (incl. 5 safeguarding iter)
Optimal objective: -0.07786
Runtime: 0.083s (83.12ms)

Let's look at the expected return and the total transaction cost:

x_opt = JuMP.value.(x);
y_opt = JuMP.value.(y);
s_opt = JuMP.value.(s);
t_opt = JuMP.value.(t);
expected_return = dot(μ, x_opt)

0.09665905111704713

expected_risk = dot(y_opt, y_opt)

0.00016794910539706204

transaction_cost = a * sum(s_opt) + b * sum( t_opt)

0.02732439918367407

References

[1] Mosek Case Studies

This page was generated using Literate.jl.