The model is a bit difficult to decipher in this R code, but it is a simplified traditional MV model: Mean-Variance Portfolio Optimization Model \[\begin{align} \min\> & \color{darkred}x'\color{darkblue}Q\color{darkred}x - \color{darkblue}\lambda \bar{\color{darkblue}r}'\color{darkred}x\\ & \sum_i \color{darkred}x_i = 1 \\ & \color{darkred}x_i \ge 0\end{align}\] Instead of having the expected returns in the objective, we also see versions of this model with a linear constraint that puts a minimum value on the expected return of the portfolio. In the example, we have \(\lambda=0\), so the objective only has a quadratic term to minimize the risk. We don't care about returns. Note that the data only has 10 observations but 100 instruments. This situation can often lead to some problems. Indeed when we run this, we see: solQP <- solve.QP(cov(mData), zeros, t(aMat), bVec, meq = 1) As a side note: the data is just completely random, so we cannot expect useful solutions. However, at least, we should be able to produce some solutions. Solution1The covariance matrix is in theory positive-semi definite. We can see this from using the definition of the covariance [2]. However, due to floating-point inaccuracies, we can end up with a covariance matrix that is just slightly non-positive semi-definite. We can see this by inspecting the eigenvalues: Q <- cov(mData) eigen(Q)$values [1] 3.977635e-02 3.789523e-02 3.258164e-02 3.093782e-02 2.605389e-02 2.466463e-02 2.127899e-02 [8] 1.897506e-02 1.627900e-02 1.871945e-17 1.454315e-17 8.114384e-18 5.109196e-18 4.392427e-18 [15] 3.680766e-18 2.728485e-18 2.649302e-18 2.192073e-18 1.937634e-18 1.914487e-18 1.856811e-18 [22] 1.449321e-18 1.391434e-18 1.211216e-18 1.196721e-18 1.127119e-18 9.776146e-19 9.670197e-19 [29] 8.257462e-19 7.572554e-19 6.283522e-19 6.030820e-19 4.056541e-19 3.610128e-19 3.604800e-19 [36] 2.770106e-19 2.530797e-19 2.263314e-19 2.133854e-19 2.089700e-19 2.024130e-19 2.006560e-19 [43] 1.688435e-19 9.412675e-20 7.033666e-20 7.012531e-20 5.942729e-20 5.673501e-20 4.444257e-20 [50] 4.259787e-20 3.124233e-21 -2.017011e-20 -2.269857e-20 -3.034668e-20 -3.075493e-20 -3.262642e-20 [57] -4.345710e-20 -6.351674e-20 -7.347496e-20 -7.739535e-20 -8.065131e-20 -9.318415e-20 -1.004645e-19 [64] -1.048992e-19 -1.185271e-19 -1.315539e-19 -1.362383e-19 -1.376705e-19 -1.456396e-19 -1.697860e-19 [71] -1.700317e-19 -1.745068e-19 -1.794480e-19 -1.835296e-19 -1.960933e-19 -2.017457e-19 -2.332259e-19 [78] -2.377082e-19 -3.232830e-19 -3.468225e-19 -3.844707e-19 -4.310383e-19 -4.467493e-19 -4.540707e-19 [85] -4.977900e-19 -6.159103e-19 -6.338311e-19 -6.484161e-19 -6.489508e-19 -6.551306e-19 -6.664641e-19 [92] -7.276688e-19 -7.422977e-19 -7.927766e-19 -1.106209e-18 -1.707077e-18 -1.742870e-18 -2.465438e-18 [99] -1.859346e-17 -3.473638e-17 Most of the eigenvalues are basically zero with quite a few just a bit negative. A positive-definite matrix would have only positive values. A positive semi-definite (PSD) matrix would allow positive and zero eigenvalues but no negative ones. There are algorithms to find a close-by positive-definite matrix. In R this is available as nearPD() in the Matrix package. Let's try this out: Q <- nearPD(cov(mData))$mat solQP <- solve.QP(Q, zeros, t(aMat), bVec, meq = 1) solQP$solution [1] 0.007861808 0.004649581 0.012677143 0.009135964 0.011926439 0.013624409 0.007916747 0.004307772 [9] 0.011808486 0.011085302 0.010229029 0.013080240 0.010085504 0.010998879 0.011099784 0.008520520 [17] 0.007451730 0.012424371 0.012677330 0.007013515 0.009070013 0.010422999 0.011041407 0.008768456 [25] 0.007594051 0.010845828 0.008091472 0.010219700 0.017100764 0.005583947 0.011981802 0.007366402 [33] 0.013770522 0.009063828 0.010744898 0.007836608 0.009580356 0.007065588 0.008219370 0.010592597 [41] 0.011775598 0.011112051 0.005555992 0.014596684 0.016019481 0.008870919 0.010940604 0.005662200 [49] 0.007056997 0.012274299 0.006510091 0.012764160 0.008258651 0.009945869 0.009891680 0.007318763 [57] 0.014960558 0.008255708 0.009655097 0.010054760 0.007309219 0.008761360 0.011616297 0.010300254 [65] 0.012043733 0.011963495 0.009991699 0.010353212 0.013896332 0.011619403 0.007160011 0.007202327 [73] 0.008530819 0.010918658 0.004754215 0.008877062 0.008526937 0.010251250 0.008190446 0.006754629 [81] 0.013209816 0.014849748 0.012213860 0.008888433 0.018188014 0.009756690 0.014786487 0.014893706 [89] 0.009737739 0.013256532 0.009990737 -0.001155889 0.009264992 0.007733719 0.014505131 0.011040198 [97] 0.012859850 0.008987104 0.008458713 0.004497735 Much better. When printing the eigenvalues of the new Q matrix, we see: eigen(Q)$values [1] 3.977635e-02 3.789522e-02 3.258164e-02 3.093782e-02 2.605389e-02 2.466462e-02 2.127899e-02 1.897506e-02 [9] 1.627900e-02 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 [17] 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 [25] 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 [33] 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 [41] 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 [49] 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 [57] 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 [65] 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 [73] 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977634e-10 3.977634e-10 3.977634e-10 3.977634e-10 [81] 3.977634e-10 3.977634e-10 3.977634e-10 3.977634e-10 3.977634e-10 3.977634e-10 3.977634e-10 3.977634e-10 [89] 3.977634e-10 3.977634e-10 3.977634e-10 3.977634e-10 3.977634e-10 3.977634e-10 3.977634e-10 3.977634e-10 [97] 3.977634e-10 3.977634e-10 3.977633e-10 3.977633e-10 All small (and negative) eigenvalues are now mapped to some tolerance. This means the Q matrix is now numerically positive definite. Solution2A simple DIY approach would be to fix the Q matrix by adding a small number to the diagonal. The following code shows a small perturbation is enough: Q <- cov(mData) + diag(1e-6,nA) Some solvers do something like this automatically. Fixing the modelLet's first fix our model a bit:
The portfolio model becomes: Portfolio Model algebraic formulation matrix notation \[\begin{align} \min & \sum_{i,j} \color{darkblue}q_{i,j} \color{darkred}x_i \color{darkred}x_j \\ & \sum_i \color{darkred}x_i = 1 \\ & \sum_i \bar{\color{darkblue}r}_i \color{darkred}x_i \ge 0.035\\ & \color{darkred}x_i \ge 0\end{align}\]\[\begin{align} \min \>& \color{darkred}x' \color{darkblue}Q \color{darkred}x \\ & \color{darkblue}e'\color{darkred}x = 1 \\ & \bar{\color{darkblue}r}' \color{darkred}x \ge 0.035\\ & \color{darkred}x \ge 0\end{align}\] QuadProg does not have facilities for bounds on the variables so we need to form explicit constraints. Its input-model looks like: \[\begin{align} \min \> & -d'b + \frac{1}{2} b'Db \\ & A'b\ge b_0\end{align}\] where the first meq constraints are equalities. The factor 0.5 is a relic from quadratic optimization (it makes the gradient a bit easier). If the model includes linear objective coefficients, it is important to take this 0.5 factor into account. Our model has only quadratic objective coefficients, so we don't care in this case. All in all, this is often a bit of an inconvenient input format. R trick. In the examples below we want to display a sparse vector. A simple trick using names can help here. v <- c(0,0,3,0,0,0,1,0,0,0,2,0,0,0,0) v [1] 0 0 3 0 0 0 1 0 0 0 2 0 0 0 0 names(v) <- paste("v",1:length(v),sep="") v[v>0] v3 v7 v11 3 1 2 The R code for our optimization model can look like: Q <- nearPD(cov(mData))$mat Another issue with QuadProg is that requires a strictly positive definite matrix. This also means it does not allow linear variables: all variables need appear to quadratically in the objective (and without all zero quadratic coefficients). After all, if the quadratic coefficient matrix would look like \[D = \left[\begin{array}{c|c}Q & 0 \\ \hline 0 & 0 \end{array} \right]\] we would end up with a matrix that is not strictly positive-definite. Often, for these types of models, CVXR is a more agreeable tool. Here is the same model in CVXR, and solved with the OSQP QP solver: library(CVXR) Here we use our covariance matrix directly. Both CVXR and OSQP accept it, even though it has tiny negative numbers. However, the OSQP documentation warns that this may be a bit dangerous [3]: We recommend the user to check the convexity of their problem before passing it to OSQP! If the user passes a non-convex problem we do not assure the solver will be able to detect it. Notice that CVXR will complain if the Q matrix is really non-positive semi-definite. Solution3In [2] a different model is proposed for these situations where we have more instruments than observations. Means-adjusted-returns Portfolio Optimization Model \[\begin{align} \min\> & \sum_t \color{darkred}w_t^2 \\ &\color{darkred}w_t = \sum_i (\color{darkblue}r_{t,i} - \bar{\color{darkblue}r}_i)\color{darkred}x_i \\ & \sum_i \color{darkred}x_i = 1 \\ & \sum_i \bar{\color{darkblue}r}_i \color{darkred}x_i \ge 0.035\\ & \color{darkred}x_i \ge 0\end{align}\] Some of the advantages of this model are:
This model is easily implemented in CVXR: library(CVXR) Using QuadProg this is a much more difficult exercise. The quadratic coefficient matrix \(D\) will look like: \[D = \left[\begin{array}{c|c}Q & 0 \\ \hline 0 & \varepsilon I \end{array} \right]\] The R code can look like: solQP <- solve.QP(cov(mData), zeros, t(aMat), bVec, meq = 1) 0 ConclusionQuadProg is a well-known QP solver for R. Unfortunately it has some drawbacks:
For quick experimentation with these types of models, CVXR is likely a better environment. This exercise showed: |
