quanttoolbox.copula

copula.families

Python alternatives

Hybrid: statsmodels.distributions.copula.api (already a dependency of this package) implements Gaussian, Student-t, Clayton, Frank, and Gumbel copulas and was used to numerically verify clayton_cdf/clayton_pdf, frank_cdf/frank_pdf, gumbel_cdf/gumbel_pdf, and gaussian_copula_cdf/gaussian_copula_pdf to full float64 precision. It is not wrapped directly: statsmodels' classes fix one theta/corr per instance, while this port (like the original MATLAB) vectorizes theta per observation; and StudentTCopula.cdf() raises NotImplementedError in statsmodels (no closed form is implemented there), so student_copula_cdf had to be self-implemented via scipy.stats.multivariate_t/a bivariate wrapper regardless. gaussian_copula_cdf/_pdf and student_copula_cdf/_pdf reuse quanttoolbox.stats.multivariate.bvn_cdf/bvn_pdf/bvt_cdf for the bivariate case and quanttoolbox.stats.distributions.mvn_cdf for the general Gaussian case, instead of a third from-scratch implementation of bivariate/multivariate-normal-CDF machinery already present twice elsewhere in this package. The remaining 13 families (AMH, Gumbel-Barnett, Galambos, Husler-Reiss, Plackett, FGM, Cubic, logistic-Gumbel, Marshall-Olkin, Sloane, nested Gumbel, and the Fréchet-Hoeffding/independence bounds) have no equivalent in any Python copula library surveyed (statsmodels, copulas, pyvinecopulib) -- keep.

The original MATLAB source's own cdfCopulaGumbel3.m/pdfCopulaGumbel3.m pair (a 3-variable nested/hierarchical Gumbel copula) has a genuine bug: the shipped PDF formula does not match d³C/du1du2du3 of its own CDF, verified independently via both a 3-D central finite difference and exact sympy symbolic differentiation (both agree with each other at ≈1.2096 for theta1=1.5, theta2=3.0, u=(0.3,0.5,0.7); the original formula gives ≈1.0064). nested_gumbel_cdf is ported; no PDF is shipped for it, rather than risk a second, differently-wrong hand derivation. Separately, several families (AMH, Husler-Reiss, Marshall-Olkin, FGM, Sloane) never had a PDF in the original at all -- no PDF is fabricated for those either.

quanttoolbox.copula.families

CDFs and PDFs for the HSF toolbox's copula families: the Fréchet-Hoeffding bounds and independence copula, the Gaussian and Student-t (elliptical) copulas, five one-parameter Archimedean copulas (Clayton, Frank, Gumbel, AMH, Gumbel-Barnett), two extreme-value copulas (Galambos, Husler-Reiss), and five other named families (Plackett, FGM, Cubic, logistic-Gumbel, Marshall-Olkin, Sloane, nested/hierarchical Gumbel).

Ported from HSF toolbox copula/{cdfCopula*,pdfCopula*, cdfConditionalCopula*,cdfmvn,cdfSloaneCopula,contourCopula*, singularCopula*}.m (23 families across ~50 files).

Architecture -- avoiding the original's per-family boilerplate:

