"""Stats module."""
# Copyright (c) 2023-2026
# Author: Hugo Delatte <hugo.delatte@skfoliolabs.com>
# SPDX-License-Identifier: BSD-3-Clause
# Implementation derived from:
# Precise, Copyright (c) 2021, Peter Cotton.
# Riskfolio-Lib, Copyright (c) 2020-2023, Dany Cajas, Licensed under BSD 3 clause.
# Statsmodels, Copyright (C) 2006, Jonathan E. Taylor, Licensed under BSD 3 clause.
from __future__ import annotations
import math
import random
import warnings
from enum import auto
import cvxpy as cp
import numpy as np
import scipy.cluster.hierarchy as sch
import scipy.linalg as sla
import scipy.optimize as sco
import scipy.sparse.linalg as scl
import scipy.spatial.distance as scd
import scipy.special as scs
from scipy.sparse import csr_matrix
from skfolio.typing import ArrayLike, BoolArray, FloatArray, IntArray
from skfolio.utils.tools import (
AutoEnum,
_validate_non_negative_integer,
_validate_positive_integer,
)
__all__ = [
"CSWeighting",
"CorrelationMethod",
"NBinsMethod",
"assert_is_distance",
"assert_is_square",
"assert_is_symmetric",
"combination_by_index",
"commutation_matrix",
"compute_optimal_n_clusters",
"corr_to_cov",
"cov_nearest",
"cov_to_corr",
"cs_pearson_correlation",
"cs_rank",
"cs_spearman_correlation",
"inverse_multiply",
"inverse_volatility_weights",
"is_cholesky_dec",
"minimize_relative_weight_deviation",
"multiply_by_inverse",
"n_bins_freedman",
"n_bins_knuth",
"rand_weights",
"rand_weights_dirichlet",
"safe_cholesky",
"safe_divide",
"sample_unique_subsets",
"squared_mahalanobis_dist",
"squared_standardized_euclidean_dist",
"symmetric_step_up_matrix",
"symmetrize",
]
_NUMERICAL_THRESHOLD = 1e-12
_MAX_RIDGE_TRIES = 3
_RIDGE_ESCALATION_FACTOR = 10.0
_DISTANCE_ATOL = 1e-5
_DISTANCE_RTOL = 0.0
[docs]
class CorrelationMethod(AutoEnum):
"""Correlation method."""
PEARSON = auto()
SPEARMAN = auto()
[docs]
class CSWeighting(AutoEnum):
"""Cross-sectional weighting."""
BENCHMARK = auto()
REGRESSION = auto()
INVERSE_IDIO_VARIANCE = auto()
IDENTITY = auto()
[docs]
class NBinsMethod(AutoEnum):
"""Enumeration of the Number of Bins Methods.
Parameters
----------
FREEDMAN : str
Freedman method
KNUTH : str
Knuth method
"""
FREEDMAN = auto()
KNUTH = auto()
def safe_divide(
numerator: float | FloatArray,
denominator: float | FloatArray,
fill_value: float = 0.0,
*,
atol: float = 0.0,
) -> float | FloatArray:
"""Safely divide arrays or scalars.
Division is performed elementwise where `abs(denominator) > atol`. Elsewhere and for
any non-finite output produced by the division,`fill_value` is returned.
Parameters
----------
numerator : float or ndarray
Numerator values.
denominator : float or ndarray
Denominator values.
fill_value : float, default=0.0
Replacement value where safe division cannot produce a valid finite result.
atol : float, default=0.0
Absolute tolerance below which the denominator is treated as zero.
Returns
-------
out : float or ndarray
Safe division result.
"""
if atol < 0:
raise ValueError("`atol` must be non-negative.")
fill_value = float(fill_value)
numerator = np.asarray(numerator)
denominator = np.asarray(denominator)
out = np.full(np.broadcast(numerator, denominator).shape, fill_value, dtype=float)
valid = np.abs(denominator) > atol
np.divide(numerator, denominator, out=out, where=valid)
np.nan_to_num(out, copy=False, nan=fill_value, posinf=fill_value, neginf=fill_value)
return float(out) if out.ndim == 0 else out
[docs]
def n_bins_freedman(x: FloatArray) -> int:
"""Compute the optimal histogram bin size using the Freedman-Diaconis rule [1]_.
Parameters
----------
x : ndarray of shape (n_observations,)
The input array.
Returns
-------
n_bins : int
The optimal bin size.
References
----------
.. [1] "On the histogram as a density estimator: L2 theory".
Freedman & Diaconis (1981).
"""
if x.ndim != 1:
raise ValueError("`x` must be a 1d-array")
n = len(x)
p_25, p_75 = np.percentile(x, [25, 75])
d = 2 * (p_75 - p_25) / (n ** (1 / 3))
if d == 0:
return 5
n_bins = max(1, np.ceil((np.max(x) - np.min(x)) / d))
return round(n_bins)
[docs]
def n_bins_knuth(x: FloatArray) -> int:
"""Compute the optimal histogram bin size using Knuth's rule [1]_.
Parameters
----------
x : ndarray of shape (n_observations,)
The input array.
Returns
-------
n_bins : int
The optimal bin size.
References
----------
.. [1] "Optimal Data-Based Binning for Histograms".
Knuth.
"""
x = np.sort(x)
n = len(x)
def func(y: FloatArray) -> float:
y = y[0]
if y <= 0:
return np.inf
bin_edges = np.linspace(x[0], x[-1], int(y) + 1)
hist, _ = np.histogram(x, bin_edges)
return -(
n * np.log(y)
+ scs.gammaln(0.5 * y)
- y * scs.gammaln(0.5)
- scs.gammaln(n + 0.5 * y)
+ np.sum(scs.gammaln(hist + 0.5))
)
n_bins_init = n_bins_freedman(x)
n_bins = sco.fmin(func, n_bins_init, disp=0)[0]
return round(n_bins)
[docs]
def rand_weights_dirichlet(n: int) -> np.array:
"""Produces n random weights that sum to one from a Dirichlet distribution
(uniform distribution over a simplex).
Parameters
----------
n : int
Number of weights.
Returns
-------
weights : ndarray of shape (n, )
The vector of weights.
"""
return np.random.dirichlet(np.ones(n))
[docs]
def rand_weights(n: int, zeros: int = 0, seed: int | None = None) -> FloatArray:
"""Produces n random weights that sum to one from a uniform distribution
(non-uniform distribution over a simplex).
Parameters
----------
n : int
Number of weights.
zeros : int, default=0
The number of weights to randomly set to zeros.
seed : int, optional
Seed for reproducibility. If None, use an unseeded generator.
Returns
-------
weights : ndarray of shape (n, )
The vector of weights.
"""
rng = np.random.default_rng(seed)
k = rng.random(n)
if zeros > 0:
zeros_idx = rng.choice(n, zeros, replace=False)
k[zeros_idx] = 0
return k / k.sum()
[docs]
def is_cholesky_dec(x: FloatArray) -> bool:
"""Returns True if Cholesky decomposition can be computed.
The matrix must be Hermitian (symmetric if real-valued) and positive-definite.
No checking is performed to verify whether the matrix is Hermitian or not.
Parameters
----------
x : ndarray of shape (n, m)
The matrix.
Returns
-------
value : bool
True if Cholesky decomposition can be applied to the matrix, False otherwise.
"""
# Around 100 times faster than checking for positive eigenvalues with np.linalg.eigh
try:
np.linalg.cholesky(x)
return True
except np.linalg.LinAlgError:
return False
def is_positive_definite(x: FloatArray) -> bool:
"""Returns True if the matrix is positive definite.
Parameters
----------
x : ndarray of shape (n, m)
The matrix.
Returns
-------
value : bool
True if the matrix is positive definite, False otherwise.
"""
return np.all(np.linalg.eigvals(x) > 0)
[docs]
def assert_is_square(x: FloatArray) -> None:
"""Raises an error if the matrix is not square.
Parameters
----------
x : ndarray of shape (n, n)
The matrix.
Raises
------
ValueError: if the matrix is not square.
"""
if x.ndim != 2 or x.shape[0] != x.shape[1]:
raise ValueError("The matrix must be square")
[docs]
def assert_is_symmetric(
x: FloatArray, *, rtol: float = 1e-5, atol: float = 1e-8
) -> None:
"""Raises an error if the matrix is not symmetric.
Parameters
----------
x : ndarray of shape (n, m)
The matrix.
rtol : float, default=1e-5
Relative tolerance for `numpy.allclose`.
atol : float, default=1e-8
Absolute tolerance for `numpy.allclose`.
Raises
------
ValueError: if the matrix is not symmetric.
"""
assert_is_square(x)
if not np.allclose(x, x.T, rtol=rtol, atol=atol):
raise ValueError("The matrix must be symmetric")
[docs]
def assert_is_distance(x: FloatArray) -> None:
"""Raises an error if the matrix is not a distance matrix.
Parameters
----------
x : ndarray of shape (n, n)
The matrix.
Raises
------
ValueError: if the matrix is a distance matrix.
"""
assert_is_symmetric(x, rtol=_DISTANCE_RTOL, atol=_DISTANCE_ATOL)
if not np.allclose(np.diag(x), np.zeros(x.shape[0]), atol=_DISTANCE_ATOL):
raise ValueError(
"The distance matrix must have diagonal elements close to zeros"
)
[docs]
def cov_to_corr(cov: FloatArray) -> tuple[FloatArray, FloatArray]:
"""Convert a covariance matrix to a correlation matrix.
Parameters
----------
cov : ndarray of shape (n, n)
Covariance matrix.
Returns
-------
corr, std : tuple[ndarray of shape (n, n), ndarray of shape (n, )]
Correlation matrix and standard-deviation vector
"""
if cov.ndim != 2:
raise ValueError(f"`cov` must be a 2D array, got a {cov.ndim}D array")
std = np.sqrt(np.diag(cov))
corr = safe_divide(cov, np.outer(std, std), fill_value=np.nan)
np.fill_diagonal(corr, 1.0)
return corr, std
[docs]
def corr_to_cov(corr: FloatArray, std: FloatArray):
"""Convert a correlation matrix to a covariance matrix given its
standard-deviation vector.
Parameters
----------
corr : ndarray of shape (n, n)
Correlation matrix.
std : ndarray of shape (n, )
Standard-deviation vector.
Returns
-------
cov : ndarray of shape (n, n)
Covariance matrix
"""
if std.ndim != 1:
raise ValueError(f"`std` must be a 1D array, got a {std.ndim}D array")
if corr.ndim != 2:
raise ValueError(f"`corr` must be a 2D array, got a {corr.ndim}D array")
cov = corr * std * std[:, None]
return cov
_CLIPPING_VALUE = 1e-13
[docs]
def cov_nearest(
cov: FloatArray,
higham: bool = False,
higham_max_iteration: int = 100,
warn: bool = False,
):
"""Compute the nearest covariance matrix that is positive definite and with a
cholesky decomposition that can be computed. The variance is left unchanged.
A covariance matrix that is not positive definite often occurs in high
dimensional problems. It can be due to multicollinearity, floating-point
inaccuracies, or when the number of observations is smaller than the number of
assets.
First, it converts the covariance matrix to a correlation matrix.
Then, it finds the nearest correlation matrix and converts it back to a covariance
matrix using the initial standard deviation.
Cholesky decomposition can fail for symmetric positive definite (SPD) matrix due
to floating point error and inversely, Cholesky decomposition can succeed for
non-SPD matrix. Therefore, we need to test for both. We always start by testing
for Cholesky decomposition which is significantly faster than checking for positive
eigenvalues.
Parameters
----------
cov : ndarray of shape (n, n)
Covariance matrix.
higham : bool, default=False
If this is set to True, the Higham (2002) algorithm [1]_ is used,
otherwise the eigenvalues are clipped to threshold above zeros (1e-13).
The default (`False`) is to use the clipping method as the Higham
algorithm can be slow for large datasets.
higham_max_iteration : int, default=100
Maximum number of iterations of the Higham (2002) algorithm.
The default value is `100`.
warn : bool, default=False
If this is set to True, a user warning is emitted when the covariance matrix
is not positive definite and replaced by the nearest. The default is False.
Returns
-------
cov : ndarray
The nearest covariance matrix.
References
----------
.. [1] "Computing the nearest correlation matrix - a problem from finance"
IMA Journal of Numerical Analysis
Higham (2002)
"""
assert_is_square(cov)
assert_is_symmetric(cov)
# Around 100 times faster than checking eigenvalues with np.linalg.eigh
if is_cholesky_dec(cov) and is_positive_definite(cov):
return cov
if warn:
warnings.warn(
"The covariance matrix is not positive definite. "
f"The {'Higham' if higham else 'Clipping'} algorithm will be used to find "
"the nearest positive definite covariance.",
stacklevel=2,
)
corr, std = cov_to_corr(cov)
if higham:
eps = np.finfo(np.float64).eps * 5
diff = np.zeros(corr.shape)
x = corr.copy()
for _ in range(higham_max_iteration):
x_adj = x - diff
eig_vals, eig_vecs = np.linalg.eigh(x_adj)
x = eig_vecs * np.maximum(eig_vals, eps) @ eig_vecs.T
diff = x - x_adj
np.fill_diagonal(x, 1)
cov = corr_to_cov(x, std)
if is_cholesky_dec(cov) and is_positive_definite(cov):
break
else:
raise ValueError("Unable to find the nearest positive definite matrix")
else:
eig_vals, eig_vecs = np.linalg.eigh(corr)
# Clipping the eigenvalues with a value smaller than 1e-13 can cause scipy to
# consider the matrix non-psd is some corner cases (see test/test_stats.py)
x = eig_vecs * np.maximum(eig_vals, _CLIPPING_VALUE) @ eig_vecs.T
x, _ = cov_to_corr(x)
cov = corr_to_cov(x, std)
return cov
[docs]
def commutation_matrix(x):
"""Compute the commutation matrix.
Parameters
----------
x : ndarray of shape (n, m)
The matrix.
Returns
-------
K : ndarray of shape (m * n, m * n)
The commutation matrix.
"""
(m, n) = x.shape
row = np.arange(m * n)
col = row.reshape((m, n), order="F").ravel()
data = np.ones(m * n, dtype=np.int8)
k = csr_matrix((data, (row, col)), shape=(m * n, m * n))
return k
[docs]
def compute_optimal_n_clusters(distance: FloatArray, linkage_matrix: FloatArray) -> int:
r"""Compute the optimal number of clusters based on Two-Order Difference to Gap
Statistic [1]_.
The Two-Order Difference to Gap Statistic has been developed to improve the
performance and stability of the Tibshiranis Gap statistic.
It applies the two-order difference of the within-cluster dispersion to replace the
reference null distribution in the Gap statistic.
The number of cluster :math:`k` is determined by:
.. math:: \begin{cases}
\begin{aligned}
&\max_{k} & & W_{k+2} + W_{k} - 2 W_{k+1} \\
&\text{s.t.} & & 1 \ge c \ge max\bigl(8, \sqrt{n}\bigr) \\
\end{aligned}
\end{cases}
with :math:`n` the sample size and :math:`W_{k}` the within-cluster dispersions
defined as:
.. math:: W_{k} = \sum_{i=1}^{k} \frac{D_{i}}{2|C_{i}|}
where :math:`|C_{i}|` is the cardinality of cluster :math:`i` and :math:`D_{i}` its
density defined as:
.. math:: D_{i} = \sum_{u \in C_{i}} \sum_{v \in C_{i}} d(u,v)
with :math:`d(u,v)` the distance between u and v.
Parameters
----------
distance : ndarray of shape (n, n)
Distance matrix.
linkage_matrix : ndarray of shape (n - 1, 4)
Linkage matrix.
Returns
-------
value : int
Optimal number of clusters.
References
----------
.. [1] "Application of two-order difference to gap statistic".
Yue, Wang & Wei (2009)
"""
cut_tree = sch.cut_tree(linkage_matrix)
n = cut_tree.shape[1]
max_clusters = min(n, max(8, round(np.sqrt(n))))
dispersion = []
for k in range(max_clusters):
level = cut_tree[:, n - k - 1]
cluster_density = []
for i in range(np.max(level) + 1):
cluster_idx = np.argwhere(level == i).flatten()
cluster_dists = scd.squareform(
distance[cluster_idx, :][:, cluster_idx], checks=False
)
if cluster_dists.shape[0] != 0:
cluster_density.append(np.nan_to_num(cluster_dists.mean()))
dispersion.append(np.sum(cluster_density))
dispersion = np.array(dispersion)
gaps = np.roll(dispersion, -2) + dispersion - 2 * np.roll(dispersion, -1)
gaps = gaps[:-2]
# k=0 represents one cluster
k = np.argmax(gaps) + 2
return k
[docs]
def minimize_relative_weight_deviation(
weights: FloatArray,
min_weights: FloatArray,
max_weights: FloatArray,
solver: str = "CLARABEL",
solver_params: dict | None = None,
) -> FloatArray:
r"""
Apply weight constraints to an initial array of weights by minimizing the relative
weight deviation of the final weights from the initial weights.
.. math::
\begin{cases}
\begin{aligned}
&\min_{w} & & \Vert \frac{w - w_{init}}{w_{init}} \Vert_{2}^{2} \\
&\text{s.t.} & & \sum_{i=1}^{N} w_{i} = 1 \\
& & & w_{min} \leq w_i \leq w_{max}, \quad \forall i
\end{aligned}
\end{cases}
Parameters
----------
weights : ndarray of shape (n_assets,)
Strictly positive initial weights summing to one.
min_weights : ndarray of shape (n_assets,)
Minimum assets weights (weights lower bounds).
max_weights : ndarray of shape (n_assets,)
Maximum assets weights (weights upper bounds).
solver : str, default="CLARABEL"
The solver to use. The default is "CLARABEL" which is written in Rust and has
better numerical stability and performance than ECOS and SCS.
For more details about available solvers, check the CVXPY documentation:
https://www.cvxpy.org/tutorial/advanced/index.html#choosing-a-solver
solver_params : dict, optional
Solver parameters. For example, `solver_params=dict(verbose=True)`.
The default (`None`) is to use the CVXPY default.
For more details about solver arguments, check the CVXPY documentation:
https://www.cvxpy.org/tutorial/advanced/index.html#setting-solver-options
"""
if not (weights.shape == min_weights.shape == max_weights.shape):
raise ValueError("`min_weights` and `max_weights` must have same size")
if np.any(weights <= 0):
raise ValueError("Initial weights must be strictly positive")
if not np.isclose(np.sum(weights), 1.0):
raise ValueError("Initial weights must sum to one")
if np.any(max_weights < min_weights):
raise ValueError("`min_weights` must be lower or equal to `max_weights`")
if np.all((weights >= min_weights) & (weights <= max_weights)):
return weights
if solver_params is None:
solver_params = {}
n = len(weights)
w = cp.Variable(n)
objective = cp.Minimize(cp.norm(w / weights - 1))
constraints = [cp.sum(w) == 1, w >= min_weights, w <= max_weights]
problem = cp.Problem(objective, constraints)
try:
problem.solve(solver=solver, **solver_params)
if w.value is None:
raise cp.SolverError("No solution found")
except (cp.SolverError, scl.ArpackNoConvergence):
raise cp.SolverError(
f"Solver '{solver}' failed. Try another"
" solver, or solve with solver_params=dict(verbose=True) for more"
" information"
) from None
return w.value
[docs]
def combination_by_index(idx: int, n: int, k: int) -> IntArray:
"""
Retrieve the k-combination at a given lexicographic position without enumerating
all combinations.
This function implements the *unranking* algorithm (also known as the combinatorial
number system or "combinadic") to retrieve the specific k-combination corresponding
to a given lexicographic `idx` without generating all C(n, k) possible subsets.
Given a universe of size `n`, there are M = C(n, k) possible subsets of size k.
This function returns the subset corresponding to the `idx` in lex order.
This approach is crucial when M = C(n, k) is too large to generate or store all
combinations, and you need to draw random subsets uniformly (sampling k=5 from n=100
gives M ≈ 7.5e7).
Time complexity: O(k)
Space complexity: O(k)
Parameters
----------
idx : int
Index (rank) of the desired combination in lex order. Must satisfy
0 <= idx < C(n, k).
n : int
Size of the universe.
k : int
Size of each combination (0 <= k <= n).
Returns
-------
combination : ndarray of shape (k,)
1D integer array of length k containing the sorted k-combination.
Raises
------
ValueError
If parameters are out of valid range.
References
----------
..[1] "The Art of Computer Programming", Vol. 4A: Combinatorial Algorithms,
Section 7.2.1.3. Knuth, D. E. (1998).
"""
total = math.comb(n, k)
if idx < 0 or idx >= total:
raise ValueError(
f"Index {idx} out of range for C({n},{k})={total} combinations."
)
combination = np.empty(k, dtype=int)
remaining_rank = idx
next_element = 0
for pos in range(k):
remaining_slots = k - pos
x = next_element
block_size = math.comb(n - x - 1, remaining_slots - 1)
while block_size <= remaining_rank:
remaining_rank -= block_size
x += 1
block_size = math.comb(n - x - 1, remaining_slots - 1)
combination[pos] = x
next_element = x + 1
return combination
[docs]
def sample_unique_subsets(
n: int, k: int, n_subsets: int, random_state: int | None = None
) -> IntArray:
"""
Generate unique k-element subsets from a universe of size n using combinatorial
unranking.
Each subset is drawn without replacement (elements within subset are distinct) and
no subset is repeated across draws. Ranks are sampled uniformly without replacement
over [0, C(n, k)).
Time complexity: O(n_subsets * k)
Space complexity: O(n_subsets * k)
Parameters
----------
n : int
Universe size.
k : int
Subset size (0 <= k <= n).
n_subsets : int
Number of distinct subsets to generate (0 <= n_subsets <= C(n, k)).
random_state : int, RandomState instance, optional
Seed or random state to ensure reproducibility.
Returns
-------
subsets : ndarray of shape (n_subsets, k)
2D integer array of shape (n_subsets, k) where each row is a sorted
k-combination.
Raises
------
ValueError
If any parameters are out of valid ranges.
"""
if n < 0:
raise ValueError(f"n must be non-negative, got {n}.")
if k < 0 or k > n:
raise ValueError(f"k={k} must satisfy 0 <= k <= n={n}.")
total = math.comb(n, k)
if n_subsets < 0 or n_subsets > total:
raise ValueError(
f"n_subsets={n_subsets} must satisfy 0 <= n_subsets <= C({n},{k})={total}."
)
rng = random.Random(random_state)
ranks = rng.sample(range(total), k=n_subsets)
# random.sample has a special-case for range objects that avoids building a list of
# length M=C(n,k) and runs in O(n_subsets) time and space as opposed to
# `choice(total, size=n_subsets, replace=False)` which run in O(M) space and raises
# ArrayMemoryError for very big M.
subsets = np.empty((n_subsets, k), dtype=int)
for i, rank in enumerate(ranks):
subsets[i, :] = combination_by_index(rank, n, k)
return subsets
[docs]
def inverse_multiply(a: FloatArray, b: FloatArray) -> FloatArray:
"""Multiply the inverse of matrix a by matrix b.
We use np.linalg.solve as it tends to produce more accurate results than
np.linalg.inv.
Parameters
----------
a : ndarray of shape (n, n)
Square matrix.
b : ndarray of shape (n, m)
Matrix.
Returns
-------
m : ndarray of shape (n, m)
The inverse of matrix a multiplied by matrix b.
"""
assert_is_square(a)
if a.shape[1] != b.shape[0]:
raise ValueError("Wrong dimension")
return np.linalg.solve(a, b)
[docs]
def multiply_by_inverse(a: FloatArray, b: FloatArray) -> FloatArray:
"""Multiply matrix a by the inverse of matrix b.
We use np.linalg.solve as it tends to produce more accurate results than
np.linalg.inv.
Parameters
----------
a : ndarray of shape (n, m)
Matrix.
b : ndarray of shape (n, n)
Square matrix.
Returns
-------
m : ndarray of shape (n, m)
The matrix a multiplied by the inverse of matrix b.
"""
return inverse_multiply(b.T, a.T).T
[docs]
def symmetric_step_up_matrix(n1: int, n2: int) -> FloatArray:
"""Compute the Symmetric step-up matrix M such that `M @ np.ones(n2) = np.ones(n1)`.
Parameters
----------
n1 : int
First dimension.
n2 : int
Second dimension.
Returns
-------
m : ndarray of shape (n1, n2)
The Symmetric step-up matrix.
"""
assert abs(n1 - n2) <= 1
if n1 == n2:
return np.eye(n1)
if n1 < n2:
return symmetric_step_up_matrix(n2, n1).T * n1 / n2
m = np.zeros((n1, n2))
j_row = np.ones((1, n2)) / n2
e = np.eye(n2)
for j in range(n1):
mj = np.concatenate([e[:j], j_row, e[j:]], axis=0)
m += mj / n1
return m
[docs]
def symmetrize(matrix: FloatArray, where: BoolArray | None = None) -> None:
r"""In-place symmetrization: :math:`M \leftarrow (M + M^T) / 2`.
When `where` is provided, only the sub-block indexed by the mask is
symmetrized, leaving the rest of the matrix untouched. This is useful for
matrices that contain NaN rows/columns where a full transpose would
propagate NaNs into the finite block.
Parameters
----------
matrix : ndarray of shape (n, n)
Square matrix to symmetrize in-place.
where : ndarray of shape (n,), optional
Boolean mask indicating which rows/columns to include. If `None`,
the full matrix is symmetrized.
"""
if where is None or np.all(where):
matrix[:] = 0.5 * (matrix + matrix.T)
elif np.any(where):
idx = np.where(where)[0]
ix = np.ix_(idx, idx)
matrix[ix] = 0.5 * (matrix[ix] + matrix[ix].T)
[docs]
def safe_cholesky(
covariance: FloatArray,
ridge_scale: float = _NUMERICAL_THRESHOLD,
max_tries: int = _MAX_RIDGE_TRIES,
) -> FloatArray:
r"""Compute a Cholesky factor :math:`L` from covariance :math:`\Sigma`.
Fast path: try plain Cholesky on the input as-is.
Fallback: symmetrize and add ridge :math:`\lambda I` with escalation until SPD:
.. math:: \Sigma_{reg} = (\Sigma + \Sigma^T)/2 + \lambda I \approx L L^T
Parameters
----------
covariance : ndarray of shape (n_assets, n_assets)
Covariance matrix :math:`\Sigma`.
ridge_scale : float, default=1e-12
Relative ridge size, as a fraction of the average absolute covariance
diagonal. If that scale is zero, a positive numerical floor is used.
max_tries : int, default=3
Maximum number of ridge escalations before raising an error.
Returns
-------
chol : ndarray of shape (n_assets, n_assets)
Lower triangular Cholesky factor :math:`L` such that
:math:`\Sigma \approx L L^T`.
Raises
------
ValueError
If Cholesky decomposition fails after all retry attempts.
"""
covariance = np.asarray(covariance)
# Fast path: try plain Cholesky on the raw covariance
try:
return sla.cholesky(covariance, lower=True, check_finite=False)
except sla.LinAlgError:
pass
# Symmetrize and apply ridge escalation
cov = 0.5 * (covariance + covariance.T)
base_diag_scale = float(np.mean(np.abs(np.diag(cov))))
if base_diag_scale > 0.0 and np.isfinite(base_diag_scale):
scale = base_diag_scale
else:
scale = max(float(np.max(np.abs(cov))), 1.0)
ridge = max(ridge_scale * scale, np.finfo(float).eps * scale)
for _ in range(max_tries):
cov_reg = cov.copy()
cov_reg[np.diag_indices_from(cov_reg)] += ridge
try:
return sla.cholesky(cov_reg, lower=True, check_finite=False)
except sla.LinAlgError:
ridge *= _RIDGE_ESCALATION_FACTOR
raise ValueError(
f"Cholesky failed after {max_tries} attempts; "
f"last ridge={ridge / _RIDGE_ESCALATION_FACTOR:.3e}."
)
[docs]
def squared_standardized_euclidean_dist(
returns: FloatArray, covariance: FloatArray
) -> float:
r"""Squared standardized Euclidean distance.
.. math:: d^2 = \sum_i (r_i\,/\,\sigma_i)^2
This is the squared Mahalanobis distance using only the diagonal of the
covariance matrix (ignoring correlations).
Parameters
----------
returns : ndarray of shape (n_assets,)
Asset return vector.
covariance : ndarray of shape (n_assets, n_assets)
Covariance matrix.
Returns
-------
float
Sum of squared standardized returns (non-negative).
Under correct calibration: :math:`\mathbb{E}[d^2] = n_{\text{assets}}`.
"""
std = np.sqrt(np.maximum(np.diag(covariance), _NUMERICAL_THRESHOLD))
return float(np.sum((returns / std) ** 2))
def _squared_mahalanobis_dist_from_cholesky(
X: FloatArray, cholesky: FloatArray, mean: FloatArray | None = None
) -> FloatArray | float:
r"""Squared Mahalanobis distance from a pre-computed Cholesky factor.
.. math:: d^2 = (r - \mu)^\top \Sigma^{-1} (r - \mu)
Parameters
----------
X : ndarray of shape (n_observations, n_assets) or (n_assets,)
Price returns of the assets. If 1-D, treated as a single observation and a
scalar is returned.
cholesky : ndarray of shape (n_assets, n_assets)
Lower triangular Cholesky factor :math:`L` such that
:math:`\Sigma = L L^\top`.
mean : ndarray of shape (n_assets,), optional
Mean vector :math:`\mu` subtracted from each row. If `None`, data
are assumed already centred.
Returns
-------
d2 : ndarray of shape (n_observations,) or float
Squared Mahalanobis distances (non-negative).
"""
X = np.asarray(X)
is_1d = X.ndim == 1
X = np.atleast_2d(X)
_, n_assets = X.shape
if cholesky.shape != (n_assets, n_assets):
raise ValueError(
f"cholesky shape {cholesky.shape} is incompatible with "
f"returns shape {X.shape}."
)
if mean is not None:
mean = np.asarray(mean)
if mean.ndim != 1 or mean.shape[0] != n_assets:
raise ValueError(
f"mean must be 1D of length {n_assets}, got shape {mean.shape}."
)
X = X - mean
y = sla.solve_triangular(cholesky, X.T, lower=True, check_finite=False)
d2 = np.maximum(np.sum(y * y, axis=0), 0.0)
if is_1d:
return float(d2[0])
return d2
[docs]
def squared_mahalanobis_dist(
X: FloatArray,
covariance: FloatArray,
mean: FloatArray | None = None,
ridge_scale: float = _NUMERICAL_THRESHOLD,
max_tries: int = _MAX_RIDGE_TRIES,
) -> FloatArray | float:
r"""Squared Mahalanobis distance via Cholesky decomposition.
.. math:: d^2 = (r - \mu)^\top \Sigma^{-1} (r - \mu)
Parameters
----------
X : ndarray of shape (n_observations, n_assets) or (n_assets,)
Price returns of the assets. If 1-D, treated as a single observation and a
scalar is returned.
covariance : ndarray of shape (n_assets, n_assets)
Covariance matrix :math:`\Sigma`.
mean : ndarray of shape (n_assets,), optional
Mean vector :math:`\mu` subtracted from each row. If `None`, data
are assumed already centred.
ridge_scale : float, default=1e-12
Relative ridge size, as a fraction of the average covariance diagonal.
max_tries : int, default=3
Maximum number of ridge escalations before raising an error.
Returns
-------
d2 : ndarray of shape (n_observations,) or float
Squared Mahalanobis distances (non-negative).
"""
chol = safe_cholesky(covariance, ridge_scale=ridge_scale, max_tries=max_tries)
return _squared_mahalanobis_dist_from_cholesky(X, cholesky=chol, mean=mean)
def _cs_pearson_correlation_2d_3d(
a: FloatArray,
b: FloatArray,
weights: FloatArray | None,
min_count: int,
eps: float,
) -> FloatArray:
r"""Weighted Pearson correlation between 2D and 3D cross-sections.
`a` has shape `(n_observations, n_assets)` and `b` has shape
`(n_observations, n_assets, n_series)`. Correlations are computed over the
asset axis and returned with shape `(n_observations, n_series)`.
"""
a_finite = np.isfinite(a)
b_finite = np.isfinite(b)
if weights is None:
a_weight = a_finite.astype(float)
n_valid = np.einsum("tn,tnk->tk", a_weight, b_finite, optimize=True)
weight_sum = n_valid
else:
w = np.asarray(weights, dtype=float)
if np.any(np.isfinite(w) & (w < 0.0)):
raise ValueError("`weights` must be non-negative.")
if w.ndim == 1:
if w.shape[0] != a.shape[1]:
raise ValueError(
f"`weights` length must match the cross-sectional axis, got "
f"{w.shape[0]} and {a.shape[1]}."
)
w = np.broadcast_to(w[np.newaxis, :], a.shape)
elif w.shape != a.shape:
raise ValueError(
f"`weights` must have shape {(a.shape[1],)} or {a.shape}, "
f"got {w.shape}."
)
effective = np.isfinite(w) & (w > 0.0) & a_finite
a_weight = np.where(effective, w, 0.0)
n_valid = np.einsum(
"tn,tnk->tk", effective.astype(float), b_finite, optimize=True
)
weight_sum = np.einsum("tn,tnk->tk", a_weight, b_finite, optimize=True)
a_weight_sum = a_weight.sum(axis=1)
a_center = safe_divide(
np.sum(np.where(a_finite, a_weight * a, 0.0), axis=1),
a_weight_sum,
fill_value=0.0,
)
a_centered = np.where(a_finite, a - a_center[:, np.newaxis], 0.0)
b_centered = np.where(b_finite, b, 0.0)
b_center = safe_divide(
np.einsum("tn,tnk->tk", a_weight, b_centered, optimize=True),
weight_sum,
fill_value=0.0,
)
b_centered -= b_center[:, np.newaxis, :]
b_centered[~b_finite] = 0.0
weighted_a = a_weight * a_centered
a_sum = np.einsum("tn,tnk->tk", weighted_a, b_finite, optimize=True)
a_square_sum = np.einsum(
"tn,tnk->tk", weighted_a * a_centered, b_finite, optimize=True
)
b_sum = np.einsum("tn,tnk->tk", a_weight, b_centered, optimize=True)
b_square_sum = np.einsum(
"tn,tnk,tnk->tk", a_weight, b_centered, b_centered, optimize=True
)
cross_sum = np.einsum("tn,tnk->tk", weighted_a, b_centered, optimize=True)
covariance = cross_sum - safe_divide(a_sum * b_sum, weight_sum, fill_value=np.nan)
a_variance = a_square_sum - safe_divide(a_sum**2, weight_sum, fill_value=np.nan)
b_variance = b_square_sum - safe_divide(b_sum**2, weight_sum, fill_value=np.nan)
a_variance = np.maximum(a_variance, 0.0)
b_variance = np.maximum(b_variance, 0.0)
denom = np.sqrt(a_variance * b_variance)
corr = safe_divide(covariance, denom, fill_value=np.nan, atol=eps)
corr[
(n_valid < min_count)
| (weight_sum <= eps)
| (a_variance <= eps)
| (b_variance <= eps)
] = np.nan
return corr
def _use_cs_pearson_correlation_2d_3d(
a: FloatArray, b: FloatArray, axis: int
) -> tuple[FloatArray, FloatArray] | None:
"""Return ordered 2D/3D inputs for the optimized cross-sectional path."""
if axis not in (1, -2):
return None
if a.ndim == 2 and b.ndim == 3 and a.shape == b.shape[:2]:
return a, b
if a.ndim == 3 and b.ndim == 2 and b.shape == a.shape[:2]:
return b, a
return None
[docs]
def cs_pearson_correlation(
a: FloatArray,
b: FloatArray,
weights: FloatArray | None = None,
axis: int = 0,
min_count: int = 3,
eps: float = 1e-12,
) -> float | FloatArray:
r"""Weighted cross-sectional Pearson correlation.
Computes the weighted Pearson correlation between *a* and *b* along `axis`.
All other dimensions are treated as independent batch dimensions over which the
computation is vectorized.
For vectors :math:`a` and :math:`b` with weights :math:`w`:
.. math::
\rho = \frac{
\sum_n w_n \,(a_n - \bar a)\,(b_n - \bar b)
}{
\sqrt{\sum_n w_n \,(a_n - \bar a)^2}\;
\sqrt{\sum_n w_n \,(b_n - \bar b)^2}
}
where :math:`\bar a = \sum_n w_n a_n / \sum_n w_n` (and likewise for :math:`\bar b`).
Parameters
----------
a : ndarray
First array.
b : ndarray
Second array, broadcastable to the same shape as *a*.
weights : ndarray, optional
Non-negative weights, broadcastable to *a* along `axis`.
`None` uses equal weights. Non-finite and zero weights are excluded
from weighted correlations.
axis : int, default=0
The cross-sectional axis along which correlation is computed.
min_count : int, default=3
Minimum number of effective observations along `axis`. Without
`weights`, this is the number of jointly finite observations. With
`weights`, this is the number of jointly finite observations with
finite strictly positive weight.
eps : float, default=1e-12
Denominator threshold below which `NaN` is returned to guard against
near-constant vectors.
Returns
-------
corr : float or ndarray
Scalar when inputs are 1D, otherwise an array with `axis`
removed.
"""
a = np.asarray(a, dtype=float)
b = np.asarray(b, dtype=float)
fast_path = _use_cs_pearson_correlation_2d_3d(a, b, axis)
if fast_path is not None:
return _cs_pearson_correlation_2d_3d(
fast_path[0], fast_path[1], weights=weights, min_count=min_count, eps=eps
)
try:
shape = np.broadcast_shapes(a.shape, b.shape)
except ValueError:
raise ValueError(
f"`a` and `b` must be broadcastable, got shapes {a.shape} and {b.shape}."
) from None
a = np.broadcast_to(a, shape)
b = np.broadcast_to(b, shape)
if axis < -a.ndim or axis >= a.ndim:
# AxisError moved under np.exceptions in NumPy 2.
raise getattr(np, "exceptions", np).AxisError(axis, a.ndim)
axis = axis % a.ndim
valid = np.isfinite(a) & np.isfinite(b)
if weights is not None:
weight = np.asarray(weights, dtype=float)
if np.any(np.isfinite(weight) & (weight < 0.0)):
raise ValueError("`weights` must be non-negative.")
if weight.ndim == 1:
if weight.shape[0] != shape[axis]:
raise ValueError(
f"`weights` length must match the size of `axis`, got "
f"{weight.shape[0]} and {shape[axis]}."
)
weight_shape = [1] * a.ndim
weight_shape[axis] = shape[axis]
weight = weight.reshape(weight_shape)
else:
while weight.ndim < a.ndim:
weight = weight[..., np.newaxis]
weight = np.broadcast_to(weight, shape).copy()
effective = valid & np.isfinite(weight) & (weight > 0.0)
n_valid = effective.sum(axis=axis)
weight[~effective] = 0.0
weight_sum = weight.sum(axis=axis)
else:
n_valid = valid.sum(axis=axis)
weight_sum = n_valid.astype(float)
a = a.copy()
a[~valid] = 0.0
b = b.copy()
b[~valid] = 0.0
with np.errstate(divide="ignore", invalid="ignore"):
if weights is not None:
a_mean = (weight * a).sum(axis=axis) / weight_sum
b_mean = (weight * b).sum(axis=axis) / weight_sum
else:
a_mean = a.sum(axis=axis) / weight_sum
b_mean = b.sum(axis=axis) / weight_sum
a -= np.expand_dims(a_mean, axis=axis)
b -= np.expand_dims(b_mean, axis=axis)
a[~valid] = 0.0
b[~valid] = 0.0
if weights is not None:
covariance = (weight * a * b).sum(axis=axis)
a_variance = (weight * a**2).sum(axis=axis)
b_variance = (weight * b**2).sum(axis=axis)
else:
covariance = (a * b).sum(axis=axis)
a_variance = (a**2).sum(axis=axis)
b_variance = (b**2).sum(axis=axis)
denom = np.sqrt(a_variance * b_variance)
with np.errstate(divide="ignore", invalid="ignore"):
corr = np.where(denom > eps, covariance / denom, np.nan)
corr = np.where(n_valid >= min_count, corr, np.nan)
return float(corr) if corr.ndim == 0 else corr
[docs]
def cs_rank(a: FloatArray, axis: int = 0) -> FloatArray:
"""Cross-sectional rank along an axis.
Ranks are 1-based. `NaN` values remain `NaN` and are excluded from the ranking
(i.e. only finite values receive ranks).
Parameters
----------
a : ndarray
Input array.
axis : int, default=0
Axis along which to rank.
Returns
-------
ranks : ndarray
Same shape as *a*, dtype `float64`.
"""
a = np.asarray(a, dtype=float)
valid = np.isfinite(a)
a_filled = np.where(valid, a, np.inf)
order = np.argsort(np.argsort(a_filled, axis=axis), axis=axis)
ranks = order.astype(float) + 1.0
ranks[~valid] = np.nan
return ranks
[docs]
def cs_spearman_correlation(
a: FloatArray,
b: FloatArray,
axis: int = 0,
min_count: int = 3,
eps: float = 1e-12,
) -> float | FloatArray:
r"""Cross-sectional Spearman rank correlation.
Ranks the jointly finite values of *a* and *b* along `axis` with :func:`cs_rank`,
then computes their Pearson correlation via :func:`cs_pearson_correlation`
(unweighted).
Parameters
----------
a : ndarray
First array.
b : ndarray
Second array, same shape as *a*.
axis : int, default=0
The cross-sectional axis along which correlation is computed.
min_count : int, default=3
Minimum number of jointly finite observations along `axis`.
eps : float, default=1e-12
Denominator threshold below which `NaN` is returned.
Returns
-------
corr : float or ndarray
Scalar when inputs are 1D, otherwise an array with `axis`
removed.
"""
a = np.asarray(a, dtype=float)
b = np.asarray(b, dtype=float)
try:
shape = np.broadcast_shapes(a.shape, b.shape)
except ValueError:
raise ValueError(
f"`a` and `b` must be broadcastable, got shapes {a.shape} and {b.shape}."
) from None
a = np.broadcast_to(a, shape)
b = np.broadcast_to(b, shape)
valid = np.isfinite(a) & np.isfinite(b)
if np.all(valid):
a_rank = cs_rank(a, axis=axis)
b_rank = cs_rank(b, axis=axis)
else:
a_rank = cs_rank(np.where(valid, a, np.nan), axis=axis)
b_rank = cs_rank(np.where(valid, b, np.nan), axis=axis)
return cs_pearson_correlation(
a_rank,
b_rank,
axis=axis,
min_count=min_count,
eps=eps,
)
[docs]
def inverse_volatility_weights(covariance: FloatArray) -> FloatArray:
r"""Inverse-volatility portfolio weights from a covariance matrix.
Computes weights proportional to the inverse standard deviation:
:math:`w_i \propto 1/\sigma_i`, normalized to sum to 1.
Parameters
----------
covariance : ndarray of shape (n, n)
Covariance matrix.
Returns
-------
w : ndarray of shape (n,)
Normalized weights summing to 1.
"""
std = np.sqrt(np.maximum(np.diag(covariance), _NUMERICAL_THRESHOLD))
w = 1.0 / std
return w / w.sum()
def _forward_mean_return(X: ArrayLike, horizon: int = 5, lag: int = 1) -> FloatArray:
r"""Compute a forward H-period mean return target.
Computes the mean of the next H returns for each observation. This is used
as the regression target in alpha estimation models.
Under the as-of indexing convention, the target for observation :math:`t` is:
.. math::
y_t = \frac{1}{H}\sum_{s=\ell}^{\ell + H - 1} X_{t+s}
where :math:`\ell` is the signal lag. The default `lag=1` means scores at
:math:`t` predict returns starting at :math:`t{+}1`. An internal `lag=0`
is supported for arrays that are already shifted into a forecast/realization
alignment. The last
`lag + horizon - 1` rows are NaN due to incomplete forward windows.
Parameters
----------
horizon : int, default=5
Number of forward periods to average. Must be >= 1.
* `horizon=1`: Returns next-period values.
* `horizon>1`: Returns mean of next H periods.
lag : int, default=1
Number of periods between the score date and the first return in the target
window. Must be >= 0. Public alpha estimators validate `signal_lag >= 1` to
respect the as-of time-indexing convention.
Returns
-------
y : ndarray of shape (T, N)
Forward mean returns. Last `lag + horizon - 1` rows are NaN.
Notes
-----
The computation uses an O(n) cumsum algorithm and handles NaN values by
excluding them from both sum and count (i.e., `nanmean` semantics).
"""
_validate_positive_integer(horizon, "horizon")
_validate_non_negative_integer(lag, "lag")
horizon = int(horizon)
lag = int(lag)
X = np.asarray(X)
if X.ndim != 2:
raise ValueError(f"X must be 2D (T, N), got {X.shape}")
n_observations, n_assets = X.shape
if n_observations == 0:
return np.empty((0, n_assets), dtype=np.float64)
# NaN-safe sums: replace NaN with 0, and count valid observations
X0 = np.nan_to_num(X, nan=0.0)
V0 = np.isfinite(X).astype(np.float64)
target_gap = lag + horizon - 1
# Pad target_gap rows at the end so incomplete forward windows are treated as missing.
# (n_obs + target_gap, n_assets)
X0 = np.vstack([X0, np.zeros((target_gap, n_assets))])
V0 = np.vstack([V0, np.zeros((target_gap, n_assets))])
# Prefix a zero row so csum has length (n_obs + target_gap + 1, n_assets)
csum = np.vstack([np.zeros((1, n_assets)), np.cumsum(X0, axis=0)])
ccnt = np.vstack([np.zeros((1, n_assets)), np.cumsum(V0, axis=0)])
# For each t, window is [t + lag, t + lag + horizon) in the padded array.
start = np.arange(n_observations) + lag
end = np.arange(n_observations) + lag + horizon
# (n_obs, n_assets)
sum_fwd = csum[end] - csum[start]
cnt_fwd = ccnt[end] - ccnt[start]
y = safe_divide(sum_fwd, cnt_fwd, fill_value=np.nan)
if target_gap > 0:
y[-target_gap:] = np.nan
return y
def _market_returns(
asset_returns: ArrayLike,
weights: ArrayLike,
estimation_mask: ArrayLike | None = None,
) -> FloatArray:
"""Compute market returns on the estimation universe.
Parameters
----------
asset_returns : array-like of shape (n_observations, n_assets)
Asset returns at each observation.
weights : array-like of shape (n_observations, n_assets)
Asset weights per observation, typically market capitalizations.
estimation_mask : array-like of shape (n_observations, n_assets), optional
Boolean mask indicating which entries are eligible for market-return
construction. If `None`, all entries are eligible.
Returns
-------
market_ret : ndarray of shape (n_observations,)
Market (cap-weighted) return at each date.
Raises
------
ValueError
If inputs do not have matching 2D shapes.
ValueError
If no eligible asset has finite returns and finite positive total
weight at any observation.
"""
asset_returns = np.asarray(asset_returns, dtype=float)
weights = np.asarray(weights, dtype=float)
if asset_returns.ndim != 2:
raise ValueError(
"asset_returns must be a 2D array of shape (n_observations, n_assets)."
)
if weights.shape != asset_returns.shape:
raise ValueError(
"weights must have the same shape as asset_returns; "
f"got weights.shape={weights.shape} and "
f"asset_returns.shape={asset_returns.shape}."
)
if estimation_mask is None:
estimation_mask = np.ones(asset_returns.shape, dtype=bool)
else:
estimation_mask = np.asarray(estimation_mask, dtype=bool)
if estimation_mask.shape != asset_returns.shape:
raise ValueError(
"estimation_mask must have the same shape as asset_returns; "
f"got estimation_mask.shape={estimation_mask.shape} and "
f"asset_returns.shape={asset_returns.shape}."
)
valid = estimation_mask & np.isfinite(asset_returns) & np.isfinite(weights)
weights = np.where(valid, weights, 0.0)
asset_returns = np.where(valid, asset_returns, 0.0)
w_sum = weights.sum(axis=1, keepdims=True)
valid_rows = w_sum[:, 0] > 0
if not np.all(valid_rows):
bad_obs = int(np.where(~valid_rows)[0][0])
raise ValueError(
"Market return is undefined because no estimable asset has finite "
f"returns and finite positive total weight at observation index {bad_obs}."
)
norm_w = np.divide(
weights, w_sum, out=np.zeros_like(weights, dtype=float), where=w_sum > 0
)
return np.sum(norm_w * asset_returns, axis=1)
