This notebook is based off a tutorial given at MIT on Feb 20th, 2015. That notebook is also in this folder, and covers more content than this notebook, including ellipsoidal uncertainty sets and cutting-edge data-driven sets.
Here we use Interact.jl to provide interactivity, and Gadfly.jl for visualizations. These features will not work if you looking at a static copy of this notebook (via, e.g. nbviewer.ipython.org - you must be running the notebook.
We're going to be looking at a fairly simple problem today, portfolio optimization. The general goal is to select a portfolio of assets to maximize return, but we also take into account the "risk" of the portfolio. In general, the higher the return, the higher the risk. Our goal is to find portfolios that trade off return and risk in such a way as to make sure that, for a given risk appetite/tolerance, we find the highest return.
In [ ]:
# You will need some packages
using Distributions # For data generation
using JuMP # a modeling language for optimization
using JuMPeR # adds robust optimization support to JuMP
using Interact # add interactivity to IJulia
using Gadfly # visualization
# You will also need a solver - any LP solver should
# work. See juliaopt.org for a list of available solvers.
# I'll be using Gurobi
using Gurobi
# Packages take quite a while to load in Julia 0.3
# Gadfly, in particular, takes many seconds
# Fortunately, this should be a one-time penalty!
In [ ]:
# We'll define some helper functions to generate data
# and to simulate the performance of our portfolios.
# These returns are normally distributed, which isn't
# very ideal/at all realistic, but it will do.
function generate_data(n, num_obs)
data = zeros(num_obs, n)
# Linking factors
beta = [(i-1.0)/n for i = 1:n]
for obs_ind = 1:num_obs
# Common market factor, mean 3%, sd 2%, truncate at +- 3 sd
z = rand(Normal(0.03, 0.02))
z = max(z, 0.03 - 3*0.02)
z = min(z, 0.03 + 3*0.02)
for asset_ind = 1:n
# Idiosyncratic contribution, mean 0%, sd 3%, truncated at +- 3 sd
asset = rand(Normal(0.00, 0.03))
asset = max(asset, 0.00 - 3*0.03)
asset = min(asset, 0.00 + 3*0.03)
data[obs_ind, asset_ind] = beta[asset_ind] * z + asset
end
end
return 100*data
end
# We'll also write a function to simulate the
# performance of a portfolio, and display the results
# graphically using Gadfly
function eval_portfolio(returns, x)
port_returns = returns * x
μ = round(mean(port_returns),2)
σ = round(std(port_returns),2)
hstack(
# Return density plot
plot(x=port_returns,
Geom.density,
Scale.x_continuous(minvalue=-10,maxvalue=15),
Scale.y_continuous(minvalue=0.00,maxvalue=0.25),
xintercept=[mean(port_returns)],Geom.vline,
layer(x=[-10.0],y=[0.20,0.15],label=["μ=$μ","σ=$σ"],
Geom.label(position=:right))
),
# Portfolio bar plot
plot(x=1:length(x), y=x,
Geom.bar,
Scale.y_continuous(minvalue=0.00,maxvalue=1.00))
)
end
# Lets generate some data that we'll use throughout
srand(2000)
# 5 assets, 10000 observations
returns = generate_data(5, 10000);
To get warmed up, we'll model the classic Markowitz portfolio optmization problem. Given $n$ assets, with average returns $\mathbf{r}$ and covariance matrix of returns $\mathbf{\Sigma}$, we can express the problem of finding the minimum risk portfolio for a given return as $$ \begin{alignat}{2} \min_{\mathbf{x}\geq 0} \quad & \mathbf{x}^T \mathbf{\Sigma} \mathbf{x} \\ \text{subject to} \quad & \mathbf{e}^T \cdot \mathbf{x} = 1 \\ & \mathbf{r}^T \cdot \mathbf{x} \geq R \end{alignat} $$ You can imagine a flipped problem where we maximize return subject to a risk tolerance. We also note that we've defined risk to be the variance of the resulting portfolio - this isn't probably a good idea in practice, but this is an older-style model.
In [ ]:
function solve_det_port(returns, min_return)
n = size(returns,2)
# Create a JuMP Model
# If you aren't using Gurobi, you'll need to change this
detport = Model(solver=GurobiSolver(OutputFlag=0))
# The allocation
@defVar(detport, 0 <= x[1:n] <= 1)
# Must be a valid allocation
@addConstraint(detport, sum{x[i], i=1:n} == 1)
# Must make min return
r = mean(returns,1)
@addConstraint(detport, sum{r[i]*x[i], i=1:n} >= min_return)
# Objective: minimize risk (variance)
Σ = cov(returns)
@setObjective(detport, Min, sum{Σ[i,j]*x[i]*x[j], i=1:n, j=1:n})
# Solve it!
solve(detport)
# Return the portfolio
getValue(x)[:]
end
# We'll let the user modify two parameters:
# Unoptimized portfolio: take a convex combination of an
# equal-share portfolio and a maximal risk portfolio.
# Minimum return of optimized portfolio
# Fun thing to try: adjust sliders until mean return
# is the same, then compare standard deviation and the
# allocations.
@manipulate for λ in 0.0:0.1:1.0, minR in 0.0:0.1:2.4
# Construct mixed portfolio
mix_port = λ*0.2*ones(5) + (1-λ)*[0,0,0,0,1]
# Optimized portfolio
opt_port = solve_det_port(returns, minR)
# Plot
vstack(
eval_portfolio(returns, mix_port),
eval_portfolio(returns, opt_port)
)
end
# This will take a while the first time. This time is
# the compilation for Gadfly - precompilation is really
# needed, and is coming soon(TM)
In a robust portfolio optmization problem, we still have our $n$ assets and our goal is to balance risk and return. Instead of using the mean and variance of the portfolio, we'll instead try to maximize the worst-case return over our uncertainty set, i.e. $$ \begin{alignat}{2} \min_{\mathbf{x}\geq 0} \quad & \min_{r \in U} \mathbf{r}^T \cdot \mathbf{x} \\ \text{subject to} \quad & \mathbf{e}^T \cdot \mathbf{x} = 1 \end{alignat} $$ The tradeoff between risk and return is baked into the uncertainty set $U$.
Polyhedral uncertainty sets are pretty popular in the literature because they don't turn (MI)LPs into (MI)SOCPs, even though an ellipsoidal set would seem to be more appropriate. The set we will do is inspired by a definition of a multivariate normal distribution. Given average returns $\bar{\mathbf{r}}$, a covariance matrix of returns $\Sigma$, we can say (assuming returns are normal) $$ \begin{alignat}{2} \mathbf{r} & = \Sigma^\frac{1}{2} \mathbf{z} + \bar{\mathbf{r}} \end{alignat} $$ We'll use that to inspire the uncertainty set: $$ \begin{alignat}{2} U & = \left\{ \mathbf{r} = \Sigma^{0.5} \mathbf{z} + \bar{\mathbf{r}}, \quad \| \mathbf{z} \|_1 \leq \Gamma \right\} \end{alignat} $$ where $\Gamma$ controls the conservativeness.
In [ ]:
function solve_poly_port(returns, Γ)
n = size(returns,2)
polyport = RobustModel(solver=GurobiSolver(OutputFlag=0))
@defVar(polyport, 0 <= x[1:n] <= 1)
@addConstraint(polyport, sum{x[i], i=1:n} == 1)
# Uncertainty
# Calculate the data we need
avg_r = mean(returns,1)
Σ = cov(returns)
Σhalf = chol(Σ)
# Create our Uncertains
# The returns
@defUnc(polyport, r[i=1:n])
# The factors
@defUnc(polyport, z[i=1:n])
# The absolute value of the factors
@defUnc(polyport, zabs[i=1:n] >= 0)
# Link all these together
# Constraints that are JUST on Uncertains are added
# to the uncertainty set automatically, not to the problem
for i in 1:n
@addConstraint(polyport, r[i] == sum{Σhalf[i,j]*z[j], j=1:n} + avg_r[i])
@addConstraint(polyport, zabs[i] >= z[i])
@addConstraint(polyport, zabs[i] >= -z[i])
end
@addConstraint(polyport, sum{zabs[i], i=1:n} <= Γ)
# Objective - same as before
@defVar(polyport, obj <= 10)
@setObjective(polyport, Max, obj)
@addConstraint(polyport, obj <= sum{r[i]*x[i], i=1:n})
# Print it to have a look
#print(polyport)
# Solve it!
solve(polyport)
#solve(polyport, prefer_cuts=true)
#solve(boxport, debug_printfinal=true)
# Return the portfolio
getValue(x)[:]
end
# We'll let the user modify two parameters:
# Unoptimized portfolio: take a convex combination of an
# equal-share portfolio and a maximal risk portfolio.
# Γ in the robust model
# Fun thing to try: adjust sliders until mean return
# is the same, then compare standard deviation and the
# allocations.
# e.g. λ=0.3, Γ=0.5
@manipulate for λ in 0.0:0.1:1.0, Γ in 0:0.25:5
# Construct mixed portfolio
mix_port = λ*0.2*ones(5) + (1-λ)*[0,0,0,0,1]
# Optimized portfolio
opt_port = solve_poly_port(returns, Γ)
# Plot
vstack(
eval_portfolio(returns, mix_port),
eval_portfolio(returns, opt_port)
)
end
# Will take some time the first time - thats JuMPeR
# compiling, again just a one of cost.
# Polyhedral uncertainty sets aren't very good for
# this use case, it seems!