Portfolio optimization with risk as constraint?

The standard mean-variance portfolio optimization models have the form:

Model M1
\[\begin{align}\min\>&\color{darkred}{\mathit{Risk}}=\color{darkred}x^T\color{darkblue}Q\color{darkred}x\\ & \color{darkblue}r^T\color{darkred}x \ge \color{darkblue}{\mathit{MinimumReturn}} \\  &\sum_i \color{darkred}x_i = 1\\ & \color{darkred}x \ge 0\end{align}\]

or may be:

Model M2
\[\begin{align}\min\>&\color{darkred}z=\color{darkred}x^T\color{darkblue}Q\color{darkred}x - \color{darkblue}\lambda \cdot\color{darkblue}r^T\color{darkred}x\\ &\sum_i \color{darkred}x_i = 1\\ & \color{darkred}x \ge 0\end{align}\]

where

• \(\color{darkblue}Q\) is a variance-covariance matrix (we assume here it is positive semi-definite [1]),
• \(\color{darkblue}\lambda\) is an exogenous constant that sometimes is varied to draw an efficient frontier, and
• \(\color{darkblue}r\) are the returns.

The last variant is:

Model M3
\[\begin{align}\max\>&\color{darkred}{\mathit{Return}}=\color{darkblue}r^T\color{darkred}x\\ & \color{darkred}x^T\color{darkblue}Q\color{darkred}x \le\color{darkblue}{\mathit{MaximumRisk}}\\ &\sum_i \color{darkred}x_i = 1 \\ & \color{darkred}x \ge 0\end{align}\]

I almost never encounter this model M3. This is a convex quadratically constrained problem, which can be solved quite easily using modern solvers. Note that the minimum-risk constraint must be an inequality in order for the model to stay convex.  (In model M1 we could have replaced the minimum return constraint with an equality constraint. But for quadratic constraints this is not the case: any nonlinear equality constraint is non-convex.)

I can invent two reasons why this model is not very popular:

• Traditionally, models like M1 and M2 have been solved with QP (Quadratic Programming) solvers. QP models have linear constraints and a (convex) quadratic objective. Solvers for quadratically constrained problems are newer. Nowadays, QP problems maybe even reformulated into a second-order cone representation. (All models discussed here are SOCP representable [2].) Essentially things have turned around: moderns solvers may like to have the risk in the constraints.
• It may be not so easy to know a reasonable numeric value for the risk. So having the risk as a constraint may not be very intuitive.

My question: is a QP formulation (model M1) still faster than a QCP (Quadratically Constrained Problem) model (i.e. model M3)? I have no clue. Let's make a small experiment.

Experiment 1

Here we solve a single model M1 and model M3 and compare results. They look like:

----    140 PARAMETER results  

                    M1          M3

MinReturn        0.500
MaxRisk                      0.098
num assets     350.000     350.000
vars           351.000     351.000
equs             3.000       3.000
status         Optimal     Optimal
obj              0.098       0.500
time             0.875       0.469

The exogenous value for minReturn is 0.5 (this was return over 4 years). The value for MaxRisk in model M3 was set to optimal risk in model M1. (So we should get as optimal objective in M3 our initial MinReturn=0.5). 

A bug in GAMS does not allow me to retrieve the iteration counts, but they were 16 and 15 interior-point iterations (and not 0 as reported by GAMS). This comparison does not say much: the problem is too easy. 

 

Experiment 2

Here we add a cardinality constraint to the problem: allow just 10 assets in our portfolio. This makes the problem an integer programming problem: I added binary variables to the model. 

The models are:

Model M1mip	Model M3mip
\[\begin{align}\min\>&\color{darkred}{\mathit{Risk}}=\color{darkred}x^T\color{darkblue}Q\color{darkred}x\\ & \color{darkblue}r^T\color{darkred}x \ge \color{darkblue}{\mathit{MinimumReturn}} \\  &\sum_i \color{darkred}x_i = 1\\ &\color{darkred}x \le \color{darkred}\delta \\& \sum_i \color{darkred}\delta_i = 10 \\ & \color{darkred}x \ge 0\\ & \color{darkred}\delta \in \{0,1\}\end{align}\]	\[\begin{align}\max\>&\color{darkred}{\mathit{Return}}=\color{darkblue}r^T\color{darkred}x\\ & \color{darkred}x^T\color{darkblue}Q\color{darkred}x \le\color{darkblue}{\mathit{MaximumRisk}}\\ &\sum_i \color{darkred}x_i = 1 \\ &\color{darkred}x \le \color{darkred}\delta \\& \sum_i \color{darkred}\delta_i = 10 \\ & \color{darkred}x \ge 0\\ & \color{darkred}\delta \in \{0,1\}\end{align}\]

Here are the results:

 

----    160 PARAMETER results  

                    M1          M3       M1mip       M3mip

MinReturn        0.500                   0.500
MaxRisk                      0.098                   0.099
num assets     350.000     350.000     350.000     350.000
vars           351.000     351.000     701.000     701.000
  discr                                350.000     350.000
equs             3.000       3.000     354.000     354.000
status         Optimal     Optimal     Optimal     Optimal
obj              0.098       0.500       0.099       0.500
time             0.875       0.469       3.172      15.203
nodes                                  559.000    6845.000
iterations                            1588.000   84589.000

