skfolio.moments.GeodesicShrinkageCovariance#

class skfolio.moments.GeodesicShrinkageCovariance(covariance_estimator=None, shrinkage=0.1, target=SCALED_IDENTITY, nearest=True, higham=False, higham_max_iteration=100)[source]#

Covariance estimator with geodesic shrinkage.

Classical shrinkage estimators (e.g. ShrunkCovariance, LedoitWolf) regularize the sample covariance matrix S towards a well-conditioned target T using a linear (Euclidean) combination:

\[(1 - \alpha) \cdot S + \alpha \cdot T\]

GeodesicShrinkageCovariance instead moves S towards T along the geodesic of the manifold of Symmetric Positive Definite (SPD) matrices under the affine-invariant Riemannian metric (AIRM) [1] [2]:

\[\Sigma(\alpha) = S^{1/2} \left(S^{-1/2} \, T \, S^{-1/2}\right)^{\alpha} S^{1/2}\]

At shrinkage $\alpha = 0$, the estimate is the starting covariance from covariance_estimator, after any repair requested by nearest. At $\alpha = 1$, it is target, subject to the final nearest repair.

The estimate remains SPD for every shrinkage in [0, 1] whenever S and target are themselves positive definite.

When target commutes with S, their shared eigenvectors are preserved and their eigenvalues are interpolated geometrically rather than arithmetically. For the default SCALED_IDENTITY target $\mu I$, with $\mu = \mathrm{trace}(S) / n$, the interpolated eigenvalues are $\lambda_i^{1-\alpha} \mu^{\alpha}$, whereas linear shrinkage gives $(1 - \alpha) \lambda_i + \alpha \mu$. The condition number of the geodesic estimate is exactly $\kappa(S)^{1-\alpha}$. Unlike linear shrinkage, the trace is generally not preserved at intermediate values.

The general interpolation is computed with a generalized eigendecomposition relative to the target, avoiding an explicit inverse square root of S.

Parameters:

covariance_estimatorBaseCovariance, optional

Covariance estimator used to compute the starting covariance matrix S (the estimate returned when shrinkage=0). The default (None) is to use EmpiricalCovariance.

shrinkagefloat, default=0.1

Shrinkage intensity alpha in the geodesic interpolation. Must be between 0 (no shrinkage) and 1 (fully shrunk to target) inclusive. The default value is 0.1.

targetGeodesicShrinkageTarget or array-like of shape (n_assets, n_assets), default=GeodesicShrinkageTarget.SCALED_IDENTITY

The shrinkage target T. The string values "scaled_identity" and "diagonal" are also accepted:

SCALED_IDENTITY: $\mu \cdot I$, with $\mu = \mathrm{trace}(S) / n$. This is the same default target used by scikit-learn’s ShrunkCovariance. Its condition number is exactly 1, so shrinkage controls how far the estimate moves from S towards a perfectly well-conditioned matrix.

DIAGONAL: $\mathrm{diag}(S)$. At shrinkage=1 the variances equal those of S and all correlations are zero. At intermediate values the geodesic also changes the variances.

array-like: a user-provided SPD target matrix of shape (n_assets, n_assets). A well-conditioned target is recommended for numerical accuracy.

nearestbool, default=True

If this is set to True, the starting covariance is repaired before interpolation and the resulting covariance is repaired afterwards. Custom targets must be positive definite and are never repaired. The covariance is replaced by 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. For more details, see cov_nearest. The default is True.

highambool, default=False

If this is set to True, the Higham (2002) algorithm is used to find the nearest PD covariance, otherwise the eigenvalues are clipped to a threshold above zeros (1e-13). The default is False and uses the clipping method as the Higham algorithm can be slow for large datasets.

higham_max_iterationint, default=100

Maximum number of iterations of the Higham (2002) algorithm. The default value is 100.

Attributes:

covariance_ndarray of shape (n_assets, n_assets): Estimated covariance.
covariance_estimator_BaseCovariance: Fitted covariance_estimator.
location_ndarray of shape (n_assets,): Estimated mean, available when the fitted covariance estimator exposes it.
n_features_in_int: Number of assets seen during fit.
feature_names_in_ndarray of shape (n_features_in_,): Names of assets seen during fit. Defined only when X has assets names that are all strings.

Methods

`fit`(X[, y])	Fit the Geodesic Shrinkage Covariance estimator.
`get_metadata_routing`()	Get metadata routing of this object.
`get_params`([deep])	Get parameters for this estimator.
`mahalanobis`(X_test)	Compute the squared Mahalanobis distance of observations.
`score`(X_test[, y])	Compute the mean log-likelihood of observations under the estimated model.
`set_params`(**params)	Set the parameters of this estimator.
`set_score_request`(*[, X_test])	Configure whether metadata should be requested to be passed to the `score` method.

References

[1]

“Positive Definite Matrices”. Bhatia, R. (2007). Princeton University Press.

[2]

“Geodesically parameterized covariance estimation”. Musolas, A., Smith, S.T. & Marzouk, Y. (2021). SIAM Journal on Matrix Analysis and Applications.

fit(X, y=None, **fit_params)[source]#

Fit the Geodesic Shrinkage Covariance estimator.

Parameters:

Xarray-like of shape (n_observations, n_assets): Price returns of the assets.
yIgnored: Not used, present for API consistency by convention.
**fit_paramsdict: Parameters to pass to the underlying estimators. Only available if enable_metadata_routing=True, which can be set by using sklearn.set_config(enable_metadata_routing=True). See Metadata Routing User Guide for more details.

