add (Choi 2022) method

Author: mirca

riskparityportfolio/tests/test_vanilla_rpp.py

@@ -1,19 +1,22 @@
 import riskparityportfolio as rpp
 import numpy as np
+import pytest
 
-
-def test():
-    n = 1000
+@pytest.mark.parametrize("method", ["choi", "spinu"])
+def test(method):
+    n = 100
     t = 1000 * n
     x = np.random.normal(size=t).reshape((n, -1))
-    corr = x @ x.T / t
-    b = np.ones(len(corr)) / len(corr)
-    w = rpp.vanilla.design(corr, b)
-    rc = w @ (corr * w)
+    cov = np.cov(x)
+    print(cov.shape)
+    b = np.ones(len(cov)) / len(cov)
+    w = rpp.vanilla.design(cov, b, maxiter=1000, method=method)
+    print(w)
+    rc = w @ (cov * w)
     rc = rc / np.sum(rc)
     # assert that the portfolio respect the budget constraint
     np.testing.assert_almost_equal(np.sum(w), 1)
-    # assert that the portfolio repect the no-shortselling constraint
-    np.testing.assert_equal(all(w > 0), True)
+    # assert that the portfolio respect the no-shortselling constraint
+    np.testing.assert_equal(all(w >= 0), True)
     # assert that the desired risk contributions are attained
-    np.testing.assert_allclose(rc, b, rtol=1e-4)
+    np.testing.assert_allclose(rc, b, atol=1/(10*n))

riskparityportfolio/vanilla.cpp

@@ -1,19 +1,26 @@
+#include <cstring>
+#include <Eigen/Dense>
 #include <pybind11/pybind11.h>
 #include <pybind11/eigen.h>
 
 using namespace Eigen;
 using namespace std;
 
+typedef const Eigen::Matrix<double, Eigen::Dynamic, 1> c_vector_t;
+typedef const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> c_matrix_t;
+typedef Eigen::Matrix<double, Eigen::Dynamic, 1> vector_t;
+typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> matrix_t;
+
 // Cyclical coordinate descent for Spinu's formulation
 // of the risk parity portfolio problem
-Eigen::VectorXd risk_parity_portfolio_ccd_spinu(const Eigen::MatrixXd& Sigma,
-                                                const Eigen::VectorXd& b,
-                                                const double tol = 1E-8,
-                                                const unsigned int maxiter = 50) {
+vector_t risk_parity_portfolio_ccd_spinu(c_matrix_t& Sigma,
+                                         c_vector_t& b,
+                                         const double tol = 1E-4,
+                                         const unsigned int maxiter = 100) {
   double aux, x_diff, xk_sum;
   const unsigned int n = b.size();
-  Eigen::VectorXd xk = (1 / Sigma.diagonal().array().sqrt()).matrix();
-  Eigen::VectorXd x_star(n), Sigma_xk(n), rc(n);
+  vector_t xk = (1 / Sigma.diagonal().array().sqrt()).matrix();
+  vector_t x_star(n), Sigma_xk(n), rc(n);
   Sigma_xk = Sigma * xk;
   for (unsigned int k = 0; k < maxiter; ++k) {
     for (unsigned int i = 0; i < n; ++i) {
@@ -34,14 +41,56 @@ Eigen::VectorXd risk_parity_portfolio_ccd_spinu(const Eigen::MatrixXd& Sigma,
 }
 
 
+// cyclical coordinate descent algo by Choi & Chen 2022
+// ref: https://arxiv.org/pdf/2203.00148.pdf
+vector_t risk_parity_portfolio_ccd_choi(c_matrix_t& cov,
+                                        c_vector_t& b,
+                                        const double tol = 1E-4,
+                                        const unsigned int maxiter = 100) {
+  const unsigned int n = b.size();
+  vector_t a(n);
+  vector_t vol = cov.diagonal().array().sqrt();
+  matrix_t invvol_mat = (1 / vol.array()).matrix().asDiagonal();
+  matrix_t corr = invvol_mat * cov * invvol_mat;
+  matrix_t adj = corr;
+  adj.diagonal().array() = 0;
+  vector_t wk = vector_t::Ones(n);
+  wk = (wk.array() / std::sqrt(corr.sum())).matrix();
+  for (unsigned int k = 0; k < maxiter; ++k) {
+    // compute portfolio weights
+    a = 0.5 * adj * wk;
+    wk = ((a.array() * a.array() + b.array()).sqrt() - a.array()).matrix();
+    if ((wk.array() * (corr * wk).array() - b.array()).abs().maxCoeff() < tol)
+      break;
+  }
+  vector_t w = wk.array() / vol.array();
+  return (w / w.sum()).matrix();
+}
+
+
+vector_t rpp_design(c_matrix_t& cov,
+                    c_vector_t& b,
+                    const double tol = 1E-4,
+                    const unsigned int maxiter = 100,
+                    const char* method = "spinu"
+                    ){
+  if (strcmp(method, "spinu"))
+    return risk_parity_portfolio_ccd_spinu(cov, b, tol, maxiter);
+  else
+    return risk_parity_portfolio_ccd_choi(cov, b, tol, maxiter);
+}
+
 namespace py = pybind11;
 
 
 PYBIND11_MODULE(vanilla, m) {
   m.doc() = "design of risk parity portfolios";
-  m.def("design", &risk_parity_portfolio_ccd_spinu,
-        py::arg("Sigma"), py::arg("b"),
-        py::arg("tol") = 1E-8, py::arg("maxiter") = 50,
+  m.def("design", &rpp_design,
+        py::arg("Sigma"),
+        py::arg("b"),
+        py::arg("tol") = 1E-4,
+        py::arg("maxiter") = 50,
+        py::arg("method") = "spinu",
         R"pbdoc(
           A function to design vanilla risk parity (budgeting) portfolios.
 
@@ -63,6 +112,9 @@ PYBIND11_MODULE(vanilla, m) {
           maxiter : int
             maximum number of iterations. Default value is 50
 
+          method : str
+            which method to use. Available: "spinu" and "choi"
+
 
           Example
           -------