add risk parity loss function & risk budgeting allocator (#98)

Co-authored-by: Jan Krepl <[email protected]>

deepdow/layers/__init__.py

@@ -2,7 +2,8 @@
 
 from .collapse import (AttentionCollapse, AverageCollapse, ElementCollapse, ExponentialCollapse,
                        MaxCollapse, SumCollapse)
-from .allocate import (AnalyticalMarkowitz, NCO, NumericalMarkowitz, Resample, SoftmaxAllocator,
+from .allocate import (AnalyticalMarkowitz, NCO, NumericalMarkowitz,
+                       NumericalRiskBudgeting, Resample, SoftmaxAllocator,
                        SparsemaxAllocator, WeightNorm)
 from .misc import Cov2Corr, CovarianceMatrix, KMeans, MultiplyByConstant
 from .transform import Conv, RNN, Warp, Zoom
@@ -20,6 +21,7 @@
            'MultiplyByConstant',
            'NCO',
            'NumericalMarkowitz',
+           'NumericalRiskBudgeting',
            'Resample',
            'RNN',
            'SoftmaxAllocator',

deepdow/layers/allocate.py

@@ -545,3 +545,64 @@ def forward(self, x):
         normalized = clamped / clamped.sum()
 
         return torch.stack(n_samples * [normalized], dim=0)
+
+
+class NumericalRiskBudgeting(nn.Module):
+    """Convex optimization layer stylized into portfolio optimization problem.
+
+    Parameters
+    ----------
+    n_assets : int
+        Number of assets.
+
+    Attributes
+    ----------
+    cvxpylayer : CvxpyLayer
+        Custom layer used by a third party package called cvxpylayers.
+
+    References
+    ----------
+    [1] https://github.com/cvxgrp/cvxpylayers
+    [2] https://papers.ssrn.com/sol3/papers.cfm?abstract_id=2297383
+    [3] https://mpra.ub.uni-muenchen.de/37749/2/MPRA_paper_37749.pdf
+    """
+
+    def __init__(self, n_assets, max_weight=1):
+        """Construct."""
+        super().__init__()
+        covmat_sqrt = cp.Parameter((n_assets, n_assets))
+        b = cp.Parameter(n_assets, nonneg=True)
+
+        w = cp.Variable(n_assets)
+
+        term_1 = 0.5 * cp.sum_squares(covmat_sqrt @ w)
+        term_2 = b @ cp.log(w)
+
+        objective = cp.Minimize(term_1 - term_2)  # refer [2]
+        constraint = [cp.sum(w) == 1, w >= 0, w <= max_weight]  # refer [2]
+
+        prob = cp.Problem(objective, constraint)
+
+        assert prob.is_dpp()
+
+        self.cvxpylayer = CvxpyLayer(prob, parameters=[covmat_sqrt, b], variables=[w])
+
+    def forward(self, covmat_sqrt, b):
+        """Perform forward pass.
+
+        Parameters
+        ----------
+        covmat : torch.Tensor
+            Of shape (n_samples, n_assets, n_assets) representing the covariance matrix.
+
+        b : torch.Tensor
+            Of shape (n_samples, n_assets) representing the budget,
+            risk contribution from each component (asset) is equal to the budget, refer [3]
+
+        Returns
+        -------
+        weights : torch.Tensor
+            Of shape (n_samples, n_assets) representing the optimal weights as determined by the convex optimizer.
+
+        """
+        return self.cvxpylayer(covmat_sqrt, b)[0]

deepdow/losses.py

@@ -3,6 +3,7 @@
 All losses are designed for minimization.
 """
 from types import MethodType
+from .layers import CovarianceMatrix
 
 import torch
 
@@ -1038,3 +1039,64 @@ def __repr__(self):
                                                                                   self.returns_channel,
                                                                                   self.input_type,
                                                                                   self.output_type)
+
+
+class RiskParity(Loss):
+    """Risk Parity Portfolio.
+
+    Parameters
+    ----------
+    returns_channel : int
+        Which channel of the `y` target represents returns.
+
+    Attributes
+    ----------
+    covariance_layer : deepdow.layers.CoverianceMatrix
+        Covarioance matrix layer.
+
+    References
+    ----------
+    https://papers.ssrn.com/sol3/papers.cfm?abstract_id=2297383
+    """
+
+    def __init__(self, returns_channel=0):
+        self.returns_channel = returns_channel
+        self.covariance_layer = CovarianceMatrix(sqrt=False)
+
+    def __call__(self, weights, y):
+        """Compute loss.
+
+        Parameters
+        ----------
+        weights : torch.Tensor
+            Tensor of shape `(n_samples, n_assets)` representing the predicted
+            weights by our portfolio optimizer.
+
+        y : torch.Tensor
+            Tensor of shape `(n_samples, n_channels, horizon, n_assets)`
+            representing the evolution over the next `horizon` timesteps.
+
+        Returns
+        -------
+         torch.Tensor
+            Tensor of shape `(n_samples,)` representing the per sample risk parity.
+        """
+        n_assets = weights.shape[-1]
+        covar = self.covariance_layer(y[:, self.returns_channel, ...])  # (n_samples, n_assets, n_assets)
+
+        weights = weights.unsqueeze(dim=1)
+        volatility = torch.sqrt(torch.matmul(weights,
+                                             torch.matmul(covar,
+                                                          weights.permute((0, 2, 1)))))  # (n_samples, 1, 1)
+        c = (covar * weights) / volatility  # (n_samples, n_assets, n_assets)
+        risk = volatility / n_assets  # (n_samples, 1, 1)
+
+        budget = torch.matmul(weights, c)  # (n_samples, n_assets, n_assets)
+        rp = torch.sum((risk - budget)**2, dim=-1).view(-1)  # (n_samples,)
+
+        return rp
+
+    def __repr__(self):
+        """Generate representation string."""
+        return "{}(returns_channel={})".format(self.__class__.__name__,
+                                               self.returns_channel)

docs/source/layers.rst

@@ -262,6 +262,32 @@ tensors (batched along the sample dimension):
 
 
 
+.. warning::
+
+    The major downside of using this allocator is a significant decrease in speed.
+
+NumericalRiskBudgeting
+**********************
+Proposed in [Spinu2013]_.
+
+.. math::
+
+    \begin{aligned}
+    \min_{\textbf{w}} \quad & \frac{1}{2}{\textbf{w}}^{T}  \boldsymbol{\Sigma} \textbf{w} - \sum_{i=1}^{N} b_i  \log(w_i) \\
+    \textrm{s.t.} \quad & \sum_{i=1}^{N}w_i = 1 \\
+    \quad & w_i >= 0, i \in \{1,...,N\}\\
+    \quad & w_i <= w_{\text{max}}, i \in \{1,...,N\}\\
+    \end{aligned}
+
+where the :math:`b_i, i=1,..,N` are the risk budgets. The user needs to provide
+:code:`n_assets` (:math:`N` in the above formulation) and :code:`max_weight`
+(:math:`w_{\text{max}}`) when constructing this layer. To perform a forward pass one passes the following
+tensors (batched along the sample dimension):
+
+- :code:`covmat_sqrt` - Corresponds to a (matrix) square root of the covariance matrix :math:`\boldsymbol{\Sigma}`
+- :code:`b` - Risk budgets
+
+
 .. warning::
 
     The major downside of using this allocator is a significant decrease in speed.
@@ -494,6 +520,9 @@ References
 .. [Prado2019]
    Lopez de Prado, M. (2019). A Robust Estimator of the Efficient Frontier. Available at SSRN 3469961.
 
+.. [Spinu2013]
+   Spinu, Florin, An Algorithm for Computing Risk Parity Weights (July 30, 2013). Available at SSRN: https://ssrn.com/abstract=2297383 or http://dx.doi.org/10.2139/ssrn.2297383
+
 .. [Jiang2017]
    Jiang, Zhengyao, and Jinjun Liang. "Cryptocurrency portfolio management with deep reinforcement learning." 2017 Intelligent Systems Conference (IntelliSys). IEEE, 2017
 
@@ -522,4 +551,4 @@ References
    Bodnar, Taras, Nestor Parolya, and Wolfgang Schmid. "On the equivalence of quadratic optimization problems commonly used in portfolio theory." European Journal of Operational Research 229.3 (2013): 637-644.
 
 .. [Jaderberg2015]
-   Jaderberg, Max, Karen Simonyan, and Andrew Zisserman. "Spatial transformer networks." Advances in neural information processing systems. 2015.
+   Jaderberg, Max, Karen Simonyan, and Andrew Zisserman. "Spatial transformer networks." Advances in neural information processing systems. 2015.

docs/source/losses.rst

@@ -168,6 +168,16 @@ The **negative** of mean portfolio returns over the :code:`horizon` time steps.
 
     {\mu}^{\textbf{w}} = \frac{\sum_{i}^{\text{horizon}} r^{\textbf{w}}_{i} }{\text{horizon}}
 
+RiskParity
+**********
+
+.. math::
+
+   \sum_{i=1}^{N}\Big(\frac{\sigma}{N} - w_i \big(\frac{\Sigma\textbf{w}}{\sigma}\big)_i\Big) ^ 2
+
+where :math:`\sigma=\sqrt{\textbf{w}^T\Sigma\textbf{w}}` and :math:`\Sigma` is
+the covariance matrix of asset returns.
+
 Quantile (Value at Risk)
 ************************
 The **negative** of the :code:`p`-quantile of portfolio returns. Note that in the background it solved via

tests/test_layers.py

@@ -4,7 +4,8 @@
 from deepdow.layers import (AverageCollapse, AttentionCollapse, ElementCollapse,
                             ExponentialCollapse, MaxCollapse,
                             SumCollapse)
-from deepdow.layers import (AnalyticalMarkowitz, NCO, NumericalMarkowitz, Resample,
+from deepdow.layers import (AnalyticalMarkowitz, NCO, NumericalMarkowitz,
+                            NumericalRiskBudgeting, Resample,
                             SoftmaxAllocator, SparsemaxAllocator, WeightNorm)
 from deepdow.layers import Cov2Corr, CovarianceMatrix, KMeans, MultiplyByConstant
 from deepdow.layers import Conv, RNN, Warp, Zoom
@@ -326,6 +327,9 @@ def test_basic(self, Xy_dummy):
         weights = popt(rets, covmat_sqrt, gamma, alpha)
 
         assert weights.shape == (n_samples, n_assets)
+        assert torch.allclose(weights.sum(dim=-1), torch.ones(n_samples,
+                                                              device=device,
+                                                              dtype=dtype))
         assert weights.dtype == X.dtype
         assert weights.device == X.device
 
@@ -787,3 +791,30 @@ def test_equality_with_warp(self):
         x_warped = layer_warp(X, tform)
 
         assert torch.allclose(x_zoomed, x_warped)
+
+
+class TestNumericalRiskBudgeting:
+    def test_basic(self, Xy_dummy):
+        X, _, _, _ = Xy_dummy
+        device, dtype = X.device, X.dtype
+        n_samples, n_channels, lookback, n_assets = X.shape
+
+        popt = NumericalRiskBudgeting(n_assets)
+
+        covmat_sqrt__ = torch.rand((n_assets, n_assets)).to(device=X.device, dtype=X.dtype)
+        covmat_sqrt_ = covmat_sqrt__ @ covmat_sqrt__
+        covmat_sqrt_.add_(torch.eye(n_assets, dtype=dtype, device=device))
+
+        covmat_sqrt = torch.stack(n_samples * [covmat_sqrt_])
+
+        budgets = torch.rand((n_samples, n_assets), dtype=dtype, device=device)
+        budgets /= budgets.sum(dim=-1, keepdim=True)
+
+        weights = popt(covmat_sqrt, budgets)
+
+        assert weights.shape == (n_samples, n_assets)
+        assert torch.allclose(weights.sum(dim=-1), torch.ones(n_samples,
+                                                              device=device,
+                                                              dtype=dtype))
+        assert weights.dtype == X.dtype
+        assert weights.device == X.device

tests/test_losses.py

@@ -3,13 +3,18 @@
 import pytest
 import torch
 
-from deepdow.losses import (Alpha, CumulativeReturn, LargestWeight, Loss, MaximumDrawdown, MeanReturns, Quantile,
-                            SharpeRatio, Softmax, SortinoRatio, SquaredWeights, StandardDeviation, TargetMeanReturn,
-                            TargetStandardDeviation, WorstReturn, log2simple, portfolio_returns,
-                            portfolio_cumulative_returns, simple2log)
+from deepdow.losses import (Alpha, CumulativeReturn, LargestWeight, Loss,
+                            MaximumDrawdown, MeanReturns, RiskParity, Quantile,
+                            SharpeRatio, Softmax, SortinoRatio, SquaredWeights,
+                            StandardDeviation, TargetMeanReturn,
+                            TargetStandardDeviation, WorstReturn, log2simple,
+                            portfolio_returns, portfolio_cumulative_returns,
+                            simple2log)
 
-ALL_LOSSES = [Alpha, CumulativeReturn, LargestWeight, MaximumDrawdown, MeanReturns, Quantile, SharpeRatio, Softmax,
-              SortinoRatio, SquaredWeights, StandardDeviation, TargetMeanReturn, TargetStandardDeviation, WorstReturn]
+ALL_LOSSES = [Alpha, CumulativeReturn, LargestWeight, MaximumDrawdown,
+              MeanReturns, RiskParity, Quantile, SharpeRatio, Softmax,
+              SortinoRatio, SquaredWeights, StandardDeviation,
+              TargetMeanReturn, TargetStandardDeviation, WorstReturn]
 
 
 class TestHelpers:
@@ -351,3 +356,44 @@ def test_no_drawdowns(self):
         loss = loss_inst(w, y)[0]
 
         assert loss == 0
+
+
+class TestRiskParity:
+    @staticmethod
+    def stupid_compute(w, y):
+        """Straightforward implementation with list comprehensions."""
+        from deepdow.layers import CovarianceMatrix
+
+        n_samples, n_assets = w.shape
+        covar = CovarianceMatrix(sqrt=False)(y[:, 0, ...])  # returns_channel=0
+
+        var = torch.cat([(w[[i]] @ covar[i]) @ w[[i]].permute(1, 0) for i in range(n_samples)], dim=0)
+        vol = torch.sqrt(var)
+
+        lhs = vol / n_assets
+        rhs = torch.cat([(1 / vol[i]) * w[[i]] * (w[[i]] @ covar[i]) for i in range(n_samples)], dim=0)
+
+        res = torch.tensor([((lhs[i] - rhs[i]) ** 2).sum() for i in range(n_samples)])
+
+        return res
+
+    def test_correct_fpass(self):
+        device, dtype = torch.device("cpu"), torch.float32
+        n_samples, n_channels, horizon, n_assets = 2, 3, 10, 5
+
+        # Generate weights and targets
+        torch.manual_seed(2)
+        y = torch.rand((n_samples, n_channels, horizon, n_assets),
+                       dtype=dtype,
+                       device=device)
+
+        weights = torch.rand((n_samples, n_assets),
+                             dtype=dtype,
+                             device=device)
+
+        weights /= weights.sum(dim=-1, keepdim=True)
+
+        res_stupid = self.stupid_compute(weights, y)
+        res_actual = RiskParity()(weights, y)
+
+        assert torch.allclose(res_stupid, res_actual)