Returns:

selfGeodesicShrinkageCovariance: Fitted estimator.

get_metadata_routing()[source]#

Get metadata routing of this object.

Please check User Guide on how the routing mechanism works.

Returns:

routingMetadataRequest: A MetadataRequest encapsulating routing information.

get_params(deep=True)#

Get parameters for this estimator.

Parameters:

deepbool, default=True: If True, will return the parameters for this estimator and contained subobjects that are estimators.

Returns:

paramsdict: Parameter names mapped to their values.

mahalanobis(X_test)#

Compute the squared Mahalanobis distance of observations.

The squared Mahalanobis distance of an observation $r$ is defined as:

\[d^2 = (r - \mu)^T \Sigma^{-1} (r - \mu)\]

where $\Sigma$ is the estimated covariance matrix (self.covariance_) and $\mu$ is the estimated mean (self.location_ if available, otherwise zero).

This distance measure accounts for correlations between assets and is useful for:

Outlier detection in portfolio returns
Risk-adjusted distance calculations
Identifying unusual market regimes

Parameters:

X_testarray-like of shape (n_observations, n_assets) or (n_assets,): Observations for which to compute the squared Mahalanobis distance. Each row represents one observation. If 1D, treated as a single observation. Assets with non-finite fitted variance are excluded from inference. After this asset-level filtering, each row is evaluated using the remaining available values only, covering row-level missing values such as market holidays or pre/post-listing. When rows have different observation patterns, the returned distances follow $\chi^2$ distributions with different degrees of freedom. Rows with no finite retained observation return NaN.

Returns:

distancesndarray of shape (n_observations,) or float: Squared Mahalanobis distance for each observation. Returns a scalar if input is 1D.

Examples

>>> import numpy as np
>>> from skfolio.moments import EmpiricalCovariance
>>> rng = np.random.default_rng(0)
>>> X = rng.standard_normal((100, 3))
>>> model = EmpiricalCovariance()
>>> model.fit(X)
EmpiricalCovariance()
>>> distances = model.mahalanobis(X)
>>> # The mean squared distance should be close to the number of assets (3).
>>> print(distances.mean())
2.9...

score(X_test, y=None)#

Compute the mean log-likelihood of observations under the estimated model.

Evaluates how well the fitted covariance matrix explains new observations, assuming a multivariate Gaussian distribution. This is useful for:

Model selection (comparing different covariance estimators)
Cross-validation of covariance estimation methods
Assessing goodness-of-fit

The log-likelihood for a single observation $r$ is:

\[\log p(r | \mu, \Sigma) = -\frac{1}{2} \left[ n \log(2\pi) + \log|\Sigma| + (r - \mu)^T \Sigma^{-1} (r - \mu) \right]\]

where $n$ is the number of assets, $\Sigma$ is the estimated covariance matrix (self.covariance_), and $\mu$ is the estimated mean (self.location_ if available, otherwise zero).

Parameters:

X_testarray-like of shape (n_observations, n_assets): Observations for which to compute the log-likelihood. Typically held-out test data not used during fitting. Assets with non-finite fitted variance are excluded from inference. This typically happens when the fitted covariance cannot be estimated for an asset, for example before listing, after delisting, or during a warmup period. After this asset-level filtering, each row of X_test is scored using the remaining available values only. This covers row-level missing values in X_test, such as market holidays or pre/post-listing.
yIgnored: Not used, present for scikit-learn API consistency.

Returns:

scorefloat: Mean log-likelihood of the observations. Higher values indicate better fit. The score is averaged over all observations.

Examples

>>> import numpy as np
>>> from skfolio.moments import EmpiricalCovariance, LedoitWolf
>>> rng = np.random.default_rng(0)
>>> X_train = rng.standard_normal((100, 5))
>>> X_test = rng.standard_normal((50, 5))
>>> emp = EmpiricalCovariance().fit(X_train)
>>> lw = LedoitWolf().fit(X_train)
>>> # Compare models on held-out data
>>> print("Empirical:", emp.score(X_test))
Empirical: -6.97...
>>> print("LedoitWolf:", lw.score(X_test))
LedoitWolf: -6.88...

set_params(**params)#

Set the parameters of this estimator.

The method works on simple estimators as well as on nested objects (such as Pipeline). The latter have parameters of the form <component>__<parameter> so that it’s possible to update each component of a nested object.

Parameters:

**paramsdict: Estimator parameters.

Returns:

selfestimator instance: Estimator instance.

set_score_request(*, X_test='$UNCHANGED$')#

Configure whether metadata should be requested to be passed to the score method.

Note that this method is only relevant when this estimator is used as a sub-estimator within a meta-estimator and metadata routing is enabled with enable_metadata_routing=True (see sklearn.set_config). Please check the User Guide on how the routing mechanism works.

The options for each parameter are:

True: metadata is requested, and passed to score if provided. The request is ignored if metadata is not provided.
False: metadata is not requested and the meta-estimator will not pass it to score.
None: metadata is not requested, and the meta-estimator will raise an error if the user provides it.
str: metadata should be passed to the meta-estimator with this given alias instead of the original name.

The default (sklearn.utils.metadata_routing.UNCHANGED) retains the existing request. This allows you to change the request for some parameters and not others.

Added in version 1.3.

Parameters:

X_teststr, True, False, or None, default=sklearn.utils.metadata_routing.UNCHANGED: Metadata routing for X_test parameter in score.

Returns:

selfobject: The updated object.