# Portfolio optimization

This tutorial was generated using Literate.jl. Download the source as a .jl file.

This tutorial was originally contributed by Arpit Bhatia.

Optimization models play an increasingly important role in financial decisions. Many computational finance problems can be solved efficiently using modern optimization techniques.

This tutorial solves the famous Markowitz Portfolio Optimization problem with data from lecture notes from a course taught at Georgia Tech by Shabbir Ahmed.

## Required packages

This tutorial uses the following packages

using JuMP
import DataFrames
import Ipopt
import MultiObjectiveAlgorithms as MOA
import Plots
import Statistics
import StatsPlots

## Formulation

Suppose we are considering investing 1000 dollars in three non-dividend paying stocks, IBM (IBM), Walmart (WMT), and Southern Electric (SEHI), for a one-month period.

We will use the initial money to buy shares of the three stocks at the current market prices, hold these for one month, and sell the shares off at the prevailing market prices at the end of the month.

As a rational investor, we hope to make some profit out of this endeavor, that is, the return on our investment should be positive.

Suppose we bought a stock at $p$ dollars per share in the beginning of the month, and sold it off at $s$ dollars per share at the end of the month. Then the one-month return on a share of the stock is $\frac{s-p}{p}$.

Since the stock prices are quite uncertain, so is the end-of-month return on our investment. Our goal is to invest in such a way that the expected end-of-month return is at least $50 or 5%. Furthermore, we want to make sure that the “risk” of not achieving our desired return is minimum. Note that we are solving the problem under the following assumptions: 1. We can trade any continuum of shares. 2. No short-selling is allowed. 3. There are no transaction costs. We model this problem by taking decision variables$x_{i}, i=1,2,3,$denoting the dollars invested in each of the 3 stocks. Let us denote by$\tilde{r}_{i}$the random variable corresponding to the monthly return (increase in the stock price) per dollar for stock$i$. Then, the return (or profit) on$x_{i}$dollars invested in stock$i$is$\tilde{r}_{i} x_{i},$and the total (random) return on our investment is$\sum_{i=1}^{3} \tilde{r}_{i} x_{i}.$The expected return on our investment is then$\mathbb{E}\left[\sum_{i=1}^{3} \tilde{r}_{i} x_{i}\right]=\sum_{i=1}^{3} \overline{r}_{i} x_{i},$where$\overline{r}_{i}$is the expected value of the$\tilde{r}_{i}.$Now we need to quantify the notion of “risk” in our investment. Markowitz, in his Nobel prize winning work, showed that a rational investor’s notion of minimizing risk can be closely approximated by minimizing the variance of the return of the investment portfolio. This variance is given by: $$$\operatorname{Var}\left[\sum_{i=1}^{3} \tilde{r}_{i} x_{i}\right] = \sum_{i=1}^{3} \sum_{j=1}^{3} x_{i} x_{j} \sigma_{i j}$$$ where$\sigma_{i j}$is the covariance of the return of stock$i$with stock$j$. Note that the right hand side of the equation is the most reduced form of the expression and we have not shown the intermediate steps involved in getting to this form. We can also write this equation as: $$$\operatorname{Var}\left[\sum_{i=1}^{3} \tilde{r}_{i} x_{i}\right] =x^\top Q x$$$ Where$Q$is the covariance matrix for the random vector$\tilde{r}. Finally, we can write the model as: \begin{aligned} \min x^\top Q x \\ \text { s.t. } \sum_{i=1}^{3} x_{i} \leq 1000 \\ \overline{r}^\top x \geq 50 \\ x \geq 0 \end{aligned} ## Data For the data in our problem, we use the stock prices given below, in monthly values from November 2000, through November 2001. df = DataFrames.DataFrame( [ 93.043 51.826 1.063 84.585 52.823 0.938 111.453 56.477 1.000 99.525 49.805 0.938 95.819 50.287 1.438 114.708 51.521 1.700 111.515 51.531 2.540 113.211 48.664 2.390 104.942 55.744 3.120 99.827 47.916 2.980 91.607 49.438 1.900 107.937 51.336 1.750 115.590 55.081 1.800 ], [:IBM, :WMT, :SEHI], ) 13×3 DataFrame RowIBMWMTSEHI Float64Float64Float64 193.04351.8261.063 284.58552.8230.938 3111.45356.4771.0 499.52549.8050.938 595.81950.2871.438 6114.70851.5211.7 7111.51551.5312.54 8113.21148.6642.39 9104.94255.7443.12 1099.82747.9162.98 1191.60749.4381.9 12107.93751.3361.75 13115.5955.0811.8 Next, we compute the percentage return for the stock in each month: returns = diff(Matrix(df); dims = 1) ./ Matrix(df[1:end-1, :]) 12×3 Matrix{Float64}: -0.0909042 0.0192374 -0.117592 0.317645 0.0691744 0.0660981 -0.107023 -0.118137 -0.062 -0.0372369 0.00967774 0.533049 0.197132 0.0245391 0.182197 -0.0278359 0.000194096 0.494118 0.0152087 -0.0556364 -0.0590551 -0.0730406 0.145487 0.305439 -0.0487412 -0.140428 -0.0448718 -0.0823425 0.0317639 -0.362416 0.178261 0.0383915 -0.0789474 0.0709025 0.0729508 0.0285714 The expected monthly return is: r = vec(Statistics.mean(returns; dims = 1)) 3-element Vector{Float64}: 0.026002150277777348 0.008101316405671459 0.07371590949198982 and the covariance matrix is: Q = Statistics.cov(returns) 3×3 Matrix{Float64}: 0.018641 0.00359853 0.00130976 0.00359853 0.00643694 0.00488727 0.00130976 0.00488727 0.0686828 ## JuMP formulation model = Model(Ipopt.Optimizer) set_silent(model) @variable(model, x[1:3] >= 0) @objective(model, Min, x' * Q * x) @constraint(model, sum(x) <= 1000) @constraint(model, r' * x >= 50) optimize!(model) solution_summary(model) * Solver : Ipopt * Status Result count : 1 Termination status : LOCALLY_SOLVED Message from the solver: "Solve_Succeeded" * Candidate solution (result #1) Primal status : FEASIBLE_POINT Dual status : FEASIBLE_POINT Objective value : 2.26344e+04 Dual objective value : 4.52688e+04 * Work counters Solve time (sec) : 8.84104e-03  The optimal allocation of our assets is: value.(x) 3-element Vector{Float64}: 497.045529849864 -9.670578501816873e-9 502.9544801594809 So we spend497 on IBM, and $503 on SEHI. This results in a variance of: scalar_variance = value(x' * Q * x) 22634.417849884147 and an expected return of: scalar_return = value(r' * x) 49.999999500002374 ## Multi-objective portfolio optimization The previous model returned a single solution that minimized the variance, ensuring that our expected return was at least$50. In practice, we might be willing to accept a slightly higher variance if it meant a much increased expected return. To explore this problem space, we can instead formulate our portfolio optimization problem with two objectives:

