In this notebook we show how to use the Second Order Cone (SOC) constraint in the variance portfolio optimization problem.
The minimization of portfolio variance is a quadratic optimization problem that can be posed as:
$$ \begin{equation} \begin{aligned} & \underset{x}{\text{min}} & & x^{\tau} \Sigma x \\ & \text{s.t.} & & \mu x^{\tau} \geq \bar{\mu} \\ & & & \sum_{i=1}^{N} x_i = 1 \\ & & & x_i \geq 0 \; ; \; \forall \; i =1, \ldots, N \\ \end{aligned} \end{equation} $$Where $x$ are the weights of assets, $\mu$ is the mean vector of expected returns and $\bar{\mu}$ the minimum expected return of portfolio.
####################################
# Downloading Data
####################################
!pip install --quiet yfinance
import numpy as np
import pandas as pd
import yfinance as yf
import warnings
warnings.filterwarnings("ignore")
yf.pdr_override()
pd.options.display.float_format = '{:.4%}'.format
# Date range
start = '2016-01-01'
end = '2019-12-30'
# Tickers of assets
assets = ['JCI', 'TGT', 'CMCSA', 'CPB', 'MO', 'APA', 'MMC', 'JPM',
'ZION', 'PSA', 'BAX', 'BMY', 'LUV', 'PCAR', 'TXT', 'TMO',
'DE', 'MSFT', 'HPQ', 'SEE', 'VZ', 'CNP', 'NI', 'T', 'BA']
assets.sort()
# Downloading data
data = yf.download(assets, start = start, end = end)
data = data.loc[:,('Adj Close', slice(None))]
data.columns = assets
# Calculating returns
Y = data[assets].pct_change().dropna()
display(Y.head())
[*********************100%***********************] 25 of 25 completed
APA | BA | BAX | BMY | CMCSA | CNP | CPB | DE | HPQ | JCI | ... | NI | PCAR | PSA | SEE | T | TGT | TMO | TXT | VZ | ZION | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
Date | |||||||||||||||||||||
2016-01-05 | -2.0257% | 0.4057% | 0.4036% | 1.9693% | 0.0180% | 0.9305% | 0.3678% | 0.5783% | 0.9483% | -1.1953% | ... | 1.5881% | 0.0212% | 2.8236% | 0.9758% | 0.6987% | 1.7539% | -0.1730% | 0.2410% | 1.3735% | -1.0857% |
2016-01-06 | -11.4863% | -1.5878% | 0.2412% | -1.7557% | -0.7727% | -1.2473% | -0.1736% | -1.1239% | -3.5867% | -0.9551% | ... | 0.5547% | 0.0212% | 0.1592% | -1.5647% | -0.1466% | -1.0155% | -0.7653% | -3.0048% | -0.9034% | -2.9145% |
2016-01-07 | -5.1388% | -4.1922% | -1.6573% | -2.7699% | -1.1047% | -1.9769% | -1.2207% | -0.8855% | -4.6059% | -2.5394% | ... | -2.2066% | -3.0309% | -1.0410% | -3.1557% | -1.6148% | -0.2700% | -2.2845% | -2.0570% | -0.5492% | -3.0020% |
2016-01-08 | 0.2736% | -2.2705% | -1.6037% | -2.5425% | 0.1099% | -0.2241% | 0.5707% | -1.6402% | -1.7641% | -0.1649% | ... | -0.1538% | -1.1366% | -0.7308% | -0.1449% | 0.0895% | -3.3838% | -0.1117% | -1.1387% | -0.9720% | -1.1254% |
2016-01-11 | -4.3383% | 0.1692% | -1.6851% | -1.0215% | 0.0914% | -1.1792% | 0.5674% | 0.5288% | 0.6616% | 0.0330% | ... | 1.6436% | 0.0000% | 0.9869% | -0.1450% | 1.2225% | 1.4570% | 0.5367% | -0.4607% | 0.5800% | -1.9919% |
5 rows × 25 columns
####################################
# Minimizing Portfolio Variance
####################################
import cvxpy as cp
from timeit import default_timer as timer
# Defining initial inputs
mu = Y.mean().to_numpy().reshape(1,-1)
sigma = Y.cov().to_numpy()
# Defining initial variables
x = cp.Variable((mu.shape[1], 1))
# Budget and weights constraints
constraints = [cp.sum(x) == 1,
x <= 1,
x >= 0]
# Defining risk objective
risk = cp.quad_form(x, sigma)
objective = cp.Minimize(risk)
weights = pd.DataFrame([])
# Solving the problem with several solvers
prob = cp.Problem(objective, constraints)
solvers = ['ECOS', 'SCS']
for i in solvers:
prob.solve(solver=i)
# Showing Optimal Weights
weights_1 = pd.DataFrame(x.value, index=assets, columns=[i])
weights = pd.concat([weights, weights_1], axis=1)
display(weights)
ECOS | SCS | |
---|---|---|
APA | 0.0002% | -0.0007% |
BA | 0.0012% | 0.0006% |
BAX | 5.2647% | 5.2662% |
BMY | 4.3893% | 4.3904% |
CMCSA | 2.1705% | 2.1772% |
CNP | 6.9881% | 6.9878% |
CPB | 3.2388% | 3.2403% |
DE | 0.0707% | 0.0874% |
HPQ | 0.0001% | -0.0005% |
JCI | 2.8381% | 2.8435% |
JPM | 6.9786% | 6.9421% |
LUV | 2.8524% | 2.8590% |
MMC | 12.5818% | 12.6038% |
MO | 7.2325% | 7.2291% |
MSFT | 0.0002% | -0.0006% |
NI | 11.4513% | 11.4451% |
PCAR | 0.0003% | -0.0019% |
PSA | 14.9188% | 14.9124% |
SEE | 0.1563% | 0.1671% |
T | 6.4205% | 6.4208% |
TGT | 4.0908% | 4.0965% |
TMO | 0.0004% | 0.0005% |
TXT | 0.0007% | -0.0021% |
VZ | 8.3448% | 8.3447% |
ZION | 0.0087% | -0.0074% |
As we can see the use of CVXPY's quad_form in portfolio optimization can give small negative values to weights that must be zero.
# Calculating Annualized Portfolio Stats
var = weights * (Y.cov() @ weights) * 252
var = var.sum().to_frame().T
std = np.sqrt(var)
ret = Y.mean().to_frame().T @ weights * 252
stats = pd.concat([ret, std, var], axis=0)
stats.index = ['Return', 'Std. Dev.', 'Variance']
display(stats)
ECOS | SCS | |
---|---|---|
Return | 12.8507% | 12.8498% |
Std. Dev. | 10.3737% | 10.3738% |
Variance | 1.0761% | 1.0761% |
The maximization of portfolio return is a problem with a quadratic constraint that can be posed as:
$$ \begin{equation} \begin{aligned} & \underset{x}{\text{max}} & & \mu x^{\tau} \\ & \text{s.t.} & & x^{\tau} \Sigma x \leq \bar{\sigma}^{2} \\ & & & \sum_{i=1}^{N} x_i = 1 \\ & & & x_i \geq 0 \; ; \; \forall \; i =1, \ldots, N \\ \end{aligned} \end{equation} $$Where $x$ are the weights of assets and $\bar{\sigma}$ is the maximum expected standard deviation of portfolio..
#########################################
# Maximizing Portfolio Return with
# Variance Constraint
#########################################
import cvxpy as cp
from timeit import default_timer as timer
# Defining initial inputs
mu = Y.mean().to_numpy().reshape(1,-1)
sigma = Y.cov().to_numpy()
# Defining initial variables
x = cp.Variable((mu.shape[1], 1))
sigma_hat = 15 / (252**0.5 * 100)
ret = mu @ x
# Budget and weights constraints
constraints = [cp.sum(x) == 1,
x <= 1,
x >= 0]
# Defining risk constraint and objective
risk = cp.quad_form(x, sigma)
constraints += [risk <= sigma_hat**2] # variance constraint
objective = cp.Maximize(ret)
weights = pd.DataFrame([])
# Solving the problem with several solvers
prob = cp.Problem(objective, constraints)
solvers = ['ECOS', 'SCS']
for i in solvers:
prob.solve(solver=i)
# Showing Optimal Weights
weights_1 = pd.DataFrame(x.value, index=assets, columns=[i])
weights = pd.concat([weights, weights_1], axis=1)
display(weights)
ECOS | SCS | |
---|---|---|
APA | 0.0000% | 0.0006% |
BA | 9.5892% | 9.6764% |
BAX | 12.6176% | 12.5937% |
BMY | 0.0000% | 0.0051% |
CMCSA | 0.0000% | 0.0072% |
CNP | 3.7021% | 3.2811% |
CPB | 0.0000% | 0.0089% |
DE | 6.2565% | 6.3273% |
HPQ | 0.0000% | -0.0003% |
JCI | 0.0000% | 0.0053% |
JPM | 6.5739% | 6.5254% |
LUV | 0.0000% | 0.0048% |
MMC | 18.7461% | 18.5184% |
MO | 0.0000% | 0.0161% |
MSFT | 35.7121% | 36.2356% |
NI | 0.0278% | 0.0099% |
PCAR | 0.0000% | 0.0008% |
PSA | 0.0000% | 0.0179% |
SEE | 0.0000% | 0.0057% |
T | 0.0000% | 0.0126% |
TGT | 6.7741% | 6.7077% |
TMO | 0.0002% | 0.0011% |
TXT | 0.0000% | 0.0028% |
VZ | 0.0001% | 0.0117% |
ZION | 0.0001% | -0.0002% |
The small negative values also appear when we use CVXPY's quad_form in constraints.
# Calculating Annualized Portfolio Stats
var = weights * (Y.cov() @ weights) * 252
var = var.sum().to_frame().T
std = np.sqrt(var)
ret = Y.mean().to_frame().T @ weights * 252
stats = pd.concat([ret, std, var], axis=0)
stats.index = ['Return', 'Std. Dev.', 'Variance']
display(stats)
ECOS | SCS | |
---|---|---|
Return | 26.0352% | 26.1001% |
Std. Dev. | 15.0000% | 15.0672% |
Variance | 2.2500% | 2.2702% |
An alternative problem is to minimize the standard deviation (square root of variance). To do this we can use the SOC constraint. The minimization of portfolio standard deviation can be posed as:
$$ \begin{equation} \begin{aligned} & \underset{x}{\text{min}} & & g \\ & \text{s.t.} & & \mu x^{\tau} \geq \bar{\mu} \\ & & & \sum_{i=1}^{N} x_i = 1 \\ & & & \left\|\Sigma^{1/2} x\right\| \leq g \\ & & & x_i \geq 0 \; ; \; \forall \; i =1, \ldots, N \\ \end{aligned} \end{equation} $$Where $\left\|\Sigma^{1/2} x\right\| \leq g$ is the SOC constraint, $x$ are the weights of assets, $\mu$ is the mean vector of expected returns, $\bar{\mu}$ the minimum expected return of portfolio and $r$ is the matrix of observed returns.
Note: the SOC constraint can be expressed as $(g,\Sigma^{1/2} x) \in Q^{n+1}$, this notation is used to model the SOC constraint in CVXPY.
#########################################
# Minimizing Portfolio Standard Deviation
#########################################
from scipy.linalg import sqrtm
# Defining initial inputs
mu = Y.mean().to_numpy().reshape(1,-1)
sigma = Y.cov().to_numpy()
G = sqrtm(sigma)
# Defining initial variables
x = cp.Variable((mu.shape[1], 1))
g = cp.Variable(nonneg=True)
# Budget and weights constraints
constraints = [cp.sum(x) == 1,
x >= 0]
# Defining risk objective
risk = g
constraints += [cp.SOC(g, G @ x)] # SOC constraint
constraints += [risk <= sigma_hat] # variance constraint
objective = cp.Minimize(risk)
weights = pd.DataFrame([])
# Solving the problem with several solvers
prob = cp.Problem(objective, constraints)
solvers = ['ECOS', 'SCS']
for i in solvers:
prob.solve(solver=i)
# Showing Optimal Weights
weights_1 = pd.DataFrame(x.value, index=assets, columns=[i])
weights = pd.concat([weights, weights_1], axis=1)
display(weights)
ECOS | SCS | |
---|---|---|
APA | 0.0000% | 0.0000% |
BA | 0.0000% | 0.0001% |
BAX | 5.2634% | 5.2632% |
BMY | 4.3887% | 4.3886% |
CMCSA | 2.1705% | 2.1705% |
CNP | 6.9872% | 6.9870% |
CPB | 3.2390% | 3.2391% |
DE | 0.0789% | 0.0791% |
HPQ | 0.0000% | 0.0002% |
JCI | 2.8376% | 2.8375% |
JPM | 6.9842% | 6.9839% |
LUV | 2.8523% | 2.8522% |
MMC | 12.5805% | 12.5803% |
MO | 7.2314% | 7.2313% |
MSFT | 0.0000% | 0.0003% |
NI | 11.4509% | 11.4508% |
PCAR | 0.0000% | 0.0001% |
PSA | 14.9186% | 14.9183% |
SEE | 0.1621% | 0.1628% |
T | 6.4196% | 6.4196% |
TGT | 4.0904% | 4.0904% |
TMO | 0.0000% | 0.0002% |
TXT | 0.0000% | 0.0001% |
VZ | 8.3447% | 8.3446% |
ZION | 0.0000% | 0.0001% |
As we can see the use of CVXPY's SOC constraint in portfolio optimization solves the error that we see when we use quad_form.
# Calculating Annualized Portfolio Stats
var = weights * (Y.cov() @ weights) * 252
var = var.sum().to_frame().T
std = np.sqrt(var)
ret = Y.mean().to_frame().T @ weights * 252
stats = pd.concat([ret, std, var], axis=0)
stats.index = ['Return', 'Std. Dev.', 'Variance']
display(stats)
ECOS | SCS | |
---|---|---|
Return | 12.8508% | 12.8509% |
Std. Dev. | 10.3737% | 10.3737% |
Variance | 1.0761% | 1.0761% |
The maximization of portfolio return using SOC constraints can be posed as:
$$ \begin{equation} \begin{aligned} & \underset{x}{\text{max}} & & \mu x^{\tau} \\ & \text{s.t.} & & g \leq \bar{\sigma} \\ & & & \left\|\Sigma^{1/2} x\right\| \leq g \\ & & & \sum_{i=1}^{N} x_i = 1 \\ & & & x_i \geq 0 \; ; \; \forall \; i =1, \ldots, N \\ \end{aligned} \end{equation} $$#########################################
# Maximizing Portfolio Return with
# Standard Deviation Constraint
#########################################
from scipy.linalg import sqrtm
# Defining initial inputs
mu = Y.mean().to_numpy().reshape(1,-1)
sigma = Y.cov().to_numpy()
G = sqrtm(sigma)
# Defining initial variables
x = cp.Variable((mu.shape[1], 1))
g = cp.Variable(nonneg=True)
sigma_hat = 15 / (252**0.5 * 100)
ret = mu @ x
# Budget and weights constraints
constraints = [cp.sum(x) == 1,
x <= 1,
x >= 0]
# Defining risk constraint and objective
risk = g
constraints += [cp.SOC(g, G @ x)] # SOC constraint
constraints += [risk <= sigma_hat] # standard deviation constraint
objective = cp.Maximize(ret)
weights = pd.DataFrame([])
# Solving the problem with several solvers
prob = cp.Problem(objective, constraints)
solvers = ['ECOS', 'SCS']
for i in solvers:
prob.solve(solver=i)
# Showing Optimal Weights
weights_1 = pd.DataFrame(x.value, index=assets, columns=[i])
weights = pd.concat([weights, weights_1], axis=1)
display(weights)
ECOS | SCS | |
---|---|---|
APA | 0.0000% | 0.0003% |
BA | 9.5885% | 9.5835% |
BAX | 12.6190% | 12.6236% |
BMY | 0.0000% | -0.0011% |
CMCSA | 0.0000% | -0.0010% |
CNP | 3.7154% | 3.7281% |
CPB | 0.0000% | -0.0014% |
DE | 6.2555% | 6.2519% |
HPQ | 0.0000% | 0.0006% |
JCI | 0.0000% | -0.0006% |
JPM | 6.5725% | 6.5909% |
LUV | 0.0000% | -0.0008% |
MMC | 18.7512% | 18.7542% |
MO | 0.0000% | -0.0017% |
MSFT | 35.7087% | 35.6702% |
NI | 0.0135% | 0.0331% |
PCAR | 0.0000% | -0.0001% |
PSA | 0.0000% | -0.0018% |
SEE | 0.0000% | -0.0007% |
T | 0.0000% | -0.0015% |
TGT | 6.7755% | 6.7784% |
TMO | 0.0001% | 0.0001% |
TXT | 0.0000% | 0.0001% |
VZ | 0.0000% | -0.0013% |
ZION | 0.0000% | 0.0009% |
CVXPY's SOC constraint also solves the error that we see when we use quad_form in constraints.
# Calculating Annualized Portfolio Stats
var = weights * (Y.cov() @ weights) * 252
var = var.sum().to_frame().T
std = np.sqrt(var)
ret = Y.mean().to_frame().T @ weights * 252
stats = pd.concat([ret, std, var], axis=0)
stats.index = ['Return', 'Std. Dev.', 'Variance']
display(stats)
ECOS | SCS | |
---|---|---|
Return | 26.0352% | 26.0316% |
Std. Dev. | 15.0000% | 14.9958% |
Variance | 2.2500% | 2.2487% |
For more portfolio optimization models and applications, you can see the CVXPY based library Riskfolio-Lib.