"""Bivariate Student's t Copula Estimation."""
# Copyright (c) 2023-2026
# Author: Hugo Delatte <hugo.delatte@skfoliolabs.com>
# Credits: Matteo Manzi, Vincent Maladière, Carlo Nicolini
# SPDX-License-Identifier: BSD-3-Clause
from __future__ import annotations
import math
import numpy as np
import scipy.integrate as si
import scipy.optimize as so
import scipy.special as sp
import scipy.stats as st
import sklearn.utils.validation as skv
from skfolio.distribution.copula._base import _RHO_BOUNDS, BaseBivariateCopula
from skfolio.distribution.copula._utils import _apply_margin_swap
from skfolio.typing import ArrayLike, FloatArray
# Student's t copula with dof less than 2.0 is so extremely heavy-tailed that even the
# mean (and many moments) of the distribution do not exist. So Impractical in practice,
# and dof above 50 tends to a Gaussian Copula so we limit it to the interval [2, 50] for
# improved stability and robustness.
_DOF_BOUNDS = (2.0, 50.0)
[docs]
class StudentTCopula(BaseBivariateCopula):
r"""Bivariate Student's t Copula Estimation.
The bivariate Student's t copula density is defined as:
.. math::
C_{\nu, \rho}(u, v) = T_{\nu, \rho} \Bigl(t_{\nu}^{-1}(u),\;t_{\nu}^{-1}(v)\Bigr)
where:
- :math:`\nu > 0` is the degrees of freedom.
- :math:`\rho \in (-1, 1)` is the correlation coefficient.
- :math:`T_{\nu, \rho}(x, y)` is the CDF of the bivariate t-distribution.
- :math:`t_{\nu}^{-1}(p)` is the quantile function (inverse CDF) of the
univariate t-distribution.
Student's t copula with degrees of freedom (dof) less than 2.0 is extremely
heavy-tailed, to the extent that even the mean (and many moments) do not exist,
rendering it impractical. Conversely, for dof above 50 the t copula behaves
similarly to a Gaussian copula. Thus, for improved stability and robustness,
the dof is limited to the interval [2, 50].
.. note::
Rotations are not needed for elliptical copula (e.g., Gaussian or Student-t)
because its correlation parameter :math:`\rho \in (-1, 1)` naturally covers
both positive and negative dependence, and they exhibit symmetric tail behavior.
Parameters
----------
itau : bool, default=True
itau : bool, default=True
If True, :math:`\rho` is estimated using the Kendall's tau inversion method;
otherwise, we use the MLE (Maximum Likelihood Estimation) method. The MLE is
slower but more accurate.
kendall_tau : float, optional
If `itau` is True and `kendall_tau` is provided, this
value is used; otherwise, it is computed.
tolerance : float, default=1e-4
Convergence tolerance for the MLE optimization.
random_state : int, RandomState instance or None, default=None
Seed or random state to ensure reproducibility.
Attributes
----------
rho_ : float
Fitted correlation coefficient (:math:`\rho`) in [-1, 1].
dof_ : float
Fitted degrees of freedom (:math:`\nu`) > 2.
Examples
--------
>>> from skfolio.datasets import load_sp500_dataset
>>> from skfolio.preprocessing import prices_to_returns
>>> from skfolio.distribution import StudentTCopula, compute_pseudo_observations
>>>
>>> # Load historical prices and convert them to returns
>>> prices = load_sp500_dataset()
>>> X = prices_to_returns(prices).tail(504)
>>> X = X[["AAPL", "JPM"]]
>>>
>>> # Convert returns to pseudo observation in the interval [0,1]
>>> X = compute_pseudo_observations(X)
>>>
>>> # Initialize the Copula estimator
>>> model = StudentTCopula()
>>>
>>> # Fit the model to the data.
>>> model.fit(X)
StudentTCopula(...)
>>>
>>> # Display the fitted parameter and tail dependence coefficients
>>> print(model.fitted_repr)
StudentTCopula(rho=0.385, dof=4.90)
>>> print(model.lower_tail_dependence)
0.157...
>>> print(model.upper_tail_dependence)
0.157...
>>>
>>> # Compute the log-likelihood, total log-likelihood, CDF, Partial Derivative,
>>> # Inverse Partial Derivative, AIC, and BIC
>>> log_likelihood = model.score_samples(X)
>>> score = model.score(X)
>>> cdf = model.cdf(X)
>>> p = model.partial_derivative(X)
>>> u = model.inverse_partial_derivative(X)
>>> aic = model.aic(X)
>>> bic = model.bic(X)
>>>
>>> # Generate 5 new samples
>>> samples = model.sample(n_samples=5)
>>>
>>> # Plot the tail concentration function.
>>> fig = model.plot_tail_concentration()
>>>
>>> # Plot a 2D contour of the estimated PDF.
>>> fig = model.plot_pdf_2d()
>>>
>>> # Plot a 3D surface of the estimated PDF.
>>> fig = model.plot_pdf_3d()
References
----------
.. [1] "An Introduction to Copulas (2nd ed.)",
Nelsen (2006)
.. [2] "Multivariate Models and Dependence Concepts",
Joe, Chapman & Hall (1997)
.. [3] "Quantitative Risk Management: Concepts, Techniques and Tools",
McNeil, Frey & Embrechts (2005)
.. [4] "The t Copula and Related Copulas",
Demarta & McNeil (2005)
.. [5] "Copula Methods in Finance",
Cherubini, Luciano & Vecchiato (2004)
"""
rho_: float
dof_: float
_n_params = 2
def __init__(
self,
itau: bool = True,
kendall_tau: float | None = None,
tolerance: float = 1e-4,
random_state: int | None = None,
) -> None:
super().__init__(random_state=random_state)
self.itau = itau
self.kendall_tau = kendall_tau
self.tolerance = tolerance
[docs]
def fit(self, X: ArrayLike, y: None = None) -> StudentTCopula:
r"""Fit the Bivariate Student's t Copula.
If `itau` is True, it uses a Kendall-based two-step method:
- Estimates the correlation parameter (:math:`\rho`) from Kendall's
tau inversion.
- Optimizes the degrees of freedom (:math:`\nu`) by maximizing the
log-likelihood.
Otherwise, it uses the full MLE method: optimizes both :math:`\rho` and
:math:`\nu` by maximizing the log-likelihood.
Parameters
----------
X : array-like of shape (n_observations, 2)
An array of bivariate inputs `(u, v)` where each row represents a
bivariate observation. Both `u` and `v` must be in the interval [0, 1],
having been transformed to uniform marginals.
y : None
Ignored. Provided for compatibility with scikit-learn's API.
Returns
-------
self : StudentTCopula
Returns the instance itself.
"""
X = self._validate_X(X, reset=True)
if self.kendall_tau is None:
kendall_tau = st.kendalltau(X[:, 0], X[:, 1]).statistic
else:
kendall_tau = self.kendall_tau
# Either used directly or for initial guess
rho_from_tau = np.clip(
np.sin((np.pi * kendall_tau) / 2.0),
a_min=_RHO_BOUNDS[0],
a_max=_RHO_BOUNDS[1],
)
if self.itau:
res = so.minimize_scalar(
_neg_log_likelihood,
args=(
rho_from_tau,
X,
),
bounds=_DOF_BOUNDS,
method="bounded",
options={"xatol": self.tolerance},
)
if not res.success:
raise RuntimeError(f"Optimization failed: {res.message}")
self.dof_ = res.x
self.rho_ = rho_from_tau
else:
# We'll use L-BFGS-B for the optimization because:
# 1) The bivariate Student-t copula's negative log-likelihood is smooth,
# making gradient-based methods more efficient than derivative-free
# methods.
# 2) L-BFGS-B directly supports simple box bounds (e.g., -1 < rho < 1,
# 0 < nu < 50).
# 3) It's typically faster and more stable for small-dimensional problems
# than more general constraint solvers (like trust-constr or SLSQP)
result = so.minimize(
fun=lambda x: _neg_log_likelihood(dof=x[0], rho=x[1], X=X),
x0=np.array([3.0, rho_from_tau]),
bounds=(_DOF_BOUNDS, _RHO_BOUNDS),
method="L-BFGS-B",
tol=self.tolerance,
)
if not result.success:
raise RuntimeError(f"Optimization failed: {result.message}")
self.dof_, self.rho_ = result.x
return self
[docs]
def cdf(self, X: ArrayLike) -> FloatArray:
r"""Compute the CDF of the bivariate Student-t copula.
Parameters
----------
X : array-like of shape (n_observations, 2)
An array of bivariate inputs `(u, v)` where each row represents a
bivariate observation. Both `u` and `v` must be in the interval `[0, 1]`,
having been transformed to uniform marginals.
Returns
-------
cdf : ndarray of shape (n_observations,)
CDF values for each observation in X.
Raises
------
RuntimeError
If the numerical integration does not converge.
Notes
-----
The CDF is computed by deterministic adaptive quadrature over a finite
angular interval. For interior inputs, let :math:`x=t_\nu^{-1}(u)`,
:math:`y=t_\nu^{-1}(v)`, :math:`a=(x+y)/(2\sqrt{\nu})`, and
:math:`b=(x-y)/(2\sqrt{\nu})`. Then
.. math::
C(u,v) = \max(u+v-1, 0) + \frac{1}{2\pi}
\int_0^{\arccos(-\rho)}
\left(1 + \frac{a^2}{\sin^2(\theta/2)}
+ \frac{b^2}{\cos^2(\theta/2)}\right)^{-\nu/2} d\theta.
This follows by integrating the CDF's derivative with respect to
correlation from -1 to :math:`\rho`, with the substitution
:math:`r=-\cos(\theta)`.
"""
skv.check_is_fitted(self)
X = self._validate_X(X, reset=False, clip=False)
return _student_t_copula_cdf(X, rho=self.rho_, dof=self.dof_)
[docs]
def survival_probability(self, X: ArrayLike) -> FloatArray:
r"""Return the probability that both variables exceed their respective
thresholds.
Parameters
----------
X : array-like of shape (n_observations, 2)
Uniform marginal probability levels `(u, v)` in `[0, 1]`. Each row
specifies the thresholds for the two variables.
Returns
-------
probabilities : ndarray of shape (n_observations,)
Joint probabilities :math:`P(U > u, V > v)` for each row in X.
Raises
------
RuntimeError
If the numerical integration does not converge.
See Also
--------
cdf : Probability that both variables are below their respective thresholds.
tail_concentration : Conditional probabilities at a common quantile.
Notes
-----
Radial symmetry gives :math:`S(u, v) = C(1-u, 1-v)`. Evaluating the
reflected CDF avoids subtracting nearly equal probabilities in the upper
tail. It uses the same deterministic quadrature as :meth:`cdf`.
"""
skv.check_is_fitted(self)
X = self._validate_X(X, reset=False, clip=False)
return _student_t_copula_cdf(1.0 - X, rho=self.rho_, dof=self.dof_)
[docs]
def partial_derivative(
self, X: ArrayLike, first_margin: bool = False
) -> FloatArray:
r"""Compute the h-function (partial derivative) for the bivariate Student's t
copula.
The h-function with respect to the second margin represents the conditional
distribution function of :math:`u` given :math:`v`:
.. math:: \begin{aligned}
h(u \mid v) &= \frac{\partial C(u,v)}{\partial v} \\
&= t_{\nu+1}\!\left(\frac{t_\nu^{-1}(u) - \rho\,t_\nu^{-1}(v)}
{\sqrt{\frac{(1-\rho^2)\left(\nu + \left(t_\nu^{-1}(v)\right)^2\right)}{\nu+1}}}\right).
\end{aligned}
where:
- :math:`\nu > 0` is the degrees of freedom.
- :math:`\rho \in (-1, 1)` is the correlation coefficient.
- :math:`t_{\nu}^{-1}(p)` is the quantile function (inverse CDF) of the
univariate \(t\)-distribution.
Parameters
----------
X : array-like of shape (n_observations, 2)
An array of bivariate inputs `(u, v)` where each row represents a
bivariate observation. Both `u` and `v` must be in the interval `[0, 1]`,
having been transformed to uniform marginals.
first_margin : bool, default=False
If True, compute the partial derivative with respect to the first
margin `u`; otherwise, compute the partial derivative with respect to the
second margin `v`.
Returns
-------
p : ndarray of shape (n_observations,)
h-function values :math:`h(u \mid v) \;=\; p` for each observation in X.
"""
skv.check_is_fitted(self)
X = self._validate_X(X, reset=False)
X = _apply_margin_swap(X, first_margin=first_margin)
# Compute the inverse CDF (percent point function) using stdtrit for better
# performance
u_inv, v_inv = sp.stdtrit(self.dof_, X).T
# Compute the denominator: sqrt((1 - rho^2) * (nu + y^2) / (nu + 1))
z = (u_inv - self.rho_ * v_inv) / (
np.sqrt((1 - self.rho_**2) * (self.dof_ + v_inv**2) / (self.dof_ + 1))
)
# Student's t CDF with (nu+1) degrees of freedom using stdtr for better
# performance
p = sp.stdtr(self.dof_ + 1, z)
return p
[docs]
def inverse_partial_derivative(
self, X: ArrayLike, first_margin: bool = False
) -> FloatArray:
r"""Compute the inverse of the bivariate copula's partial derivative, commonly
known as the inverse h-function [1]_.
Let :math:`C(u, v)` be a bivariate copula. The h-function with respect to the
second margin is defined by
.. math::
h(u \mid v) \;=\; \frac{\partial\,C(u, v)}{\partial\,v},
which is the conditional distribution of :math:`U` given :math:`V = v`.
The **inverse h-function**, denoted :math:`h^{-1}(p \mid v)`, is the unique
value :math:`u \in [0,1]` such that
.. math::
h(u \mid v) \;=\; p,
\quad \text{where } p \in [0,1].
In practical terms, given :math:`(p, v)` in :math:`[0, 1]^2`,
:math:`h^{-1}(p \mid v)` solves for the :math:`u` satisfying
:math:`p = \partial C(u, v)/\partial v`.
Parameters
----------
X : array-like of shape (n_observations, 2)
An array of bivariate inputs `(p, v)`, each in the interval `[0, 1]`.
- The first column `p` corresponds to the value of the h-function.
- The second column `v` is the conditioning variable.
first_margin : bool, default=False
If True, compute the inverse partial derivative with respect to the first
margin `u`; otherwise, compute the inverse partial derivative with respect
to the second margin `v`.
Returns
-------
u : ndarray of shape (n_observations,)
A 1D-array of length `n_observations`, where each element is the computed
:math:`u = h^{-1}(p \mid v)` for the corresponding pair in `X`.
References
----------
.. [1] "Multivariate Models and Dependence Concepts", Joe, H. (1997)
.. [2] "An Introduction to Copulas", Nelsen, R. B. (2006)
"""
skv.check_is_fitted(self)
X = self._validate_X(X, reset=False)
X = _apply_margin_swap(X, first_margin=first_margin)
p_inv = sp.stdtrit(self.dof_ + 1, X[:, 0])
v_inv = sp.stdtrit(self.dof_, X[:, 1])
u_inv = (
p_inv
* np.sqrt((self.dof_ + v_inv**2) / (self.dof_ + 1) * (1 - self.rho_**2))
+ self.rho_ * v_inv
)
u = sp.stdtr(self.dof_, u_inv)
return u
[docs]
def score_samples(self, X: ArrayLike) -> FloatArray:
"""Compute the log-likelihood of each sample (log-pdf) under the model.
Parameters
----------
X : array-like of shape (n_observations, 2)
An array of bivariate inputs `(u, v)` where each row represents a
bivariate observation. Both `u` and `v` must be in the interval `[0, 1]`,
having been transformed to uniform marginals.
Returns
-------
density : ndarray of shape (n_observations,)
The log-likelihood of each sample under the fitted copula.
"""
skv.check_is_fitted(self)
X = self._validate_X(X, reset=False)
log_density = _sample_scores(X=X, rho=self.rho_, dof=self.dof_)
return log_density
@property
def lower_tail_dependence(self) -> float:
"""Theoretical lower tail dependence coefficient."""
skv.check_is_fitted(self)
arg = -np.sqrt((self.dof_ + 1) * (1 - self.rho_) / (1 + self.rho_))
return 2 * sp.stdtr(self.dof_ + 1, arg)
@property
def upper_tail_dependence(self) -> float:
"""Theoretical upper tail dependence coefficient."""
return self.lower_tail_dependence
@property
def fitted_repr(self) -> str:
"""String representation of the fitted copula."""
return f"{self.__class__.__name__}(rho={self.rho_:0.3f}, dof={self.dof_:0.2f})"
def _student_t_copula_cdf(X: FloatArray, rho: float, dof: float) -> FloatArray:
"""Evaluate the Student-t copula CDF using the angular integral in `cdf`."""
smaller = np.min(X, axis=1)
larger = np.max(X, axis=1)
lower = np.maximum(smaller - (1 - larger), 0)
cdf = smaller.copy()
interior = np.flatnonzero((smaller > 0) & (larger < 1))
quantiles = sp.stdtrit(dof, X[interior])
upper_angle = math.acos(-rho)
scale = 2 * math.sqrt(dof)
exponent = -float(dof) / 2
def integrand(angle: float, a: float, b: float) -> float:
"""Evaluate the nonnegative angular integrand in logarithmic form."""
sin_term = a / math.sin(angle / 2)
cos_term = b / math.cos(angle / 2)
return math.exp(
exponent * math.log1p(sin_term * sin_term + cos_term * cos_term)
) / (2 * math.pi)
for i, (x, y) in zip(interior, quantiles, strict=True):
# Use Python scalars to avoid NumPy overhead in the quadrature callback.
a = float((x + y) / scale)
b = float((x - y) / scale)
# Near x = -y, the integrand changes on a narrow scale at angle zero.
# Geometric breakpoints keep adaptive quadrature from missing this layer.
width = 2 * math.atan2(abs(a), math.hypot(1, b))
points = []
if width > 0:
while width < upper_angle:
points.append(width)
width *= 8
# Scaling the absolute tolerance preserves accuracy when C(q, q) is
# divided by q to compute tail concentration.
result = si.quad(
integrand,
0,
upper_angle,
args=(a, b),
points=points,
epsabs=1e-8 * smaller[i],
epsrel=1e-8,
limit=200,
full_output=True,
)
if len(result) > 3:
raise RuntimeError(f"Student-t copula CDF integration failed: {result[3]}")
# Enforce the Frechet bounds against floating-point rounding only after
# checking that the integration converged.
cdf[i] = min(max(lower[i] + result[0], lower[i]), smaller[i])
return cdf
def _neg_log_likelihood(dof: float, rho: float, X: FloatArray) -> float:
"""Negative log-likelihood function for optimization.
Parameters
----------
X : array-like of shape (n_observations, 2)
An array of bivariate inputs `(u, v)` where each row represents a
bivariate observation. Both `u` and `v` must be in the interval `[0, 1]`,
having been transformed to uniform marginals.
rho : float
Correlation copula parameter.
dof : float
Degree of freedom copula parameter.
Returns
-------
value : float
The negative log-likelihood value.
"""
return -np.sum(_sample_scores(X=X, rho=rho, dof=dof))
def _sample_scores(X: FloatArray, rho: float, dof: float) -> FloatArray:
"""Compute the log-likelihood of each sample (log-pdf) under the bivariate
Gaussian copula model.
Parameters
----------
X : array-like of shape (n_observations, 2)
An array of bivariate inputs `(u, v)` where each row represents a
bivariate observation. Both `u` and `v` must be in the interval `[0, 1]`,
having been transformed to uniform marginals.
rho : float
Gaussian copula parameter.
Returns
-------
density : ndarray of shape (n_observations,)
The log-likelihood of each sample under the fitted copula.
Raises
------
ValueError
If rho is not in (-1, 1) or dof is not positive.
"""
if not (-1.0 <= rho <= 1.0):
raise ValueError("rho must be between -1 and 1.")
if not 1.0 <= dof <= 50:
raise ValueError("Degrees of freedom `dof` must be between 1 and 50.")
# Inverse CDF (ppf) using stdtrit for better performance
x, y = sp.stdtrit(dof, X).T
a = 1.0 - rho**2
log_density = (
sp.gammaln((dof + 2.0) / 2.0)
+ sp.gammaln(dof / 2.0)
- 2.0 * sp.gammaln((dof + 1.0) / 2.0)
- np.log(a) / 2
+ (dof + 1.0) / 2.0 * (np.log1p(x**2 / dof) + np.log1p(y**2 / dof))
- (dof + 2.0) / 2.0 * np.log1p((x**2 - 2 * rho * x * y + y**2) / a / dof)
)
return log_density
