This notebook makes use of the latex_envs Jupyter extension for equations and references. If the LaTeX is not rendering properly in your local installation of Jupyter , it may be because you have not installed this extension. Details at https://jupyter-contrib-nbextensions.readthedocs.io/en/latest/nbextensions/latex_envs/README.html
We consider a classic portfolio problem with loss risk constraints:
maximizex∈ℜnˉrTxsubject toProb(rTx≤α)≤β,∑ni=1xi=1,x≥0. Here we suppose the return of aseets r∈ℜn is Gaussian with mean ˉr∈ℜn and covariance V∈ℜn×n, α is a given unwanted return level (e.g. an excessive lose) and β is a given maximum probability.We will demonstrate that problem (1) can be transformed into a second-order cone programming (SOCP) of the following form and solved easily using the SOCP solver introduced at Mark 27 in the nAG Library:
minimizex∈ℜncTxsubject tolA≤Ax≤uA,lx≤x≤ux,x∈K, where A∈ℜm×n, lA,uA∈ℜm, c,lx,ux∈ℜn are the problem data, and K=Kn1×⋯×Knr×ℜnl where Kni is either a quadratic cone or a rotated quadratic cone defined as follows:Let u=rTx with σ=xTVx denoting its variance, the probability constraint in problem (1) can be written as
Prob(u−ˉu√σ≤α−ˉu√σ)≤β. Note that (u−ˉu)/√σ is a standard Gaussian random variable, the probability above is simply Φ((α−ˉu)/√σ), where Φ(z)=1√2π∫z−∞e−t2/2dt is the CDF of the standard Gaussian random variable. Thus the probability constraint in problem (1) can be written as α−ˉu√σ≤Φ−1(β) or, equivalently, ˉu+Φ−1(β)√σ≥α. From ˉu=ˉrTx and σ=xTVx=‖Fx‖2 by factorizing V=FTF, (8) is equivalent to ˉrTx+Φ−1(β)‖Fx‖≥α. Depending on the value of β, (9) could be convex or concave. By setting β≤0.5 (which is reasonable for risk control), (9) is convex and can be written as ˉrTx+Φ−1(β)t=α and ‖Fx‖≤t. Note that by letting Fx=y, ‖Fx‖=‖y‖≤t in (10) fits exactly the quadratic cone constraint (3), thus can be solved by second-order cone programming. The final equivalent SOCP model to problem (1) is maximizex∈ℜnˉrTxsubject to∑ni=1xi=1,Fx=y,ˉrTx+Φ−1(β)t=α,(t;y)∈Kn+1q,x≥0.We demonstrate how to use the nAG SOCP solver to model and solve problem (1) by solving the equivalent SOCP problem (11).
# Import utility libraries and the nAG Library
import numpy as np
import math as mt
from scipy.sparse import coo_matrix
from scipy.stats import norm
from naginterfaces.library import opt
from naginterfaces.library import lapackeig
As an example, we set α=0.0001, β=0.05 and randomly generate ˉr and V for 8 assets as follows.
# Set alpha and beta
alpha = 0.0001
beta = 0.05
# Fix random seed
np.random.seed(9)
# Number of assets
n_assets = 8
# Vector of expected returns
r = np.ones(n_assets)*.02 + np.random.rand(n_assets)*.15
# Covariance matrix
V = np.matrix(np.random.randn(n_assets, n_assets))
V = V.T * V
V = V / np.max(np.abs(np.diag(V))) * .2
Now factorize V=FTF using eigenvalue decomposition from the nAG Library.
# Note one could use sparse factorization if V is input as sparse matrix
U, lamda = lapackeig.dsyevd('V', 'L', V)
# Find positive eigenvalues and corresponding eigenvectors
i = 0
k = 0
F = []
while i<len(lamda):
if lamda[i] > 0:
F = np.append(F, mt.sqrt(lamda[i])*U[:,i])
k += 1
i += 1
F = F.reshape((k, n_assets))
For modelling, nAG SOCP solver requires several input arguments for objective and constraints. Now we initialize the data that will be used to feed nAG SOCP solver.
# Number of variables
n = n_assets
# Number of constraints
m = 0
# Objective coefficient c
c = np.full(n, 0.0, float)
# Bounds on variables
blx = np.full(n, -1.e20, float)
bux = np.full(n, 1.e20, float)
# Linear constraint bu <= Ax <= bu
# A in coordinate list format (COO)
irowa = np.empty(0, int)
icola = np.empty(0, int)
a = np.empty(0, float)
# Bounds on Ax
bl = np.empty(0, float)
bu = np.empty(0, float)
# Cone constraints
ctype = []
group = []
Because we will add auxiliary variables and constraints during the process, it is necessary we keep tracking the number of variables and constraints in the model by maintaining the up-to-date problem size.
# Initialize the up-to-date problem size
n_up = n
m_up = m
Now we keep modifying the above data during adding objective and constraints one by one. First is the objective coefficient c.
# Add objective function min -r'x
c = -r
Then we add the long-only constraint.
# Number of linear constraints will increase by 1
m += 1
# Set lower bound on x to 0
blx[0:n] = np.zeros(n)
# Add sum(x) = 1
irowa = np.append(irowa, np.full(n, m_up+1, dtype=int))
icola = np.append(icola, np.arange(1, n+1))
a = np.append(a, np.full(n, 1.0, dtype=float))
bl = np.append(bl, np.full(1, 1.0, dtype=float))
bu = np.append(bu, np.full(1, 1.0, dtype=float))
Now add the probability constraint by adding Fx=y, ˉrTx+Φ−1(β)t=α and (t;y)∈Kn+1q.
# Get quantile function of beta
quantile = norm.ppf(beta)
# Up-to-date problem size
m_up = m
n_up = n
# Then k + 1 more variables need to be added together with
# k + 1 linear constraints and a second-order cone contraint
# Enlarge the model
n = n + k + 1
m = m + k + 1
# All the added auxiliary variables do not take part in obj
c = np.append(c, np.zeros(k+1))
# Enlarge bounds on x, add inf bounds on the new added k+1 variables
blx = np.append(blx, np.full(k+1, -1.e20, dtype=float))
bux = np.append(bux, np.full(k+1, 1.e20, dtype=float))
# Enlarge linear constraints
# Sparsity pattern of F (COO)
row, col = np.nonzero(F)
val = F[row, col]
# Convert to 1-based and move row down by m
# Add Fx = y and Phi^-1(beta)*t + r_bar'x - alpha = 0
# [x,t,y]
row = row + 1 + m_up
col = col + 1
row = np.append(row, np.arange(m_up+1, m_up+k+1+1))
col = np.append(np.append(col, np.arange(n_up+2, n_up+k+1+1)), n_up+1)
val = np.append(val, np.append(np.full(k, -1.0, dtype=float), quantile))
irowa = np.append(irowa, row)
icola = np.append(icola, col)
a = np.append(a, val)
bl = np.append(bl, np.append(np.zeros(k), alpha))
bu = np.append(bu, np.append(np.zeros(k), alpha))
# coeffient of x in Phi^-1(beta)*t + r'x - alpha = 0
irowa = np.append(irowa, np.full(n_assets, m_up+k+1, dtype=int))
icola = np.append(icola, np.arange(1, n_assets+1))
a = np.append(a, r)
# Enlarge cone constraints
ctype.extend('Q')
group_temp = np.arange(n_up+1, n_up+k+1+1)
group.append(group_temp)
By now, all the data we need is ready. Feed them to the nAG SOCP solver and solve.
# Create problem handle
handle = opt.handle_init(n)
# Set objective function
opt.handle_set_linobj(handle, c)
# Set box constraints
opt.handle_set_simplebounds(handle, blx, bux)
# Set linear constraints
opt.handle_set_linconstr(handle, bl, bu, irowa, icola, a)
# Set cone constraints
i = 0
while i<len(ctype):
opt.handle_set_group(handle, ctype[i], 0, group[i])
i += 1
# Set options
for option in [
'Print Options = NO',
'Print File = 1',
'SOCP Scaling = A'
]:
opt.handle_opt_set(handle, option)
# Call socp interior point solver
slt = opt.handle_solve_socp_ipm(handle)
Now we can print the optimal portfolio and the corresponding return.
# Optimal portfolio
slt.x[0:n_assets]
# Optimal expected return
r.dot(slt.x[0:n_assets])
In this notebook, we demonstrated how to use the nAG Library to model and solve portfolio optimization with probability constraint via SOCP. It is worth noting that SOCP is widely used in portfolio optimization due to its flexibility and versatility to handle a large variety of problems with different kinds of constraints, not only the probability constraint mentioned above, e.g. leverage constraint, turnover constraint, max position constraint and tracking-error constraint. We refer the readers to \cite{AG03, LVBL98, nAGDOC} for more details.
[1] Alizadeh Farid and Goldfarb Donald, ``Second-order cone programming'', Mathematical programming, vol. 95, number 1, pp. 3--51, 2003.
[2] Lobo Miguel Sousa, Vandenberghe Lieven, Boyd Stephen et al., ``Applications of second-order cone programming'', Linear algebra and its applications, vol. 284, number 1-3, pp. 193--228, 1998.
[3] Numerical Algorithms Group, ``nAG documentation'', 2019. online