Source code for statsmodels.distributions.copula.archimedean

# -*- coding: utf-8 -*-
"""
Created on Fri Jan 29 19:19:45 2021

Author: Josef Perktold
License: BSD-3

"""
import sys

import numpy as np
from scipy import stats, integrate, optimize

from . import transforms
from .copulas import Copula
from statsmodels.tools.rng_qrng import check_random_state


def _debye(alpha):
    EPSILON = np.finfo(np.float32).eps

    def integrand(t):
        return t / (np.exp(t) - 1)

    debye_value = integrate.quad(integrand, EPSILON, alpha)[0] / alpha
    return debye_value


[docs]class ArchimedeanCopula(Copula): """Base class for Archimedean copulas Parameters ---------- transform : instance of transformation class Archimedean generator with required methods including first and second derivatives args : tuple Optional copula parameters. Copula parameters can be either provided when creating the instance or as arguments when calling methods. k_dim : int Dimension, number of components in the multivariate random variable. Currently only bivariate copulas are verified. Support for more than 2 dimension is incomplete. """ def __init__(self, transform, args=(), k_dim=2): super().__init__(k_dim=k_dim) self.args = args self.transform = transform self.k_args = 1 def _handle_args(self, args): # TODO: how to we handle non-tuple args? two we allow single values? # Model fit might give an args that can be empty if isinstance(args, np.ndarray): args = tuple(args) # handles empty arrays, unpacks otherwise if not isinstance(args, tuple): # could still be a scalar or numpy scalar args = (args,) if len(args) == 0 or args == (None,): # second condition because we converted None to tuple args = self.args return args
[docs] def cdf(self, u, args=()): """Evaluate cdf of Archimedean copula.""" args = self._handle_args(args) axis = -1 phi = self.transform.evaluate phi_inv = self.transform.inverse cdfv = phi_inv(phi(u, *args).sum(axis), *args) # clip numerical noise out = cdfv if isinstance(cdfv, np.ndarray) else None cdfv = np.clip(cdfv, 0., 1., out=out) # inplace if possible return cdfv
[docs] def pdf(self, u, args=()): """Evaluate pdf of Archimedean copula.""" args = self._handle_args(args) axis = -1 u = np.asarray(u) if u.shape[-1] > 2: msg = "pdf is currently only available for bivariate copula" raise ValueError(msg) # phi = self.transform.evaluate # phi_inv = self.transform.inverse phi_d1 = self.transform.deriv phi_d2 = self.transform.deriv2 cdfv = self.cdf(u, args=args) pdfv = - np.product(phi_d1(u, *args), axis) pdfv *= phi_d2(cdfv, *args) pdfv /= phi_d1(cdfv, *args)**3 return pdfv
[docs] def logpdf(self, u, args=()): """Evaluate log pdf of multivariate Archimedean copula.""" # TODO: replace by formulas, and exp in pdf args = self._handle_args(args) axis = -1 u = np.asarray(u) if u.shape[-1] > 2: msg = "pdf is currently only available for bivariate copula" raise ValueError(msg) phi_d1 = self.transform.deriv phi_d2 = self.transform.deriv2 cdfv = self.cdf(u, args=args) # I need np.abs because derivatives are negative, # is this correct for mv? logpdfv = np.sum(np.log(np.abs(phi_d1(u, *args))), axis) logpdfv += np.log(np.abs(phi_d2(cdfv, *args) / phi_d1(cdfv, *args)**3)) return logpdfv
def _arg_from_tau(self, tau): # for generic compat return self.theta_from_tau(tau)
[docs]class ClaytonCopula(ArchimedeanCopula): r"""Clayton copula. Dependence is greater in the negative tail than in the positive. .. math:: C_\theta(u,v) = \left[ \max\left\{ u^{-\theta} + v^{-\theta} -1 ; 0 \right\} \right]^{-1/\theta} with :math:`\theta\in[-1,\infty)\backslash\{0\}`. """ def __init__(self, theta=None, k_dim=2): if theta is not None: args = (theta,) else: args = () super().__init__(transforms.TransfClayton(), args=args, k_dim=k_dim) if theta is not None: if theta <= -1 or theta == 0: raise ValueError('Theta must be > -1 and !=0') self.theta = theta
[docs] def rvs(self, nobs=1, args=(), random_state=None): if self.k_dim != 2: msg = "rvs is only available for bivariate copula" raise NotImplementedError(msg) rng = check_random_state(random_state) th, = self._handle_args(args) x = rng.random((nobs, self.k_dim)) v = stats.gamma(1. / th).rvs(size=(nobs, 1), random_state=rng) return (1 - np.log(x) / v) ** (-1. / th)
[docs] def pdf(self, u, args=()): u = np.atleast_2d(u) th, = self._handle_args(args) a = (th + 1) * np.prod(u, axis=1) ** -(th + 1) b = np.sum(u ** -th, axis=1) - 1 c = -(2 * th + 1) / th return a * b ** c
[docs] def logpdf(self, u, args=()): # we skip Archimedean logpdf, that uses numdiff return super(ArchimedeanCopula, self).logpdf(u, args=args)
[docs] def cdf(self, u, args=()): u = np.atleast_2d(u) th, = self._handle_args(args) return (np.sum(u ** (-th), axis=1) - 1) ** (-1.0 / th)
[docs] def tau(self, theta=None): # Joe 2014 p. 168 if theta is None: theta = self.theta return theta / (theta + 2)
[docs] def theta_from_tau(self, tau): return 2 * tau / (1 - tau)
[docs]class FrankCopula(ArchimedeanCopula): r"""Frank copula. Dependence is symmetric. .. math:: C_\theta(\mathbf{u}) = -\frac{1}{\theta} \log \left[ 1- \frac{ \prod_j (1-\exp(- \theta u_j)) }{ (1 - \exp(-\theta)-1)^{d - 1} } \right] with :math:`\theta\in \mathbb{R}\backslash\{0\}, \mathbf{u} \in [0, 1]^d`. """ def __init__(self, theta=None, k_dim=2): if theta is not None: args = (theta,) else: args = () super().__init__(transforms.TransfFrank(), args=args, k_dim=k_dim) if theta is not None: if theta == 0: raise ValueError('Theta must be !=0') self.theta = theta
[docs] def rvs(self, nobs=1, args=(), random_state=None): if self.k_dim != 2: msg = "rvs is only available for bivariate copula" raise NotImplementedError(msg) rng = check_random_state(random_state) th, = self._handle_args(args) x = rng.random((nobs, self.k_dim)) v = stats.logser.rvs(1. - np.exp(-th), size=(nobs, 1), random_state=rng) return -1. / th * np.log(1. + np.exp(-(-np.log(x) / v)) * (np.exp(-th) - 1.))
# explicit BV formulas copied from Joe 1997 p. 141 # todo: check expm1 and log1p for improved numerical precision
[docs] def pdf(self, u, args=()): u = np.atleast_2d(u) th, = self._handle_args(args) if u.shape[-1] != 2: return super().pdf(u) g_ = np.exp(-th * np.sum(u, axis=1)) - 1 g1 = np.exp(-th) - 1 num = -th * g1 * (1 + g_) aux = np.prod(np.exp(-th * u) - 1, axis=1) + g1 den = aux ** 2 return num / den
[docs] def cdf(self, u, args=()): u = np.atleast_2d(u) th, = self._handle_args(args) dim = u.shape[-1] if dim != 2: return super().cdf(u) num = np.prod(1 - np.exp(- th * u), axis=1) den = (1 - np.exp(-th)) ** (dim - 1) return -1.0 / th * np.log(1 - num / den)
[docs] def logpdf(self, u, args=()): u = np.atleast_2d(u) th, = self._handle_args(args) if u.shape[-1] == 2: # bivariate case u1, u2 = u[..., 0], u[..., 1] b = 1 - np.exp(-th) pdf = np.log(th * b) - th * (u1 + u2) pdf -= 2 * np.log(b - (1 - np.exp(- th * u1)) * (1 - np.exp(- th * u2))) return pdf else: # for now use generic from base Copula class, log(self.pdf(...)) # we skip Archimedean logpdf, that uses numdiff super(ArchimedeanCopula, self).logpdf(u, args)
[docs] def cdfcond_2g1(self, u, args=()): """Conditional cdf of second component given the value of first. """ u = np.atleast_2d(u) th, = self._handle_args(args) if u.shape[-1] == 2: # bivariate case u1, u2 = u[..., 0], u[..., 1] cdfc = np.exp(- th * u1) cdfc /= np.expm1(-th) / np.expm1(- th * u2) + np.expm1(- th * u1) return cdfc else: raise NotImplementedError("u needs to be bivariate (2 columns)")
[docs] def ppfcond_2g1(self, q, u1, args=()): """Conditional pdf of second component given the value of first. """ u1 = np.asarray(u1) th, = self._handle_args(args) if u1.shape[-1] == 1: # bivariate case, conditional on value of first variable ppfc = - np.log(1 + np.expm1(- th) / ((1 / q - 1) * np.exp(-th * u1) + 1)) / th return ppfc else: raise NotImplementedError("u needs to be bivariate (2 columns)")
[docs] def tau(self, theta=None): # Joe 2014 p. 166 if theta is None: theta = self.theta debye_value = _debye(theta) return 1 + 4 * (debye_value - 1) / theta
[docs] def theta_from_tau(self, tau): MIN_FLOAT_LOG = np.log(sys.float_info.min) MAX_FLOAT_LOG = np.log(sys.float_info.max) def _theta_from_tau(alpha): return self.tau(theta=alpha) - tau result = optimize.least_squares(_theta_from_tau, 1, bounds=( MIN_FLOAT_LOG, MAX_FLOAT_LOG)) theta = result.x[0] return theta
[docs]class GumbelCopula(ArchimedeanCopula): r"""Gumbel copula. Dependence is greater in the positive tail than in the negative. .. math:: C_\theta(u,v) = \exp\!\left[ -\left( (-\log(u))^\theta + (-\log(v))^\theta \right)^{1/\theta} \right] with :math:`\theta\in[1,\infty)`. """ def __init__(self, theta=None, k_dim=2): if theta is not None: args = (theta,) else: args = () super().__init__(transforms.TransfGumbel(), args=args, k_dim=k_dim) if theta is not None: if theta <= 1: raise ValueError('Theta must be > 1') self.theta = theta
[docs] def rvs(self, nobs=1, args=(), random_state=None): if self.k_dim != 2: msg = "rvs is only available for bivariate copula" raise NotImplementedError(msg) rng = check_random_state(random_state) th, = self._handle_args(args) x = rng.random((nobs, self.k_dim)) v = stats.levy_stable.rvs( 1. / th, 1., 0, np.cos(np.pi / (2 * th)) ** th, size=(nobs, 1), random_state=rng ) return np.exp(-(-np.log(x) / v) ** (1. / th))
[docs] def pdf(self, u, args=()): u = np.atleast_2d(u) th, = self._handle_args(args) xy = -np.log(u) xy_theta = xy ** th sum_xy_theta = np.sum(xy_theta, axis=1) sum_xy_theta_theta = sum_xy_theta ** (1.0 / th) a = np.exp(-sum_xy_theta_theta) b = sum_xy_theta_theta + th - 1.0 c = sum_xy_theta ** (1.0 / th - 2) d = np.prod(xy, axis=1) ** (th - 1.0) e = np.prod(u, axis=1) ** (- 1.0) return a * b * c * d * e
[docs] def cdf(self, u, args=()): u = np.atleast_2d(u) th, = self._handle_args(args) h = np.sum((-np.log(u)) ** th, axis=1) cdf = np.exp(-h ** (1.0 / th)) return cdf
[docs] def logpdf(self, u, args=()): # we skip Archimedean logpdf, that uses numdiff return super(ArchimedeanCopula, self).logpdf(u, args=args)
[docs] def tau(self, theta=None): # Joe 2014 p. 172 if theta is None: theta = self.theta return (theta - 1) / theta
[docs] def theta_from_tau(self, tau): return 1 / (1 - tau)