Here we see that the M1mip model is doing a bit better than the M3mip model. Just one observation, so not sure whether we can draw conclusions from this. It could be that the quadratic objective helps the branch-and-bound process a bit. 

Conclusions

• Modern quadratic solvers have no problem with convex quadratically constrained models. When coming from non-linear programming this is a bit different: for general NLP models, we often prefer to have the nonlinearities to be in the objective.
• But the devil is in the details. Our little MIP example seems to prefer a quadratic objective above a quadratic constraint.

References

1. Covariance Matrix not positive definite in portfolio models, https://yetanothermathprogrammingconsultant.blogspot.com/2018/04/covariance-matrix-not-positive-definite.html
2. Second-order Cone Programming, https://en.wikipedia.org/wiki/Second-order_cone_programming

Appendix: GAMS and R code

$ontext   portfolio models   min risk, s.t. returns >= MINRET   max returns, s.t. risk <= MAXRISK $offtext $set script     script.R $set gdx        data.gdx $set rcommand   '"c:\Program Files\R\R-4.0.3\bin\Rscript.exe" --vanilla %script%' $onecho > %script% options(echo=TRUE) library(pprobeData) library(Matrix) library(reshape2) library(gdxrrw) # data in pprobeData: #    xassetPrices #    xassetLogReturns # I rescale the risk/return a bit r <- xassetLogReturns[1:1000,] mu <- colSums(r) summary(mu) # make sure the Covariance matrix is PD covmat <- 1000.0*as.matrix(nearPD(cov(r))$mat) summary(as.vector(covmat)) # # organize data to be exported to GAMS # mudf <- data.frame(     as.factor(names(mu)),     mu ) names(mudf) <- c("i","mu") attr(mudf,"symName") = "mu" attr(mudf,"domains") = c("i") attr(mudf,"domInfo") = "relaxed" covdf <- melt(covmat) covdf2 <- data.frame(     as.factor(covdf$Var1),     as.factor(covdf$Var2),     covdf$value ) names(covdf2) <- c("i","j","cov") attr(covdf2,"symName") = "cov" attr(covdf2,"domains") = c("i","i") attr(covdf2,"domInfo") = "relaxed" wgdx.lst("%gdx%",mudf,covdf2) $offecho $call '%rcommand%' set dummy 'for ordering' /MinReturn,MaxRisk/; set i; parameters     mu(i)      'returns'     cov(i,i)   'covariance matrix' ; $gdxin  %gdx% $loaddc i<mu.dim1 mu cov display i; alias(i,j); *---------------------------------------------- * models *---------------------------------------------- scalars    minret  'minimum return' /0.5/    maxrisk  'maximum allowed risk'    num    'number of assets we can choose' /10/ ; positive variable x(i) 'fraction to invest'; binary variable delta(i); variable z 'objective'; equation     qobj         'quadratic objective'     lobj         'linear objective'     budget       'allocations sum up to 1'     minreturn    'minimum return constraint (linear)'     emaxrisk     'maximum risk constraint (quadratic)'     notused      'delta=0 => x=0'     numassets    'number of assets in portfolio' ; * objectives qobj..    z =e= sum((i,j),x(i)*cov(i,j)*x(j)); lobj..    z =e= sum(i,mu(i)*x(i)); * return/risk constraints minreturn..  sum(i, mu(i)*x(i)) =g= minret; emaxrisk..  sum((i,j),x(i)*cov(i,j)*x(j)) =l= maxrisk; * allocations are fractions budget.. sum(i, x(i)) =e= 1; * cardinality constraints notused(i).. x(i) =l= delta(i); numassets.. sum(i, delta(i)) =e= num; * continuous models model m1 /qobj,budget,minreturn/; model m3 /lobj,budget,emaxrisk/; * cardinality constrained models model m1mip /qobj,budget,minreturn,notused,numassets/; model m3mip /lobj,budget,emaxrisk,notused,numassets/; *---------------------------------------------- * reporting *---------------------------------------------- acronym Optimal,NonOptimal; parameter results(*,*); $macro report(m,label,ex,exval)  \     results(ex,label) = exval;   \     results('num assets',label) = card(i); \     results('vars',label) = m.numVar; \     results('  discr',label) = m.numDVar; \     results('equs',label) = m.numEqu; \     results('status',label)$(m.solvestat=1) = Optimal; \     results('status',label)$(m.solvestat<>1) = NonOptimal; \     results('obj',label) = z.l; \     results('time',label) = m.resusd; \     results('nodes',label)$(m.nodusd<>NA) = m.nodusd; \     results('iterations',label) = m.iterusd; *---------------------------------------------- * solve *---------------------------------------------- option qcp=cplex,miqcp=cplex,optcr=0,threads=8, bratio=1; solve m1 minimizing z using qcp; report(m1,'M1','MinReturn',minret) maxrisk = z.l; solve m3 maximizing z using qcp; report(m3,'M3','MaxRisk',maxrisk) display results; solve m1mip minimizing z using miqcp; display x.l,delta.l,z.l; report(m1mip,'M1mip','MinReturn',minret) maxrisk = z.l; solve m3mip maximizing z using miqcp; display x.l,delta.l,z.l; report(m3mip,'M3mip','MaxRisk',maxrisk) display results;