• The four Fréchet-Hoeffding-bound / independence functions each existed twice in the original (cdfCopulaUpper.m for n dimensions vs. cdfCopulaUpper2.m for 2; same for Lower/Product) purely because MATLAB has no convenient "2-column array or two scalars" polymorphism. comonotonicity_cdf/independence_cdf here take a single u array of shape (n_obs, n_dim) and work for any n_dim including 2, collapsing each pair of .m files into one function. countermonotonicity_cdf stays bivariate-only (it's the only one of the four that is a genuine copula, i.e. has uniform margins, exactly when n_dim = 2; frechet_lower_bound is the general-n non-copula bound cdfCopulaLower.m computes).
• gaussian_copula_cdf/_pdf and student_copula_cdf/_pdf likewise collapse the n-dimensional (cdfCopulaNormal.m/cdfCopulaStudent.m) and bivariate (cdfCopulaNormal2.m/cdfCopulaStudent2.m) originals into one function each, special-casing n_dim = 2 to call quanttoolbox.stats.multivariate.bvn_cdf/bvt_cdf -- already a tested, scipy-backed, near-singular-correlation-robust implementation -- rather than rederiving bivariate-normal/Student CDF logic a second time inside this module. The general-n path calls quanttoolbox.stats.distributions.mvn_cdf (Gaussian) or scipy.stats.multivariate_t.cdf (Student, no existing n-dimensional wrapper to reuse).
• Archimedean/extreme-value/other named families (Clayton, Frank, Gumbel, AMH, Gumbel-Barnett, Galambos, Husler-Reiss, Plackett, FGM, Cubic, logistic-Gumbel, Marshall-Olkin, Sloane) are each only a few lines of closed-form numpy, so are transliterated directly rather than forced through a generic Archimedean-generator abstraction that would save little code while obscuring each family's actual (and structurally quite different) closed form. Clayton/Frank/Gumbel's CDF, PDF, and Kendall's tau were verified to match statsmodels.distributions.copula exactly (statsmodels is already a dependency of this package); they are kept as plain vectorized numpy functions here rather than thin wrappers around the statsmodels class-based API because the original (and the rest of this port) vectorizes theta itself per observation, which statsmodels' one-theta-per-instance design does not support.
• Every family the original ported a CDF for but never ported a PDF for (AMH, Husler-Reiss, Marshall-Olkin, FGM, Sloane) stays that way here -- no PDF is fabricated for families the original author evidently chose not to derive one for.
• cdfCopulaGumbel3.m (a 3-variable nested/hierarchical Gumbel copula) is ported as nested_gumbel_cdf, but its companion pdfCopulaGumbel3.m is not ported: independently re-deriving d^3 C / du1 du2 du3 both symbolically (via sympy) and via a 3-D central finite difference disagreed with the original formula (verified at theta1=1.5, theta2=3.0, u=(0.3, 0.5, 0.7): symbolic/numeric both give ~1.2096, the original formula gives ~1.0064) -- the original's PDF has a genuine bug. Given how easy it would be to introduce a different subtle error in a hand-derived replacement for this rarely-used trivariate extension, the safer choice is to ship the (verified-correct) CDF only.

amh_cdf(u1, u2, theta)

The Ali-Mikhail-Haq (AMH) copula CDF. No PDF is provided in the original -- see module docstring; amh_conditional_cdf (the partial derivative w.r.t. u1) is available for simulation.

Original: copula/cdfCopulaAMH.m

Source code in src/quanttoolbox/copula/families.py

def amh_cdf(
    u1: np.ndarray | float, u2: np.ndarray | float, theta: np.ndarray | float
) -> np.ndarray:
    """The Ali-Mikhail-Haq (AMH) copula CDF. No PDF is provided in the
    original -- see module docstring; `amh_conditional_cdf` (the partial
    derivative w.r.t. `u1`) is available for simulation.

    Original: copula/cdfCopulaAMH.m
    """
    u1 = np.asarray(u1, dtype=float)
    u2 = np.asarray(u2, dtype=float)
    theta = np.asarray(theta, dtype=float)
    return u1 * u2 / (1.0 - theta * (1.0 - u1) * (1.0 - u2))

amh_conditional_cdf(u1, u2, theta)

The AMH copula's conditional CDF Pr(U2 <= u2 | U1 = u1) (dC/du1), used for conditional simulation.

Original: copula/cdfConditionalCopulaAMH.m

Source code in src/quanttoolbox/copula/families.py

def amh_conditional_cdf(
    u1: np.ndarray | float, u2: np.ndarray | float, theta: np.ndarray | float
) -> np.ndarray:
    """The AMH copula's conditional CDF ``Pr(U2 <= u2 | U1 = u1)``
    (``dC/du1``), used for conditional simulation.

    Original: copula/cdfConditionalCopulaAMH.m
    """
    u1 = np.asarray(u1, dtype=float)
    u2 = np.asarray(u2, dtype=float)
    theta = np.asarray(theta, dtype=float)
    numerator = (1.0 - theta) * u2 + theta * u2**2
    denominator = (1.0 - theta * (1.0 - u1) * (1.0 - u2)) ** 2
    return numerator / denominator

clayton_cdf(u1, u2, theta)

The Clayton copula CDF (verified to match statsmodels.distributions.copula.archimedean.ClaytonCopula).

Original: copula/cdfCopulaClayton.m

Source code in src/quanttoolbox/copula/families.py

def clayton_cdf(
    u1: np.ndarray | float, u2: np.ndarray | float, theta: np.ndarray | float
) -> np.ndarray:
    """The Clayton copula CDF (verified to match
    `statsmodels.distributions.copula.archimedean.ClaytonCopula`).

    Original: copula/cdfCopulaClayton.m
    """
    u1 = np.asarray(u1, dtype=float)
    u2 = np.asarray(u2, dtype=float)
    theta = np.asarray(theta, dtype=float)
    return np.maximum(u1 ** (-theta) + u2 ** (-theta) - 1.0, 0.0) ** (-1.0 / theta)

clayton_pdf(u1, u2, theta)

The Clayton copula density (verified to match statsmodels).

Original: copula/pdfCopulaClayton.m

Source code in src/quanttoolbox/copula/families.py

def clayton_pdf(
    u1: np.ndarray | float, u2: np.ndarray | float, theta: np.ndarray | float
) -> np.ndarray:
    """The Clayton copula density (verified to match `statsmodels`).

    Original: copula/pdfCopulaClayton.m
    """
    u1 = np.asarray(u1, dtype=float)
    u2 = np.asarray(u2, dtype=float)
    theta = np.asarray(theta, dtype=float)
    base = np.maximum(u1 ** (-theta) + u2 ** (-theta) - 1.0, 0.0)
    pdf = (1.0 + theta) * ((u1 * u2) ** (-theta - 1.0)) * base ** (-(2.0 * theta + 1.0) / theta)
    return np.where(np.isnan(pdf), 0.0, pdf)

comonotonicity_cdf(u)

The comonotonicity (Fréchet-Hoeffding upper-bound, "M") copula: C(u) = min(u_1, ..., u_n). u has shape (n_obs, n_dim).

Original: copula/cdfCopulaUpper.m (n-dim) and cdfCopulaUpper2.m (bivariate) -- merged, see module docstring.

Source code in src/quanttoolbox/copula/families.py

def comonotonicity_cdf(u: np.ndarray) -> np.ndarray:
    """The comonotonicity (Fréchet-Hoeffding upper-bound, "M") copula:
    ``C(u) = min(u_1, ..., u_n)``. `u` has shape ``(n_obs, n_dim)``.

    Original: copula/cdfCopulaUpper.m (n-dim) and cdfCopulaUpper2.m
    (bivariate) -- merged, see module docstring.
    """
    u = np.asarray(u, dtype=float)
    return np.min(u, axis=1)

comonotonicity_support(u1)

The comonotonicity copula's singular support: u2 = u1.

Original: copula/singularCopulaUpper2.m

Source code in src/quanttoolbox/copula/families.py

def comonotonicity_support(u1: np.ndarray | float) -> np.ndarray:
    """The comonotonicity copula's singular support: ``u2 = u1``.

    Original: copula/singularCopulaUpper2.m
    """
    return np.asarray(u1, dtype=float)

countermonotonicity_cdf(u1, u2)

The countermonotonicity (Fréchet-Hoeffding lower-bound, "W") copula: C(u1, u2) = max(u1 + u2 - 1, 0). Only a genuine copula (uniform margins) in 2 dimensions -- see frechet_lower_bound for the general -n bound.

Original: copula/cdfCopulaLower2.m

Source code in src/quanttoolbox/copula/families.py

def countermonotonicity_cdf(u1: np.ndarray | float, u2: np.ndarray | float) -> np.ndarray:
    """The countermonotonicity (Fréchet-Hoeffding lower-bound, "W") copula:
    ``C(u1, u2) = max(u1 + u2 - 1, 0)``. Only a genuine copula (uniform
    margins) in 2 dimensions -- see `frechet_lower_bound` for the general
    -`n` bound.

    Original: copula/cdfCopulaLower2.m
    """
    u1 = np.asarray(u1, dtype=float)
    u2 = np.asarray(u2, dtype=float)
    return np.maximum(u1 + u2 - 1.0, 0.0)

countermonotonicity_support(u1)

The countermonotonicity copula's singular support: u2 = 1 - u1.

Original: copula/singularCopulaLower2.m

Source code in src/quanttoolbox/copula/families.py

def countermonotonicity_support(u1: np.ndarray | float) -> np.ndarray:
    """The countermonotonicity copula's singular support: ``u2 = 1 - u1``.

    Original: copula/singularCopulaLower2.m
    """
    return 1.0 - np.asarray(u1, dtype=float)

cubic_cdf(u1, u2, theta)

The cubic copula CDF.

Original: copula/cdfCopulaCubic.m

Source code in src/quanttoolbox/copula/families.py

def cubic_cdf(
    u1: np.ndarray | float, u2: np.ndarray | float, theta: np.ndarray | float
) -> np.ndarray:
    """The cubic copula CDF.

    Original: copula/cdfCopulaCubic.m
    """
    u1 = np.asarray(u1, dtype=float)
    u2 = np.asarray(u2, dtype=float)
    theta = np.asarray(theta, dtype=float)
    return u1 * u2 + theta * u1 * (u1 - 1.0) * (2.0 * u1 - 1.0) * u2 * (u2 - 1.0) * (2.0 * u2 - 1.0)

cubic_pdf(u1, u2, theta)

The cubic copula density.

Original: copula/pdfCopulaCubic.m

Source code in src/quanttoolbox/copula/families.py

def cubic_pdf(
    u1: np.ndarray | float, u2: np.ndarray | float, theta: np.ndarray | float
) -> np.ndarray:
    """The cubic copula density.

    Original: copula/pdfCopulaCubic.m
    """
    u1 = np.asarray(u1, dtype=float)
    u2 = np.asarray(u2, dtype=float)
    theta = np.asarray(theta, dtype=float)
    return 1.0 + theta * (6.0 * u1**2 - 6.0 * u1 + 1.0) * (6.0 * u2**2 - 6.0 * u2 + 1.0)

fgm_cdf(u1, u2, theta)

The Farlie-Gumbel-Morgenstern (FGM) copula CDF. No PDF is provided in the original -- see module docstring.

Original: copula/cdfCopulaFGM.m

Source code in src/quanttoolbox/copula/families.py

def fgm_cdf(
    u1: np.ndarray | float, u2: np.ndarray | float, theta: np.ndarray | float
) -> np.ndarray:
    """The Farlie-Gumbel-Morgenstern (FGM) copula CDF. No PDF is provided
    in the original -- see module docstring.

    Original: copula/cdfCopulaFGM.m
    """
    u1 = np.asarray(u1, dtype=float)
    u2 = np.asarray(u2, dtype=float)
    theta = np.asarray(theta, dtype=float)
    return u1 * u2 * (1.0 + theta * (1.0 - u1) * (1.0 - u2))

frank_cdf(u1, u2, theta)

The Frank copula CDF (verified to match statsmodels).

Original: copula/cdfCopulaFrank.m

Source code in src/quanttoolbox/copula/families.py

def frank_cdf(
    u1: np.ndarray | float, u2: np.ndarray | float, theta: np.ndarray | float
) -> np.ndarray:
    """The Frank copula CDF (verified to match `statsmodels`).

    Original: copula/cdfCopulaFrank.m
    """
    u1 = np.asarray(u1, dtype=float)
    u2 = np.asarray(u2, dtype=float)
    theta = np.asarray(theta, dtype=float)
    a = np.exp(-theta * u1) - 1.0
    b = np.exp(-theta * u2) - 1.0
    return -np.log(1.0 + a * b / (np.exp(-theta) - 1.0)) / theta

frank_contour(u1, alpha, theta)

u2 such that the Frank copula's CDF equals alpha along the curve through u1 -- a level-set/contour curve used for plotting. Returns nan where the contour has no (real-valued, in-[0,1]) solution.

Original: copula/contourCopulaFrank.m

Source code in src/quanttoolbox/copula/families.py

def frank_contour(
    u1: np.ndarray | float, alpha: np.ndarray | float, theta: np.ndarray | float
) -> np.ndarray:
    """`u2` such that the Frank copula's CDF equals `alpha` along the
    curve through `u1` -- a level-set/contour curve used for plotting.
    Returns `nan` where the contour has no (real-valued, in-``[0,1]``)
    solution.

    Original: copula/contourCopulaFrank.m
    """
    u1 = np.asarray(u1, dtype=float)
    alpha = np.asarray(alpha, dtype=float)
    theta = np.asarray(theta, dtype=float)

    u2 = (
        -np.log(
            1.0
            + (np.exp(-alpha * theta) - 1.0) * (np.exp(-theta) - 1.0) / (np.exp(-theta * u1) - 1.0)
        )
        / theta
    )
    u2 = np.where(np.abs(np.imag(u2)) > 0.01, np.nan, np.real(u2))
    return np.where((u2 < 0.0) | (u2 > 1.0), np.nan, u2)

frank_level_curve(u1, t, theta)

u2(u1) such that the Frank copula C(u1, u2; theta) = t, from the closed-form inversion of the Frank copula CDF, evaluated on the outer-product grid of u1 (rows) against theta (columns) -- theta may be a vector of scenarios, broadcasting against u1 in one call, unlike frank_contour above which broadcasts u1 and theta elementwise. Since any copula satisfies C(u1, u2) <= min(u1, u2) <= u1 (its value at u2 = 1 is exactly u1, its maximum over u2), a solution u2 in [0, 1] only exists for u1 >= t; below that the log argument goes negative and this returns nan for theta > 0 (matching the level curve's true domain) -- but, unlike frank_contour, does not additionally clamp results outside [0, 1] to nan, so for theta < 0 (or theta near 0) and u1 < t it can return an out-of-domain value greater than 1 rather than nan; this is carried over unchanged from the source notebook, which only exercises u1 >= t.

This is the same closed-form as frank_contour immediately above (copula/contourCopulaFrank.m) -- not re-derived here -- reshaped for grid evaluation over a vector theta and without frank_contour's extra domain clamp, matching the source notebook's own implementation exactly. Promoted from HSF-Notebooks chapter 15a.

Source code in src/quanttoolbox/copula/families.py

def frank_level_curve(u1: np.ndarray, t: np.ndarray | float, theta: np.ndarray) -> np.ndarray:
    """``u2(u1)`` such that the Frank copula ``C(u1, u2; theta) = t``,
    from the closed-form inversion of the Frank copula CDF, evaluated on
    the outer-product grid of `u1` (rows) against `theta` (columns) --
    `theta` may be a vector of scenarios, broadcasting against `u1` in
    one call, unlike `frank_contour` above which broadcasts `u1` and
    `theta` elementwise. Since any copula satisfies ``C(u1, u2) <=
    min(u1, u2) <= u1`` (its value at ``u2 = 1`` is exactly `u1`, its
    maximum over `u2`), a solution ``u2`` in ``[0, 1]`` only exists for
    ``u1 >= t``; below that the log argument goes negative and this
    returns `nan` for `theta > 0` (matching the level curve's true
    domain) -- but, unlike `frank_contour`, does *not* additionally
    clamp results outside ``[0, 1]`` to `nan`, so for `theta < 0` (or
    `theta` near ``0``) and ``u1 < t`` it can return an out-of-domain
    value greater than ``1`` rather than `nan`; this is carried over
    unchanged from the source notebook, which only exercises `u1 >= t`.

    This is the same closed-form as `frank_contour` immediately above
    (`copula/contourCopulaFrank.m`) -- not re-derived here -- reshaped
    for grid evaluation over a vector `theta` and without `frank_contour`'s
    extra domain clamp, matching the source notebook's own
    implementation exactly. Promoted from HSF-Notebooks chapter 15a.
    """
    u1 = np.asarray(u1, dtype=float)[:, None]
    theta = np.asarray(theta, dtype=float)[None, :]
    num = (np.exp(-theta * t) - 1.0) * (np.exp(-theta) - 1.0)
    denom = np.exp(-theta * u1) - 1.0
    with np.errstate(invalid="ignore"):
        return -1.0 / theta * np.log(1.0 + num / denom)

frank_pdf(u1, u2, theta)

The Frank copula density (verified to match statsmodels).

Original: copula/pdfCopulaFrank.m

Source code in src/quanttoolbox/copula/families.py

def frank_pdf(
    u1: np.ndarray | float, u2: np.ndarray | float, theta: np.ndarray | float
) -> np.ndarray:
    """The Frank copula density (verified to match `statsmodels`).

    Original: copula/pdfCopulaFrank.m
    """
    u1 = np.asarray(u1, dtype=float)
    u2 = np.asarray(u2, dtype=float)
    theta = np.asarray(theta, dtype=float)
    eta = 1.0 - np.exp(-theta)
    v1 = 1.0 - np.exp(-theta * u1)
    v2 = 1.0 - np.exp(-theta * u2)
    return np.exp(-theta * (u1 + u2)) * theta * eta / (eta - v1 * v2) ** 2

frechet_lower_bound(u)

The Fréchet-Hoeffding lower bound for n_dim dimensions: max(u_1 + ... + u_n - n_dim + 1, 0). Not itself a copula for n_dim > 2 (no distribution attains it as its copula), but every copula is bounded below by it. u has shape (n_obs, n_dim).

Original: copula/cdfCopulaLower.m

Source code in src/quanttoolbox/copula/families.py

def frechet_lower_bound(u: np.ndarray) -> np.ndarray:
    """The Fréchet-Hoeffding lower bound for `n_dim` dimensions:
    ``max(u_1 + ... + u_n - n_dim + 1, 0)``. Not itself a copula for
    `n_dim > 2` (no distribution attains it as its copula), but every
    copula is bounded below by it. `u` has shape ``(n_obs, n_dim)``.

    Original: copula/cdfCopulaLower.m
    """
    u = np.asarray(u, dtype=float)
    n_dim = u.shape[1]
    return np.maximum(np.sum(u, axis=1) - n_dim + 1.0, 0.0)

galambos_cdf(u1, u2, theta)

The Galambos copula CDF.

Original: copula/cdfCopulaGalambos.m

Source code in src/quanttoolbox/copula/families.py

def galambos_cdf(
    u1: np.ndarray | float, u2: np.ndarray | float, theta: np.ndarray | float
) -> np.ndarray:
    """The Galambos copula CDF.

    Original: copula/cdfCopulaGalambos.m
    """
    u1 = np.asarray(u1, dtype=float)
    u2 = np.asarray(u2, dtype=float)
    theta = np.asarray(theta, dtype=float)
    u1t = -np.log(u1)
    u2t = -np.log(u2)
    return u1 * u2 * np.exp((u1t ** (-theta) + u2t ** (-theta)) ** (-1.0 / theta))

galambos_pdf(u1, u2, theta)

The Galambos copula density.

Original: copula/pdfCopulaGalambos.m

Source code in src/quanttoolbox/copula/families.py

def galambos_pdf(
    u1: np.ndarray | float, u2: np.ndarray | float, theta: np.ndarray | float
) -> np.ndarray:
    """The Galambos copula density.

    Original: copula/pdfCopulaGalambos.m
    """
    u1 = np.asarray(u1, dtype=float)
    u2 = np.asarray(u2, dtype=float)
    theta = np.asarray(theta, dtype=float)
    u1t = -np.log(u1)
    u2t = -np.log(u2)
    cdf = galambos_cdf(u1, u2, theta)
    p = (
        1.0
        - (u1t ** (-theta) + u2t ** (-theta)) ** (-1.0 / theta - 1.0)
        * (u1t ** (-theta - 1.0) + u2t ** (-theta - 1.0))
        + (u1t ** (-theta) + u2t ** (-theta)) ** (-1.0 / theta - 2.0)
        * (u1t * u2t) ** (-1.0 - theta)
        * (1.0 + theta + (u1t ** (-theta) + u2t ** (-theta)) ** (-1.0 / theta))
    )
    return cdf / (u1 * u2) * p

gaussian_copula_cdf(u, corr)

The Gaussian copula CDF: C(u) = Phi_corr(Phi^-1(u_1), ..., Phi^-1(u_n)), Phi_corr the correlation-corr multivariate normal CDF. u has shape (n_obs, n_dim); for n_dim = 2, delegates to quanttoolbox.stats.multivariate.bvn_cdf (see module docstring).

Original: copula/cdfCopulaNormal.m (n-dim) and cdfCopulaNormal2.m (bivariate) -- merged.

Source code in src/quanttoolbox/copula/families.py

def gaussian_copula_cdf(u: np.ndarray, corr: np.ndarray) -> np.ndarray:
    """The Gaussian copula CDF: ``C(u) = Phi_corr(Phi^-1(u_1), ...,
    Phi^-1(u_n))``, `Phi_corr` the correlation-`corr` multivariate normal
    CDF. `u` has shape ``(n_obs, n_dim)``; for `n_dim = 2`, delegates to
    `quanttoolbox.stats.multivariate.bvn_cdf` (see module docstring).

    Original: copula/cdfCopulaNormal.m (n-dim) and cdfCopulaNormal2.m
    (bivariate) -- merged.
    """
    u = np.asarray(u, dtype=float)
    corr = np.asarray(corr, dtype=float)
    x = normdist.ppf(u)
    if u.shape[1] == 2:
        return bvn_cdf(x[:, 0], x[:, 1], corr[0, 1])
    return mvn_cdf(x, np.zeros(u.shape[1]), corr)

gaussian_copula_conditional_cdf(u1, u2, rho)

The bivariate Gaussian copula's conditional CDF Pr(U1 <= u1 | U2 = u2), used for conditional simulation.

Original: copula/cdfConditionalCopulaNormal2.m

Source code in src/quanttoolbox/copula/families.py

def gaussian_copula_conditional_cdf(
    u1: np.ndarray | float, u2: np.ndarray | float, rho: np.ndarray | float
) -> np.ndarray:
    """The bivariate Gaussian copula's conditional CDF ``Pr(U1 <= u1 |
    U2 = u2)``, used for conditional simulation.

    Original: copula/cdfConditionalCopulaNormal2.m
    """
    u1 = np.asarray(u1, dtype=float)
    u2 = np.asarray(u2, dtype=float)
    rho = np.asarray(rho, dtype=float)
    return normdist.cdf((normdist.ppf(u1) - rho * normdist.ppf(u2)) / np.sqrt(1.0 - rho**2))

gaussian_copula_pdf(u, corr)

The Gaussian copula density.

Original: copula/pdfCopulaNormal.m (n-dim) and pdfCopulaNormal2.m (bivariate) -- merged.

Source code in src/quanttoolbox/copula/families.py

def gaussian_copula_pdf(u: np.ndarray, corr: np.ndarray) -> np.ndarray:
    """The Gaussian copula density.

    Original: copula/pdfCopulaNormal.m (n-dim) and pdfCopulaNormal2.m
    (bivariate) -- merged.
    """
    u = np.asarray(u, dtype=float)
    corr = np.asarray(corr, dtype=float)
    x = normdist.ppf(u)
    if u.shape[1] == 2:
        rho = corr[0, 1]
        return bvn_pdf(x[:, 0], x[:, 1], 0.0, 0.0, 1.0, 1.0, rho) / (
            normdist.pdf(x[:, 0]) * normdist.pdf(x[:, 1])
        )
    corr_inv = np.linalg.inv(corr) - np.eye(u.shape[1])
    quad_form = np.einsum("ij,jk,ik->i", x, corr_inv, x)
    return np.exp(-0.5 * quad_form) / np.sqrt(np.linalg.det(corr))

gaussian_copula_pdf_grid(u1, u2, rho)

Evaluate the bivariate Gaussian copula density on the full outer-product grid u1 x u2. gaussian_copula_pdf above takes paired (n_obs, 2) columns, not a broadcasting row/column pair, so this builds the grid explicitly with numpy.meshgrid, flattens it into pairs for gaussian_copula_pdf, and reshapes the result back to the grid's shape.

Not ported from the MATLAB HSF toolbox -- no .m file in hfs-archive implements this as a function (the toolbox instead relies on MATLAB's implicit column/row broadcasting, e.g. chap15_copula1.m's pdfCopulaNormal2(u1, u2, rho) with u1 a column and u2 a row, which numpy does not do for gaussian_copula_pdf's paired-columns API). Promoted from HSF-Notebooks chapter 15a.

Source code in src/quanttoolbox/copula/families.py

def gaussian_copula_pdf_grid(u1: np.ndarray, u2: np.ndarray, rho: float) -> np.ndarray:
    """Evaluate the bivariate Gaussian copula density on the full
    outer-product grid ``u1 x u2``. `gaussian_copula_pdf` above takes
    paired ``(n_obs, 2)`` columns, not a broadcasting row/column pair,
    so this builds the grid explicitly with `numpy.meshgrid`, flattens
    it into pairs for `gaussian_copula_pdf`, and reshapes the result
    back to the grid's shape.

    Not ported from the MATLAB HSF toolbox -- no `.m` file in
    `hfs-archive` implements this as a function (the toolbox instead
    relies on MATLAB's implicit column/row broadcasting, e.g.
    `chap15_copula1.m`'s ``pdfCopulaNormal2(u1, u2, rho)`` with `u1` a
    column and `u2` a row, which `numpy` does not do for
    `gaussian_copula_pdf`'s paired-columns API). Promoted from
    HSF-Notebooks chapter 15a.
    """
    u1 = np.asarray(u1, dtype=float)
    u2 = np.asarray(u2, dtype=float)
    grid_u1, grid_u2 = np.meshgrid(u1, u2, indexing="ij")
    u_pairs = np.column_stack([grid_u1.ravel(), grid_u2.ravel()])
    corr = np.array([[1.0, rho], [rho, 1.0]])
    z = gaussian_copula_pdf(u_pairs, corr)
    return z.reshape(grid_u1.shape)

gumbel_barnett_cdf(u1, u2, theta)

The Gumbel-Barnett copula CDF.

Original: copula/cdfCopulaGumbelBarnett.m

Source code in src/quanttoolbox/copula/families.py

def gumbel_barnett_cdf(
    u1: np.ndarray | float, u2: np.ndarray | float, theta: np.ndarray | float
) -> np.ndarray:
    """The Gumbel-Barnett copula CDF.

    Original: copula/cdfCopulaGumbelBarnett.m
    """
    u1 = np.asarray(u1, dtype=float)
    u2 = np.asarray(u2, dtype=float)
    theta = np.asarray(theta, dtype=float)
    return u1 * u2 * np.exp(-theta * np.log(u1) * np.log(u2))

gumbel_barnett_pdf(u1, u2, theta)

The Gumbel-Barnett copula density.

Original: copula/pdfCopulaGumbelBarnett.m

Source code in src/quanttoolbox/copula/families.py

def gumbel_barnett_pdf(
    u1: np.ndarray | float, u2: np.ndarray | float, theta: np.ndarray | float
) -> np.ndarray:
    """The Gumbel-Barnett copula density.

    Original: copula/pdfCopulaGumbelBarnett.m
    """
    u1 = np.asarray(u1, dtype=float)
    u2 = np.asarray(u2, dtype=float)
    theta = np.asarray(theta, dtype=float)
    log_u1 = np.log(u1)
    log_u2 = np.log(u2)
    return (1.0 - theta - theta * (log_u1 + log_u2) + theta**2 * log_u1 * log_u2) * np.exp(
        -theta * log_u1 * log_u2
    )

gumbel_cdf(u1, u2, theta)

The Gumbel copula CDF (verified to match statsmodels).

Original: copula/cdfCopulaGumbel.m

Source code in src/quanttoolbox/copula/families.py

def gumbel_cdf(
    u1: np.ndarray | float, u2: np.ndarray | float, theta: np.ndarray | float
) -> np.ndarray:
    """The Gumbel copula CDF (verified to match `statsmodels`).

    Original: copula/cdfCopulaGumbel.m
    """
    u1 = np.asarray(u1, dtype=float)
    u2 = np.asarray(u2, dtype=float)
    theta = np.asarray(theta, dtype=float)
    return np.exp(-(((-np.log(u1)) ** theta + (-np.log(u2)) ** theta) ** (1.0 / theta)))

gumbel_pdf(u1, u2, theta)

The Gumbel copula density (verified to match statsmodels).

Original: copula/pdfCopulaGumbel.m

Source code in src/quanttoolbox/copula/families.py

def gumbel_pdf(
    u1: np.ndarray | float, u2: np.ndarray | float, theta: np.ndarray | float
) -> np.ndarray:
    """The Gumbel copula density (verified to match `statsmodels`).

    Original: copula/pdfCopulaGumbel.m
    """
    u1 = np.asarray(u1, dtype=float)
    u2 = np.asarray(u2, dtype=float)
    theta = np.asarray(theta, dtype=float)
    u1t = -np.log(u1)
    u2t = -np.log(u2)
    w = u1t**theta + u2t**theta
    pdf = (
        ((u1t * u2t) ** (theta - 1.0))
        * (w ** (1.0 / theta) + theta - 1.0)
        / (w ** (2.0 - 1.0 / theta))
        / (u1 * u2)
    )
    return pdf * gumbel_cdf(u1, u2, theta)

husler_reiss_cdf(u1, u2, theta)

The Husler-Reiss copula CDF. No PDF is provided in the original -- see module docstring.

Original: copula/cdfCopulaHuslerReiss.m

Source code in src/quanttoolbox/copula/families.py

def husler_reiss_cdf(
    u1: np.ndarray | float, u2: np.ndarray | float, theta: np.ndarray | float
) -> np.ndarray:
    """The Husler-Reiss copula CDF. No PDF is provided in the original --
    see module docstring.

    Original: copula/cdfCopulaHuslerReiss.m
    """
    u1 = np.asarray(u1, dtype=float)
    u2 = np.asarray(u2, dtype=float)
    theta = np.asarray(theta, dtype=float)
    # `+ 0.0` normalizes the -0.0 that `-np.log(1.0)` produces (IEEE 754
    # negation of +0.0) back to +0.0 -- left as -0.0, `u2t / u1t` below
    # would evaluate to -inf instead of +inf at the u1 == 1 boundary, and
    # `log(-inf)` is NaN instead of the correct +inf.
    u1t = -np.log(u1) + 0.0
    u2t = -np.log(u2) + 0.0
    phi1 = 1.0 / theta + 0.5 * theta * np.log(u1t / u2t)
    phi2 = 1.0 / theta + 0.5 * theta * np.log(u2t / u1t)
    return np.exp(-u1t * normdist.cdf(phi1) - u2t * normdist.cdf(phi2))

independence_cdf(u)

The independence (product) copula: C(u) = u_1 * u_2 * ... * u_n. u has shape (n_obs, n_dim).

Original: copula/cdfCopulaProduct.m (n-dim) and cdfCopulaProduct2.m (bivariate) -- merged, see module docstring.

Source code in src/quanttoolbox/copula/families.py

def independence_cdf(u: np.ndarray) -> np.ndarray:
    """The independence (product) copula: ``C(u) = u_1 * u_2 * ... *
    u_n``. `u` has shape ``(n_obs, n_dim)``.

    Original: copula/cdfCopulaProduct.m (n-dim) and cdfCopulaProduct2.m
    (bivariate) -- merged, see module docstring.
    """
    u = np.asarray(u, dtype=float)
    return np.prod(u, axis=1)

logistic_gumbel_cdf(u1, u2)

The (parameter-free) logistic-Gumbel copula CDF.

Original: copula/cdfCopulaLogisticGumbel.m

Source code in src/quanttoolbox/copula/families.py

def logistic_gumbel_cdf(u1: np.ndarray | float, u2: np.ndarray | float) -> np.ndarray:
    """The (parameter-free) logistic-Gumbel copula CDF.

    Original: copula/cdfCopulaLogisticGumbel.m
    """
    u1 = np.asarray(u1, dtype=float)
    u2 = np.asarray(u2, dtype=float)
    return u1 * u2 / (u1 + u2 - u1 * u2)

logistic_gumbel_contour(u1, alpha)

u2 such that the logistic-Gumbel copula's CDF equals alpha along the curve through u1, i.e. a level-set/contour curve used for plotting: values of u1 < alpha fall outside the valid contour and return nan.

Original: copula/contourCopulaLogisticGumbel.m

Source code in src/quanttoolbox/copula/families.py

def logistic_gumbel_contour(u1: np.ndarray | float, alpha: np.ndarray | float) -> np.ndarray:
    """`u2` such that the logistic-Gumbel copula's CDF equals `alpha`
    along the curve through `u1`, i.e. a level-set/contour curve used for
    plotting: values of `u1 < alpha` fall outside the valid contour and
    return `nan`.

    Original: copula/contourCopulaLogisticGumbel.m
    """
    u1 = np.asarray(u1, dtype=float)
    alpha = np.asarray(alpha, dtype=float)
    u1 = np.where(u1 < alpha, np.nan, u1)
    return alpha * u1 / (u1 + alpha * u1 - alpha)

logistic_gumbel_pdf(u1, u2)

The (parameter-free) logistic-Gumbel copula density.

Original: copula/pdfCopulaLogisticGumbel.m

Source code in src/quanttoolbox/copula/families.py

def logistic_gumbel_pdf(u1: np.ndarray | float, u2: np.ndarray | float) -> np.ndarray:
    """The (parameter-free) logistic-Gumbel copula density.

    Original: copula/pdfCopulaLogisticGumbel.m
    """
    u1 = np.asarray(u1, dtype=float)
    u2 = np.asarray(u2, dtype=float)
    return 2.0 * u1 * u2 / (u1 + u2 - u1 * u2) ** 3

marshall_olkin_cdf(u1, u2, theta1, theta2)

The Marshall-Olkin copula CDF. Has a singular component along u1**theta1 == u2**theta2 (see marshall_olkin_singular_support); no PDF is provided in the original.

Original: copula/cdfCopulaMarshallOlkin.m

Source code in src/quanttoolbox/copula/families.py

def marshall_olkin_cdf(
    u1: np.ndarray | float,
    u2: np.ndarray | float,
    theta1: np.ndarray | float,
    theta2: np.ndarray | float,
) -> np.ndarray:
    """The Marshall-Olkin copula CDF. Has a singular component along
    ``u1**theta1 == u2**theta2`` (see `marshall_olkin_singular_support`);
    no PDF is provided in the original.

    Original: copula/cdfCopulaMarshallOlkin.m
    """
    u1 = np.asarray(u1, dtype=float)
    u2 = np.asarray(u2, dtype=float)
    theta1 = np.asarray(theta1, dtype=float)
    theta2 = np.asarray(theta2, dtype=float)
    cond = u1**theta1 > u2**theta2
    return (
        u1 ** (1.0 - theta1)
        * u2 ** (1.0 - theta2)
        * (u1**theta1 + np.where(cond, u2**theta2 - u1**theta1, 0.0))
    )

marshall_olkin_singular_support(u1, theta1, theta2)

The Marshall-Olkin copula's singular support: u2 = u1 ** (theta1 / theta2).

Original: copula/singularCopulaMarshallOlkin.m

Source code in src/quanttoolbox/copula/families.py

def marshall_olkin_singular_support(
    u1: np.ndarray | float, theta1: float, theta2: float
) -> np.ndarray:
    """The Marshall-Olkin copula's singular support: ``u2 = u1 **
    (theta1 / theta2)``.

    Original: copula/singularCopulaMarshallOlkin.m
    """
    u1 = np.asarray(u1, dtype=float)
    return u1 ** (theta1 / theta2)

nested_gumbel_cdf(u1, u2, u3, theta1, theta2)

A 3-variable nested (hierarchical) Gumbel copula: (u1, u2) are coupled at strength theta2 first, then combined with u3 at strength theta1 (requires theta2 >= theta1 >= 1 for a valid copula). No PDF is provided -- see module docstring.

Original: copula/cdfCopulaGumbel3.m

Source code in src/quanttoolbox/copula/families.py

def nested_gumbel_cdf(
    u1: np.ndarray | float,
    u2: np.ndarray | float,
    u3: np.ndarray | float,
    theta1: np.ndarray | float,
    theta2: np.ndarray | float,
) -> np.ndarray:
    """A 3-variable nested (hierarchical) Gumbel copula: (u1, u2) are
    coupled at strength `theta2` first, then combined with `u3` at
    strength `theta1` (requires ``theta2 >= theta1 >= 1`` for a valid
    copula). No PDF is provided -- see module docstring.

    Original: copula/cdfCopulaGumbel3.m
    """
    u1t = -np.log(np.asarray(u1, dtype=float))
    u2t = -np.log(np.asarray(u2, dtype=float))
    u3t = -np.log(np.asarray(u3, dtype=float))
    theta1 = np.asarray(theta1, dtype=float)
    theta2 = np.asarray(theta2, dtype=float)

    inner = u1t**theta2 + u2t**theta2
    return np.exp(-((inner ** (theta1 / theta2) + u3t**theta1) ** (1.0 / theta1)))

plackett_cdf(u1, u2, theta)

The Plackett copula CDF.

Original: copula/cdfCopulaPlackett.m

Source code in src/quanttoolbox/copula/families.py

def plackett_cdf(
    u1: np.ndarray | float, u2: np.ndarray | float, theta: np.ndarray | float
) -> np.ndarray:
    """The Plackett copula CDF.

    Original: copula/cdfCopulaPlackett.m
    """
    u1 = np.asarray(u1, dtype=float)
    u2 = np.asarray(u2, dtype=float)
    theta = np.asarray(theta, dtype=float)
    eta = theta - 1.0
    w = 1.0 + eta * (u1 + u2)
    return 0.5 * (w - np.sqrt(w**2 - 4.0 * theta * eta * u1 * u2)) / eta

plackett_pdf(u1, u2, theta)

The Plackett copula density.

Original: copula/pdfCopulaPlackett.m

Source code in src/quanttoolbox/copula/families.py

def plackett_pdf(
    u1: np.ndarray | float, u2: np.ndarray | float, theta: np.ndarray | float
) -> np.ndarray:
    """The Plackett copula density.

    Original: copula/pdfCopulaPlackett.m
    """
    u1 = np.asarray(u1, dtype=float)
    u2 = np.asarray(u2, dtype=float)
    theta = np.asarray(theta, dtype=float)
    eta = theta - 1.0
    w = 1.0 + eta * (u1 + u2)
    return theta * (w - 2.0 * eta * u1 * u2) / (w**2 - 4.0 * theta * eta * u1 * u2) ** 1.5

sloane_cdf(u1, u2, rho)

The Sloane copula CDF. No PDF is provided in the original.

Original: copula/cdfSloaneCopula.m

Source code in src/quanttoolbox/copula/families.py

def sloane_cdf(
    u1: np.ndarray | float, u2: np.ndarray | float, rho: np.ndarray | float
) -> np.ndarray:
    """The Sloane copula CDF. No PDF is provided in the original.

    Original: copula/cdfSloaneCopula.m
    """
    u1 = np.asarray(u1, dtype=float)
    u2 = np.asarray(u2, dtype=float)
    rho = np.asarray(rho, dtype=float)

    u1v = np.arccosh(1.0 / u1**2)
    u2v = np.arccosh(1.0 / u2**2)
    is_u1_smaller = u1v <= u2v
    zeta = np.where(is_u1_smaller, u1v, u2v)
    xi = np.abs(u1v - u2v)

    c1 = np.cosh(xi) * np.cosh(zeta * np.sqrt(1.0 + rho)) * np.cosh(zeta * np.sqrt(1.0 - rho))
    c2 = np.sinh(xi) * np.cosh(zeta * np.sqrt(1.0 - rho)) * np.sinh(zeta * np.sqrt(1.0 + rho))
    c3 = np.sinh(xi) * np.cosh(zeta * np.sqrt(1.0 + rho)) * np.sinh(zeta * np.sqrt(1.0 - rho))

    c = c1 + 0.5 * np.sqrt(1.0 + rho) * c2 + 0.5 * np.sqrt(1.0 - rho) * c3
    c = 1.0 / np.sqrt(c)
    return np.where(np.isreal(c), np.real(c), 0.0)

student_copula_cdf(u, corr, nu)

The Student-t copula CDF with nu degrees of freedom. u has shape (n_obs, n_dim); for n_dim = 2, delegates to quanttoolbox.stats.multivariate.bvt_cdf.

Original: copula/cdfCopulaStudent.m (n-dim) and cdfCopulaStudent2.m (bivariate, looping over a scenario grid of rho/nu pairs) -- merged: sweep scenarios with a plain Python loop over corr/nu at the call site instead of baking scenario-looping into the function.

Source code in src/quanttoolbox/copula/families.py

def student_copula_cdf(u: np.ndarray, corr: np.ndarray, nu: float) -> np.ndarray:
    """The Student-t copula CDF with `nu` degrees of freedom. `u` has
    shape ``(n_obs, n_dim)``; for `n_dim = 2`, delegates to
    `quanttoolbox.stats.multivariate.bvt_cdf`.

    Original: copula/cdfCopulaStudent.m (n-dim) and cdfCopulaStudent2.m
    (bivariate, looping over a scenario grid of `rho`/`nu` pairs) --
    merged: sweep scenarios with a plain Python loop over `corr`/`nu` at
    the call site instead of baking scenario-looping into the function.
    """
    u = np.asarray(u, dtype=float)
    corr = np.asarray(corr, dtype=float)
    from scipy.stats import t as tdist

    x = tdist.ppf(u, nu)
    if u.shape[1] == 2:
        return bvt_cdf(x[:, 0], x[:, 1], corr[0, 1], nu)
    return multivariate_t.cdf(x, loc=np.zeros(u.shape[1]), shape=corr, df=float(nu))

student_copula_pdf(u, corr, nu)

The Student-t copula density with nu degrees of freedom.

Original: copula/pdfCopulaStudent.m (n-dim) and pdfCopulaStudent2.m (bivariate) -- merged.

Source code in src/quanttoolbox/copula/families.py

def student_copula_pdf(u: np.ndarray, corr: np.ndarray, nu: float) -> np.ndarray:
    """The Student-t copula density with `nu` degrees of freedom.

    Original: copula/pdfCopulaStudent.m (n-dim) and pdfCopulaStudent2.m
    (bivariate) -- merged.
    """
    from scipy.special import gammaln
    from scipy.stats import t as tdist

    u = np.asarray(u, dtype=float)
    corr = np.asarray(corr, dtype=float)
    n_dim = u.shape[1]
    x = tdist.ppf(u, nu)

    corr_inv = np.linalg.inv(corr)
    corr_det = np.linalg.det(corr)
    quad_form = np.einsum("ij,jk,ik->i", x, corr_inv, x)

    log_num = gammaln((nu + n_dim) / 2.0) + (n_dim - 1) * gammaln(nu / 2.0)
    log_den = n_dim * gammaln((nu + 1.0) / 2.0) + 0.5 * np.log(corr_det)
    log_pdf = (
        -0.5 * (nu + n_dim) * np.log1p(quad_form / nu)
        + 0.5 * (nu + 1.0) * np.sum(np.log1p(x**2 / nu), axis=1)
        + log_num
        - log_den
    )
    pdf = np.exp(log_pdf)
    return np.where(np.isnan(pdf), 0.0, pdf)

copula.dependence

Python alternatives

Keep -- Kendall's tau and Spearman's rho for the families above, plus the Debye/dilogarithm special functions their closed forms need and the empirical dependogram (pseudo-observations via marginal ranks). clayton_tau/frank_tau/gumbel_tau/gaussian_tau were verified against statsmodels' .tau() methods. The original SpearmanCopula.m is itself a generic double-integral estimator (rho = 12 * integral2(C,0,1,0,1) - 3, for any copula CDF C) that the original only ever plugged Clayton and Gumbel into (both lack closed-form Spearman's rho); spearman_rho_numeric here (via scipy.integrate.dblquad) keeps that genericity explicit and usable for any family in families.py, not just the two the original wired up -- clayton_rho/gumbel_rho are one-line instantiations of it.

quanttoolbox.copula.dependence

Kendall's tau and Spearman's rho for the copula families in families.py, the two special functions (Debye, dilogarithm) their closed forms need, and the empirical dependogram (pseudo-observations via marginal ranks).

Ported from HSF toolbox copula/{KendallCopula*,SpearmanCopula*, DebyeFunction,diLogFunction,dependogram}.m.

Architecture -- avoiding the original's per-family boilerplate:

• SpearmanCopula.m is already, in the original, a generic estimator: rho = 12 * integral2(C, 0, 1, 0, 1) - 3 for any copula CDF C. SpearmanCopulaClayton.m/SpearmanCopulaGumbel.m are just this generic function called with cdfCopulaClayton/cdfCopulaGumbel plugged in (looping over a theta array) -- because Clayton and Gumbel have no closed-form Spearman's rho. That generic estimator is spearman_rho_ numeric here (scipy.integrate.dblquad in place of integral2); clayton_rho/gumbel_rho are one-line calls to it, and it works unchanged for any other family in families.py that lacks a closed form (AMH, Plackett, ... -- pass e.g. lambda u1, u2: families. plackett_cdf(u1, u2, theta)), not just the two the original wired it up for.
• Frank and Gaussian do have closed forms (SpearmanCopulaFrank.m, SpearmanCopulaNormal.m), so frank_rho/gaussian_rho use those directly instead of the (much slower) numeric integral.

clayton_rho(theta)

Spearman's rho of the Clayton copula (no closed form -- numeric integration via spearman_rho_numeric).

Original: copula/SpearmanCopulaClayton.m

Source code in src/quanttoolbox/copula/dependence.py

def clayton_rho(theta: float) -> float:
    """Spearman's rho of the Clayton copula (no closed form -- numeric
    integration via `spearman_rho_numeric`).

    Original: copula/SpearmanCopulaClayton.m
    """
    from quanttoolbox.copula.families import clayton_cdf

    return spearman_rho_numeric(lambda u1, u2: float(clayton_cdf(u1, u2, theta)))

clayton_tau(theta)

Kendall's tau of the Clayton copula: theta / (theta + 2) (verified to match statsmodels).

Original: copula/KendallCopulaClayton.m

Source code in src/quanttoolbox/copula/dependence.py

def clayton_tau(theta: np.ndarray | float) -> np.ndarray:
    """Kendall's tau of the Clayton copula: ``theta / (theta + 2)``
    (verified to match `statsmodels`).

    Original: copula/KendallCopulaClayton.m
    """
    theta = np.asarray(theta, dtype=float)
    return theta / (theta + 2.0)

debye_function(x, k)

The Debye function D_k(x) = (k / x^k) * integral_0^x t^k / (e^t - 1) dt, extended to x < 0 via the reflection identity D_k(-x) = D_k(x) + k*x / (k+1). Used by frank_tau/frank_rho.

Original: copula/DebyeFunction.m

Source code in src/quanttoolbox/copula/dependence.py

def debye_function(x: np.ndarray | float, k: int) -> np.ndarray:
    """The Debye function ``D_k(x) = (k / x^k) * integral_0^x t^k / (e^t
    - 1) dt``, extended to `x < 0` via the reflection identity ``D_k(-x)
    = D_k(x) + k*x / (k+1)``. Used by `frank_tau`/`frank_rho`.

    Original: copula/DebyeFunction.m
    """
    x_arr = np.atleast_1d(np.asarray(x, dtype=float))
    out = np.empty_like(x_arr)
    for i, xx in enumerate(x_arr.ravel()):
        is_nonneg = xx >= 0.0
        xa = abs(xx)
        if xa == 0.0:
            out.ravel()[i] = 1.0 / (k + 1.0)
            continue
        integral, _ = quad(lambda t: t**k / np.expm1(t), 0.0, xa)
        value = (k / xa**k) * integral
        if not is_nonneg:
            value = value + k * xa / (1.0 + k)
        out.ravel()[i] = value
    return out if np.asarray(x).ndim > 0 else out[0]

dependogram(data)

Pseudo-observations of the empirical copula: each column of data (shape (n_obs, n_dim)) replaced by its marginal ranks (ties averaged), scaled to (0, 1) by rank / (n_obs + 1).

Original: copula/dependogram.m

Source code in src/quanttoolbox/copula/dependence.py

def dependogram(data: np.ndarray) -> np.ndarray:
    """Pseudo-observations of the empirical copula: each column of
    `data` (shape ``(n_obs, n_dim)``) replaced by its marginal ranks
    (ties averaged), scaled to ``(0, 1)`` by ``rank / (n_obs + 1)``.

    Original: copula/dependogram.m
    """
    data = np.asarray(data, dtype=float)
    n = data.shape[0]
    ranks = np.apply_along_axis(rankdata, 0, data)
    return ranks / (n + 1.0)

dilog_function(x)

integral_x^1 log(t) / (1 - t) dt (related to the dilogarithm Li2). Not used by any closed form in this module (the original likewise never calls it from KendallCopula*.m/SpearmanCopula*.m) -- ported as a standalone special function for completeness.

Original: copula/diLogFunction.m

Source code in src/quanttoolbox/copula/dependence.py

def dilog_function(x: np.ndarray | float) -> np.ndarray:
    """``integral_x^1 log(t) / (1 - t) dt`` (related to the dilogarithm
    ``Li2``). Not used by any closed form in this module (the original
    likewise never calls it from `KendallCopula*.m`/`SpearmanCopula*.m`)
    -- ported as a standalone special function for completeness.

    Original: copula/diLogFunction.m
    """
    x_arr = np.atleast_1d(np.asarray(x, dtype=float))
    out = np.empty_like(x_arr)
    for i, xx in enumerate(x_arr.ravel()):
        integral, _ = quad(lambda t: np.log(t) / (1.0 - t), xx, 1.0)
        out.ravel()[i] = integral
    return out if np.asarray(x).ndim > 0 else out[0]

frank_rho(theta)

Spearman's rho of the Frank copula: 1 - 12*(D_1(theta) - D_2(theta)) / theta.

Original: copula/SpearmanCopulaFrank.m

Source code in src/quanttoolbox/copula/dependence.py

def frank_rho(theta: np.ndarray | float) -> np.ndarray:
    """Spearman's rho of the Frank copula: ``1 - 12*(D_1(theta) -
    D_2(theta)) / theta``.

    Original: copula/SpearmanCopulaFrank.m
    """
    theta = np.asarray(theta, dtype=float)
    d1 = debye_function(theta, 1)
    d2 = debye_function(theta, 2)
    return 1.0 - 12.0 * (d1 - d2) / theta

frank_tau(theta)

Kendall's tau of the Frank copula: 1 - 4*(1 - D_1(theta)) / theta (verified to match statsmodels).

Original: copula/KendallCopulaFrank.m

Source code in src/quanttoolbox/copula/dependence.py

def frank_tau(theta: np.ndarray | float) -> np.ndarray:
    """Kendall's tau of the Frank copula: ``1 - 4*(1 - D_1(theta)) /
    theta`` (verified to match `statsmodels`).

    Original: copula/KendallCopulaFrank.m
    """
    theta = np.asarray(theta, dtype=float)
    d1 = debye_function(theta, 1)
    return 1.0 - 4.0 * (1.0 - d1) / theta

gaussian_rho(rho)

Spearman's rho of the Gaussian copula: 6 * asin(rho/2) / pi.

Original: copula/SpearmanCopulaNormal.m

Source code in src/quanttoolbox/copula/dependence.py

def gaussian_rho(rho: np.ndarray | float) -> np.ndarray:
    """Spearman's rho of the Gaussian copula: ``6 * asin(rho/2) / pi``.

    Original: copula/SpearmanCopulaNormal.m
    """
    rho = np.asarray(rho, dtype=float)
    return 6.0 * np.arcsin(rho / 2.0) / np.pi

gaussian_tau(rho)

Kendall's tau of the Gaussian copula: 2 * asin(rho) / pi (verified to match statsmodels).

Original: copula/KendallCopulaNormal.m

Source code in src/quanttoolbox/copula/dependence.py

def gaussian_tau(rho: np.ndarray | float) -> np.ndarray:
    """Kendall's tau of the Gaussian copula: ``2 * asin(rho) / pi``
    (verified to match `statsmodels`).

    Original: copula/KendallCopulaNormal.m
    """
    rho = np.asarray(rho, dtype=float)
    return 2.0 * np.arcsin(rho) / np.pi

gumbel_rho(theta)

Spearman's rho of the Gumbel copula (no closed form -- numeric integration via spearman_rho_numeric).

Original: copula/SpearmanCopulaGumbel.m

Source code in src/quanttoolbox/copula/dependence.py

def gumbel_rho(theta: float) -> float:
    """Spearman's rho of the Gumbel copula (no closed form -- numeric
    integration via `spearman_rho_numeric`).

    Original: copula/SpearmanCopulaGumbel.m
    """
    from quanttoolbox.copula.families import gumbel_cdf

    return spearman_rho_numeric(lambda u1, u2: float(gumbel_cdf(u1, u2, theta)))

gumbel_tau(theta)

Kendall's tau of the Gumbel copula: 1 - 1/theta (verified to match statsmodels).

Original: copula/KendallCopulaGumbel.m

Source code in src/quanttoolbox/copula/dependence.py

def gumbel_tau(theta: np.ndarray | float) -> np.ndarray:
    """Kendall's tau of the Gumbel copula: ``1 - 1/theta`` (verified to
    match `statsmodels`).

    Original: copula/KendallCopulaGumbel.m
    """
    theta = np.asarray(theta, dtype=float)
    return 1.0 - 1.0 / theta

spearman_rho_numeric(cdf_fn)

Spearman's rho of an arbitrary bivariate copula cdf_fn(u1, u2), via the generic identity rho = 12 * integral_0^1 integral_0^1 C(u1, u2) du1 du2 - 3. Works for any copula CDF, not just the ones with a Python wrapper below -- see module docstring.

Original: copula/SpearmanCopula.m

Source code in src/quanttoolbox/copula/dependence.py

def spearman_rho_numeric(cdf_fn: Callable[[float, float], float]) -> float:
    """Spearman's rho of an arbitrary bivariate copula `cdf_fn(u1, u2)`,
    via the generic identity ``rho = 12 * integral_0^1 integral_0^1
    C(u1, u2) du1 du2 - 3``. Works for any copula CDF, not just the ones
    with a Python wrapper below -- see module docstring.

    Original: copula/SpearmanCopula.m
    """
    integral, _ = dblquad(lambda u2, u1: cdf_fn(u1, u2), 0.0, 1.0, 0.0, 1.0)
    return 12.0 * integral - 3.0

copula.simulate

Python alternatives

Keep -- no general-purpose Python library exposes generic conditional-CDF-inversion copula simulation the way this module and its MATLAB source do. The original rndCopula2.m is itself a generic bivariate engine (draw u1, v2 ~ U(0,1), solve conditional_cdf(u1, u2) = v2 for u2 by bisection) that the original only ever plugged Gumbel into; simulate_from_conditional_cdf here (built on this package's already-vectorized quanttoolbox.optim.bisection.bisection) keeps that genericity explicit -- it works for any family exposing a conditional CDF (e.g. families.amh_conditional_cdf), not just Gumbel. AMH and Frank have closed-form quantile inversions in the original and use those directly instead (simulate_amh, simulate_frank) -- faster, and exact rather than iterative. simulate_gaussian_copula/simulate_student_copula merge the n-dimensional and bivariate-special-case MATLAB originals into one Cholesky-based function each.

quanttoolbox.copula.simulate

Simulation for the copula families in families.py: closed-form quantile inversion where the original derived one (AMH, Frank), a generic root-finding conditional-inversion engine otherwise (Gumbel, and any other bivariate family), direct Cholesky-based simulation for the Gaussian/Student copulas, and an empirical-quantile transform for turning uniform copula draws into draws matching an observed marginal sample (NORTA-style simulation).

Ported from HSF toolbox copula/{rndCopula2,rndCopulaAMH,rndCopulaFrank, rndCopulaGumbel,rndCopulaNormal,rndCopulaNormal2,rndCopulaStudent, rndCopulaStudent2,rndnCopula,rndCopulaEmpiricalQuantile}.m.

Architecture:

• rndCopula2.m is already, in the original, a generic bivariate copula simulator: given any conditional-CDF function handle cndCopula2(u1, u2) (= dC(u1,u2)/du1), it draws u1, v2 ~ U(0,1) and solves cndCopula2(u1, u2) = v2 for u2 by bisection. rndCopulaGumbel.m is exactly this generic engine, called with Gumbel's own conditional CDF plugged in, because Gumbel has no closed-form quantile inversion. That generic engine is simulate_from_conditional_cdf here, built on this package's own quanttoolbox.optim.bisection.bisection (already vectorized) rather than a second hand-rolled bisection loop -- and, like the original, it works for any bivariate family with a conditional-CDF function, not just Gumbel: pass e.g. families.amh_conditional_cdf to simulate from AMH this way instead of simulate_amh's closed-form inversion.
• AMH and Frank do have closed-form conditional-CDF inversions (rndCopulaAMH.m's quadratic-formula root, rndCopulaFrank.m's direct log-formula), so simulate_amh/simulate_frank use those directly -- faster and exact, no root-finding needed.
• rndCopulaNormal.m (n-dim) and rndCopulaNormal2.m (bivariate, a rho in (-1, 1) special case of the same Cholesky construction) are merged into one simulate_gaussian_copula; same for rndCopulaStudent.m/rndCopulaStudent2.m -> simulate_student_copula.

empirical_quantile_transform(u, x)

Transform uniform copula draws u (shape (n_samples, n_dim)) into draws from the empirical (linearly-interpolated) marginal distribution of an observed sample x -- either shape (n_ref,) (the same reference sample reused for every column of u) or shape (n_ref, n_dim) (one reference sample per column). NORTA-style simulation: pairs a copula's dependence structure with an arbitrary observed marginal instead of a named distribution family.

Original: copula/rndCopulaEmpiricalQuantile.m

Source code in src/quanttoolbox/copula/simulate.py

def empirical_quantile_transform(u: np.ndarray, x: np.ndarray) -> np.ndarray:
    """Transform uniform copula draws `u` (shape ``(n_samples, n_dim)``)
    into draws from the empirical (linearly-interpolated) marginal
    distribution of an observed sample `x` -- either shape ``(n_ref,)``
    (the same reference sample reused for every column of `u`) or shape
    ``(n_ref, n_dim)`` (one reference sample per column). NORTA-style
    simulation: pairs a copula's dependence structure with an arbitrary
    observed marginal instead of a named distribution family.

    Original: copula/rndCopulaEmpiricalQuantile.m
    """
    u = np.asarray(u, dtype=float)
    n_dim = u.shape[1]
    x = np.asarray(x, dtype=float)
    if x.ndim == 1:
        x = np.tile(x[:, None], (1, n_dim))

    n_ref = x.shape[0]
    y = np.zeros_like(u)
    for i in range(n_dim):
        z = np.sort(x[:, i])
        z_padded = np.concatenate(([z[0]], z, [z[-1]]))
        w = n_ref * u[:, i]
        wt = np.floor(w).astype(int)
        frac = w - wt
        wt = wt + 1  # 1-based index into z_padded, matching the original
        y[:, i] = z_padded[wt] + frac * (z_padded[wt + 1] - z_padded[wt])
    return y

simulate_amh(theta, n_samples, random_state=None)

Simulate from the AMH copula via its closed-form conditional-CDF inversion (the smaller root of a quadratic in u2).

Original: copula/rndCopulaAMH.m

Source code in src/quanttoolbox/copula/simulate.py

def simulate_amh(
    theta: float, n_samples: int, random_state: object = None
) -> tuple[np.ndarray, np.ndarray]:
    """Simulate from the AMH copula via its closed-form conditional-CDF
    inversion (the smaller root of a quadratic in `u2`).

    Original: copula/rndCopulaAMH.m
    """
    rng = np.random.default_rng(random_state)
    u1 = rng.random(n_samples)
    v = rng.random(n_samples)

    a = v * theta**2 * (1.0 - u1) ** 2 - theta
    b = 2.0 * theta * v * (1.0 - theta + theta * u1) * (1.0 - u1) - (1.0 - theta)
    c = v * (1.0 - theta + theta * u1) ** 2
    delta = b**2 - 4.0 * a * c

    u2 = (-b - np.sqrt(delta)) / (2.0 * a)
    return u1, u2

simulate_frank(theta, n_samples, random_state=None)

Simulate from the Frank copula via its closed-form conditional-CDF inversion.

Original: copula/rndCopulaFrank.m

Source code in src/quanttoolbox/copula/simulate.py

def simulate_frank(
    theta: float, n_samples: int, random_state: object = None
) -> tuple[np.ndarray, np.ndarray]:
    """Simulate from the Frank copula via its closed-form conditional-CDF
    inversion.

    Original: copula/rndCopulaFrank.m
    """
    rng = np.random.default_rng(random_state)
    u1 = rng.random(n_samples)
    v2 = rng.random(n_samples)

    u2 = (
        -np.log(1.0 + v2 * (np.exp(-theta) - 1.0) / (v2 + (1.0 - v2) * np.exp(-theta * u1))) / theta
    )
    return u1, u2

simulate_from_conditional_cdf(conditional_cdf, n_samples, random_state=None, config=None)

Simulate n_samples draws (u1, u2) from a bivariate copula given only its conditional CDF conditional_cdf(u1, u2) = dC(u1, u2) / du1: draws u1, v2 ~ U(0,1) and solves conditional_cdf(u1, u2) = v2 for u2 by bisection. Works for any family exposing a conditional CDF (families.amh_conditional_cdf, families.gaussian_copula_conditional_cdf, a hand-written one for a family not in families.py at all, ...), not just Gumbel -- see module docstring.

Original: copula/rndCopula2.m (generic engine) and rndCopulaGumbel.m (Gumbel plugged into it)

Source code in src/quanttoolbox/copula/simulate.py

def simulate_from_conditional_cdf(
    conditional_cdf: Callable[[np.ndarray, np.ndarray], np.ndarray],
    n_samples: int,
    random_state: object = None,
    config: BisectionConfig | None = None,
) -> tuple[np.ndarray, np.ndarray]:
    """Simulate `n_samples` draws `(u1, u2)` from a bivariate copula given
    only its conditional CDF ``conditional_cdf(u1, u2) = dC(u1, u2) /
    du1``: draws `u1, v2 ~ U(0,1)` and solves ``conditional_cdf(u1, u2) =
    v2`` for `u2` by bisection. Works for any family exposing a
    conditional CDF (`families.amh_conditional_cdf`,
    `families.gaussian_copula_conditional_cdf`, a hand-written one for a
    family not in `families.py` at all, ...), not just Gumbel -- see
    module docstring.

    Original: copula/rndCopula2.m (generic engine) and rndCopulaGumbel.m
    (Gumbel plugged into it)
    """
    if config is None:
        config = BisectionConfig()
    rng = np.random.default_rng(random_state)
    u1 = rng.random(n_samples)
    v2 = rng.random(n_samples)

    eps = np.finfo(float).eps
    a = np.full(n_samples, eps)
    b = np.full(n_samples, 1.0 - eps)
    u2 = bisection(lambda u: conditional_cdf(u1, u) - v2, a, b, config)
    return u1, np.asarray(u2)

simulate_gaussian_copula(corr, n_samples, random_state=None)

Simulate from the Gaussian copula with correlation matrix corr (k_dim x k_dim, k_dim >= 2): Cholesky-correlate standard normal draws, then map through the standard normal CDF.

Original: copula/rndCopulaNormal.m (n-dim) and rndCopulaNormal2.m (bivariate) -- merged, see module docstring.

Source code in src/quanttoolbox/copula/simulate.py

def simulate_gaussian_copula(
    corr: np.ndarray, n_samples: int, random_state: object = None
) -> np.ndarray:
    """Simulate from the Gaussian copula with correlation matrix `corr`
    (``k_dim x k_dim``, ``k_dim >= 2``): Cholesky-correlate standard
    normal draws, then map through the standard normal CDF.

    Original: copula/rndCopulaNormal.m (n-dim) and rndCopulaNormal2.m
    (bivariate) -- merged, see module docstring.
    """
    corr = np.asarray(corr, dtype=float)
    rng = np.random.default_rng(random_state)
    z = rng.standard_normal((n_samples, corr.shape[0]))
    correlated = z @ np.linalg.cholesky(corr).T
    return normdist.cdf(correlated)

simulate_gumbel(theta, n_samples, random_state=None, config=None)

Simulate from the Gumbel copula via simulate_from_conditional_cdf (no closed-form quantile inversion exists for Gumbel).

Original: copula/rndCopulaGumbel.m

Source code in src/quanttoolbox/copula/simulate.py

def simulate_gumbel(
    theta: float, n_samples: int, random_state: object = None, config: BisectionConfig | None = None
) -> tuple[np.ndarray, np.ndarray]:
    """Simulate from the Gumbel copula via `simulate_from_conditional_cdf`
    (no closed-form quantile inversion exists for Gumbel).

    Original: copula/rndCopulaGumbel.m
    """

    def conditional_cdf(u1: np.ndarray, u2: np.ndarray) -> np.ndarray:
        u1t = -np.log(u1)
        u2t = -np.log(u2)
        w = u1t**theta + u2t**theta
        beta = 1.0 / theta
        return np.exp(-(w**beta)) * (1.0 + (u2t / u1t) ** theta) ** (beta - 1.0) / u1

    return simulate_from_conditional_cdf(conditional_cdf, n_samples, random_state, config)

simulate_student_copula(corr, nu, n_samples, random_state=None)

Simulate from the Student-t copula with correlation matrix corr and nu degrees of freedom: Cholesky-correlate standard normal draws, scale by an independent chi-squared draw, then map through the Student-t CDF.

Original: copula/rndCopulaStudent.m (n-dim) and rndCopulaStudent2.m (bivariate) -- merged, see module docstring.

Source code in src/quanttoolbox/copula/simulate.py

def simulate_student_copula(
    corr: np.ndarray, nu: float, n_samples: int, random_state: object = None
) -> np.ndarray:
    """Simulate from the Student-t copula with correlation matrix `corr`
    and `nu` degrees of freedom: Cholesky-correlate standard normal
    draws, scale by an independent chi-squared draw, then map through
    the Student-t CDF.

    Original: copula/rndCopulaStudent.m (n-dim) and rndCopulaStudent2.m
    (bivariate) -- merged, see module docstring.
    """
    corr = np.asarray(corr, dtype=float)
    rng = np.random.default_rng(random_state)
    n = rng.standard_normal((n_samples, corr.shape[0]))
    chi2 = rng.chisquare(nu, size=n_samples)

    correlated = n @ np.linalg.cholesky(corr).T
    scaled = correlated / np.sqrt(chi2 / nu)[:, None]
    return tdist.cdf(scaled, nu)