#!/usr/bin/env python3
# coding: utf-8
# @Author: ArthurBernard
# @Email: arthur.bernard.92@gmail.com
# @Date: 2019-09-12 14:52:08
# @Last modified by: ArthurBernard
# @Last modified time: 2019-11-05 17:20:52
""" Algorithms of portfolio allocation.
Risk-based and mean-variance methods that compute portfolio weights
from a sample of asset returns. Each function takes a 2-D array where
columns are asset price/return series and returns the optimal weights
as a 1-D array summing to one. All of them accept an optional ``cov=``
callable that replaces the sample-covariance estimation — see
:mod:`fynance.portfolio.covariance` for ready-made estimators.
A walk-forward wrapper, :func:`rolling_allocation`, applies any of
these methods on a rolling training window — useful for backtesting
allocation strategies without lookahead bias; a ``cov=`` estimator
reaches the allocator unchanged through its ``**kwargs`` forwarding.
Main entry points
-----------------
- :func:`ERC` — Equal Risk Contribution (risk-parity).
- :func:`RBP` — Risk Budgeting Portfolio (ERC generalized to arbitrary
per-asset risk budgets).
- :func:`HRP` — Hierarchical Risk Parity.
- :func:`IVP` — Inverse Variance Portfolio.
- :func:`MDP` — Maximum Diversified Portfolio.
- :func:`MVP`, :func:`MVP_uc` — Minimum Variance Portfolio (constrained
/ unconstrained).
- :func:`rolling_allocation` — walk-forward wrapper.
- :mod:`fynance.portfolio.covariance` — covariance estimators (sample,
Ledoit-Wolf shrinkage, EWMA, factor model, Marchenko-Pastur denoising)
usable as the ``cov=`` parameter of the allocators above.
"""
from __future__ import annotations
# Built-in packages
import warnings
from typing import Callable
# Third-party
import numpy as np
import scipy.cluster.hierarchy as sch
from numpy.typing import NDArray
from scipy.optimize import Bounds, LinearConstraint, minimize
# Local packages
from fynance.metrics import diversified_ratio
__all__ = ['ERC', 'HRP', 'IVP', 'MDP', 'MVP', 'MVP_uc', 'RBP', 'rolling_allocation']
# =========================================================================== #
# covariance estimation seam #
# =========================================================================== #
def _estimate_cov(
X: NDArray[np.float64],
cov: Callable[[NDArray[np.float64]], NDArray[np.float64]] | None,
) -> NDArray[np.float64]:
r""" Estimate the covariance matrix used by an allocator.
Opt-in seam shared by all the allocators below: when ``cov`` is
``None`` this is a thin, behavior-preserving wrapper around
:func:`numpy.cov` (the historical estimator); when ``cov`` is a
callable it is called on ``X`` and the returned matrix is validated
(square, matching ``X``'s number of columns, symmetric within a tight
tolerance) before being handed back to the caller.
Parameters
----------
X : array_like
Returns panel, shape ``(T, N)``.
cov : callable or None
Callable mapping the ``(T, N)`` array ``X`` to an ``(N, N)``
covariance matrix, e.g.
:func:`fynance.portfolio.covariance.ledoit_wolf`. If ``None``
(default), the sample covariance (:func:`numpy.cov`) is used.
Returns
-------
np.ndarray
Symmetric ``(N, N)`` covariance matrix.
Raises
------
ValueError
If ``cov`` is a callable and its return value is not square, does
not match ``X``'s number of columns, or is not symmetric within
``1e-8``.
"""
if cov is None:
return np.atleast_2d(np.cov(X, rowvar=False))
N = np.asarray(X).shape[1]
sigma = np.atleast_2d(np.asarray(cov(X), dtype=np.float64))
if sigma.shape != (N, N):
raise ValueError(
f"cov callable must return an ({N}, {N}) covariance matrix "
f"matching X's {N} columns, got shape {sigma.shape}."
)
if not np.allclose(sigma, sigma.T, atol=1e-8):
raise ValueError(
f"cov callable returned a non-symmetric {sigma.shape} matrix "
f"(max |sigma - sigma.T| = {np.max(np.abs(sigma - sigma.T)):.3e})."
)
return sigma
# =========================================================================== #
# Equal Risk Contribution #
# =========================================================================== #
[docs]
def ERC(
X: NDArray[np.float64],
w0: NDArray[np.float64] | None = None,
up_bound: float = 1.,
low_bound: float = 0.,
cov: Callable[[NDArray[np.float64]], NDArray[np.float64]] | None = None,
) -> NDArray[np.float64]:
r""" Get weights of Equal Risk Contribution portfolio allocation.
Risk-parity allocation: each asset contributes the same amount of
total portfolio variance. ERC sits between the naive 1/N and the
minimum-variance portfolio — it requires only the covariance
matrix and is robust to noisy expected returns, which makes it a
common choice when return forecasts are unreliable.
The optimizer (SLSQP) minimizes a smooth surrogate of the
risk-contribution dispersion under sum-to-one and box constraints.
Notes
-----
Weights of Equal Risk Contribution, as described by S. Maillard, T.
Roncalli and J. Teiletche [1]_, verify the following problem:
.. math::
w = \text{arg min } f(w) \\
u.c. \begin{cases}w'e = 1 \\
0 \leq w_i \leq 1 \\
\end{cases}
With:
.. math::
f(w) = N \sum_{i=1}^{N}w_i^2 (\Omega w)_i^2
- \sum_{i,j=1}^{N} w_i w_j (\Omega w)_i (\Omega w)_j
Where :math:`\Omega` is the variance-covariance matrix of `X` and :math:`N`
the number of assets.
Parameters
----------
X : array_like
Each column is a series of price or return's asset.
w0 : array_like, optional
Initial weights for the optimizer.
up_bound, low_bound : float, optional
Respectively maximum and minimum values of weights, such that low_bound
:math:`\leq w_i \leq` up_bound :math:`\forall i`. Default is 0 and 1.
cov : callable, optional
Callable mapping the ``(T, N)`` training array to an ``(N, N)``
covariance matrix, e.g.
:func:`fynance.portfolio.covariance.ledoit_wolf`; default `None`
keeps the sample covariance.
Returns
-------
array_like
Weights that minimize the Equal Risk Contribution portfolio.
References
----------
.. [1] http://thierry-roncalli.com/download/erc-slides.pdf
Examples
--------
>>> import numpy as np
>>> from fynance.portfolio.covariance import ledoit_wolf
>>> rng = np.random.default_rng(0)
>>> X = rng.normal(0.0, 0.01, size=(200, 4))
>>> w = ERC(X, cov=ledoit_wolf)
>>> w.shape
(4, 1)
>>> bool(np.isclose(w.sum(), 1.0))
True
"""
T, N = X.shape
SIGMA = _estimate_cov(X, cov)
if N == 1:
return np.ones([1, 1])
up_bound = max(up_bound, 1 / N)
# Clamp low_bound so the box stays compatible with the sum-to-one
# constraint: with low_bound > 1/N every feasible weight already exceeds
# its share, the constraint is infeasible and SLSQP would silently return
# weights summing to more than one (as HRP/IVP already guard against).
low_bound = min(low_bound, 1 / N)
# The risk-contribution surrogate is quartic in the covariance, so on
# return-scale inputs (values ~1e-4) it collapses to ~1e-16, far below
# SLSQP's default ftol and the optimizer stops at the 1/N start. Rescale
# the covariance to unit trace so the objective is O(1) regardless of the
# input scale; this leaves the argmin (the weights) unchanged.
scale = np.trace(SIGMA) / N
if scale <= 0:
return np.ones([N, 1]) / N
SIGMA = SIGMA / scale
def f_ERC(w):
w = w.reshape([N, 1])
arg = N * np.sum(w ** 2 * (SIGMA @ w) ** 2)
return arg - np.sum(w * (SIGMA @ w) * np.sum(w * (SIGMA @ w)))
# Set inital weights
if w0 is None:
w0 = np.ones([N]) / N
const_sum = LinearConstraint(np.ones([1, N]), [1], [1])
const_ind = Bounds(low_bound * np.ones([N]), up_bound * np.ones([N]))
result = minimize(
f_ERC,
w0,
method='SLSQP',
constraints=[const_sum],
bounds=const_ind,
options={'ftol': 1e-12, 'maxiter': 1000},
)
return result.x.reshape([N, 1])
# =========================================================================== #
# Risk Budgeting Portfolio #
# =========================================================================== #
def _validate_budgets(budgets: NDArray[np.float64] | None, N: int) -> NDArray[np.float64]:
""" Validate and normalize a risk-budget vector.
Parameters
----------
budgets : array_like or None
Candidate risk budgets, expected length ``N``. ``None`` means the
equal-budget default ``1 / N``.
N : int
Number of assets.
Returns
-------
np.ndarray
Length-``N`` budget vector, strictly positive and summing to
exactly 1.
Raises
------
ValueError
If `budgets` does not have length `N`, contains a non-positive
entry, or its sum deviates from 1 by more than ``1e-8``.
"""
if budgets is None:
return np.full(N, 1.0 / N)
b = np.asarray(budgets, dtype=np.float64).ravel()
if b.shape != (N,):
raise ValueError(
f"budgets must be a length-{N} vector, got shape {b.shape}."
)
if np.any(b <= 0):
raise ValueError("budgets entries must be strictly positive.")
total = b.sum()
if abs(total - 1.0) > 1e-8:
raise ValueError(
f"budgets must sum to 1 (within 1e-8), got sum={total!r}."
)
# Silently renormalize away any float rounding within the tolerance.
return b / total
[docs]
def RBP(
X: NDArray[np.float64],
budgets: NDArray[np.float64] | None = None,
w0: NDArray[np.float64] | None = None,
up_bound: float = 1.,
low_bound: float = 0.,
cov: Callable[[NDArray[np.float64]], NDArray[np.float64]] | None = None,
) -> NDArray[np.float64]:
r""" Get weights of a Risk Budgeting Portfolio allocation.
Generalizes :func:`ERC` to arbitrary per-asset risk budgets: instead of
equalizing risk contributions, each asset ``i`` is allocated a target
share ``b_i`` of total portfolio variance, with :math:`\sum_i b_i = 1`.
Passing ``budgets=None`` falls back to the equal-budget case
:math:`b_i = 1/N \, \forall i`, which reproduces :func:`ERC`.
The optimizer (SLSQP) minimizes a smooth least-squares surrogate of the
gap between each asset's risk contribution and its target budget, under
sum-to-one and box constraints.
Notes
-----
Weights of the Risk Budgeting Portfolio, as described by T. Roncalli
[6]_, verify the following problem:
.. math::
w = \text{arg min } f(w) \\
u.c. \begin{cases}w'e = 1 \\
0 \leq w_i \leq 1 \\
\end{cases}
With:
.. math::
f(w) = \sum_{i=1}^{N} \left( w_i (\Omega w)_i
- b_i \, w' \Omega w \right)^2
Where :math:`\Omega` is the variance-covariance matrix of `X`, :math:`N`
the number of assets and :math:`b` the target risk-budget vector (with
:math:`\sum_{i=1}^{N} b_i = 1`). With :math:`b_i = 1/N \, \forall i` this
objective shares the same minimizers as :func:`ERC`'s (both vanish
exactly when every asset's risk contribution matches its budget).
Parameters
----------
X : array_like
Each column is a series of price or return's asset.
budgets : array_like, optional
Target risk budget per asset, length ``N``, strictly positive
entries summing to 1 (within ``1e-8``, silently renormalized inside
that tolerance). Default `None` spreads the budget equally
(:math:`b_i = 1/N`), reproducing :func:`ERC`.
w0 : array_like, optional
Initial weights for the optimizer.
up_bound, low_bound : float, optional
Respectively maximum and minimum values of weights, such that low_bound
:math:`\leq w_i \leq` up_bound :math:`\forall i`. Default is 0 and 1.
cov : callable, optional
Callable mapping the ``(T, N)`` training array to an ``(N, N)``
covariance matrix, e.g.
:func:`fynance.portfolio.covariance.ledoit_wolf`; default `None`
keeps the sample covariance.
Returns
-------
array_like
Weights whose risk contributions match `budgets`.
Raises
------
ValueError
If `budgets` does not have length `N`, contains a non-positive
entry, or its sum deviates from 1 by more than ``1e-8``.
References
----------
.. [6] T. Roncalli, "Introduction to Risk Parity and Budgeting", 2013,
https://arxiv.org/abs/1403.1889
See Also
--------
ERC
fynance.portfolio.attribution.risk_contribution
Examples
--------
Two independent assets with volatilities ``sigma_1=0.01`` and
``sigma_2=0.03``: for a diagonal covariance the risk-budgeting
first-order condition reduces to :math:`w_i \sigma_i \propto
\sqrt{b_i}`, i.e. :math:`w_i \propto \sqrt{b_i} / \sigma_i`.
>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> sigma1, sigma2 = 0.01, 0.03
>>> X = np.column_stack([
... rng.normal(0.0, sigma1, 2000),
... rng.normal(0.0, sigma2, 2000),
... ])
>>> b = np.array([0.8, 0.2])
>>> w = RBP(X, budgets=b)
>>> w.shape
(2, 1)
>>> bool(np.isclose(w.sum(), 1.0))
True
>>> expected = np.sqrt(b) / np.array([sigma1, sigma2])
>>> expected = expected / expected.sum()
>>> bool(np.allclose(w.flatten(), expected, atol=1e-2))
True
"""
T, N = X.shape
b = _validate_budgets(budgets, N)
SIGMA = _estimate_cov(X, cov)
if N == 1:
return np.ones([1, 1])
up_bound = max(up_bound, 1 / N)
# Clamp low_bound so the box stays compatible with the sum-to-one
# constraint: with low_bound > 1/N every feasible weight already exceeds
# its share, the constraint is infeasible and SLSQP would silently return
# weights summing to more than one (as HRP/IVP already guard against).
low_bound = min(low_bound, 1 / N)
# Same scale-safety rationale as ERC: rescale the covariance to unit
# trace so the quartic objective is O(1) regardless of the input scale;
# this leaves the argmin (the weights) unchanged.
scale = np.trace(SIGMA) / N
if scale <= 0:
return np.ones([N, 1]) / N
SIGMA = SIGMA / scale
b_col = b.reshape([N, 1])
def f_RBP(w):
w = w.reshape([N, 1])
sigma_w = SIGMA @ w
port_var = w.T @ sigma_w
return np.sum((w * sigma_w - b_col * port_var) ** 2)
# Set inital weights
if w0 is None:
w0 = np.ones([N]) / N
const_sum = LinearConstraint(np.ones([1, N]), [1], [1])
const_ind = Bounds(low_bound * np.ones([N]), up_bound * np.ones([N]))
result = minimize(
f_RBP,
w0,
method='SLSQP',
constraints=[const_sum],
bounds=const_ind,
options={'ftol': 1e-12, 'maxiter': 1000},
)
return result.x.reshape([N, 1])
# =========================================================================== #
# HRP developed by Marcos Lopez de Prado #
# =========================================================================== #
def _get_quasi_diag(link: NDArray[np.float64]) -> list[int]:
""" Compute quasi diagonal matrix.
Parameter
---------
link: list of N lists
Linkage matrix, N list (cluster) of 4-tuple such that the two first
elements are the costituents, third report the distance between the two
first, and fourth is the number of element (<= N) in this cluster.
Returns
-------
sortIx: list
Sorted list of items.
"""
link = link.astype(int)
numItems = link[-1, 3] # number of original leaf items
# Use a plain Python list for cheap appends/inserts
items = [int(link[-1, 0]), int(link[-1, 1])]
while max(items) >= numItems:
expanded = []
for item in items:
if item >= numItems:
cluster_idx = item - numItems
expanded.append(int(link[cluster_idx, 0]))
expanded.append(int(link[cluster_idx, 1]))
else:
expanded.append(item)
items = expanded
return items
def _get_rec_bisec(mat_cov: NDArray[np.float64], sortIx: list[int]) -> NDArray[np.float64]:
""" Compute weights via recursive bisection.
Parameters
----------
mat_cov: np.ndarray
Matrix variance-covariance (N x N).
sortIx: list or np.ndarray of int
Sorted list of asset indices (0..N-1).
Returns
-------
np.ndarray
Weight vector of shape (N,) indexed by sortIx order.
"""
n = len(sortIx)
w = np.ones(n)
cItems = [list(range(n))]
while len(cItems) > 0:
cItems = [i[j: k] for i in cItems for j, k in (
(0, int(len(i) / 2)),
(int(len(i) / 2), len(i))
) if len(i) > 1]
for i in range(0, len(cItems), 2):
cItems0_idx = cItems[i]
cItems1_idx = cItems[i + 1]
cItems0 = [sortIx[j] for j in cItems0_idx]
cItems1 = [sortIx[j] for j in cItems1_idx]
cVar0 = _get_cluster(mat_cov, cItems0)
cVar1 = _get_cluster(mat_cov, cItems1)
alpha = 1 - cVar0 / (cVar0 + cVar1)
w[cItems0_idx] *= alpha
w[cItems1_idx] *= 1 - alpha
return w
def _get_cluster(mat_cov: NDArray[np.float64], cItems: list[int]) -> float:
""" Compute cluster for variance.
Parameters
----------
mat_cov: np.ndarray
Covariance matrix.
cItems: list or np.ndarray of int
Cluster asset indices.
Returns
-------
cVar: float
Cluster variance
"""
cov_ = mat_cov[np.ix_(cItems, cItems)]
w_ = _get_IVP(cov_).reshape(-1, 1)
cVar = (w_.T @ cov_) @ w_
return float(cVar.item())
def _get_IVP(mat_cov: NDArray[np.float64]) -> NDArray[np.float64]:
""" Compute the inverse-variance matrix.
Parameters
----------
mat_cov : NDArray[np.float64]
Variance-covariance matrix.
Returns
-------
NDArray[np.float64]
Inverse-variance weights.
"""
ivp = 1. / np.diag(mat_cov)
ivp /= np.sum(ivp)
return ivp
[docs]
def HRP(
X: NDArray[np.float64],
method: str = 'single',
metric: str = 'euclidean',
low_bound: float = 0.,
up_bound: float = 1.0,
cov: Callable[[NDArray[np.float64]], NDArray[np.float64]] | None = None,
) -> NDArray[np.float64]:
r""" Get weights of the Hierarchical Risk Parity allocation.
Cluster-based allocation that avoids inverting the full
covariance matrix. Compared with classical Markowitz solutions,
HRP is far more stable when ``N`` is large or assets are highly
correlated, because it groups similar assets first and only
allocates risk *within* and *between* clusters.
Three steps: (1) build a correlation distance and run hierarchical
linkage, (2) reorder the covariance matrix into quasi-diagonal
form, (3) recursively bisect the ordered tree, allocating weights
by inverse variance to each branch.
Notes
-----
Hierarchical Risk Parity algorithm is developed by Marco Lopez de Prado
[2]_. First step is clustering and second step is allocating weights.
Parameters
----------
X : array_like
Each column is a price or return's asset series. Some errors will
happen if one or more series are constant.
method, metric: str
Parameters for linkage algorithm, default ``method='single'`` and
``metric='euclidean'``.
low_bound, up_bound : float
Respectively minimum and maximum value of weights, such that low_bound
:math:`\leq w_i \leq` up_bound :math:`\forall i`. Default is 0 and 1.
cov : callable, optional
Callable mapping the ``(T, N)`` training array to an ``(N, N)``
covariance matrix, e.g.
:func:`fynance.portfolio.covariance.ledoit_wolf`; default `None`
keeps the sample covariance.
Returns
-------
np.ndarray
Vector of weights computed by HRP algorithm.
References
----------
.. [2] https://ssrn.com/abstract=2708678
"""
X = np.asarray(X, dtype=np.float64)
T, N = X.shape
if N == 1:
return np.ones([1, 1])
up_bound = max(up_bound, 1.0 / N)
low_bound = min(low_bound, 1.0 / N)
mat_cov = _estimate_cov(X, cov)
diag_cov = np.sqrt(np.diag(mat_cov))
outer_diag = np.outer(diag_cov, diag_cov)
with np.errstate(invalid='ignore', divide='ignore'):
mat_corr = np.divide(mat_cov, outer_diag)
mat_corr = np.clip(np.nan_to_num(mat_corr, nan=0.0, posinf=0.0, neginf=0.0), -1.0, 1.0)
mat_dist = ((1.0 - mat_corr) / 2.0) ** 0.5
mat_dist_upper = mat_dist[np.triu_indices(N, k=1)]
link = sch.linkage(mat_dist_upper, method=method, metric=metric)
sortIx = _get_quasi_diag(link)
w_sorted = _get_rec_bisec(mat_cov, sortIx)
# Scatter the cluster-ordered weights back to original asset order.
w = np.empty(N)
w[np.asarray(sortIx)] = w_sorted
return _normalize(w.reshape([N, 1]), up_bound=up_bound, low_bound=low_bound)
# =========================================================================== #
# Inverse Variance Portfolio #
# =========================================================================== #
[docs]
def IVP(
X: NDArray[np.float64],
normalize: bool = False,
low_bound: float = 0.,
up_bound: float = 1.0,
cov: Callable[[NDArray[np.float64]], NDArray[np.float64]] | None = None,
) -> NDArray[np.float64]:
r""" Get weights of the Inverse Variance Portfolio allocation.
Notes
-----
w are computed by the inverse of the asset's variance [3]_ such that:
.. math::
w_i = \frac{1}{\sigma_k^2} (\sum_{i} \frac{1}{\sigma_i^2})^{-1}
With :math:`\sigma_i^2` is the variance of asset i.
Parameters
----------
X : array_like
Each column is a price or return's asset series.
normalize : bool, optional
If True clip the weights to the ``[low_bound, up_bound]`` box and
renormalize them to sum to one (via :func:`_normalize`), preserving
the inverse-variance ordering. If False (default) return the raw
inverse-variance weights :math:`w_i \propto 1 / \sigma_i^2`, which
already sum to one and lie in :math:`[0, 1]`.
low_bound, up_bound : float, optional
Respectively minimum and maximum values of weights, such that low_bound
:math:`\leq w_i \leq` up_bound :math:`\forall i`. Default is 0 and 1.
cov : callable, optional
Callable mapping the ``(T, N)`` training array to an ``(N, N)``
covariance matrix, e.g.
:func:`fynance.portfolio.covariance.ledoit_wolf`; default `None`
keeps the sample covariance.
Returns
-------
np.ndarray
Vector of weights computed by the IVP algorithm.
References
----------
.. [3] https://en.wikipedia.org/wiki/Inverse-variance_weighting
"""
mat_cov = _estimate_cov(X, cov)
N = mat_cov.shape[0]
if N == 1:
return np.ones([1, 1])
w = _get_IVP(mat_cov)
if normalize:
up_bound = max(up_bound, 1 / N)
low_bound = min(low_bound, 1 / N)
w = _normalize(w, up_bound=up_bound, low_bound=low_bound)
return w.reshape([N, 1])
# =========================================================================== #
# Minimum Variance Portfolio #
# =========================================================================== #
[docs]
def MVP(
X: NDArray[np.float64],
normalize: bool = False,
cov: Callable[[NDArray[np.float64]], NDArray[np.float64]] | None = None,
) -> NDArray[np.float64]:
r""" Get weights of the Minimum Variance Portfolio allocation.
Closed-form Markowitz allocation that minimizes the portfolio
variance subject only to a sum-to-one constraint. Weights are not
constrained to be positive — short positions are allowed and the
weights returned by the analytical formula can be negative or
larger than one. Use :func:`MVP_uc` for a constrained variant.
The covariance matrix must be invertible; pseudo-inverse is used
as a fallback when ``X`` has linearly dependent columns.
Notes
-----
The vector of weights noted :math:`w` that minimize the portfolio variance
[4]_ is define as below:
.. math:: w = \frac{\Omega^{-1} e}{e' \Omega^{-1} e} \\
.. math:: \text{With } \sum_{i=1}^{N} w_i = 1
Where :math:`\Omega` is the asset's variance-covariance matrix and
:math:`e` is a vector of ones.
Parameters
----------
X : array_like
Each column is a time-series of price or return's asset.
normalize : boolean, optional
If True normalize the weigths such that :math:`0 \leq w_i \leq 1` and
:math:`\sum_{i=1}^{N} w_i = 1`, :math:`\forall i`. Default is False.
cov : callable, optional
Callable mapping the ``(T, N)`` training array to an ``(N, N)``
covariance matrix, e.g.
:func:`fynance.portfolio.covariance.ledoit_wolf`; default `None`
keeps the sample covariance.
Returns
-------
array_like
Vector of weights to apply to the assets.
References
----------
.. [4] https://breakingdownfinance.com/finance-topics/modern-portfolio-theory/minimum-variance-portfolio/
See Also
--------
HRP
"""
mat_cov = _estimate_cov(X, cov)
if mat_cov.shape[0] == 1:
return np.ones([1, 1])
# Inverse variance matrix
try:
iv = np.linalg.inv(mat_cov)
except np.linalg.LinAlgError:
iv = np.linalg.pinv(mat_cov)
e = np.ones([iv.shape[0], 1])
w = (iv @ e) / (e.T @ iv @ e)
if normalize:
w = w - np.min(w)
return w / np.sum(w)
return w
[docs]
def MVP_uc(
X: NDArray[np.float64],
w0: NDArray[np.float64] | None = None,
up_bound: float = 1.,
low_bound: float = 0.,
cov: Callable[[NDArray[np.float64]], NDArray[np.float64]] | None = None,
) -> NDArray[np.float64]:
r""" Get weights of the Minimum Variance Portfolio under constraints.
Numerical (SLSQP) Markowitz allocation that minimizes portfolio
variance subject to box constraints on each weight, in addition
to the sum-to-one constraint. Use this variant whenever short
selling must be excluded or per-asset caps need to be enforced;
use :func:`MVP` for the unconstrained closed-form solution.
Notes
-----
Weights of Minimum Variance Portfolio verify the following problem:
.. math::
w = \text{arg min } w' \Omega w \\
u.c. \begin{cases}w'e = 1 \\
0 \leq w_i \leq 1 \\
\end{cases}
Where :math:`\Omega` is the variance-covariance matrix of `X` and :math:`e`
a vector of ones.
Parameters
----------
X : array_like
Each column is a series of price or return's asset.
w0 : array_like, optional
Initial weights for the optimizer.
up_bound, low_bound : float, optional
Respectively maximum and minimum values of weights, such that low_bound
:math:`\leq w_i \leq` up_bound :math:`\forall i`. Default is 0 and 1.
cov : callable, optional
Callable mapping the ``(T, N)`` training array to an ``(N, N)``
covariance matrix, e.g.
:func:`fynance.portfolio.covariance.ledoit_wolf`; default `None`
keeps the sample covariance.
Returns
-------
array_like
Weights that minimize the variance of the portfolio.
"""
mat_cov = _estimate_cov(X, cov)
N = X.shape[1]
if N == 1:
return np.ones([1, 1])
up_bound = max(up_bound, 1 / N)
# Clamp low_bound so the box stays compatible with the sum-to-one
# constraint: with low_bound > 1/N the feasible set is empty and SLSQP
# would silently return weights summing to more than one (as HRP/IVP
# already guard against).
low_bound = min(low_bound, 1 / N)
# On return-scale inputs the portfolio variance w'Sigma w is ~1e-4 or
# smaller, which can sit below SLSQP's default ftol so the optimizer stops
# at the 1/N start. Rescale the covariance to unit trace so the objective
# is O(1); the variance-minimizing weights are invariant to this scaling.
scale = np.trace(mat_cov) / N
if scale <= 0:
return np.ones([N, 1]) / N
mat_cov = mat_cov / scale
def f_MVP(w):
w = w.reshape([N, 1])
return w.T @ mat_cov @ w
# Set inital weights
if w0 is None:
w0 = np.ones([N]) / N
# Set constraints and minimze
const_sum = LinearConstraint(np.ones([1, N]), [1], [1])
const_ind = Bounds(low_bound * np.ones([N]), up_bound * np.ones([N]))
result = minimize(
f_MVP,
w0,
method='SLSQP',
constraints=[const_sum],
bounds=const_ind,
options={'ftol': 1e-12, 'maxiter': 1000},
)
return result.x.reshape([N, 1])
# =========================================================================== #
# Maximum Diversification Portfolio developed by Choueifaty and Coignard #
# =========================================================================== #
def _diversified_ratio_from_cov(
w: NDArray[np.float64],
sigma: NDArray[np.float64],
) -> float:
r""" Diversification ratio of weights `w` evaluated from a fixed covariance.
Same ratio as :func:`fynance.metrics.diversified_ratio`,
:math:`D(w) = (w' s) / \sqrt{w' \Sigma w}` with :math:`s =
\sqrt{\text{diag}(\Sigma)}`, but computed from a caller-supplied
:math:`\Sigma` instead of recomputing the sample covariance of the
returns panel on every call — the ``cov=`` path of :func:`MDP`.
Parameters
----------
w : array_like
Portfolio weights, shape ``(N,)`` or ``(N, 1)``.
sigma : array_like
Covariance matrix, shape ``(N, N)``.
Returns
-------
float
Diversification ratio of `w` under `sigma`.
"""
w = np.asarray(w, dtype=np.float64).reshape(-1, 1)
s = np.sqrt(np.diag(sigma)).reshape(-1, 1)
return ((w.T @ s) / np.sqrt(w.T @ sigma @ w)).item()
[docs]
def MDP(
X: NDArray[np.float64],
w0: NDArray[np.float64] | None = None,
up_bound: float = 1.,
low_bound: float = 0.,
cov: Callable[[NDArray[np.float64]], NDArray[np.float64]] | None = None,
) -> NDArray[np.float64]:
r""" Get weights of Maximum Diversified Portfolio allocation.
Notes
-----
Weights of Maximum Diversification Portfolio, as described by Y. Choueifaty
and Y. Coignard [5]_, verify the following problem:
.. math::
w = \text{arg max } D(w) \\
u.c. \begin{cases}w'e = 1 \\ 0 \leq w_i \leq 1 \\ \end{cases}
Where :math:`D(w)` is the diversified ratio of portfolio weighted by `w`.
Parameters
----------
X : array_like
Each column is a series of price or return's asset.
w0 : array_like, optional
Initial weights for the optimizer.
up_bound, low_bound : float, optional
Respectively maximum and minimum values of weights, such that low_bound
:math:`\leq w_i \leq` up_bound :math:`\forall i`. Default is 0 and 1.
cov : callable, optional
Callable mapping the ``(T, N)`` training array to an ``(N, N)``
covariance matrix, e.g.
:func:`fynance.portfolio.covariance.ledoit_wolf`; default `None`
keeps the sample covariance. When `None` (default) the objective
calls :func:`diversified_ratio`, which recomputes the (biased)
sample covariance of `X` at every optimizer iteration. When a
callable is given, the covariance is instead estimated **once**
via that callable and the diversification ratio is evaluated from
this fixed matrix at every iteration (see
:func:`_diversified_ratio_from_cov`) — the two evaluation paths
can therefore differ slightly even when `cov` reduces to the
sample covariance, because :func:`diversified_ratio` uses the
biased (``ddof=0``) sample covariance recomputed from `X` while
`cov` is called once upfront.
Returns
-------
array_like
Weights that maximize the diversified ratio of the portfolio.
See Also
--------
diversified_ratio
References
----------
.. [5] `Choueifaty, Y., and Coignard, Y., 2008, Toward Maximum
Diversification. <https://www.tobam.fr/wp-content/uploads/2014/12/
TOBAM-JoPM-Maximum-Div-2008.pdf>`_
"""
T, N = X.shape
if N == 1:
return np.ones([1, 1])
up_bound = max(up_bound, 1 / N)
if cov is None:
# Set function to minimize: diversified_ratio recomputes the
# (biased) sample covariance of X at every call, unchanged.
def f_max_divers_weights(w):
return -diversified_ratio(X, W=w)
else:
# Estimate the covariance once and maximize the diversification
# ratio evaluated from this fixed matrix at every iteration.
sigma = _estimate_cov(X, cov)
def f_max_divers_weights(w):
return -_diversified_ratio_from_cov(w, sigma)
# Set inital weights
if w0 is None:
w0 = np.ones([N]) / N
# Set constraints and minimze
const_sum = LinearConstraint(np.ones([1, N]), [1], [1])
const_ind = Bounds(low_bound * np.ones([N]), up_bound * np.ones([N]))
result = minimize(
f_max_divers_weights,
w0,
method='SLSQP',
constraints=[const_sum],
bounds=const_ind
)
return result.x.reshape([N, 1])
# =========================================================================== #
# Rolling allocation #
# =========================================================================== #
[docs]
def rolling_allocation(
f: Callable[..., NDArray[np.float64]],
X: NDArray[np.float64],
n: int = 252,
s: int = 63,
ret: bool = True,
drift: bool = True,
**kwargs,
) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
r""" Roll an algorithm of portfolio allocation.
Generic walk-forward backtester for any allocation function ``f``
(e.g. :func:`ERC`, :func:`HRP`, :func:`MVP`). At each step the
weights are estimated on a training window of length ``n`` and held
for the next ``s`` periods. By construction this respects strict
temporal ordering — no future data leaks into the weights — which
is the same no-lookahead pattern used in
:class:`fynance.models.rolling._RollingBasis` for ML models.
Assets that are constant on more than 50% of the training window
are dropped from that step's allocation; the remaining weight mass
is redistributed across the active assets.
Notes
-----
Weights are computed on the past data from ``t - n`` to ``t`` and are
applied to backtest on data from ``t`` to ``t + s``.
.. math::
\forall t \in [n, T], w_{t:t+s} = f(X_{t-n:t})
A ``cov=`` covariance-estimator callable (see
:mod:`fynance.portfolio.covariance`) reaches the allocator ``f``
unchanged through ``**kwargs``, e.g.
``rolling_allocation(MVP_uc, X, n=252, s=63, cov=ledoit_wolf)``.
Parameters
----------
f : callable
Allocation algorithm that take as parameters a subarray of ``X``
and ``**kwargs``, and return a vector (as ``np.ndarray``) of weights.
X : array_like
Data matrix, each columns is a series of prices, indexes or
performances, each row is a observation at time ``t``.
n, s : int
Respectively the number of observations to compute weights and the
number of observations to roll. Default is ``n=252`` and ``s=63``.
ret : bool, optional
If True (default) pass to ``f`` the returns of ``X``. Otherwise pass
``X`` to ``f``.
drift : bool, optional
If False performance of the portfolio is computed as if we rebalance
the weights of asset at each timeframe. Otherwise we let to drift the
weights. Default is True.
**kwargs
Any keyword arguments to pass to ``f``.
Returns
-------
numpy.ndarray
Performance of the portfolio allocated following ``f``, shape
``(T,)``. The first ``n`` observations are held at the initial
value ``100.``.
numpy.ndarray
Weights of the portfolio per period, shape ``(T, n_assets)``.
"""
X = np.asarray(X, dtype=np.float64)
if X.ndim == 1:
X = X.reshape(-1, 1)
X = _ffill(X)
T, K = X.shape
w_mat = np.full((T, K), np.nan)
portfolio = np.full(T, 100.)
if ret:
X_ = np.empty_like(X)
X_[0] = np.nan
X_[1:] = X[1:] / X[:-1] - 1.
else:
X_ = X
# Walk-forward: weights computed on X[t - n:t], held over X[t:t + s].
t = n
while t < T - 1:
# Training window: rows [t - n, t - 1].
sub_X = X_[max(t - n, 0):t]
# Drop assets that are constant on more than 50% of the window.
assets = np.flatnonzero(_active_assets(sub_X, n))
sub_X = _bfill(sub_X)
# Compute weights
if assets.size == 1:
w = np.array([[1.]])
else:
w = np.asarray(f(sub_X[:, assets], **kwargs), dtype=np.float64)
row = np.zeros(K)
row[assets] = w.flatten()
w_mat[t] = row
# Compute portfolio performance over the test window [t, t + s].
t_end = min(t + s, T - 1)
perf = _perf_alloc(_bfill(X[t:t_end + 1][:, assets]), w=w, drift=drift)
portfolio[t:t_end + 1] = portfolio[t] * perf.flatten()
t += s
w_mat = _ffill(w_mat)
w_mat[np.isnan(w_mat)] = 0.
return portfolio, w_mat
# =========================================================================== #
# Tools #
# =========================================================================== #
def _ffill(a: NDArray[np.float64]) -> NDArray[np.float64]:
# Forward-fill NaNs along axis 0, column by column (no pandas needed).
a = a.copy()
for k in range(a.shape[1]):
col = a[:, k]
mask = np.isnan(col)
if not mask.any():
continue
idx = np.where(~mask, np.arange(col.size), 0)
np.maximum.accumulate(idx, out=idx)
col[mask] = col[idx[mask]]
return a
def _bfill(a: NDArray[np.float64]) -> NDArray[np.float64]:
# Backward-fill NaNs along axis 0 (mirror of :func:`_ffill`).
return _ffill(a[::-1])[::-1]
def _active_assets(sub_X: NDArray[np.float64], n: int) -> NDArray[np.bool_]:
# True for columns with less than 50% of identical observations (NaN
# counted as a value), matching the original value-counts filter.
thr = 0.5 * n
K = sub_X.shape[1]
keep = np.empty(K, dtype=bool)
for k in range(K):
col = sub_X[:, k]
nan_count = int(np.isnan(col).sum())
vals = col[~np.isnan(col)]
max_freq = int(np.unique(vals, return_counts=True)[1].max()) if vals.size else 0
keep[k] = max(max_freq, nan_count) < thr
return keep
def _perf_alloc(X: NDArray[np.float64], w: NDArray[np.float64], drift: bool = True) -> NDArray[np.float64]:
# Compute portfolio performance following specified weights
if w.ndim == 1:
w = w.reshape([w.size, 1])
if drift:
return (X / X[0, :]) @ w
perf = np.zeros(X.shape)
perf[1:] = (X[1:] / X[:-1] - 1)
return np.cumprod(perf @ w + 1) # type: ignore[return-value]
def _normalize(w: NDArray[np.float64], low_bound: float = 0., up_bound: float = 1., sum_w: float = 1., max_iter: int = 1000) -> NDArray[np.float64]:
# Iterative algorithm to set bounds. Copy the input so the caller's array
# is never mutated in place.
w = np.array(w, dtype=np.float64)
if up_bound < sum_w / w.size or low_bound > sum_w / w.size:
raise ValueError('Low or up bound exceeded sum weight constraint.')
j = 0
while (min(w) < low_bound or max(w) > up_bound) and j < max_iter:
for i in range(w.size):
w[i] = min(w[i], up_bound)
w[i] = max(w[i], low_bound)
w = sum_w * (w / sum(w))
j += 1
if j >= max_iter:
warnings.warn(
'Iterative normalize algorithm exceeded max iterations',
stacklevel=2,
)
return w
