Source code for arch.univariate.distribution

"""
Distributions to use in ARCH models.  All distributions must inherit from
:class:`Distribution` and provide the same methods with the same inputs.
"""

from abc import ABCMeta, abstractmethod
from collections.abc import Sequence
from typing import Callable, Optional, Union
import warnings

import numpy as np
from numpy import (
    abs,
    array,
    asarray,
    empty,
    exp,
    float64,
    int64,
    integer,
    isscalar,
    log,
    nan,
    ndarray,
    ones_like,
    pi,
    sign,
    sqrt,
    sum,
)
from numpy.random import Generator, RandomState, default_rng
from scipy.special import comb, gamma, gammainc, gammaincc, gammaln
import scipy.stats as stats

from arch.typing import ArrayLike, ArrayLike1D, Float64Array, Float64Array1D
from arch.utility.array import AbstractDocStringInheritor, ensure1d, to_array_1d

__all__ = ["Distribution", "Normal", "StudentsT", "SkewStudent", "GeneralizedError"]


[docs] class Distribution(metaclass=ABCMeta): """ Template for subclassing only """ def __init__( self, random_state: Optional[RandomState] = None, *, seed: Union[int, RandomState, Generator, None] = None, ) -> None: self._name = "Distribution" self.num_params: int = 0 self._parameters: Optional[Float64Array1D] = None if random_state is not None: if seed is not None: raise ValueError( "seed cannot be simultaneously used with random_state. Use " "seed to future proof your code." ) if not isinstance(random_state, RandomState): raise TypeError( "random_state must contain a RandomState instance when not None." ) warnings.warn( "random_state is deprecated and will be removed in a future update. " "Use seed instead.", FutureWarning, ) _seed: Union[int, RandomState, Generator, None] = random_state else: _seed = seed if _seed is None: self._generator: Union[Generator, RandomState] = default_rng() elif isinstance(_seed, (int, integer)): self._generator = default_rng(_seed) elif isinstance(_seed, (RandomState, Generator)): self._generator = _seed else: raise TypeError( "seed must by a NumPy Generator or RandomState or int, if not None." ) @property def name(self) -> str: """The name of the distribution""" return self._name def _check_constraints( self, parameters: Union[Sequence[float], ArrayLike1D, None] ) -> Float64Array1D: bounds = self.bounds(empty(0)) if parameters is not None: params = asarray(ensure1d(parameters, "parameters", False)) nparams = len(params) else: nparams = 0 params = empty(0) if nparams != len(bounds): raise ValueError(f"parameters must have {len(bounds)} elements") if len(bounds) == 0: return empty(0) for p, n, b in zip(params, self.name, bounds): if not (b[0] <= p <= b[1]): raise ValueError( f"{n} does not satisfy the bounds requirement of ({b[0]}, {b[1]})" ) return asarray(params) @property def generator(self) -> Union[RandomState, Generator]: """The NumPy Generator or RandomState attached to the distribution""" return self._generator @property def random_state(self) -> Union[RandomState, Generator]: """ The NumPy RandomState attached to the distribution .. deprecated:: 5.0 random_state is deprecated. Use generator instead. """ warnings.warn( "random_state is deprecated and will be removed in a future " "update. Use generator instead.", FutureWarning, ) return self._generator @abstractmethod def _simulator(self, size: Union[int, tuple[int, ...]]) -> Float64Array: """ Simulate i.i.d. draws from the distribution Parameters ---------- size : int or tuple Shape of the draws from the distribution Returns ------- rvs : ndarray Simulated pseudo random variables Notes ----- Must call `simulate` before using `_simulator` """
[docs] @abstractmethod def simulate( self, parameters: Union[int, float, Sequence[Union[float, int]], ArrayLike1D] ) -> Callable[[Union[int, tuple[int, ...]]], Float64Array]: """ Simulates i.i.d. draws from the distribution Parameters ---------- parameters : ndarray Distribution parameters Returns ------- simulator : callable Callable that take a single output size argument and returns i.i.d. draws from the distribution """
[docs] @abstractmethod def constraints(self) -> tuple[Float64Array, Float64Array]: """ Construct arrays to use in constrained optimization. Returns ------- A : ndarray Constraint loadings b : ndarray Constraint values Notes ----- Parameters satisfy the constraints A.dot(parameters)-b >= 0 """
[docs] @abstractmethod def bounds(self, resids: ArrayLike1D) -> list[tuple[float, float]]: """ Parameter bounds for use in optimization. Parameters ---------- resids : ndarray Residuals to use when computing the bounds Returns ------- bounds : list List containing a single tuple with (lower, upper) bounds """
[docs] @abstractmethod def loglikelihood( self, parameters: Union[Sequence[float], ArrayLike1D], resids: ArrayLike, sigma2: ArrayLike, individual: bool = False, ) -> Union[float, Float64Array1D]: """ Loglikelihood evaluation. Parameters ---------- parameters : ndarray Distribution shape parameters resids : ndarray nobs array of model residuals sigma2 : ndarray nobs array of conditional variances individual : bool, optional Flag indicating whether to return the vector of individual log likelihoods (True) or the sum (False) Notes ----- Returns the loglikelihood where resids are the "data", and parameters and sigma2 are inputs. """
[docs] @abstractmethod def starting_values(self, std_resid: ArrayLike1D) -> Float64Array1D: """ Construct starting values for use in optimization. Parameters ---------- std_resid : ndarray Estimated standardized residuals to use in computing starting values for the shape parameter Returns ------- sv : ndarray The estimated shape parameters for the distribution Notes ----- Size of sv depends on the distribution """
[docs] @abstractmethod def parameter_names(self) -> list[str]: """ Names of distribution shape parameters Returns ------- names : list (str) Parameter names """
[docs] @abstractmethod def ppf( self, pits: Union[float, Sequence[float], ArrayLike1D], parameters: Union[Sequence[float], ArrayLike1D, None] = None, ) -> Union[float, Float64Array]: """ Inverse cumulative density function (ICDF) Parameters ---------- pits : {float, ndarray} Probability-integral-transformed values in the interval (0, 1). parameters : ndarray, optional Distribution parameters. Use ``None`` for parameterless distributions. Returns ------- i : {float, ndarray} Inverse CDF values """
[docs] @abstractmethod def cdf( self, resids: Union[Sequence[float], ArrayLike1D], parameters: Union[Sequence[float], ArrayLike1D, None] = None, ) -> Float64Array: """ Cumulative distribution function Parameters ---------- resids : ndarray Values at which to evaluate the cdf parameters : ndarray Distribution parameters. Use ``None`` for parameterless distributions. Returns ------- f : ndarray CDF values """
[docs] @abstractmethod def moment( self, n: int, parameters: Union[Sequence[float], ArrayLike1D, None] = None ) -> float: """ Moment of order n Parameters ---------- n : int Order of moment parameters : ndarray, optional Distribution parameters. Use None for parameterless distributions. Returns ------- float Calculated moment """
[docs] @abstractmethod def partial_moment( self, n: int, z: float = 0.0, parameters: Union[Sequence[float], ArrayLike1D, None] = None, ) -> float: r""" Order n lower partial moment from -inf to z Parameters ---------- n : int Order of partial moment z : float, optional Upper bound for partial moment integral parameters : ndarray, optional Distribution parameters. Use None for parameterless distributions. Returns ------- float Partial moment References ---------- .. [1] Winkler et al. (1972) "The Determination of Partial Moments" *Management Science* Vol. 19 No. 3 Notes ----- The order n lower partial moment to z is .. math:: \int_{-\infty}^{z}x^{n}f(x)dx See [1]_ for more details. """
def __str__(self) -> str: return self._description() def __repr__(self) -> str: return self.__str__() + ", id: " + hex(id(self)) + "" def _description(self) -> str: return self.name + " distribution"
[docs] class Normal(Distribution, metaclass=AbstractDocStringInheritor): """ Standard normal distribution for use with ARCH models Parameters ---------- random_state : RandomState, optional .. deprecated:: 5.0 random_state is deprecated. Use seed instead. seed : {int, Generator, RandomState}, optional Random number generator instance or int to use. Set to ensure reproducibility. If using an int, the argument is passed to ``np.random.default_rng``. If not provided, ``default_rng`` is used with system-provided entropy. """ def __init__( self, random_state: Optional[RandomState] = None, *, seed: Union[None, int, RandomState, Generator] = None, ) -> None: super().__init__(random_state=random_state, seed=seed) self._name = "Normal"
[docs] def constraints(self) -> tuple[Float64Array, Float64Array]: return empty(0), empty(0)
[docs] def bounds(self, resids: ArrayLike1D) -> list[tuple[float, float]]: return []
[docs] def loglikelihood( self, parameters: Union[Sequence[float], ArrayLike1D], resids: ArrayLike, sigma2: ArrayLike, individual: bool = False, ) -> Union[float, Float64Array1D]: r"""Computes the log-likelihood of assuming residuals are normally distributed, conditional on the variance Parameters ---------- parameters : ndarray The normal likelihood has no shape parameters. Empty since the standard normal has no shape parameters. resids : ndarray The residuals to use in the log-likelihood calculation sigma2 : ndarray Conditional variances of resids individual : bool, optional Flag indicating whether to return the vector of individual log likelihoods (True) or the sum (False) Returns ------- ll : float The log-likelihood Notes ----- The log-likelihood of a single data point x is .. math:: \ln f\left(x\right)=-\frac{1}{2}\left(\ln2\pi+\ln\sigma^{2} +\frac{x^{2}}{\sigma^{2}}\right) """ lls = -0.5 * (log(2 * pi) + log(sigma2) + resids**2.0 / sigma2) if individual: return lls else: return sum(lls)
[docs] def starting_values(self, std_resid: ArrayLike1D) -> Float64Array1D: return empty(0)
def _simulator(self, size: Union[int, tuple[int, ...]]) -> Float64Array: return self._generator.standard_normal(size)
[docs] def simulate( self, parameters: Union[int, float, Sequence[Union[float, int]], ArrayLike1D] ) -> Callable[[Union[int, tuple[int, ...]]], Float64Array]: return self._simulator
[docs] def parameter_names(self) -> list[str]: return []
[docs] def cdf( self, resids: Union[Sequence[float], ArrayLike1D], parameters: Union[Sequence[float], ArrayLike1D, None] = None, ) -> Float64Array: self._check_constraints(parameters) return stats.norm.cdf(asarray(resids))
[docs] def ppf( self, pits: Union[float, Sequence[float], ArrayLike1D], parameters: Union[Sequence[float], ArrayLike1D, None] = None, ) -> Float64Array: self._check_constraints(parameters) scalar = isscalar(pits) if scalar: assert isinstance(pits, float) _pits = array([pits], dtype=np.float64) else: _pits = asarray(pits, dtype=float) ppf = stats.norm.ppf(_pits) if scalar: return ppf[0] else: return ppf
[docs] def moment( self, n: int, parameters: Union[Sequence[float], ArrayLike1D, None] = None ) -> float: if n < 0: return nan return stats.norm.moment(n)
[docs] def partial_moment( self, n: int, z: float = 0.0, parameters: Union[Sequence[float], ArrayLike1D, None] = None, ) -> float: if n < 0: return nan elif n == 0: return stats.norm.cdf(z) elif n == 1: return -stats.norm.pdf(z) else: return -(z ** (n - 1)) * stats.norm.pdf(z) + (n - 1) * self.partial_moment( n - 2, z, parameters )
[docs] class StudentsT(Distribution, metaclass=AbstractDocStringInheritor): """ Standardized Student's distribution for use with ARCH models Parameters ---------- random_state : RandomState, optional .. deprecated:: 5.0 random_state is deprecated. Use seed instead. seed : {int, Generator, RandomState}, optional Random number generator instance or int to use. Set to ensure reproducibility. If using an int, the argument is passed to ``np.random.default_rng``. If not provided, ``default_rng`` is used with system-provided entropy. """ def __init__( self, random_state: Optional[RandomState] = None, *, seed: Union[None, int, RandomState, Generator] = None, ) -> None: super().__init__(random_state=random_state, seed=seed) self._name = "Standardized Student's t" self.num_params: int = 1
[docs] def constraints(self) -> tuple[Float64Array, Float64Array]: return array([[1], [-1]]), array([2.05, -500.0])
[docs] def bounds(self, resids: ArrayLike1D) -> list[tuple[float, float]]: return [(2.05, 500.0)]
[docs] def loglikelihood( self, parameters: Union[Sequence[float], ArrayLike1D], resids: ArrayLike, sigma2: ArrayLike, individual: bool = False, ) -> Union[float, Float64Array1D]: r"""Computes the log-likelihood of assuming residuals are have a standardized (to have unit variance) Student's t distribution, conditional on the variance. Parameters ---------- parameters : ndarray Shape parameter of the t distribution resids : ndarray The residuals to use in the log-likelihood calculation sigma2 : ndarray Conditional variances of resids individual : bool, optional Flag indicating whether to return the vector of individual log likelihoods (True) or the sum (False) Returns ------- ll : float The log-likelihood Notes ----- The log-likelihood of a single data point x is .. math:: \ln\Gamma\left(\frac{\nu+1}{2}\right) -\ln\Gamma\left(\frac{\nu}{2}\right) -\frac{1}{2}\ln(\pi\left(\nu-2\right)\sigma^{2}) -\frac{\nu+1}{2}\ln(1+x^{2}/(\sigma^{2}(\nu-2))) where :math:`\Gamma` is the gamma function. """ nu = parameters[0] lls = gammaln((nu + 1) / 2) - gammaln(nu / 2) - log(pi * (nu - 2)) / 2 lls -= 0.5 * (log(sigma2)) lls -= ((nu + 1) / 2) * (log(1 + (resids**2.0) / (sigma2 * (nu - 2)))) if individual: return lls else: return sum(lls)
[docs] def starting_values(self, std_resid: ArrayLike1D) -> Float64Array1D: """ Construct starting values for use in optimization. Parameters ---------- std_resid : ndarray Estimated standardized residuals to use in computing starting values for the shape parameter Returns ------- sv : ndarray Array containing starting valuer for shape parameter Notes ----- Uses relationship between kurtosis and degree of freedom parameter to produce a moment-based estimator for the starting values. """ k = stats.kurtosis(std_resid, fisher=False) sv = max((4.0 * k - 6.0) / (k - 3.0) if k > 3.75 else 12.0, 4.0) return array([sv])
def _simulator(self, size: Union[int, tuple[int, ...]]) -> Float64Array: assert self._parameters is not None parameters = self._parameters std_dev = sqrt(parameters[0] / (parameters[0] - 2)) return self._generator.standard_t(self._parameters[0], size=size) / std_dev
[docs] def simulate( self, parameters: Union[int, float, Sequence[Union[float, int]], ArrayLike1D] ) -> Callable[[Union[int, tuple[int, ...]]], Float64Array]: parameters = ensure1d(parameters, "parameters", False) if parameters[0] <= 2.0: raise ValueError("The shape parameter must be larger than 2") self._parameters = ensure1d(parameters, "parameters", series=False).astype( float ) return self._simulator
[docs] def parameter_names(self) -> list[str]: return ["nu"]
[docs] def cdf( self, resids: Union[Sequence[float], ArrayLike1D], parameters: Union[Sequence[float], ArrayLike1D, None] = None, ) -> Float64Array: parameters = self._check_constraints(parameters) nu = parameters[0] var = nu / (nu - 2) return stats.t(nu, scale=1.0 / sqrt(var)).cdf(asarray(resids))
[docs] def ppf( self, pits: Union[float, Sequence[float], ArrayLike1D], parameters: Union[Sequence[float], ArrayLike1D, None] = None, ) -> Float64Array: parameters = self._check_constraints(parameters) _pits = asarray(pits, dtype=float) nu = parameters[0] var = nu / (nu - 2) return stats.t(nu, scale=1.0 / sqrt(var)).ppf(_pits)
[docs] def moment( self, n: int, parameters: Union[Sequence[float], ArrayLike1D, None] = None ) -> float: if n < 0: return nan parameters = self._check_constraints(parameters) nu = parameters[0] var = nu / (nu - 2) return stats.t.moment(n, nu, scale=1.0 / sqrt(var))
[docs] def partial_moment( self, n: int, z: float = 0.0, parameters: Union[Sequence[float], ArrayLike1D, None] = None, ) -> float: parameters = self._check_constraints(parameters) nu = parameters[0] var = nu / (nu - 2) scale = 1.0 / sqrt(var) moment = (scale**n) * self._ord_t_partial_moment(n, z / scale, nu) return moment
@staticmethod def _ord_t_partial_moment(n: int, z: float, nu: float) -> float: r""" Partial moments for ordinary parameterization of Students t df=nu Parameters ---------- n : int Order of partial moment z : float Upper bound for partial moment integral nu : float Degrees of freedom Returns ------- float Calculated moment References ---------- .. [1] Winkler et al. (1972) "The Determination of Partial Moments" *Management Science* Vol. 19 No. 3 Notes ----- The order n lower partial moment to z is .. math:: \int_{-\infty}^{z}x^{n}f(x)dx See [1]_ for more details. """ if n < 0 or n >= nu: return nan elif n == 0: moment = stats.t.cdf(z, nu) elif n == 1: c = gamma(0.5 * (nu + 1)) / (sqrt(nu * pi) * gamma(0.5 * nu)) e = 0.5 * (nu + 1) moment = (0.5 * (c * nu) / (1 - e)) * ((1 + (z**2) / nu) ** (1 - e)) else: t1 = (z ** (n - 1)) * (nu + z**2) * stats.t.pdf(z, nu) t2 = (n - 1) * nu * StudentsT._ord_t_partial_moment(n - 2, z, nu) moment = (1 / (n - nu)) * (t1 - t2) return moment
[docs] class SkewStudent(Distribution, metaclass=AbstractDocStringInheritor): r""" Standardized Skewed Student's distribution for use with ARCH models Parameters ---------- random_state : RandomState, optional .. deprecated:: 5.0 random_state is deprecated. Use seed instead. seed : {int, Generator, RandomState}, optional Random number generator instance or int to use. Set to ensure reproducibility. If using an int, the argument is passed to ``np.random.default_rng``. If not provided, ``default_rng`` is used with system-provided entropy. Notes ----- The Standardized Skewed Student's distribution ([1]_) takes two parameters, :math:`\eta` and :math:`\lambda`. :math:`\eta` controls the tail shape and is similar to the shape parameter in a Standardized Student's t. :math:`\lambda` controls the skewness. When :math:`\lambda=0` the distribution is identical to a standardized Student's t. References ---------- .. [1] Hansen, B. E. (1994). Autoregressive conditional density estimation. *International Economic Review*, 35(3), 705–730. <https://www.ssc.wisc.edu/~bhansen/papers/ier_94.pdf> """ def __init__( self, random_state: Optional[RandomState] = None, *, seed: Union[None, int, RandomState, Generator] = None, ) -> None: super().__init__(random_state=random_state, seed=seed) self._name = "Standardized Skew Student's t" self.num_params: int = 2
[docs] def constraints(self) -> tuple[Float64Array, Float64Array]: return array([[1, 0], [-1, 0], [0, 1], [0, -1]]), array([2.05, -300.0, -1, -1])
[docs] def bounds(self, resids: ArrayLike1D) -> list[tuple[float, float]]: return [(2.05, 300.0), (-1, 1)]
[docs] def loglikelihood( self, parameters: Union[Sequence[float], ArrayLike1D], resids: ArrayLike, sigma2: ArrayLike, individual: bool = False, ) -> Float64Array1D: r""" Computes the log-likelihood of assuming residuals are have a standardized (to have unit variance) Skew Student's t distribution, conditional on the variance. Parameters ---------- parameters : ndarray Shape parameter of the skew-t distribution resids : ndarray The residuals to use in the log-likelihood calculation sigma2 : ndarray Conditional variances of resids individual : bool, optional Flag indicating whether to return the vector of individual log likelihoods (True) or the sum (False) Returns ------- ll : float The log-likelihood Notes ----- The log-likelihood of a single data point x is .. math:: \ln\left[\frac{bc}{\sigma}\left(1+\frac{1}{\eta-2} \left(\frac{a+bx/\sigma} {1+sgn(x/\sigma+a/b)\lambda}\right)^{2}\right) ^{-\left(\eta+1\right)/2}\right], where :math:`2<\eta<\infty`, and :math:`-1<\lambda<1`. The constants :math:`a`, :math:`b`, and :math:`c` are given by .. math:: a=4\lambda c\frac{\eta-2}{\eta-1}, \quad b^{2}=1+3\lambda^{2}-a^{2}, \quad c=\frac{\Gamma\left(\frac{\eta+1}{2}\right)} {\sqrt{\pi\left(\eta-2\right)} \Gamma\left(\frac{\eta}{2}\right)}, and :math:`\Gamma` is the gamma function. """ parameters_arr: Float64Array1D parameters_arr = np.ravel(asarray(parameters, dtype=float)) eta, lam = parameters_arr const_c = self.__const_c(parameters_arr) const_a = self.__const_a(parameters_arr) const_b = self.__const_b(parameters_arr) resids = resids / sigma2**0.5 lls = log(const_b) + const_c - log(sigma2) / 2 if abs(lam) >= 1.0: lam = sign(lam) * (1.0 - 1e-6) llf_resid = ( (const_b * resids + const_a) / (1 + sign(resids + const_a / const_b) * lam) ) ** 2 lls -= (eta + 1) / 2 * log(1 + llf_resid / (eta - 2)) if individual: return lls else: return sum(lls)
[docs] def starting_values(self, std_resid: ArrayLike1D) -> Float64Array1D: """ Construct starting values for use in optimization. Parameters ---------- std_resid : ndarray Estimated standardized residuals to use in computing starting values for the shape parameter Returns ------- sv : ndarray Array containing starting valuer for shape parameter Notes ----- Uses relationship between kurtosis and degree of freedom parameter to produce a moment-based estimator for the starting values. """ k = stats.kurtosis(std_resid, fisher=False) sv = max((4.0 * k - 6.0) / (k - 3.0) if k > 3.75 else 12.0, 4.0) return array([sv, 0.0])
def _simulator(self, size: Union[int, tuple[int, ...]]) -> Float64Array: # No need to normalize since it is already done in parameterization assert self._parameters is not None if isinstance(self._generator, Generator): uniforms = self._generator.random(size=size) else: uniforms = self._generator.random_sample(size=size) ppf = self.ppf(uniforms, self._parameters) assert isinstance(ppf, ndarray) return ppf
[docs] def simulate( self, parameters: Union[int, float, Sequence[Union[float, int]], ArrayLike1D] ) -> Callable[[Union[int, tuple[int, ...]]], Float64Array]: parameters = ensure1d(parameters, "parameters", False) if parameters[0] <= 2.0: raise ValueError("The shape parameter must be larger than 2") if abs(parameters[1]) > 1.0: raise ValueError( "The skew parameter must be smaller than 1 in absolute value" ) self._parameters = ensure1d(parameters, "parameters", series=False).astype( float ) return self._simulator
[docs] def parameter_names(self) -> list[str]: return ["eta", "lambda"]
def __const_a(self, parameters: Float64Array1D) -> float: """ Compute a constant. Parameters ---------- parameters : ndarray Shape parameters of the skew-t distribution Returns ------- a : float Constant used in the distribution """ eta, lam = parameters c = self.__const_c(parameters) return float(4 * lam * exp(c) * (eta - 2) / (eta - 1)) def __const_b(self, parameters: Float64Array1D) -> float: """ Compute b constant. Parameters ---------- parameters : ndarray Shape parameters of the skew-t distribution Returns ------- b : float Constant used in the distribution """ lam = float(parameters[1]) a = self.__const_a(parameters) return (1 + 3 * lam**2 - a**2) ** 0.5 @staticmethod def __const_c(parameters: Float64Array1D) -> float: """ Compute c constant. Parameters ---------- parameters : ndarray Shape parameters of the skew-t distribution Returns ------- c : float Log of the constant used in loglikelihood """ eta = parameters[0] # return gamma((eta+1)/2) / ((pi*(eta-2))**.5 * gamma(eta/2)) return float( gammaln((eta + 1) / 2) - gammaln(eta / 2) - log(pi * (eta - 2)) / 2 )
[docs] def cdf( self, resids: Union[Sequence[float], ArrayLike1D], parameters: Union[Sequence[float], ArrayLike1D, None] = None, ) -> Float64Array: parameters = self._check_constraints(parameters) scalar = isscalar(resids) _resids = ensure1d(resids, "resids").astype(float) eta, lam = parameters a = self.__const_a(parameters) b = self.__const_b(parameters) var = eta / (eta - 2) y1 = (b * _resids + a) / (1 - lam) * sqrt(var) y2 = (b * _resids + a) / (1 + lam) * sqrt(var) tcdf = stats.t(eta).cdf p = (1 - lam) * tcdf(y1) * (_resids < (-a / b)) p += (_resids >= (-a / b)) * ((1 - lam) / 2 + (1 + lam) * (tcdf(y2) - 0.5)) if scalar: p = p[0] return p
[docs] def ppf( self, pits: Union[float, Sequence[float], Float64Array, ArrayLike], parameters: Union[Sequence[float], ArrayLike1D, None] = None, ) -> Union[float, Float64Array]: parameters = self._check_constraints(parameters) scalar = isscalar(pits) if scalar: pits = array([pits]) pits = asarray(pits, dtype=float) eta, lam = parameters a = self.__const_a(parameters) b = self.__const_b(parameters) cond = pits < (1 - lam) / 2 icdf1 = stats.t.ppf(pits[cond] / (1 - lam), eta) icdf2 = stats.t.ppf(0.5 + (pits[~cond] - (1 - lam) / 2) / (1 + lam), eta) icdf = -999.99 * ones_like(pits) assert isinstance(icdf, ndarray) icdf[cond] = icdf1 icdf[~cond] = icdf2 icdf = icdf * (1 + sign(pits - (1 - lam) / 2) * lam) * (1 - 2 / eta) ** 0.5 - a icdf = icdf / float64(b) if scalar: return float(icdf[0]) assert isinstance(icdf, ndarray) return icdf
[docs] def moment( self, n: int, parameters: Union[Sequence[float], ArrayLike1D, None] = None ) -> float: parameters = self._check_constraints(parameters) eta, lam = parameters if n < 0 or n >= eta: return nan a = self.__const_a(parameters) b = self.__const_b(parameters) loc = -a / b lscale = sqrt(1 - 2 / eta) * (1 - lam) / b rscale = sqrt(1 - 2 / eta) * (1 + lam) / b moment = 0.0 for k in range(n + 1): # binomial expansion around loc # 0->inf right partial moment for ordinary t(eta) r_pmom = ( 0.5 * (gamma(0.5 * (k + 1)) * gamma(0.5 * (eta - k)) * eta ** (0.5 * k)) / (sqrt(pi) * gamma(0.5 * eta)) ) l_pmom = ((-1) ** k) * r_pmom lhs = (1 - lam) * (lscale**k) * (loc ** (n - k)) * l_pmom rhs = (1 + lam) * (rscale**k) * (loc ** (n - k)) * r_pmom moment += comb(n, k) * (lhs + rhs) return moment
[docs] def partial_moment( self, n: int, z: float = 0.0, parameters: Union[Sequence[float], ArrayLike1D, None] = None, ) -> float: parameters_arr = self._check_constraints(parameters) eta, lam = parameters_arr if n < 0 or n >= eta: return nan a = self.__const_a(parameters_arr) b = self.__const_b(parameters_arr) loc = -a / b lscale = sqrt(1 - 2 / eta) * (1 - lam) / b rscale = sqrt(1 - 2 / eta) * (1 + lam) / b moment = 0.0 for k in range(n + 1): # binomial expansion around loc lbound = min(z, loc) lhs = ( (1 - lam) * (loc ** (n - k)) * (lscale**k) * StudentsT._ord_t_partial_moment( k, z=(lbound - loc) / lscale, nu=float(eta) ) ) if z > loc: rhs = ( (1 + lam) * (loc ** (n - k)) * (rscale**k) * ( StudentsT._ord_t_partial_moment( k, z=(z - loc) / rscale, nu=float(eta) ) - StudentsT._ord_t_partial_moment(k, z=0.0, nu=float(eta)) ) ) else: rhs = 0.0 moment += comb(n, k) * (lhs + rhs) return moment
[docs] class GeneralizedError(Distribution, metaclass=AbstractDocStringInheritor): """ Generalized Error distribution for use with ARCH models Parameters ---------- random_state : RandomState, optional .. deprecated:: 5.0 random_state is deprecated. Use seed instead. seed : {int, Generator, RandomState}, optional Random number generator instance or int to use. Set to ensure reproducibility. If using an int, the argument is passed to ``np.random.default_rng``. If not provided, ``default_rng`` is used with system-provided entropy. """ def __init__( self, random_state: Optional[RandomState] = None, *, seed: Union[None, int, RandomState, Generator] = None, ) -> None: super().__init__(random_state=random_state, seed=seed) self._name = "Generalized Error Distribution" self.num_params: int = 1
[docs] def constraints(self) -> tuple[Float64Array, Float64Array]: return array([[1], [-1]]), array([1.01, -500.0])
[docs] def bounds(self, resids: ArrayLike1D) -> list[tuple[float, float]]: return [(1.01, 500.0)]
[docs] def loglikelihood( self, parameters: Union[Sequence[float], ArrayLike1D], resids: ArrayLike, sigma2: ArrayLike, individual: bool = False, ) -> Float64Array1D: r""" Computes the log-likelihood of assuming residuals are have a Generalized Error Distribution, conditional on the variance. Parameters ---------- parameters : ndarray Shape parameter of the GED distribution resids : ndarray The residuals to use in the log-likelihood calculation sigma2 : ndarray Conditional variances of resids individual : bool, optional Flag indicating whether to return the vector of individual log likelihoods (True) or the sum (False) Returns ------- ll : float The log-likelihood Notes ----- The log-likelihood of a single data point x is .. math:: \ln\nu-\ln c-\ln\Gamma(\frac{1}{\nu})-(1+\frac{1}{\nu})\ln2 -\frac{1}{2}\ln\sigma^{2} -\frac{1}{2}\left|\frac{x}{c\sigma}\right|^{\nu} where :math:`\Gamma` is the gamma function and :math:`\ln c` is .. math:: \ln c=\frac{1}{2}\left(\frac{-2}{\nu}\ln2+\ln\Gamma(\frac{1}{\nu}) -\ln\Gamma(\frac{3}{\nu})\right). """ nu = parameters[0] log_c = 0.5 * (-2 / nu * log(2) + gammaln(1 / nu) - gammaln(3 / nu)) c = exp(log_c) lls = log(nu) - log_c - gammaln(1 / nu) - (1 + 1 / nu) * log(2) lls -= 0.5 * log(sigma2) lls -= 0.5 * abs(resids / (sqrt(sigma2) * c)) ** nu if individual: return lls else: return sum(lls)
[docs] def starting_values(self, std_resid: ArrayLike1D) -> Float64Array1D: """ Construct starting values for use in optimization. Parameters ---------- std_resid : ndarray Estimated standardized residuals to use in computing starting values for the shape parameter Returns ------- sv : ndarray Array containing starting valuer for shape parameter Notes ----- Defaults to 1.5 which is implies heavier tails than a normal """ return to_array_1d(array([1.5]))
def _simulator(self, size: Union[int, tuple[int, ...]]) -> Float64Array: assert self._parameters is not None parameters = self._parameters nu = parameters[0] randoms = self._generator.standard_gamma(1 / nu, size) ** (1.0 / nu) if isinstance(self._generator, Generator): random_ints = self._generator.integers(0, 2, size, dtype=int64) else: random_ints = self._generator.randint(0, 2, size, dtype=int64) randoms *= 2.0 * random_ints - 1.0 scale = sqrt(gamma(3.0 / nu) / gamma(1.0 / nu)) return randoms / scale
[docs] def simulate( self, parameters: Union[int, float, Sequence[Union[float, int]], ArrayLike1D] ) -> Callable[[Union[int, tuple[int, ...]]], Float64Array]: parameters = ensure1d(parameters, "parameters", False) if parameters[0] <= 1.0: raise ValueError("The shape parameter must be larger than 1") self._parameters = ensure1d(parameters, "parameters", series=False).astype( float ) return self._simulator
[docs] def parameter_names(self) -> list[str]: return ["nu"]
[docs] def ppf( self, pits: Union[float, Sequence[float], ArrayLike1D], parameters: Union[Sequence[float], ArrayLike1D, None] = None, ) -> Float64Array: parameters = self._check_constraints(parameters) _pits = asarray(pits, dtype=float) nu = parameters[0] var = stats.gennorm(nu).var() return stats.gennorm(nu, scale=1.0 / sqrt(var)).ppf(_pits)
[docs] def cdf( self, resids: Union[Sequence[float], ArrayLike1D], parameters: Union[Sequence[float], ArrayLike1D, None] = None, ) -> Float64Array: parameters = self._check_constraints(parameters) nu = parameters[0] var = stats.gennorm(nu).var() _resids = ensure1d(resids, "resids", series=False).astype(float) return stats.gennorm(nu, scale=1.0 / sqrt(var)).cdf(_resids)
[docs] def moment( self, n: int, parameters: Union[Sequence[float], ArrayLike1D, None] = None ) -> float: if n < 0: return nan parameters = self._check_constraints(parameters) nu = parameters[0] var = stats.gennorm(nu).var() return stats.gennorm.moment(n, nu, scale=1.0 / sqrt(var))
[docs] def partial_moment( self, n: int, z: float = 0.0, parameters: Union[Sequence[float], ArrayLike1D, None] = None, ) -> float: parameters = self._check_constraints(parameters) nu = parameters[0] scale = 1.0 / sqrt(stats.gennorm(nu).var()) moment = (scale**n) * self._ord_gennorm_partial_moment(n, z / scale, nu) return moment
@staticmethod def _ord_gennorm_partial_moment(n: int, z: int, beta: float) -> float: r""" Partial moment for ordinary generalized normal parameterization. Parameters ---------- n : int Order of partial moment z : float Upper bound for partial moment integral beta : float Parameter of generalized normal Returns ------- float Partial moment Notes ----- The standard parameterization follows: .. math:: f(x)=\frac{\beta}{2\Gamma(\beta^{-1})}\text{exp}\left(-|x|^{\beta}\right) """ if n < 0: return nan w = 0.5 * beta / gamma(1 / beta) # integral over (-inf, min(z,0)) lz = abs(min(z, 0)) ** beta lterm = ( w * ((-1) ** n) * (1 / beta) * gamma((n + 1) / beta) * gammaincc((n + 1) / beta, lz) ) # remaining integral rz = max(0, z) ** beta rterm = w * (1 / beta) * gamma((n + 1) / beta) * gammainc((n + 1) / beta, rz) moment = lterm + rterm return moment