Merge pull request #31 from convexfi/add_dummy_variables

Changes: Finer control of the iterative algorithm and dummy variables allowed in the constraints.

.gitignore

@@ -5,7 +5,8 @@ build/
 riskparityportfolio.egg-info/
 var/
 .DS_Store
-.coverage
+.idea/
+.coverage*
 docs/source/tutorials/.ipynb_checkpoints/*
 .eggs/
 .ipynb_checkpoints/*

setup.py

@@ -4,7 +4,7 @@
 import setuptools
 import os
 
-__version__ = "0.5.1"
+__version__ = "0.6.0"
 
 # Prepare and send a new release to PyPI
 if "release" in sys.argv[-1]:

src/riskparityportfolio/rpp.py

@@ -18,14 +18,22 @@ def __init__(self):
 
 
 class RiskParityPortfolio:
-    """Designs risk parity portfolios by solving the following optimization problem
+    """Designs risk parity portfolios by solving the following optimization problem:
 
-    minimize R(w) - alpha * mu.T * w + lambda * w.T Sigma w
-    subject to Cw = c, Dw <= d
+      minimize     R(w) - alpha * mu.T w + lambda * w.T Sigma w
+      subject to   Cw = c, Dw <= d
 
-    where R is a risk concentration function, and alpha and lambda are trade-off
+    where R(w) is a risk concentration function, and alpha and lambda are trade-off
     parameters for the expected return and the variance, respectively.
 
+    The risk concentration R(w) is computed as the squared l2-norm of the risk concentration vector,
+    which by default is obtained from RiskContribOverVarianceMinusBudget():
+
+      R(w) = sum_i (MR_i/sum(MR_i) - budget_i)^2
+
+    where MR_i = w_i * (Sigma @ w)_i are the marginal risks (sum(MR_i) = Var(w)), and
+          budget_i are the risk budgets (by default 1/n).
+
     Parameters
     ----------
     covariance : array, shape=(n, n)
@@ -36,6 +44,20 @@ class RiskParityPortfolio:
         weights of the portfolio
     risk_concentration : class
         any valid child class of RiskConcentrationFunction
+
+    Examples:
+    # Set up:
+    >>> import numpy as np
+    >>> import riskparityportfolio as rpp
+    >>> n = 10
+    >>> U = np.random.multivariate_normal(mean=np.zeros(n), cov=0.1 * np.eye(n), size=round(.7 * n)).T
+    >>> Sigma = U @ U.T + np.eye(n)
+    # Basic usage with default constraints sum(w) = 1 and w >= 0:
+    >>> my_portfolio = rpp.RiskParityPortfolio(Sigma)
+    >>> my_portfolio.design()
+    >>> my_portfolio.weights
+    # Basic usage with equality and inequality constraints:
+    >>> my_portfolio.design(Cmat=Cmat, cvec=cvec, Dmat=Dmat, dvec=dvec)
     """
 
     def __init__(
@@ -172,13 +194,15 @@ def add_variance(self, lmd):
         self.lmd = lmd
         self.has_variance = True
 
-    def design(self, **kwargs):
+    def design(self, verbose=True, **kwargs):
         """Optimize the portfolio.
 
         Parameters
         ----------
+        verbose : boolean
+            Whether to print the optimization process.
         kwargs : dict
-            Dictionary of parameters to be passed to SuccessiveConvexOptimizer.
+            Dictionary of parameters to be passed to SuccessiveConvexOptimizer().
         """
         self.sca = SuccessiveConvexOptimizer(self, **kwargs)
-        self.sca.solve()
+        self.sca.solve(verbose)

src/riskparityportfolio/sca.py

@@ -1,4 +1,5 @@
 import numpy as np
+import warnings
 from tqdm import tqdm
 
 try:
@@ -101,9 +102,25 @@ def maxiter(self, value):
 
 class SuccessiveConvexOptimizer:
     """
-    Successive Convex Approximation optimizer tailored for the risk parity problem.
+    Successive Convex Approximation optimizer tailored for the risk parity problem including the linear constraints:
+       Cmat @ w  = cvec
+       Dmat @ w <= dvec,
+    where matrices Cmat and Dmat have n columns (n being the number of assets). Based on the paper:
+
+    Feng, Y., and Palomar, D. P. (2015). SCRIP: Successive convex optimization methods for risk parity portfolios design.
+    IEEE Trans. Signal Processing, 63(19), 5285–5300.
+
+    By default, the constraints are set to sum(w) = 1 and w >= 0, i.e.,
+       Cmat = np.ones((1, n))
+       cvec = np.array([1.0])
+       Dmat = -np.eye(n)
+       dvec = np.zeros(n).
+
+    Notes:
+        1) If equality constraints are not needed, set Cmat = np.empty((0, n)) and cvec = [].
+        2) If the matrices Cmat and Dmat have more than n columns, it is assumed that the additional columns
+        (same number for both matrices) correspond to dummy variables (which do not appear in the objective function).
     """
-
     def __init__(
         self,
         portfolio,
@@ -119,28 +136,27 @@ def __init__(
         dvec=None,
     ):
         self.portfolio = portfolio
-        self.tau = tau or 0.05 * np.sum(np.diag(self.portfolio.covariance)) / (
-            2 * self.portfolio.number_of_assets
-        )
+        self.tau = tau or 1e-4  # 0.05 * np.trace(self.portfolio.covariance) / (2 * self.portfolio.number_of_assets)
         sca_validator = SuccessiveConvexOptimizerValidator()
         self.gamma = sca_validator.gamma = gamma
         self.zeta = sca_validator.zeta = zeta
         self.funtol = sca_validator.funtol = funtol
         self.wtol = sca_validator.wtol = wtol
         self.maxiter = sca_validator.maxiter = maxiter
         self.Cmat = Cmat
-        self.Dmat = Dmat
+        self.Dmat = Dmat  # Dmat @ w <= dvec
         self.cvec = cvec
         self.dvec = dvec
-        self.CCmat = np.vstack((self.Cmat, self.Dmat)).T
-        self.bvec = np.concatenate((self.cvec, self.dvec))
+        self.number_of_vars = self.Cmat.shape[1]
+        self.number_of_dummy_vars = self.number_of_vars - self.portfolio.number_of_assets
+        self.dummy_vars = np.zeros(self.number_of_dummy_vars)
+        self.CCmat = np.vstack((self.Cmat, -self.Dmat)).T  # CCmat.T @ w >= bvec
+        self.bvec = np.concatenate((self.cvec, -self.dvec))
         self.meq = self.Cmat.shape[0]
         self._funk = self.get_objective_function_value()
         self.objective_function = [self._funk]
         self._tauI = self.tau * np.eye(self.portfolio.number_of_assets)
-        self.Amat = (
-            self.portfolio.risk_concentration.jacobian_risk_concentration_vector()
-        )
+        self.Amat = self.portfolio.risk_concentration.jacobian_risk_concentration_vector()
         self.gvec = self.portfolio.risk_concentration.risk_concentration_vector
 
     @property
@@ -151,7 +167,7 @@ def Cmat(self):
     def Cmat(self, value):
         if value is None:
             self._Cmat = np.atleast_2d(np.ones(self.portfolio.number_of_assets))
-        elif np.atleast_2d(value).shape[1] == self.portfolio.number_of_assets:
+        elif np.atleast_2d(value).shape[1] >= self.portfolio.number_of_assets:
             self._Cmat = np.atleast_2d(value)
         else:
             raise ValueError(
@@ -166,9 +182,9 @@ def Dmat(self):
     @Dmat.setter
     def Dmat(self, value):
         if value is None:
-            self._Dmat = np.eye(self.portfolio.number_of_assets)
-        elif np.atleast_2d(value).shape[1] == self.portfolio.number_of_assets:
-            self._Dmat = -np.atleast_2d(value)
+            self._Dmat = -np.eye(self.portfolio.number_of_assets)
+        elif np.atleast_2d(value).shape[1] == self.Cmat.shape[1]:
+            self._Dmat = np.atleast_2d(value)
         else:
             raise ValueError(
                 "Dmat shape {} doesnt agree with the number of"
@@ -200,7 +216,7 @@ def dvec(self, value):
         if value is None:
             self._dvec = np.zeros(self.portfolio.number_of_assets)
         elif len(value) == self.Dmat.shape[0]:
-            self._dvec = -np.atleast_1d(value)
+            self._dvec = np.atleast_1d(value)
         else:
             raise ValueError(
                 "dvec shape {} doesnt agree with Dmat shape"
@@ -215,42 +231,74 @@ def get_objective_function_value(self):
             obj += self.portfolio.lmd * self.portfolio.volatility ** 2
         return obj
 
-    def iterate(self):
+    def iterate(self, verbose=True):
         wk = self.portfolio.weights
         g = self.gvec(wk)
         A = np.ascontiguousarray(self.Amat(wk))
         At = np.transpose(A)
         Q = 2 * At @ A + self._tauI
-        q = 2 * np.matmul(At, g) - np.matmul(Q, wk)
+        q = 2 * np.matmul(At, g) - Q @ wk  # np.matmul() is necessary here since g is not a numpy array
         if self.portfolio.has_variance:
             Q += self.portfolio.lmd * self.portfolio.covariance
         if self.portfolio.has_mean_return:
             q -= self.portfolio.alpha * self.portfolio.mean
-        w_hat = quadprog.solve_qp(Q, -q, C=self.CCmat, b=self.bvec, meq=self.meq)[0]
-        self.portfolio.weights = wk + self.gamma * (w_hat - wk)
+        if self.number_of_dummy_vars > 0:
+            Q = np.vstack([np.hstack([Q, np.zeros((self.portfolio.number_of_assets, self.number_of_dummy_vars))]),
+                           np.hstack([np.zeros((self.number_of_dummy_vars, self.portfolio.number_of_assets)),
+                                      self.tau * np.eye(self.portfolio.number_of_assets)])])
+            q = np.concatenate([q, -self.tau * self.dummy_vars])
+        # Call QP solver (min 0.5*x.T G x + a.T x  s.t.  C.T x >= b) controlling for ill-conditioning:
+        try:
+            w_hat = quadprog.solve_qp(Q, -q, C=self.CCmat, b=self.bvec, meq=self.meq)[0]
+        except ValueError as e:
+            if str(e) == "matrix G is not positive definite":
+                warnings.warn(
+                    "Matrix Q is not positive definite: adding regularization term and then calling QP solver again.")
+                # eigvals = np.linalg.eigvals(Q)
+                # print("    - before regularization: cond. number = {:,.0f}".format(max(eigvals) / min(eigvals)))
+                # print("    - after regularization: cond. number = {:,.0f}".format(max(eigvals + np.trace(Q)/1e7) / min(eigvals + np.trace(Q)/1e7)))
+                Q += np.eye(Q.shape[0]) * np.trace(Q)/1e7
+                w_hat = quadprog.solve_qp(Q, -q, C=self.CCmat, b=self.bvec, meq=self.meq)[0]
+            else:
+                # If the error is different, re-raise it
+                raise
+        self.portfolio.weights = wk + self.gamma * (w_hat[:self.portfolio.number_of_assets] - wk)
         fun_next = self.get_objective_function_value()
         self.objective_function.append(fun_next)
         has_w_converged = (
-            np.abs(self.portfolio.weights - wk)
-            <= 0.5 * self.wtol * (np.abs(self.portfolio.weights) + np.abs(wk))
+                (np.abs(self.portfolio.weights - wk) <= self.wtol * 0.5 * (np.abs(self.portfolio.weights) + np.abs(wk)))
+                | ((np.abs(self.portfolio.weights) < 1e-6) & (np.abs(wk) < 1e-6))
         ).all()
         has_fun_converged = (
-            np.abs(self._funk - fun_next)
-            <= 0.5 * self.funtol * (np.abs(self._funk) + np.abs(fun_next))
-        ).all()
-        if has_w_converged or has_fun_converged:
+                (np.abs(self._funk - fun_next) <= self.funtol * 0.5 * (np.abs(self._funk) + np.abs(fun_next)))
+                | ((np.abs(self._funk) <= 1e-10) & (np.abs(fun_next) <= 1e-10))
+        )
+        if self.number_of_dummy_vars > 0:
+            have_dummies_converged = (
+                    (np.abs(w_hat[self.portfolio.number_of_assets:] - self.dummy_vars) <= self.wtol * 0.5 *
+                     (np.abs(w_hat[self.portfolio.number_of_assets:]) + np.abs(self.dummy_vars)))
+                    | ((np.abs(w_hat[self.portfolio.number_of_assets:]) < 1e-6) & (np.abs(self.dummy_vars) < 1e-6))
+            ).all()
+            self.dummy_vars = w_hat[self.portfolio.number_of_assets:]
+        else:
+            have_dummies_converged = True
+        if (has_w_converged and have_dummies_converged) or has_fun_converged:
+            # if verbose:
+            #     print(f"  Has func. converged: {has_fun_converged}; has w converged: {has_w_converged}")
             return False
         self.gamma = self.gamma * (1 - self.zeta * self.gamma)
         self._funk = fun_next
         return True
 
-    def solve(self):
+    def solve(self, verbose=True):
         i = 0
-        with tqdm(total=self.maxiter) as pbar:
-            while self.iterate() and i < self.maxiter:
-                i += 1
-                pbar.update()
-
+        iterator = range(self.maxiter)
+        if verbose:
+            iterator = tqdm(iterator)
+        for _ in iterator:
+            if not self.iterate(verbose=verbose):
+                break
+            i += 1
 
 def project_line_and_box(weights, lower_bound, upper_bound):
     def objective_function(variable, weights):