1. to minimize the variance
2. to maximize the expected return

The solution to this bi-objective problem is the efficient frontier of modern portfolio theory, and each point in the solution is a point with the best return for a fixed level of risk.

model = Model(() -> MOA.Optimizer(Ipopt.Optimizer))
set_silent(model)

We also need to choose a solution algorithm for MOA. For our problem, the efficient frontier will have an infinite number of solutions. Since we cannot find all of the solutions, we choose an approximation algorithm and limit the number of solution points that are returned:

set_optimizer_attribute(model, MOA.Algorithm(), MOA.EpsilonConstraint())
set_optimizer_attribute(model, MOA.SolutionLimit(), 50)

Now we can define the rest of the model:

@variable(model, x[1:3] >= 0)
@constraint(model, sum(x) <= 1000)
@expression(model, variance, x' * Q * x)
@expression(model, expected_return, r' * x)
# We want to minimize variance and maximize expected return, but we must pick
# a single objective sense Min, and negate any Max objectives:
@objective(model, Min, [variance, -expected_return])
optimize!(model)
solution_summary(model)
* Solver : MOA[algorithm=MultiObjectiveAlgorithms.EpsilonConstraint, optimizer=Ipopt]

* Status
Result count       : 50
Termination status : OPTIMAL
Message from the solver:
"Solve complete. Found 50 solution(s)"

* Candidate solution (result #1)
Primal status      : FEASIBLE_POINT
Dual status        : NO_SOLUTION
Objective value    : [2.58170e-08,-5.36777e-05]
Objective bound    : [5.78303e-09,-7.37159e+01]
Dual objective value : -7.37159e+01

* Work counters
Solve time (sec)   : 6.36807e-01


The algorithm found 50 different solutions. Let's plot them to see how they differ:

objective_space = Plots.hline(
[scalar_return];
label = "Single-objective solution",
linecolor = :red,
)
Plots.vline!(objective_space, [scalar_variance]; label = "", linecolor = :red)
Plots.scatter!(
objective_space,
[value(variance; result = i) for i in 1:result_count(model)],
[value(expected_return; result = i) for i in 1:result_count(model)];
xlabel = "Variance",
ylabel = "Expected Return",
label = "",
title = "Objective space",
markercolor = "white",
markersize = 5,
legend = :bottomright,
)
for i in 1:result_count(model)
y = objective_value(model; result = i)
Plots.annotate!(objective_space, y, -y, (i, 3))
end

decision_space = StatsPlots.groupedbar(
vcat([value.(x; result = i)' for i in 1:result_count(model)]...);
bar_position = :stack,
label = ["IBM" "WMT" "SEHI"],
xlabel = "Solution #",
ylabel = "Investment (\\$)",
title = "Decision space",
)
Plots.plot(objective_space, decision_space; layout = (2, 1), size = (600, 600))

Perhaps our trade-off wasn't so bad after all. Our original solution corresponded to picking a solution #17. If we buy more SEHI, we can increase the return, but the variance also increases. If we buy less SEHI, such as a solution like #5 or #6, then we can achieve the corresponding return without deploying all of our capital. We should also note that at no point should we buy WMT.